一種快速簡便優秀的全域性曲線調整與區域性資訊想結合的非線性彩色增強演算法（多圖深度分析和探索）
本文的分析基於《Adaptive and integrated neighborhood-dependent approach for nonlinear enhancement of color images》一文相關內容，但對其進行了深度的改良。
我們首先解讀或者說翻譯下這篇論文。
論文公佈的時間是2005年了，已經是比較久遠的了，我第一次看到該論文大概是在2014年，後面在2016年左右實現了該演算法，這裡還有當時開發留下的記錄，居然是除夕左右做的。佩服自己。
以前沒有特別注意到該演算法的效果，覺得也就那樣吧，所以沒有怎麼去發揮它，但是，最近再次審視，發現他除了實現實現簡單、速度快，而且還具有效果佳、適應性廣、不破壞本身就光照好的位置等等眾多優點，似乎比目前我看到的低照度增強演算法都好。
演算法本身的步驟分為三步，第一步是根據的影象亮度分佈建立一個自適應的全域性對映函式，這一步大大的提高了影象中暗部畫素的畫素值，同時也壓縮了影象的動態範圍。第二步是所謂的自適應對比度增強，根據畫素領域內的平均值和畫素值本身的比例，做一個對映，提高整體的對比度。後續還有一個步驟是顏色恢復的過程。
第一步的全域性曲線調整如下所示：
首先計算出彩色影象的亮度值，這個其實可以有很多方式，包括常用的YUV空間的Y通道，HSL空間的L分量，甚至可以使用我提到的對比度保留去色那種方式獲取，論文裡使用的是最普通的計算公式：
我們把亮度歸一化得到：
然後我們用下面這個傳輸函式（對映函式）來將亮度影象進行一個線性的增強：
其中的引數Z值由影象本身的內容決定，如下所示：
式中的L表示亮度影象的累計直方圖（CDF）達到0.1時的色階值，也就是說如果亮度影象中90%的畫素值都大於150，則Z=1，如果10%或者更多的畫素值都小於50，則Z取值為0，否則其他情況Z則根據L的值線性插值。
稍作分析下，如果Z=0，說明影象中存在大量的偏暗畫素，影象有必要變亮一些，如果Z=1，則說明影象已經很亮了，則此時影象無需繼續加亮處理。介於兩者之間時，我們也就做中和處理。
那麼我們採用的傳輸函式是否能達到這種需求呢，我們來看下公式（3）的組成和曲線。
以Z=0為例，如上圖2所示，曲線6為公式（3）最後對應的結果。我們看看其組成。公式（3）的第一和第二部分對應圖2中的曲線2和3，他們的和得到曲線4。曲線2中在低亮度區域的增強很顯著，在高亮度區域逐漸變緩，曲線3是一個線性函式，隨亮度增加線性的減少。公式（3）的第三部分對應曲線5，曲線4和曲線5累加後得到曲線6。可見最終的曲線在低亮度區域顯著增強亮度，在高亮度區域緩慢增加亮度。當然這裡也是可以採用其他具有類似作用的曲線的。
對於不同的Z值，影象給出了最終的相應曲線，可見，隨著Z值的增加，曲線逐漸變為一條直線（不變化），而我們前面的自適應Z值計算的過程也隨著影象亮度的增加而增加，也就是說如果圖形原本就很亮，則Z值就越大，基本上影象就沒什麼調整了，這基本上就是自適應了。
第二步：自適應對比度增強。
經過第一步處理後，影象的亮度自適應的得到了增強，但是影象的對比度明顯減少了。普通的全域性對比度增強演算法的過程是使亮的畫素更亮，暗的畫素更暗(似乎和第一步有點相反的感覺），這樣影象的動態方範圍就更廣了，同時由於這種方法在處理不考慮領域的資訊，對於那些領域和他只有細微的差異的畫素，其細節很難得到有效提升。因此，我們需要考慮區域性的對比度增強，這種增強下，不同位置相同畫素值得畫素在增強後可以得到不同的結果，因為他們一般會有不同的領域畫素值。噹噹前畫素值比周邊畫素的平均值大時，我們增大當前畫素值，而當前畫素值比周邊畫素平均值小時，我們減少它的值。這樣的處理過程結果值和當前畫素值的絕對值沒有關係，至於領域資訊有關了。這樣，影象的的對比度和細節都能得到有效的提升，同時影象的動態範圍也有得到有效的壓縮。
通常，一個比較好的領域計算方式是高斯模糊，我們採用的計算公式如下所示：
式中，S(x,y)為對比度增強的結果。指數E(x,y)如下所示：
下標Conv表示卷積，注意這裡是對原始的亮度資料進行卷積。P是一個和影象有關的引數，如果原始影象的對比度比較差，P應該是一個較大的值，來提高影象的整體對比度，我們通過求原始亮度圖的全域性均方差來決定P值的大小。
當全域性均方差小於3時（說明影象大部分地方基本是同一個顏色了，對比度很差），此時P值取大值，當均方差大於10時，說明原圖的對比度還是可以的，減少增強的程度，均方差介於3和10之間則適當線性增強。
這個對比度增強的過程分析如下，
當卷積值小於原始值時，也就是說中心點的亮度大於周邊的亮度，此時E（x,y）必然小於1，由於I’（x,y）在前面已經歸一化，他是小於或等於1的，此時式（7）
的值必然大於原始亮度值，也就是亮的更亮（這裡的亮不是說全域性的亮，而是區域性的亮）。如果卷積值大於原始值，說明中心點的亮度比周邊的暗，此時E(x,y)大於1，導致處式（7）處理後的結果值更暗。
為了得到更好的對比度增強效果，我們一般都使用多尺度的卷積和增強，因為各個不同的尺度能帶來不同的全域性資訊。一般來說，尺度較小時，能提高區域性的對比度，但是可能整體看起來不是很協調，尺度較大時，能獲得整體影象的更多資訊，但是細節增強的力度稍差。中等程度的尺度在細節和協調方面做了協調。所以，一般類似於MSRCR，我們用不同尺度的資料來混合得到更為合理的結果。
尺度的大小，我們可以設定為固定值，比如5,20,120等，也可以根據影象的大小進行一個自適應的調整。
第三步：顏色恢復。
此步比較簡單，一般就是使用下式：
Lamda通常取值就為1，這樣可以保證影象整體沒有色彩偏移。
以上就是論文的主要步驟，按照這個步驟去寫程式碼也不是一件非常困難的事情：
int IM_BacklightRepair(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride)
{
int Channel = Stride / Width;
if (Channel != 3) return IM_STATUS_INVALIDPARAMETER;
if ((Src == NULL) || (Dest == NULL)) return IM_STATUS_NULLREFRENCE;
if ((Width <= 0) || (Height <= 0)) return IM_STATUS_INVALIDPARAMETER;
int Status = IM_STATUS_OK;
int RadiusS = 5, RadiusM = 20, RadiusL = 120;
const int LowLevel = 50, HighLevel = 150;
const float MinCDF = 0.1f;
int *Histgram = (int *)calloc(256, sizeof(int));
unsigned char *Table = (unsigned char *)malloc(256 * 256 * sizeof(unsigned char)); // 各尺度的模糊
unsigned char *BlurS = (unsigned char *)malloc(Width * Height * sizeof(unsigned char)); // 各尺度的模糊
unsigned char *BlurM = (unsigned char *)malloc(Width * Height * sizeof(unsigned char));
unsigned char *BlurL = (unsigned char *)malloc(Width * Height * sizeof(unsigned char));
unsigned char *Luminance = (unsigned char *)malloc(Width * Height * sizeof(unsigned char));
if ((Histgram == NULL) || (Table == NULL) || (BlurS == NULL) || (BlurM == NULL) || (BlurL == NULL) || (Luminance == NULL))
{
Status = IM_STATUS_OUTOFMEMORY;
goto FreeMemory;
}
float Z = 0, P = 0;
Status = IM_GetLuminance(Src, Luminance, Width, Height, Stride, false); // 得到亮度分量
if (Status != IM_STATUS_OK) goto FreeMemory;
for (int Y = 0; Y < Height * Width; Y++) Histgram[Luminance[Y]]++; // 統計亮度分量的直方圖
float Sum = 0, Mean = 0, StdDev = 0;
for (int Y = 0; Y < 256; Y++) Sum += Histgram[Y] * Y; // 畫素的總和，注意用float型別儲存
Mean = Sum / (Width * Height); // 平均值
for (int Y = 0; Y < 256; Y++) StdDev += Histgram[Y] * (Y - Mean) * (Y - Mean);
StdDev = sqrtf(StdDev / (Width * Height)); // 全域性影象的均方差
int CDF = 0, L = 0;
for (L = 0; L < 256; L++)
{
CDF += Histgram[L];
if (CDF >= Width * Height * MinCDF) break; // where L is the intensity level corresponding to a cumulative distribution function CDF of 0.1.
}
if (L <= LowLevel)
Z = 0;
else if (L <= HighLevel)
Z = (L - LowLevel) * 1.0f / (HighLevel - LowLevel); // 計算Z值
else
Z = 1;
if (StdDev <= 3) // 計算P值，Also, P is determined by the globaln standard deviation of the input intensity image Ix, y as
P = 3;
else if (StdDev <= 10)
P = (27 - 2 * StdDev) / 7.0f;
else
P = 1;
for (int Y = 0; Y < 256; Y++) // Y表示的是I的卷積值
{
for (int X = 0; X < 256; X++) // X表示的I(原始亮度值)
{
float I = X * IM_INV255; // 公式2
I = (powf(I, 0.75f * Z + 0.25f) + (1 - I) * 0.4f * (1 - Z) + powf(I, 2 - Z)) * 0.5f; // 公式3
Table[Y * 256 + X] = IM_ClampToByte(255 * powf(I, powf((Y + 1.0f) / (X + 1.0f), P)) + 0.5f); // 公式7及8
}
}
Status = IM_GaussBlur(Luminance, BlurS, Width, Height, Width, RadiusS);
if (Status != IM_STATUS_OK) goto FreeMemory;
Status = IM_GaussBlur(Luminance, BlurM, Width, Height, Width, RadiusM);
if (Status != IM_STATUS_OK) goto FreeMemory;
Status = IM_GaussBlur(Luminance, BlurL, Width, Height, Width, RadiusL);
if (Status != IM_STATUS_OK) goto FreeMemory;
for (int Y = 0; Y < Height; Y++)
{
unsigned char *LinePS = Src + Y * Stride;
unsigned char *LinePD = Dest + Y * Stride;
int Index = Y * Width;
for (int X = 0; X < Width; X++, Index++, LinePS += 3, LinePD += 3)
{
int L = Luminance[Index];
if (L == 0)
{
LinePD[0] = 0; LinePD[1] = 0; LinePD[2] = 0;
}
else
{
int Value = ((Table[L + (BlurS[Index] << 8)] + Table[L + (BlurM[Index] << 8)] + Table[L + (BlurL[Index] << 8)]) / 3); // 公式13
LinePD[0] = IM_ClampToByte(LinePS[0] * Value / L);
LinePD[1] = IM_ClampToByte(LinePS[1] * Value / L);
LinePD[2] = IM_ClampToByte(LinePS[2] * Value / L);
}
}
}
FreeMemory:
if (Histgram != NULL) free(Histgram);
if (Table != NULL) free(Table);
if (BlurS != NULL) free(BlurS);
if (BlurM != NULL) free(BlurM);
if (BlurL != NULL) free(BlurL);
if (Luminance != NULL) free(Luminance);
}
為了提高速度，我們構建了一個2維（256*256）大小的查詢表，這也是一種很常規的加速方法。
演算法的程式碼部分似乎需要的解釋的部分不多，都是一些常規的處理，也基本就是一步一步按照流程來書寫的。我們來看演算法的效果。
由這兩幅圖的結果我們初步得到這樣的結論：一是該演算法很好的保護了原本對比度和亮度就非常不錯的部分（比如兩幅圖的天空部分基本上沒有什麼變化），這比一些其他的基於Log空間的演算法，比如本人部落格裡的MSRCR，全域性Gamma校正等演算法要好很多，那些演算法處理後原本細節很不錯的地方會發生較大的變化。這不利於影象的整體和諧性。第二，就是暗部的增強效果確實不錯，很多細節和紋理都能更為清晰的表達出來。
為了更為清晰的表達出步驟1和步驟2的處理能力，我們把亮度圖、單獨的步驟1的結果圖和單獨的步驟的結果圖，以及步驟1和步驟2在一起的結果圖比較如下：
亮度圖 全域性曲線調整（步驟1）
區域性對比度增強（步驟2） 步驟1和步驟2綜合
至於如何得到這些中間結果，我想看看程式碼稍作修改就應該沒有什麼大問題吧。
可以看到，步驟1的結果圖中有一部分不是很和諧，有塊狀出現，這個在後續的步驟我們會提及如何處理。
其實我把這四個圖放在一起，我想說的就是經過這麼久的閱讀論文，我覺得所有的這類演算法都是應該是這種框架，全域性亮度調整+區域性對比度調節。不同的演算法只是體現在這兩個步驟使用不同的過程，而想MSRCR這類演算法則把他們在一個過程內同時實現，這樣可能還是沒有本演算法這樣靈活。
經過嘗試，這裡的第一步全域性亮度調整如果使用自動色階或者直方圖均衡化後得到的結果並不是很理想，至於為什麼，暫時沒有仔細的去想。
接著我們討論下這個演算法的一些問題及其改進方法。
問題一：問題我們注意到在上面的圖中全域性亮度調整後的圖中的一些明顯的視覺瑕疵在經過和區域性對比度增強混合後在最終的合成圖中似乎表現得並不是那麼誇張，但是這並不表明這個問題可以忽視，我們看一下下面這張圖的結果：
雖然原圖的亮度比較低，但是在視覺上原圖的可接受程度要比處理後的圖更為好，這主要是因為處理後的圖在暗處顯示出了很多色塊和色斑，而這些色斑在原圖中是無法直接看到的，經過增強後他們變得非常的突兀，也就是說他們增強的程度過於強烈，這個原因是核心在於步驟的1的全域性調整，在圖2中我們看到低亮度處的調整曲線十分的陡峭，這也就意味著增強的程度特別高，會出現的一個現象就是哪怕原始的畫素只差1個值，在處理後的結果中會相差幾十以上，視覺中表現為色塊的現象。我們從上述圖的（39，335）座標處取以20*20的放塊來觀察：
一種簡單的處理方式就是對放大的幅度進行限制，比如一般我們認為，前後處理的結果不應該超過4倍或者其他值，也可以根據影象的內容去自適應設定這個值，這樣就能有效的避免原文出現這樣的問題，修正的程式碼如下所示：
int Value = ((Table[L + (BlurS[Index] << 8)] + Table[L + (BlurM[Index] << 8)] + Table[L + (BlurL[Index] << 8)]) / 3); // 公式13 float Cof = IM_Min(Value * 1.0f / L, 4); LinePD[0] = IM_ClampToByte(LinePS[0] * Cof); LinePD[1] = IM_ClampToByte(LinePS[1] * Cof); LinePD[2] = IM_ClampToByte(LinePS[2] * Cof);
當Cof上限位不同值時，效果也有所區別，如下面兩圖所示：
Cof上限為4 Cof上限為8
根據個人的經驗，Cof設定為4基本上能在增強效果和瑕疵之間達到一個平衡。
問題二：邊緣問題，我們來看下面兩幅測試圖及其效果：
原圖 演算法一次處理 演算法二次處理後
原圖 演算法二次處理後
我們注意到，對於這兩幅圖，大部分的增強效果都是非常不錯的，特別是經過二次增強後，演算法的細節和飽和度等都比較不錯，但是注意到在邊緣處，比如小孩的帽子處、聖誕樹和天空的邊緣處等等明顯發黑，也就是說他沒有得到增強。這是怎麼回事呢，我們以第一幅圖為例，我們檢視下他的亮度圖單獨經過步驟1和步驟2後處理結果：
亮度圖 全域性增強圖 區域性對比度增強圖
可以看出，全域性增強圖在邊緣處未發現有任何問題，而區域性對比度圖在邊緣處變得特別黑，我們將亮度圖減去區域性對比度增強後的圖得到下圖：
在頭髮邊緣我們看到了明顯的白邊。我們分析下，在頭髮邊緣處，畫素比較暗，進行卷積時，周邊是天空或者窗戶等亮區域，這樣進行卷積時，無論卷積的半徑大小是多少，得到的結果必然是平均值大於中心畫素值（而且偏離的比較遠），根據前述的對比度增強的原則，這個時候頭髮處就應該變得更黑了，而在遠離邊緣區，卷積的值不會和中心畫素值有如此大的差異，應該對比度增強的程度也不會如此誇張，就出現上述最終的結果，我們在本文第一個貼圖的地面影子增強處也用紅色線框標註出了連這種現象。
此現象在很多具有強邊緣的影象中出現的比較明顯，而對於普通的自然照片一般難以發現，在論文作者提供的素材中似乎未有該現象發生。
解決這個問題的方案很明顯，需要使用那些能夠保證邊緣不受濾波器影響的或影響的比較少的卷積演算法，這當然就是邊緣保留濾波器的範疇了，雖然現在邊緣濾波器有很多種型別，但是最為廣泛的還是雙邊濾波器、導向濾波等等。我們這裡使用導向濾波來試驗下是否能對結果有所改進。
在導向濾波中，導向的半徑和Eps是影響濾波器最為核心的兩個引數，當Eps固定時，半徑很小時，影象有一種毛絨絨的感覺，稍大一點半徑，則影象能顯示出較好的保邊效果，在非邊緣區則出現模糊效果，而當半徑進一步增大時，整個影象的變換比較小，整體有一種淡淡的朦膿感。當半徑固定時，Eps較小時，影象較為清晰，隨著Eps增加，整個影象就越模糊，到一定程度就和同半徑模糊沒有什麼區別了。經過摸索和多次測試，個人認為半徑引數配置為IM_Min(Width, Height) / 50，Eps引數為20時，能取得非常不錯的保邊效果。此時，我們將前面的三次模糊直接用一個保邊代替，能得到下面的處理效果：
可以明顯的看出邊緣部分得到了完美的解決，無任何瑕疵出現了。
使用單個保邊濾波代替多尺度的高斯模糊，我偶然在測試一幅圖中又發現了另外一個問題，如下所示（只是從原圖中截取了部分顯示）。
原圖區域性 使用單個導向濾波後處理的結果
看到處理後的圖，感覺到非常的失望，這個是怎麼回事呢，後面我單獨測試這個圖後面亮度圖對應的導向濾波的結果，發現也是帶有明顯的紋路感覺的結果。而且嘗試了雙邊、表面模糊等其他EPF濾波器也得到了同樣的結果，但是缺意外的發現作者原始的使用多尺度的高斯模糊的結果卻相當的不錯:
原始多尺度的處理結果 高斯單尺度處理結果
於是我們又重新做了個試驗，把原始的多尺度也改成單尺度的，並且尺度大小和導向濾波用的相同，得到的結果如上圖所示。
這種分析表明對於這個影象並不是因為我用了導向濾波才導致這個紋路出現的，使用高斯濾波，在單尺度時也一樣有問題。但為什麼多尺度時這個問題就消失了呢。
針對這一現象，我做了幾個方面的分析，第一，確實存在這一類影象，在正常時我們看不出這種塊狀，但是當進行模糊或卷積時，這種塊狀就相當的明顯了，哪怕最為平滑的高斯模糊也會出現明顯的分塊現象，下面是對上述區域性原圖分別進行了半徑10和20的高斯並適當加亮（以便顯示）:
所以單尺度的高斯模糊處理後的結果也必然會帶有塊狀，但是當我們使用多尺度時，我們注意到尺度不同時這個色塊的邊界是不同的，而我們多個尺度之間時相互求平均值，此時原本在一個尺度上分界很明顯的邊界線就變得較為模糊了，下圖是上面兩個圖求平均後的結果：
相對來說色塊邊界要比前面兩個圖減弱了不少。
那麼對於使用EPF濾波器的過程，我們如果也使用多尺度的方法，能否對結果進行改善呢，我們這樣處理，也採用三個尺度的保邊濾波，一個比一個大，但是保持Eps的值不變，經過測試，得到的結果如上所示，雖然仔細看還是能看出色塊的存在，但是比原來的要好了很多。
最後，我們還是來簡單的談下演算法的優化問題。
第一個可優化的地方是2維查詢表的建立過程，開始以為只有65536個元素的計算，所以查詢表順序是沒有怎麼仔細考慮的，但是實測，這一塊佔用的時間還是蠻可觀的，有好幾毫秒，主要是因為這裡的powf是個很耗時的過程，所以我們主要稍微把迴圈的位置調換一下，就可以減少大量的計算了，如下:
for (int Y = 0; Y < 256; Y++) // Y表示的I(原始亮度值)
{
float InvY = 1.0f / (Y + 0.01f);
float I = Y * IM_INV255; // 公式2
I = (powf(I, 0.75f * Z + 0.25f) + (1 - I) * 0.5f * (1 - Z) + powf(I, 2 - Z)) * 0.5f; // 公式3
for (int X = 0; X < 256; X++) // X表示的是I的卷積值
{
Table[Y * 256 + X] = IM_ClampToByte(255 * powf(I, powf((X + 0.01f) * InvY, P)) + 0.5f); // 公式7及8
}
}
然後耗時的部分就是LinePD[0] = IM_ClampToByte(LinePS[0] * Cof);這樣的程式碼了，這個也可以通過定點化處理，並配合SSE優化做處理。
優化完成後的程式處理1080P的24點陣圖像大概需要50ms。
最後我們分享幾組利用該演算法處理影象的結果。
當然，有些圖在本演算法處理後還可以加上自動色階或直方圖均衡等操作進一步提升影象的質量。
至此，關於低照度影象的增強演算法，我想我應該不會再怎麼去碰他了，也該休息了。
本文演算法的測試例程見 : http://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，位於選單Enhace->BackLightRepair選單中。
&n