基礎版本的基礎版本 直方圖均衡化系列

我們最終所需要的算法名字是:Contrast Limited Adaptive Histogram Equalization限制對比度自適應直方圖均衡化,每個單詞取首字母可以縮寫為:CLAHE

直方圖的原理以及作用可以參考這里,總之就是,算最終的均衡化是經歷過計算直方圖以及累積分布直方圖兩個過程,其理論基礎在參考的文獻里解釋的非常清楚.但是從原始的HE到CLAHE中間經歷了多個優化算法,網上很少有直觀的代碼實現,好不容易在這篇文章里找到.下面的介紹主要是結合代碼,進行部分解釋說明.

直方圖均衡化,HE:

(資料圖片僅供參考)

這個是最基礎版本,代碼表現的很清楚

Mat eaualizeHist_GO(Mat src){    int width = src.cols;    int height= src.rows;    Mat HT_GO = src.clone();    int tmp[256] ={0};    float C[256] = {0.0};    int total = width*height;      for (int i=0 ;i(i,j);            tmp[index] ++;        }    }    //計算累積函數      for(int i = 0;i < 256 ; i++){          if(i == 0)              C[i] = 1.0f * tmp[i] / total;          else              C[i] = C[i-1] + 1.0f * tmp[i] / total;      }      //這里的累積函數分配的方法非常直觀高效    for(int i = 0;i < src.rows;i++){          for(int j = 0;j < src.cols;j++){                  int index = src.at(i,j);            HT_GO.at(i,j) = C[index] * 255  ;        }      }      return HT_GO;}

自適應直方圖均衡化,AHE:

把整個大圖分成8X8=64份,當然,這里的舉例圖像的長寬正好是8的倍數.寫入到存儲直方圖的數組的順序是按照沒一行的從左到右, 一共八行.

