现代图像处理算法教程-一-

布客飞龙
132 阅读58分钟

现代图像处理算法教程（一）

原文：Modern Algorithms for Image Processing

协议：CC BY-NC-SA 4.0

一、简介

这本书包含图像处理算法的描述，如降噪，包括减少脉冲噪声，对比度增强，阴影校正，边缘检测，等等。该书的配套网站上包含了用 C#编程语言实现算法的项目源代码。源代码是 Windows 窗体项目，而不是 Microsoft 基础类库(MFC)项目。这些项目中的控件和图形是通过简单易懂的方法实现的。我选择了这种实现控件和图形服务的方式，而不是基于 MFC 的方式，因为使用 MFC 的集成开发环境(IDE)非常昂贵。此外，使用 MFC 的软件相当复杂。它在一个项目中包括许多不同的文件，用户很大程度上无法理解这些文件的意义和用途。相反，Windows 窗体及其实用工具是免费的，更容易理解。它们提供了类似于 MFC 的控件和图形。

为了提供快速的图像处理，我们将 Windows 窗体中的标准类Bitmap的对象转换为我们的类CImage的对象，它们的方法很快，因为它们直接访问像素集，而使用Bitmap的标准方式包括实现相对较慢的方法GetPixel和SetPixel或使用LockBits的方法，它们很快，但不能用于索引图像。

类CImage相当简单:它包含图像的属性width、height和nBits以及实际项目中使用的方法。我的方法在专门讨论项目的章节中有描述。下面是我们这个阶层CImage的定义。

  class CImage
  { public Byte[] Grid;
     public int width, height, nBits;

     public CImage() { } // default constructor

     public CImage(int nx, int ny, int nbits) // constructor
     {
        width = nx;
        height = ny;
        nBits = nbits;
        Grid = new byte[width * height * (nBits / 8)];
    }

    public CImage(int nx, int ny, int nbits, byte[] img) // constructor
    {
        width = nx;
        height = ny;
        nBits = nbits;
        Grid = new byte[width * height * (nBits / 8)];
        for (int i = 0; i < width * height * nBits / 8; i++) Grid[i] = img[i];
    }
  } //*********************** end of class CImage *****************

项目描述中描述了类CImage的方法。

二、降噪

数字图像经常因随机误差而失真，随机误差通常被称为*噪声。*主要有两种噪声:高斯噪声和脉冲噪声(见图 2-1 )。高斯噪声是具有类似于高斯分布的概率分布的统计噪声。它主要出现在采集期间(例如，在传感器中)。这可能是由于照明不足或传感器温度过高造成的。它来自许多自然来源，如导体中原子的热振动，称为热噪声。它影响图像的所有像素。

脉冲噪声也称为椒盐噪声，表现为稀疏出现的亮像素和暗像素。它来自于脉冲失真，比如那些来自于电子设备附近的电焊所产生的失真，或者是由于旧照片的不正确存储所导致的失真。它影响像素组的一小部分。

img/474294_1_En_2_Fig1_HTML.jpg

图 2-1

噪声的例子:(a)高斯噪声；脉冲噪声

最简单的过滤器

我们首先考虑降低高斯噪声强度的方法。本章稍后将讨论降低脉冲噪声强度的算法。

降低高斯噪声强度的最有效的方法是用在 P 附近的一小组像素的亮度的平均值来代替像素 P 的亮度。这种方法基于随机值理论的事实:N 个均匀分布的随机值的平均值的标准偏差比单个值的标准偏差小√ N 倍。这种方法使用遮罩对图像进行二维卷积，遮罩是排列在一个 W × W 像素的正方形中的权重数组，实际像素 P 位于正方形的中间。本章稍后将介绍该过滤器的源代码。这种方法有两个缺点:它相当慢，因为它对图像的每个像素使用 W × W 加法，并且它使图像模糊。它将近似同质区域边界处的精细边缘转换为斜坡，其中亮度在宽度等于 W 像素的条纹中线性变化。快速均值滤波器克服了第一个缺点。然而，第二个缺点非常重要，它妨碍了均值滤波器用于消除噪声。然而，平均滤波器对于改善具有阴影的图像(即，表示非均匀照明物体的图像)是重要的，我们将在后面看到。我建议另一个过滤器的目的是消除噪音。

首先，让我们描述一下最简单的均值滤波器。我给出了这个简单方法的源代码，读者可以在自己的程序中使用它。在这段代码以及许多其他代码示例中，我们使用了某些类，这些类将在下一节中定义。

最简单的均值滤波器

非加权平均滤波器计算一个滑动方形窗口中的平均灰度值，该窗口为 W × W 像素，其中 W 是方形滑动窗口的宽度和高度。窗口大小 W 越大，对高斯噪声的抑制越强:滤波器以因子 W 降低噪声。为了对称，数值 W 通常为奇数 W = 2 × h + 1。窗口内像素的坐标( x + xx，y + yy )以当前中心像素( x，y )为中心对称变化，间隔:h≤xx≤+h和h≤YY≤+h**h= 1，2，3

在图像边界附近，窗口部分位于图像之外。在这种情况下，计算失去了其自然的对称性，因为只有图像内部的像素可以被平均。解决边界问题的合理方法是控制坐标( x + xx，y + yy )，如果它们指向图像之外。如果这些坐标在图像之外，灰度值的求和必须暂停，并且除数 nS 不应增加。

这里给出了一个颜色平均算法的例子。在我们的代码中，我们经常使用由某些符号行表示的注释:标有=的行标记循环的开始和结束，标有减号的行标记if指令，等等。这使得代码的结构更加清晰可见。

该算法最简单的慢速版本有四个嵌套的for循环。

public void CImage::Averaging(CImage input, int HalfWidth)
{ int nS, sum;
  for (int y=0; y<height; y++) //================================
  { for (int x=0; x<width; x++) //==============================
    { nS=sum=0;
      for (int yy=-HalfWidth; yy<=HalfWidth; yy++) //=======
      {  if (y+yy>=0 && y+yy<input.height )
          for (int xx=-HalfWidth; xx<=HalfWidth; xx++) //===
          { if (x+xx>=0 && x+xx<input.width)
            { sum+=input.Grid[x+xx+width*(y+yy)];
               nS++;
            }
          } //====== end for (xx... ===================
      } //======== end for (yy... =======================
      Grid[x+width*y]=(sum+nS/2)/nS; //+nS/2 is for rounding
    } //============ end for (x... ========================
  } //============== end for (y... ========================
} //****************** end Averaging ****************************

这是源代码:读者可以复制它并放入它的 C#源代码中，它就可以工作了。

参数HalfWidth是滑动窗口宽度的一半。窗口的宽度和高度都等于2*HalfWidth+1。前面代码中的变量x和y是Grid,中像素的索引，xx和yy是滑动平均窗口内像素的索引。

最内层for循环中sum的计算需要 W × W 加法和 W × W 对输入图像的每个像素进行图像访问，相当耗时。

让我们再一次注意到，虽然平均滤波器在降低高斯噪声的强度方面非常有效，但是它们强烈地模糊了图像。因此，它们不应用于降噪。为此，我建议使用本章稍后描述的 sigma 滤波器。平均用于阴影校正(参见第四章)。因此，它将主要用于相当大的滑动窗口，大到图像宽度的一半。那么最简单的平均程序变得非常耗时，以致于实际上不可能使用它。例如，在 1000 × 1000 像素的灰度图像和 400 × 400 像素的滑动窗口的情况下，这对于阴影校正是典型的，函数Averaging在标准 PC 上的运行时间可能需要大约 20 分钟。

快速平均滤波器

最简单的平均滤波器对图像的每个像素进行加法运算和除法运算。例如，在滑动窗口为 5 × 5 = 25 个像素的情况下，它对每个像素进行 5 5 + 1 = 26 次操作。有些应用(例如阴影校正)使用更大的滑动窗口；比如 400 × 400 的一个= 16 万像素。如果应用使用最简单的均值滤波器和如此大的滑动窗口，它可以运行几分钟。

可以使用以下基本思想来加速平均滤波器:可以首先计算 1 × W 像素的小一维窗口中的灰度值之和。我们称这些总和的数组为*列。*在以下对基本思想的描述中，我们使用笛卡尔坐标系，图像的列的索引是横坐标 x ，行的索引是纵坐标 y 。我们认为纵坐标轴是向下的，从上到下，这在图像处理中是常见的，并且不同于数学中纵坐标轴的方向。

快速均值滤波器解决的问题与最简单的滤波器相同，但速度要快得多。快速滤镜以与最简单滤镜相同的方式模糊图像。因此，它的主要应用是阴影校正，而不是噪声抑制。

在使用刚才提到的基本思想时，可以将每个像素的运算次数从 W × W 减少到≈4:快速滤波器计算并保存每列 W 像素的灰度值之和，每列的中间像素位于图像的实际行中(图 2-2 )。然后，滤波器直接计算其中心像素位于一行开始处的窗口上的总和；也就是说，将各列中保存的总数相加。然后，窗口沿该行移动一个像素，筛选器通过将窗口右边框的列总和的值相加，并减去左边框的列总和的值，来计算下一个位置的总和。需要检查要增加或减少的列是否在图像内部。如果不是，则必须跳过相应的加法或减法。

img/474294_1_En_2_Fig2_HTML.png

图 2-2

快速平均滤波器的功能说明

由于对列和的计算应用了类似的过程，每个像素的平均加法或减法次数减少到≈2 + 2 = 4。窗口内的总和必须仅针对每行开头的像素直接计算(即，通过添加HalfWidth + 1 列总和)。必须仅对图像第一行的像素直接计算列的总和。

当前进到图像的下一行时，过滤器通过增加每列下端以下的灰度值并减去每列上端的灰度值来更新列的值(图 2-2 )。在这种情况下，还需要检查要加或减的灰度值是否在图像中。一旦计算出窗口中灰度值的总和，该滤波器就将该总和除以(舍入)窗口与图像相交处的像素数，并将结果保存在输出图像的相应像素中。

这是为过滤灰度图像而设计的快速平均滤波器的最简单版本的源代码。

 public int FastAverageM(CImage Inp, int hWind, Form1 fm1)
// Filters the grayscale image "Inp" and returns the result in 'Grid' of the
// calling image.
{
    width  = Inp.width ; height = Inp.height; // elements of class CImage
    Grid = new byte[width  * height];
    int[] SumColmn; int[] nPixColmn;
    SumColmn = new int[width];
    nPixColmn = new int[width];
    for (int i = 0; i < width ; i++) SumColmn[i] = nPixColmn[i] = 0;

    int nPixWind = 0, SumWind = 0;
    for (int y = 0; y < height + hWind; y++) //=============================
    {
      int yout = y - hWind, ysub = y - 2 * hWind - 1;
      SumWind = 0; nPixWind = 0;
      int y1 = 1 + (height + hWind) / 100;
      for (int x = 0; x < width + hWind; x++) //===========================
      {
        int xout = x - hWind, xsub = x - 2 * hWind - 1;    // 1\. and 2\. addition
        if (y < height && x < width )
        {
          SumColmn[x] += Inp.Grid[x + width  * y];
          nPixColmn[x]++;     // 3\. and 4\. addition
        }

        if (ysub >= 0 && x < width )
        {
          SumColmn[x] -= Inp.Grid[x + width * ysub];
          nPixColmn[x]--;
        }
        if (yout >= 0 && x < width )
        {
          SumWind += SumColmn[x];
          nPixWind += nPixColmn[x];
        }
        if (yout >= 0 && xsub >= 0)
        {
          SumWind -= SumColmn[xsub];
          nPixWind -= nPixColmn[xsub];
        }
        if (xout >= 0 && yout >= 0)
          Grid[xout + width  * yout] = (byte)((SumWind + nPixWind / 2) / nPixWind);
      } //===================== end for (int x = 0; =====================
    } //====================== end for (int y = 0; ======================
    return 1;

    } //************************* end FastAverageM ***********************

接下来我将展示为彩色和灰度图像设计的快速平均滤波器的通用源代码。它使用变量int nbyte，对于彩色图像设置为 3，对于灰度图像设置为 1。我们为(2*hWind + 1) 2 像素的滑动窗口中的颜色强度之和定义了一个三元素的数组SumWind[3]，用于红色、绿色和蓝色强度之和。在灰度图像的情况下，只有元素SumWind[0]被使用。我们使用下面描述的变量。

具有坐标(c+nbyte*x, nbyte*y)的位置是颜色通道的位置，红色、绿色或蓝色通道之一，其强度被添加到数组SumColmn的相应元素。位置(c+nbyte*x, nbyte*ysub)是一个颜色通道的位置，其强度将从SumColmn中减去。变量c+nbyte*x是短列的横坐标，其内容将被添加到SumWind[c]。变量c+nbyte*xsub是短列的横坐标，其内容将从SumWind[c]中减去。

public int FastAverageUni(CImage Inp, int hWind, Form1 fm1)
// Filters the color or grayscale image "Inp" and returns the result in
// 'Grid' of calling image.
{ int c = 0, nByte = 0;
   if (Inp.N_Bits == 8) nByte = 1;
   else nByte = 3;
   width  = Inp.width ; height = Inp.height; // elements of the class "Cimage"
   Grid = new byte[nByte * width  * height];
   int[] nPixColmn;
   nPixColmn = new int[width];
   for (int i = 0; i < width ; i++) nPixColmn[i] = 0;
   int[,] SumColmn;
   SumColmn = new int[width, 3];
   int nPixWind = 0, xout = 0, xsub = 0;
   int[] SumWind = new int[3];
   for (int y = 0; y < height + hWind; y++) //=============================
  { int yout = y - hWind, ysub = y - 2 * hWind - 1;
     nPixWind = 0;
     for (c = 0; c < nByte; c++) SumWind[c] = 0;
     int y1 = 1 + (height + hWind) / 100;
     for (int x = 0; x < width  + hWind; x++) //============================
     { xout = x - hWind;
        xsub = x - 2 * hWind - 1;    // 1\. and 2\. addition
        if (y < height && x < width) // 3\. and 4\. addition
       { for (c=0; c< nByte; c++)
            SumColmn[x, c] += Inp.Grid[c + nByte*(x + width*y)];
         nPixColmn[x]++;
      }
      if (ysub >= 0 && x < width )
      { for (c=0; c<nByte; c++)
            SumColmn[x, c] -=Inp.Grid[c+nByte*(x+ width*ysub)];
        nPixColmn[x]--;
      }
      if (yout >= 0 && x < width )
      { for (c = 0; c < nByte; c++) SumWind[c] += SumColmn[x, c];
         nPixWind += nPixColmn[x];
      }
      if (yout >= 0 && xsub >= 0)
      { for (c = 0; c < nByte; c++) SumWind[c] -= SumColmn[xsub, c];
         nPixWind -= nPixColmn[xsub];
      }
      if (xout >= 0 && yout >= 0)
        for (c = 0; c < nByte; c++)
          Grid[c+nByte*(xout+width*yout)]=(byte)( SumWind[c] / nPixWind);
    } //============= end for (int x = 0;  =============================
  } //============== end for (int y = 0;  ==============================
  return 1;
} //***************** end FastAverageUni ********************************

该源代码可以在相应的 Windows 窗体项目中使用。它不是最快的版本；通过去除内部循环中的乘法运算，速度可以提高 50%。一些乘法可以在开始循环之前执行；其他一些可以通过添加来代替。更快的版本可以包含以下九个循环，而不是在FastAverageM或FastAverageUni中具有索引 y 和 x 的两个循环:

for (int yOut=0; yOut<=hWind; yOut++)
{ for(int xOut=0; xOut<=hWind; xOut++){...} //Loop 1
   for(xOut=hWind+1; xOut<=width-hWind-1; xOut++) {...} //Loop 2
   for(xOut=width-hWind; xOut<width; xOut++) {...} //Loop 3
}
for(yOut=hWind+1; yOut<=height-hWind-1; yOut++)
{ for(xOut=0; xOut<=hWind; xOut++) {...} //Loop 4
   for(xOut=hWind+1; xOut<=width-hWind-1; xOut++) {...} //Loop 5
   for(xOut=width-hWind; xOut<width; xOut++) {...} //Loop 6
}
for(yOut=height-hWind; yOut<height; yOut++)
{ for(xOut=0; xOut<=hWind; xOut++) {...} /Loop 7
   for(xOut=hWind+1; xOut<=width-hWind-1; xOut++) {...} //Loop 8
   for(xOut=width-hWind; xOut<width; xOut++) {...} //Loop 9
}

