本文的分析基於《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菜單中。
