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

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

相关文章

来自专栏AI科技大本营的专栏

朴素贝叶斯分类器详解及中文文本舆情分析(附代码实践)

作者 | 杨秀璋(笔名:Eastmount),贵州财经大学信息学院老师,硕士毕业于北京理工大学,主要研究方向是Web数据挖掘、知识图谱、Python数据分析、图...

643
来自专栏YoungGy

机器为什么可以学习

机器学习、人工智能炙手可热,但是机器到底为什么可以学习呢?本文将从霍夫丁不等式讲到VC维,探究机器学习的原因所在。 机器怎么可能学习 机器确实可能会学习 机器为...

1995
来自专栏数据科学与人工智能

【陆勤践行】SVM之核函数

前段时间热映的《星际穿越》想必大家都看过,在这部烧脑大片中,主角库珀进入到了高维度空间,在那里,时间这个维度变成实体存在,人们可以像散步一样沿着时间这个维度来回...

2058
来自专栏大数据挖掘DT机器学习

用R语言写个贝叶斯模型 预测我的妻子是否怀孕

在2015年的二月21日,我的妻子已经33天没有来月经了,她怀孕了,这真是天大的好消息! 通常月经的周期是大约一个月,如果你们夫妇打算怀孕,那么月经没来或许是一...

2949
来自专栏人工智能

决策树及ID3算法学习

决策树是一种用树形结构来辅助行为研究、决策分析以及机器学习的方式,是机器学习中的一种基本的分类方法。

1.5K16
来自专栏企鹅号快讯

Python数据分析与实战挖掘

基础篇 书推荐:《用python做科学计算》 ? 扩展库 简介 Numpy数组支持,以及相应的高效处理函数 Scipy矩阵支持,以及相应的矩阵数值计算模块 Ma...

3525
来自专栏机器学习之旅

Kaggle&TianChi分类问题相关纯算法理论剖析导读Bias-Variance-TradeofGBDT 理论剖析GBDT 实战剖析

17/12/30-update :很多朋友私密我想要代码,甚至利用金钱诱惑我,好吧,我沦陷了。因为原始代码涉及到公司的特征工程及一些利益trick,所以我构造了...

753
来自专栏人工智能

卷积神经网络学习笔记

1.卷积神经网络的图像识别原理: 通过过滤函数 来描绘出图像的边界: 过滤函数和图像相同区域的数值进行相乘,得到新的图像, 新图像则只剩下边图像。 cros...

21610
来自专栏大数据挖掘DT机器学习

逻辑回归之基础知识及手写数字识别实例

作者:daniel-D 原文:http://www.cnblogs.com/daniel-D/   这学期 Pattern Recognition 课程的 p...

5499
来自专栏红色石头的机器学习之路

Coursera吴恩达《构建机器学习项目》课程笔记(1)-- 机器学习策略(上)

《Structuring Machine Learning Projects》(构建机器学习项目)这门课是Andrw Ng深度学习专项课程中的第三门课。这门课主...

2170

扫码关注云+社区