【算法隨記二】線卷積積分及其在圖像加強和特效方面的應用（一）。

　　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左右。下面我稍微對個人代碼作個解釋。

　　第1、參考源做者的代碼，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