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，還是有點慢的，那么我們看下他的代碼怎么寫的：
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像素，結果如下圖：
可以看到此時產生了一個和原始矢量場趨勢一致的圖，可以認為他就是沿矢量場方向進行了運動模糊的一種特效，我們設計不同的矢量場，就能得到不同的運動模糊效果，那么特別有意義的是，如果這個矢量場時由圖像本身的內容生成的，那將使得算法的效果更有意義，此部分我們將在下一章里給出討論，這里先發布幾個圖供參考：
簡易的油畫效果
指紋的增強顯示
人物特效顯示
本文相關代碼下載：https://files.cnblogs.com/files/Imageshop/Line_Integral_Convolution.rar
