基於《Combining Sketch and Tone for Pencil Drawing Production》的影象鉛筆畫演算法的實現

阿新 • • 發佈：2019-01-06

一，借鑑：

本文借鑑了CSDN博主風吹夏天對此論文演算法的理解：風吹夏天的影象鉛筆畫演算法，以及香港中文大學Cewu Lu等人寫的該論文的主頁。原文作者和博主風吹夏天都給過程式碼，但是程式碼不全。我仔細看了原論文和該博主的文章後，基本上，大致演算法思想就理通了。本文對於具體演算法細節就不細述了，個人建議還是先進行原論文的研讀，對該論文所進行的步驟先有個大致的瞭解。

二，演算法思路

1，首先需要產生筆畫結構

大致思路：對原影象進行梯度運算，得出大致輪廓 -----> 設計8個方向的卷積核，並依據論文中的公式對畫素依據最大值所在子影象進行方向分類 -----> 再依據論文中的公式對每個方向的畫素分別和對應的卷積核進行卷積

幾個注意的地方：

對於卷積核的大小，論文中寫的是原圖片高和寬的三十分之一，論文作者設計卷積核的目的是希望產生鉛筆畫中一筆一筆勾畫出的痕跡，但如果太大，會導致輪廓圖不清晰，卷積核大小應該設定在3~13間效果好些；

根據輸入圖片的差異，原圖片可能需要去噪，若最後輪廓不清晰，可能需要直方圖均衡化使輪廓相對清晰。

2，其次產生色調圖

論文作者根據對大量鉛筆圖片的研究，發現它們的直方圖如下，而我們要做的就是將原圖片的色階的直方圖轉換為以下的直方圖，而以下的直方圖可以由高斯分佈函式，均勻分佈函式和拉布拉斯分佈函式帶權值相加後近似得到(三個權值相加為1)，它們的權值比為w1:w2:w3，此後，進行SML對映將原圖片的直方圖轉換為以下的直方圖就可以了。

注意的地方：一般採用w1:w2:w3 = 2:22:76，因為鉛筆畫較談，w3越大，圖片越亮，但對於有些圖片，如果太亮，有些較暗的細節就會消失，此時可以選擇w1:w2:w3 = 11:37:52，此外也可以根據最後效果調整三個引數大小，一般，w3介於52~76間，w1最好在10以下。

3，由鉛筆畫背景和色調圖生成色調渲染圖

我使用的是如下的鉛筆背景圖：

這一步主要是為了解下面這個方程組中的β矩陣，其中H(x)為上面這張鉛筆背景圖，J(x)為第二步生成的色調圖，現在主要求解β(x)，最後色調渲染圖T(x) = H(x).^β(x)，這樣做的目的是使色調圖有鉛筆反覆抹畫的痕跡，且保持直方圖與鉛筆畫的直方圖一致：

下面是求解公式，風吹夏天博主的文章有詳細的求解過程，大致思路是先把影象矩陣變成向量，再進行求解：

注意的地方：根據最後效果的考慮，需要對背景鉛筆圖片和生成的色調渲染圖片進行gamma矯正。

4，最後將第一步生成的結構圖片與第二三步生成的色調渲染圖片直接點乘就可以了

注意的地方：由於實現過程中，有許多地方用到矩陣點乘和矩陣相乘，因此在計算時應該先用im2double函式將影象矩陣中每個畫素對映到0~1間進行計算。

5，輸入/得到彩色影象 VS 輸入/得到灰度影象

輸入灰度影象/彩色影象：輸入時進行判斷就維數判斷就可以了，但在呼叫函式grayPencilDrawing時，輸入的一定要是灰度影象；

輸出灰度影象/彩色影象：輸出灰度鉛筆畫影象，呼叫函式grayPencilDrawing，函式返回時，返回的就是灰度鉛筆影象；輸出彩色鉛筆畫影象，這時輸入的肯定是彩色影象，需要先將原彩色影象由rgb空間轉成yuv空間，對於yuv格式，其中y分量是色階，uv分量是色調，只需要對y分量呼叫函式grayPencilDrawing，最後再由yuv空間轉成rgb空間。以下是yuv和rgb轉換公式：

注意的地方：需要對輸入進行iif else條件判斷。

三，具體程式碼

具體需要改的引數在第一個程式碼段中都有標識，可以先讀懂程式碼，然後修改相關引數，此外，對於每一張圖片，所需要的引數有可能不同，主要依據最後效果修改引數：

程式碼1：輸入圖片和背景圖片：該程式碼可能作為一種介面程式碼，將論文中需要修改的引數都提取了出來，只需修改4~17行的相關引數，直接函式呼叫就可以檢視結果了。此外，對最後rgb空間和yuv空間進行了考慮：

