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

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

  可见,这种操作实际上把图像整体调亮了。

  由于映射后的色阶已经是浮点数了,因此,直方图信息的统计就必须换种方式了,论文给出的公式为:

      公式很复杂, 其实就是有点类似线性插值那种意思,不认识了那两个数学符号了,就是向上取整和向下取整。

      这样的对数熵直方图信息会由于巨大的色阶调整,导致很多色阶是没有直方图信息的,一般可以对这个直方图信息进行下高斯平滑的,得到新的直方图。

  最后图像的对数熵,计算方法如下:

  其中:    

第二:论文根据资料《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.

      我稍微看了下这个算法,确实是个不错的想法,不过我并没有去实践,我采用了另外一种粗略的优化方式。

      首先,很明显,为了计算这些最优参数,我们没有必要直接在原图大小上直接计算,这点在原论文也有说明,我们即使把他们的宽高分别缩小到原图的1/5甚至1/10计算出来的结果也不会有太大的差异,而这些参数的差异对最终的的结果影响也不大,但是计算量就能减少到原来的1/25和1/100。

     接着,我们观察到a、b以及c的最优结果范围一般都在-2和2之间,并且从g的计算公式中知道,由于r是属于0和1之间的正数,r^2, r^4, r^6在数值递减的非常快,比如r=0.8,则三者对应的结果就分别为0.64、0.4096、0.2621,因此,a和b及c在公式中的对最后结果的影响也依次越来越小。

     那么,我们可以参考以前的对比度保留之彩色图像去色算法---基础算法也可以上档次一文中的优化方式,把a, b ,c 三个参数分别在[-2,2]之间离散化,考虑到参数稍微差异不会对结果有太大的影响,以及a、b、c的重要性,我们可以设置a、b、c三者的离散间隔分别为0.2、0.3、0.4,然后综合上述判断a、b、c是否为合理组合的函数,离散取样的计算量组合大概有300种可能,对小图计算着300种可能性的耗时是完全可以接受的,甚至考虑极端一点,把c的计算放到循环外侧,即C取固定值0,然后计算出优选的a和b值后,在计算C值。

      上述计算过程并未考虑暗角中心点的范围,我们是固定把暗角的中心点放置在图像的正中心位置的,即 (Width/2, Height /2),实际上,对于大部分拍摄的图来说,暗角就是位于中心位置的,因此这种假设也无可厚非,因为暗角中心计算的增加必然会严重增加计算量, 为了求出暗角中心的合理位置,我们在计算出上述a、b、c后,在小图中以一定步长按照公式计算出粗略的中心位置,再放大到原图中去。

      计算出上述a、b、c以及中心点后,就可以再次按照校正公式来进行校正了,注意暗角的影响对每个通道都是等同的,因此,每个通道都应该乘以相同的值。

  下面贴出一些用论文中的算法处理的结果图:

     注意到上面最后一副图的结果,那个女的婚纱以及衣服那些地方已经严重的过曝了,我不清楚理论上造成整个原因的是什么,但是如果把计算i(L)的公式中的N修改为小一点的值,比如64,则可以避免这个结果。

     github上的那个代码则对这个对数熵的过程做了一点改造,这个改造相当暴力,就是什么呢,他把原来的[0,255]直接量化为8个等级,量化的依据是整形LOG2函数,即0->0,[1, 2]->1,[3, 4]->3,[5,8]->4,[9,16]->5,[17,32]->6,[33,64]->7,[65,128]->8,[129,255]->9, 原来的一条曲线映射函数变成了阶跃函数了。这样直方图实际上只有9个值了,那么也不需要什么直方图插值和高斯模糊了,直方图则可以用整形表示,相对来说速度也能有很大的提升,并且也能克服上述最后一张图片出现的那个瑕疵,其结果如下:

      最后我们贴一些代码对上述过程予以解释:

   第一个是判断a、b、c是否为合理值的函数:

//    按论文中公式18得条件判断是否是合理的参数
bool IsValidABC_Original(float A, float B, float C)
{
    const int MAX_BRIGHTNESS_MULTIPLICATION = 3;
    if ((1 + A + B + C) > MAX_BRIGHTNESS_MULTIPLICATION)    return false;        //    当r==1时,出现最大的亮度调整
    if (C == 0)
    {
        if (A <= 0 || (A + 2 * B <= 0))            //    如果C==0,则根据公式(15)知,当r==0时,A必须大于0,而当r==1时,A+2B必须大于0
            return false;
    }
    else
    {
        float D = 4 * B * B - 12 * A * C;
        float qMins = (-2 * B - sqrtf(D)) / (6 * C);        //    公式(17)
        float qPlus = (-2 * B + sqrtf(D)) / (6 * C);
        if (C < 0)
        {
            if (D >= 0)
            {
                if (qMins > 0 || qPlus < 1)
                    return false;
            }
            else
                return false;
        }
        else
        {
            if (D >= 0)
            {
                if (!((qMins <= 0 && qPlus <= 0) || (qMins >= 1 && qPlus >= 1)))
                    return false;
            }
        }
    }
    return true;
}

 可见除了文中公式的一些限制,我们还增加了几个额外的小限制,比如最大亮度调节比例为3等等。

    第二个是计算指定参数下计算对数熵的过程:

