13行代碼實現最快速最高效的積分圖像算法。


  研究圖像到一定程度的人,應該都對積分圖像有所了解,大家在百度或者google中都可以搜索到大量的相關博客,我這里不做多介紹。用積分圖也確實能解決很多實際的問題,比如我博客中的基於局部均方差相關信息的圖像去噪及其在實時磨皮美容算法中的應用 一文我就在網上看到很多人用累計積分圖和乘積積分圖來實現了。不過我瀏覽了很多人的博客,覺得很多人哪怕是圖像方面的似乎還比較牛的人都對這個積分圖的理解還不到位或者說理解有誤,這里有必要借助我微弱的力量再次校正下。

     首先一個普遍的問題就是:積分圖像的大小。你可以看到,90%自己寫的積分圖的作者都是把積分圖像定位為W * H,但是不知道他們有沒有注意到OpenCV里積分圖的相關文檔,我這里貼出OpenCV的官方文檔:

  C++: void integral(InputArray image, OutputArray sum, int sdepth=-1

Parameters:
  • image – Source image as W X H  , 8-bit or floating-point (32f or 64f).
  • sum – Integral image as (W + 1) X (H + 1) , 32-bit integer or floating-point (32f or 64f).
  • sdepth – Desired depth of the integral and the tilted integral images, CV_32SCV_32F, or CV_64F.

 

 

 

     注意到sum這一項的大小了嗎, (W + 1) X (H + 1)而非W X H,為什么呢,我們知道,某個點的積分圖反應的是原圖中此位置左上角所有像素之和,這里是的累加和是不包括這個點像素本身的,那么這樣,對於原圖的第一行和第一列的所有像素,其對應位置的積分圖就應該是0, 這樣考慮到所有的像素,為了能容納最后一列和最后一行的情況,最終的積分圖就應該是 (W + 1) X (H + 1)大小。

     如果你還是希望定義成W X H大小,那么就必須每次判斷你訪問的積分圖的位置,作為寫程序來說,這樣做對程序的性能和代碼的簡潔性都是不好的,並且你稍微不注意就會把代碼寫錯。

     第二,就是積分圖的計算的優化,很多博客也都描述了他們的優化方式,雖然他們都是描述的同一個算法,比如百度上比較靠前的博文: 【圖像處理】快速計算積分圖  中就用下述前兩幅圖描述了他的優化過程:

          --------->    --------->  

                        (1)  最原始方案                 (2)網絡上的優化方案                  (3)我的優化方案

 

     不錯,這樣做已經很不錯了,但是有一個問題就是我們需要多一個大小為W的內存來更新保存每一列相應的累計和,但是我們如果換成行方向呢,如上面的(3)所示,則只需要一個變量來累計行方向上的值,而我們知道單個變量CPU可能會把他放置到寄存器中或者我們可以強制把他聲明為register變量的,而列方向上的內存數據就無法做到這一點了。

     另外,我們申請了積分圖后,是沒有必要先給他來一個清零操作的,因為后續我們會對每一個位置的值都填充0的。

     好,下面我們貼出了優化后的代碼:

 1 void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
 2 {
 3     memset(Integral, 0, (Width + 1) * sizeof(int));                    //    第一行都為0
 4     for (int Y = 0; Y < Height; Y++)
 5     {
 6         unsigned char *LinePS = Src + Y * Stride;
 7         int *LinePL = Integral + Y * (Width + 1) + 1;                  //    上一行位置            
 8         int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;            //    當前位置,注意每行的第一列的值都為0
 9         LinePD[-1] = 0;                                                //    第一列的值為0
10         for (int X = 0, Sum = 0; X < Width; X++)
11         {
12             Sum += LinePS[X];                                        //    行方向累加
13             LinePD[X] = LinePL[X] + Sum;                            //    更新積分圖
14         }
15     }
16 }

     代碼很簡單,但是集成了很多細節信息,只有有心人才能理解。

     在PC上,也許還是可以考慮到SSE優化,要使用SSE,就必須對Sum這個做改造,又要用一個W長度的int型內存記錄,然后對 LinePD[X] = LinePL[X] + Sum; 用SSE優化,一次性處理4個數據據的加法,但是這么做SSE的提速和前面sum那句的處理的耗時那個更厲害,我沒有去驗證了。

     那么上面方法(2)的進一步的語法優化的代碼如下(比原作者那個應該效率會好很多的)

 1 void GetGrayIntegralImage(unsigned char *Src, int *Integral, int Width, int Height, int Stride)
 2 {
 3     int *ColSum = (int *)calloc(Width, sizeof(int));        //    用的calloc函數哦,自動內存清0
 4     memset(Integral, 0, (Width + 1) * sizeof(int));
 5     for (int Y = 0; Y < Height; Y++)
 6     {
 7         unsigned char *LinePS = Src + Y * Stride;
 8         int *LinePL = Integral + Y * (Width + 1) + 1;
 9         int *LinePD = Integral + (Y + 1) * (Width + 1) + 1;
10         LinePD[-1] = 0;
11         for (int X = 0; X < Width; X++)
12         {
13             ColSum[X] += LinePS[X];
14             LinePD[X] = LinePD[X - 1] + ColSum[X];
15         }
16     }
17     free(ColSum);
18 }

  可以看到,在最內層的循環里多了幾個內存變量的訪問,朋友們可以自己測試,我的優化方案在這一塊的耗時只有方案(2)的一半左右,而且代碼更簡潔。

      當我們需要訪問中心點為(x, y),半徑為r的范圍內的矩形像素內的累積值,相應的坐標計算就應該為:

          Integral(x - r, y - r) + Integral(x + r + 1, y + r + 1) - Integral(x - r, y + r + 1) - Integral(x + r + 1, y - r)

  注意上式計算中的加1,這主要是因為積分圖是計算左上角像素的累計值這個特性決定的。

     那么下面我貼出用積分圖實現0(1)的均值模糊的部分代碼:

void BoxBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int Radius)
{
    int *Integral = (int *)malloc((Width + 1) * (Height + 1) * sizeof(int));
    GetGrayIntegralImage(Src, Integral, Width, Height, Stride);
    //#pragma omp parallel for
    for (int Y = 0; Y < Height; Y++)
    {
        int Y1 = max(Y - Radius, 0);
        int Y2 = min(Y + Radius + 1, Height - 1);
        /*int Y1 = Y - Radius;
        int Y2 = Y + Radius + 1;
        if (Y1 < 0) Y1 = 0;
        if (Y2 > Height) Y2 = Height;*/
        int *LineP1 = Integral + Y1 * (Width + 1);
        int *LineP2 = Integral + Y2 * (Width + 1);
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            int X1 = max(X - Radius, 0);
            int X2 = min(X + Radius + 1, Width);
            //    int X1 = X - Radius;
            //    if (X1 < 0) X1 = 0;
            //    int X2 = X + Radius + 1;
            //    if (X2 >= Width) X2 = Width - 1;
            int Sum = LineP2[X2] - LineP1[X2] - LineP2[X1] + LineP1[X1];
            int PixelCount = (X2 - X1) * (Y2 - Y1);
            LinePD[X] = (Sum + (PixelCount >> 1)) / PixelCount;
        }
    }
    free(Integral);
}

  代碼也很簡潔明了,注意下我注釋掉的幾個部分。

  第一:

    //#pragma omp parallel for

  由於進行的積分圖的操作,每個像素點的周邊半徑為r區域內的像素之和的計算就是前后無關的了,因此像素和像素之間的計算就是獨立的了,這樣就可以並行執行,這里用了OpenMP作為簡易並行實現的一個例子而已。

     第二:我原來編程習慣怎是不怎么喜歡用max和min這樣的函數,我總覺得用if這種判斷效率應該會比max或者min高,而實際上卻是max厲害一些,我看了下反匯編,max和min充分利用了條件按傳送指令比如cmovg,cmovl,而一般的if語句很大程度上會調用jmp指令的(當然和常數比較一般編譯也會用條件傳送替代的,比如和0比較),jmp指令時比較耗時的。因此.......

     還有一個注意點就是積分圖的最大X坐標是Width,最大Y坐標是Height, 這點和圖像有所不同,正如上述代碼所示。

     通過積分圖技術實現的均值模糊和之前我在文章解析opencv中Box Filter的實現並提出進一步加速的方案(源碼共享) 中介紹的方式(非SSE優化的代碼)耗時基本差不多,內存占用也差不多,但是積分圖技術有個優勢,就是如果某個算法需要計算同一個圖像的多個半徑的模糊值,則積分圖只需要計算一次,只在眾多的基於多尺度模糊的算法中也是能提速的方案之一。

     上述是普通的累計積分圖,我們常用的還有平方積分圖,還有兩幅圖對應的乘積積分圖,本質上和本例是沒有啥,只需將 Sum += LinePS[X]; 這一句稍作修改。但是有一點,一直是我不太願意用積分圖的主要原因,就是保存積分圖數據的內存區域或者數據類型問題,普通的積分圖還好,就算是每個像素都為255,也只有當圖像大於2800*3000左右時,才會超出int類型所能表達的范圍,如果用uint類型表示,能容納的圖像大小又能提高一倍,一般來說夠用了,但是如果是平方積分圖,int類型在極端情況下只能處理不大於 256*256大小的圖像,這樣在實際中基本上是無用的,因此,可能我們就需要float類型或者int64為來表示。但是,第一,float類型會引來計算速度下降,特別是非PC的環境下。第二,int64會占用更大的內存,大約是8倍的原圖像內存,同時在32位系統上也會帶來速度的下降。

     對於彩色的圖像,處理方式也類似,這里就不贅述了。

     本文對應的代碼下載地址:http://files.cnblogs.com/files/Imageshop/IntegralImage.rar

     

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM