图像去暗角
??????????????????????????????????????????????????? 圖像去暗角
暗角圖像是一種在現(xiàn)實(shí)中較為常見的圖像，其主要特征就是在圖像四個角有較為顯著的亮度下降，比如下面兩幅圖。根據(jù)其形成的成因，主要有3種：natural vignetting, pixel vignetting, 以及mechanic?vignetting，當(dāng)然，不管他的成因如何，如果能夠把暗角消除或者局部消除，則就有很好的工程意義。
??
?
? ? ? 這方面的資料和論文也不是很多，我最早是在2014年看到Y(jié). Zheng等人的論文《Single image vignetting correction》以及同樣有他們撰寫的論文《Single image vignetting correction using radial gradient symmetry》有講這方面的算法，不過其實(shí)現(xiàn)的復(fù)雜度較高，即使能編程實(shí)現(xiàn)，速度估計也很慢，其實(shí)用性就不高了。
? ? ? 前不久，偶爾的機(jī)會看到一篇名為《Single-Image Vignetting Correction by Constrained Minimization of log-Intensity Entropy》的論文，并且在github上找到了相關(guān)的一些參考代碼，雖然那個代碼寫的實(shí)在是惡心加無聊，但是對于我來說這并不重要，只要稍有參考，在結(jié)合論文那自己來實(shí)現(xiàn)就不是難事了。
? ? ?論文里的算法核心其實(shí)說起來也沒啥難的，我就我的理解來簡單的描述下：
? ? ?第一：去暗角可以說是陰影校正的一種特例，而將整副圖像的熵最小化也被證明為進(jìn)行陰影校正的一種有效方法，但是普通的熵在優(yōu)化過程中會優(yōu)化到局部最優(yōu)的。因此論文中提出了一種對數(shù)熵的概念（Log-Intensity Entropy），論文中用數(shù)據(jù)做了說明，假設(shè)一副普通正常的圖像其直方圖是單峰分布，那么如果這幅圖像有暗角，其直方圖必然會存在另外一個低明度的分布，如下圖所示：
?
我們校正暗角的過程就是使低明度的分布向原來的正常明度靠近，由上圖第一行的數(shù)據(jù)可以看到，普通的熵計算直到兩個直方圖有部分重疊的時候熵才會下降，之前熵一直都是增加的，而對數(shù)熵則在沒有重疊前至少是保持不增的，因此能夠更好的獲取全局最優(yōu)解。
? ? ?那么論文提出的對數(shù)熵的計算公式為：
? ? ?首先先將亮度進(jìn)行對數(shù)映射，映射公式為： ? ??
? ? ? ? ? ? ? ? ? ? ??
也就是將[0,255]內(nèi)的像素值映射到[0, N-1]內(nèi)，但不是線性映射，而是基于對數(shù)關(guān)系的映射，通常N就是取256，這樣映射后的像素范圍還是[0,255]，但是注意這里的i(L)已經(jīng)是浮點(diǎn)數(shù)了。我們繪制出N等于256時上式的曲線：
? ? ? ? ? ? ? ? ? ??
可見，這種操作實(shí)際上把圖像整體調(diào)亮了。
由于映射后的色階已經(jīng)是浮點(diǎn)數(shù)了，因此，直方圖信息的統(tǒng)計就必須換種方式了，論文給出的公式為：
? ? ? ? ? ? ?
? ?
假設(shè)對于0-255灰度范圍的像素，進(jìn)行如下映射，可以把灰度值映射至[0,N-1]
i(L) = (N ? 1) log(1 + L)/ log 256
假設(shè)N=256，由于映射之后，圖像灰度值為雙精度，而Hist[K]中表示的K必須是整型數(shù)據(jù)。對于每個像素的灰度值大小都是雙精度的圖像，為了統(tǒng)計雙精度的灰度直方圖信息，我們必須進(jìn)行相應(yīng)的加權(quán)處理。、、分別是向下取整，和向上取整，的取值范圍為[0,1]。下圖中橫坐標(biāo)表示的直方圖的K值，縱坐標(biāo)表示的是像素數(shù)量值，那么向上取整和向下取整兩者之差必為1。而恰好是在[0,1]。假設(shè)，，其，那么，對于介于、某一雙精度灰度值而言，我們可以根據(jù)其小數(shù)部分確定該雙精度值在、處的大小。由三角形相似原理可知：
?
???? 參考論文中，
????? 注意取整之后的大小關(guān)系 。向下取整，k小于，向上取整K要大于。
???? 當(dāng)某一雙精度灰度值位于上下取整數(shù)值之間的時候，對向下取整的所對應(yīng)的直方圖bin像素數(shù)量的貢獻(xiàn)量為，對向上取整所對應(yīng)的直方圖bin像素數(shù)量的貢獻(xiàn)量為。相當(dāng)于把位于上下取整之間的x所插值生成的y大小，按照其權(quán)重分別分配給前一個直方圖bin和后一個直方圖bin .
? ? ? 公式很復(fù)雜， 其實(shí)就是有點(diǎn)類似線性插值那種意思，不認(rèn)識了那兩個數(shù)學(xué)符號了，就是向上取整和向下取整。
? ? ? 這樣的對數(shù)熵直方圖信息會由于巨大的色階調(diào)整，導(dǎo)致很多色階是沒有直方圖信息的，一般可以對這個直方圖信息進(jìn)行下高斯平滑的，得到新的直方圖。
最后圖像的對數(shù)熵，計算方法如下：
? ? ? ? ? ? ? ? ? ? ? ? ? ? ??
其中：? ??
??????????????????????????????
? ? ? ? ?
?第二：論文根據(jù)資料《Radiometric alignment and vignetting calibration》提出了一個暗角圖像亮度下降的關(guān)系式，而我去看《Radiometric alignment and vignetting calibration》這篇論文時，他的公式又是從《Vignette and exposure calibration and compensation》中獲取的，所以這個論文的作者寫得文章還不夠嚴(yán)謹(jǐn)。這個公式是一個擁有五個未知參數(shù)的算式，如下所示：
? ? ? ? ? ? ? ? ???
其中：
? ? ? ? ? ? ??
? ? ?其中，x和y是圖像每一點(diǎn)的坐標(biāo)，而則表示暗角的中心位置，他們和a、b、c均為未知量。
?我們可以看到，當(dāng)r=0時，校正系數(shù)為1，即無需校正。當(dāng)r=1時，校正系數(shù)為1+a+b+c。
? ? ?那么經(jīng)過暗角校正后的圖像就為：
? ? ? ? ? ? ?
按照我們的常識，暗角圖像從暗角的中心點(diǎn)想四周應(yīng)該是逐漸變暗的，根據(jù)上式函數(shù)g應(yīng)該是隨著r單調(diào)遞增的（因為我們是校正暗角圖像，所以越靠近邊緣上式的乘法中g(shù)值也就應(yīng)該越大)，因此函數(shù)g的一階導(dǎo)數(shù)應(yīng)該大于0，即：
? ? ? ? ? ? ?
? ? ?同時，我們注意到參數(shù)r的范圍很明顯應(yīng)該在[0,1]之間，這樣上式則可以轉(zhuǎn)換為：
? ? ? ? ? ? ? ??
如果令，則上式變?yōu)?#xff1a;
? ? ? ? ? ?
根據(jù)二次不等式相關(guān)知識，令：
? ? ? ? ?
? ? ?則論文總結(jié)了滿足下述關(guān)系式的a,b,c就能滿足上述要求了：
? ? ? ??
這個我也沒有去驗證了。
? ? ? 第三: 上面描述了校正暗角圖像的公式（帶參數(shù)）以及評價一副圖像是否有暗角的指標(biāo)，那么最后一步就是用這個指標(biāo)來確定公式的參數(shù)。我們未知的參數(shù)有5個，即a、b、c以及暗角的中心點(diǎn)。解這種受限的最優(yōu)問題是有專門的算法的，也是非常計算耗時的。因此，作者提出了一種快速的算法：Hill climbing with rejection of invalid solutions.
? ? ? Hill climbing with rejection of invalid solutions：簡而言之，我們把a(bǔ),b,c三個陰影校正參數(shù)初始化為0，每個參數(shù)獨(dú)立的按照一定步長增加并減小，并計算對應(yīng)的圖像熵。經(jīng)過六次計算之后，我們獲取到了一組使得圖像熵最小的amin,bmin,cmin參數(shù)，然后從這開始再次向之前方法一樣計算，直到找到圖像最小熵對應(yīng)的三個參數(shù)。
? ? ? 首先，很明顯，為了計算這些最優(yōu)參數(shù)，我們沒有必要直接在原圖大小上直接計算，這點(diǎn)在原論文也有說明，我們即使把他們的寬高分別縮小到原圖的1/5甚至1/10計算出來的結(jié)果也不會有太大的差異，而這些參數(shù)的差異對最終的的結(jié)果影響也不大，但是計算量就能減少到原來的1/25和1/100。
? ? ? 計算出上述a、b、c以及中心點(diǎn)后，就可以再次按照校正公式來進(jìn)行校正了，注意暗角的影響對每個通道都是等同的，因此，每個通道都應(yīng)該乘以相同的值。
下面貼出一些用論文中的算法處理的結(jié)果圖：
? ?
??
? ?
? ? ??
代碼實(shí)現(xiàn)：
#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; }//檢查參數(shù)的合法性 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; }?結(jié)論：參考相關(guān)文獻(xiàn)和代碼，最終計算出來的圖像總是不理想，存在數(shù)據(jù)溢出的情況。不知道哪出錯了。
參考文獻(xiàn)：
1.圖像增強(qiáng)系列之圖像自動去暗角算法。
2.https://github.com/HJCYFY/Vignetting-Correction? 論文對應(yīng)的代碼
?
總結(jié)
- 上一篇： 002_2 gtsam/unstable
- 下一篇： 时间复杂度（算法的渐进时间复杂度）