欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

图像去暗角

程序员文章站 2024-03-03 19:17:10
...

                                                    图像去暗角

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

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

参考文献:

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

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

 

相关标签: 阴影校正