赛璐璐风格图像线稿提取，算法优化中
一、导读
在日系插画领域，赛璐璐风格不会过多地塑造出复杂的光影和小调子。以线条造型为主，清晰的光影表达，以及明快的配色是这个风格的画师普遍追求的目标。基于这个前提，赛璐璐风格就是做自动线稿提取算法的很好切入点。能否详尽地提取出原始图像线稿中所有的线条，并且结果能让人在视觉基础上感到满意，应该是这个算法应该追求的主要目标。与此同时，这样的线稿提取算法也能成为未来自动草稿细化的基础。为了更好地实现机器辅助设计，我们应该深刻地从图像处理的立场上研究这种风格。
二、基本算法构架的解释以及测试集的选择
算法主要分为三个部分：
1.首先对图像进行降噪，并且平滑整个图像，抹去图像中的噪音和一些小调子，以更好地进行边缘提取.其中，我选择的降噪算法是经典的ROF降噪模型：其优点在于经历多次迭代后仍然可以保持图像的边缘，并且同时可以删去噪音和小调子。相关知识请参考如下网站：
https://blog.csdn.net/cyh706510441/article/details/45194223
2.使用Canny边缘提取算法，提取出图像中的主要细节。
有关Canny边缘提取算法的前置知识请参考：
https://blog.csdn.net/Jamence/article/details/84645554
3.结果的后续美观性上的处理。
测试集的选择为：1.画师Anmi；2.画师Rei；3.画师Kirarin；4.动画《JOJO的奇妙冒险~黄金之风》
步骤1对于一些有小调子并且存在噪音的图，例如插画师rei以及Anmi的插画，就很有必要进行这一步，其可以防止Canny算法过渡地提取边缘。另外值得注意的是：不论是JPG格式还是PNG格式的图像，其在MATLAB中都是以红绿蓝三个通道（也就是一个三维数组），储存了三个不同的灰度值，最后叠加形成的彩色图。我们提取边缘的时候，应该对三个通道都进行提取，最后叠加成图。在图像降噪的环节中，我们的ROF降噪模型的迭代步数都是100，下面是原图以及实验结果：
可以看见，整体的形体被提取了出来，而且也没有损失细节。后续算法的改进会更加注重线条的优化，以及噪点的删去。同时，我们可以看到对红绿蓝三层通道进行处理后的的线稿提取结果，还可以加强外轮廓（这在线稿细化中是非常重要的）。虽然仍然会有小调子和噪音的出现。但不论如何，Canny边缘提取算法和ROF降噪模型的完美配合，显示出了广阔的应用前景。
对于步骤二，为何使用Canny边缘提取算法，是因为例如传统的Sobel梯度边缘提取算法，提取出的是图像梯度值较高的“强边缘”，但会忽视图像中的弱边缘。下图中咦Kirarin老师的图为例，左起至右依次是原图—Sobel边缘提取算法—Canny算子的结果图：
对于线稿提取，我们自然是不想遗漏任何的线条信息的。但Canny算子此时也会显出它的弊端：我们也不是毫无保留地提取所有线条。如何更好地改进Canny算法，使得其可以提取出赛璐璐图片的“主要信息”而忽略其“次要信息”呢？首先我们需要定义何种线条信息在图片中是“主要的”，何种是“次要的”。紧接着我们可以考虑设置一个阈值——即当某一像素点的梯度大于一定取值的时候，这个边缘信息才被保留下来，否则便舍弃。或许在未来，这一领域会和神经网络紧密联系起来，而且值得一提的是，日本的很多学者也已经开始用神经网络解决此类问题了。
对于步骤三，对提取结果的进一步优化，使其符合我们的视觉审美，这一步我目前仍然在优化和进行中：目前我的思路如下：一个好的线稿无非要满足一下两个特整：1.外轮廓粗而内轮廓细（例如上图的向日葵的花瓣外轮廓很粗很实在，而内轮廓很细或者很淡。）2.对一些闭塞进行加深（通俗一点讲就是两条线相交的位置，临近交点的线条都要加深） 而除此之外，便是要靠艺术家自身的设计功力去塑造漂亮的线稿了。请记住我们的目的永远不是用机器取代人来设计，而是使用机器解决艺术创作中的重复性工作，解放创造力。
三、代码的公布以及使用范例
首先是ROF降噪算法的代码（注意是函数，本体要调用这个代码）
function J = tv(I,iter,dt,ep,lam,I0,C)
%% Total Variation denoising.
%% Example: J=tv(I,iter,dt,ep,lam,I0)
%% Input: I - 样本图像 (二进制浮点数 array 灰度图 level 1-256),
% iter - 迭代次数
% dt - 时间步长 [0.2],
% ep - epsilon (正规化因子) [1],
% lam - 松弛因子 [0],
% I0 - 无噪声图像 [I0=I]
% ([]内为默认值)
% Output: 像素矩阵
%% 判断是否输入相关参数
if ~exist(\'ep\')
ep=1;
end
if ~exist(\'dt\')
dt=ep/5; % dt需满足CFL条件
end
if ~exist(\'lam\')
lam=0;
end
if ~exist(\'I0\')
I0=I;
end
if ~exist(\'C\')
C=0;
end
[ny,nx]=size(I); ep2=ep^2;
%% 循环迭代
for i=1:iter,
% 构建数值微分
I_x = (I(:,[2:nx nx])-I(:,[1 1:nx-1]))/2;
I_y = (I([2:ny ny],:)-I([1 1:ny-1],:))/2;
I_xx = I(:,[2:nx nx])+I(:,[1 1:nx-1])-2*I;
I_yy = I([2:ny ny],:)+I([1 1:ny-1],:)-2*I;
Dp = I([2:ny ny],[2:nx nx])+I([1 1:ny-1],[1 1:nx-1]);
Dm = I([1 1:ny-1],[2:nx nx])+I([2:ny ny],[1 1:nx-1]);
I_xy = (Dp-Dm)/4;
% 计算递推关系式
Num = I_xx.*(ep2+I_y.^2)-2*I_x.*I_y.*I_xy+I_yy.*(ep2+I_x.^2);
Den = (ep2+I_x.^2+I_y.^2).^(3/2);
I_t = Num./Den + lam.*(I0-I+C);
I=I+dt*I_t; %% 更新图像
J= I;
end 然后是Canny边缘提取算法：
function m = canny1step( src, lowTh)
%canny函数第一步，求去x，y方向的偏导，模板如下：
% Gx
% 1 -1
% 1 -1
% Gy
% -1 -1
% 1 1
%------------------------------------
% 输入：
% src：图像，如果不是灰度图转成灰度图
% lowTh：低阈值
% 输出：
% m： 两个偏导的平方差，反映了边缘的强度
% theta：反映了边缘的方向
% sector：将方向分为3个区域，具体如下
% 2 1 0
% 3 X 3
% 0 1 2
% canny1：非极大值
% canny2：双阈值抑制
% bin ： 二值化
%---------------------------------------
[Ay Ax dim ] = size(src);
%转换为灰度图
if dim>1
src = rgb2gray(src);
end
src = double(src);
m = zeros(Ay, Ax);
theta = zeros(Ay, Ax);
sector = zeros(Ay, Ax);
canny1 = zeros(Ay, Ax);%非极大值抑制
canny2 = zeros(Ay, Ax);%双阈值检测和连接
bin = zeros(Ay, Ax);
for y = 1:(Ay-1)
for x = 1:(Ax-1)
gx = src(y, x) + src(y+1, x) - src(y, x+1) - src(y+1, x+1);
gy = -src(y, x) + src(y+1, x) - src(y, x+1) + src(y+1, x+1);
m(y,x) = (gx^2+gy^2)^0.5 ;
%--------------------------------
theta(y,x) = atand(gx/gy) ;
tem = theta(y,x);
%--------------------------------
if (tem<67.5)&&(tem>22.5)
sector(y,x) = 0;
elseif (tem<22.5)&&(tem>-22.5)
sector(y,x) = 3;
elseif (tem<-22.5)&&(tem>-67.5)
sector(y,x) = 2;
else
sector(y,x) = 1;
end
%--------------------------------
end
end
%-------------------------
%非极大值抑制
%------> x
% 2 1 0
% 3 X 3
%y 0 1 2
for y = 2:(Ay-1)
for x = 2:(Ax-1)
if 0 == sector(y,x) %右上 - 左下
if ( m(y,x)>m(y-1,x+1) )&&( m(y,x)>m(y+1,x-1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 1 == sector(y,x) %竖直方向
if ( m(y,x)>m(y-1,x) )&&( m(y,x)>m(y+1,x) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 2 == sector(y,x) %左上 - 右下
if ( m(y,x)>m(y-1,x-1) )&&( m(y,x)>m(y+1,x+1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
elseif 3 == sector(y,x) %横方向
if ( m(y,x)>m(y,x+1) )&&( m(y,x)>m(y,x-1) )
canny1(y,x) = m(y,x);
else
canny1(y,x) = 0;
end
end
end%end for x
end%end for y
%---------------------------------
%双阈值检测
ratio = 2;
for y = 2:(Ay-1)
for x = 2:(Ax-1)
if canny1(y,x)<lowTh %低阈值处理
canny2(y,x) = 0;
bin(y,x) = 0;
continue;
elseif canny1(y,x)>ratio*lowTh %高阈值处理
canny2(y,x) = canny1(y,x);
bin(y,x) = 1;
continue;
else %介于之间的看其8领域有没有高于高阈值的，有则可以为边缘
tem =[canny1(y-1,x-1), canny1(y-1,x), canny1(y-1,x+1);
canny1(y,x-1), canny1(y,x), canny1(y,x+1);
canny1(y+1,x-1), canny1(y+1,x), canny1(y+1,x+1)];
temMax = max(tem);
if temMax(1) > ratio*lowTh
canny2(y,x) = temMax(1);
bin(y,x) = 1;
continue;
else
canny2(y,x) = 0;
bin(y,x) = 0;
continue;
end
end
end%end for x
end%end for y
end%end of function最后是我们的主体代码
close all; clear all; clc; mI = double( imread( \'C:/Matlab/IMG_3228.png\' ) );%此处输入你的图像 sigma1 = 15; %高斯正态分布标准差 gausFilter = fspecial(\'gaussian\',[3 3],sigma1); %高斯滤波器 %设置高斯滤波是提取边缘后的进一步去除噪音点 mI1 = mI(:,:,1); mI2 = mI(:,:,1); mI3 = mI(:,:,1);%三个图层分别对应红绿蓝三个通道的灰度图 % mI1 = TVSample(mI( :, :, 1 ),100); % mI2 = TVSample(mI( :, :, 2 ),100); % mI3 = TVSample(mI( :, :, 3 ),100); [ny,nx]=size(mI1); A = canny(mI1,graythresh(mI1)); %将图像进行TV迭代， TV函数的格式为：TV（图像，迭代次数，（后续参数）） % A = edge(mI1);%为了让你方便和其他边缘提取算法比较，这里直接用MAT内置的edge来做，使用的时候去除 % 就行 A = abs(A);%取像素绝对值 A = imcomplement(A);%二值化 A = mat2gray(A);%转为标准灰度（0-255） A = imfilter(A,gausFilter,\'replicate\'); %高斯滤波除一下结果图的噪音 % A = imbinarize(A,\'global\');%像素取反，更清楚地看见边缘 B = canny(mI2,graythresh(mI1)); %将图像进行TV迭代， TV函数的格式为：TV（图像，迭代次数，（后续参数）） % B = edge(mI2); B = abs(B); B = imcomplement(B); B = mat2gray(B); B = imfilter(B,gausFilter,\'replicate\'); % B = imbinarize(B,\'global\'); C = canny(mI3,graythresh(mI1)); %将图像进行TV迭代， TV函数的格式为：TV（图像，迭代次数，（后续参数）） % C = edge(mI3); C = abs(C); C = imcomplement(C); C = mat2gray(C); C = imfilter(C,gausFilter,\'replicate\'); % C = imbinarize(C,\'global\'); D = A+B+C; % D = 1 * D.^5 % tresh = min(min(D)) % D(find(tresh+1<D<tresh+3))=D(find(tresh+1<D<tresh+3))-1; % [ny,nx]=size(A); % E = A % for i = 3:ny-2 % for j = 3:nx-2 % if A(i,j) == 0 % E(i+1,j) = 0; E(i-1,j) = 0;E(i,j+1) = 0;E(i,j-1) = 0; % E(i+1,j+1) = 0;E(i+1,j-1) = 0;E(i-1,j+1) = 0;E(i-1,j-1) = 0; % % end % end % end %这里的一段代码是为后续优化设计的，还没完善。 figure(1); imagesc( A ); colormap( gray ); axis off;%输出提取后1通道图像 figure(2); imagesc( B ); colormap( gray ); axis off;%输出提取后的2通道图像 figure(3); imagesc( C ); colormap( gray ); axis off;%输出提取后的3通道图像 figure(4); imagesc( D ); colormap( gray ); axis off; %输出提取后的复合图像
四、写在最后
我曾经一直在迷茫：其实我数学并不好，学什么都要比别人吃力很多，我喜欢画画，我为什么不能选自己喜欢的东西。但不管如何抱怨，命运还是把我带到了国内数学系很牛的一个大学。我从自暴自弃到现在，思考了很多东西：大概这一切的确冥冥之中有天意，我的很多画家朋友用他们自己的努力，去推动着商业美术的发展，那我也能用我自己的方式去实现一些东西，不管会绕多远的路，遇到多少挫折也好，就像老谢的lesson3教的：最快的捷径就是绕远路。
我很喜欢布加拉提说的：“我不后悔，虽然是身处于这样的世界，可我还是想走在“自己所坚信的道路”上，虽然现在我只能逃走，但只要找到弱点，就一定能打败老板，我一定会找出他的弱点！”即使滚石早就预言着他不论如何都会死，但他用他的觉悟诠释了什么叫做“知其不可而为之”。
ARRIVEDERC！