clc;
clear;
%==============================引數修改====================================
Input = imread('input1.jpg');%輸入原始圖片
dirNum = 8;%表示具有dirNum個方向的卷積核
[m, n, src] = size(Input);
ks = floor(min(m/50, n/50) + 0.5);%論文中要求為高或寬的1/30，實際上不能這麼大
ks = 10;%卷積核大小，建築類圖片為10~13，風景花卉類圖片為3~6
strokeDepth = 2;%對生成的輪廓圖多描幾遍

w1 = 0.57;%高斯分佈權重
w2 = 0.37;%均勻分佈權重
w3 = 0.06;%拉布拉斯分佈權重

Back = rgb2gray(imread('pencil.jpg'));%輸入背景鉛筆圖
backDepth = 0.7;%對鉛筆背景圖pencil.jpg多描幾遍
renderDepth = 0.7;%對生成的紋理渲染圖多描幾遍

%如果輸入的是灰度影象，對其進行處理
if src == 2
    gray_out = rgb2gray(Input);
    gray_out = grayPencilDrawing(gray_out, dirNum, ks, strokeDepth, w1, w2, w3, Back, backDepth, renderDepth);
    rgb_out = Input;
end
%如果輸入的是彩色影象，對其進行處理
if src == 3
    rgb_out = mat2gray(Input);
    R = rgb_out(:, :, 1);
    G = rgb_out(:, :, 2);
    B = rgb_out(:, :, 3);
    Y = 0.299*R + 0.587*G + 0.114*B;
    U = -0.147*R- 0.289*G + 0.436*B;
    V = 0.615*R - 0.515*G - 0.100*B;
    Y = uint8(Y .* 255);
    Y = grayPencilDrawing(Y, dirNum, ks, strokeDepth, w1, w2, w3, Back, backDepth, renderDepth);
    %yuv_out = cat(3, Y, U, V);%按照第三維組合
    rgb_out(:,:,1) = Y + 1.14 * V;
    rgb_out(:,:,2) = Y - 0.39 * U - 0.58 * V;
    rgb_out(:,:,3) = Y + 2.03 * U;
    figure, imshow(rgb_out), title('彩色鉛筆圖');
    imwrite(rgb_out, 'coloredPencil.jpg');
end

程式碼2：具體呼叫函式的實現，這一步也是上述演算法思路前四個的實現

function Out = grayPencilDrawing(Input, dirNum, ks, strokeDepth, w1, w2, w3, Back, backDepth, renderDepth)
%===============================1：獲取輪廓圖===============================
%讀入單通道的圖片
    [h, w] = size(Input);
    I_gray = Input;
    I_gray1 = I_gray;
    %Input = histeq(Input, 256); %針對有些圖片層次不清晰，若最後輪廓圖不清晰，加上這兩句
    %figure, imshow(Input);
    I_gray = im2double(I_gray);
    imX = [abs(I_gray(:,1:(end-1)) - I_gray(:,2:end)), zeros(h, 1)];%最右邊增加一列0
    imY = [abs(I_gray(1:(end-1),:) - I_gray(2:end,:)); zeros(1, w)];%最下邊增加一行0
    imX = imX .^ 2;
    imY = imY .^ 2;
    imEdge = sqrt(imX + imY);%論文中邊緣公式(1)實現
    imEdge = immultiply(imEdge, 5);    
%%%%%figure, imshow(imEdge, []), title('邊緣提取');
    imwrite(imEdge, 'stroke_tmp1.jpg');
%水平方向為1的卷積核
    kerRef = zeros(ks*2+1);
    kerRef(ks+1,:) = 1;
%對卷積核依次旋轉180/dirNum度並使用論文中的公式先計算8個方向卷積，計算出邊緣
    response = zeros(h, w, dirNum);
    for n = 1: dirNum
        ker = imrotate(kerRef, (n-1)*180/dirNum, 'bilinear', 'crop');
        response(:, :, n) = conv2(imEdge, ker, 'same');%對梯度圖進行八個方向卷積核，論文中邊緣公式(2)實現
    end
    [~ , index] = max(response, [], 3); %index為一個m*n矩陣，每個元素為dirNum個矩形中對應元素的最大縱座標號(範圍為1~dirNum)，最大值即為該點方向
    C = zeros(h, w, dirNum);
    for n=1:dirNum
        C(:,:,n) = imEdge .* (index == n); %index == 2也是一個矩陣，元素值為0或1。之後與imEdge對應元素相乘，論文中邊緣公式(3)實現
    end
    Spn = zeros(h, w, dirNum);
    for n = 1: dirNum
        ker = imrotate(kerRef, (n-1)*180/dirNum, 'bilinear', 'crop');
        Spn(:, :, n) = conv2(C(:,:,n), ker, 'same'); %對得到的各個方面的響應再次進行方向卷積，論文中邊緣公式(4)實現
    end
    Sp = sum(Spn, 3); %按照第三維對這dirNum個矩陣的元素各自累加
    Sp = (Sp - min(Sp(:))) / (max(Sp(:)) - min(Sp(:)));%min(sp(:))和max(sp(:))分別表示矩陣sp中元素的最小值和最大值，即灰度拉伸
    Stroke = 1 - Sp;%由於梯度圖背景黑線白，現在使其背景白線黑
    Stroke = Stroke .^ strokeDepth;%輪廓圖描深
    figure, imshow(Stroke), title('Stroke圖');
    imwrite(Stroke, 'stroke.jpg');
    
