??????????????????????????????????????????????????? 圖像去暗角

暗角圖像是一種在現實中較為常見的圖像，其主要特征就是在圖像四個角有較為顯著的亮度下降，比如下面兩幅圖。根據其形成的成因，主要有3種：natural vignetting, pixel vignetting, 以及mechanic?vignetting，當然，不管他的成因如何，如果能夠把暗角消除或者局部消除，則就有很好的工程意義。

??

?

　　

? ? ? 這方面的資料和論文也不是很多，我最早是在2014年看到Y. Zheng等人的論文《Single image vignetting correction》以及同樣有他們撰寫的論文《Single image vignetting correction using radial gradient symmetry》有講這方面的算法，不過其實現的復雜度較高，即使能編程實現，速度估計也很慢，其實用性就不高了。

? ? ? 前不久，偶爾的機會看到一篇名為《Single-Image Vignetting Correction by Constrained Minimization of log-Intensity Entropy》的論文，并且在github上找到了相關的一些參考代碼，雖然那個代碼寫的實在是惡心加無聊，但是對于我來說這并不重要，只要稍有參考，在結合論文那自己來實現就不是難事了。

? ? ?論文里的算法核心其實說起來也沒啥難的，我就我的理解來簡單的描述下：

? ? ?第一：去暗角可以說是陰影校正的一種特例，而將整副圖像的熵最小化也被證明為進行陰影校正的一種有效方法，但是普通的熵在優化過程中會優化到局部最優的。因此論文中提出了一種對數熵的概念（Log-Intensity Entropy），論文中用數據做了說明，假設一副普通正常的圖像其直方圖是單峰分布，那么如果這幅圖像有暗角，其直方圖必然會存在另外一個低明度的分布，如下圖所示：

?

　　我們校正暗角的過程就是使低明度的分布向原來的正常明度靠近，由上圖第一行的數據可以看到，普通的熵計算直到兩個直方圖有部分重疊的時候熵才會下降，之前熵一直都是增加的，而對數熵則在沒有重疊前至少是保持不增的，因此能夠更好的獲取全局最優解。

? ? ?那么論文提出的對數熵的計算公式為：

? ? ?首先先將亮度進行對數映射，映射公式為： ? ??

? ? ? ? ? ? ? ? ? ? ??

　　也就是將[0,255]內的像素值映射到[0, N-1]內，但不是線性映射，而是基于對數關系的映射，通常N就是取256，這樣映射后的像素范圍還是[0,255]，但是注意這里的i(L)已經是浮點數了。我們繪制出N等于256時上式的曲線：

? ? ? ? ? ? ? ? ? ??

　　可見，這種操作實際上把圖像整體調亮了。

　　由于映射后的色階已經是浮點數了，因此，直方圖信息的統計就必須換種方式了，論文給出的公式為：

? ? ? ? ? ? ?

? ?

假設對于0-255灰度范圍的像素，進行如下映射，可以把灰度值映射至[0,N-1]

i(L) = (N ? 1) log(1 + L)/ log 256

假設N=256，由于映射之后，圖像灰度值為雙精度，而Hist[K]中表示的K必須是整型數據。對于每個像素的灰度值大小都是雙精度的圖像，為了統計雙精度的灰度直方圖信息，我們必須進行相應的加權處理。、、分別是向下取整，和向上取整，的取值范圍為[0,1]。下圖中橫坐標表示的直方圖的K值，縱坐標表示的是像素數量值，那么向上取整和向下取整兩者之差必為1。而恰好是在[0,1]。假設，，其，那么，對于介于、某一雙精度灰度值而言，我們可以根據其小數部分確定該雙精度值在、處的大小。由三角形相似原理可知：

?

???? 參考論文中，

????? 注意取整之后的大小關系 。向下取整，k小于，向上取整K要大于。

???? 當某一雙精度灰度值位于上下取整數值之間的時候，對向下取整的所對應的直方圖bin像素數量的貢獻量為，對向上取整所對應的直方圖bin像素數量的貢獻量為。相當于把位于上下取整之間的x所插值生成的y大小，按照其權重分別分配給前一個直方圖bin和后一個直方圖bin .

? ? ? 公式很復雜， 其實就是有點類似線性插值那種意思，不認識了那兩個數學符號了，就是向上取整和向下取整。

? ? ? 這樣的對數熵直方圖信息會由于巨大的色階調整，導致很多色階是沒有直方圖信息的，一般可以對這個直方圖信息進行下高斯平滑的，得到新的直方圖。

　　最后圖像的對數熵，計算方法如下：

? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

　　其中：? ??

??????????????????????????????

? ? ? ? ?

?第二：論文根據資料《Radiometric alignment and vignetting calibration》提出了一個暗角圖像亮度下降的關系式，而我去看《Radiometric alignment and vignetting calibration》這篇論文時，他的公式又是從《Vignette and exposure calibration and compensation》中獲取的，所以這個論文的作者寫得文章還不夠嚴謹。這個公式是一個擁有五個未知參數的算式，如下所示：

