# 【算法随记二】线卷积积分及其在图像增强和特效方面的应用（一）。

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).

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

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

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

```  1 ///        flow imaging (visualization) through Line Integral Convolution     ///
2 void FlowImagingLIC(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage,
3     float*  p_LUT0, float*  p_LUT1, float   krnlen)
4 {
5     int        vec_id;                        ///ID in the VECtor buffer (for the input flow field)
7     int        advcts;                        ///number of ADVeCTion stepS per direction (a step counter)
8     int        ADVCTS = int(krnlen * 3);    ///MAXIMUM number of advection steps per direction to break dead loops
9
10     float    vctr_x;                        ///x-component  of the VeCToR at the forefront point
11     float    vctr_y;                        ///y-component  of the VeCToR at the forefront point
12     float    clp0_x;                        ///x-coordinate of CLiP point 0 (current)
13     float    clp0_y;                        ///y-coordinate of CLiP point 0    (current)
14     float    clp1_x;                        ///x-coordinate of CLiP point 1 (next   )
15     float    clp1_y;                        ///y-coordinate of CLiP point 1 (next   )
16     float    samp_x;                        ///x-coordinate of the SAMPle in the current pixel
17     float    samp_y;                        ///y-coordinate of the SAMPle in the current pixel
18     float    tmpLen;                        ///TeMPorary LENgth of a trial clipped-segment
19     float    segLen;                        ///SEGment   LENgth
20     float    curLen;                        ///CURrent   LENgth of the streamline
21     float    prvLen;                        ///PReVious  LENgth of the streamline
22     float    W_ACUM;                        ///ACcuMulated Weight from the seed to the current streamline forefront
23     float    texVal;                        ///TEXture VALue
24     float    smpWgt;                        ///WeiGhT of the current SaMPle
25     float    t_acum[2];                    ///two ACcUMulated composite Textures for the two directions, perspectively
26     float    w_acum[2];                    ///two ACcUMulated Weighting values   for the two directions, perspectively
27     float*    wgtLUT = NULL;                ///WeiGhT Look Up Table pointing to the target filter LUT
28     float    len2ID = (DISCRETE_FILTER_SIZE - 1) / krnlen;    ///map a curve LENgth TO an ID in the LUT
29     ///for each pixel in the 2D output LIC image///
30     for (int j = 0; j < Height; j++)
31     for (int i = 0; i < Width; i++)
32     {
33         ///init the composite texture accumulators and the weight accumulators///
34         t_acum[0] = t_acum[1] = w_acum[0] = w_acum[1] = 0.0f;
37         {
38             ///init the step counter, curve-length measurer, and streamline seed///
40             curLen = 0.0f;
41             clp0_x = i + 0.5f;
42             clp0_y = j + 0.5f;
43
44             ///access the target filter LUT///
45             wgtLUT = (advDir == 0) ? p_LUT0 : p_LUT1;
46
47             ///until the streamline is advected long enough or a tightly  spiralling center / focus is encountered///
49             {
50                 ///access the vector at the sample///
51                 vec_id = (int(clp0_y) * Width + int(clp0_x)) << 1;
52                 vctr_x = pVectr[vec_id];
53                 vctr_y = pVectr[vec_id + 1];
54
55                 ///in case of a critical point///
56                 if (vctr_x == 0.0f && vctr_y == 0.0f)
57                 {
60                     break;
61                 }
62
63                 ///negate the vector for the backward-advection case///
64                 vctr_x = (advDir == 0) ? vctr_x : -vctr_x;
65                 vctr_y = (advDir == 0) ? vctr_y : -vctr_y;
66
67                 ///clip the segment against the pixel boundaries --- find the shorter from the two clipped segments///
68                 ///replace  all  if-statements  whenever  possible  as  they  might  affect the computational speed///
69                 segLen = LINE_SQUARE_CLIP_MAX;
70                 segLen = (vctr_x < -VECTOR_COMPONENT_MIN) ? (int(clp0_x) - clp0_x) / vctr_x : segLen;
71                 segLen = (vctr_x >  VECTOR_COMPONENT_MIN) ? (int(int(clp0_x) + 1.5f) - clp0_x) / vctr_x : segLen;
72                 segLen = (vctr_y < -VECTOR_COMPONENT_MIN) ?
73                     (((tmpLen = (int(clp0_y) - clp0_y) / vctr_y)  <  segLen) ? tmpLen : segLen)
74                     : segLen;
75                 segLen = (vctr_y >  VECTOR_COMPONENT_MIN) ?
76                     (((tmpLen = (int(int(clp0_y) + 1.5f) - clp0_y) / vctr_y)  <  segLen) ? tmpLen : segLen)
77                     : segLen;
78
79                 ///update the curve-length measurers///
80                 prvLen = curLen;
81                 curLen += segLen;
82                 segLen += 0.0004f;
83
84                 ///check if the filter has reached either end///
85                 segLen = (curLen > krnlen) ? ((curLen = krnlen) - prvLen) : segLen;
86
87                 ///obtain the next clip point///
88                 clp1_x = clp0_x + vctr_x * segLen;
89                 clp1_y = clp0_y + vctr_y * segLen;
90
91                 ///obtain the middle point of the segment as the texture-contributing sample///
92                 samp_x = (clp0_x + clp1_x) * 0.5f;
93                 samp_y = (clp0_y + clp1_y) * 0.5f;
94
95                 ///obtain the texture value of the sample///
96                 texVal = pNoise[int(samp_y) * Width + int(samp_x)];
97
98                 ///update the accumulated weight and the accumulated composite texture (texture x weight)///
99                 W_ACUM = wgtLUT[int(curLen * len2ID)];
100                 smpWgt = W_ACUM - w_acum[advDir];
102                 t_acum[advDir] += texVal * smpWgt;
103
104                 ///update the step counter and the "current" clip point///
106                 clp0_x = clp1_x;
107                 clp0_y = clp1_y;
108
109                 ///check if the streamline has gone beyond the flow field///
110                 if (clp0_x < 0.0f || clp0_x >= Width || clp0_y < 0.0f || clp0_y >= Height)  break;
111             }
112         }
113
114         ///normalize the accumulated composite texture///
115         texVal = (t_acum[0] + t_acum[1]) / (w_acum[0] + w_acum[1]);
116
117         ///clamp the texture value against the displayable intensity range [0, 255]
118         texVal = (texVal <   0.0f) ? 0.0f : texVal;
119         texVal = (texVal > 255.0f) ? 255.0f : texVal;
120         pImage[j * Width + i] = (unsigned char)texVal;
121     }
122 }```

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

