任意半徑局部直方圖類算法在PC中快速實現的框架。

     在圖像處理中，局部算法一般來說，在很大程度上會獲得比全局算法更為好的效果，因為他考慮到了圖像領域像素的信息，而很多局部算法可以借助於直方圖獲得加速。同時，一些常規的算法，比如中值濾波、最大值濾波、最小值濾波、表面模糊等等都可以通過局部直方圖進行加速。而傳統的獲取局部直方圖計算量很大，特別是半徑增加時，耗時會成平方關系增加。一些局部算法只有在半徑較大時才會獲得很好的效果，因此，必須找到一種合適的加速計算局部直方圖的方式。

      在參考Median Filter in Constant Time.pdf一文附帶的C的代碼的基礎上，本文提出了基於SSE加速的恆長任意半徑局部直方圖獲取技術，可以大大加速算法的計算時間，特別是大半徑時的提速更為明顯。

      主要的優化思路是，沿着列方向一行一行的更行整行的列直方圖，新的一行對應的列直方圖更新時只需要減去已經不再范圍內的那個像素同時加入新進入的像素的直方圖信息。之后，對於一行中的第一個像素點，累加半徑輻射范圍內的列直方圖，得到改點的局部直方圖，對於行中的其他的像素，則類似於更新行直方圖，先減去不在范圍內那列的列直方圖，然后加上移入范圍內的列直方圖。由於采用了基於SSE函數的加速過程，直方圖想加和相減的速度較普通的加減法有了10倍以上的提速，因此大大的提高了整體的實用性。

       具體的過程我用代碼加以說明：

  　　1、一些公用的內存分配過程

    TMatrix *Row = NULL, *Col = NULL; unsigned char *LinePS, *LinePD; int  X, Y, K, Width = Src->Width, Height = Src->Height; int *RowOffset, *ColOffSet; unsigned short *ColHist    = (unsigned short *)IS_AllocMemory(256 * (Width + 2 * Radius) * sizeof(unsigned short), true); if (ColHist == NULL) {Ret = IS_RET_ERR_OUTOFMEMORY; goto Done8;} unsigned short *Hist    = (unsigned short *)IS_AllocMemory(256 * sizeof(unsigned short), true); if (Hist == NULL) {Ret = IS_RET_ERR_OUTOFMEMORY; goto Done8;} Ret = GetValidCoordinate(Width, Height, Radius, Radius, Radius, Radius, Edge, &Row, &Col);        // 獲取坐標偏移量
    if (Ret != IS_RET_OK) goto Done8;

　　其中的ColHist用於保存一行像素對應的列直方圖 ，注意這里的行是用的擴展后的行的大小即：Width + 2 * Radius。IS_AllocMemory是個內部使用了_mm_malloc定義的內存分配函數，主要是考慮SSE函數的16字節對齊問題。

      Hist變量用於保存每個像素點的局部直方圖數據，任何基於局部直方圖技術的函數最終都演變為對於該函數進行各種各樣的計算。    

      GetValidCoordinate是一個用於輔助邊界處像素點處理的函數，具體可詳見附件中給出的代碼。

      2、更新一行像素的列直方圖

for (Y = 0; Y < Height; Y++) { if (Y == 0)                                            // 第一行的列直方圖,要重頭計算
 { for (K = -Radius; K <= Radius; K++) { LinePS = Src->Data + ColOffSet[K] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++) { ColHist[X * 256 + LinePS[RowOffset[X]]]++; } } } else                                                // 其他行的列直方圖，更新就可以了
 { LinePS = Src->Data + ColOffSet[Y - Radius - 1] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++)        // 刪除移出范圍內的那一行的直方圖數據
 { ColHist[X * 256 + LinePS[RowOffset[X]]]--; } LinePS = Src->Data + ColOffSet[Y + Radius] * Src->WidthStep; for (X = -Radius; X < Width + Radius; X++)        // 增加進入范圍內的那一行的直方圖數據
 { ColHist[X * 256 + LinePS[RowOffset[X]]]++; } }
　　//　　依次獲取一行每個像素的局部直方圖
 //　　根據局部直方圖獲的結果
}

