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

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

??

?

　　

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

? ? ? 前不久，偶爾的機(jī)會(huì)看到一篇名為《Single-Image Vignetting Correction by Constrained Minimization of log-Intensity Entropy》的論文，并且在github上找到了相關(guān)的一些參考代碼，雖然那個(gè)代碼寫(xiě)的實(shí)在是惡心加無(wú)聊，但是對(duì)于我來(lái)說(shuō)這并不重要，只要稍有參考，在結(jié)合論文那自己來(lái)實(shí)現(xiàn)就不是難事了。

? ? ?論文里的算法核心其實(shí)說(shuō)起來(lái)也沒(méi)啥難的，我就我的理解來(lái)簡(jiǎn)單的描述下：

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

?

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

? ? ?那么論文提出的對(duì)數(shù)熵的計(jì)算公式為：

? ? ?首先先將亮度進(jìn)行對(duì)數(shù)映射，映射公式為： ? ??

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

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

? ? ? ? ? ? ? ? ? ??

　　可見(jiàn)，這種操作實(shí)際上把圖像整體調(diào)亮了。

　　由于映射后的色階已經(jīng)是浮點(diǎn)數(shù)了，因此，直方圖信息的統(tǒng)計(jì)就必須換種方式了，論文給出的公式為：

? ? ? ? ? ? ?

? ?

假設(shè)對(duì)于0-255灰度范圍的像素，進(jìn)行如下映射，可以把灰度值映射至[0,N-1]

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

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

?

???? 參考論文中，

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

???? 當(dāng)某一雙精度灰度值位于上下取整數(shù)值之間的時(shí)候，對(duì)向下取整的所對(duì)應(yīng)的直方圖bin像素?cái)?shù)量的貢獻(xiàn)量為，對(duì)向上取整所對(duì)應(yīng)的直方圖bin像素?cái)?shù)量的貢獻(xiàn)量為。相當(dāng)于把位于上下取整之間的x所插值生成的y大小，按照其權(quán)重分別分配給前一個(gè)直方圖bin和后一個(gè)直方圖bin .

? ? ? 公式很復(fù)雜， 其實(shí)就是有點(diǎn)類(lèi)似線(xiàn)性插值那種意思，不認(rèn)識(shí)了那兩個(gè)數(shù)學(xué)符號(hào)了，就是向上取整和向下取整。

? ? ? 這樣的對(duì)數(shù)熵直方圖信息會(huì)由于巨大的色階調(diào)整，導(dǎo)致很多色階是沒(méi)有直方圖信息的，一般可以對(duì)這個(gè)直方圖信息進(jìn)行下高斯平滑的，得到新的直方圖。

　　最后圖像的對(duì)數(shù)熵，計(jì)算方法如下：

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

　　其中：? ??

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

? ? ? ? ?

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

? ? ? ? ? ? ? ? ???

　　其中：

? ? ? ? ? ? ??

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

　 ?我們可以看到，當(dāng)r=0時(shí)，校正系數(shù)為1，即無(wú)需校正。當(dāng)r=1時(shí)，校正系數(shù)為1+a+b+c。

? ? ?那么經(jīng)過(guò)暗角校正后的圖像就為：

? ? ? ? ? ? ?

　　按照我們的常識(shí)，暗角圖像從暗角的中心點(diǎn)想四周應(yīng)該是逐漸變暗的，根據(jù)上式函數(shù)g應(yīng)該是隨著r單調(diào)遞增的（因?yàn)槲覀兪切Ｕ到菆D像，所以越靠近邊緣上式的乘法中g(shù)值也就應(yīng)該越大)，因此函數(shù)g的一階導(dǎo)數(shù)應(yīng)該大于0，即：

? ? ? ? ? ? ?

? ? ?同時(shí)，我們注意到參數(shù)r的范圍很明顯應(yīng)該在[0,1]之間，這樣上式則可以轉(zhuǎn)換為：

? ? ? ? ? ? ? ??

　　如果令，則上式變?yōu)?#xff1a;

? ? ? ? ? ?

　　根據(jù)二次不等式相關(guān)知識(shí)，令：

? ? ? ? ?

? ? ?則論文總結(jié)了滿(mǎn)足下述關(guān)系式的a,b,c就能滿(mǎn)足上述要求了：

? ? ? ??

　　這個(gè)我也沒(méi)有去驗(yàn)證了。

? ? ? 第三: 上面描述了校正暗角圖像的公式（帶參數(shù)）以及評(píng)價(jià)一副圖像是否有暗角的指標(biāo)，那么最后一步就是用這個(gè)指標(biāo)來(lái)確定公式的參數(shù)。我們未知的參數(shù)有5個(gè)，即a、b、c以及暗角的中心點(diǎn)。解這種受限的最優(yōu)問(wèn)題是有專(zhuān)門(mén)的算法的，也是非常計(jì)算耗時(shí)的。因此，作者提出了一種快速的算法：Hill climbing with rejection of invalid solutions.

? ? ? Hill climbing with rejection of invalid solutions：簡(jiǎn)而言之，我們把a(bǔ),b,c三個(gè)陰影校正參數(shù)初始化為0，每個(gè)參數(shù)獨(dú)立的按照一定步長(zhǎng)增加并減小，并計(jì)算對(duì)應(yīng)的圖像熵。經(jīng)過(guò)六次計(jì)算之后，我們獲取到了一組使得圖像熵最小的amin,bmin,cmin參數(shù)，然后從這開(kāi)始再次向之前方法一樣計(jì)算，直到找到圖像最小熵對(duì)應(yīng)的三個(gè)參數(shù)。

? ? ? 首先，很明顯，為了計(jì)算這些最優(yōu)參數(shù)，我們沒(méi)有必要直接在原圖大小上直接計(jì)算，這點(diǎn)在原論文也有說(shuō)明，我們即使把他們的寬高分別縮小到原圖的1/5甚至1/10計(jì)算出來(lái)的結(jié)果也不會(huì)有太大的差異，而這些參數(shù)的差異對(duì)最終的的結(jié)果影響也不大，但是計(jì)算量就能減少到原來(lái)的1/25和1/100。

? ? ? 計(jì)算出上述a、b、c以及中心點(diǎn)后，就可以再次按照校正公式來(lái)進(jìn)行校正了，注意暗角的影響對(duì)每個(gè)通道都是等同的，因此，每個(gè)通道都應(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; }//計(jì)算圖像熵 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)和代碼，最終計(jì)算出來(lái)的圖像總是不理想，存在數(shù)據(jù)溢出的情況。不知道哪出錯(cuò)了。

參考文獻(xiàn)：

1.圖像增強(qiáng)系列之圖像自動(dòng)去暗角算法。

2.https://github.com/HJCYFY/Vignetting-Correction? 論文對(duì)應(yīng)的代碼

?

總結(jié)

以上是生活随笔為你收集整理的图像去暗角的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。