相关论文的链接：Combining Sketch and Tone for Pencil Drawing Production
第一次看《Combining Sketch and Tone for Pencil Drawing Production》一文是在两年前，随意看了一下，觉得论文里的公式比较多，以为实现有一定的难度，没有去细究，最近在作者主页上看到有 [code of direction classification] 部分代码，下载后觉得还是有自己实现的可能，下面记录下自己实现过程中的一些体会和心得。
铅笔画其实一直是一个比较难以获得较为理想效果的算法，我看到的论文里这篇文章应该是说相当优秀的。总的来说，其算法分为两个步骤：
1、Line Drawing with Strokes 得到一幅图 S。
2、Tone Mapping 得到另外一副图T。
3、得到最终结果 R = S * T;
应该说第一步决定了最终的效果，作者通过以下四个步骤得到S图。
（1）、对原图进行边缘检测，作者论文给出的公式是：
按照这个公式实现的效果实际上检测的效果很弱，我认为作者真正意义上可能不是使用的改公式，因为这一步对最终效果的影响很大， 我采用了一些其他能够更好的检测出效果的边缘检测算法，如果Sobel或者PS里FindEges之类的算法。
作者认为这个公式得到的结果含有太多的噪音并且边缘部分的线条很多不连续，因此提出继续一下几个步骤得到更稳定的效果。
（2) 对得到的G进行8个方向的卷积，卷积核为沿指定的方向为1，其他的值均为0（实际上考虑抗锯齿问题，用了双线性插值得到卷积核的），卷积核的大小论文提出为图像宽度或高度的1/30（这个我觉得有点不行，当太大时，会有明显的不合理线条出现)，具体公式如下：
在论文给出的相关代码中，有如下部分：
%% convolution kernel with horizontal direction
kerRef = zeros(ks*2+1);
kerRef(ks+1,:) = 1;
%% classification
response = zeros(H,W,dirNum);
for ii = 0 : (dirNum-1)
ker = imrotate(kerRef, ii*180/dirNum, 'bilinear', 'crop');
response(:,:,ii+1) = conv2(imEdge, ker, 'same');
end
其实这里的卷积就是按照指定的角度的运动模糊，我们可以用matlab的代码来进行验证：
比如当角度为22.5（ii = 1)，核的大小为5（ks = 5)时，按照上面的代码得到的ker变量为：
归一化后的结果为：
而对应的运动模糊的卷积矩阵为（对应的matlab代码为： H = fspecial('motion', 2*5+1, 22.5）：
可见只有很小的差别。
（3） 得到各个方向的卷积结果后，对每一个像素点，具有最大的卷积值的那一个方向的响应设置为G，而其他方向的响应设置为0，原文的话语是：
The classification is performed by selecting the maximum value among the responses in all directions and is written as
我觉得有点不可思议的是上面语句说的是最大值，下面给出的公式确是最小值，这难道是大家的笔误。（实际上应该是最大值的）。
还有注意的是这里的G(p)是指公式（1）中的G。
（4）对得到的各个方面的响应再次进行方向卷积，即：
论文中提出要对这个结果进行卷积并进行反相处理得到结果S，这个其实就看你自己的编码方式了。
这里给出以上4个步骤一些中间结果：
原图 边缘检测图 22.5度的卷积图 22.5度的响应图 中间结果图S
第二部的Tone Mapping 实际上也是由两个步骤来实现的。
（1）直方图匹配。
论文中并没有这样起小标题，而是用了较大的篇幅说了一堆事情，我总结一句话就是，根据对大量的手工绘制的铅笔画图像数据的观察和分析，其直方图的分布和我们拍摄的图像有很大的不同，但是都成一定的规律，这个规律可以用一定的经验公式来表达。因此，我们可以设定一个固定的直方图，然后将图像自身的直方图映射到这个直方图，作为结果。
简单的阐述下过程吧。借作者论文中的一副图像来做说明。在下图中，(a)是一副手工绘制的铅笔画，(b)图是硬性划分的阴影、中间调及高光图像，(c)图分别对应阴影、中间调及高光的直方图分布。可以看出，阴影部分基本成正态分布，中间调大致为均匀分布，而高光部分成拉普拉斯分布。因此，作者对三个部分分别构建了三个函数来模拟曲线。
a、高光部分。
在人手工绘制的铅笔画中，由于纸张一般都是白色，因此，高光占有的比例实际上肯定是非常大的，这在直方图中反应就是在接近色阶255时分布曲线越陡峭。作者提出以下函数作为这部分的分布曲线。
而关于中间调，则用一截水平分布的线条来模拟：
而暗调部分则表明了图像深度的变化，用一个高斯曲线来模拟：
以上公式对于不同的铅笔画来说，各部分的权重是不一样，因此作者提出了一个综合的公式来获取最终的铅笔画对应的直方图：
根据不同的需要，可以调节不同的权重系数来达到不同的混合效果。
作者根据经验，提出了一组数据：
如果权重都相同，三部分的曲线如下图所示：
可见暗调部分的比例和人工绘制的不协调，如果按照上表中论文给出的数据，得到的最终混合直方图效果如下图：
似乎也和人工的结果不一致，因此我认为此处也是论文的一个笔误或者错误，w1和W3的值很明显弄反了，即W1应该为52，W3为11，修改后的直方图为：
基本和人工绘制的一致了，同时注意到上述曲线有两个巨变之处，实际处理时需要对曲线进行一定程度的平滑最好。
附绘制曲线的matlab代码：
a=1:255;
p1 = 1/9*exp(-(255-a)/9);
p3=1/sqrt(2*pi*11) * exp(-(a-80).*(a-80) / (2.0*11*11));
p2= 1:255;
p2(:)=1/(225-105);
p2(1:105)=0;
p2(225:255)=0;
p = 0.52 * p1 + 0.37 * p2 + 0.11 * p3;
plot(a,p)
接下来的工作就是进行直方图匹配，这方面的资料可以参考何斌那本VC++的数字图像处理数，里面有SML和GML两种匹配方式。
（2）纹理渲染
这一部分论文说的也很简单，就是求解一个方程：
这不是我擅长的东西，有兴趣的朋友可能要自己研究下，我也没有实现它，由这一步得到中间结果T。
那最后一步就是将S 和 T相乘，就类似于PS中的正片叠底 混合算法。
以上的一些操作都是针对灰度图像，对于彩色图像，如果直接将三通道分开，然后分别调用灰度算法，再合成这样处理， 是有问题的，出来的结果很不理想，这主要是由于各通道在进行Line Drawing with Strokes时所获得的线条方向不太可能完全一致，导致合成后偏色，解决的方式有多种，比如将RGB转换到HSL空间，然后对L分量进行处理，然后在转换到RGB空间，或者借用LAB空间作为中转平台也是可以的。
下面贴出我的C++部分的核心代码（并不能直接运行，大概体现了算法的思路）:
extern IS_RET EdgeCoarse(TMatrix *Src, TMatrix *Dest); extern IS_RET MotionBlur(TMatrix *Src, TMatrix *Dest, int Length, float Angle, EdgeMode Edge); extern IS_RET SMLHistgramMaping(TMatrix *Src, TMatrix *Dest, int* HistgramB, int *HistgramG, int *HistgramR); extern IS_RET BlendImage(TMatrix *Base, TMatrix *Mixed, TMatrix *Result, BlendMode BlendOp, int Opacity); extern IS_RET GuassBlur(TMatrix *Src, TMatrix *Dest, float Radius); extern IS_RET DecolorizationWithContrastPreserved(TMatrix *Src, TMatrix *Dest, int Level = 64, float Sigma = 0.05); IS_RET __stdcall PencilDrawing(TMatrix *Src, TMatrix *Dest, int LineLength) { if (Src == NULL || Dest == NULL) return IS_RET_ERR_NULLREFERENCE; if (Src->Data == NULL || Dest->Data == NULL) return IS_RET_ERR_NULLREFERENCE; if (Src->Width != Dest->Width || Src->Height != Dest->Height || Src->Channel != Dest->Channel || Src->Depth != Dest->Depth || Src->WidthStep != Dest->WidthStep) return IS_RET_ERR_PARAMISMATCH; if (Src->Depth != IS_DEPTH_8U || Dest->Depth != IS_DEPTH_8U) return IS_RET_ERR_NOTSUPPORTED; IS_RET Ret = IS_RET_OK; if (Src->Data == Dest->Data) { TMatrix *Clone = NULL; Ret = IS_CloneMatrix(Src, &Clone); if (Ret != IS_RET_OK) return Ret; Ret = PencilDrawing(Clone, Dest, LineLength); IS_FreeMatrix(&Clone); return Ret; } if (Src->Channel == 1) { int Amount = 2 * LineLength + 1; int Width = Src->Width, Height = Src->Height; int X, Y, Z, Index, MaxValue, Sum; float Value; unsigned char *LinePS, *LinePD; TMatrix **Response = (TMatrix **)malloc(8 * sizeof(TMatrix)); TMatrix *ImageEdge = NULL; Ret = IS_CreateMatrix(Width, Height, IS_DEPTH_8U, 1, &ImageEdge); TMatrix *LineShape = NULL; Ret = IS_CreateMatrix(Width, Height, IS_DEPTH_8U, 1, &LineShape); if (Ret != IS_RET_OK) goto Done8; Ret = GuassBlur(Src, ImageEdge, 1); // 去除点噪音 if (Ret != IS_RET_OK) goto Done8; Ret = EdgeCoarse(ImageEdge, ImageEdge); // 两个参数相同对速度无影响，对应论文公式1 if (Ret != IS_RET_OK) goto Done8; for (Z = 0; Z < 8; Z++) { Ret = IS_CreateMatrix(Width, Height, IS_DEPTH_8U, 1, &Response[Z]); if (Ret != IS_RET_OK) goto Done8; Ret = MotionBlur(ImageEdge, Response[Z], Amount, Z * 180.0 / 8, EdgeMode::Smear); // 对应论文公式2 if (Ret != IS_RET_OK) goto Done8; // 各个方向卷积 } for (Y = 0; Y < Height; Y++) { LinePS = ImageEdge->Data + Y * ImageEdge->WidthStep; for (X = 0; X < Width; X++) { MaxValue = 0; for (Z = 0; Z < 8; Z++) { LinePD = (Response[Z]->Data + Y * Response[Z]->WidthStep + X); if (MaxValue < LinePD[0]) { Index = Z; MaxValue = LinePD[0]; } } for (Z = 0; Z < 8; Z++) { LinePD = (Response[Z]->Data + Y * Response[Z]->WidthStep); // 对应公式3 if (Z == Index) LinePD[X] = LinePS[X]; else LinePD[X] = 0; } } } for (Z = 0; Z < 8; Z++) { MotionBlur(Response[Z], Response[Z], Amount, Z * 180.0 / 8, EdgeMode::Smear); // 对应公式S' } for (Y = 0; Y < Height; Y++) { LinePD = LineShape->Data + Y * LineShape->WidthStep; for (X = 0; X < Width; X++) { Sum = 0; for (Z = 0; Z < 8; Z++) { LinePS = (Response[Z]->Data + Y * Response[Z]->WidthStep); Sum += LinePS[X]; } LinePD[X] = (255 - ClampToByte(Sum) * 0.5); // The final pencil stroke map S is obtained by inverting pixel values and mapping them to [0,1]. } } float *HistgramF = (float *)IS_AllocMemory(256 * sizeof(float)); float *HistgramFC = (float *)IS_AllocMemory(256 * sizeof(float)); int *Histgram = (int *)IS_AllocMemory(256 * sizeof(int)); int Ua = 105, Ub = 225, Mud = 90, DeltaB = 9, DeltaD = 11, Omega1 = 76, Omega2 = 22, Omega3 = 2, Iter = 5; for (Y = 0; Y < 256; Y++) { if (Y < Ua || Y > Ub) // 表1中的参数 Value = 0; else Value = 1.0 / (Ub - Ua); HistgramF[Y] = (Omega2 * Value + 1.0 / DeltaB * exp(-(255.0 - Y) / DeltaB) * Omega1 + 1.0 /sqrt(2 * PI * 11) * exp(-(Y - Mud) * (Y - Mud) / (2.0 * DeltaD * DeltaD)) * Omega3) * 0.01; HistgramFC[Y] = HistgramF[Y]; // 拷贝一个备份 } for (Z = 0; Z < Iter; Z++) // 这样的直方图并不平滑，做一点平滑处理 { HistgramFC[0] = (HistgramF[0] + HistgramF[1]) / 2; // 第一点 for (Y = 1; Y < 255; Y++) HistgramFC[Y] = (HistgramF[Y - 1] + HistgramF[Y] + HistgramF[Y + 1]) / 3; // 中间的点 HistgramFC[255] = (HistgramF[254] + HistgramF[255]) / 2; // 最后一点 memcpy(HistgramF, HistgramFC, 256 * sizeof(float)); } for (Y = 0; Y < 256; Y++) Histgram[Y] = HistgramF[Y] * Width * Height; TMatrix *ToneMap = NULL; Ret = IS_CreateMatrix(Width, Height, IS_DEPTH_8U, 1, &ToneMap); if (Ret != IS_RET_OK) goto Done8; Ret = GuassBlur(Src, ToneMap, 1); // Initially, the grayscale input I is slightly Gaussian smoothed. if (Ret != IS_RET_OK) goto Done8; SMLHistgramMaping(ToneMap, ToneMap, Histgram, Histgram, Histgram); // we adjust the tone maps using simple histogram matching in all the three layers and superpose them again. BlendImage(ToneMap, LineShape, Dest, BlendMode::Multiply, 255); // We combine the pencil stroke S and tonal texture T by multiplying the stroke and texture values for each pixel to accentuate important contours Done8: IS_FreeMatrix(&ImageEdge); IS_FreeMatrix(&ToneMap); IS_FreeMatrix(&LineShape); IS_FreeMemory(HistgramF); IS_FreeMemory(HistgramFC); IS_FreeMemory(Histgram); for (Z = 0; Z < 8; Z++)IS_FreeMatrix(&Response[Z]); IS_FreeMemory(*Response); } else { unsigned char *LinePS, *LinePD; int X, Y, Z, Width = Src->Width, Height = Src->Height; TMatrix *Gray = NULL, *GrayC = NULL; IS_RET Ret = IS_CreateMatrix(Src->Width, Src->Height, IS_DEPTH_8U, 1, &Gray); if (Ret != IS_RET_OK) goto Done; Ret = IS_CreateMatrix(Src->Width, Src->Height, IS_DEPTH_8U, Src->Channel, &GrayC); Ret = DecolorizationWithContrastPreserved(Src, Gray); if (Ret != IS_RET_OK) goto Done; PencilDrawing(Gray, Gray, LineLength); for (Y = 0; Y < Height; Y++) { LinePS = Gray->Data + Y * Gray->WidthStep; LinePD = GrayC->Data + Y * GrayC->WidthStep; for (X = 0; X < Width; X++) // 恢复V分量 { LinePD[0] = LinePS[X]; LinePD[1] = LinePS[X]; LinePD[2] = LinePS[X]; LinePD += 3; } } BlendImage(Dest, GrayC, Dest, BlendMode::Luminosity, 255); Done: IS_FreeMatrix(&Gray); IS_FreeMatrix(&GrayC); } }
贴一些处理的效果：
原图 处理结果图
和论文里处理的效果还有不小的差距，这主要是由于最后一步纹理没有做，然后就是还有细节有问题，不过也算有点收获。
提供一个测试小工具： http://files.cnblogs.com/files/Imageshop/PencilDrawing.rar
****************************作者： laviewpbt 时间： 2015.2.21 联系QQ: 33184777 转载请保留本行信息**********************