首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >线性卷积积分及其在图像增强和特效方面的应用

线性卷积积分及其在图像增强和特效方面的应用

作者头像
OpenCV学堂
发布2019-05-14 10:22:58
1K0
发布2019-05-14 10:22:58
举报

 LIC (Line Integral Convolution) is a well-known texture synthesis technique proposed by Cabral and Leedom [33] at Lawrence Livermore National Laboratory in ACM SigGraph 93. It employs a low-pass filter to convolve an input noise texture along pixel-centered symmetrically bi-directional streamlines to exploit spatial correlation in the flow direction. LIC provides a global dense representation of the flow, emulating what happens when a rectangular area of massless fine sand is blown by strong wind (Figure 1).

以上内容和图片摘自

http://www.zhanpingliu.org/Research/FlowVis/LIC/LIC.htm,

这是一个中国人的网站,并且还是一个很老的网站。

关于LIC的实现代码,网络上流传的最为广泛的也是这里的代码,详见:http://www.zhanpingliu.org/Research/FlowVis/LIC/LIC_Source.htm

按照个人的理解,LIC算法他就是把一幅矢量场数据用图像的方式可视化出来,那么对于某一个固定位置的点,他其实是只有当前点的矢量值的,一个带有方向信息的Vector(一般都是归一化后的数据,即矢量的长度为1)。要把该点转换为一个像素值,那么我们需要首先有个基点,对于全局来说就是一幅图,通常情况下,我们会产生一幅和矢量数据大小一样的白噪声图来作为基点图,之后,我们采用卷积的方法,沿着当前点的矢量方向前进一定的距离D,得到新的坐标位置,记录下这个位置在基点图中的像素值,并累加,之后,这个新位置也有他对应的矢量方向,沿着这个矢量方向继续前进,并执行相同的累加操作,直到前进了指定的距离后,再计算累加后的平均值最为可视化后的像素值。

整个流程其实看起来就是沿着某一条线,对线上相关离散点进行卷积,所以严格的说可以叫做折线卷积。同时,为了保证对称性,我们会在卷积起点时也会沿着矢量相反的方向进行卷积,很明显,这个反响卷积的路线和正向卷积的一半来说不会时对称的。

那么谈到代码,我把上面zhanpingliu的方案整理成了一个可运行的C++方案,并未做任何的优化,对一幅512*512的灰度图做测试,由矢量数据生成图像的耗时约为145ms,还是有点慢的,那么我们看下他的代码怎么写的:

///        flow imaging (visualization) through Line Integral Convolution     ///
void FlowImagingLIC(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage,
    float*  p_LUT0, float*  p_LUT1, float   krnlen)
{
    int        vec_id;                        ///ID in the VECtor buffer (for the input flow field)
    int        advDir;                        ///ADVection DIRection (0: positive;  1: negative)
    int        advcts;                        ///number of ADVeCTion stepS per direction (a step counter)
    int        ADVCTS = int(krnlen * 3);    ///MAXIMUM number of advection steps per direction to break dead loops    

    float    vctr_x;                        ///x-component  of the VeCToR at the forefront point
    float    vctr_y;                        ///y-component  of the VeCToR at the forefront point
    float    clp0_x;                        ///x-coordinate of CLiP point 0 (current)
    float    clp0_y;                        ///y-coordinate of CLiP point 0    (current)
    float    clp1_x;                        ///x-coordinate of CLiP point 1 (next   )
    float    clp1_y;                        ///y-coordinate of CLiP point 1 (next   )
    float    samp_x;                        ///x-coordinate of the SAMPle in the current pixel
    float    samp_y;                        ///y-coordinate of the SAMPle in the current pixel
    float    tmpLen;                        ///TeMPorary LENgth of a trial clipped-segment
    float    segLen;                        ///SEGment   LENgth
    float    curLen;                        ///CURrent   LENgth of the streamline
    float    prvLen;                        ///PReVious  LENgth of the streamline        
    float    W_ACUM;                        ///ACcuMulated Weight from the seed to the current streamline forefront
    float    texVal;                        ///TEXture VALue
    float    smpWgt;                        ///WeiGhT of the current SaMPle
    float    t_acum[2];                    ///two ACcUMulated composite Textures for the two directions, perspectively
    float    w_acum[2];                    ///two ACcUMulated Weighting values   for the two directions, perspectively
    float*    wgtLUT = NULL;                ///WeiGhT Look Up Table pointing to the target filter LUT
    float    len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen;    ///map a curve LENgth TO an ID in the LUT
    ///for each pixel in the 2D output LIC image///
    for (int j = 0; j < Height; j++)
    for (int i = 0; i < Width; i++)
    {
        ///init the composite texture accumulators and the weight accumulators///
        t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;
        ///for either advection direction///
        for (advDir = 0; advDir < 2; advDir++)
        {
            ///init the step counter, curve-length measurer, and streamline seed///
            advcts = 0;
            curLen = 0.0f;
            clp0_x = i + 0.5f;
            clp0_y = j + 0.5f;

            ///access the target filter LUT///
            wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;

            ///until the streamline is advected long enough or a tightly  spiralling center / focus is encountered///
            while (curLen < krnlen && advcts < ADVCTS)
            {
                ///access the vector at the sample///
                vec_id = (int(clp0_y) * Width + int(clp0_x)) << 1;
                vctr_x = pVectr[vec_id];
                vctr_y = pVectr[vec_id + 1];

                ///in case of a critical point///
                if (vctr_x == 0.0f && vctr_y == 0.0f)
                {
                    t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir];           ///this line is indeed unnecessary
                    w_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir];
                    break;
                }

                ///negate the vector for the backward-advection case///
                vctr_x = (advDir == 0) ? vctr_x : -vctr_x;
                vctr_y = (advDir == 0) ? vctr_y : -vctr_y;

                ///clip the segment against the pixel boundaries --- find the shorter from the two clipped segments///
                ///replace  all  if-statements  whenever  possible  as  they  might  affect the computational speed///
                segLen = LINE_SQUARE_CLIP_MAX;
                segLen = (vctr_x < -VECTOR_COMPONENT_MIN) ? (int(clp0_x) - clp0_x) / vctr_x : segLen;
                segLen = (vctr_x >  VECTOR_COMPONENT_MIN) ? (int(int(clp0_x) + 1.5f) - clp0_x) / vctr_x : segLen;
                segLen = (vctr_y < -VECTOR_COMPONENT_MIN) ?
                    (((tmpLen = (int(clp0_y) - clp0_y) / vctr_y)  <  segLen) ? tmpLen : segLen)
                    : segLen;
                segLen = (vctr_y >  VECTOR_COMPONENT_MIN) ?
                    (((tmpLen = (int(int(clp0_y) + 1.5f) - clp0_y) / vctr_y)  <  segLen) ? tmpLen : segLen)
                    : segLen;

                ///update the curve-length measurers///
                prvLen = curLen;
                curLen += segLen;
                segLen += 0.0004f;

                ///check if the filter has reached either end///
                segLen = (curLen > krnlen) ? ((curLen = krnlen) - prvLen) : segLen;

                ///obtain the next clip point///
                clp1_x = clp0_x + vctr_x * segLen;
                clp1_y = clp0_y + vctr_y * segLen;

                ///obtain the middle point of the segment as the texture-contributing sample///
                samp_x = (clp0_x + clp1_x) * 0.5f;
                samp_y = (clp0_y + clp1_y) * 0.5f;

                ///obtain the texture value of the sample///
                texVal = pNoise[int(samp_y) * Width + int(samp_x)];

                ///update the accumulated weight and the accumulated composite texture (texture x weight)///
                W_ACUM = wgtLUT[int(curLen * len2ID)];
                smpWgt = W_ACUM - w_acum[advDir];
                w_acum[advDir] = W_ACUM;
                t_acum[advDir] += texVal * smpWgt;

                ///update the step counter and the "current" clip point///
                advcts++;
                clp0_x = clp1_x;
                clp0_y = clp1_y;

                ///check if the streamline has gone beyond the flow field///
                if (clp0_x < 0.0f || clp0_x >= Width || clp0_y < 0.0f || clp0_y >= Height)  break;
            }
        }

        ///normalize the accumulated composite texture///
        texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);

        ///clamp the texture value against the displayable intensity range [0, 255]
        texVal = (texVal <   0.0f) ? 0.0f : texVal;
        texVal = (texVal > 255.0f) ? 255.0f : texVal;
        pImage[j * Width + i] = (unsigned char)texVal;
    }
}

