颜色协调模型Color Harmoniztion

前言

最近做换脸，在肤色调整的那一块，看到一个有意思的文章，复现一波玩玩。不过最后一步掉链子了，有兴趣的可以一起讨论把链子补上。

主要是github上大佬的那个复现代码和原文有点差异，而且代码复杂度过高，阅读费劲，这里为了清晰理解理论知识，就一步一步按照论文的每个章节走，不过有很大一部分代码都借鉴大佬了。

国际惯例，参考文献：

论文《Color Harmonization》

opencv超像素分割

大佬的实现

简介

这篇论文主要干啥呢，看下图

注意小姑凉的衣服的颜色从左边的原图变成了右边的变化了的图片。

那么论文的工作就是：在尽量保持原色的基础上，调整图片颜色，使得图片整体颜色看起来更加协调。

基本原理就是将RGB图像转成HSV，然后调整色相(hue)直方图。下图中左边的彩色圆环代表原始图片的色相直方图，右边为论文调整后的色相直方图。

感觉这个论文用来去除杂色应该挺有用，不过更好玩的是能交互式自动调色。

【注】要知道色度hue的范围是$0°\sim360°$0°∼360° ，opencv里面提取出来的是$0°\sim 180°$0°∼180°

流程复现

环境python3、opencv-python==3.4.2.16、opencv-contrib-python==3.4.2.16

注意如果找不到ximgproc，手动卸载opencv-contrib-python再重新安装，自动卸载安装无法解决问题。

颜色模板

首先，文章介绍了8种比较适合用于调色的模板，注意灰色区域的大小固定，但是位置不是固定的，可以绕着圆心旋转。

当图像的色相都落在灰色区域，就是颜色处理的比较好的图像。

文末介绍了每个模板的灰色区域面积比：

• 模板中大的扇形区域：V,Y,X占26%的面积，角度是$93.6°$93.6°

• 模板中小的扇形区域：i,L,I,Y占5%的面积，角度是$18°$18°

• 模板L占22%面积，角度是$79.0°$79.0°

• 模板T占50%面积，角度是$180°$180°

• 模板I,X,Y的两个扇形区的夹角是$180°$180°

• 模板L两个扇形区的夹角是$90°$90°

不过大佬已经把这部分列出来了，两个数，第一个是扇形区域的中心，第二个是扇形区域的宽度：

#定义模板,分别定义的中心与边界偏转角度
HueTemplates = {
    "i"       : [( 0.00, 0.05)],
    "V"       : [( 0.00, 0.26)],
    "L"       : [( 0.00, 0.05), ( 0.25, 0.22)],
    "mirror_L": [( 0.00, 0.05), (-0.25, 0.22)],
    "I"       : [( 0.00, 0.05), ( 0.50, 0.05)],
    "T"       : [( 0.25, 0.50)],
    "Y"       : [( 0.00, 0.26), ( 0.50, 0.05)],
    "X"       : [( 0.00, 0.26), ( 0.50, 0.26)],
}

随便显示一个瞅瞅

#预览模板
def show_temp(template_name,template_alpha):
    canvas = np.zeros((canvas_h, canvas_w, 3))  #画布
    cv2.circle(canvas, (yc, xc), circle_r, (255,255,255), -1) #画圆
    for t in HueTemplates[template_name]:
        center = t[0]*360 + template_alpha
        width = t[1]*360
        start = center - width/2
        end = center + width/2
        cv2.ellipse(canvas,(yc,xc),(circle_r,circle_r),0,start,end,(0,0,0),-1,cv2.LINE_AA)
    cv2.circle(canvas, (yc, xc), 10, (0,0,0), -1,cv2.LINE_AA) #画中心点
    canvas = np.array(canvas,np.uint8)    
    canvas = cv2.cvtColor(canvas,cv2.COLOR_BGR2RGB)
    return canvas

plt.imsave('./temp/show_temp.png',show_temp('X',90),cmap='gray')

色相直方图

就是单纯的将图像的H通道提取出来，看看像素个数。

def count_hue_histogram(X):
    N = 360
    H = X[:, :, 0].astype(np.int32) * 2    
    H_flat = H.flatten()
    
    histo = np.zeros(N)
    for i in range(N):
        histo[i] = np.sum(H_flat==i);
    return histo

然后把它画到一个环上