%===============================2：獲取色調圖===============================    
%直方圖匹配
    p = zeros(1, 256);
    v = zeros(1, 256);
    i = 0: 255;
    p1 = (1/9) * exp(-(255 - i)/9);
    %axis([0 255 0 1]);
    p2 = zeros(1, 256);
    p2(105: 225) = 1/(225-105);
    p3 = (1/sqrt(2*pi*11)) * exp(-(i-90).*(i-90)/(2.0*11*11));
    p = p1*w1 + p2*w2 + p3*w3;
    t = p;
    %figure, plot(imhist(I_gray)/(h*w)), title('原始直方圖');
    for i = 1: h
        for j = 1: w
            v(I_gray1(i, j)+1) = v(I_gray1(i, j)+1) + 1;
        end
    end
    v = v / (h*w);
    for i = 2: 256
        p(i) = p(i-1) + p(i);
        v(i) = v(i-1) + v(i);
    end
    p = p /p(256);
%計算源影象到目標影象累積直方圖中各灰度級的差的絕對值
    scrMin = zeros(256, 256);
    for x = 1: 256
        for y = 1: 256
            scrMin(x, y) = abs(v(x) - p(y));
        end
    end
    for x = 1: 256
        minX = 1;
        minValue = scrMin(x, 1);
        for y = 2: 256
            if(minValue > scrMin(x, y))
                minValue = scrMin(x, y);
                minX = y;
            end
        end
        HistogramMapping(x) = minX;
    end
%對原圖根據HistogramMapping進行灰度級對映
    I_tone = I_gray1;
    for x = 1: h
        for y = 1: w
            I_tone(x, y) = HistogramMapping(I_gray1(x, y) + 1) - 1;
        end
    end
%%%%%figure, subplot(211), plot(t), title('目標直方圖');
%%%%%subplot(212), plot(imhist(I_tone)/(h*w)), title('色調對映圖的直方圖');%因為進行對映時，目標影象中有些灰度值不存在，因此曲線中接近0的曲折部分表示某些灰度值不存在
    figure;
    subplot(121), imshow(I_gray1), title('原始灰度圖');
    subplot(122), imshow(I_tone), title('色調對映圖');
    imwrite(I_tone, 'tone.jpg');    

%===============================3：紋理渲染圖=============================== 
    %subplot(251), imshow(I);%執行這9行可以檢視背景鉛筆圖效果
    a = 1;
    for i = 0.2: 0.2: 1.8
        %a = floor(i / 0.2 + 0.5) + 1;
        a = a + 1;
        I0 = ((double(Back)/255).^i);
        %subplot(2, 5, a), imshow(I0);
    end
    I0 = ((double(Back)/255) .^ backDepth);
    %figure, imshow(I0);
    I0 = uint8(I0 * 255);
    imwrite(I0, 'New_texture.jpg');
%開始計算P .^ β = J中的β矩陣
    theta = 0.2;
    P = im2double(imread('New_texture.jpg'));
    J = im2double(imread('tone.jpg'));
%初始化為向量方便計算
    P = imresize(P, [h, w]);
    P = reshape(P, h*w, 1);
    logP = log(P);
    logP = spdiags(logP, 0, h*w, h*w);
    
    J = imresize(J, [h, w]);
    J = reshape(J, h*w, 1);
    logJ = log(J);
%兩個梯度
    e = ones(h*w, 1);
    Dx = spdiags([-e, e], [0, h], h*w, h*w);
    Dy = spdiags([-e, e], [0, 1], h*w, h*w);
%帶入求解公式計算
    A = theta * (Dx * Dx' + Dy * Dy') + (logP)' * logP;
    b = (logP)' * logJ;
%計算向量形式的beta
    beta = pcg(A, b, 1e-6, 60);
%轉化成矩陣形式並拉伸
    beta = reshape(beta, h, w);
    beta = (beta - min(beta(:))) / (max(beta(:)) - min(beta(:)))*5;
    
    P = reshape(P, h, w);
    T = P .^ beta;
    A = T .^ renderDepth;
%%%%%imshow(A), title('色調渲染圖');;
    imwrite(A, 'new_tone.jpg');
    
%===============================3：結構圖與色調圖整合===============================
    I = im2double(imread('stroke.jpg'));
    S = A .* I;
    figure, imshow(S), title('灰色鉛筆畫');
    imwrite(S, 'grayPencil.jpg');
    Out = S;
end

四，一些實現效果

原圖灰度鉛筆畫彩色鉛筆畫