九个循环中的每一个循环都处理图像的一部分(见图 2-3 )，该部分要么是hWind + 1 像素宽，要么是hWind + 1 像素高。

img/474294_1_En_2_Fig3_HTML.png

图 2-3

图像的九个部分对应于九个循环

仅当满足条件hWind ≤ min(宽度，高度)/2 - 1 时，才能使用此版本的快速平均滤波器。在这个版本的例程中，变量为xOut的内部循环不包含乘法运算和if指令。该例程比之前描述的FastAverageM快 60%。然而，它要长得多，调试起来也困难得多。速度的提高并不重要:这段代码用 0.7 秒处理一个 2448 × 3264 像素的大彩色图像，滑动窗口为 1200 × 1200 像素，而FastAverageM需要 1.16 秒。这些计算时间几乎与滑动窗口的大小无关:对于 5 × 5 像素的滑动窗口，它们相应地为 0.68 和 1.12 秒。

快速高斯滤波器

平均滤波器产生平滑的图像，其中可以看到原始图像中不存在的一些矩形形状。这些形状的出现是因为平均滤镜将每个亮像素转换为滑动窗口大小的均匀亮矩形。一旦一个像素具有明显不同于相邻像素值的亮度，矩形就变得可见。这是不必要的失真。当使用高斯滤波器时，可以避免这种情况，该高斯滤波器将要相加的灰度值乘以根据二维高斯定律随着距窗口中心的距离而衰减的值。此外，高斯滤波器能更好地抑制噪声。重量示例如图 2-4 所示。

img/474294_1_En_2_Fig4_HTML.png

图 2-4

经典高斯滤波器滑动窗口中的权重示例

这些值称为过滤器的权重。对应于二维高斯定律的权重是小于 1 的浮点数:

w(、和=(2π】和

**它们可以预先计算出来，保存在一个二维数组中，数组的大小与滑动窗口的大小相对应(图 2-4 )。然后将滤波图像的灰度值或颜色通道乘以权重，并计算乘积的和。这个过程需要对要滤波的灰度图像的每个像素进行W2浮点乘法和 W 2 加法，其中 W 是滑动窗口的宽度。在彩色图像的情况下，这个数字是 3W2。

使用统计学知识，即许多等效概率分布的卷积趋向于高斯分布，有可能获得大致相同的结果。这一过程的收敛如此之快，以至于仅计算三个矩形分布的卷积就足以获得良好的近似。矩形分布在区间内密度恒定，在区间外密度为零。如果一个图像用一个滤波器处理三次，其结果相当于用三个矩形权重的卷积进行加权滤波。因此，为了近似执行图像的高斯滤波，用快速平均滤波器对图像滤波三次就足够了。该过程要求独立于窗口大小的像素的每个颜色通道进行 4 × 3 = 12 次整数加法。我们已经计算出等效高斯分布的标准与平均滤波器的滑动窗口的半宽度成比例。在表 2-1 hWind中，是平均窗口的半宽度，Sigma是随机变量的标准，其分布对应于通过快速滤波器的三重滤波计算的权重。

表 2-1。随着 hWind 的增加，适马/hWind 的关系趋于 1

|

hWind

|

one

|

Two

|

three

|

four

|

five

| | --- | --- | --- | --- | --- | --- | | Sigma | One point four one four | Two point four five six | Three point four five nine | Four point four five six | Five point four nine | | Sigma/hWind | One point four one four | One point two one eight | One point one five three | One point one one four | One point zero nine eight |

可以看到，当窗口宽度增加时，关系Sigma/hWind趋于 1。

图 2-5 显示了近似高斯滤波器的权重与真实高斯权重的差异。

img/474294_1_En_2_Fig5_HTML.png

图 2-5

近似高斯滤波器和真实高斯滤波器的权重

中值滤波器

平均和高斯滤波器可最有效地抑制高斯噪声。具有宽度为 W = 2 h + 1 像素的滑动窗口的平均滤波器将同质区域的陡峭边缘转换为宽度为 W 的斜坡。在高斯滤波器的情况下，斜坡更陡。但是，这两种滤镜都会模糊图像，因此不应用于噪声抑制。

大多数关于图像处理的教科书都推荐使用中值滤波器来抑制噪声。中值滤波器对滑动窗口中的颜色强度进行排序，并用位于排序序列中间的强度替换滑动窗口中间的强度。中值滤波器也可用于抑制脉冲噪声或椒盐噪声。

然而，几乎没有哪本教科书能让读者注意到中值滤波器的重大缺陷。首先，它严重扭曲了形象。具有(2 * h + 1) 2 像素的滑动窗口的中值滤波器删除宽度小于 h 像素的每个条纹。它还删除了矩形每个角上大约 2 个 h 像素的三角形部分。更重要的是，如果条纹之间的空间宽度也等于 h ，它会反转包含宽度为 h 的一些平行条纹的图像的一部分(比较图 2-6 和图 2-7 )。如果读者注意到 median 根据多数做出决策，这就很容易理解了:如果滑动窗口中的大多数像素都是暗的，那么中心像素就会变暗。

img/474294_1_En_2_Fig7_HTML.png

图 2-7

中值为 5 × 5 像素的滤波后的相同图像

img/474294_1_En_2_Fig6_HTML.png

图 2-6

原始图像和 5 × 5 像素的滑动窗口

也不推荐使用中值来抑制脉冲噪声，因为这将删除与噪声无关的细线形状的对象。我在后面的章节中建议了一个有效的方法。

适马滤波器:最有效的一种

西格玛滤波器以与平均滤波器相同的方式减少噪声:通过平均许多灰度值或颜色。西格玛滤波器的思想是仅对滑动窗口中与中心像素的强度相差不超过固定参数容差的那些强度(即颜色通道的灰度值或强度)进行平均。根据这一思想，西格玛滤波器减少了高斯噪声，并保留了图像中的边缘不模糊。

西格玛过滤器是由 John-Sen Lee (1983)提出的。然而，直到最近，它仍然鲜为人知:据我所知，没有任何图像处理教科书提到过它。在一篇专业论文中只提到过一次 Chochia (1984)。

托马西和曼杜奇(1998)提出了一种类似于西格玛滤波器的滤波器，他们称之为双边滤波器 *。*他们建议为被平均的颜色分配两种权重:域权重随着平均像素与滑动窗口的中心像素的距离增加而变小，范围权重随着被平均的像素和中心像素的颜色强度之间的差异增加而变小。两个权重都可以定义为高斯分布的密度。该滤波器工作良好:它减少了高斯噪声，并保留了边缘的清晰度。然而，它本质上比 sigma 滤波器慢；例如，处理 2500 × 3500 像素的彩色图像，双边滤波器需要 30 秒，而最简单的 sigma 滤波器只需要 7 秒。因此，双边滤波器大约比 sigma 滤波器慢四倍。双边过滤器的作者在参考文献中没有提到 sigma 过滤器。

为了解释为什么 sigma 滤波器工作得如此好，它可以被视为众所周知的期望最大化(EM)算法的第一次迭代。将滑动窗口中的颜色细分成两个子集，接近中心像素颜色的子集和远离中心像素颜色的子集，可以被认为是期望步骤，而接近颜色的平均是最大化步骤。众所周知，EM 算法收敛相当快。因此，单次第一次迭代带来的结果接近滑动窗口中像素子集的平均值，该像素子集的颜色很可能接近中间像素的颜色。

当比较 sigma 滤波器和双边滤波器时，可以注意到 sigma 滤波器使用类似于双边滤波器的算法，其中高斯分布被可以更快计算的更简单的矩形分布代替。

让我们首先展示使用 sigma 过滤器过滤灰度图像的伪代码。我们用Input表示输入图像，用Output表示输出图像。让一个坐标为(X, Y)的像素滑过输入图像。设M为(X, Y)处的灰度值。对每个位置(X, Y)做:

sum=0; number=0; M=Input(X,Y);
Let another pixel (x, y) run through the gliding window with the center of window at (X, Y).
for each pixel Input(x,y) do:
if (abs(Input(x,y) - M) < tolerance)
{ sum+=Input(x,y);
   number++;
}
Output(X,Y)=Round(sum / number);

现在让我们展示这个简单解决方案的源代码。我们用加法代替了经常重复的乘法，使这个方法更快(13%)。

public int SigmaSimpleColor (CImage Inp, int hWind, int Toleranz)
// The sigma filter for 3 byte color images.
{
   int[] gvMin = new int[3], gvMax = new int[3], nPixel = new int[3];
   int [] Sum = new int[3];
   int c, hWind3 = hWind * 3, NX3 = width * 3, x3, xyOut;
   N_Bits = Inp.N_Bits; // "N_Bits is an element of the class CImage
  for (int y = xyOut = 0; y < height; y++) // ==============================
  {
    int y1, yStart = Math.Max(y - hWind, 0) * NX3,
              yEnd = Math.Min(y + hWind, height - 1) * NX3;
    for (int x = x3 = 0; x < width; x++, x3+=3, xyOut+=3) //=================
    {
       int x1, xStart = Math.Max(x3 - hWind3, 0),
                 xEnd = Math.Min(x3 + hWind3, NX3 - 3);
      for (c=0; c<3; c++)
      {
         Sum[c] = 0;
         nPixel[c] = 0;
        gvMin[c] = Math.Max(0, Inp.Grid[c + xyOut] - Toleranz);
        gvMax[c] = Math.Min(255, Inp.Grid[c + xyOut]+Toleranz);
      }
      for (y1 = yStart; y1 <= yEnd; y1 += NX3)
        for (x1 = xStart; x1 <= xEnd; x1 += 3)
          for (c = 0; c < 3; c++)
          {
             if (Inp.Grid[c+x1+y1]>=gvMin[c] && Inp.Grid[c+x1+y1]<=gvMax[c])
            { Sum[c]+=Inp.Grid[c+x1+y1];
               nPixel[c]++;
            }
         }
      for (c = 0; c < 3; c++) //======================================
      { if (nPixel[c] > 0) Grid[c+xyOut] = (byte)((Sum[c] + nPixel[c]/2)/nPixel[c]);
        else Grid[c+xyOut]=0;
      } //================= end for (c... ===========================
    } //================== end for (int x... ==========================
  } //=================== end for (int y... ===========================
  return 1;
} //********************** end SigmaSimpleColor *************************

这种解决方案工作得很好，但是如果滑动窗口的大小大于 5 × 5 像素，则它相当慢:对于灰度图像，它需要大约每像素 OPP = 4 **W2 操作，或者对于彩色图像，它需要 9 **W2 操作。在大多数实际情况下，使用窗口大小为 3 × 3 或 5 × 5 像素的 sigma 滤波器就足够了。因此SigmaSimpleColor几乎可以在任何地方使用。

稍后我们将看到更快版本的 sigma 滤波器。不幸的是，在这种情况下不可能应用快速平均滤波器中使用的方法，因为该过程是非线性的。由于使用了局部直方图，该过程可以变得更快。直方图是一个数组，其中每个元素包含窗口中相应灰度值或颜色强度的出现次数。西格玛滤波器通过更新过程计算窗口每个位置的直方图:窗口右边界的垂直列中的灰度值或颜色强度用于增加直方图的相应值，而左边界的值用于减少它们。

设 OPP 是每个像素的运算次数。23 W 是实现直方图所需的运算次数，23(2tol+1)是计算直方图的 3(2tol+1)个值和相应的像素数所需的运算次数。这样总的 OPP = 23* W +23(2* tol +1)。

这是西格玛过滤器的源代码，带有彩色图像的局部直方图。

public int SigmaColor(CImage Inp, int hWind, int Toleranz, Form1 fm1)
// The sigma filter for color images with 3 bytes per pixel.
{
  int gv, y1, yEnd, yStart;
  int[] gvMin = new int[3], gvMax = new int[3], nPixel = new int[3];
  int Sum = new int[3];
  int[,] hist = new int[256, 3];
  int c;
  for (y1 = 0; y1 < 256; y1++) for (c = 0; c < 3; c++) hist[y1, c] = 0;
  N_Bits = Inp.N_Bits;
  fm1.progressBar1.Value = 0;
  int yy = 5 + nLoop * height / denomProg;
  for (int y = 0; y < height; y++) //=====================================
  {
    if ((y % yy) == 0) fm1.progressBar1.PerformStep();
    yStart = Math.Max(y - hWind, 0);
    yEnd = Math.Min(y + hWind, height - 1);
    for (int x = 0; x < width; x++) //====================================
    {
      for (c = 0; c < 3; c++)  //========================================
      {
        if (x == 0) //-------------------------------------------
        {
          for (gv = 0; gv < 256; gv++) hist[gv,c] = 0;
          for (y1 = yStart; y1 <= yEnd; y1++)
            for (int xx = 0; xx <= hWind; xx++)
              hist[Inp.Grid[c + 3 * (xx + y1 * width)], c]++;
        }
        else
        {
          int x1 = x + hWind, x2 = x - hWind - 1;
          if (x1 < width - 1)
            for (y1 = yStart; y1 <= yEnd; y1++)
               hist[Inp.Grid[c + 3 * (x1 + y1 * width)], c]++;
          if (x2 >= 0)
            for (y1 = yStart; y1 <= yEnd; y1++)
            {
              hist[Inp.Grid[c + 3 * (x2 + y1 * width)], c]--;
              if (hist[Inp.Grid[c + 3 * (x2 + y1 * width)], c] < 0) return -1;
            }
        } //---------------- end if (x==0) ------------------------

        Sum[c] = 0; nPixel[c] = 0;
        gvMin[c] = Math.Max(0, Inp.Grid[c + 3 * x + 3 * width * y] - Toleranz);
        gvMax[c] = Math.Min(255, Inp.Grid[c + 3 * x + 3 * width * y] + Toleranz);
        for (gv = gvMin[c]; gv <= gvMax[c]; gv++)
        {
          Sum[c] += gv * hist[gv, c]; nPixel[c] += hist[gv, c];

        }
        if (nPixel[c] > 0)
          Grid[c + 3 * x + 3 * width * y] = (byte)((Sum[c] + nPixel[c] / 2) / nPixel[c]);
        else Grid[c + 3 * x + 3 * width * y] = Inp.Grid[c + 3 * x + 3 * width * y];
      } //================ end for (c... ============================
    } //================= end for (int x... ===========================
  } //================== end for (int y... ============================
  return 1;
} //********************** end SigmaColor ******************************

如果滑动窗口的宽度小于 7(hWind小于 3)，方法SigmaSimpleColor比SigmaColor快。宽度值越大SigmaColor越快。SigmaColor的工作时间随滑动窗口的宽度变化相当缓慢。参见表 2-2 。

表 2-2

1200 × 1600 像素彩色图像的工作时间(秒)

|   |

hWind

| | --- | --- | | 方法 | one | Two | three | six | | SigmaSimpleColor | Zero point seven one | One point six eight | Three point one one | Ten point two six | | SigmaColor | Two point four three | Two point five five | Two point six three | Three point zero two |

本书中描述的几乎所有项目都使用了方法SigmaSimpleColor和灰度图像的类似方法。

脉冲噪声的抑制

正如已经指出的，脉冲噪声影响单个的、偶尔选择的像素或相邻像素的小组，而不是图像的所有像素。后者是高斯噪声的特征。我们区分暗脉冲噪声和亮脉冲噪声。暗噪声包含亮度低于其邻域亮度的像素或像素组，而在亮噪声的情况下，像素的亮度高于其邻域的亮度。

抑制灰度或彩色图像中的脉冲噪声的问题相当复杂:在暗噪声的情况下，需要自动检测满足以下条件的像素的所有子集 S :

  1. 子集 S 的所有像素必须具有低于或等于阈值 T 的亮度。

  2. 子集 S 被连接。

  3. 像素数(即 S 的面积)低于预定值 M 。

  4. 不属于 S 但是与 S 的像素相邻的所有像素必须具有高于 T 的亮度。

  5. 对于不同的子集 S ，阈值 T 可以不同。

这个问题很难，因为阈值 T 是未知的。理论上，在一个接一个地测试 T 的所有可能值的同时解决问题是可能的。然而，这样的算法会非常慢。

米（meter 的缩写））I. Schlesinger (1997 年)提出了以下快速解决方案的想法:

我们提出了这样一个过程，即在二进制图像中，利用后续阈值处理的“噪声去除”给出了与利用后续常用噪声去除的阈值处理相同的结果。这种等价适用于每个阈值的值。

项目PulseNoise的算法由四个步骤组成:(1)按照像素亮度的递增顺序对像素进行排序，(2)去除白噪声，(3)按照像素亮度的递减顺序对像素进行排序，以及(4)去除黑噪声。

使用长度等于不同亮度数量的附加阵列，在像素的双重扫描期间完成像素的排序。

在有序像素的单倍扫描期间进行白噪声去除。设 t i 为亮度为 v 的像素(t i )。像素的处理包括检查是否存在包含该像素的白噪声区域 G。最初该区域只包含 t i 。然后该区域增长，使得如果(1)t’是区域 G 的邻居，则某个像素 t’被包括在其中；以及(2)v(t ')≥v(tI)。区域增长，直到满足以下条件之一:

1。所获得的区域 G 的尺寸超过了预定尺寸。在这种情况下，区域 G 不被认为是噪声，并且下一个像素 t (i+1) 被处理。

2。没有像素可以包含在 G 中。在这种情况下，G 是一个噪声，其所有像素的亮度等于相邻像素的最大亮度。

去除黑噪声类似于去除白噪声。"

下面是我们在 Windows 窗体项目PulseNoiseWF中实现该算法的详细描述和源代码。我们做了一些无关紧要的改变:我们使用亮度直方图histo[256]并定义了一个二维数组Index[256][*]，而不是根据亮度对像素进行排序，这太慢了。该阵列包含每个亮度值light的不同数量的像素索引(x + Width*y)的Index[light][histo[light]]。值(x + Width*y)是坐标为(x，y)的像素的索引，通过该索引可以在图像中找到该像素。因此，例如，值Index[light][10]是亮度等于light的所有像素的集合中的第十个像素的索引。当从小亮度开始读取数组Index时，我们获得以亮度增加的顺序排序的像素的索引，但是当以相反的方向读取时，我们获得以亮度减少的顺序排序的像素的索引。

该项目包含一个文件CImage.cs，该文件定义了包含用于预处理图像的方法的类CImage。这些方法处理CImage类的对象，这比在 Windows 窗体中处理位图更容易也更快。

该项目还包含文件CPnoise.cs，其类CPnoise包含实现脉冲噪声抑制的方法。还有一个包含类Queue的小文件Queue.cs，它定义了实现众所周知的结构Queue的简单方法。队列是先进先出(FIFO)的数据结构。在 FIFO 数据结构中，添加到队列中的第一个元素将是第一个被移除的元素。这相当于这样的要求:一旦添加了新元素，在删除新元素之前，必须删除在它之前添加的所有元素。我在下面描述了这个项目的方法。

我们一起考虑灰度和彩色图像的情况，因为这两种情况之间的唯一区别是，在彩色图像的情况下，像素的三个通道(红色、绿色和蓝色)的共同亮度被认为与灰度图像的灰度值相同。我们使用第三章中描述的方法MaxC(byte R, byte G, byte B)来计算彩色像素的亮度。在这个项目中，我们只在给找到的集合 S 的像素分配颜色的指令中直接使用颜色信息。

项目PulseNoiseWF显示如图 2-8 所示的表格。当用户单击“打开图像”时，“打开文件对话 1”对话框打开，用户可以选择图像。图像打开，并根据打开图像的像素格式，通过方法BitmapToGrid或BitmapToGridOld转换成对象CImage.Orig。BitmapToGrid使用BitmapData类的方法LockBits，该方法允许直接评估位图的像素；这比使用GetPixel的标准评估要快得多。这里显示的是BitmapToGrid的代码。

img/474294_1_En_2_Fig8_HTML.jpg

图 2-8

项目的形式PulseNoiseWF

private void BitmapToGrid(Bitmap bmp, byte[] Grid)
{
  Rectangle rect = new Rectangle(0, 0, bmp.Width, bmp.Height);
  BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadWrite,
                                                         bmp.PixelFormat);

  int nbyte;
  switch (bmp.PixelFormat)
  {
    case PixelFormat.Format24bppRgb:        nbyte = 3; break;
    case PixelFormat.Format8bppIndexed:    nbyte = 1; break;
    default: MessageBox.Show("Not suitable pixel format=" + bmp.PixelFormat);
                  return;
  }
  IntPtr ptr = bmpData.Scan0;
  int length = Math.Abs(bmpData.Stride) * bmp.Height;
  byte[] rgbValues = new byte[length];
  System.Runtime.InteropServices.Marshal.Copy(ptr, rgbValues, 0, length);

  progressBar1.Visible = true;
  for (int y = 0; y < bmp.Height; y++)
  {
    int y1 = 1 + bmp.Height / 100;
    if (y % y1 == 1) progressBar1.PerformStep();

    for (int x = 0; x < bmp.Width; x++)
    {
      if (nbyte == 1)  // nbyte is global according to the PixelFormat of "bmp"
      {

        Color color = bmp.Palette.Entries[rgbValues[x+Math.Abs(bmpData.Stride)*y]];
        Grid[3 * (x + bmp.Width * y) + 0] = color.B;
        Grid[3 * (x + bmp.Width * y) + 1] = color.G;
        Grid[3 * (x + bmp.Width * y) + 2] = color.R;
      }
      else
        for (int c = 0; c < nbyte; c++)
           Grid[c + nbyte * (x + bmp.Width * y)] =
                           rgbValues[c + nbyte * x + Math.Abs(bmpData.Stride) * y];
    }
  }
  bmp.UnlockBits(bmpData);
  progressBar1.Visible = false;
} //*********************** end BitmapToGrid ********************

方法BitmapToGridOld使用GetPixel。在 8 位图像的情况下，必须使用颜色表Palette来定义每个像素的颜色。这在GetPixel中比在Palette.Entries中进行得更快，因为在BitmapToGrid中就是这种情况。因此，我们在 8 位原始图像的情况下使用BitmapToGridOld。

表格Form1包含以下定义:

    private Bitmap origBmp;
    private Bitmap Result; // result of processing
    CImage Orig;  // copy of original image
    CImage Work;  // work image
    public Point[] v = new Point[20];
    // "v" are corners of excluded rectangles
    int Number, // number of defined elements "v"
      maxNumber = 8;
    bool Drawn = false;

图像显示在表单左侧的pictureBox1中。用户现在可以改变要消除的暗斑和亮斑的最大尺寸(像素数)的值。这可以用带有标签Delete dark和Delete light的工具numericUpDown来完成。

当用户点击脉冲噪声时，会启动相应的程序残端。这是树桩的代码。

private void button3_Click(object sender, EventArgs e)     // Impulse noise
{
  Work.Copy(Orig);  // "Work" is the work image defined outside
  int nbyte = 3;
  Drawn = true;
  Work.DeleteBit0(nbyte);
  int maxLight, minLight;
  int[] histo = new int[256];
  for (int i = 0; i < 256; i++) histo[i] = 0;
  int light, index;
  int y1 = Work.height / 100;
  for (int y = 0; y < Work.height; y++) //==========================
  {
    if (y % y1 == y1 - 1) progressBar1.PerformStep();
    for (int x = 0; x < Work.width; x++) //=======================
    {
      index = x + y * Work.width; // Index of the pixel (x, y)
      light = MaxC(Work.Grid[3 * index+2] & 254,
                    Work.Grid[3 * index + 1] & 254,
                    Work.Grid[3 * index + 0] & 254);
      if (light < 0) light = 0;
      if (light > 255) light = 255;
      histo[light]++;
    } //=============== end for (int x=1; .. ===================
  } //================ end for (int y=1; .. =====================
  progressBar1.Visible = false;
  for (maxLight = 255; maxLight > 0; maxLight--) if (histo[maxLight] != 0) break;
  for (minLight = 0; minLight < 256; minLight++) if (histo[minLight] != 0) break;
  CPnoise PN = new CPnoise(histo, 1000, 4000);
  PN.Sort(Work, histo, Number, pictureBox1.Width, pictureBox1.Height, this);
  progressBar1.Visible = false;
  int maxSizeD = 0;
  if (textBox1.Text != "") maxSizeD = int.Parse(textBox1.Text);
  int maxSizeL = 0;
  if (textBox2.Text != "") maxSizeL = int.Parse(textBox2.Text);
  PN.DarkNoise(ref Work, minLight, maxLight, maxSizeD, this);
  progressBar1.Visible = false;
  Work.DeleteBit0(nbyte);
  PN.LightNoise(ref Work, minLight, maxLight, maxSizeL, this);
  progressBar1.Visible = false;
  Result = new Bitmap(origBmp.Width, origBmp.Height,
                                               PixelFormat.Format24bppRgb);
  progressBar1.Visible = true;
  int i1 = 1 + nbyte * origBmp.Width * origBmp.Height / 100;
  for (int i = 0; i < nbyte * origBmp.Width * origBmp.Height; i++)

  {
    if (i % i1 == 1) progressBar1.PerformStep();
    if (Work.Grid[i] == 252 || Work.Grid[i] == 254) Work.Grid[i] = 255;
  }
  progressBar1.Visible = false;
  GridToBitmap(Result, Work.Grid);
  pictureBox2.Image = Result;
  Graphics g = pictureBox1.CreateGraphics();
  Pen myPen = new Pen(Color.Blue);

  // Drawing the excluded rectangles
  for (int n = 0; n < Number; n += 2)
  {
    g.DrawLine(myPen, v[n + 1].X, v[n + 0].Y, v[n + 1].X, v[n + 1].Y);
    g.DrawLine(myPen, v[n + 0].X, v[n + 0].Y, v[n + 1].X, v[n + 0].Y);
    g.DrawLine(myPen, v[n + 0].X, v[n + 0].Y, v[n + 0].X, v[n + 1].Y);
    g.DrawLine(myPen, v[n + 0].X, v[n + 1].Y, v[n + 1].X, v[n + 1].Y);
  }

  progressBar1.Visible = false;
} //**************************** end Impulse noise ***************

stump 调用方法DeleteBit0，删除图像中所有像素的第 0 位，使得标记已经使用的像素成为可能。然后计算图像亮度的直方图，定义最小亮度MinLight和最大亮度MaxLight。

然后调用根据亮度对像素进行分类的方法Sort。没有必要通过已知的分类算法对像素进行分类，因为这将花费太多时间。相反，我们的方法Sort获得图像亮度的直方图histo[256]，并填充索引的二维数组Index[256][ ]。Index[light]的论点light是像素的明度。子数组Index[light][ ]是一个索引数组，其中nPixel[light]是具有亮度light的像素的数量。值nPixel[light]等于histo[light]。子阵列Index[light][ ]的每个元素是一个坐标为(x, y)的像素的索引，该像素的亮度等于light。这里x是一列的编号，y是一行的编号，width是图像中的列数。方法Sort保存第i个像素的索引ind=x+ width*y，其亮度light和坐标(x, y)在Index[light][i]中。因此，图像的所有像素按照亮度的递增顺序进行排序。如果需要降序，可以反方向读取数组Index。下面是定义类CPnoise及其方法的源代码。我们在类CPnoise的定义之外描述了这些方法并展示了它们的代码，尽管是根据 C#的规定。他们应该在里面。这种方法更方便，因为我们可以在方法代码附近找到方法的描述。

class CPnoise
{ unsafe
  public int[][] Index; // saving indexes of all pixels ordered by lightness
  int[] Comp; // contains indexes of pixels of the connected subset
  int[] nPixel; // number of pixels with certain lightness
  int MaxSize;// admissible size of the connected subset
  Queue Q1;
  unsafe
  public CPnoise(int[] Histo, int Qlength, int Size)  // Constructor
  {
    this.maxSize = Size;
    this.Q1 = new Queue(Qlength);//necessary to find connected subsets
    this.Comp = new int[MaxSize];
    this.nPixel=new int[256];//256 is the number of lightness values
    for (int light = 0; light < 256; light++) nPixel[light] = 0;
    Index = new int[256][];
    for (int light = 0; light < 256; light++)
      Index[light] = new int[Histo[light] + 1];
  } //************************ end of Constructor ******************

正如我们在前面的代码button3_Click中看到的，方法Sort在直方图计算之后被调用。它将图像所有像素的索引插入到索引的二维数组Index[256][ ]中。应该消除脉冲噪声的图像有时包含不是噪声并且不应该被消除的小斑点。例如，图像中人的眼睛有时是小而暗的斑点，其大小类似于暗噪声的斑点的大小。为了防止它们被消除，我们开发了一个简单的程序:用户应该用鼠标在输入图像中画出一些矩形，包括不应该被消除的位置。方法Sort获得这些矩形的参数v[ ]，并跳过将这些矩形内的像素插入到数组Index[256][ ]中。为此，采用简单的方法getCond，该方法使用Form1中定义的Points的数组v[ ]。这里没有给出代码。下面是方法Sort的源代码。

public int Sort(CImage Image, int[] histo, int Number, int picBox1Width,
                                              int picBox1Height, Form1 fm1)
{
int light, i;
double ScaleX = (double)picBox1Width / (double)Image.width;
double ScaleY = (double)picBox1Height / (double)Image.height;
double Scale; // Scale of the presentation of the image in "pictureBox1"
if (ScaleX < ScaleY) Scale = ScaleX;
else Scale = ScaleY;
bool COLOR;
if (Image.nBits == 24) COLOR = true;
else COLOR = false;
double marginX = (double)(picBox1Width - Scale *Image.width)*0.5; // left of image
double marginY=(double)(picBox1Height - Scale*Image.height)*0.5; // above image
bool Condition = false; // Skipping pixel (x, y) if it lies in a global rectangles "fm1.v"
fm1.progressBar1.Value = 0;
fm1.progressBar1.Step = 1;
fm1.progressBar1.Visible = true;
fm1.progressBar1.Maximum = 100;
for (light = 0; light < 256; light++) nPixel[light] = 0;
for (light = 0; light < 256; light++)
  for (int light1 = 0; light1 < histo[light]; light1++) Index[light][light1] = 0;

int y1 =  Image.height / 100;
for (int y = 1; y < Image.height; y++) //=================================
{
  if (y % y1 == y1 - 1) fm1.progressBar1.PerformStep();
  for (int x = 1; x < Image.width; x++) //================================
  {
    Condition = false;
    for (int k = 0; k < Number; k += 2)
      Condition = Condition || getCond(k, x, y, marginX, marginY, Scale, fm1);
    if (Condition)  continue;
    i = x + y * Image.width; // Index of the pixel (x, y)
    if (COLOR)
      light = MaxC(Image.Grid[3 * i+2] & 254, Image.Grid[3 * i + 1] & 254,
                                              Image.Grid[3 * i + 0] & 254);
    else light = Image.Grid[i] & 252;
    if (light < 0) light = 0;
    if (light > 255) light = 255;
    Index[light][nPixel[light]] = i; // record of index "i" of a pixel with lightness "light"
    if (nPixel[light] < histo[light]) nPixel[light]++;
  } //============== end for (int x=1; .. =============================
} //=============== end for (int y=1; .. ==============================
fm1.progressBar1.Visible = false;
return 1;

} //***************** end Sort *****************************************

方法Sort读取图像的所有像素，并用像素的索引填充二维数组Index[256][ ]。这个数组由方法DarkNoise和LightNoise使用。

我们首先考虑暗噪声的情况。方法DarkNoise包含一个带有变量light的for循环，它从最大亮度减 2 开始扫描light的值，直到最小亮度minLight。这个循环包含另一个带有变量i的循环，它从 0 到nPixel[light]扫描二维数组Index[light][i]的第二个索引的值。它从数组Index[light][i]中读取像素的索引，并测试具有该索引的像素是否具有等于light的亮度以及其标签是否等于零。如果两个条件都满足，那么对索引为Index[light][i]的像素调用方法BreadthFirst_D。该方法构建亮度小于或等于light的像素的连通子集，包含索引为Index[light][i]的起始像素。

下面是DarkNoise的源代码:

public int DarkNoise(ref CImage Image, int minLight, int maxLight, int maxSize, Form1 fm1)
{
  bool COLOR = (Image.nBits == 24);
  int ind3 = 0, // index multiplied with 3
    Label2, Lum, rv = 0;
  if (maxSize == 0) return 0;
  fm1.progressBar1.Maximum = 100;
  fm1.progressBar1.Step = 1;
  fm1.progressBar1.Value = 0;
  fm1.progressBar1.Visible = false;
  int bri1 = 2;
  fm1.progressBar1.Visible = true;
  for (int light = maxLight - 2; light >= minLight; light--) //==============
  {
    if ((light % bri1) == 1) fm1.progressBar1.PerformStep();
    for (int i = 0; i < nPixel[light]; i++) //=============================
    {
      ind3 = 3 * Index[light][i];
      if (COLOR)
      {
        Label2 = Image.Grid[2 + ind3] & 1;
        Lum = MaxC(Image.Grid[2 + ind3] & 254, Image.Grid[1 + ind3] & 254,
                                               Image.Grid[0 + ind3] & 254);
      }
      else
      {
        Label2 = Image.Grid[Index[light][i]] & 2;
        Lum = Image.Grid[Index[light][i]] & 252;
      }
      if (Lum == light && Label2 == 0)
      {
        rv = BreadthFirst_D(ref Image, i, light, maxSize);
      }
    } //================= end for (int i.. ============================
  } //================== end for (int light.. ==========================
  fm1.progressBar1.Visible = false;
  return rv;
} //*********************** end DarkNoise ******************************

方法DarkNoise调用方法BreadthFirst_D( )，将值Image、i、光和maxSize作为其参数。maxSize的值是所寻找的亮度小于或等于light的连通像素组中的最大允许像素数。该方法实现了已知的宽度优先搜索算法，该算法被设计来标记树结构的所有顶点。亮度小于或等于light的所有像素的连接集合是一棵树。下面是BreadthFirst_D( )的代码:

private int BreadthFirst_D(ref CImage Image, int i, int light, int maxSize)
{
   int lightNeib,   lightness of the neighbor
       index, Label1, Label2,
       maxNeib, // maxNeib is the maximum number of neighbors of a pixel
       Neib, // the index of a neighbor
       nextIndex, // index of the next pixel in the queue
       numbPix; // number of pixel indexes in "Comp"
  bool small; // equals "true" if the subset is less than "maxSize"
  bool COLOR = (Image.nBits == 24);
  index = Index[light][i];
  int[] MinBound = new int[3]; // color of pixel with min. lightness near the subset
  for (int c = 0; c < 3; c++) MinBound[c] = 300; // a value out of [0, 255]
  for (int p = 0; p < MaxSize; p++) Comp[p] = -1; // MaxSize is element of the class
  numbPix = 0;
  maxNeib = 8; // maximum number of neighbors
  small = true;
  Comp[numbPix] = index;
  numbPix++;
  if (COLOR)
    Image.Grid[1 + 3 * index] |= 1; // Labeling as in Comp (Label1)
  else
    Image.Grid[index] |= 1; // Labeling as in Comp
  Q1.input = Q1.output = 0;
  Q1.Put(index); // putting index into the queue
  while (Q1.Empty() == 0) //=  loop running while queue not empty ============
  {

    nextIndex = Q1.Get();
    for (int n = 0; n <= maxNeib; n++) // == all neighbors of "nextIndex"
    {
      Neib = Neighb(Image, nextIndex, n); // the index of the nth neighbor of nextIndex
      if (Neib < 0) continue; // Neib < 0 means outside the image
      if (COLOR)
      {
        Label1 = Image.Grid[1 + 3 * Neib] & 1;
        Label2 = Image.Grid[2 + 3 * Neib] & 1;
        lightNeib = MaxC(Image.Grid[2 + 3 * Comp[m]],
                         Image.Grid[1 + 3 * Comp[m]], Image.Grid[0 + 3 * Comp[m]]) & 254;
      }
      else
      {
        Label1 = Image.Grid[Neib] & 1;
        Label2 = Image.Grid[Neib] & 2;
        lightNeib = Image.Grid[Neib] & 252; // MaskGV;
      }
      if (lightNeib == light && Label2 > 0) small = false;
      if (lightNeib <= light) //-------------------------------------------
      {
        if (Label1 > 0) continue;
        Comp[numbPix] = Neib; // putting the element with index Neib into Comp
        numbPix++;
        if (COLOR)
          Image.Grid[1 + 3 * Neib] |= 1; // Labeling with "1" as in Comp
        else
          Image.Grid[Neib] |= 1; // Labeling with "1" as in Comp
        if (numbPix > maxSize)  // Very important condition
        {
          small = false;
          break;
        }
        Q1.Put(Neib);
      }

      else // lightNeib<light
      {
        if (Neib != index) //----------------------------------------------
        {
          if (COLOR)
          {
            if (lightNeib < (MinBound[0] + MinBound[1] + MinBound[2]) / 3)
              for (int c = 0; c < 3; c++) MinBound[c] = Image.Grid[c + 3 * Neib];
          }
          else
            if (lightNeib < MinBound[0]) MinBound[0] = lightNeib;
        } //---------------- end if (Neib!=index) -------------------------
      } //----------------  end if (lightNeib<=light) and else ------------
    } // ===================== end for (n=0; .. ========================
    if (small == false) break;
  } // ===================== end while =============================
  int lightComp; // lightness of a pixel whose index is contained in "Comp"
  for (int m = 0; m < numbPix; m++) //==================================
  {

    if (small != false && MinBound[0] < 300) //"300" means MinBound not calculated ---
    {
      if (COLOR)
        for (int c = 0; c < 3; c++) Image.Grid[c + 3 * Comp[m]] = (byte)MinBound[c];
      else
        Image.Grid[Comp[m]] = (byte)MinBound[0];
    }
    else
    {
      if (COLOR)
        lightComp = MaxC(Image.Grid[3 * Comp[m]],
                        Image.Grid[1 + 3 * Comp[m]], Image.Grid[2 + 3 * Comp[m]]) & 254;
      else
        lightComp = Image.Grid[Comp[m]] & 252; // MaskGV;
      if (lightComp == light) //-------------------------------------------
      {
        if (COLOR) Image.Grid[2 + 3 * Comp[m]] |= 1;
        else Image.Grid[Comp[m]] |= 2;
      }
      else // lightComp!=light
      {
        if (COLOR)
        {
          Image.Grid[1 + 3 * Comp[m]] &= (byte)254; // deleting label 1
          Image.Grid[2 + 3 * Comp[m]] &= (byte)254; // deleting label 2
        }
        else
          Image.Grid[Comp[m]] &= 252; // deleting the labels
      } //-------------- end if (lightComp == light) and else--------------
    } //-----------------end if (small != false) and else -----------------
  } //==================== end for (int m=0 .. ========================
  return numbPix;
} //*********************** end BreadthFirst_D ***************************

方法BreadthFirst_D获得阵列Index[light][i]中像素 P 的亮度light和数量i作为参数。该方法的任务是找到亮度小于或等于light并包含像素 p 的像素的连通集合 S ，如果该集合包含小于或等于参数maxSize的像素数量numbPix，则 S 的像素获得较亮的颜色，并作为暗像素变得不可见。否则，S 的所有像素获得一个标签，该标签指示这些像素具有亮度light并且属于多于maxSize个像素的连通子集。

该方法从Index[light][i]中取出index，将其放入队列Q中，并且只要Q不为空，就开始一个while循环运行。在循环中，从队列Q中取出值nextIndex，并且测试其索引为Neib的八个邻居中的每一个的亮度是否小于或等于light。除此之外Neib的标签也测试过。

我们使用Image.Grid中像素的最低有效位来标记像素。在灰度图像的情况下，我们使用两个最低有效位:位 0 是LabelQ1，位 1 是LabelBig2。在彩色图像的情况下，我们使用绿色通道的一位作为LabelQ1，红色通道的一位作为LabelBig2。

如果索引已经放入队列，则设置LabelQ1。如果像素属于大于maxSize的像素的连通大集合 S 并且像素的亮度等于light，则LabelBig2被设置。

如果索引为Neib的像素的LabelQ1没有设置，并且其亮度小于或等于light，那么索引Neib被保存在包含所有亮度小于或等于light的像素的索引的数组Comp中。然后索引Neib被放入队列，像素获得标签LabelQ1。然而，如果Neib的亮度大于light，该亮度将用于计算颜色MinBound，该颜色用于填充要变亮的一小组 S 的所有像素。

这两个标签对于使方法BreadthFirst_D更快都很重要:一个设置了LabelQ1的像素不能被再次输入到队列中，一个设置了LabelBig2的像素表示包含该像素的集合 S 不小。方法DarkNoise不为标有LabelBig2的像素调用方法BrightFirst_D。

在BrightFirst_D的while循环中，从队列中取出索引nextIndex作为它的下一个元素。nextIndex的所有八个邻居的索引Neib通过方法Neighb(Image, nextIndex, n) ( n从 0 到 7)依次传递。如果邻居Neib位于图像之外(nextIndex位于图像的边界)，则方法Neighb返回负值，并且Neib将被忽略。否则Neib的亮度lightNeb取自Image。如果lightNeb大于light，则Neib不属于组成比light更暗或与light同样暗的像素的连通集合(分量)的像素集合 S 。值lightNeb然后用于计算颜色MinBound，该颜色是与 S 相邻的像素中亮度最小的颜色。

如果lightNeb等于light并且像素Neib设置了LabelBig2，则逻辑变量small设置为false。否则，如果lightNeb大于等于light，有两种情况:在Neib设置了LabelQ1的情况下，忽略像素Neib；否则，将已经包含在搜索集合 S 中的像素数量numbPix与maxSize进行比较。如果大于maxSize，那么变量small被设置为false并且while循环被中断。然而，如果numbPix小于或等于maxSize，那么Neib被包含到集合 S ( Neib被保存在数组Comp中)，像素Neib得到LabelQ1集合，Neib被放入队列。

在while循环结束后，正在检查small的值。如果为真，那么索引保存在数组Comp中的 S 的所有像素得到颜色MinBound。否则亮度等于light的Comp的像素得到标签LabelBig2。 S 的所有其他像素没有标签。该方法返回值numbPix。

方法BreadthFirst_L与BreadthFirst_D类似。区别主要在于互换了一些>和<操作符(而不是指令if (numbPix > maxSise) {...})。方法BreadthFirst_L被方法LightNoise调用，类似于DarkNoise。方法LightNoise从最低亮度开始读取数组Index并调用方法BreadthFirst_L。

所有方法都适用于灰度和彩色图像。这是由于采用了方法MaxC(Red, Green, Blue)(参见第三章中的解释)，该方法根据以下简单公式计算颜色为(Red, Green, Blue)的像素的亮度:

明度= max(0.713* R，G ，0.527* B )。

读者可以在图 2-9 中看到一个将该项目应用于德国著名女画家 Louise Seidler 的一幅绘画作品的老照片的例子。

img/474294_1_En_2_Fig9_HTML.jpg

图 2-9

将项目PulseNoiseWF应用于 Louise Seidler 的一幅旧照片的例子(左)。小于 380 像素的暗点被移除(右图)。

如第四章和第五章所述，当我们使用阴影校正后的脉冲噪声抑制时，代表图纸旧照片的灰度图像的结果要好得多。图 2-9 中的图像也将在应用阴影校正方法后看起来更好。

图 2-10 是另一个有轻噪点的图像的例子。这是路易丝·塞德勒的油画《耶稣和孩子们》的照片片段。

img/474294_1_En_2_Fig10_HTML.jpg

图 2-10

将项目PulseNoiseWF应用于路易丝·塞德勒(左)的一张照片片段的示例。小于 30 像素的光斑被去除(右图)。**

三、对比度增强

我们将数字图像的对比度定义如下:

CIP=(李 max - 李min)/域

其中李 max 为最大值李 min 为图像中的最小明度域为明度的域或差值的最大可能值李max-李 min 。

刚刚定义的对比度度量显然取决于单个像素的亮度。为了获得更鲁棒的测量，有必要知道亮度值的频率:出现在少量像素中的值并不重要。频率包含在图像的直方图中。可以找到这样的亮度值MinInp和MaxInp，亮度低于MinInp和亮度高于MaxInp的像素数量小于预定值，我们称之为丢弃区域(因为像素数量是面积的度量)。则鲁棒对比度度量为

cr=(MaxInp–MinInp)/Domain。

自动线性对比度增强

大多数对比度增强程序，如 IrfanView，通过中间部分导数大于 1 的函数，使用图像中一组亮度值到另一组亮度值的映射(见图 3-1 )。

img/474294_1_En_3_Fig1_HTML.png

图 3-1

在其他程序中用于增加对比度的函数的图形示例

这里提出了另一种基于研究图像直方图的方法。我们算法的思想在于将图像中最低和最高亮度之间的间隔映射到整个间隔上；比如到[0, 255]上。重要的是要注意，具有极低或极高值的少量像素可以从本质上影响映射的结果。

因此，有必要从映射中排除具有极值的像素的预定部分(例如，1%)。为此，需要找到值MinInp和MaxInp，使得亮度小于MinInp和亮度大于MaxInp的像素数量小于图像中像素总数的 1%。

为了增加图像的对比度，我们必须增加最大和最小亮度之间的差异。这可以通过减少图像中亮度的小值并增加较大的值来实现。为了获得 1 的最大可能对比度，有必要以这样的方式变换亮度值，即用 0 代替MinInp和用 255 代替MaxInp，并且中间值单调变化；即，它们的顺序保持不变。

可以通过自动程序增加对比度。用户必须仅以图像大小的百分比来定义丢弃区域的值。建议使用查找表(LUT)来加快输出亮度的计算。LUT 是一个数组，包含尽可能多的输入值；例如 256 或从 0 到 255。LUT 的每个元素都包含预先计算的输出值。因此，只计算 256 次输出值就足够了，而不是宽度×高度的倍数，尽管宽度×高度的数值可能高达数千甚至数百万。从数组中获取一个值比计算这个值花费的时间要少得多。

也容易降低对比度。当计算 LUT 时，设置以下值就足够了:MinInp=0、MaxInp=255、MinOut>0和MaxOut<255。大幅降低对比度，然后自动提高，这是一个有趣的实验！

增加对比度的最简单方法(见图 3-2 )是通过以下分段线性变换:

img/474294_1_En_3_Fig2_HTML.png

图 3-2

最简单的分段线性对比度增强

if (Linp <= MinInp) Lout = 0;
else if (Linp >= MaxInp) Lout = 255;
        else Lout = 255*(Linp - MinInp)/(MaxInp – MinInp);

这种非常简单的方法产生了非常好的结果，但对于包含大面积黑暗或大面积明亮区域的图像来说却不是这样。在这种情况下，当采用直方图均衡化方法时，有可能获得更好的结果，这将在下一节中解释。

直方图均衡

一些图像包含较大的低对比度暗区(参见图 3-3a )。先前描述的最简单的分段线性方法均匀地拉伸亮度值，而对于这样的图像，暗区域中的低值必须比其他区域中的值更强烈地分开。这种变换不能用最简单的分段线性方法来完成。

img/474294_1_En_3_Fig3_HTML.jpg

图 3-3

对比度变化的例子:(a)带直方图的原图；(b)直方图对比度降低；(c)线性对比度增强；直方图均衡化

增加这种图像的对比度的一个众所周知的想法是直方图均衡化 : 有必要以这样一种方式转换亮度值，使得所有值获得几乎相同的频率。为了实现这一点，具有相对较高频率的灰度值 inp 1 的像素(图 3-4a )必须被许多具有较低频率的不同灰度值的像素替换。这些灰度值 out 1 、 out 2 等等必须组成一个相邻值的序列(图 3-4b )。没有已知的方法来决定这些像素中的哪些必须获得某个值。幸运的是，这不是必须的:以平均频率保持不变的方式，用具有更大差值的其他值替换原始灰度值就足够了(图 3-4c )。这可以通过直方图均衡化来实现。

img/474294_1_En_3_Fig4_HTML.png

图 3-4

变换直方图以使频率或平均频率恒定:(a)具有大频率的直方图部分；(b)具有高频率的灰度值被许多具有较低频率的值代替；(c)具有使平均频率恒定的较大差值的值

为了实现直方图均衡，有必要计算累积直方图。这是一个有 256 个元素的整数数组。对应于亮度值 L 的元素包含亮度值小于或等于 L 的像素数量。为了计算 L 的累积直方图的值，将小于或等于 L 的所有值的常用直方图的值相加就足够了:

Cum\left[L\right]=\sum \limits_{i\le L} Hist\left[i\right];

直方图均衡化的 LUT 可以通过累积直方图来计算。大多数关于图像处理的教科书建议对 LUT 的元素使用以下公式:

eq[l= 255cum[l]/max cum；(3.1)

其中 MaxCum 是累积直方图的最大值，显然等于图像的尺寸宽度*高度。

然而，当使用该公式时，如果原始图像包含大量灰度值等于 gMin，的像素，这是直方图大于 0 的最小灰度值，则会出现困难。在这种情况下，通过根据等式 3.1 计算的 LUT 变换的图像可能非常差，因为许多像素具有与比 0 大得多的和 [ gMin ]成比例的大灰度值。这种转换的一个例子如图 3-5b 所示。

img/474294_1_En_3_Fig5_HTML.jpg

图 3-5

直方图均衡的例子:(a)原始图像；(b)根据等式 3.1 用 LUT 变换；(c)根据等式 3.2 用 LUT 转换

为了解决这个问题，设置 LUT 等式的第 L 个元素的值与差值Cum[L–MinHist成比例就足够了，而不是与 Cum [ L 成比例，其中MinHist=Cum[gMin]=LUT 的元素必须根据以下公式计算:

如果 L 小于或等于 gMin ，则等式[L]= 0；其他

eq[l]= 255(cum**l-minorist)/max cum

其中 MaxCum= width*height 为 Cum [ L 的最大值。

图 [3-5c 显示了根据等式 3.2 计算的 LUT 的变换示例。

直方图均衡化的结果并不总是令人满意的。因此，合理的做法是使用分段线性方法与均衡方法的组合，称为混合方法(见图 3-6 ):

作战[【l】]=【w】**【lut】[]+(3.3)

其中 LUT L 为分段线性 LUT，权重 W 1 和 W 2 可以选择为和等于 1 的非负刹车。图 [3-6c 显示了使用 CombiLUT 与 W 1 = 0.2 和 W 2 = 0.8 的结果。

img/474294_1_En_3_Fig6_HTML.jpg

图 3-6

组合对比度增强:(a)原始图像；分段线性的；(c)直方图均衡；混合方法(20%)

测量彩色图像的亮度

为了定义彩色图像的对比度，需要计算每个像素的亮度或明度。在文献中，彩色像素的亮度有许多不同的定义。维基百科在论文“HSL 和 HSV”中提出了“明度”的不同定义；例如:

I = (R+G+B)/3 或

V = max(R，G，B)或

L = (max(R，G，B)+min(R，G，B))/2 或

Y = 0.30R + 0.59G + 0.11B

其中 R 、 *G、*和 B 分别是红色、绿色和蓝色通道的强度。作者开发了一个项目WFluminance来实验测试不同的亮度定义。

该项目在表单的左侧显示了八个不同颜色的矩形(图 3-7 ，并提供了通过相应的工具numericUpDown改变每个矩形的亮度而不改变其色调的可能性。每个彩色矩形的亮度已经改变，直到所有矩形看起来都具有相同的视觉上可接受的亮度，等于左边灰色矩形的亮度。众所周知，灰色具有彼此相等的所有三个颜色通道的强度:R = G = B。因此，可以假设灰色的亮度等于该恒定强度。

img/474294_1_En_3_Fig7_HTML.jpg

图 3-7

视觉亮度恒定的矩形

该项目根据不同的定义计算每个矩形的亮度值。结果如表 3-1 所示。从表 3-1 中可以看出，对于图 3-7 的所有八个矩形，MC 的值几乎是恒定的，看起来它们具有恒定的亮度。如图 3-7 右侧所示，其他颜色的 MC 值也测试成功，代表视觉亮度相等的矩形对。

表 3-1

根据不同定义的亮度值

|

矩形

|

颜色

|

稀有

|

G

|

B

|

M

|

L

|

Y

|

主持人

| | --- | --- | --- | --- | --- | --- | --- | --- | --- | | one | 蓝色 | Zero | Zero | Two hundred and forty-two | Two hundred and forty-two | One hundred and twenty-one | Seventeen | One hundred and twenty-seven | | Two | 品红 | One hundred and sixty | Zero | One hundred and sixty | One hundred and sixty | Eighty | Forty-five | One hundred and fourteen | | three | 红色 | One hundred and seventy-four | Zero | Zero | One hundred and seventy-four | Eighty-seven | Thirty-six | One hundred and twenty-four | | four | 黄色 | One hundred and twenty-two | One hundred and twenty-two | Zero | One hundred and twenty-two | Sixty-one | One hundred and thirteen | One hundred and twenty-two | | five | 格林（姓氏）；绿色的 | Zero | One hundred and twenty-four | Zero | One hundred and twenty-four | Sixty-two | Eighty-eight | One hundred and twenty-four | | six | 蓝绿色 | Zero | One hundred and twenty-two | One hundred and twenty-two | One hundred and twenty-two | Sixty-one | Ninety-six | One hundred and twenty-two | | seven | 蓝色 | Zero | Zero | Two hundred and forty-two | Two hundred and forty-two | One hundred and twenty-one | Seventeen | One hundred and twenty-seven | | eight | 灰色的 | One hundred and twenty-one | One hundred and twenty-one | One hundred and twenty-one | One hundred and twenty-one | One hundred and twenty-one | One hundred and twenty-one | One hundred and twenty-one |

注:M = max(R，G，B)；L = (M + min((R，G，B))；y = 0.2126 * R+0.7152 * G+0.0722 * R；MC = max(0.713R，G，0.527B)。

右半部分颜色矩形左侧的文本显示矩形中颜色的 MC 值。项目的用户可以通过改变相应numericUpDown工具的值来改变相邻灰色矩形的灰度值，直到灰色矩形的视觉亮度等于彩色矩形的视觉亮度。如果numericUpDown的值近似等于 MC 值，则这证明 MC 方法给出了该颜色亮度的良好估计。如您所见，所有矩形的值大致相等。

有一个问题是彩色像素亮度的精确定义是否重要。此定义用于将彩色图像转换为灰度图像以及阴影校正。我们已经尝试使用不同的方法来计算彩色像素的亮度，同时将不同的彩色图像转换成灰度图像。在大多数情况下，结果并不取决于方法的选择。然而，在图 3-8 所示的例子中，当使用方法 Y(表 3-1 )时，彩色图像的相当大一部分消失了。

我们还测试了阴影校正的结果如何依赖于亮度测量方法的选择。与将彩色图像转换成灰度图像的情况类似，阴影校正的结果大多不取决于该选择。但是，有些图像的阴影校正强烈依赖于该选择。例如，参见图 3-9 中的图像。

img/474294_1_En_3_Fig9_HTML.png

图 3-9

阴影校正强烈依赖于亮度测量方法的图像示例；方法 Y 和方法 I

img/474294_1_En_3_Fig8_HTML.png

图 3-8

方法 Y 下零件消失的图像示例

考虑到刚才描述的结果，我们决定在任何地方都使用 MC 方法(作为方法MaxC)。

彩色图像的对比度

为了执行彩色图像的对比度增强，有必要计算亮度的直方图，然后像在灰度图像的情况下那样处理它。用于亮度变换的 LUT 的计算可以以与之前相同的方式进行，或者利用具有丢弃区域的分段线性规则，或者利用直方图均衡化。然而，为了改变颜色通道的强度，应该使用新的过程:需要为每个像素计算转换后的亮度与原始亮度的关系，然后将三个强度 R 、 *G、*和 B 中的每一个乘以该关系。例如，像素 P 的变换红色通道为:

rt=【ro】**【lit/lio】

其中 RT 和 LiT 是变换后的红色强度和亮度，而 RO 和 LiO 是像素 P 的红色强度和亮度的原始值。以这种方式，所有三种强度可以成比例地改变，并且色调保持不变，由此对比度增加。

手动控制对比度增强

一些照片具有相当不均匀的亮度，这是由被摄物体的不均匀照明引起的。照片的这一缺点可以通过直方图均衡化方法得到很大程度的纠正，直方图均衡化方法通常是一种很好的对比度增强方法。但是，有时需要在查看结果的同时更正 LUT 的内容。

一些建议的图像处理程序(例如，Photoshop)包含允许将 LUT 的图形制作成贝塞尔曲线的工具，贝塞尔曲线是具有连续导数的平滑曲线。然而，这是一个坏主意，因为由 LUT 实现的函数一定不是平滑的。它必须是连续的，更重要的是，单调递增的。贝塞尔曲线并不总是满足最后一个条件。如果 LUT 的函数不是单调递增的，那么得到的图像几乎不会失真。

将 LUT 函数定义为如图 3-10 所示的分段线性单调递增的直线要简单得多。

img/474294_1_En_3_Fig10_HTML.png

图 3-10

LUT 的分段线性函数

断点 V1 和 V2 可以通过鼠标点击来定义。转换后的图像必须在断点定义后立即显示。

我们开发了 Windows 窗体项目WFpiecewiseLinear来实现这个过程。运行项目的形式如图 3-11 所示。

img/474294_1_En_3_Fig11_HTML.jpg

图 3-11

运行项目的形式示例WFpiecewiseLinear

该表单包含两个按钮。当用户单击打开图像时，他或她可以选择一个图像。图像将显示在左边的图片框中。同时，图像的副本出现在右边的图片框中。图像的直方图出现在表单的左下角，蓝色线条显示 LUT 的函数图。用户现在可以点击线附近的两个点。线以这样一种方式改变，即它的断点位于被单击的点上。LUT 相应变化，用户可以观察到右边图片框中图像的选择。一旦用户对图像的外观感到满意，他或她就可以单击保存结果。然后，用户可以选择名称并保存结果图像。保存图像的类型可以是Bmp或Jpg。

如前一部分所解释的，为了计算颜色通道的新强度，有必要通过 LUT 计算像素的新亮度，然后将每个旧的颜色强度乘以新亮度，并除以旧亮度。必须对图像的每个像素执行该过程。由于乘法和除法，这花费了太多的时间，并且变换后的图像出现延迟。

为了使变换图像的计算更快，我们开发了一种使用bigLUT的方法。这是一个颜色通道的新强度的 2 16 = 65，536 个值的 LUT。作为该表的一个参数，我们使用一个整数，其较低值字节是旧的颜色强度，第二个字节是亮度的旧值(不使用字节 3 和 4 ),取自先前准备的原始图像的灰度值版本。新色彩强度的 65，536 个值中的每一个都等于旧色彩强度与 LUT 值的乘积除以旧亮度。这些值在定义bigLUT时只计算一次。在项目运行期间，只需为每个像素计算自变量。新的强度仅仅是从bigLUT中读出的。下面是计算bigLUT的方法代码:

    int[] bigLUT = new int[65536];

    private void makeBigLUT(int[] LUT)
    { for (int colOld = 0; colOld < 256; colOld++)  // old color intensity
       { for (int lightOld = 1; lightOld < 256; lightOld++) //old lightness
          { int colNew = colOld * LUT[lightOld] / lightOld;
             int arg = (lightOld << 8) | colOld;
             bigLUT[arg] = colNew;
          }
       }
    }

接下来是项目最重要部分的源代码。当单击“打开图像”按钮时，第一部分开始。这个部分打开图像，生成位图origBmp，定义类CImage的三个对象，并用origBmp的内容填充对象origIm。软件WindowsForms建议通过GetPixels的方法访问位图的像素，这个方法比较慢。我们开发了一个更快的方法BitmapToGrid，它使用了类Bitmap的方法LockBits。然而，BitmapToGrid不应该用于每像素 8 位的索引位图，因为这些位图包含一个调色板，像素值总是由调色板值定义。因此，对于 8 位索引位图，方法BitmapToGrid变得太慢。对于这样的位图，我们使用带有GetPixel的标准方法。

此外，我们项目的第一部分使用我们的方法MaxC计算彩色像素的亮度，用方法ColorToGrayMC将原始图像origIm转换成灰度图像origGrayIm。它还计算origGrayIm的直方图，绘制pictureBox3中的直方图，计算包含分段线性曲线(没有对比度增强)的值的标准 LUT，并绘制代表标准 LUT 的函数的线。它还将图像contrastIm定义为结果图像，通过类似于前面提到的方法BitmapToGrid的快速方法GridToBitmap在位图ContrastBmp中制作其副本，并在pictureBox2中显示ContrastBmp。

private void button1_Click(object sender, EventArgs e) // Open image
{
  OpenFileDialog openFileDialog1 = new OpenFileDialog();
  if (openFileDialog1.ShowDialog() == DialogResult.OK)
  {
    try
    {
      origBmp = new Bitmap(openFileDialog1.FileName);
      pictureBox1.Image = origBmp;
    }
    catch (Exception ex)
    {
      MessageBox.Show("Error: Could not read file from disk. Original error: " +
        ex.Message);
    }
  }
  width = origBmp.Width;
  height = origBmp.Height;

  origIm = new CImage(width, height, 24);
  origGrayIm = new CImage(width, height, 8);  // grayscale version of origIm
  contrastIm = new CImage(width, height, 24);

  if (origBmp.PixelFormat == PixelFormat.Format8bppIndexed)
  {
    progressBar1.Visible = true;

    Color color;
    int y2 = height / 100 + 1, nbyte = 3;
    for (int y = 0; y < height; y++) //=====================
    {
      if (y % y2 == 1) progressBar1.PerformStep();
      for (int x = 0; x < width; x++) //=================
      {
        int i = x + width * y;
        color = origBmp.GetPixel(i % width, i / width);
        for (int c = 0; c < nbyte; c++)
        {
          if (c == 0) origIm.Grid[nbyte * i] = color.B;
          if (c == 1) origIm.Grid[nbyte * i + 1] = color.G;
          if (c == 2) origIm.Grid[nbyte * i + 2] = color.R;
        }
      } //========= end for ( int x... =================
    } //========== end for ( int y... ==================
    progressBar1.Visible = false;
  }
  else BitmapToGrid(origBmp, origIm.Grid);

  // Calculating the histogram:
  origGrayIm.ColorToGrayMC(origIm, this);
  for (int gv = 0; gv < 256; gv++) histo[gv] = 0;
  for (int i = 0; i < width * height; i++) histo[origGrayIm.Grid[i]]++;
  MaxHist = 0;
  for (int gv = 0; gv < 256; gv++) if (histo[gv] > MaxHist)

MaxHist = histo[gv];
  MinGV = 255;
  MaxGV = 0;
  for (MinGV = 0; MinGV < 256; MinGV++)
    if (histo[MinGV] > 0) break;
  for (MaxGV = 255; MaxGV >= 0; MaxGV--)
    if (histo[MaxGV] > 0) break;

  // Drawing the histogram:
  Graphics g = pictureBox3.CreateGraphics();
  SolidBrush myBrush = new SolidBrush(Color.LightGray);
  Rectangle rect = new Rectangle(0, 0, 256, 256);
  g.FillRectangle(myBrush, rect);

  Pen myPen = new Pen(Color.Red);
  Pen greenPen = new Pen(Color.Green);
  for (int gv = 0; gv < 256; gv++)
  {
    int hh = histo[gv]*255/MaxHist;
    if (histo[gv] > 0 && hh < 1) hh = 1;
    g.DrawLine(myPen, gv, 255, gv, 255 - hh);
    if (gv == MinGV || gv == MaxGV)
      g.DrawLine(greenPen, gv, 255, gv, 255 - hh);
  }

  // Calculating the standard LUT:
  int[] LUT = new int[256];
  int X = (MinGV + MaxGV) / 2;
  int Y = 128;
  for (int gv = 0; gv < 256; gv++)
  {
    if (gv <= MinGV) LUT[gv] = 0;
    if (gv > MinGV && gv <= X) LUT[gv] = (gv - MinGV) * Y / (X - MinGV);
    if (gv > X && gv <= MaxGV) LUT[gv] = Y + (gv - X) * (255 - Y) / (MaxGV - X);
    if (gv >= MaxGV) LUT[gv] = 255;
  }

  // Drawing the standard curve:

  int yy = 255;
  Pen myPen1 = new Pen(Color.Blue);
  g.DrawLine(myPen1, 0, yy, MinGV, yy);
  g.DrawLine(myPen1, MinGV, yy, X, yy - Y);
  g.DrawLine(myPen1, X, yy - Y, MaxGV, 0);
  g.DrawLine(myPen1, MaxGV, 0, yy, 0);

  // nbyte = 3;  origIm and contrastIm are both 24-bit images
  for (int i = 0; i < nbyte * width * height; i++) contrastIm.Grid[i] = (byte)LUT[(int)origIm.Grid[i]];

  ContrastBmp = new Bitmap(width, height, PixelFormat.Format24bppRgb);
  progressBar1.Visible = true;
  GridToBitmap(ContrastBmp, contrastIm.Grid);
  pictureBox2.Image = ContrastBmp;
  progressBar1.Visible = false;
} //***************** end Open image ***********************************

在第二章中描述了将位图内容快速复制到CImage类图像的方法BitmapToGrid。方法GridToBitmap类似，所以我不在这里介绍。

项目中用于控制LUT形状的部分允许用户在包含直方图图像的pictureBox3中点击两次，从而定义 LUT 分段线性曲线的尼克斯点。该部分使用了Form1中定义的整数变量cntClick、X1、Y1、X2和Y2。变量cntClick对点击进行计数，但是它的值3立即被1替换，因此它只能取值1和2。

当cntClick等于1时，则第一个拐点V1的坐标X1和Y1(图 3-11 )由点击点定义。从直方图的起点到计算点(X1, Y1)的分段线性曲线的第一个分支。

当cntClick等于2时，则定义第二个尼克斯点V2的坐标X2和Y2(图 3-11 )。用户应设置该点，使X2大于X1且Y2大于Y1。如果碰巧X2小于X1，那么X2被设置为X1 + 1；类似地Y2。因此分段线性曲线总是保持单调。计算从点(X1, Y1)到点(X2, Y2)的分段线性曲线的第二个分支和从点(X2, Y2)到直方图终点的第三个分支。然后如前所述计算 LUT bigLUT，并在pictureBox2中计算和显示对比图像。

下面是计算 LUT 并同时显示结果图像的部分代码。

private void  pictureBox3_MouseDown(object sender, MouseEventArgs e)
// making new LUT and the resulting image
{
  int nbyte = 3;
  Graphics g = g2Bmp;  //pictureBox3.CreateGraphics();
  SolidBrush myBrush = new SolidBrush(Color.LightGray);
  Rectangle rect = new Rectangle(0, 0, 256, 256);

  cntClick++;
  if (cntClick == 3) cntClick = 1;
  int oldX = -1, oldY = -1;

  if (cntClick == 1)
  {
    X1 = e.X;
    if (X1 < MinGV) X1 = MinGV;
    if (X1 > MaxGV) X1 = MaxGV;
    Y1 = 255 - e.Y; // (X, Y) is the clicked point in the graph of the LUT
    if (X1 != oldX || Y1 != oldY) //---------------------------------------
    {
      // Calculating the LUT for X1 and Y1:
      for (int gv = 0; gv <= X1; gv++)
      {
        if (gv <= MinGV) LUT[gv] = 0;
        if (gv > MinGV && gv <= X1) LUT[gv] = (gv - MinGV) * Y1 / (X1 - MinGV);
        if (LUT[gv] > 255) LUT[gv] = 255;

      }
    }
    oldX = X1;
    oldY = Y1;

  }

  if (cntClick == 2)
  {
    X2 = e.X;
    if (X2 < MinGV) X2 = MinGV;
    if (X2 > MaxGV) X2 = MaxGV;
    if (X2 < X1) X2 = X1 + 1;
    Y2 = 255 - e.Y; // (X2, Y2) is the second clicked point in the graph of the LUT
    if (Y2 < Y1) Y2 = Y1 + 1;

    if (X2 != oldX || Y2 != oldY) //---------------------------------------
    {
      // Calculating the LUT for X2 and Y2:
      for (int gv = X1 + 1; gv < 256; gv++)
      {
        //if (gv <= MinGV) LUT[gv] = 0;
        if (gv > X1 && gv <= X2) LUT[gv] = Y1 + (gv - X1) * (Y2 - Y1) / (X2 - X1);
        if (gv > X2 && gv <= MaxGV)
                                    LUT[gv] = Y2 + (gv - X2) * (255 - Y2) / (MaxGV - X2);
        if (LUT[gv] > 255) LUT[gv] = 255;

        if (gv >= MaxGV) LUT[gv] = 255;
      }
    }
    oldX = X2;
    oldY = Y2;
  }

  Pen myPen1 = new Pen(Color.Blue);
  makeBigLUT(LUT);

  // Drawing the curve for LUT:
  g.FillRectangle(myBrush, rect);

  Pen myPen = new Pen(Color.Red);
  g.DrawLine(myPen, MinGV, 255, MinGV, 250);
  for (int gv = 0; gv < 256; gv++)
  {
    int hh = histo[gv] * 255 / MaxHist;
    if (histo[gv] > 0 && hh < 1) hh = 1;
    g.DrawLine(myPen, gv, 255, gv, 255 - hh);
  }

  int yy = 255;
  g.DrawLine(myPen1, 0, yy-LUT[0], MinGV, yy-LUT[0]);
  g.DrawLine(myPen1, MinGV, yy - LUT[0], X1, yy - LUT[X1]);
  g.DrawLine(myPen1, X1, yy - LUT[X1], X2, yy - LUT[X2]);
  g.DrawLine(myPen1, X2, yy - LUT[X2], MaxGV, 0);
  g.DrawLine(myPen1, MaxGV, 0, 255, 0);

  // Calculating 'contrastIm':
  int[] GV = new int[3];
  int arg, colOld, colNew;
  for (int i = 0; i < nbyte * width * height; i++) contrastIm.Grid[i] = 0;
  progressBar1.Visible = true;
  int y3 = 1 + nbyte * width * height / 100;

  for (int y = 0, yn=0; y < width * height; y += width, yn+=nbyte * width) //======
  {
    if (y % y3 == 1) progressBar1.PerformStep();
    for (int x = 0, xn=0; x < width; x++, xn+=nbyte)

    {
      int lum = origGrayIm.Grid[x + y];
      for (int c = 0; c < nbyte; c++)
      {
        colOld = origIm.Grid[c + xn + yn]; // xn + yn = nbyte*(x + width * y);
        arg=(lum << 8) | colOld;
        colNew = bigLUT[arg];
        if (colNew > 255) colNew = 255;
        contrastIm.Grid[c + xn + yn] = (byte)colNew;
      }
    }
  } //================ end for (int y = 0 ... ======================
  progressBar1.Visible = false;

  // Calculating "ContrastBmp":
  GridToBitmap(ContrastBmp, contrastIm.Grid);
  pictureBox2.Image = ContrastBmp;
  pictureBox3.Image = BmpPictBox3;
  progressBar1.Visible = false;
} //****************** end pictureBox3_MouseDown ***********************

还有一小部分代码用于保存结果图像。它使用对话框SaveFileDialog。如果用户想要将结果图像保存到输入原始图像的文件中，一些附加操作是必要的。输入文件的名称OpenImageFile必须确定并保存在程序部分Open image中。此外，必须定义一个名为tmpFileName的临时文件。生成的图像保存在tmpFileName文件中。那么必须调用下面的方法:

File.Replace(tmpFileName, OpenImageFile,
       OpenImageFile.Insert(OpenImageFile.IndexOf("."), "BuP"));

该方法将tmpFileName的内容复制到输入文件OpenImageFile，删除tmpFileName，并将输入文件的旧内容保存到新文件中，新文件的名称是通过在.字符前的OpenImageFile中插入字符串BuP产生的。以下是项目这一部分的代码。

private void button2_Click(object sender, EventArgs e)  // Save result
{
  SaveFileDialog dialog = new SaveFileDialog();
  dialog.Filter =
       "Image Files(*.BMP;*.JPG;*.GIF)|*.BMP;*.JPG;*.GIF|All files (*.*)|*.*";
  if (dialog.ShowDialog() == DialogResult.OK)
  {
    string tmpFileName;
    if (dialog.FileName == OpenImageFile)
    {
      tmpFileName = OpenImageFile.Insert(OpenImageFile.IndexOf("."), "$$$");
      if (dialog.FileName.Contains("jpg"))
            ContrastBmp.Save(tmpFileName, ImageFormat.Jpeg); // saving tmpFile
      else ContrastBmp.Save(tmpFileName, ImageFormat.Bmp);

      origBmp.Dispose();
      File.Replace(tmpFileName, OpenImageFile,
                     OpenImageFile.Insert(OpenImageFile.IndexOf("."), "BuP"));
      origBmp = new Bitmap(OpenImageFile);
      pictureBox1.Image = origBmp;
    }
    else
    {
      if (dialog.FileName.Contains("jpg"))
             ContrastBmp.Save(dialog.FileName, ImageFormat.Jpeg);
      else ContrastBmp.Save(dialog.FileName, ImageFormat.Bmp);
    }
    MessageBox.Show("The result image saved under " + dialog.FileName);
  }
} //****************** end Save result ******************************

图 3-12 显示了使用所述项目获得的对比度增强的示例。该项目可用于增强彩色和灰度图像的对比度。

img/474294_1_En_3_Fig12_HTML.jpg

图 3-12

项目WFpiecewiseLinear的对比度增强示例

项目WFpiecewiseLinear可以在本书的源代码中找到。

四、阈值阴影校正

有些照片的光线不均匀，因为拍摄的对象没有得到均匀的照明。这种照片的暗区中的细微细节很难辨认。我们建议一种方法来提高这种图像的质量。

一种用于校正由安装在其他星球上的空间探测器上的自动设备拍摄的数字照片的已知方法使用从原始图像的像素亮度中减去局部平均亮度:

GVnew [ x，y ] = GVold [ x，y–表示 [ x，y ] + Const

其中 GVold x，y 是坐标为( x，y )的像素的亮度， Mean [ x，y 是像素( x，y )的邻域中的平均亮度， Const 是不依赖于坐标的亮度值， GVnew 是像素的亮度

这种方法给出的结果大多是充分的。但是，我建议一个稍微不同的方法。考虑从具有反射率 R ( x，y )的表面反射的光的亮度 L ( x，y )，该表面由产生照明 I (x，y)的光源照明:

L ( x，y ) = R ( x，y)；I(x，y))*

其中φ是落光方向和表面法线之间的角度。如果我们有表面照度的估计值 E ( x，y ，并且我们对知道表面的反射率 R ( x，y )感兴趣，并且我们假设 cosφ = 1，那么我们可以通过将观察到的亮度 L ( x， y )除以 E ( x，y )，而不是从 L ( x，y )减去 E ( x，y )。

照度 I ( x，y )主要是随着坐标 x 和 y 的变化而缓慢改变其值的函数，而表面的反射率 R ( x，y )快速变化。因此，可以通过计算每个点( x，y )附近的 L ( x，y )的平均值来获得照度的估计值。在数码摄影的情况下，数码照片像素的亮度与 L ( x，y )成正比，我们感兴趣的是函数 R ( x，y )。所以要得到一个与 R ( x，y )成正比的函数，就要计算出 E ( x，y )作为照片亮度的局部均值，再用 L ( x，y )除以 E ( x，y )。从 L ( x，y )减去 E ( x，y )也可以得到 R ( x，y )的估计值。如果所有的值 L ( x，y )， R ( x，y )， I ( x，y )，以及 E ( x，y )都与相应物理值的对数成正比，这是正确的。但是，使用通常的值 L ( x，y )、 R ( x，y )、 I ( x，y )和 E ( x，y )有时会产生良好的结果。因此，我们在项目中使用减法和除法。用户可以比较结果并选择最好的一个。

我在这里提出的项目WFshadingBin执行阴影校正(见我们在互联网上的软件)和阴影校正图像的阈值。

给定图像的局部平均值可通过第 [2 章中描述的方法Averaging计算。然而，现代数码照片很大。大多数照片包含超过 1000 个⋅1000 = 106像素。用于计算典型图像阴影校正平均值的窗口大小必须至少为 100 ⋅ 100 = 10，000 像素。对于这样的大小，Averaging方法会太慢。最好使用提供快速平均的方法，如第二章所述的FastAverageM。平均窗口的宽度由带有标签Window的工具numericUpDown1指定，作为图像宽度的一部分。这对于通常不知道已处理图像大小的用户来说更方便。

另一方面，对于绘图照片的阴影校正，我们需要一个局部平均窗口，其宽度略大于绘图中线条的宽度。后者的宽度通常在 4 到 9 个像素之间。这可能小于图像宽度的 1 %,图像宽度可能高达 2，000 像素。因此，我们决定将平均窗口的宽度指定为图像宽度的千分之一。

即使对于彩色图像，我们只需要亮度的局部平均值，而不需要颜色通道的平均值。因此，我们必须将彩色图像转换成灰度图像，其像素中包含三色亮度( R，G，B )。方法ColorToGrayMC在为彩色图像的每个像素计算值的同时执行这种转换

亮度 =最大值(R0.713，G，B0.527)

用我们的方法MaxC计算(见第三章)。然后可以使用方法FastAverageM来计算灰度图像的局部平均值。

我接下来描述这个项目WFshadingBin。本项目的形式如图 4-1 所示。

img/474294_1_En_4_Fig1_HTML.jpg

图 4-1

阴影校正后项目WFshadingBin的形式

它包含五个图片框、四个按钮和两个类型为numericUpDown的工具。第一个按钮“打开图像”启动项目的第一部分，与项目的第一部分WFpiecewiseLinear类似，不同之处在于这里定义了七个图像— OrigIm、SigmaIm、SubIm、DivIm、GrayIm、MeanIm和BinIm。图像SigmaIm包含通过西格玛滤波器对原始图像进行滤波的结果。图像用于通过减法和除法以及图像SubIm和DivIm的阈值处理来执行阴影校正。阈值将在下一节中描述。

阴影校正按钮启动项目的相应部分。它包含这里介绍的方法CorrectShading。

public void CorrectShading()
{
  int c, i, x, y;
  bool ROWSECTION = false; // used to show the values of images in one row.
  int[] color = {0, 0, 0};
  int[] color1 = {0, 0, 0};
  int Lightness=(int)numericUpDown2.Value;
  int hWind = (int)(numericUpDown1.Value * width / 2000);
  MeanIm.FastAverageM(GrayIm, hWind, this); // uses numericUpDown1
  progressBar1.Visible = true;
  progressBar1.Value = 0;

  int[] histoSub = new int[256];
  int[] histoDiv = new int[256];
  for (i = 0; i < 256; i++) histoSub[i] = histoDiv[i] = 0;
  byte lum =  0;
  byte lum1 = 0;
  int jump = height / 17; // width and height are properties of Form1
  for (y = 0; y < height; y++) //=======================================
  {
    if (y % jump == jump - 1) progressBar1.PerformStep();
    for (x = 0; x < width; x++)
    {                               // nbyteIm is member of 'Form1'
      for (c = 0; c < nbyteIm; c++) //===================================
      {
        color[c] = Round(SigmaIm.Grid[c + nbyteIm * (x + width * y)] * Lightness / (double)MeanIm.Grid[x + width * y]); // Division
        if (color[c] < 0) color[c] = 0;
        if (color[c] > 255) color[c] = 255;
        DivIm.Grid[c + nbyteIm * (x + width * y)] = (byte)color[c];

        color1[c] = SigmaIm.Grid[c + nbyteIm * (x + width * y)] + Lightness - MeanIm.Grid[x + width * y]; // Subtraction
        if (color1[c] < 0) color1[c] = 0;
        if (color1[c] > 255) color1[c] = 255;
        SubIm.Grid[c + nbyteIm * (x + width * y)] = (byte)color1[c];
      } //=============== end for (c... ===============================
      if (nbyteIm == 1)
      {
        lum = (byte)color[0];
        lum1 = (byte)color1[0];
      }
      else
      {
        lum = SigmaIm.MaxC((byte)color[2], (byte)color[1], (byte)color[0]);
        lum1 = SigmaIm.MaxC((byte)color1[2], (byte)color1[1], (byte)color1[0]);
      }
      histoDiv[lum]++;
      histoSub[lum1]++;
    }
  } //============== end for (y... ===================================

  // Calculating  MinLight and MaxLight for 'Div':
  int MaxLightDiv, MaxLightSub, MinLightDiv, MinLightSub, Sum = 0;
  for (MinLightDiv = 0; MinLightDiv < 256; MinLightDiv++)
  {
    Sum += histoDiv[MinLightDiv];
    if (Sum > width * height / 100) break;
  }
  Sum = 0;
  for (MaxLightDiv = 255; MaxLightDiv >= 0; MaxLightDiv--)
  {
    Sum += histoDiv[MaxLightDiv];
    if (Sum > width * height / 100) break;
  }

  // Calculating  MinLight and MaxLight for 'Sub':
  Sum = 0;
  for (MinLightSub = 0; MinLightSub < 256; MinLightSub++)
  {
    Sum += histoSub[MinLightSub];
    if (Sum > width * height / 100) break;
  }
  Sum = 0;
  for (MaxLightSub = 255; MaxLightSub >= 0; MaxLightSub--)
  {
    Sum += histoSub[MaxLightSub];
    if (Sum > width * height / 100) break;
  }

  // Calculating  MinLight and MaxLight for 'Sub':
  Sum = 0;
  for (MinLightSub = 0; MinLightSub < 256; MinLightSub++)
  {
    Sum += histoSub[MinLightSub];
    if (Sum > width * height / 100) break;
  }
  Sum = 0;
  for (MaxLightSub = 255; MaxLightSub >= 0; MaxLightSub--)
  {
    Sum += histoSub[MaxLightSub];
    if (Sum > width * height / 100) break;
  }

  // Calculating LUT for 'Div':
  byte[] LUT = new byte[256];
  for (i = 0; i < 256; i++)
    if (i <= MinLightDiv) LUT[i] = 0;
    else
      if (i > MinLightDiv && i <= MaxLightDiv)
        LUT[i] = (byte)(255 * (i - MinLightDiv) / (MaxLightDiv - MinLightDiv));
      else LUT[i] = 255;

  // Calculating LUTsub for 'Sub':
  byte[] LUTsub = new byte[256];
  for (i = 0; i < 256; i++)
    if (i <= MinLightSub) LUTsub[i] = 0;
    else
      if (i > MinLightSub && i <= MaxLightSub)
        LUTsub[i] = (byte)(255 * (i - MinLightSub) / (MaxLightSub - MinLightSub));
      else LUTsub[i] = 255;

  // Calculating contrasted "Div" and "Sub":
  for (i = 0; i < 256; i++) histoDiv[i] = histoSub[i] = 0;
  jump = width * height / 17;
  for (i = 0; i < width * height; i++) //==================================
  {
    if (i % jump == jump - 1) progressBar1.PerformStep();

    for (c = 0; c < nbyteIm; c++)
    {
      DivIm.Grid[c + nbyteIm * i] = LUT[DivIm.Grid[c + nbyteIm * i]];
      SubIm.Grid[c + nbyteIm * i] = LUTsub[SubIm.Grid[c + nbyteIm * i]];
    }

    if (nbyteIm == 1)
    {
      lum = DivIm.Grid[0 + nbyteIm * i];
      lum1 = SubIm.Grid[0 + nbyteIm * i];
    }
    else
    {
      lum = SigmaIm.MaxC(DivIm.Grid[2 + nbyteIm * i], DivIm.Grid[1 +  nbyteIm * i],
        DivIm.Grid[0 + nbyteIm * i]);
      lum1 = SigmaIm.MaxC(SubIm.Grid[2 + nbyteIm * i], SubIm.Grid[1 +  nbyteIm * i],
        SubIm.Grid[0 + nbyteIm * i]);
    }
    histoDiv[lum]++;
    histoSub[lum1]++;
  } //=============== end for (i = 0; ... ==============================

  // Displaying the histograms and the row sections:
  Graphics g1 = pictureBox4.CreateGraphics();
  Graphics g = pictureBox5.CreateGraphics();
  Graphics g0 = pictureBox1.CreateGraphics();
  int MaxHisto1 = 0, SecondMax1 = 0;
  int MaxHisto = 0, SecondMax = 0;
  for (i = 0; i < 256; i++)
  {
    if (histoSub[i] > MaxHisto1) MaxHisto1 = histoSub[i];
    if (histoDiv[i] > MaxHisto) MaxHisto = histoDiv[i];
  }
  for (i = 0; i < 256; i++) if (histoSub[i] != MaxHisto1 && histoSub[i] > SecondMax1) SecondMax1 = histoSub[i];
  MaxHisto1 = SecondMax1 * 4 / 3;
  for (i = 0; i < 256; i++) if (histoDiv[i] != MaxHisto && histoDiv[i] > SecondMax) SecondMax = histoDiv[i];
  MaxHisto = SecondMax * 4 / 3;

  Pen redPen = new Pen(Color.Red), yellowPen = new Pen(Color.Yellow),
                      bluePen = new Pen(Color.Blue), greenPen = new Pen(Color.Green);
  SolidBrush whiteBrush = new SolidBrush(Color.White);
  Rectangle Rect1 = new Rectangle(0, 0, 256, 200);
  g1.FillRectangle(whiteBrush, Rect1);
  Rectangle Rect = new Rectangle(0, 0, 256, 200);
  g.FillRectangle(whiteBrush, Rect);
  for (i = 0; i < 256; i++)
  {
    g1.DrawLine(redPen, i, pictureBox4.Height - histoSub[i] * 200 / MaxHisto1, i, pictureBox4.Height);
    g.DrawLine(redPen, i, pictureBox5.Height - histoDiv[i] * 200 / MaxHisto, i, pictureBox5.Height);
  }
  for (i = 0; i < 256; i += 50)
  {
    g1.DrawLine(greenPen, i, pictureBox4.Height - 200, i, pictureBox4.Height);
    g.DrawLine(greenPen, i, pictureBox5.Height - 200, i, pictureBox5.Height);
  }
  if (ROWSECTION)
  {
    y = height * 10 / 19;
    g0.DrawLine(redPen, marginX, marginY + (int)(y * Scale1),
      marginX + (int)(width * Scale1), marginY + (int)(y * Scale1));
    int xold = marginX, xs = 0;
    int yd = 0, ydOld = 256 - DivIm.Grid[0 + width * y];
    int ys = 0, ysOld = 256 - SubIm.Grid[0 + width * y];
    int yg = 0, ygOld = 256 - GrayIm.Grid[0 + width * y];
    int ym = 0, ymOld = 256 - MeanIm.Grid[0 + width * y];
    for (x = 1; x < width; x++)
    {
      xs = marginX + (int)(x * Scale1);
      yd = 256 - DivIm.Grid[x + width * y];
      ys = 256 - SubIm.Grid[x + width * y];
      yg = 256 - GrayIm.Grid[x + width * y];
      ym = 256 - MeanIm.Grid[x + width * y];
      g0.DrawLine(redPen, xold, ymOld, xs, ym);
      g0.DrawLine(yellowPen, xold, ydOld, xs, yd);
      //g0.DrawLine(bluePen, xold, ysOld, xs, ys);
      xold = xs;
      ydOld = yd;
      ysOld = ys;
      ygOld = yg;
      ymOld = ym;
    }
    g0.DrawLine(bluePen, marginX, 256 - Threshold,
    marginX + (int)(width * Scale1, 256 - Threshold);
  }

} //************** end CorrectShading ***********************************

方法FastAverageM从标签为Window in per mille to Width的工具numericUpDown1中获取参数hWind。该方法计算灰度图像的局部平均值GrayIm。请注意，为了在图片照片的情况下获得良好的效果，有时需要选择相当小的Window值。利用稍后描述的项目WFshadBinImpulse可以获得更好的结果。

在具有变量y和x的循环中计算阴影校正图像SubIm和DivIm的内容。

正如你所看到的，方法CorrectShading以两种方式计算带有校正阴影的图像。第一种方法包括将原始图像OrigIm的值除以保存在图像MeanIm中的局部平均亮度，并乘以通过工具numericUpDown2手动指定的值Lightness。第二种方法包括从原始图像OrigIm中减去MeanIm中保存的局部平均亮度，并加上Lightness。使用这两种方法很重要，因为实验表明除法的结果并不总是比减法好。方法CorrectShading以刚刚描述的两种方式计算两幅图像DivIm和SubIm。两个结果都显示给用户，用户可以决定应该保存两个结果中的哪一个。

校正图像的质量取决于两个参数:Window和Lightness。可以通过带有相应标签的工具来更改这些参数，并立即看到结果。参数Lightness指定图像的平均亮度。参数Window定义了计算原始图像平均亮度的滑动窗口的大小。必须选择这个参数，使得它近似等于由于不均匀亮度而出现的原始图像中的斑点的尺寸。例如，前面图 4-1 所示图像的最佳窗口值约为 200 像素。该值应以图像宽度的千分之一为单位进行设置。在这种情况下，从宽度= 400 开始，200 个像素构成千分之 500。用户必须定义Window=500。阴影校正后的表单视图如图 4-1 所示。

当用户对阴影校正的结果满意时，他或她可以保存图像SubIm或图像DivIm，或者对这些图像进行阈值处理。接下来描述阈值处理的过程。

对图像进行阈值处理

为灰度图像选择最佳阈值的常用方法包括使用对应于直方图最小值的灰度值作为阈值。理由如下:具有两个区域的图像的直方图看起来像两列，每个区域具有恒定的灰度级。如果存在噪声，直方图看起来像两座小山，中间有一个山谷(图 4-2 )。

img/474294_1_En_4_Fig2_HTML.png

图 4-2

包含两个区域的图像直方图示例

最佳阈值对应于直方图的局部最小值，但不总是对应于全局最小值。因此，有必要将搜索区域限制为最小值，以保证图像的至少某个预定部分将是黑色或白色。例如，如果您希望图像的至少 5%为黑色，您必须找到这样一个灰度值 minGV ，使得直方图区域的 5%位于 minGV 的左侧。只需考虑 minGV 右侧的最小值。类似地，选择 maxGV 可以保证至少该区域的期望部分是白色的。

有时直方图有不止一个局部最小值，最深的一个并不总是对应于你想要的结果。在这种情况下，对于对应于所有局部最小值的阈值，显示小的二进制图像并选择最佳图像是合适的。

另一种可能性是产生多级图像，每一级对应于两个阈值之间的空间。还必须提供对应于小于第一阈值且大于最后阈值的灰度值的等级。这意味着等于 0 的亮度和等于 255 的亮度也应该被认为是阈值。产生多级图像意味着图像被量化:亮度在阈值TI和TI+1之间的所有像素获得值(TI+TI+1)/2。为了使对局部最小值的搜索更加确定，有必要平滑直方图，因为由于噪声，直方图可能具有许多非常小的局部最小值。图 4-3 显示了一个具有四个等级的图像示例，以及对应于直方图局部最小值的三个不同阈值的阈值处理结果。

img/474294_1_En_4_Fig3_HTML.png

图 4-3

在三个候选阈值中进行选择

人们提出了许多选择最佳阈值的不同方法，但我建议另一种方法:我们的项目提出了手动改变阈值并立即看到结果的可能性。因此，用户可以选择看起来最佳的结果。

一些要进行阈值处理的图像具有强烈的阴影，这意味着局部平均亮度从图像的一部分到另一部分逐渐变化。如果较暗区域的亮侧比较亮区域的暗侧亮，则不存在分隔这两个区域的恒定阈值。在这种情况下，阴影校正非常有用。阴影校正的过程已经描述过了。正如你已经知道的，有两种阴影校正的方法:可以使用减法或除以局部平均值。用户可以在用减法或除法对阴影校正的结果进行阈值处理之间进行选择。我们的项目WFshadingBin在阴影校正后执行阈值处理。它显示了阴影校正图像的直方图。用户可以通过点击直方图(例如，在局部最小值中)来选择阈值，并立即看到结果。所选阈值在直方图中显示为一条蓝色垂直线。图 4-4 再次显示了项目WFshadingBin的形式。

img/474294_1_En_4_Fig4_HTML.jpg

图 4-4

项目的形式WFshadingBin

阈值图像有时在均匀的白色区域中包含暗点。如果发生这种情况，必须增加Window参数。

该项目可用于校正光线不足的图片。在这种情况下，平均窗口的最佳尺寸可以非常小，例如 6 或 9 个像素。这不到宽度的 1 %,宽度通常大于 1，000 像素。因此，我们决定在这个项目中定义平均窗口的大小，单位是千分之一，而不是图像宽度的百分比。

使用WFshadingBin项目进行阈值处理的结果相当令人满意。图 4-5a 显示了 Wikipedia.en 上的文章“Otsu 的方法”中提供的图像的结果，该文章描述了日本科学家 Nobuyuki Otsu 开发的一种相当复杂的方法，用于选择最佳阈值。图 4-5b 显示了通过大津法获得的结果。

img/474294_1_En_4_Fig5_HTML.jpg

图 4-5

维基百科图像的阈值化采用(a) WFshadingBin,和(b) Otsu 的方法

如你所见，结果是相似的。然而，我们的方法具有这样的优点，即用户可以容易且快速地改变阈值并立即看到结果图像。

我们项目中使用的方法CorrectShading在前面的章节中已经描述过了。这里我们展示了源代码的一部分，它画出了显示阈值的蓝线并执行阈值处理。

private void pictureBox5_MouseClick(object sender, MouseEventArgs e)
// Thresholding DivIm
{
  Threshold = e.X;
  Graphics g = pictureBox5.CreateGraphics();
  Pen bluePen = new Pen(Color.Blue);
  g.DrawLine(bluePen, Threshold, 0, Threshold, pictureBox5.Height);
  progressBar1.Visible = true;
  progressBar1.Value = 0;
  int nbyte = DivIm.N_Bits / 8;
  int jump = height / 100;
  for (int y = 0; y < height; y++)
  {
    if (y % jump == jump - 1) progressBar1.PerformStep();
    for (int x = 0; x < width; x++)
    {
      int i = x + width * y;
      if (nbyte == 1)
      {
        if (DivIm.Grid[i] > Threshold) BinIm.Grid[i] = 255;
        else BinIm.Grid[i] = 0;
      }
      else
      {
        if (DivIm.MaxC(DivIm.Grid[2 + 3*i], DivIm.Grid[1 + 3*i],
                                      DivIm.Grid[0 + 3*i]) > Threshold) BinIm.Grid[i] = 255;
        else BinIm.Grid[i] = 0;
      }
      Div_Bitmap.SetPixel(x, y, Color.FromArgb(BinIm.Grid[i], BinIm.Grid[i],
                                               BinIm.Grid[i]));
    }
  }

  pictureBox3.Image = Div_Bitmap;
  Threshold = -1;
} //****************** end pictureBox5_MouseClick ***********************

项目WFshadingBin在改善历史图像的照片方面特别有效。图 4-6 是处理历史图像片段照片的又一个例子。

img/474294_1_En_4_Fig6_HTML.jpg

图 4-6

(a)原始图像，和(b)处理过的图像

五、WFshadBinImpulse项目

我们还开发了项目WFshadBinImpulse，其中结合了阴影校正、阈值处理和脉冲噪声抑制等程序。这种组合对于处理旧图纸照片的图像特别有用(参见图 2-8 )。本项目的形式如图 5-1 所示。

img/474294_1_En_5_Fig1_HTML.jpg

图 5-1

项目的形式WFshadBinImpulse

点击打开图像，用对话框OpenFileDialog启动项目的常规部分。项目的这一部分还定义了七个类别的图像CImage — OrigIm、SigmaIm、GrayIm、MeanIm、ShadIm、BinIm和ImpulseIm。

用户应该通过选择四个选项中的一个来指定图像的种类:用深色线条绘制，或用浅色线条绘制，无绘制 Div，或无绘制 Sub。参数Window in per mille of Width、Lightness、Threshold的初始值以及抑制脉冲噪声的参数将根据这些选项的选择自动设置。用户可以通过相应的numericUpDown工具修正这些值。Window将被指定为原始图像宽度的一部分，单位为千分之一，如前所述。以图像宽度的千分之一而不是直接以像素来指定是必要的，因为窗口的最佳宽度取决于图像的大小，但是用户通常不知道图像有多大。

在用户指定了图像类型后，他或她应该单击底纹。对于阴影校正，通过前面描述的方法FastAverageM计算图像的局部平均亮度。然后通过第四章中描述的CorrectShading方法计算图像SubIm和DivIm。阴影校正后的图像显示在右侧的图片框中。用户可以校正参数的建议值，以获得可能的最佳校正图像。

CorrectShading方法还绘制阴影校正图像的直方图。用户可以通过单击直方图来指定阈值。阈值图像会立即显示。

用于抑制脉冲噪声的方法的参数的建议值也被自动设置。然而，用户应该测试这些参数的一些值。

如果图像是一幅图画，用户应在pictureBox1(原始图像)中围绕图像中不应消除小斑点(如人的眼睛)的部分绘制小矩形。用户应该用鼠标点击矩形的左上角和右下角。矩形以蓝色显示。在脉冲噪声抑制运行后，可以重新定义这些矩形。最多可以画六个矩形。称为maxNumber的矩形的加倍最大数量在Form1的开头定义。

如果用户对阴影校正和阈值处理的结果感到满意，那么他或她可以通过点击脉冲噪声来开始抑制脉冲噪声，该脉冲噪声在阈值图像中被视为小的黑色和白色斑点。如果用户对获得的结果不满意，他或她可以通过相应的numericUpDown工具尝试更改删除暗和删除亮的值。这些值指定了应该删除的斑点中的最大像素数。

如果用户对最终结果满意，那么他或她可以保存结果。需要点击保存结果，选择正确的目录，用扩展名.bmp或.jpg指定结果图像的文件名。

六、边缘检测

在已知的边缘检测方法中，简单的梯度滤波器包含某种图像平滑。索贝尔滤波器是一个常见的例子。它由两个权重矩阵定义，如图 6-1 所示。

img/474294_1_En_6_Fig1_HTML.png

图 6-1

Sobel 滤波器的权重

该滤波器计算实际像素的滑动 3 × 3 邻域中的灰度值与图 6-1 所示的相应权重 WH 和 WV 的乘积的两个和 SH 和 SV ，并返回绝对值| SH | + | SV |。该值大于给定阈值的像素属于边缘。由于边缘太厚，结果一般不令人满意。

其他已知的边缘检测方法是拉普拉斯 Kimmel (2003)和 Canny (1986)滤波器的零交叉。下面我将把这些方法与我们新的相当简单有效的方法进行比较。

拉普拉斯算子

一种有效的边缘检测方法是拉普拉斯的零交叉。拉普拉斯算子在数学上被定义为二阶偏导数的和:

Lap( F ( x 、和)=(x，【t】

*因为数字图像是不可微的，所以有必要用有限差分来代替导数。

定义 1: 表达式 D 1 ( x ，δx，F)=F(x+δx)-F(x)称为 F ( x )的第一差就

**定义二:**表达式Dn(x，δx，F)=D1(x，δx，Dn-1(x相对于 x 和 y 的第二偏差值等于

d(x， x，f=(x-【t】

d(y， y，f=(*，**

*并且δx =δy= 1 的有限拉普拉斯算子等于

2【f】(x、和=f(x–1

*可以用图 6-2 所示的滤波器在数字图像中进行计算。

img/474294_1_En_6_Fig2_HTML.png

图 6-2

用于计算有限拉普拉斯算子的滤波器

将该滤波器应用于具有灰度值 F ( x ， y )的数字图像产生以下值:

2【f】(x、和=f(x–1

*等式 6.1 可以改写如下:

2【f】(x、和=f(x–1 和+f(x、y+1)-5(x、**

**= 5 M ( x ， y ) - 5 F ( x ，y)= 5(M(x，y)-F(x，y)(6.2)

其中 M ( x ， y )为坐标为( x ， y )的像素邻域内五个像素的 F ( x ， y )的平均值。也可以通过从 P 邻域的局部均值中减去 P 的灰度值，计算出像素 P = ( x ， y )中的有限拉普拉斯算子。当使用更大的邻域时，最好使用快速平均滤波器(第二章)。

过零方法

过零方法基于检测拉普拉斯变换符号的位置。考虑图 6-3 所示的一行数字灰度图像的横截面。

img/474294_1_En_6_Fig3_HTML.png

图 6-3

一行数字图像的横截面

拉普拉斯算子在边的一侧有正值，在另一侧有负值。与远离边缘的位置上的拉普拉斯值相比，边缘附近的拉普拉斯绝对值较大。使用拉普拉斯的负值更实际——即值( F ( x ，y)——平均值)而不是值(平均值——F(x，y)——因为这些值在更大的灰度值下更大。

拉普拉斯变换其符号的位置被称为零交叉。这些位置中的一些对应于边缘的位置。零交叉总是位于两个像素之间。因此，边缘位置的这些指示总是非常微弱。

拉普拉斯曲线的过零点是闭曲线吗？

在一些图像处理的教科书中提出的拉普拉斯算子的一个优点是，拉普拉斯算子产生的边缘总是封闭的直线。然而，只有当我们设定拉普拉斯值的阈值时，这才是正确的；例如，如果我们用 1 代替拉普拉斯的正值，用 0 代替负值。然而，在这样的图像中具有值 1 的像素组的边界是闭合曲线的事实是二进制图像的优点，而不是拉普拉斯的优点。这些边缘如图 6-4b 所示。用等于 0 的阈值二值化的拉普拉斯图像中的大多数边缘是不相关的。

img/474294_1_En_6_Fig4_HTML.jpg

图 6-4

(a)原始图像，以及(b)用阈值 0 二值化的拉普拉斯算子；边缘显示为白线

如你所见，图 6-4b 中的绝大多数白线完全没有意义，不能作为边缘。这些就是所谓的无关过零点。

图 6-5 显示了拉普拉斯不相关零交叉的解释。黑线显示原始图像的灰度值。绿线代表灰度值的局部均值，红线是拉普拉斯算子。红色箭头表示不相关的过零事件。在这些位置的拉普拉斯值很小，而在相关交叉点的拉普拉斯绝对值很大。

img/474294_1_En_6_Fig5_HTML.png

图 6-5

拉普拉斯算子的不相关过零点(红色箭头)

我们将在下一节解释如何消除不相关的过零事件。

如何消除不相关的交叉

不相关的交叉可以通过两个阈值与相关的交叉区分开:只有从大于正阈值的拉普拉斯值到小于负阈值的拉普拉斯值的转变才应该被识别为相关的边缘。大多数不相关的边缘会消失，但一些相关的边缘会因图像中的噪声而中断。应对图像进行滤波以降低噪声，例如使用最简单的均值滤波器(第章第二部分)。图像应该被过滤两次:用一个小的和一个大的滑动窗口。这两个滤波图像之间的差异可以被认为是拉普拉斯算子的良好近似。

考虑图 6-6b ，图中一行的灰度值显示为一条蓝线。拉普拉斯值显示为绿线。黑线是 x 轴；红线显示用于区分相关和不相关过零事件的阈值。如果零交叉位于两个这样的像素之间，其中一个像素的拉普拉斯值大于正阈值，而另一个像素的拉普拉斯值小于负阈值，则零交叉是相关的。

img/474294_1_En_6_Fig6_HTML.jpg

图 6-6

(a)用两个阈值(白线)找到的边缘，以及(b)解释

图 6-7 解释了如何消除不相关的过零事件。

img/474294_1_En_6_Fig7_HTML.png

图 6-7

拉普拉斯值的一个例子

使用拉普拉斯算子之前的降噪

使用均值滤波器降噪时，边缘会变得模糊。结果，拉普拉斯的正值和负值之间的差异变小；拉普拉斯值逐渐变化。这使得很难区分相关和不相关的过零事件。我们建议使用 sigma 滤波器(第章第二部分)，而不是使用小窗口求平均值。那么边缘保持陡峭；拉普拉斯算子的绝对量变得更大。

图 6-8a 中的黑线显示了通过小窗口平均过滤的灰度值。在图 6-8b 中，黑线代表用西格玛滤波器过滤的灰度值。对于两者，绿线代表具有较大窗口的局部平均值，红线显示拉普拉斯值。

img/474294_1_En_6_Fig8_HTML.png

图 6-8

使用(a)平均和(b) sigma 滤波器滤波后的灰度和拉普拉斯值

数字化和极值滤波过程中的模糊

在图像的数字化过程中，即使是理想的边缘也会变得有些模糊。其原因是，投影到一组光敏元件(例如，CCD、电荷耦合器件、各种设备中使用的电子光传感器，包括数码相机)的图像中的亮区和暗区之间的边界偶尔会落在元件的中心附近，如图 6-9 所示。

img/474294_1_En_6_Fig9_HTML.jpg

图 6-9

CCD 矩阵由一个暗区和一个亮区照明

这个元素获得一个光量，我们称之为*中间值，*位于适合亮区和暗区的值之间。相应的像素位于边缘的中间，并且拉普拉斯算子获得小的值或者甚至零值。检测到的边缘出现间隙，如图 6-10 所示。

img/474294_1_En_6_Fig10_HTML.png

图 6-10

灰度值(黑线)、局部平均值(绿线)和拉普拉斯算子(红线)

我们建议使用极值过滤器来避免中间值导致的拉普拉斯间隙。这个应用于灰度图像的滤波器在一个小的滑动窗口中计算最大和最小灰度值。它还计算中心像素的灰度值与最大值和最小值之间的差异，并决定这两个值中的哪一个更接近中心值。更接近的值被分配给输出图像中的中心像素。边缘变得清晰，拉普拉斯间隙消失。这里是灰度图像的极端过滤器的源代码。

int CImage::ExtremVar(CImage &Inp, int hWind)
{ N_Bits=8; width=Inp.width; height=Inp.height;
  Grid=new unsigned char[width*height];
  int hist[256];
  for (int y=0; y<height; y++) // ================================
  { int gv, y1, yStart=__max(y-hWind,0), yEnd=__min(y+hWind,height-1);
     for (int x=0; x<width; x++) //============================
    { if (x==0) //-------------------------------------------------------
      { for (gv=0; gv<256; gv++) hist[gv]=0;
         for (y1=yStart; y1<=yEnd; y1++)
         for (int xx=0; xx<=hWind; xx++)
                                  hist[Inp.Grid[xx+y1*width]]++;
      }
      else
      { int x1=x+hWind, x2=x-hWind-1;
         if (x1<width)for (y1=yStart; y1<=yEnd; y1++)
             hist[Inp.Grid[x1+y1*width]]++;
         if (x2>=0)
          for (y1=yStart; y1<=yEnd; y1++)
          { hist[Inp.Grid[x2+y1*width]]--;
             if (deb && hist[Inp.Grid[x2+y1*width]]<0) return -1;
          }
      } //---------------- end if (x==0) ---------------------------------
      int gMin=0, gMax=255;
      for (gv=gMin; gv<=gMax; gv++)
        if (hist[gv]>0) { gMin=gv; break; }
      for (gv=gMax; gv>=0; gv--)
        if (hist[gv]>0) { gMax=gv; break; }
      if (Inp.Grid[x+width*y]-gMin<gMax-Inp.Grid[x+width*y])
         Grid[x+width*y]=gMin;
      else Grid[x+width*y]=gMax;
    } //=============== end for (int x... ================
  } //================ end for (int y... =================
  return 1;
} //******************* end ExtremVar *********************

这里是彩色和灰度图像的通用极端方法的源代码。

public int ExtremLightUni(CImage Inp, int hWind,Form1 fm1)
/* The extreme filter for color or grayscale images with variable hWind. The filter finds in the (2*hWind+1)-neighborhood of the actual pixel (x,y) the color "Color1" with minimum and the color "Color2" with the maximum lightness. "Color1" is assigned to the output pixel if its lightness is closer to the lightness of the central pixel than the lightness of "Color2". --*/
{
 byte[] CenterColor = new byte[3], Color = new byte[3], Color1 =
  new byte[3], Color2 = new byte[3];
  int c, k, nbyte = 3, x;
  if (Inp.N_Bits == 8) nbyte = 1;
  for (int y = 0; y < height; y++) // ====================================
  {
    if (y % jump == jump - 1) fm1.progressBar1.PerformStep();
    for (x = 0; x < width; x++) //======================================
    {
      for (c = 0; c < nbyte; c++) Color2[c] = Color1[c] = Color[c] =
                                 CenterColor[c] = Inp.Grid[c + nbyte * (x + y * width)];
      int MinLight = 1000, MaxLight = 0;
      for (k = -hWind; k <= hWind; k++) //==============================
      {
        if (y + k >= 0 && y + k < height)
          for (int i = -hWind; i <= hWind; i++) //==========================
          {
            if (x + i >= 0 && x + i < width) // && (i > 0 || k > 0))
            {
              for (c = 0; c < nbyte; c++)
                           Color[c] = Inp.Grid[c + nbyte * (x + i + (y + k) * width)];
              int light;
              if (nbyte == 3) light= MaxC(Color[2], Color[1], Color[0]);
              else light = Color[0];
              if (light < MinLight)
              {
                MinLight = light;
                for (c = 0; c < nbyte; c++) Color1[c] = Color[c];
              }
              if (light > MaxLight)
              {
                MaxLight = light;
                for (c = 0; c < nbyte; c++) Color2[c] = Color[c];
              }
            }
          } //=============== end for (int i... =======================
      } //=================== end for (int k... ======================
      int CenterLight = MaxC(CenterColor[2], CenterColor[1], CenterColor[0]);
      int dist1 = 0, dist2 = 0;
      dist1 = CenterLight - MinLight;
      dist2 = MaxLight - CenterLight;
      if (dist1 < dist2)
        for (c = 0; c < nbyte; c++)  Grid[c + 3 * x + y * width * 3] = Color1[c]; // Min

      else
        for (c = 0; c < nbyte; c++)  Grid[c + 3 * x + y * width * 3] = Color2[c]; // Max

    } //================== end for (int x... ==========================
  } //==================== end for (int y... ==========================
  //fm1.progressBar1.Visible = false;
  return 1;
} //********************** end ExtremLightUni ***************************

考虑在使用 sigma 滤波器和极值滤波器之后通过拉普拉斯过零方法(阈值= 10)检测的边缘的例子。图 6-11 显示了原始图像的灰度值(图 6-11a )以及用滤波器对图像进行连续步骤处理后的图像:用西格玛滤波器(图 6-11b )、用极值滤波器(图 6-11c )和带边缘的拉普拉斯算子(图 6-11d )处理后的图像。拉普拉斯正值显示为红色，负值显示为蓝色，过零点显示为白色线条。

img/474294_1_En_6_Fig11_HTML.jpg

图 6-11

具有拉普拉斯零交叉的边缘检测的例子:(a)原始图像，(b) sigma 滤波的，(c)极端滤波的，以及(d)具有边缘的拉普拉斯算子(白线)

如你所见，边缘检测相当成功。然而，在下一节中解释了边缘中的一些间隙。

拉普拉斯零交叉法的基本误差

拉普拉斯算子具有一个基本性质，即它的过零点序列在三个或四个序列相交的点的邻域中一定有间隙。图 6-12 提供了一个例子。图 6-12 显示了人工图像(a)、通过过零方法检测到的边缘(b)、理想边缘(c)以及拉普拉斯符号的值(d)。红色像素是拉普拉斯值为正值的像素，蓝色像素是负值的像素。过零点在图 6-12d 中显示为白线。

img/474294_1_En_6_Fig12_HTML.png

图 6-12

(a)带有噪声的人工图像，σ= 5；(b)在 sigma 和极端滤波之后图像上的过零点；理想的边缘；以及(d)拉普拉斯算子上的过零点(红色是正值)

注意图 6-12 中由绿色箭头标记的基本误差导致的边缘缺失。这些被认为是基本误差，因为不可能消除它们。出现这些误差是因为拉普拉斯算子是局部平均值和实际灰度值之间的差。零交叉位于拉普拉斯算子的正值像素和负值像素之间。这些像素中的一个像素的灰度值大于局部平均值，而另一个像素的灰度值小于局部平均值。然而，如果在滑动窗口中有两个以上不同的灰度值，则局部平均值只能位于两个灰度值之间。第三灰度值不能有零交叉。因此，边缘有一个间隙。由于这些间隙，通过拉普拉斯零交叉产生的边缘不包含分叉。

这些错误非常重要，因为边通常用于生成表示具有恒定颜色的区域边界的多边形。缺少分叉会使这些多边形不完整，因为它们没有正确描述区域的边界。*********

avatar
文章
阅读
粉丝