canvas_h = 600 #画布高度
canvas_w = 600 #画布宽度
yc = int(canvas_h/2) #圆心位置y
xc = int(canvas_w/2) #圆心位置x
circle_r = 250 #半径
def draw_polar_histogram(histo):
    N = 360
    histo = histo.astype(float)
    histo /= np.max(histo)
    histo *= circle_r
    canvas = np.zeros((canvas_h, canvas_w, 3))  #画布
    cv2.circle(canvas, (yc, xc), circle_r, (255,255,255), -1) #画圆
    for i in range(N):
        theta = -i * np.pi / 180 #各个hue的弧度
        count = histo[i] #各个hue的数目
        #当前hue的柱子
        y1 = yc - int(circle_r * np.sin(theta))
        x1 = xc + int(circle_r * np.cos(theta))
        y2 = yc - int((circle_r-histo[i]) * np.sin(theta))
        x2 = xc + int((circle_r-histo[i]) * np.cos(theta))

        color_HSV = np.zeros((1,1,3), dtype=np.uint8)
        color_HSV[0,0,:] = [int(i/2),255,255] #每个角度的H
        color_BGR = cv2.cvtColor(color_HSV, cv2.COLOR_HSV2BGR) #将HSV转换为BGR
        B = int(color_BGR[0,0,0])
        G = int(color_BGR[0,0,1])
        R = int(color_BGR[0,0,2])
        cv2.line(canvas, (x1,y1), (x2,y2), (B,G,R), 3,cv2.LINE_AA) #画柱子
    canvas = cv2.circle(canvas, (yc, xc), 5, (0,0,0), -1) #圆心
    canvas = np.array(canvas,np.uint8)
    canvas = cv2.cvtColor(canvas,cv2.COLOR_BGR2RGB)
    return canvas

测试代码

img_rgb = cv2.imread('peacock.png')
img_hsv = cv2.cvtColor(img_rgb,cv2.COLOR_BGR2HSV)
h_hist = count_hue_histogram(img_hsv)
img_hue = img_hsv[...,0].copy()*2.0
img_sat = img_hsv[...,1].copy()/2.55
hist_img = draw_polar_histogram(h_hist)
plt.imsave('./temp/hist_img.png',hist_img)

图像评分

正式进入论文第3章节，根据图像的色相与每个模板扇形区域的边界的最小距离，以及每个像素的饱和度，来计算得分，公式为
$F(x,(m,\alpha))=\sum_{p\in X}||H(p)-E_{T_m(\alpha)}(p)||\cdot S(p)$F(x,(m,α))=p∈X∑​∣∣H(p)−ETm​(α)​(p)∣∣⋅S(p)
X代表图像，p代表图像每个像素位置的色度值(Hue)，$E_{T_m(\alpha)}(p)$ETm​(α)​(p)代表色度值p与扇形区域最接近的边界，$S(p)$S(p)代表HSV中的饱和度S

接下来先看某个模板中，$0°\sim360°$0°∼360°每个色相与边界的最短距离(角度值差)，如果色相在模板的扇形区域内，距离为$0$0：

#圆弧距：在hue圆上的角度差
def deg_distance(a, b):
    d1 = np.abs(a - b)
    d2 = np.abs(360-d1)
    d = np.minimum(d1, d2)
    return d
#是否在区域内
def in_border(h,center,width):
    return deg_distance(h,center)<width/2
#区域外的h，计算最近边界的圆弧距
def dist_to_border(h,border):
    H1=deg_distance(h,border[0])   
    H2=deg_distance(h,border[1])   
    H_dist2bdr = np.minimum(H1,H2)
    return H_dist2bdr
#计算当前模板的每个hue值的加权
def hue_weight(temp_name,alpha):
    hweight = []
    h = np.arange(360)
    for t in HueTemplates[temp_name]:
        center = t[0]*360 + alpha #中心位置
        width = t[1]*360 #宽度
        border = [center - width/2,center + width/2]  #起止位置
        temp_dist = dist_to_border(h,border) #色相与当前边界的距离
        temp_dist[in_border(h,center,width)]=0        
        hweight.append(temp_dist)
    hweight = np.array(hweight)
    hweight = hweight.min(axis=0)
    return hweight

接下来计算第一项：
$||H(p)-E_{T_m(\alpha)}(p)||$∣∣H(p)−ETm​(α)​(p)∣∣
直接将图像H值换成对应的最近边界距离即可：