//    计算不同参数修复后的图像的整体对数熵
float CalculateEntropyFromImage(unsigned char *Src, int Width, int Height, float A, float B, float C, int CenterX, int CenterY)
{
    float Histgram[256] = { 0 };
    float Invert = 1.0f / (CenterX * CenterX + CenterY * CenterY + 1);
    float Mul_Factor = 256.f / log(256.0f);
    for (int Y = 0; Y < Height; Y++)
    {
        unsigned char *LinePS = Src + Y * Width * 4;
        int SquareY = (Y - CenterY) * (Y - CenterY);
        for (int X = 0; X < Width; X++)
        {
            int Intensity = (LinePS[0] + (LinePS[1] << 1) + LinePS[2]) >> 2;         //    公式(2)
            int RadiusSqua2 = (X - CenterX) * (X - CenterX) + SquareY;
            float R = RadiusSqua2 * Invert;                        //    公式(12)
            float Gain = 1 + (R * (A + R * (B + R * C)));            //    gain = 1 + a * r^2 + b * r^4 + c * r^6 ,公式(11)
            float Correction = Gain * Intensity;                    //    直接校正后的结果值
            if (Correction >= 255)
            {
                Correction = 255;                //    It is possible that, due to local intensity increases applied by devignetting, the corrected image intensity range exceeds 255.    
                Histgram[255]++;                //    In this case the algorithm simply adds histogram bins at the upper end without rescaling the distribution,
            }
            else
            {
                float Pos = Mul_Factor * log(Correction + 1);    //    公式(6)
                int Int = int(Pos);
                Histgram[Int] += 1 - (Pos - Int);    //    公式(7)
                Histgram[Int + 1] += Pos - Int;
            }
            LinePS += 4;
        }
    }
    float TempHist[256 + 2 * 4];            //    SmoothRadius = 4
    TempHist[0] = Histgram[4];                TempHist[1] = Histgram[3];    
    TempHist[2] = Histgram[2];                TempHist[3] = Histgram[1];
    TempHist[260] = Histgram[254];            TempHist[261] = Histgram[253];
    TempHist[262] = Histgram[252];            TempHist[263] = Histgram[251];
    memcpy(TempHist + 4, Histgram, 256 * sizeof(float));
    
    for (int X = 0; X < 256; X++)            //    公式(8),进行一个平滑操作
        Histgram[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;

    return CalculateEntropyFromHistgram_Original(Histgram, 256);
}

  其中计算熵的函数为:

//    从直方图中计算熵值,原论文中直方图肯定是浮点数
float CalculateEntropyFromHistgram_Original(float Histgram[], int Length)
{
    float Sum = 0;
    for (int X = 0; X < Length; X++)
    {
        Sum += Histgram[X];
    }
    float Entropy = 0;
    for (int X = 0; X < Length; X++)
    {
        if (Histgram[X] == 0) continue;
        float p = (float)Histgram[X] / Sum;
        Entropy += p * logf(p);
    }
    return -Entropy;
}

  其中

    int Int = int(Pos);
    Histgram[Int] += 1 - (Pos - Int);    //    公式(7)
    Histgram[Int + 1] += Pos - Int;

 就是公式7所描述的过程的实现。

  论文中的高斯模糊,我这里只是借助了一个简单的线性模糊来代替,这个不会对结果造成本质的区别。

     最后图像的校正代码大概如下:

int Devignetting_Original(unsigned char *Src, unsigned char *Dest, int Width, int Height)
{
    if ((Src == NULL) || (Dest == NULL))        return STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))            return STATUS_INVALIDPARAMETER;

    const float Step = 0.2f;        //`    粗选A\B\C三个变量的步长

    float SmallestEntropy = 1000000000000.0f;
    float A = 0, B = 0, C = 0;            
    int CenterX = Width / 2, CenterY = Height / 2;        //    中心就默认为图片中心

    for (int X = -10; X <= 10; X++)        //    多次测试,表面最优的A\B\C的范围均在[-2,2]之间
    {
        for (int Y = -10; Y <= 10; Y++)
        {
            for (int Z = -10; Z <= 10; Z++)
            {
                if (IsValidABC_Original(X * Step, Y * Step, Z * Step) == true)    //    判断这个组合时候有效
                {
                    float Entropy = CalculateEntropyFromImage(Src, Width, Height, X * Step, Y * Step, Z * Step, CenterX, CenterY);
                    if (Entropy < SmallestEntropy)                                    //    取熵值最小的
                    {
                        A = X * Step;
                        B = Y * Step;
                        C = Z * Step;
                        SmallestEntropy = Entropy;
                    }
                }
            }
        }
    }
    float Invert = 1.0 / (CenterX * CenterX + CenterY * CenterY + 1);
    for (int Y = 0; Y < Height; Y++)
    {
        byte *LinePS = Src + Y * Width * 4;
        byte *LinePD = Dest + Y * Width * 4;
        int SquareY = (Y - CenterY) * (Y - CenterY);
        for (int X = 0; X < Width; X++)
        {
            int RadiusSqua2 = (X - CenterX) * (X - CenterX) + SquareY;
            float R2 = RadiusSqua2 * Invert;                                    //    公式12
            float Gain = 1 + (R2 * (A + R2 * (B + R2 * C)));                    //    公式11                
            LinePD[0] = ClampToByte(LinePS[0] * Gain);            //    增益
            LinePD[1] = ClampToByte(LinePS[1] * Gain);
            LinePD[2] = ClampToByte(LinePS[2] * Gain);
            LinePD += 4;
            LinePS += 4;
        }
    }
    return STATUS_OK;
}