? ? ? ? ? ? ? ? ???

　　其中：

? ? ? ? ? ? ??

? ? ?其中，x和y是圖像每一點的坐標，而則表示暗角的中心位置，他們和a、b、c均為未知量。

　 ?我們可以看到，當r=0時，校正系數為1，即無需校正。當r=1時，校正系數為1+a+b+c。

? ? ?那么經過暗角校正后的圖像就為：

? ? ? ? ? ? ?

　　按照我們的常識，暗角圖像從暗角的中心點想四周應該是逐漸變暗的，根據上式函數g應該是隨著r單調遞增的（因為我們是校正暗角圖像，所以越靠近邊緣上式的乘法中g值也就應該越大)，因此函數g的一階導數應該大于0，即：

? ? ? ? ? ? ?

? ? ?同時，我們注意到參數r的范圍很明顯應該在[0,1]之間，這樣上式則可以轉換為：

? ? ? ? ? ? ? ??

　　如果令，則上式變為：

? ? ? ? ? ?

　　根據二次不等式相關知識，令：

? ? ? ? ?

? ? ?則論文總結了滿足下述關系式的a,b,c就能滿足上述要求了：

? ? ? ??

　　這個我也沒有去驗證了。

? ? ? 第三: 上面描述了校正暗角圖像的公式（帶參數）以及評價一副圖像是否有暗角的指標，那么最后一步就是用這個指標來確定公式的參數。我們未知的參數有5個，即a、b、c以及暗角的中心點。解這種受限的最優問題是有專門的算法的，也是非常計算耗時的。因此，作者提出了一種快速的算法：Hill climbing with rejection of invalid solutions.

? ? ? Hill climbing with rejection of invalid solutions：簡而言之，我們把a,b,c三個陰影校正參數初始化為0，每個參數獨立的按照一定步長增加并減小，并計算對應的圖像熵。經過六次計算之后，我們獲取到了一組使得圖像熵最小的amin,bmin,cmin參數，然后從這開始再次向之前方法一樣計算，直到找到圖像最小熵對應的三個參數。

? ? ? 首先，很明顯，為了計算這些最優參數，我們沒有必要直接在原圖大小上直接計算，這點在原論文也有說明，我們即使把他們的寬高分別縮小到原圖的1/5甚至1/10計算出來的結果也不會有太大的差異，而這些參數的差異對最終的的結果影響也不大，但是計算量就能減少到原來的1/25和1/100。

? ? ? 計算出上述a、b、c以及中心點后，就可以再次按照校正公式來進行校正了，注意暗角的影響對每個通道都是等同的，因此，每個通道都應該乘以相同的值。

　　下面貼出一些用論文中的算法處理的結果圖：

? ?

??

? ?

? ? ??

代碼實現：