# 利用template的hue权重对图像的hue直方图进行加权
def h_dist(img_h,temp_name,alpha):
    score = np.zeros((img_h.shape[0],img_h.shape[1]),dtype=np.float32)
    hw = hue_weight(temp_name,alpha)
    for i in range(360):
        score[img_h==i]= hw[i]
    return score

再跟饱和度一起，把图像与每个模板对应的得分都计算出来

#计算最小的模板得分
def cal_scores(img_hue,img_s):
    scores = np.zeros((len(HueTemplates.keys()),360),dtype=np.float32)
    temp_keys = list(HueTemplates.keys())
    for i in range(len(temp_keys)):#遍历模板
        print('temp_keys',temp_keys[i])
        for alpha in range(360):#每个旋转角度
            scores[i,alpha]=np.sum(np.multiply(h_dist(img_hue,temp_keys[i],alpha),img_s))
    return scores

提取最好的那个模板和角度，使得当前的图像得分最大

# 计算每个模板的得分
hue_scores = cal_scores(img_hue,img_sat)
#得到最好的template和alpha
[best_m,betst_alpha] = np.unravel_index(np.argmin(hue_scores),hue_scores.shape)
'''
temp_keys i
temp_keys V
temp_keys L
temp_keys mirror_L
temp_keys I
temp_keys T
temp_keys Y
temp_keys X
'''

把最好的模板、旋转角度与对应图像色相直方图放一起显示一波瞅瞅

#画图
hist_img = draw_polar_histogram(h_hist)
temp_keys = list(HueTemplates.keys())
print(temp_keys[best_m],betst_alpha)
temp_img = show_temp(temp_keys[best_m],betst_alpha)
overlay_img = cv2.addWeighted(temp_img,0.2,hist_img,1-0.2,0)
plt.figure(figsize=(8,8))
plt.imshow(overlay_img)
plt.imsave('./temp/overlay_img.png',overlay_img)
'''
T 46
'''

图像分割

上面有些模板有两个对称扇形，这两个扇形代表互补色，也就是颜色差距很大，但是当我们把扇形外面的像素压到扇形里面去的时候，中间的一些颜色到两个扇形的距离差不多，这时候就会出现相似的颜色被压到了完全不同的颜色。如下图所示

看(d)图，两个箭头从同一个位置出发，但是这个位置的色相到两个扇形的距离差不多，这时候如果随机分配颜色，就会出现孔雀脖子上相似的颜色被压到完全不同的两个颜色里面去了。

这时候，我们就需要对每个像素指定他所压缩的扇形，并且要保证图像比较连续的区域的颜色不要被压的不连续了，那么就可以先对图像每个像素进行二分类，当有两个扇形的时候，像素标签为$0或1$0或1 ，分别表示应该被划分到哪个区域。

文章的算法貌似比较耗时，大佬使用了opencv里面的一个超像素块划分方法，参考文档看前言。这个超像素块划分有点类似于目标分割，达到的效果如下：

原图戳这里

可以看出来，连续区域基本被划分到一块了。每一块都被标记为一个单独的类别。

我们现在要把这一对类别合并起来，变成两个类别。

分割流程大致可以分为这样

先不看颜色的连续性，单纯计算每个像素的色相到哪个扇形区域最近，先为每个像素单独打标签

#先计算每个像素到哪个扇形区域距离最短
def class_pix(h,temp_name,alpha):
    # 先计算当前模板下，每个hue对应哪个扇区
    img_label=np.zeros((h.shape[0],h.shape[1]),dtype=np.int32)
    hlabels = []
    hue_stand=np.arange(360)
    #扇区里外一起处理，某个扇区内的距离自己的扇区边界更近
    for t in HueTemplates[temp_name]:
        center = t[0]*360 + alpha #中心位置
        width = t[1]*360 #宽度
        border = [center - width/2,center + width/2]  #起止位置
        temp_dist = dist_to_border(hue_stand,border) #hue与当前扇区边界的距离，不要排除扇区内的hue的计算
        hlabels.append(temp_dist)
    hlabels = np.array(hlabels)
    hlabel = hlabels.argmin(axis=0)
    #再把图像的hue转换为对应扇区标签
    for i in range(360):
        img_label[h==i]=hlabel[i]
    return img_label

然后再统计opencv划分的每个超像素块里面的像素标签哪个多，就把当前超像素的标签置为谁：