上面的代码是未经特别的优化的,只是表达了大概的意思,并且把暗角的中心点默认为图像的中心点,如果暗角的中心点不是图像中心,要注意计算r时可能会出现r>1的情况,这个时候一定要注意重置r=1,否则结果就会完全不对了。

经过测试,对于没有暗角的图像,一般来说该算法不会对图片产生很大的影响,很多图片基本是无变换的,有些可能会有点区别,也就是整体变亮一点而已,因此,还是有比较好的普适性的。

     由于论文中的暗角强度减弱公式是根据一些光学原理获得的,其可能并不符合一些软件自己添加的暗角的规律,所以如果你用这中测试图去测试算法,可能不会获得非常满意的结果,

     算法速度上,经过我上面的描述的优化后的算法,对于800*600的彩色图,一般能在15ms左右处理完成,而未优化的上述代码可能要10000ms。 基本上算法的核心代码为已经贴出,这里不共享我优化后的快速实现,有能力的朋友自然能快速搞定他们,共享一个测试工程吧。

http://files.cnblogs.com/files/Imageshop/Devignetting_Test.rar

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏大数据文摘

利用 Scikit Learn的Python数据预处理实战指南

2276
来自专栏Python数据科学

【Python数据分析基础】: 数据缺失值处理

再好的模型,如果没有好的数据和特征质量,那训练出来的效果也不会有所提高。数据质量对于数据分析而言是至关重要的,有时候它的意义会在某种程度上会胜过模型算法。

933
来自专栏机器学习算法全栈工程师

【TPAMI重磅综述】 SIFT与CNN的碰撞:万字长文回顾图像检索任务十年探索历程(上篇)

基于内容的图像检索任务(CBIR)长期以来一直是计算机视觉领域重要的研究课题,自20世纪90年代早期以来,研究人员先后设计了图像的全局特征,局部特征,卷积特征的...

784
来自专栏ATYUN订阅号

机器学习:Python测试线性可分性的方法

? 线性和非线性分类 两个子集是线性可分的,如果存在一个超平面将每组的元素的所有元素的一组驻留在另一侧的超平面其他设置。我们可以描述它在2D绘图中通过分离线,...

3465
来自专栏机器之心

LSTM、GRU与神经图灵机:详解深度学习最热门的循环神经网络

选自MachineLearningMastery 作者:Jason Brownlee 机器之心编译 参与:熊猫 循环神经网络是当前深度学习热潮中最重要和最核心的...

3429
来自专栏瓜大三哥

YCbCr数据格式

  主要的采样格式有YCbCr 4:2:0、YCbCr 4:2:2、YCbCr4:1:1和 YCbCr 4:4:4。其中YCbCr 4:1:1 比较常用,其含义...

893
来自专栏AI科技评论

想了解递归神经网络?这里有一份入门教程

本文来自deeplearning4j,AI科技评论编辑。 递归神经网络是一类人工神经网络,用于识别诸如文本、基因组、手写字迹、语音等序列数据的模式,或用于识别传...

2933
来自专栏SIGAI学习与实践平台

场景文本检测—CTPN算法介绍

原创声明:本文为 SIGAI 原创文章,仅供个人学习使用,未经允许,不得转载,不能用于商业目的。

3693
来自专栏人工智能头条

优秀的排序算法如何成就了伟大的机器学习技术(视频+代码)

【导读】在机器学习中,支持向量机(SVM)算法是针对二分类任务设计的,可以分析数据,识别模式,用于分类和回归分析。训练算法构建一个模型,将新示例分配给一个类别或...

672
来自专栏IT派

随机森林之美

导语:随机森林和决策树相比,能更好的防止过拟合。虽然每个基分类器很弱,但最后组合的结果通常很强,这也类似于:“三个臭皮匠顶个诸葛亮”的思想。对比发现Random...

3544

扫码关注云+社区