图像去暗角

最新推荐文章于 2022-04-26 17:37:57 发布

大熊背

最新推荐文章于 2022-04-26 17:37:57 发布

阅读量4.1k

点赞数 3

分类专栏：图像基础算法处理文章标签：阴影校正

本文链接：https://blog.csdn.net/lz0499/article/details/82054651

版权

图像基础算法处理专栏收录该内容

63 篇文章 123 订阅

订阅专栏

图像去暗角

暗角图像是一种在现实中较为常见的图像，其主要特征就是在图像四个角有较为显著的亮度下降，比如下面两幅图。根据其形成的成因，主要有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必须是整型数据。对于每个像素的灰度值大小都是双精度的图像，为了统计双精度的灰度直方图信息，我们必须进行相应的加权处理。 $\left \lfloor iL(x,y)) \right \rfloor$ 、 $\left \lceil iL(x,y)) \right \rceil$ 、 $\mu$ 分别是向下取整，和向上取整， $\mu$ 的取值范围为[0,1]。下图中横坐标表示的直方图的K值，纵坐标表示的是像素数量值，那么向上取整和向下取整两者之差必为1。而 $\mu$ 恰好是在[0,1]。假设， $\left \lfloor iL(x,y)) \right \rfloor=x_{1}$ $\left \lceil iL(x,y)) \right \rceil=x_{2}$ ，其 $f(x_{1})=y_{_{1}}f(x_{2})=y_{2}$ ，那么，对于介于 $\left \lfloor iL(x,y)) \right \rfloor$ 、 $\left \lceil iL(x,y)) \right \rceil$ 某一双精度灰度值而言，我们可以根据其小数部分确定该双精度值在 $x_{1},x_{2}$ 、处的大小。由三角形相似原理可知：

参考论文中，

注意取整之后的大小关系。向下取整，k小于 iL(x,y) ，向上取整K要大于 iL(x,y) 。

当某一双精度灰度值位于上下取整数值之间的时候，对向下取整的所对应的直方图bin像素数量的贡献量为 $(1-\mu )y_{1}$ $（1-\mu）y_{1}$ ，对向上取整所对应的直方图bin像素数量的贡献量为 $\mu y_{2}$ 。相当于把位于上下取整之间的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;   
}

结论：参考相关文献和代码，最终计算出来的图像总是不理想，存在数据溢出的情况。不知道哪出错了。

参考文献：

1.图像增强系列之图像自动去暗角算法。

2.https://github.com/HJCYFY/Vignetting-Correction 论文对应的代码