铅笔特效算法

铅笔特效方面的算法很多，像photoshop这种牛逼软件具有这个功能就不足为奇了。偶尔闲来无事，收集一些这方面的资料，觉得蛮有意思的。最近看到一篇香港中文大学贾佳亚教授的paper，当然他只是挂个名，因为他不是第一作者，应该是他学生的作品。废话不多说，这篇paper的题目是“Combining Sketch and Tone for Pencil Drawing Production”。看到这篇paper的实现效果，简直让人惊呆了。这是我见过的最好的铅笔特效算法，没有之一。文章的原理也很简单。希望感兴趣的朋友去下来研究研究。你值得拥有！

下面贴出自己的代码，甚为粗糙，还望看官莫要见笑。

 1 clear all;
 2 close all;
 3 clc;
 4 ks=10;
 5 dirNum=8;
 6 [FileName,PathName]=uigetfile({'*.jpg';'*.png';'*.bmp'},'Open an Image File');
 7 Img=imread([PathName FileName]);
 8 figure,imshow(Img);
 9 [H, W, sc] = size(Img);
10     if sc == 3
11         I = rgb2gray(Img);
12     end
13 PC = phasecongmono(I, 5, 3, 2.5, 0.55, 2.0);
14 figure,imshow(PC);
15 PC=PC.^0.65;
16 
17 %% convolution kernel with horizontal direction 
18     kerRef = zeros(ks*2+1);
19     kerRef(ks+1,:) = 1;
20 
21 %% classification 
22     response = zeros(H,W,dirNum);
23     Sresponse = zeros(H,W,dirNum);
24     for ii = 0 : (dirNum-1)
25         ker = imrotate(kerRef, ii*180/dirNum, 'bilinear', 'crop');
26         response(:,:,ii+1) = conv2(PC, ker, 'same');
27     end
28     [~ , index] = sort(response, 3); 
29     Dirs = (index(:,:,end)-1)*180/dirNum;
30     maxindex=index(:,:,end);
31     C = zeros(H,W,dirNum);
32     for k=1:dirNum
33     for i=1:H
34         for j=1:W
35             if(maxindex(i,j)==k)
36                 C(i,j,k)=PC(i,j);
37                     
38             end
39         end
40     end
41     end
42 Cresponse = zeros(H,W,dirNum);
43 for jj = 0 : (dirNum-1)
44         ker = imrotate(kerRef, jj*180/dirNum, 'bilinear', 'crop');
45 
46         Cresponse(:,:,jj+1) = conv2(C(:,:,jj+1), ker, 'same');
47 end
48 S=0.079*Cresponse(:,:,1)+0.079*Cresponse(:,:,2)+0.079*Cresponse(:,:,3)+0.079*Cresponse(:,:,4)+0.079*Cresponse(:,:,5)+0.079*Cresponse(:,:,6)+0.079*Cresponse(:,:,7)+0.079*Cresponse(:,:,8);
49 figure,imshow(1-S);
50 a=0:255;
51 p1 = 1/9*exp(-(255-a)/9);
52 p3=1/sqrt(2*pi*11) * exp(-(a-90).*(a-90) / (2.0*11*11));
53 p2= 0:255;
54 p2(:)=1/(225-105);
55 p2(1:55)=0;
56 p2(225:256)=0;
57 p = 0.76 * p1 + 0.22 * p2 + 0.02 * p3;
58 figure,plot(a,p);
59 hgram=round(p*H*W);
60 
61 Inew=histeq(uint8(I),hgram);
62 H=fspecial('gaussian',13,2);
63 Inew = imfilter(Inew,H,'replicate');
64 Inew=double(Inew);
65 Sn=1-mapminmax(S);
66 out=Inew.*Sn*0.55;
67 figure,imshow(uint8(out));

需要提醒的是，本代码未能实现texture rendering的步骤，所以效果与文章还是有一定差距的。