图像去暗角

                                                    图像去暗角

暗角图像是一种在现实中较为常见的图像,其主要特征就是在图像四个角有较为显著的亮度下降,比如下面两幅图。根据其形成的成因,主要有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  论文对应的代码

 

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值