LIC (Line Integral Convolution) is a well-known texture synthesis technique proposed by Cabral and Leedom at in . 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的实现代码，网络上流传的最为广泛的也是这里的代码，详见：
按照个人的理解，LIC算法他就是把一幅矢量场数据用图像的方式可视化出来，那么对于某一个固定位置的点，他其实是只有当前点的矢量值的，一个带有方向信息的Vector（一般都是归一化后的数据，即矢量的长度为1）。要把该点转换为一个像素值，那么我们需要首先有个基点，对于全局来说就是一幅图，通常情况下，我们会产生一幅和矢量数据大小一样的白噪声图来作为基点图，之后，我们采用卷积的方法，沿着当前点的矢量方向前进一定的距离D，得到新的坐标位置，记录下这个位置在基点图中的像素值，并累加，之后，这个新位置也有他对应的矢量方向，沿着这个矢量方向继续前进，并执行相同的累加操作，直到前进了指定的距离后，再计算累加后的平均值最为可视化后的像素值。
整个流程其实看起来就是沿着某一条线，对线上相关离散点进行卷积，所以严格的说可以叫做折线卷积。同时，为了保证对称性，我们会在卷积起点时也会沿着矢量相反的方向进行卷积，很明显，这个反响卷积的路线和正向卷积的一半来说不会时对称的。
那么谈到代码，我把上面的方案整理成了一个可运行的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) 6 int advDir; ///ADVection DIRection (0: positive; 1: negative) 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; 35 ///for either advection direction/// 36 for (advDir = 0; advDir < 2; advDir++) 37 { 38 ///init the step counter, curve-length measurer, and streamline seed/// 39 advcts = 0; 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/// 48 while (curLen < krnlen && advcts < ADVCTS) 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 { 58 t_acum[advDir] = (advcts == 0) ? 0.0f : t_acum[advDir]; ///this line is indeed unnecessary 59 w_acum[advDir] = (advcts == 0) ? 1.0f : w_acum[advDir]; 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];101 w_acum[advDir] = W_ACUM;102 t_acum[advDir] += texVal * smpWgt;103 104 ///update the step counter and the "current" clip point///105 advcts++;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也就可以了。当流线移动到了新的位置后，我们取这个位置对应的像素来进行卷积。
第二，为了简便，对于流线超出了边缘的部分，我们直接使用了边缘的像素值来代替，这样就造成了边缘的值和原始的效果有所不同。
第三，当遇到某一个点处的X和Y矢量都为0时，理论上讲流线的卷积就应该停止了，而本代码没有考虑这个事情，实际上此时后续取样坐标点就不会在产生任何变化了，流线也就一直停止在这个位置处了。
那么如果要按照原始的算法来该我自己的代码，也时很简单的一件事情的，首先，沿流线和流线相反方向的计算必须要分开了，因为他们可能在不同的位置终止，第二，我们还必须要记录下实际流线中有效的取样点的数量。第三，要注意判断流线取样点的数量是不是为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时，改进的版本和上述快速版本的区别如下;
当流线长度变为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像素，结果如下图：
可以看到此时产生了一个和原始矢量场趋势一致的图，可以认为他就是沿矢量场方向进行了运动模糊的一种特效，我们设计不同的矢量场，就能得到不同的运动模糊效果，那么特别有意义的是，如果这个矢量场时由图像本身的内容生成的，那将使得算法的效果更有意义，此部分我们将在下一章里给出讨论，这里先发布几个图供参考：
简易的油画效果
指纹的增强显示
人物特效显示
本文相关代码下载：