　　可見，這部分和普通的局部優化方式類似，沒有什么特殊的地方。

　　3、依次獲取一行每個像素的局部直方圖

    for (Y = 0; Y < Height; Y++) { 　　　　　//　　更新一行像素的列直方圖 memset(Hist, 0, 256 * sizeof(unsigned short));        // 每一行直方圖數據清零先 LinePS = Src->Data + Y * Src->WidthStep; LinePD = Dest->Data + Y * Dest->WidthStep; for (X = 0; X < Width; X++) { if (X == 0) { for (K = -Radius; K <= Radius; K++)            // 行第一個像素，需要重新計算 
                    HistgramAddShort(ColHist + K * 256, Hist); } else { 　　/* 　HistgramAddShort(ColHist + RowOffset[X + Radius] * 256, Hist); 　　HistgramSubShort(ColHist + RowOffset[X - Radius - 1] * 256, Hist); 　　　　　　 */ HistgramSubAddShort(ColHist + RowOffset[X - Radius - 1] * 256, ColHist + RowOffset[X + Radius] * 256, Hist);  // 行內其他像素，依次刪除和增加就可以了
 }

　　　　　　　　//　　根據局部直方圖獲的結果

 LinePS++; LinePD++; } }

　　上面處理的過程其實和2的過程的優化道理是類似的，只不過一個是行方向，一個是列方向，聰明者自然能明白，稍微愚鈍者請自己多多斟酌，自然有豁然開朗的時刻。

　　4、 根據局部直方圖獲的結果

　　根據不同的算法需求，結合局部直方圖信息來獲取結果，比如最大值算法可以用如下方式獲得：

    for (K = 255; K >= 0; K--)
    {
        if (Hist[K] != 0)
        {
            LinePD[X] = K;
            break;
        }
    }

    　關於直方圖累加的代碼如下：