一百多行代码,说真的,我看的有点晕,为什么呢,其一就是这个代码写的真心的不是很好,里面由很多冗余的部分,虽然注释确实是写的不少,难以理解的其实也就是第69到77行,计算seglen的部分,这里的4行判断同时能满足的也就两条,他全部放在一起了,他这里计算最后取样坐标居然不是以矢量数据为主要依据,而是前一次的取样位置。我有点不明白是为什么。同时,那个第82行也非常重要,如果没有那一行,结果就完全不正确,谁知道这是为什么。

我试着按照自己的理解,并参考上面的代码最这个过程做了修改,使得代码看起来简洁而且运行速度也能快一点。具体如下:

void FastFlowImagingLIC(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage, float   krnlen)
{
    float Step = 0.333333f;                         //    每次前进0.33333像素的流线距离                                        
    int LoopAmount = int(krnlen / Step);            //    流线是左右对称的
    if (LoopAmount <= 0)    LoopAmount = 1;
    int Weight = LoopAmount * 2;
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePD = pImage + Y * Width;
        for (int X = 0; X < Width; X++)
        {
            float PosX_P = X + 0.5f, PosY_P = Y + 0.5f, PosX_N = X + 0.5f, PosY_N = Y + 0.5f;
            int Sum = 0;
            for (int Z = 0; Z < LoopAmount; Z++)
            {
                int Index_P = (IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1)) << 1;
                int Index_N = (IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1)) << 1;
                PosX_P += pVectr[Index_P] * 0.333333f;    PosY_P += pVectr[Index_P + 1] * 0.333333f;            //    X和Y都是归一化的,所以*0.333333f后流线距离也就是0.25了
                PosX_N -= pVectr[Index_N] * 0.333333f;    PosY_N -= pVectr[Index_N + 1] * 0.333333f;
                Index_P = IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1);
                Index_N = IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1);
                Sum += pNoise[Index_P] + pNoise[Index_N];
            }
            LinePD[X] = (unsigned char)(Sum / Weight);
        }
    }
}

一共也只有二十几行,效果如下:

除了边缘的部分,整体基本没有什么大的区别,速度上大概能到100ms左右。下面我稍微对我的代码做个解释。

第一、参考源作者的代码,ADVCTS 设置为int(krnlen * 3),即计算次数为流线长度的3倍,我们可以认为他一次沿着流线走1/3像素的距离,因此,我们这里取Step = 0.33333f,接着,我们的流线的起点就是要计算的当前点的坐标,按照当前点的矢量方向或反矢量方向前进1/3像素,因为这个算法中我们要求Vector变量在使用之前必须是归一化的,所以X和Y坐标各自乘以Step也就可以了。当流线移动到了新的位置后,我们取这个位置对应的像素来进行卷积。

第二,为了简便,对于流线超出了边缘的部分,我们直接使用了边缘的像素值来代替,这样就造成了边缘的值和原始的效果有所不同。

第三,当遇到某一个点处的X和Y矢量都为0时,理论上讲流线的卷积就应该停止了,而本代码没有考虑这个事情,实际上此时后续取样坐标点就不会在产生任何变化了,流线也就一直停止在这个位置处了。

那么如果要按照原始的算法来该我自己的代码,也时很简单的一件事情的,首先,沿流线和流线相反方向的计算必须要分开了,因为他们可能在不同的位置终止,第二,我们还必须要记录下实际流线中有效的取样点的数量。第三,要注意判断流线取样点的数量是不是为0。

另外,无论是原始的代码,还是改动后的,其实取样这一块都可以进一步加以改进,可以看到,取矢量值时我们得到的矢量坐标是浮点数,在基点图中取样的坐标点也是浮点数,而我们都直接把他们取整后在计算坐标的,如果不考虑耗时,而要获得更好的效果,就应该对他们插值,比如可以用双线性插值做个简单的优化,这样能获得更好的视觉效果。

