前言

最近做換臉，在膚色調整的那一塊，看到一個有意思的文章，復現一波玩玩。不過最后一步掉鏈子了，有興趣的可以一起討論把鏈子補上。

主要是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