Mat aheGO(Mat src,int _step = 8){    Mat AHE_GO = src.clone();    int block = _step;    int width = src.cols;    int height = src.rows;    int width_block = width/block; //每個小格子的長和寬    int height_block = height/block;    //存儲各個直方圖      int tmp2[8*8][256] ={0};    float C2[8*8][256] = {0.0};    //分塊    int total = width_block * height_block;     for (int i=0;i<BLOCK;I++){ for="" (int="" j="0;j<block;j++){" int="" index="src.at<uchar" start_x="i*width_block;" width_block;="" start_y="j*height_block;" +="" height_block;="" num="i+block*j;" 遍歷小塊,計算直方圖="" ii="start_x" end_x="start_x" ii++)="" for(int="" jj="start_y" (jj,ii);                    tmp2[num][index]++;                  }              }             //計算累積分布直方圖              for(int k = 0 ; k < 256 ; k++){                  if( k == 0)                      C2[num][k] = 1.0f * tmp2[num][k] / total;                  else                      C2[num][k] = C2[num][k-1] + 1.0f * tmp2[num][k] / total;              }          }    }    //將統計結果寫入    for (int i=0;i<BLOCK;I++){ for="" (int="" j="0;j<block;j++){" int="" index="src.at<uchar" start_x="i*width_block;" width_block;="" start_y="j*height_block;" +="" height_block;="" num="i+block*j;" 遍歷小塊,計算直方圖="" ii="start_x" end_x="start_x" for(int="" jj="start_y" (jj,ii);                    //結果直接寫入AHE_GO中去                    AHE_GO.at(jj,ii) = C2[num][index] * 255  ;                }              }         }    }    return AHE_GO;}

限制對比度直方圖均衡化:CLHE.

對比度指的是一副圖片中最亮的白和最暗的黑之間的反差大小=(max-min)/(max+min).對比度增強可以使用線性拉伸的方式,即,直接將(min1,max1)的范圍擴大到(min2,max2).但是對于比較復雜的圖片沒有什么效果.直方圖均衡化,也是一種拉伸對比度的方式.

HE算法在一種情況下，效果不好，如果一個圖片中有大塊的暗區或者亮區的話，效果非常不好。這個的原因，也非常好理解，因為HE其實要求一個圖片中必須有10%的最亮的像素點，必須有10%第二亮的像素點，必須有10%第三亮的像素點……假設有一張純黑的圖片，你想想經過HE處理之后，會出現什么情況？答案就是一部分黑的像素也會被強行搞成白的.來自.

因此單純的直方圖均衡化,AH比較適合一副圖像整體偏暗或者亮.如果局部有個比較亮或者暗的,就會多出很多干擾信息.因此進行對比度限制,可以減弱帶來的不利影響.

主要體現在:

假如左圖為原圖的直方圖,在計算累積函數的時候,先轉換成右圖.這樣可以達到限制對比度的功能.即進行削峰,并將削下來的量均勻的分配給每個色階(或者亮度值).

Mat clheGO(Mat src,int _step = 8){    int width = src.cols;    int height= src.rows;    Mat CLHE_GO = src.clone();    int tmp[256] ={0};    float C[256] = {0.0};    int total = width*height;      for (int i=0 ;i(i,j);            tmp[index] ++;        }    }    /限制對比度計算部分,注意這個地方average的計算不一定科學    int average = width * height / 255/64;      int LIMIT = 4 * average;      int steal = 0;      for(int k = 0 ; k < 256 ; k++){          if(tmp[k] > LIMIT){              steal += tmp[k] - LIMIT;              tmp[k] = LIMIT;          }      }      int bonus = steal/256;      //hand out the steals averagely      for(int k = 0 ; k < 256 ; k++){          tmp[k] += bonus;      }      ///    //計算累積函數      for(int i = 0;i < 256 ; i++){          if(i == 0)              C[i] = 1.0f * tmp[i] / total;          else              C[i] = C[i-1] + 1.0f * tmp[i] / total;      }      //這里的累積函數分配的方法非常直觀高效    for(int i = 0;i < src.rows;i++){          for(int j = 0;j < src.cols;j++){                  int index = src.at(i,j);            CLHE_GO.at(i,j) = C[index] * 255  ;        }      }      return CLHE_GO;}

CLAHE

將上面的AHE結合就是不帶有插值算法的CLAHE:

Mat claheGoWithoutInterpolation(Mat src, int _step = 8){    Mat CLAHE_GO = src.clone();    int block = _step;//pblock    int width = src.cols;    int height= src.rows;    int width_block = width/block; //每個小格子的長和寬    int height_block = height/block;    //存儲各個直方圖      int tmp2[8*8][256] ={0};    float C2[8*8][256] = {0.0};    //分塊    int total = width_block * height_block;     for (int i=0;i<BLOCK;I++){ for="" (int="" j="0;j<block;j++){" int="" index="src.at<uchar" start_x="i*width_block;" width_block;="" start_y="j*height_block;" +="" height_block;="" num="i+block*j;" 遍歷小塊,計算直方圖="" ii="start_x" end_x="start_x" for(int="" jj="start_y" (jj,ii);                    tmp2[num][index]++;                  }              }             //裁剪和增加操作，也就是clahe中的cl部分            //這里的參數 對應《Gem》上面 fCliplimit  = 4  , uiNrBins  = 255            int average = width_block * height_block / 255;              int LIMIT = 4 * average;              int steal = 0;              for(int k = 0 ; k < 256 ; k++){                  if(tmp2[num][k] >LIMIT){                      steal += tmp2[num][k] - LIMIT;                      tmp2[num][k] = LIMIT;                  }              }              int bonus = steal/256;              //hand out the steals averagely              for(int k = 0 ; k < 256 ; k++){                  tmp2[num][k] += bonus;              }              //計算累積分布直方圖              for(int k = 0 ; k < 256 ; k++){                  if( k == 0)                      C2[num][k] = 1.0f * tmp2[num][k] / total;                  else                      C2[num][k] = C2[num][k-1] + 1.0f * tmp2[num][k] / total;              }          }    }    //計算變換后的像素值      //將統計結果寫入    for (int i=0;i<BLOCK;I++){ for="" (int="" j="0;j<block;j++){" int="" index="src.at<uchar" start_x="i*width_block;" width_block;="" start_y="j*height_block;" +="" height_block;="" num="i+block*j;" 遍歷小塊,計算直方圖="" ii="start_x" end_x="start_x" for(int="" jj="start_y" (jj,ii);                    //結果直接寫入AHE_GO中去                    CLAHE_GO.at(jj,ii) = C2[num][index] * 255  ;                }              }         }         }      return CLAHE_GO;}

這樣的結果會有嚴重的網格感覺,需要使用插值的方式進行解決:

Mat claheGO(Mat src,int _step = 8){    Mat CLAHE_GO = src.clone();    int block = _step;//pblock    int width = src.cols;    int height= src.rows;    int width_block = width/block; //每個小格子的長和寬    int height_block = height/block;    //存儲各個直方圖      int tmp2[8*8][256] ={0};    float C2[8*8][256] = {0.0};    //分塊    int total = width_block * height_block;     for (int i=0;i<BLOCK;I++) for="" (int="" j="0;j<block;j++)" int="" index="src.at<uchar" start_x="i*width_block;" width_block;="" start_y="j*height_block;" +="" height_block;="" num="i+block*j;" 遍歷小塊,計算直方圖="" ii="start_x" end_x="start_x" ii++)="" for(int="" jj="start_y" (jj,ii);                    tmp2[num][index]++;                  }              }             //裁剪和增加操作，也就是clahe中的cl部分            //這里的參數 對應《Gem》上面 fCliplimit  = 4  , uiNrBins  = 255            int average = width_block * height_block / 255;              //關于參數如何選擇，需要進行討論。不同的結果進行討論            //關于全局的時候，這里的這個cl如何算，需要進行討論             int LIMIT = 40 * average;              int steal = 0;              for(int k = 0 ; k < 256 ; k++)              {                  if(tmp2[num][k] > LIMIT){                      steal += tmp2[num][k] - LIMIT;                      tmp2[num][k] = LIMIT;                  }              }              int bonus = steal/256;              //hand out the steals averagely              for(int k = 0 ; k < 256 ; k++)              {                  tmp2[num][k] += bonus;              }              //計算累積分布直方圖              for(int k = 0 ; k < 256 ; k++)              {                  if( k == 0)                      C2[num][k] = 1.0f * tmp2[num][k] / total;                  else                      C2[num][k] = C2[num][k-1] + 1.0f * tmp2[num][k] / total;              }          }    }    //計算變換后的像素值      //根據像素點的位置，選擇不同的計算方法      for(int  i = 0 ; i < width; i++)      {          for(int j = 0 ; j < height; j++)          {              //four coners              if(i <= width_block/2 && j < = height_block/2)              {                  int num = 0;                  CLAHE_GO.at(j,i) = (int)(C2[num][CLAHE_GO.at(j,i)] * 255);              }else if(i <= 2="" j="">= ((block-1)*height_block + height_block/2)){                  int num = block*(block-1);                  CLAHE_GO.at(j,i) = (int)(C2[num][CLAHE_GO.at(j,i)] * 255);              }else if(i >= ((block-1)*width_block+width_block/2) && j <= int="" num="block-1;" else="" i="">= ((block-1)*width_block+width_block/2) && j >= ((block-1)*height_block + height_block/2)){                  int num = block*block-1;                  CLAHE_GO.at(j,i) = (int)(C2[num][CLAHE_GO.at(j,i)] * 255);              }              //four edges except coners              else if( i <= 2="" int="" num_i="0;" num_j="(j" -="" num1="num_j*block" num2="num1" float="" p="(j" q="1-p;" else="" i="">= ((block-1)*width_block+width_block/2)){                  //線性插值                  int num_i = block-1;                  int num_j = (j - height_block/2)/height_block;                  int num1 = num_j*block + num_i;                  int num2 = num1 + block;                  float p =  (j - (num_j*height_block+height_block/2))/(1.0f*height_block);                  float q = 1-p;                  CLAHE_GO.at(j,i) = (int)((q*C2[num1][CLAHE_GO.at(j,i)]+ p*C2[num2][CLAHE_GO.at(j,i)])* 255);              }else if( j <= 2="" int="" num_i="(i" -="" num_j="0;" num1="num_j*block" num2="num1" float="" p="(i" q="1-p;" else="" j="">= ((block-1)*height_block + height_block/2) ){                  //線性插值                  int num_i = (i - width_block/2)/width_block;                  int num_j = block-1;                  int num1 = num_j*block + num_i;                  int num2 = num1 + 1;                  float p =  (i - (num_i*width_block+width_block/2))/(1.0f*width_block);                  float q = 1-p;                  CLAHE_GO.at(j,i) = (int)((q*C2[num1][CLAHE_GO.at(j,i)]+ p*C2[num2][CLAHE_GO.at(j,i)])* 255);              }              //雙線性插值            else{                  int num_i = (i - width_block/2)/width_block;                  int num_j = (j - height_block/2)/height_block;                  int num1 = num_j*block + num_i;                  int num2 = num1 + 1;                  int num3 = num1 + block;                  int num4 = num2 + block;                  float u = (i - (num_i*width_block+width_block/2))/(1.0f*width_block);                  float v = (j - (num_j*height_block+height_block/2))/(1.0f*height_block);                  CLAHE_GO.at(j,i) = (int)((u*v*C2[num4][CLAHE_GO.at(j,i)] +                       (1-v)*(1-u)*C2[num1][CLAHE_GO.at(j,i)] +                      u*(1-v)*C2[num2][CLAHE_GO.at(j,i)] +                      v*(1-u)*C2[num3][CLAHE_GO.at(j,i)]) * 255);              }              //最后這步，類似高斯平滑            CLAHE_GO.at(j,i) = CLAHE_GO.at(j,i) + (CLAHE_GO.at(j,i) << 8) + (CLAHE_GO.at(j,i) << 16);                 }      }    return CLAHE_GO;}

插值的方式就是根據一個點周圍四個點的值來確定.只是說其值為累積分布直方圖的值.分成64個塊.每個塊又分成四個小塊.整個圖像的四個頂點所在的小塊不用插值.除此之外的四個邊采用"單"線性插值.剩下的為雙線性插值.因此可以簡單的理解為,只有相鄰的小塊才會進行插值.

整體來講上面遺留了兩個問題:

int average = width_block * height_block / 255;  int LIMIT = 40 * average;  int steal = 0;

1、在進行CLAHE中CL的計算，也就是限制對比度的計算的時候，參數的選擇缺乏依據。在原始的《GEMS》中提供的參數中， fCliplimit  = 4  , uiNrBins  = 255.但是在OpenCV的默認參數中，這里是40.就本例而言，如果從結果上反推，我看10比較好。這里參數的選擇缺乏依據；

2、CLHE是可以用來進行全局直方圖增強的，那么這個時候，這個average 如何計算，肯定不是width * height/255,這樣就太大了，算出來的LIMIT根本沒有辦法獲得。

優化

該博主的圖像處理系列很值得細細品味。這里也結合他對直方圖均衡化系列的文章，提取出比較新穎的觀點進行總結。

一：三通道聯合處理：

for (Y = 0; Y < Height; Y++){        Pointer = Scan0 + Y * Stride;               // 定位到每個掃描行的第一個像素，以避免溶于數據的影響        for (X = 0; X < Width; X++){            HistGram[*Pointer]++;                   // Blue            HistGram[*(Pointer + 1)]++;             // Green            HistGram[*(Pointer + 2)]++;             // Red                Pointer += 3;                           // 移向下一個像素        }    }    Num = 0;    for (Y = 0; Y < 256; Y++){        Num = Num + HistGram[Y];        Lut[Y] = (byte)((float)Num / (Width * Height * 3) * 255);       // 計算映射表    }    for (Y = 0; Y < Height; Y++){        Pointer = Scan0 + Y * Stride;        for (X = 0; X < Width * 3; X += 3){            Pointer[X] = Lut[Pointer[X]];            Pointer[X + 1] = Lut[Pointer[X + 1]];            Pointer[X + 2] = Lut[Pointer[X + 2]];        }    }

二：強化的基于局部直方圖裁剪均衡化的對比度調節算法

比如說該篇中除了提到了局部直方圖和全局直方圖以及亮度分量的融合方法，也對累積分布直方圖中（局部均衡化后映射表）的平滑思想做了介紹：

第一種思想就是，將色階（bins，比如256）均勻分成K（比如說是32）份，每份的開始的值和索引（原始像素值）構成二維序列點。這樣根據這K個點可以擬合出一條平滑曲線。這個曲線插值成離散的后即可成為新的映射表。這樣的累積分布直方圖就會更加的平滑。

第二種思想就是，將映射表（累計分布直方圖）中的值，進行一維的均值或者高斯濾波，同樣也可以達到平滑的作用。

三：自動色階

這篇是基于自動色階的算法來實現的圖像增強，首先，自動色階是另外一種裁剪直方圖的方式，通過設置lowcut和highcut兩個參數來確定裁剪直方圖兩頭的程度，裁剪的兩端設置為極端值（比如說0或者1），中間的在重新映射到0-1（0-255）的范圍，使像素值進行拉伸（可以線性,即拉伸時的系數為1，也可以gamma曲線的方式），從而達到對比度增強的效果。

# 裁剪PixelAmount = Width * Height                     "所有像素的數目      Sum = 0      For Y = 0 To 255          Sum = Sum + HistBlue(Y)          If Sum >= PixelAmount * LowCut * 0.01 Then   "注意PS界面里的那個百分號              MinBlue = Y                              "得到藍色分量的下限              Exit For                                 "退出循環          End If      Next     Sum = 0     For Y = 255 To 0 Step -1         Sum = Sum + HistBlue(Y)         If Sum >= PixelAmount * HighCut * 0.01 Then  "注意PS界面里的那個百分號             MaxBlue = Y                              "得到藍色分量的上限             Exit For                                 "退出循環         End If     Next# 映射     For Y = 0 To 255         If Y <= 0="" minblue="" then="" elseif="" y="">= MaxBlue Then             BlueMap(Y) = 255         Else             BlueMap(Y) = (Y - MinBlue) / (MaxBlue - MinBlue) * 255      "線性映射         End If     Next

自動對比度與自動色階稍有不同的地方是自動色階各通道（對于多通道，如果是灰度圖這種單通道，兩者算法一樣），動對比度算法首先獲取三個通道下限值的最小值，以及上限值的最大值，然后以此為新的上下限，計算映射表。

上面提到，在拉伸過程中，也可以使用gamma矯正的方式，即在這里：

((Y - Min) / (Max - Min)) * 255

線性拉伸為1,gamma的方式可以為：

pow((float)(Y - Min) / (Max - Min), Gamma) * 255

gamma的值可以根據如下求得：

float Avg = 0, Mean = 0, Sum = 0;    for (int Y = 0; Y < 256; Y++)    {        Sum += Histgram[Y];        Avg += Y * Histgram[Y];    }    Mean = Avg / Sum;    float Gamma = log(0.5f) / log((float)(Mean - Min) / (Max - Min));    if (Gamma < 0.1f)        Gamma = 0.1f;    else if (Gamma > 10)        Gamma = 10;

局部自適應自動色階除了上面介紹的自動色階方法還結合了CLAHE，當然這里的自適應限制對比度直方圖均衡話中的限制就不需要了，因為限制也是一種裁剪方式。在計算最后的裁剪位置后，還可以通過設置一個參數來進行調整對比度的程度。

void MakeMapping(int* Histgram,float CutLimit=0.01,float Contrast = 1){    int I, Sum = 0, Amount = 0;    const int Level = 256;    for (I = 0; I < Level; I++) Amount += Histgram[I];    int MinB =0 ,MaxB=255;    int Min = 0,Max=255;    for (I = 0; I < Level; I++){        if  (Histgram[I]!=0){            Min = I ;            break;}    }    for(I = Level-1; I >= 0; I--){        if  (Histgram[I]!=0){            Max = I ;            break;}    }    for (I = 0; I < Level; I++){        Sum = Sum + Histgram[I];        if (Sum >= Amount * CutLimit){            MinB = I;                                          break;}    }      Sum = 0;    for(I = Level-1; I >= 0; I--){        Sum = Sum +Histgram[I];        if (Sum >= Amount * CutLimit ){            MaxB = I ;                                        break;}       }    int Delta = (Max - Min) * Contrast * 0.5  ;    Min = Min - Delta;    Max = Max +    Delta ;    if (Min    < 0) Min = 0;    if (Max > 255) Max = 255;    if (MaxB!=MinB){        for (I = 0; I < Level; I++){            if (IMaxB) Histgram[I]=Max;            else                Histgram[I] = (Max-Min)* (I - MinB) / (MaxB - MinB) + Min ;             }    }    else{        for (I = 0; I < Level; I++) Histgram[I]=MaxB;        //     必須有，不然會有一些圖像平坦的部位效果出錯    }}

圖像增強

整理上面提到博主的一些其他圖像增強的方法。

這里提到的增強平時常用的一種銳化方法，過程即是將圖像分為高低頻，然后高頻部分的乘以一個大于1的系數。這樣得到的原圖的邊緣部分就會放大，達到了銳化或對比度增強的效果。

系數的選擇可以使用動態的方法，這樣可以避免某些地方過大的銳化造成的震鈴現象。D可以使用全局平均值或者全局均方差。

這里提到的是（多尺度視網膜增強算法）MSRCR的色彩增強算法。主要是根據下面的公式來：

Log[R(x,y)] = Log[I(x,y)]-Log[L(x,y)]

I是原始圖像，R是增強后的圖像，L為原圖經過高斯或者均值模糊后的圖片。后面提到的尺度也就是模糊核的半徑。

通過上式反推出R，即增強后的圖像。文章中提到的量化算法比較新穎。本來求解中涉及到的log函數需要通過exp才能的出R，但是可以通過求解log[R(x,y)]（value）的最大（max）最?。╩in）值后，線性方法的方式得到量化結果。不知道為什么這樣做，難道只是為了提速？

R(x,y) = ( Value - Min ) / (Max - Min) * (255-0)

所謂的多尺度就是進行多個尺度的模糊，然后進行權重求和。

Log[R(x,y)] =  Log[R(x,y)] + Weight(i)* ( Log[Ii(x,y)]-Log[Li(x,y)])

其中Weight(i)表示每個尺度對應的權重，要求各尺度權重之和必須為1，經典的取值為等權重。

但是，SSR（單尺度）和MSR（多尺度）在最大尺度相同的時候效果誰好誰壞不太好說。

這種方式會導致色差，解決的方式一般可以通過一下過程：

（1）分別計算出 Log[R(x,y)]中R/G/B各通道數據的均值Mean和均方差Var（注意是均方差）。   （2）利用類似下述公式計算各通道的Min和Max值。            Min = Mean - Dynamic * Var;              Max = Mean + Dynamic * Var;    (3)  對Log[R(x,y)]的每一個值Value，進行線性映射：            R(x,y) = ( Value - Min ) / (Max - Min) * (255-0) ，同時要注意增加一個溢出判斷,即：           if (R(x,y) > 255)  R(x,y) =255; else if (R(x,y) < 0) R(x,y)=0;

Dynamic取值越小，圖像的對比度月想，一般為2-3能取得比較好的效果。

這里是一種gamma矯正的方式進習慣你對比度增強的算法。主要是進行一個動態的gamma值。使得小于128的像素使用gamma值為0-1的范圍，從而使原來的元素值變大。反之，亦然。并且距離中心128越遠，gamma值的變化會越劇烈。

其中，BFmask是雙邊濾波或者均值濾波模糊后的圖像。 α一般取2。也可以動態矯正，比如說，對于低對比度的圖像，應該需要較強烈的校正，因此α值應該偏大，而對于有較好對比度的圖，α值應該偏向于1，從而產生很少的校正量。

對于三種通道RGB的的處理方式，除了三通道分別處理外，也可以使用如下方式：

一：原圖轉換到YUV或者HSV這中帶亮度的顏色空間中，然后用新得到的luminance值代替Y通道或V通道，然后在轉換會RGB空間

二：新的luminance值和原始luminance值的比值作為三通到的增強系數，這樣三通道可以得到同樣程度的增強。

三：可以根據該公式：

第一種方法容易出現結果圖色彩偏淡，第二種每個分量易出現過飽和，第三種可能要稍微好一點，建議使用第三種。