#include <iostream> #include <iomanip> #include <fstream> #include <string> #include <opencv2/core/core.hpp> #include <opencv2/imgproc/imgproc.hpp> #include <opencv2/imgproc/imgproc_c.h> #include <opencv2/highgui/highgui.hpp>using namespace std; using namespace cv; float calH(float a,float b,float c,Mat GrayImg); bool check(float a,float b,float c); int CalculateEntropyFromImage(Mat img,Mat result) {Mat aa = img.clone(); const int rows = img.rows;const int cols = img.cols;Mat resize_img;resize(aa,resize_img,Size(cols/10,rows/10),0,0,INTER_LINEAR);float a=0,b=0,c=0;float a_min=0,b_min=0,c_min=0;float delta=2;float Hmin = calH(a,b,c,img);while(delta>1/256){float a_temp=a+delta;if(check(a_temp,b,c)){float H = calH(a_temp,b,c,resize_img);if(Hmin>H){a_min =a_temp;b_min =b;c_min =c;Hmin = H; }}a_temp=a-delta;if(check(a_temp,b,c)){float H = calH(a_temp,b,c,resize_img);if(Hmin>H){a_min =a_temp;b_min =b;c_min =c;Hmin = H; }}float b_temp=b+delta;if(check(a,b_temp,c)){float H = calH(a,b_temp,c,resize_img);if(Hmin>H){a_min =a;b_min =b_temp;c_min =c;Hmin = H; }}b_temp=b-delta;if(check(a,b_temp,c)){float H = calH(a,b_temp,c,resize_img);if(Hmin>H){a_min =a;b_min =b_temp;c_min =c;Hmin = H; }}float c_temp=c+delta;if(check(a,b,c_temp)){float H = calH(a,b,c_temp,resize_img);if(Hmin>H){a_min =a;b_min =b;c_min =c_temp;Hmin = H; }}c_temp=c-delta;if(check(a,b,c_temp)){float H = calH(a,b,c_temp,resize_img);if(Hmin>H){a_min =a;b_min =b;c_min =c_temp;Hmin = H; }}delta = delta/2.0f;}// cout<<"***************"<<endl;cout<<"amin "<<a_min<<"bmin "<<b_min<<"cmin "<<c_min<<endl;//Mat result(img.size(),CV_8UC1);int c_x = cols/2,c_y = rows/2;const float d = sqrt((float)c_x*c_x+c_y*c_y+1); for(int row=0;row<rows;++row){uchar *data = aa.ptr<uchar>(row);uchar *value = result.ptr<uchar>(row);for(int col=0;col<cols;++col){float r=sqrt((float)((row-c_y)*(row-c_y)+(col-c_x)*(col-c_x)))/d;float r2 = r*r;float r4 = r2*r2;float r6 = r2*r2*r2;float g = 1+a_min*r2+b_min*r4+c_min*r6;// this will cause overflow // ToDo: The image should be normalized to the original brightnessvalue[col] = (data[col]*g);if(value[col]>255.0f)value[col] = 255;if(value[col]<0.0f)value[col]=0;}}convertScaleAbs( result, result); return 0; }int main(int argc, char* argv[]) {if(argc<2)return 0;string path = argv[1];Mat img = imread(path,IMREAD_UNCHANGED);Mat RGB_result(img.size(),CV_8UC3);Mat result = Mat::zeros(img.size(), CV_8UC1);cvtColor(img,img,COLOR_BGR2GRAY); CalculateEntropyFromImage(img,result);imshow("HJ",result);imwrite("HJ.bmp",result);waitKey(0);return 0; }//檢查參數的合法性 bool check(float a,float b,float c) {if ((a>0) && (b==0) && (c==0))return true;if (a>=0 && b>0 && c==0)return true;if (c==0 && b<0 && -a<=2*b)return true;if (c>0 && b*b<3*a*c)return true;if (c>0 && b*b == 3*a*c && b>=0)return true;if (c>0 && b*b == 3*a*c && -b>=3*c)return true;float q_p = (-2*b+sqrt(4*b*b-12*a*c))/(6*c);if (c>0 && b*b > 3*a*c && q_p<=0)return true;float q_d = (-2*b-sqrt(4*b*b-12*a*c))/(6*c);if (c>0 && b*b > 3*a*c && q_d>=1)return true;if (c<0 && b*b > 3*a*c && q_p >=1 && q_d<=0)return true;return false; }//計算圖像熵 float calH(float a,float b,float c,Mat GrayImg) {Mat GrayFloatImg(GrayImg.size(),CV_32FC1);float histogram[256];memset(histogram,0,sizeof(float)*256);const int rows = GrayImg.rows;const int cols = GrayImg.cols; float c_x = cols/2.0,c_y = rows/2.0;const float d = sqrt(c_x*c_x+c_y*c_y+1);for(int row=0;row<rows;++row){uchar *data = GrayImg.ptr<uchar>(row);float *value = GrayFloatImg.ptr<float>(row);for(int col=0;col<cols;++col){float r=sqrt((row-c_y)*(row-c_y)+(col-c_x)*(col-c_x))/d;float r2 = r*r;float r4 = r2*r2;float r6 = r2*r2*r2;float g = 1+a*r2+b*r4+c*r6;value[col] = data[col]*g;if(value[col] >= 255){value[col] = 255.0;histogram[255]++;}else {value[col] = (255.0f)*log(1+value[col])/log(256.0f);int k_d = (int)(floor(value[col]));int k_u = (int)(ceil(value[col]));histogram[k_d] += (1+k_d-value[col]);histogram[k_u] += (k_u-value[col]);}}}float TempHist[256 + 2 * 4]; // SmoothRadius = 4TempHist[0] = histogram[4]; TempHist[1] = histogram[3]; TempHist[2] = histogram[2]; TempHist[3] = histogram[1];TempHist[260] = histogram[254]; TempHist[261] = histogram[253];TempHist[262] = histogram[252]; TempHist[263] = histogram[251];memcpy(TempHist + 4, histogram, 256 * sizeof(float));// smoothfor (int X = 0; X < 256; X++)histogram[X] = (TempHist[X] + 2 * TempHist[X + 1] + 3 * TempHist[X + 2] + 4 * TempHist[X + 3] + 5 * TempHist[X + 4] + 4 * TempHist[X + 5] + 3 * TempHist[X + 6] + 2 * TempHist[X + 7]) + TempHist[X + 8] / 25.0f;float sum =0;for(int i=0;i<256;++i){sum += histogram[i];}float H=0,pk;for(int i=0;i<256;++i){pk = histogram[i]/sum;if(pk!=0)H += pk * log(pk); }return -H; }

?結論：參考相關文獻和代碼，最終計算出來的圖像總是不理想，存在數據溢出的情況。不知道哪出錯了。

參考文獻：

1.圖像增強系列之圖像自動去暗角算法。

2.https://github.com/HJCYFY/Vignetting-Correction? 論文對應的代碼

?

總結

以上是生活随笔為你收集整理的图像去暗角的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。