还有一种近似的方法,就是我们考虑对于一个特定的点,卷积的方向就一直不改变,就以当前点的矢量方向执行严格的线卷积,当然这个时候,对于那些具有强烈矢量变换的区域,这种方法就有点效果问题了,但是如果卷积的长度不大,一半来说区别不大,如下所示:

void FastFlowImagingLIC_Appr(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage, float   krnlen)
{
    float Step = 0.333333f;                                    //    每次前进0.33333像素的流线距离                                        
    int LoopAmount = int(krnlen / Step);            //    流线是左右对称的
    if (LoopAmount <= 0)    LoopAmount = 1;
    int Weight = LoopAmount * 2;
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePD = pImage + Y * Width;
        float *LinePV = pVectr + Y * Width * 2;
        for (int X = 0; X < Width; X++, LinePV += 2)
        {
            float PosX_P = X + 0.5f, PosY_P = Y + 0.5f, PosX_N = X + 0.5f, PosY_N = Y + 0.5f;
            float VectorX = LinePV[0] * 0.333333f, VectorY = LinePV[1] * 0.333333f;
            int Sum = 0;
            for (int Z = 0; Z < LoopAmount; Z++)
            {
                PosX_P += VectorX;    PosY_P += VectorY;            //    X和Y都是归一化的,所以*0.333333f后流线距离也就是0.25了
                PosX_N -= VectorX;    PosY_N -= VectorY;
                int Index_P = IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1);
                int Index_N = IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1);
                Sum += pNoise[Index_P] + pNoise[Index_N];
            }
            LinePD[X] = (unsigned char)(Sum / Weight);
        }
    }
}

比如,当流线长度为10时,改进的版本和上述快速版本的区别如下;

当流线长度变为30时,区别如下:

我们注意他们的正中心部位,可以看到快速的版本中间出现了明显的紊乱和不同,这主要是中心部分矢量场的紊动特别频繁,相邻位置处的差异比骄大,签署的近似就会带来不太正确的结果。

但是这种近似的速度却能大幅提高,大概只需要35ms,而且还有一个好处就是这个算法可以使用SSE去进行优化了,优化后能做到18m左右。

在原始代码里,有p_LUT0及p_LUT1两个查找表,并且是线性的,所以在这里其实是毫无作用的,但是这说明作者还是想到了,这个积分可以不是普通的均值积分,也可以是类似高斯这种权重随流线距离起点距离成反比的样式的啊。这就可以自由发挥了。有兴趣的朋友可以自己尝试下。

前面一直说的基点图,在这里使用的白噪音的图像,源代码里由如下算法:

// make white noise as the LIC input texture     ///
void MakeWhiteNoise(int  Width, int  Height, unsigned char*  pNoise)
{
    srand(100000);
    for (int j = 0; j < Height; j++)
    for (int X = 0; X < Width; X++)
    {
        int  r = rand();
        r = ((r & 0xff) + ((r & 0xff00) >> 8)) & 0xff;
        pNoise[j * Width + X] = (unsigned char)r;
    }
}

rand()函数其实是个很耗时的函数,如果真的要处理大图像,上述代码的效率是很低的,工程应用上还有很多其他技巧来实现上述类似的效果的,这里不赘述。

如果我们使用另外一幅图像来代替这个白噪声图像,那么出来的结果是什么样呢,我们做个测试,输入一个lena图,流线长度设计为30像素,结果如下图:

可以看到此时产生了一个和原始矢量场趋势一致的图,可以认为他就是沿矢量场方向进行了运动模糊的一种特效,我们设计不同的矢量场,就能得到不同的运动模糊效果,那么特别有意义的是,如果这个矢量场时由图像本身的内容生成的,那将使得算法的效果更有意义,此部分我们将在下一章里给出讨论,这里先发布几个图供参考:

简易油画效果!

指纹增强!

图像增强

作者文章相关代码下载:

https://files.cnblogs.com/files/Imageshop/Line_Integral_Convolution.rar

作者是图像处理算法优化与图像增强方面的大神,其代表作品Imageshop

博客地址:

https://www.cnblogs.com/Imageshop/

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-05-06,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 OpenCV学堂 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
相关产品与服务
图像处理
图像处理基于腾讯云深度学习等人工智能技术,提供综合性的图像优化处理服务,包括图像质量评估、图像清晰度增强、图像智能裁剪等。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档