暗角图像是一种在现实中较为常见的图像，其主要特征就是在图像四个角有较为显著的亮度下降，比如下面两幅图。根据其形成的成因，主要有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; //` 粗选ABC三个变量的步长 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++) // 多次测试，表面最优的ABC的范围均在[-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
