图像去暗角
图像去暗角
暗角图像是一种在现实中较为常见的图像，其主要特征就是在图像四个角有较为显著的亮度下降，比如下面两幅图。根据其形成的成因，主要有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 brightness
value[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 = 4
TempHist[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));
// smooth
for (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;
}
结论：参考相关文献和代码，最终计算出来的图像总是不理想，存在数据溢出的情况。不知道哪出错了。
参考文献：
2.https://github.com/HJCYFY/Vignetting-Correction 论文对应的代码