/// <summary>
/// 無符號短整形直方圖數據相加，Y = X + Y， 整理時間2014.12.28; 
/// </summary>
/// <param name="X">加數。</param>
/// <param name="Y">被加數，結果保存於該數中。</param>
/// <remarks>使用了SSE優化。</remarks>
void HistgramAddShort(unsigned short *X, unsigned short *Y)
{
    *(__m128i*)(Y + 0)        =    _mm_add_epi16(*(__m128i*)&Y[0],        *(__m128i*)&X[0]);        //    不要想着用自己寫的匯編超過他的速度了，已經試過了
    *(__m128i*)(Y + 8)        =    _mm_add_epi16(*(__m128i*)&Y[8],        *(__m128i*)&X[8]);
    *(__m128i*)(Y + 16)        =    _mm_add_epi16(*(__m128i*)&Y[16],    *(__m128i*)&X[16]);
    *(__m128i*)(Y + 24)        =    _mm_add_epi16(*(__m128i*)&Y[24],    *(__m128i*)&X[24]);
    *(__m128i*)(Y + 32)        =    _mm_add_epi16(*(__m128i*)&Y[32],    *(__m128i*)&X[32]);
    *(__m128i*)(Y + 40)        =    _mm_add_epi16(*(__m128i*)&Y[40],    *(__m128i*)&X[40]);
    *(__m128i*)(Y + 48)        =    _mm_add_epi16(*(__m128i*)&Y[48],    *(__m128i*)&X[48]);
    *(__m128i*)(Y + 56)        =    _mm_add_epi16(*(__m128i*)&Y[56],    *(__m128i*)&X[56]);
    *(__m128i*)(Y + 64)        =    _mm_add_epi16(*(__m128i*)&Y[64],    *(__m128i*)&X[64]);
    *(__m128i*)(Y + 72)        =    _mm_add_epi16(*(__m128i*)&Y[72],    *(__m128i*)&X[72]);
    *(__m128i*)(Y + 80)        =    _mm_add_epi16(*(__m128i*)&Y[80],    *(__m128i*)&X[80]);
    *(__m128i*)(Y + 88)        =    _mm_add_epi16(*(__m128i*)&Y[88],    *(__m128i*)&X[88]);
    *(__m128i*)(Y + 96)        =    _mm_add_epi16(*(__m128i*)&Y[96],    *(__m128i*)&X[96]);    
    *(__m128i*)(Y + 104)    =    _mm_add_epi16(*(__m128i*)&Y[104],    *(__m128i*)&X[104]);
    *(__m128i*)(Y + 112)    =    _mm_add_epi16(*(__m128i*)&Y[112],    *(__m128i*)&X[112]);
    *(__m128i*)(Y + 120)    =    _mm_add_epi16(*(__m128i*)&Y[120],    *(__m128i*)&X[120]);
    *(__m128i*)(Y + 128)    =    _mm_add_epi16(*(__m128i*)&Y[128],    *(__m128i*)&X[128]);
    *(__m128i*)(Y + 136)    =    _mm_add_epi16(*(__m128i*)&Y[136],    *(__m128i*)&X[136]);
    *(__m128i*)(Y + 144)    =    _mm_add_epi16(*(__m128i*)&Y[144],    *(__m128i*)&X[144]);
    *(__m128i*)(Y + 152)    =    _mm_add_epi16(*(__m128i*)&Y[152],    *(__m128i*)&X[152]);
    *(__m128i*)(Y + 160)    =    _mm_add_epi16(*(__m128i*)&Y[160],    *(__m128i*)&X[160]);
    *(__m128i*)(Y + 168)    =    _mm_add_epi16(*(__m128i*)&Y[168],    *(__m128i*)&X[168]);
    *(__m128i*)(Y + 176)    =    _mm_add_epi16(*(__m128i*)&Y[176],    *(__m128i*)&X[176]);
    *(__m128i*)(Y + 184)    =    _mm_add_epi16(*(__m128i*)&Y[184],    *(__m128i*)&X[184]);
    *(__m128i*)(Y + 192)    =    _mm_add_epi16(*(__m128i*)&Y[192],    *(__m128i*)&X[192]);
    *(__m128i*)(Y + 200)    =    _mm_add_epi16(*(__m128i*)&Y[200],    *(__m128i*)&X[200]);
    *(__m128i*)(Y + 208)    =    _mm_add_epi16(*(__m128i*)&Y[208],    *(__m128i*)&X[208]);
    *(__m128i*)(Y + 216)    =    _mm_add_epi16(*(__m128i*)&Y[216],    *(__m128i*)&X[216]);
    *(__m128i*)(Y + 224)    =    _mm_add_epi16(*(__m128i*)&Y[224],    *(__m128i*)&X[224]);    
    *(__m128i*)(Y + 232)    =    _mm_add_epi16(*(__m128i*)&Y[232],    *(__m128i*)&X[232]);
    *(__m128i*)(Y + 240)    =    _mm_add_epi16(*(__m128i*)&Y[240],    *(__m128i*)&X[240]);
    *(__m128i*)(Y + 248)    =    _mm_add_epi16(*(__m128i*)&Y[248],    *(__m128i*)&X[248]);
}

　　_mm_add_epi16可以一次性完成16個short類型的數據的加法，比傳統的add指令快了很多倍。

     由於_mm_add_epi16是這對短整形數據進行的處理，因此，一般情況下改指令所能處理的半徑不能大於127，如果需要大於127，則需要修改過程序中的short類型為int,同時需要使用_mm_add_epi32指令，這樣程序的速度會有所下降。

　　經過測試，在我的I5的台式機中，1024*768圖像在直方圖更新上所需要的平均之間約為30ms，相比局部算法的核心就算部分時間（比如上述的求最大值），可能大部分耗時並不在這里。

     附件的代碼中有個完整的測試工程，並有我目前所有的TMatrix結構的完整代碼，我以后的文章都將以改結構為依托進行處理。

     代碼還共享了很多處理的函數，我很自信一定值得朋友去學習的。

     這種前后依賴的算法都有一個很致命的缺點，就是不可以並行，把圖像分段處理，也會造成過多初始化耗時。

      代碼下載地址：http://files.cnblogs.com/files/Imageshop/BaseFile.rar

 

****************************作者： laviewpbt   時間： 2015.4.20    聯系QQ:  33184777 轉載請保留本行信息**********************