正常的结果　　　　　　　　　　　　　　　　　　　　　　　　　　取消了segLen += 0.0004f的结果

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

``` 1 void FastFlowImagingLIC(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage, float   krnlen)
2 {
3     float Step = 0.333333f;                         //    每次前进0.33333像素的流线距离
4     int LoopAmount = int(krnlen / Step);            //    流线是左右对称的
5     if (LoopAmount <= 0)    LoopAmount = 1;
6     int Weight = LoopAmount * 2;
7     for (int Y = 0; Y < Height; Y++)
8     {
9         unsigned char *LinePD = pImage + Y * Width;
10         for (int X = 0; X < Width; X++)
11         {
12             float PosX_P = X + 0.5f, PosY_P = Y + 0.5f, PosX_N = X + 0.5f, PosY_N = Y + 0.5f;
13             int Sum = 0;
14             for (int Z = 0; Z < LoopAmount; Z++)
15             {
16                 int Index_P = (IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1)) << 1;
17                 int Index_N = (IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1)) << 1;
18                 PosX_P += pVectr[Index_P] * 0.333333f;    PosY_P += pVectr[Index_P + 1] * 0.333333f;            //    X和Y都是归一化的，所以*0.333333f后流线距离也就是0.25了
19                 PosX_N -= pVectr[Index_N] * 0.333333f;    PosY_N -= pVectr[Index_N + 1] * 0.333333f;
20                 Index_P = IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1);
21                 Index_N = IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1);
22                 Sum += pNoise[Index_P] + pNoise[Index_N];
23             }
24             LinePD[X] = (unsigned char)(Sum / Weight);
25         }
26     }
27 }```

一共也只有二十几行，效果如下：

原始的效果　　　　　　　　　　　　　　　　　　　　　　　　　　　改动的效果

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

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

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

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

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

``` 1 void FastFlowImagingLIC_Appr(int  Width, int  Height, float*  pVectr, unsigned char*  pNoise, unsigned char*  pImage, float   krnlen)
2 {
3     float Step = 0.333333f;                                    //    每次前进0.33333像素的流线距离
4     int LoopAmount = int(krnlen / Step);            //    流线是左右对称的
5     if (LoopAmount <= 0)    LoopAmount = 1;
6     int Weight = LoopAmount * 2;
7     for (int Y = 0; Y < Height; Y++)
8     {
9         unsigned char *LinePD = pImage + Y * Width;
10         float *LinePV = pVectr + Y * Width * 2;
11         for (int X = 0; X < Width; X++, LinePV += 2)
12         {
13             float PosX_P = X + 0.5f, PosY_P = Y + 0.5f, PosX_N = X + 0.5f, PosY_N = Y + 0.5f;
14             float VectorX = LinePV[0] * 0.333333f, VectorY = LinePV[1] * 0.333333f;
15             int Sum = 0;
16             for (int Z = 0; Z < LoopAmount; Z++)
17             {
18                 PosX_P += VectorX;    PosY_P += VectorY;            //    X和Y都是归一化的，所以*0.333333f后流线距离也就是0.25了
19                 PosX_N -= VectorX;    PosY_N -= VectorY;
20                 int Index_P = IM_ClampI(PosY_P, 0, Height - 1) * Width + IM_ClampI(PosX_P, 0, Width - 1);
21                 int Index_N = IM_ClampI(PosY_N, 0, Height - 1) * Width + IM_ClampI(PosX_N, 0, Width - 1);
22                 Sum += pNoise[Index_P] + pNoise[Index_N];
23             }
24             LinePD[X] = (unsigned char)(Sum / Weight);
25         }
26     }
27 }```

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

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

但是这种近似的速度却能大幅提高，大概只需要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像素，结果如下图：

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

简易的油画效果

指纹的增强显示

人物特效显示

0 条评论

• ### 远程办公经验为0，如何将日常工作平滑过度到线上?

我是一名创业者，我的公司（深圳市友浩达科技有限公司）在2018年8月8日开始运营，现在还属于微型公司。这个春节假期，我一直十分关注疫情动向，也非常关心其对公司带来的影响。

• ### 数据中台，概念炒作还是另有奇效？ | TVP思享

作者简介：史凯，花名凯哥，腾讯云最具价值专家TVP，ThoughtWorks数据智能业务总经理。投身于企业数字化转型工作近20年。2000年初，在IBM 研发企业级中间件，接着加入埃森哲，为大型企业提供信息化架构规划，设计，ERP，云平台，数据仓库构建等技术咨询实施服务，随后在EMC负责企业应用转型业务，为企业提供云迁移，应用现代化服务。现在专注于企业智能化转型领域，是数据驱动的数字化转型的行业布道者，数据中台的推广者，精益数据创新体系的创始人，2019年荣获全球Data IQ 100人的数据赋能者称号，创业邦卓越生态聚合赋能官TOP 5。2019年度数字化转型专家奖。打造了行业第一个数据创新的数字化转型卡牌和工作坊。创建了精益数据创新方法论体系构建数据驱动的智能企业，并在多个企业验证成功，正在向国内外推广。

• ### 扩展 Kubernetes 之 CRI

使用 cri-containerd 的调用流程更为简洁, 省去了上面的调用流程的 1，2 两步

• ### 扩展 Kubernetes 之 Kubectl Plugin

kubectl 功能非常强大, 常见的命令使用方式可以参考 kubectl --help，或者这篇文章

• ### 多种登录方式定量性能测试方案

最近接到到一个测试任务，某服务提供了两种登录方式：1、账号密码登录；2、手机号+验证码登录。要对这两种登录按照一定的比例进行压测。

• ### 线程安全类在性能测试中应用

首先验证接口参数签名是否正确，然后加锁去判断订单信息和状态，处理用户增添VIP时间事务，成功之后释放锁。锁是针对用户和订单的分布式锁，使用方案是用的redis。

• ### 使用CDN(jsdelivr) 优化博客访问速度

PS: 此篇文章适用于 使用 Github pages 或者 coding pages 的朋友,其他博客也类似.

• ### 扩展 Kubernetes 之 CNI

Network Configuration 是 CNI 输入参数中最重要当部分, 可以存储在磁盘上

• ### 聚焦【技术应变力】云加社区沙龙online重磅上线！

云加社区结合特殊时期热点，挑选备受关注的音视频流量暴增、线下业务快速转线上、紧急上线防疫IoT应用等话题，邀请众多业界专家，为大家提供连续十一天的干货分享。从视野、预判、应对等多角度，帮助大家全面提升「技术应变力」！

• ### 京东购物小程序购物车性能优化实践

它是小程序开发工具内置的一个可视化监控工具，能够在 OS 级别上实时记录系统资源的使用情况。