#制作每个像素的标签，确定每个像素往哪个扇区shift
def split_img(img_hsv,img_h,temp_name,alpha):
    num_superpixels=200
    h_cls = class_pix(img_hue,temp_name,alpha)
    SEEDS = cv2.ximgproc.createSuperpixelSEEDS(img_hsv.shape[1],img_hsv.shape[0],img_hsv.shape[2],num_superpixels,10)
    SEEDS.iterate(img_hsv,4)
    V=np.zeros((img_hue.shape))
    N=V.shape[0]
    grid_num = SEEDS.getNumberOfSuperpixels() #超像素块个数
    labels = SEEDS.getLabels() #原图每个像素的属于哪个超像素块(标签)
    for i in range(grid_num):
        P=[[],[]]
        s=np.average(h_cls[labels==i]) #当前超像素块的所有hue的标签的均值
        if(s>0.5):#像素块大部分朝哪个扇区，就确定是哪个扇区
            s=1
        else:
            s=0
        h_cls[labels==i]=s
    return h_cls

这两个函数就保证了，颜色又连续，每个颜色被压缩的距离又最小，色差变化不会太大。

压色相

每个像素往哪个扇形块压，我们都知道了，那就执行最后一步，把像素按照下式压到对应扇形区域去
$H'(p)=C(p)+\frac{w}{2}(1-G_\sigma(||H(p)-C(p)||))$H′(p)=C(p)+2w​(1−Gσ​(∣∣H(p)−C(p)∣∣))
上式中C(p)代表扇形区域的中心，$w$w代表扇形区域的宽度，$G_\sigma$Gσ​代表标准高斯函数，均值为0，方差为$\sigma$σ，作者建议这个方差$\sigma$σ 等于$\frac{w}{2}$2w​

那么先把高斯函数写出来呗：
$\frac{1}{\sqrt{(2\pi)}\sigma}e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}$(2π)​σ1​e−21​(σx−μ​)2

#标准正态分布采样
def normalized_gaussian(X,mu,sigma):
    X=np.asarray(X).astype(np.float64)
    sigma=np.asarray(sigma).astype(np.float64)
    mu=X-mu
    M2=np.multiply(mu,mu)
    S2=np.multiply(sigma,sigma)
    return (1/(sigma*np.sqrt(2*np.pi)))*np.exp(-M2/(2*S2))

最后一步，根据$H'(p)$H′(p)的公式，压色相值

# 对每个标签块分别shift到对应扇区
def hue_shift(img_hsv,temp_name,alpha):
    img_hue = img_hsv[...,0].copy()*2.0
    new_img_hue = np.zeros_like(img_hue)
    h_cls = split_img(img_hsv,img_hue,temp_name,alpha) #确定每个像素对应的扇区

    i=0
    for t in HueTemplates[temp_name]: #每个扇区分别处理
        center = t[0]*360 + alpha
        width = t[1]*360
        mask = (h_cls==i) #当前扇区所对应的所有标签
        temp_hue = img_hue*mask #排除其他扇区超像素块的干扰
        center_dist = deg_distance(temp_hue,center)
        G = normalized_gaussian(center_dist,0,width/2)
        new_h = center + (1-G)*width/2
        new_img_hue=new_img_hue+new_h*mask
        i=i+1
    new_img_hsv = img_hsv.copy()
    new_img_hsv[...,0]=new_img_hue
    result_img = cv2.cvtColor(new_img_hsv,cv2.COLOR_HSV2RGB)
    return result_img

搞定！！

可视化看看结果

temp_name=temp_keys[best_m]
alpha=betst_alpha
plt.imshow(hue_shift(img_hsv,temp_name,alpha))

嗯？真心绿。

其实还有更好玩的，自己随便指定一个template和alpha，能玩一天，比如

temp_name='X'#temp_keys[best_m]
alpha=betst_alpha
plt.imshow(hue_shift(img_hsv,temp_name,alpha))
plt.imsave('./temp/result1.png',hue_shift(img_hsv,temp_name,alpha))

艾玛~~好玩

问题

最后一步的实现应该有问题，和大佬的代码有点出入，可能我对论文的最后一步理解有遗漏，后面再来查缺补漏吧。有大神看论文以后，也可以在这个基础上完善一下。

博客代码下载：

链接：https://pan.baidu.com/s/1DavXVylYx3Sd0K-LygnBoQ
提取码：ub0n