Sauvola算法公式推导
近日微博上有一位网友针对《Sauvola 实现图像二值化》一文问我有没有看懂Sauvola算法。
实话实说，我确实是不懂Sauvola算法的原理，在此之前，我也从来有没有尝试过想去看懂它。看到这位网友的评论后，我找到并学习了一些关于Sauvola算法的资料，目前仍然没有完全理解Sauvola算法的精髓，但我找到了《Sauvola 实现图像二值化》一文中所发布的源代码公式的来源，并验证了该源代码与Sauvola算法的原理公式基本相符，在此分享我的学习心得与公式推导过程，欢迎讨论。
1.我眼中的Sauvola算法原理
参考文献1：《Efficient multiscale Sauvola’s binarization》 密码: diu8，该文献先给出Sauvola算法的计算公式，然后提出了该算法的几点局限性并分别给出了改善对策；
参考文献2：《Adaptive document image binarization》密码: fhaj，该文献的作者是J. Sauvola，它分别讲述了图像中非文字部分与文字部分的二值化方法，文字部分的二值化方法即《Sauvola 实现图像二值化》中代码所实现的Sauvola算法。该文献中提到Sauvola算法是Niblack算法的改进版本，有两个公式来源于Niblack算法。《An Introduction to Image Processing》(作者：W. Niblack，1986年出版)一书中有详细讲述Niblack算法，可惜我遍寻网络，也没有找到该书的电子版，全国图书馆咨询联盟，谷歌图书以及多个外文书库都仅仅是收录了书名，并无内容的电子版。
参考文献3：《基于改进 Otsu和 Niblack的图像二值化方法》密码: mfw4，还好我在文献3中找到了Sauvola算法中沿用的Niblack算法的公式。
我们常见的图像二值化算法大致可分为全局阈值方法与局部阈值方法这两种类型。其中OTSU算法是全局阈值的代表，而Sauvola算法则是局部阈值方法的标杆。Sauvola算法的输入是灰度图像，它以当前像素点为中心，根据当前像素点邻域内的灰度均值与标准方差来动态计算该像素点的阈值。
假定当前像素点的坐标为(x,y)，以该点为中心的领域为r*r，g(x,y)表示(x,y)处的灰度值，Sauvola算法的步骤为：
2.Sauvola算法C#代码讲解
本文中Sauvola算法的代码更新了部分注释，给出阈值计算公式的推导过程：
本文贴出的代码在图像尺寸比较大时，图像积分的计算值会超出图像积分存储量的数据类型。若您需要用本文发布的代码，还请检查图像尺寸是否会令代码中相关变量溢出，若有请自行修改为合适的变量，谢谢！
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 |
/// <summary> /// 实现Sauvola算法实现图像二值化 /// </summary> /// <param name="bin_image">用于存储二值化完成的图像</param> /// <param name="gray_image">用于存储等待二值化完成的灰度图像</param> public static void Sauvola(Byte[,] bin_image, Byte[,] gray_image) { //以当前像素点为中心的邻域的宽度 int w = 40; //使用者自定义的修正系数 double k = 0.3; //邻域边界距离中心点的距离 int whalf = w >> 1; int MAXVAL = 256; int image_width = gray_image.GetLength(0); int image_height = gray_image.GetLength(1); int[,] integral_image = new int[image_width, image_height]; int[,] integral_sqimg = new int[image_width, image_height]; int[,] rowsum_image = new int[image_width, image_height]; int[,] rowsum_sqimg = new int[image_width, image_height]; int xmin,ymin,xmax,ymax; double diagsum,idiagsum,diff,sqdiagsum,sqidiagsum,sqdiff,area; double mean,std,threshold; for (int j = 0; j < image_height; j++) { rowsum_image[0, j] = gray_image[0, j]; rowsum_sqimg[0, j] = gray_image[0, j] * gray_image[0, j]; } for (int i = 1; i < image_width; i++) { for (int j = 0; j < image_height; j++) { //计算图像范围内任意宽度窗口(邻域)的灰度值之和 rowsum_image[i, j] = rowsum_image[i - 1, j] + gray_image[i, j]; //计算图像范围内任意宽度窗口(邻域)的灰度值平方之和 rowsum_sqimg[i, j] = rowsum_sqimg[i - 1, j] + gray_image[i, j] * gray_image[i, j]; } } for (int i = 0; i < image_width; i++) { integral_image[i, 0] = rowsum_image[i, 0]; integral_sqimg[i, 0] = rowsum_sqimg[i, 0]; } for (int i = 0; i < image_width; i++) { for (int j = 1; j < image_height; j++) { //计算图像范围内任意宽度窗口(邻域)的灰度值的积分 integral_image[i, j] = integral_image[i, j - 1] + rowsum_image[i, j]; //计算图像范围内任意宽度窗口(邻域)的灰度值平方的积分 integral_sqimg[i, j] = integral_sqimg[i, j - 1] + rowsum_sqimg[i, j]; } } //Calculate the mean and standard deviation using the integral image for(int i=0; i<image_width; i++){ for(int j=0; j<image_height; j++){ xmin = Math.Max(0,i-whalf); ymin = Math.Max(0, j - whalf); xmax = Math.Min(image_width - 1, i + whalf); ymax = Math.Min(image_height - 1, j + whalf); area = (xmax-xmin+1)*(ymax-ymin+1); // area can't be 0 here // proof (assuming whalf >= 0): // we'll prove that (xmax-xmin+1) > 0, // (ymax-ymin+1) is analogous // It's the same as to prove: xmax >= xmin // image_width - 1 >= 0 since image_width > i >= 0 // i + whalf >= 0 since i >= 0, whalf >= 0 // i + whalf >= i - whalf since whalf >= 0 // image_width - 1 >= i - whalf since image_width > i // --IM if (area <= 0) throw new Exception("Binarize: area can't be 0 here"); if (xmin == 0 && ymin == 0) { // Point at origin diff = integral_image[xmax, ymax]; sqdiff = integral_sqimg[xmax, ymax]; } else if (xmin == 0 && ymin > 0) { // first column diff = integral_image[xmax, ymax] - integral_image[xmax, ymin - 1]; sqdiff = integral_sqimg[xmax, ymax] - integral_sqimg[xmax, ymin - 1]; } else if (xmin > 0 && ymin == 0) { // first row diff = integral_image[xmax, ymax] - integral_image[xmin - 1, ymax]; sqdiff = integral_sqimg[xmax, ymax] - integral_sqimg[xmin - 1, ymax]; } else { // rest of the image diagsum = integral_image[xmax, ymax] + integral_image[xmin - 1, ymin - 1]; idiagsum = integral_image[xmax, ymin - 1] + integral_image[xmin - 1, ymax]; //以(i,j)为中心点的w邻域内灰度值的积分 diff = diagsum - idiagsum; sqdiagsum = integral_sqimg[xmax, ymax] + integral_sqimg[xmin - 1, ymin - 1]; sqidiagsum = integral_sqimg[xmax, ymin - 1] + integral_sqimg[xmin - 1, ymax]; //以(i,j)为中心点的w邻域内灰度值平方的积分 sqdiff = sqdiagsum - sqidiagsum; } //以(i,j)为中心点的w邻域内的灰度均值 mean = diff/area; //以(i,j)为中心点的w邻域内的标准方差 std = Math.Sqrt((sqdiff - diff*diff/area)/(area-1)); //根据Sauvola计算公式和以(i,j)为中心点的w邻域内的灰度均值与标准方差来计算当前点(i,j)的二值化阈值 threshold = mean*(1+k*((std/128)-1)); //根据当前点的阈值对当前像素点进行二值化 if(gray_image[i,j] < threshold) bin_image[i,j] = 0; else bin_image[i,j] = (byte)(MAXVAL-1); } } } |
Comments
发表回复
此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据。
这个算法中的计算平方积分图的时候存在数据越界的情况
您好，请问您说的是本文发布的代码中存在数组越界吗？可以麻烦您指出是哪一块吗？首先，我用这段代码二值化过多个图片，并未报错；另外，我刚才有仔细查看一遍代码，计算平方积分图与计算领域灰度均值数组下标均一致，并未发现有越界的状况。可能是我眼拙，麻烦指点，谢谢！
如果图像尺寸比较大的时候，积分图的计算会发生超出数据定义的类型！是数据越界，不是数组越界！
谢谢指点。在图像尺寸比较大的时候，确实存在这个问题。
std = Math.Sqrt((sqdiff – diff*diff/area)/(area-1));
这行代码貌似跟推论不一致。应该怎么理解呢？
朋友，你的问题难倒我了，程序来自于github，公式推导来自于各相关文献，我从来没想过为啥程序里作为除数的area需要减一，刚才我重新看了一遍，也没有想到什么合适的理由，抱歉。不过我觉得减一或不减一应该都可以，有时间我会验证一下。