关于Cewu Lu等的《Combining Sketch and Tone for Pencil Drawing Production》一文铅笔画算法的理解和笔录。

相关论文的链接: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

                              原图

                   边缘检测图

                         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 转载请保留本行信息***************

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏专知

【论文推荐】最新七篇图像检索相关论文—草图、Tie-Aware、场景图解析、叠加跨注意力机制、深度哈希、人群估计

2333
来自专栏专知

【论文推荐】最新六篇行人再识别相关论文—特定视角、多目标、双注意匹配网络、联合属性-身份、迁移学习、多通道金字塔型

【导读】专知内容组整理了最近六篇行人再识别(Person Re-Identification)相关文章,为大家进行介绍,欢迎查看! 1. Learning Vi...

7715
来自专栏专知

【干货】Python大数据处理库PySpark实战——使用PySpark处理文本多分类问题

【导读】近日,多伦多数据科学家Susan Li发表一篇博文,讲解利用PySpark处理文本多分类问题的详情。我们知道,Apache Spark在处理实时数据方面...

7.5K8
来自专栏专知

【论文推荐】最新八篇视频描述生成相关论文—在线视频理解、联合定位和描述事件、生成视频、跨模态注意力机制、联合事件检测和描述

【导读】专知内容组整理近期八篇视频描述生成(Video Captioning)相关文章,为大家进行介绍,欢迎查看!

1615
来自专栏DHUtoBUAA

正六边形网格化(Hexagonal Grids)原理与实现

 在路径规划、游戏设计栅格法应用中,正六边形网格不如矩形网格直接和常见,但是正六边形具有自身的应用特点,更适用于一些特殊场景中,比如旷阔的海洋、区域或者太空。...

4075
来自专栏PPV课数据科学社区

【学习】常用的机器学习&数据挖掘知识点

Basis(基础): MSE(Mean Square Error 均方误差),LMS(LeastMean Square 最小均方),LSM(Least Squa...

34312
来自专栏CVer

[计算机视觉论文速递] 2018-03-11

通知:这篇推文有10篇论文速递信息,涉及目标检测、行人重识别Re-ID、图像检索和Zero-Shot Learning等方向 这篇文章本来是在2018-03-1...

4188
来自专栏专知

【论文推荐】最新八篇推荐系统相关论文—可解释推荐、上下文感知推荐系统、异构知识库嵌入、深度强化学习、移动推荐系统

【导读】专知内容组既昨天推出八篇推荐系统相关论文之后,今天为大家又推出八篇推荐系统(Recommendation System)相关论文,欢迎查看!

3723
来自专栏专知

【论文推荐】最新八篇图像检索相关论文—三元组、深度特征图、判别式、卷积特征聚合、视觉-关系知识图谱、大规模图像检索

2463
来自专栏木子昭的博客

K近邻(knn)算法预测电影类型案例1案例2 Facebook入住地点

K近邻思想: 根据你的"邻居们"来确定你的类别 你一觉醒来,不知道自己身在何方里,你能通过计算机定位到周围5个"最近的"邻居,其中有4个身处火星,1个身处月...

3605

扫码关注云+社区