现代图像处理算法教程（二）
七、一种新的边缘检测方法
如图 6-11 的示例所示，sigma 和极值滤波器的实施提高了拉普拉斯过零的质量，从而也提高了边缘的质量。需要强调的是，这两种滤波器极大地提高了图像的质量，使得应用最简单的边缘检测方法成为可能。这是灰度值梯度的二值化:
度**(、和)=【gv】*、和【y】*****
度和 ( x 、和 ) = GV x 、和——gv
这里 GV ( x ， y )是像素的灰度值( x ， y )。必须将| Gradx(x， y )|和| Grady(x， y )|的值与一个阈值进行比较。如果| Gradx(x， y )|的值大于阈值，则有一小块边缘位于像素( x ， y )和(x—1，y )之间。对于| Grad*(x， y )|:如果该值大于阈值，则边缘在像素( x ， y )和( x ，*y—1)之间运行。
*在彩色图像的情况下，我们使用像素的色差( x 1 、 y 1 )和( x 2 、 y 2 )。我们称这种方法为二值化梯度。这是像素的颜色通道的强度差的绝对值之和:
difcolor =∑color(I、x、和)–color(I，，
其中 color( i ， x ， y 1 )是像素的第 i 个颜色通道的强度( x 1 ， y 1 )，对于 color( i ，x2
图 7-1 显示了使用拉普拉斯零交叉(a)和二值化梯度(b)的边缘检测结果的比较示例。你可以看到这些图像几乎是一样的。
图 7-1
拉普拉斯零交叉(a)和二值化梯度(b)
比较绿色箭头指示的位置:通过二值化梯度找到的边缘没有由于过零方法的基本误差而在过零边缘出现的间隙。当使用色差时，边缘的进一步处理变得更加精确。
用于编码边缘的装置
对边缘进行编码不是一个简单的问题:例如，如果发现一对色差很大的像素( x 、 y )和(x–1、 y )，那么必须在包含边缘的图像中标记一个指示边缘位置的像素。假设我们决定用边缘来标记图像中的像素( x ， y )。在像素对( x 、 y )和( x 、y–1)的情况下，我们决定标记像素( x 、 y )。然后，在图 7-2 所示的具有不同边缘元素的三种情况下，阈值等于 40，同一像素，即阴影像素，将被标记，并且不可能在这三种情况之间进行区分。
图 7-2
边缘不同的三幅图像
如果我们决定用大的色差来标记一对相邻像素中的另一个像素，那么在其他图像的情况下也会发生相同的情况。解决这个问题的唯一可能的方法是引入另一种编码边缘的结构，即包含像素和位于两个相邻像素之间的边缘元素的不同元素的结构。这种结构被称为抽象细胞复合体(ACC 参见 Kovalevsky，1989，2008)。
抽象细胞复合体的概念
我们在这里将数字平面视为一个二维的细胞复合体，而不是一组像素(见图 7-3 )。因此，我们的数字平面除了像素之外还包含裂缝和点，它们被认为是小正方形。裂缝是这些方块的边；点是裂缝的端点，因此是像素的角。点是零维单元，裂缝是一维单元，像素是二维单元。
图 7-3
包含阴影子集裂纹边界的小型二维复合体的示例
将平面视为抽象的细胞复合体有许多优点:不再有第八章中描述的连通性悖论，子集的边界变成了面积为零的细曲线，区域的边界和它的补的边界是相同的，等等。数字曲线尤其是数字直线的定义和处理变得更加简单和清晰。从图像的经济编码和精确重建的角度来看，最重要的优点是能够通过极其简单和快速的算法(参见 Kovalevsky，(1990))来填充裂纹边界的内部，这在将边界表示为像素组时是不能应用的。
本节的其余部分包含对本演示重要的拓扑概念的简短总结。更多细节和拓扑基础请参考 Kovavlevsky (2008)。熟悉细胞复合体的读者可以跳过这一节的其余部分。
单元之间存在一种约束关系:一个较低维度的单元可能会约束一些较高维度的单元。回头参考图 7-3 所示的小型二维复合体的例子。像素被表示为正方形的内部，裂缝被表示为正方形的边，点(即 0 单元)是裂缝的端点，同时也是像素的角。
现在让我们介绍一些我们在续集中需要的概念。复合体的子集 S 的边界裂缝是将属于 S 的像素与不属于 S 的另一个像素分开的裂缝。图 7-3 中阴影子集的边界裂纹绘制为粗线。子集 S 的边界(也称为裂纹边界)是 S 的所有边界裂纹以及这些裂纹的所有端点的集合。边界不包含像素，因此是一个面积为零的稀疏集合。边界的连通子集称为边界曲线。关于连通性的概念，请参考 Kovalevsky (1989)。
我们认为数字平面是一个*笛卡尔二维复合体；*即作为平面坐标轴的两个一维复形的笛卡尔积。 x 坐标是一行的编号； y 坐标是行号。我们在这里使用计算机图形的坐标系统；即正 x 轴从左向右运行，正 y 轴从上向下运行。
为了保存detected边缘，我们需要一种特殊的 ACC:Cartesian two-dimensional ACC``.使用这种 ACC，引入细胞坐标成为可能。与数字图像的情况完全一样， x 坐标是一列的编号，而 y 坐标是一行的编号。
对于数字图像，我们使用通常的坐标，其中 x 从 0 变化到宽度–1， y 从 0 变化到高度–1。我们称这些坐标为标准。ACC 中代表宽 × 高像素图像的单元格坐标面积更大: x 坐标从 0 变为 2* 宽， y 从 0 变为 2* 高。因此，表示宽度为×高度为像素的图像的 ACC 的大小具有(2 宽度为 + 1) × (2* 高度为 + 1)个单元的大小。我们将图像中代表给定数字图像 ACC 的细胞坐标称为*组合坐标。我们称包含 ACC 的图像为组合坐标中的图像。**
*注意，像素的组合坐标 x 和 y 都是奇数，而垂直裂纹的 x 坐标是偶数，其 y 坐标是奇数。在水平裂纹的情况下，情况相反:它的 x 坐标是奇数，而它的 y 坐标是偶数。一个点的两个组合坐标都是偶数。
在我们的一些项目中，我们处理检测到的边缘。例如，为了分析对象的形状，用多边形来近似检测到的边缘是相当方便的。要找到多边形，我们必须沿着边走。与通常只包含像素的图像相比，在组合坐标中跟踪图像的边缘要简单和方便得多。在我们的项目中，组合坐标中的图像被称为CombIm。
一种简单的边缘编码方法
该方法使用包含单元复合体的图像CombIm,该单元复合体的尺寸对应于已处理图像的尺寸:宽度CombIm等于 2 *宽度+ 1，高度等于 2 *高度+ 1，其中宽度和高度是已处理图像的尺寸。添加+ 1 是必要的，以便为右边的点和CombIm的底部边界留出位置，这有时对处理很重要。
方法LabelCells一行接一行地读取由极值滤波器处理的图像，并为每个像素( x，y )测试该像素( x，y )与其相邻像素之一(x-1y)和( *x，y-*1)的颜色的绝对差值是否大于给定的阈值。阈值可以由用户定义，用户应该知道在低阈值下，边缘裂纹的数量变大，并且边缘的某些部分变厚的概率也很大。阈值的正确值应该通过实验找到。
如果像素( x，y )和(x—1，y )的色差的绝对值大于阈值，则位于这些像素之间的垂直裂纹获得标记1。类似地，如果像素( x，y )和( *x，y—*1)的色差大于阈值，则相应的水平裂缝被标记为1。
与裂纹相关的两个点的标记将同时增加。如果一个点与许多裂纹相关联(最多可以有四个关联裂纹)，则该点的标签会增加许多倍。有限地，一个点的标号的值等于与该点相关的裂纹的数量。该信息可以在编码边缘的处理过程中使用。
下面是LabelCells方法的源代码。
public void LabelCells(int th, CImage Image3)
{
int difH, difV, nbyte, NXB = Image3.width, x, y;
byte Lab = 1;
if (Image3.N_Bits == 24) nbyte = 3;
else nbyte = 1;
for (x = 0; x < width * height; x++) Grid[x] = 0;
byte[] Colorh = new byte[3];
byte[] Colorp = new byte[3];
byte[] Colorv = new byte[3];
for (y = 0; y < height; y += 2)
for (x = 0; x < width; x += 2) // through all points
Grid[x + width * y] = 0;
for (y = 1; y < height; y += 2)
for (x = 1; x < width; x += 2) // through the right and upper pixels
{
if (x >= 3) //-- vertical cracks: abs.dif{(x/2, y/2)-((x-2)/2, y/2)} -------
{
for (int c = 0; c < nbyte; c++)
{
Colorv[c] = Image3.Grid[c + nbyte *((x-2)/2)+nbyte*NXB*(y/2)];
Colorp[c] = Image3.Grid[c + nbyte * (x / 2) + nbyte *NXB*(y/2)];
}
if (nbyte == 3) difV = ColorDifAbs(Colorp, Colorv);
else difV = Math.Abs(Colorp[0] - Colorv[0]);
if (difV < 0) difV = -difV;
if (difV > th)
{
Grid[x - 1 + width * y] = Lab; // vertical crack
Grid[x - 1 + width * (y - 1)]++; // point above the crack;
Grid[x - 1 + width * (y + 1)]++; // point below the crack
}
} //------------------------ end if (x>=3) --------------------------
if (y >= 3) //--- horizontal cracks: abs.dif{(x/2, y/2)-(x/2, (y-2)/2)} ---
{
for (int c = 0; c < nbyte; c++)
{
Colorh[c] = Image3.Grid[c + nbyte *(x/2)+nbyte*NXB*((y-2)/2)];
Colorp[c] = Image3.Grid[c + nbyte * (x / 2) + nbyte *NXB*(y/2)];
}
if (nbyte == 3) difH = ColorDifAbs(Colorp, Colorh);
else difH = Math.Abs(Colorp[0] - Colorh[0]);
if (difH > th)
{
Grid[x + width * (y - 1)] = Lab; // horizontal crack
Grid[x - 1 + width * (y - 1)]++; // point left of crack
Grid[x + 1 + width * (y - 1)]++; // point right of crack
}
} //------------------------ end if (y>=3) --------------------------
} //================= end for (x=1;... =====================
} //******************* end LabelCells *************************
二值化梯度方法的改进
二值化梯度方法的最简单版本有一个缺点:对于灰度值或颜色的差异，用小阈值产生的边缘有时太厚。如果图像在某些位置模糊不清，就会出现这种情况，从而产生强度渐变。在斜坡宽度大于 4 个像素的情况下(这种情况很少发生)，通常使用的窗口为 5 × 5 像素的极值滤波器(前面源代码中的参数hWind等于 2)无法完全消除斜坡。然后边缘可以模糊。
为了消除这个问题，我们开发了一个改进版本的方法LabelCells，我们称之为LabelCellsSign。这种方法的目标类似于 Canny(1986)的无最大值抑制方法的目标，但是实现的解决方案非常不同。
让我们首先考虑应用于数字图像的梯度概念。人们不能使用梯度的经典概念(在二维情况下)定义为具有两个分量的向量，这两个分量是强度相对于坐标 x 和 *y，*的导数，因为数字图像不是欧几里德空间。需要用有限差分来代替导数:颜色的强度 I 相对于 x 的导数由 I ( x + 1，y)-I(x， y )和强度 I 相对于 y 的导数亮度 I ( x ， y )是像素的属性( x ， y )，其中 x 是列的编号， y 是行的编号。但是， I ( x + 1，y)–I(x， y )和 I ( x ，y+1)–I(在这种情况下，ACC 的概念变得非常有用:两个差异中的每一个都可以分配给位于一对像素之间的相应裂纹(一维单元)。然后，应将梯度分配给与两个裂纹相关的点(零维单元)。
首先考虑灰度图像的情况。有必要为图像中的每个位置计算灰度值的梯度，标记梯度的绝对值大于为边缘检测定义的阈值的位置，找到沿着平行于梯度的线的标记位置的连接子集，并找到该子集中梯度的绝对值达到最大值的位置。这个位置属于边缘。
很明显，在梯度的绝对值达到最大值的位置，其分量等于相邻像素中灰度值的差也达到最大值。我们感兴趣的是属于边缘的裂缝。因此，没有必要使用沿着梯度方向在连通子集中寻找位置的复杂过程。测试一行中的所有相邻像素对，以找到灰度值的差大于阈值的像素对，并在这些像素对中寻找差的最大值就足够了。这样，边缘的垂直裂缝将被发现。通过测试一列中的相邻像素对，可以以类似的方式找到边缘的水平裂缝。
然而，应该考虑梯度的方向或灰度值差异的符号。如果我们只考虑绝对值，那么图 7-4 中用红色和绿色箭头表示的情况无法区分，将会找到绝对值的单个最大值。但是，在这种情况下，必须找到两条边，因为有两条斜坡:一条斜坡在暗条上方亮度递减，另一条在暗条下方亮度递增。因此，我们必须考虑灰度值差异的符号，并寻找正差异的最大值和负差异的最小值。
图 7-4
灰度图像中亮度梯度的示例
让我们首先考虑灰度图像情况下的方法LabelCellsSign，尽管它也适用于彩色图像。
该方法实现了某种有限状态机。它使用一个变量State，其值对应于方法的不同状态。在方法开始时，State获得零值。
该方法属于类CImage并被称为CombIm.LabelCellsSign(threshold, InputImage)，其中InputImage是由极端过滤器计算的图像。CombIm是包含用于编码边缘的细胞复合体的图像。因此，其尺寸为(如之前描述的方法LabelCells的情况)2*width + 1和2*height + 1，其中width和height是InputImage的尺寸。
方法LabelCellsSign使用两对for循环，其变量x和y扫描包含在图像CombIm中的 ACC 的二维单元的所有位置，以保存边缘。
第一对线圈用于检测垂直裂缝。变量为y的外部循环遍历从1到height - 1的所有奇数值。带有变量x的内部循环从值x = 3开始，因为它使用值x - 2，该值必须保留在图像区域CombIm。在x循环的开始，计算输入图像站的两个相邻像素((x - 1) / 2, y / 2)和(x / 2, y / 2) ( x和y是图像CombIm中的坐标)的颜色差difV。在灰度图像的情况下，它是通过减法直接计算的。然而，在彩色图像的情况下，它是用方法ColorDifSign计算的，这将在后面解释。与方法LabelCells不同，这里计算的差值不是绝对值，而是一个符号。
变量Inp表示前面提到的两个像素的色差与边缘检测指定的Threshold的关系:如果差值大于Threshold，则Inp等于1，其他地方等于0。变量State和Inp组成控制变量Contr，等于3*State + Inp控制开关指令。变量xStartP包含大于Threshold的色差序列开始的x的值。类似地，变量xStartM包含x的值，其中color差小于负阈值-Threshold的序列开始。这里是switch指令的伪代码。
switch (Contr)
{
case 4: if (x > xStartP && difV > maxDif)
{ maxDiv = difV;
xopt = x;
}
break;
case 3: label the vertical crack at (xopt - 1, y) and its end points;
State = 0;
break;
case 2: label the vertical crack at (xopt - 1, y) and its end points;
minDif = difV;
xStartM = x;
State = -1;
break;
case 1: maxDiv = difV; xopt = x; xStartP = x;
State = 1;
break;
case 0: break;
case -1: maxDiv = difV; xopt = x; xStartP = x;
State = -1;
break;
case -2: label the vertical crack at (xopt - 1, y) and its end points;
maxDif = difV;
xStartP = x;
State = 1;
case -3: label the vertical crack at (xopt - 1, y) and its end points;
State = 0;
break;
case -4: if (x > xStartM && difV < minDif)
{ minDiv = difV;
xopt = x;
}
break;
} //:::::::::::::::::::::: end switch :::::::::::::::::::::::::::::::.:::::
在行首，变量State、Inp和Contr等于0。如果Inp变为值1，那么Contr也变为1，并且在case 1:中，最大值maxDif被设置为difV , xopt和xStartP被设置为x，并且State被设置为1。数值xStartP是色差大于Threshold的裂纹序列的起点坐标。如果Inp保留1，则Contr变为4，因为State == 1和difV的最大值是沿着Inp保留1的裂纹顺序计算的。
如果Inp变为0，那么Contr变为3，在最大位置xopt - 1标记一条垂直裂纹。这是像素序列中唯一一个difV大于Threshold的位置，在此标记了垂直裂缝。因此，在运行值y处，垂直裂纹序列变细。
水平裂纹及其端点的标记发生在第二对for环中，其中外环为x环，内环为y环。switch指令看起来与垂直裂纹的指令相似，但是变量difV被替换为difH，变量xopt被替换为yopt，变量xStartP被替换为yStartP，变量xStartM被替换为yStartM。下面是方法LabelCellsSign的源代码:
public int LabelCellsSign(int th, CImage Extrm)
{
int difH, difV, c, maxDif, minDif, nByte, NXB = Extrm.width, x, y, xopt, yopt;
int Inp, State, Contr, xStartP, xStartM, yStartP, yStartM;
if (Extrm.N_Bits == 24) nByte = 3;
else nByte = 1;
for (x = 0; x < width * height; x++) Grid[x] = 0;
byte[] Colorp = new byte[3], Colorh = new byte[3], Colorv = new byte[3];
maxDif = 0; minDif = 0;
for (y = 1; y < height; y += 2) //====== vertical cracks ==============
{
State = 0;
xopt = -1;
xStartP = xStartM = -1;
for (x = 3; x < width; x += 2) //==============================
{
for (c = 0; c < nByte; c++)
{
Colorv[c] = Extrm.Grid[c + nByte * ((x - 2) / 2) + nByte*NXB*(y / 2)];
Colorp[c] = Extrm.Grid[c + nByte * (x / 2) + nByte * NXB * (y / 2)];
}
if (nByte == 3) difV = ColorDifSign(Colorp, Colorv);
else difV = Colorp[0] - Colorv[0];
if (difV > th) Inp = 1;
else
if (difV > -th) Inp = 0;
else Inp = -1;
Contr = State * 3 + Inp;
switch (Contr) //:::::::::::::::::::::::::::::::::::::::::::::::::
{
case 4:
if (x > xStartP && difV > maxDif)
{
maxDif = difV;
xopt = x;
}
break;
case 3:
Grid[xopt - 1 + width * y] = 1; // vertical crack
Grid[xopt - 1 + width * (y - 1)]++; // point above
Grid[xopt - 1 + width * (y + 1)]++; // point below
State = 0;
break;
case 2:
Grid[xopt - 1 + width * y] = 1; // vertical crack
Grid[xopt - 1 + width * (y - 1)]++; // point above
Grid[xopt - 1 + width * (y + 1)]++; // point below
minDif = difV;
xopt = x;
xStartM = x;
State = -1;
break;
case 1: maxDif = difV; xopt = x; xStartP = x; State = 1; break;
case 0: break;
case -1: minDif = difV; xopt = x; xStartM = x; State = -1; break;
case -2:
Grid[xopt - 1 + width * y] = 1; // vertical crack
Grid[xopt - 1 + width * (y - 1)]++; // point above
Grid[xopt - 1 + width * (y + 1)]++; // point below
maxDif = difV;
xopt = x;
xStartP = x;
State = 1;
break;
case -3:
Grid[xopt - 1 + width * y] = 1; // vertical crack
Grid[xopt - 1 + width * (y - 1)]++; // point above
Grid[xopt - 1 + width * (y + 1)]++; // point below
State = 0;
break;
case -4:
if (x > xStartM && difV < minDif)
{
minDif = difV;
xopt = x;
}
break;
} //:::::::::::::::::::::: end switch ::::::::::::::::::::::::::::::
} //================ end for (x=3;... =====================
} //================= end for (y=1;... ======================
for (x = 1; x < width; x += 2) //=== horizontal cracks ===============
{
State = 0;
minDif = 0; yopt = -1; yStartP = yStartM = 0;
for (y = 3; y < height; y += 2) //=============================
{
for (c = 0; c < nByte; c++)
{
Colorh[c] = Extrm.Grid[c + nByte * (x / 2) + nByte*NXB*((y - 2) / 2)];
Colorp[c] = Extrm.Grid[c + nByte * (x / 2) + nByte * NXB * (y / 2)];
}
if (nByte == 3) difH = ColorDifSign(Colorp, Colorh);
else difH = Colorp[0] - Colorh[0];
if (difH > th)
Inp = 1;
else
if (difH > -th) Inp = 0;
else Inp = -1;
Contr = State * 3 + Inp;
switch (Contr) //:::::::::::::::::::::::::::::::::::::::::::::::::
{
case 4:
if (y > yStartP && difH > maxDif)
{
maxDif = difH;
yopt = y;
}
break;
case 3:
Grid[x + width * (yopt - 1)] = 1; // horizontal crack
Grid[x - 1 + width * (yopt - 1)]++; // left point
Grid[x + 1 + width * (yopt - 1)]++; // right point
State = 0;
break;
case 2:
Grid[x + width * (yopt - 1)] = 1; // horizontal crack
Grid[x - 1 + width * (yopt - 1)]++; // left point
Grid[x + 1 + width * (yopt - 1)]++; // right point
yopt = y;
State = -1;
break;
case 1: maxDif = difH; yopt = y; yStartP = y; State = 1; break;
case 0: break;
case -1: minDif = difH; yopt = y; yStartM = y; State = -1; break;
case -2:
Grid[x + width * (yopt - 1)] = 1; // horizontal crack
Grid[x - 1 + width * (yopt - 1)]++; // left point
Grid[x + 1 + width * (yopt - 1)]++; // right point
yopt = y;
State = 1;
break;
case -3:
Grid[x + width * (yopt - 1)] = 1; // horizontal crack
Grid[x - 1 + width * (yopt - 1)]++; // left point
Grid[x + 1 + width * (yopt - 1)]++; // right point
State = 0;
break;
case -4:
if (y > yStartM && difH < minDif)
{
minDif = difH;
yopt = y;
}
break;
} //:::::::::::::::::::::: end switch ::::::::::::::::::::::::::::::
} //================ end for (y=3;... =====================
} //================= end for (x=1;... ======================
return 1;
} //******************** end LabelCellsSign **********************
我们之前解释了这种方法在灰度图像的情况下是如何工作的。现在让我们解释一下它在彩色图像中的作用。
我们应该首先考虑应用于彩色数字图像的渐变概念。首先，在彩色图像的情况下，不存在单一的梯度，而是三个梯度，每个梯度对应于三原色之一:红色、绿色和蓝色。这些梯度中的每一个都由相邻像素中相应颜色的两个强度差来定义:作为 I ( x + 1，y–I(x， y )的 x 分量和作为 I ( x ，y的分量亮度 I ( x ， y )是像素的属性( x ， y )，其中 x 是列的编号， y 是行的编号。但是，一个 I ( x + 1，y)–I(x， y )或者 I ( x ，y+1)–I(x)在这种情况下，ACC 的概念是非常有用的:两个差异中的每一个可以被分配给位于一对像素之间的两个相应的裂缝(一维单元)中的一个。然后，应将梯度分配给与两个裂纹相关的点(零维单元)。
按照 Canny (1986)的正确思路，需要求梯度绝对值的最大值。在彩色图像的情况下，有必要决定应该使用三种颜色梯度值或它们的分量的哪种组合来正确地表示颜色的变化。一种可能性是使用三原色差的三个绝对值的和，而不是灰度值的差，正如我们在灰度图像的情况下所做的那样。但是，与灰度图像的情况一样，总和必须有符号。如果我们取颜色强度的有符号差的和，不同原色的差可以具有不同的符号，并且相互补偿，使得它们的和变小。我们的研究表明，彩色图像中边缘位置附近的三个颜色梯度具有非常不同的方向，因此它们实际上可以在总和中相互补偿。
我们已经决定使用分配有亮度差符号的绝对差之和。这是通过 ColorDifSign 方法实现的，如下所示。
int ColorDifSign(byte[] Colp, byte[] Colh)
{
int Dif = 0;
for (int c = 0; c < 3; c++) Dif += Math.Abs(Colp[c] - Colh[c]);
int Sign;
if (MaxC(Colp[2], Colp[1], Colp[0])-MaxC(Colh[2], Colh[1], Colh[0])>0)
Sign = 1;
else Sign = -1;
return (Sign * Dif) / 3;
}
在这段代码中，MaxC( )是第三章中描述的方法，返回像素颜色的亮度。
为了研究方法LabelCellsSign的功能，我们在项目WFdetectEdges1中开发了一种工具，用于显示图像ExtremIm的一条线的亮度变化、色差与边缘检测阈值的关系以及检测到的边缘垂直裂缝。用户可以点击右边图片框中的一个点。然后，在下面的第三个图片框中会出现从单击点开始的被单击线的放大部分。在每个上部图片框中，会出现一条水平线，显示下部框中显示的线条部分。图 7-5 以项目WFdetectEdges1的形式为例。
图 7-5 中的下方曲线显示了图像SigmaIm中一条线的亮度。绝对值大于阈值的色差在上部曲线中显示为彩色垂直线。红线表示正差异，绿线表示负差异。红线和绿线都表示检测到的边缘垂直裂纹的位置。
图 7-5
项目的形式WFdetectEdges1
图 7-6 显示了带有图 7-5 图像片段边缘的细胞复合体。这是通过稍后描述的方法DrawComb的边缘表示。图 7-5 和图 7-6 中所示的边缘是通过LabelCellsSign在色带中寻找最大绝对色差的方法产生的，色带是两个同质区域之间的一个窄条。这条带子的颜色变化很快。如你所见，这种方法效果很好。
图 7-6
在图 7-5 的图像片段中检测到的边缘
然而，使用我们的方法ExtremLightUni,滑动窗口的尺寸相对较大，必须大于斜坡的宽度，这产生了在斜坡内看起来像细线的大部分窄边缘。因此，可以成功使用前面描述的更简单的方法LabelCells。
为了证明我们的边缘检测方法的成功，我们复制了来自维基百科的文章“Canny Edge Detector”中的图像，并将我们的方法应用于该文章中的彩色图像。图 7-7 显示了这个图像的一个片段(a)，我们的方法检测到的边缘(b)，以及 Canny 边缘检测器检测到的边缘。
图 7-7
(a)图像的片段，(b)通过我们的方法检测的边缘，以及(c)通过 Canny 边缘检测器检测的边缘
二值化梯度方法的进一步改进
通过方法LabelCellsSign计算出的边缘总是很细很细。然而，由于一些图像的不均匀结构，经常存在小块边缘，这会不必要地干扰边缘的进一步处理。因此，我们开发了方法CleanCombNew，它由三部分组成。第一部分将看起来像正方形的四个裂缝的每个结构转换成由两个裂缝组成的角；第二部分删除了与分叉点关联的单个裂纹，这些裂纹有一个关联端点；第三部分删除包含的单元数量小于预定阈值的边的分量。这个阈值是CleanCombNew的一个参数。
第一部分和第二部分的代码很简单，所以我们不做解释。然而，第三部分的代码相当复杂。在解释它之前，我应该提到我们已经开发了一个方法 DrawComb，它显示了一个放大的片段CombIm,将裂缝表示为短白线，将点表示为小的彩色矩形。对于可被检测为具有单个裂纹的入射点的终点，颜色为红色。带有两条裂纹的点的颜色为绿色，带有三条裂纹的点的颜色为黄色，带有四条裂纹的点的颜色为紫色。用户可以通过用鼠标点击右边的图片框来选择片段的位置。单击点指定片段左上角的位置。片段的大小是标准的；它是以这样的方式选择的，即放大的片段适合左边图片框的窗口。图 7-8 显示了一个例子。所选择的片段在图 7-8 的右侧显示为一个白色小方块。
图 7-8
由DrawComb表示的边缘片段
方法的第三部分CleanCombNew必须找到边的连通分量并计算它们的单元数。连通分量的概念在图论中是众所周知的。请注意，表示为单元复合体的边是顶点为点、边为裂缝的图。一条边的连通分量是它的子集，其中任意两个元素存在包含这些元素的关联路径，关联路径是任意两个相邻元素关联的序列a1，a2，a3，… a n 。图的一条边与它的端点关联，而这些端点又与这条边关联。
通过众所周知的广度优先搜索方法来寻找边的连通分量。这是一个遍历图形的算法。它从图中的某个任意顶点开始，在移动到下一级邻居之前，首先探索邻居顶点。
在方法CleanCombNew中，在一对for循环中搜索起始点，其中变量x和y仅取偶数。请注意，我们的细胞复合体中的任何一点都有两个偶数坐标。
搜索以前未使用的点。使用的点获得标记的最高位7(标记 128)。当找到一个带有标签1、3或4的未使用点时，方法ComponClean开始。它使用图像CombIm的网格副本，因为ComponClean标记了它使用的网格中的一些单元格，这些标记会干扰CombIm的进一步使用。方法ComponClean通过其子例程Trace跟踪组件的所有单元，并对其进行计数。只要计数的细胞数量小于阈值Size作为CleanCombNew的参数，细胞的坐标就作为值x + width * y保存在数组Index中。扫描完一个组件的所有单元后，将它的编号与Size进行比较:如果小于Size，则删除保存的单元。图 7-9 显示了使用CleanCombNew前后的边缘示例。
图 7-9
使用参数Size = 21 的CleanCombNew之前和之后的边缘(a)和(b)
Canny 边缘检测器
Canny (1986)著名的边缘检测器计算灰度值的梯度 G ( x ， y )=|Sobel。x |+|索贝尔。Y|使用 Sobel 滤镜(第六章)，计算边缘的方向(作为 45 的倍数)，并删除边缘上不具有梯度绝对值| G ( x ， y )|的局部最大值的所有候选。这是非最大值抑制 (NMS)的过程:如果像素( x ， y )的值为| G ( x ， y )|并且其不位于边缘方向的邻居的值更大，那么该值被设置为等于零。之后应用一个特殊的程序:它使用两个阈值T1T2。扫描图像，直到找到| G ( x ， y )|的值大于 T 2 的像素。这个像素位于边缘。该边缘在两个方向上被追踪，并且所有值大于TT1的像素被标记为属于该边缘。在这最后一步之后，算法就完成了。
图 7-10 比较了 Canny 算法和我们的二值化梯度的结果。正如你所看到的，二值化梯度提供了更多的细节，并且比 Canny 算法简单。
图 7-10
(Canny 算法和(b)二值化梯度的结果
彩色图像中的边缘
拉普拉斯的过零方法和 Canny 算法都要求将彩色图像转换成灰度图像。然而，在这种转换之后，彩色图像中的一些边缘消失了。如果两个相邻区域颜色不同但亮度相同，就会出现这种情况。
二值化梯度可以检测具有相同亮度但不同颜色的两个区域之间的真实颜色边缘。为此，需要适用于彩色图像的 sigma 滤波器和极限滤波器。第二章介绍了彩色图像的西格玛滤波器。用于彩色图像的极端过滤器具有一个特殊的特征:它必须找到具有最大差异的两种颜色，而不是在小滑动窗口中找到灰度值的最大值和最小值。两种颜色之间的差异是红色、绿色和蓝色通道的差异的平方和。输出图像中的中心像素获得与输入图像的中心像素的颜色具有较小差异的两种最不同颜色的颜色。
如前所述，我们的边缘检测器使用ColorDifSign方法计算两个相邻像素 P 1 (红色 1 ，绿色 1 ，蓝色 1 和 P 2 (红色 2 ，绿色 2 ，蓝色 2 )之间的色差，作为具有亮度差符号的颜色通道的绝对差的总和:
ColorDifSign =
sign of(light2–light 1)*(|红色2–红色1|+|绿色2–绿色1|+|蓝色2–蓝色 1 |)。
如果两个像素的色差为正且大于预定阈值，或者色差为负且小于(-阈值)，则边缘位于两个像素之间。我们的方法的优点可以在具有不同颜色和相等强度的两个区域之间的边缘上看到，如图 7-11 中红色花朵和具有几乎相同亮度的绿色叶子之间的边缘。注意箭头指示的位置。如你所见，图 7-11d 中缺失了许多边缘。
图 7-11
(a)原始图像；(b)彩色图像的边缘；(c)原始转换成灰度值；(d)灰度图像的边缘
结论
我们的边缘检测器(改进的二值化梯度)提供了比众所周知的 Canny 算法更好的结果。我们的方法也检测边缘的彩色图像更精确的方法比使用不同的灰度值。**
八、一种新的图像压缩方法
提出了一种新的基于边缘检测的图像压缩方法。这个想法是基于这样的假设，即图像中的重要信息位于靠近边缘的像素中。所有其他像素仅包含颜色的均匀分布，这可以通过用有限差分法数值求解偏微分方程来重建。这里我们用上一章描述的方法进行边缘检测。我们将边缘表示为 ACC，如第七章所述。这种表示具有优势，因为边缘的元素是由两个相邻像素的颜色差异定义的。因此，边缘的元素既不属于这两个像素，而是属于位于这两个像素之间的空间元素。这种空间元素在数字图像中不存在；但是，它存在于 ACC 中。这是被称为裂缝的一维单元，而像素是二维单元。裂缝的两端是称为点的零维单元。当把像素看成小方块时，裂缝就是这些方块的边；点是他们的角，是裂缝的终点。
这些边缘构成了一个细线网络，这些细线是一系列裂纹和点，其中每条线段要么是一条闭合曲线，要么是一个序列，其起点和终点要么是三条或四条线段相交的分支点(图 8-1 )。
图 8-1
边缘的例子
这种结构可以用所谓的小区列表来描述。单元格列表的概念基于将图像表示为 ACC(参考 Kovalevsky (1989，2008))。
使用单元复合体对边界进行编码
将数字图像视为 ACC 通过对子集的边界进行编码带来了优势:表示为一系列裂缝和点的边界变成了具有零面积的细曲线；区域的边界和它的补集的边界是相同的。数字曲线尤其是数字直线的定义和处理变得更加简单和清晰。从图像的经济编码和精确重建的角度来看，最重要的优点是能够使用极其简单和快速的算法来填充裂纹边界的内部，这种算法在将边界表示为像素组时无法应用。
第七章简要总结了对本演示重要的拓扑概念。更多细节和拓扑基础请参考 Kovalevsky (2008)。
我们经常需要欧几里得坐标来讨论模拟图像的数字化问题。为了保持在细胞复合体的框架中，我们建议将欧几里得坐标视为具有相对较大分母的有理数。这对应于任何计算机模型，因为计算机中的浮点变量是具有大分母的有理数。为了对细胞复合体中的边缘进行编码，我们使用大小为(2 宽度+ 1)(2 高度+1)的图像，其中宽度高度是原始图像的大小。我们称该图像为组合坐标中的图像或CombIm，因为该图像中被视为坐标的列数和行数指定了所有维度的细胞的位置(零、一和二)。与通常只包含像素的图像相比，在该图像中找到不同尺寸的细胞和跟踪边缘更加方便。组合坐标中像素的两个坐标都是奇数；一个裂纹有一个奇数和一个偶数坐标；一个点有两个偶数坐标。对比图 8-2 。
图 8-2
零维、一维和二维细胞的组合坐标
边由一系列裂纹和点组成。边缘不包含像素。边缘的裂缝位于不同颜色的两个像素之间。入射到边缘裂缝的两个像素的颜色具有大于预定阈值的差异。
图像的边缘集合由线组成，线是裂缝和点的连接序列。一条线可以关闭；在这种情况下，它不包含端点和分支点。闭合线的每一点都与两条裂纹相关。一条非闭合线的起点和终点与一条、三条或四条裂纹相关，但不与两条裂纹相关。一条线的每个内点恰好与两条裂纹相关(图 8-3 )。
图 8-3
线条(粗体线段)、分支点和短线的示例
要完整地描述一幅图像的边缘，只需要一列线条。为了重建图像，知道与每条线的裂缝相关的像素的颜色就足够了。然而，这些颜色沿着一条线的一边几乎是不变的。如果位于一条线一侧的像素组中的颜色有很大变化，那么也会有垂直于该线的边缘裂纹；但是，没有这样的裂缝。因为线条一侧的颜色变化缓慢，所以将线条每一侧的颜色分别保存在线条开头和结尾的一个像素中就足够了。线的一侧的所有其他颜色可以通过在线的开始和结束处的这些颜色之间的插值来计算。
为了执行这种插值，需要知道线的所有元素的准确位置，这可以通过起点的坐标和线的所有裂纹的方向序列来指定。定向线的裂缝(即指定哪个点是起点)可以有四个方向之一:方向 0 向右，方向 1 向下，方向 2 向左，方向 3 向上。这样，可以非常节省地对图像的边缘线和相邻的颜色进行编码。图像边界的颜色对于重建也是必要的。只知道矩形图像四边的颜色就足够了。边界处的所有其他颜色都可以通过插值来计算。在该插值期间，也可以考虑与边界像素序列交叉的边缘裂缝处的颜色。
在对线侧的颜色和边界处的颜色进行插值之后，在线之间的区域中仍然存在大量像素，这些区域在所有插值之后仍然是空的。这些像素中的颜色可以通过沿着图像的行和列的线性插值来近似计算。这种插值必须从边界或线处的颜色开始，并且可以沿着一行继续，直到已经出现的颜色出现。可以沿着列进行类似的插值。然后每个像素获得两种颜色:一种来自沿行的插值，另一种来自沿列的插值。这两种颜色是平均的。经过这些插值后，颜色的分布不是很均匀。通过拉普拉斯偏微分方程的数字解进行平滑，可以使其更加均匀。经过这种平滑处理后，重建图像看起来与原始图像非常相似。
项目描述WFcompressPal
运行项目的形式如图 8-4 所示。
图 8-4
运行项目的形式WFcompressPal
当用户点击打开图像时，他或她可以选择一个图像，该图像将作为原始图像OrigBmp打开并显示在左侧图片框中。在项目的这个区块中定义了八个工作图像。必须将OrigBmp转换成工作图像OrigIm。软件WindowsForms建议通过GetPixels的方法访问位图的像素，这个方法相当慢。我们开发了一个更快的方法BitmapToImage，它使用了类Bitmap的方法LockBits来处理彩色图像。然而，如果原始图像是每像素 8 比特的图像，则意味着打开的图像是索引图像。方法LockBits不应该用于每像素 8 位的索引位图，因为这些位图包含一个调色板，像素值总是由调色板值定义。因此，LockBits方法对于 8 位索引位图来说太慢了。对于这样的位图，我们使用标准方法GetPixel。除了EdgeIm和MaskIm之外的所有工作图像将被定义为 24 位图像。图像EdgeIm和MaskIm始终是 8 位图像。
现在，用户可以决定是否抑制脉冲噪声。当原始图像不包含脉冲噪声时，该过程也是有用的，因为去除亮度不同于周围像素亮度的小斑点使得图像更加均匀，并且人类感知的图像质量保持不变。使图像更加均匀会提高压缩率。
如果用户决定使用脉冲噪声抑制，他或她可以改变要去除的斑点的最大尺寸，而带有删除暗和删除亮标签的工具的值。只需将“删除暗”和“删除亮”的值设置为零，就可以跳过此过程。用户可以点击脉冲噪声。第二章描述了抑制脉冲噪声的方法。
现在有必要对图像进行分割。当用户点击 Segment 时，脉冲噪声被抑制的图像将被 sigma 滤波器抑制高斯噪声(第二章)和极值滤波器增加相邻像素之间的色差(第六章)处理。极值过滤器的处理结果保存在图像ExtremIm中。然后，在彩色图像的情况下，将计算代表被处理图像的颜色的 256 种颜色的调色板Palet,并且创建具有该调色板的索引图像Pal。图像SegmentIm是通过将索引图像Pal转换成真彩色图像而获得的真彩色图像。它显示在右边的图片框中。
调色板的计算通过MakePalette方法进行。这个方法创建了一个包含 256 种颜色的数组Palette，并提供了颜色编号，称为indices。Palette的颜色近似代表原始图像中出现的所有颜色。
方法MakePalette首先在空间中指定一个三维立方体，其坐标轴对应于红色、绿色和蓝色通道。立方体有大小
[最小[c]，最大[c]]，c = 0，1，2
其中Min[c]是图像中出现的颜色通道c的最小值，Max[c]是最大值。因此，这个三维颜色立方体包含了图像中出现的所有颜色。颜色的立方体细分为 10×10×10 = 1000 个小立方体。颜色立方体中小立方体的位置由三个坐标iComp[c], c = 0, 1, 2指定。每个数字iComp[c]可以取 0 到 9 的值。根据预定值MaxComp[3]={9, 9, 9}，通过将像素颜色的三个通道值(红色、绿色和蓝色)除以值Div[c]，为图像的每个像素计算三个数字iComp[c], c = 0, 1, 2。三个数字iComp[c]指定了值index，该值可以在 0 到 999 之间。index指定像素颜色所属的小立方体的位置。index的值为:
并且在方法MakePalette中计算取决于MaxComp的值Weight[c]。方法MakePalette为每个小立方体计算颜色位于该小立方体中的原始图像的像素数。
非空小立方体的数量指定了调色板中颜色的数量。我们希望这个数字在 200 到 255 之间。如果低于 200，则MaxComp[1](绿色通道)的值增加。因此，小立方体的数量增加，并且活动调色板元素的数量也随之增加。如果提到的数字超过最大值 255(我们为索引 0 保留调色板的一个值为(0，0，0))，则MaxComp[1]减少。因此，活动小立方体的数量获得 200 和 255 之间的值。在小图像的情况下，这种控制过程不起作用。当MaxComp[1]大于 20 时中断，在这种情况下，活动调色板元素的数量可能超出所需的间隔。
方法MakePalette为每个非空的小立方体计算属于该立方体的所有颜色的平均值。还计算颜色的标准偏差(这部分代码在下面的文本版本中被省略)。感兴趣的标准偏差只是给出计算的调色板如何精确地表示原始图像的颜色的概念。
平均值用作调色板的值。对非空的小立方体进行计数，计数值是调色板中使用的颜色(或颜色号)的新索引数。调色板保存为数组。MakePalette创建调色板图像Pal。颜色索引从图像Pal复制到图像Comb的像素中。这里是MakePalette的源代码。
public int MakePalette(CImage Img, int[] Palet, Form1 fm1)
// Produces a palette of 256 elements optimally representing colors of the image "Img"
{
bool deb = false;
int MaxN = 1000, ResIndex = 10000;
int[] Sum = new int[3 * MaxN];
int[] nPix = new int[MaxN]; // "nPix[Index]" is number of pixels with index
double[] Mean = new double[3 * MaxN];
int[] Sum0 = new int[3 * ResIndex];
int[] nPix0 = new int[ResIndex];
int c, i, jump, nIndex, n, MaxIndex;
int[] iChan = new int[3], Div = new int[3], Weight = new int[3];
int[] NumbInterv = { 11, 12, 9 };
// Computing minimum and maximum of the color channels:
int[] Min = { 256, 256, 256 }, Max = { 0, 0, 0 };
fm1.progressBar1.Visible = true;
fm1.progressBar1.Step = 1;
if (Img.width * Img.height > 300) jump = Img.width * Img.height / 20;
else jump = 3;
for (i = 0; i < Img.width * Img.height; i++) //=======================
{
if (i % jump == jump - 1) fm1.progressBar1.PerformStep();
for (c = 0; c < 3; c++)
{
if (Img.Grid[3 * i + c] < Min[c]) Min[c] = Img.Grid[3 * i + c];
if (Img.Grid[3 * i + c] > Max[c]) Max[c] = Img.Grid[3 * i + c];
}
} //=============================== end for (i... =============
int nIndexMin = 228, nIndexMax = 255;
nIndex = 0;
// Changing NumbInterv[] to get nIndex between nIndexMin and nIndexMax:
do //===================================================
{
Weight[0] = 1;
Weight[1] = NumbInterv[0] + 1;
Weight[2] = NumbInterv[0] + Weight[1] * NumbInterv[1] + 1;
for (i = 0; i < ResIndex; i++) Sum0[i] = 0;
for (i = 0; i < 3 * MaxN; i++) //==================
{
Sum[i] = 0;
Mean[i] = 0.0;
nPix[i / 3] = 0;
} //================== end for (i ... ===========
for (c = 0; c < 3; c++)
{
Div[c] = (int)(0.5 + (double)(Max[c] - Min[c]) / (double)NumbInterv[c]);
if (Div[c] == 0) Div[c] = 1;
if ((Max[c] - Min[c]) / Div[c] > NumbInterv[c]) NumbInterv[c]++;
}
MaxIndex = Weight[0] * NumbInterv[0] + Weight[1] * NumbInterv[1] +
Weight[2] * NumbInterv[2];
int maxIndex = 0;
if (MaxIndex >= ResIndex)
{ MessageBox.Show("MakePalette, Overflow: MaxIndex=" +
MaxIndex + " > ResIndex=" + ResIndex + "; return -1.");
return -1;
}
for (i = 0; i < ResIndex; i++) nPix0[i] = 0;
int Index = 0;
for (i = 0; i < Img.width * Img.height; i++) //=======================
{
Index = 0;
for (c = 0; c < 3; c++)
{
iChan[c] = (Img.Grid[3 * i + c] - Min[c]) / Div[c];
Index += Weight[c] * iChan[c];
}
if (Index > maxIndex) maxIndex = Index;
if (Index > ResIndex - 1)
{
MessageBox.Show("MP, Overflow: Index=" + Index + " too great; return -1.");
return -1;
}
for (c = 0; c < 3; c++) Sum0[3 * Index + c] += Img.Grid[3 * i + c];
nPix0[Index]++;
} //================ end for (i = 0; ... ============================
nIndex = 0;
for (i = 0; i <= MaxIndex; i++) if (nPix0[i] > 0) nIndex++;
int minInd = 0, maxInd = 0;
if (nIndex < nIndexMin)
{
if (NumbInterv[0] <= NumbInterv[1] && NumbInterv[0] <= NumbInterv[2])
minInd = 0;
else
if (NumbInterv[1] <= NumbInterv[0] && NumbInterv[1] <= NumbInterv[2])
minInd = 1;
else
if (NumbInterv[2] <= NumbInterv[0] && NumbInterv[2] <= NumbInterv[1])
minInd = 2;
NumbInterv[minInd]++;
}
if (nIndex > nIndexMax)
{
if (NumbInterv[0] >= NumbInterv[1] && NumbInterv[0] >= NumbInterv[2])
maxInd = 0;
else
if (NumbInterv[1] >= NumbInterv[0] && NumbInterv[1] >= NumbInterv[2])
maxInd = 1;
else
if (NumbInterv[2] >= NumbInterv[0] && NumbInterv[2] >= NumbInterv[1])
maxInd = 2;
NumbInterv[maxInd]--;
}
if (nIndex >= nIndexMin && nIndex <= nIndexMax || NumbInterv[1] > 20)
{
if (deb)
MessageBox.Show("MakePalette: nIndex=" + nIndex + " is OK. break.");
break;
}
} while (nIndex > nIndexMax || nIndex < nIndexMin); //===================
int[] NewIndex = new int[MaxIndex];
if (MaxIndex > 300) jump = MaxIndex / 20;
else jump = 3;
for (i = n = 0; i < MaxIndex; i++) //===================================
{
if (i % jump == jump - 1) fm1.progressBar1.PerformStep();
if (nPix0[i] > 0) //---------------------------------------------------
{
n++;
if (n > MaxN - 1)
{
MessageBox.Show("MP: Overflow in Sum; n=" + n + "< MaxN=" + MaxN);.
return -1;
}
NewIndex[i] = n;
nPix[n] = nPix0[i];
for (c = 0; c < 3; c++)
{
Sum[3 * n + c] = Sum0[3 * i + c];
if (nPix[n] > 0) Mean[3 * n + c] = (double)(Sum[3 * n + c]) / (double)(nPix[n]);
else Mean[3 * n + c] = 0;
}
} //-------------------------- end if (nPix0... -----------------------
} //====================== end for (i... ===========================
int MaxNewIndex = n;
if (Img.width * Img.height > 300) jump = Img.width * Img.height / 20;
else jump = 3;
// Putting NewIndex into "this.Grid":
for (i = 0; i < Img.width * Img.height; i++) //==========================
{
if (i % jump == jump - 1) fm1.progressBar1.PerformStep();
int Index = 0;
for (c = 0; c < 3; c++)
{
iChan[c] = (Img.Grid[3 * i + c] - Min[c]) / Div[c];
Index += Weight[c] * iChan[c];
}
if (Index >= MaxIndex) Index = MaxIndex - 1;
Grid[i] = (byte)NewIndex[Index];
if (Grid[i] == 0) Grid[i] = (byte)MaxNewIndex;
} //=============== end for (i=0; ... ===============================
// Calculating "Palet" and "this.Palette":
byte R = 0, G = 0, B = 0;
jump = MaxNewIndex / 20;
for (n = 0; n <= MaxNewIndex; n++)
{
if (n % jump == jump -1) fm1.progressBar1.PerformStep();
if (n > 0)
{
if (Mean[3 * n + 2] < 255.0) R = (byte)Mean[3 * n + 2];
if (Mean[3 * n + 1] < 255.0) G = (byte)Mean[3 * n + 1];
if (Mean[3 * n + 0] < 255.0) B = (byte)Mean[3 * n + 0];
Palet[n] = RGB(R, G, B);
}
else
{
Palet[n] = 0;
}
}
return 1;
} //******************* end MakePalette ********************************
现在可以执行边缘检测。用户可以选择用于边缘检测的阈值，该阈值指定两个相邻像素的颜色之间的最小差值，在该值处，这些像素之间的裂纹将被标记为属于边缘。用户可以在标签阈值下的工具numericUpDown的小窗口中输入一个新值，或者点击该设置右侧的两个小箭头之一来增加或减少该值。用户应该考虑，在较高的阈值下，将产生较少量的边缘元素。这导致更高的压缩率，但也降低了恢复图像的质量。
设计用于检测边缘的项目部分使用组合坐标(参见第七章)和更大的尺寸(2 宽度+ 1)(2 *高度+1)定义了每像素一个字节的图像CombIm。这是必要的，因为在标准尺寸的图像中表示裂缝和点是困难的。图像CombIm将包含代表检测到的边缘的细胞复合体。在第七章中可以找到对细胞复合体的描述。方法LabelCellsSign测试图像的每对相邻像素ExtremIm。如果像素的色差大于阈值，则这些像素之间的裂缝用1标记。这意味着对应于该裂纹的图像元素CombIm获得值1。同时，裂纹端点的标记将增加，从而获得与该点相关的裂纹数量相等的标记。点的标签示例如图 8-5 所示。
图 8-5
维和一维单元格的标签
如前所述，调色板索引从由MakePalette产生的图像Pal复制到图像Comb的像素中。在灰度图像的情况下，灰度值是从Image3复制的，通过ExtremVar方法进行处理。
现在可以开始零件Encode了。当用户点击编码时，使用图像CombIm作为参数来启动方法SearchLin。该方法查看CombIm的所有点，并为相应字节的最低有效位中标记为 1、3 或 4 的每个点调用方法ComponLin。这些是边缘线的端点或分支点。由 SearchLin 调用的方法ComponLin跟踪和编码连接组件的所有行。追踪过程中遇到的所有点的标签将被删除。
在SearchLin处理完标签为 1、3 或 4 的所有点后，它寻找标签为 2 的未删除点。这些点属于没有端点和分支点的组件。这样的组件是一条闭合曲线。
这里是SearchLin的代码。
public int SearchLin(ref CImage Comb, Form1 fm1)
{ int Lab, rv, x, y;
for (x=0; x<MaxByte; x++) Byte[x]=0;
if (Step[0].Y < 0) x = Step[0].Y;
fm1.progressBar1.Value = 0;
int y1 = Comb.nLoop * CNY / Comb.denomProg;
for (y=0; y<CNY; y+=2)
{
if ((y % y1) == 0) fm1.progressBar1.PerformStep();
for (x=0; x<CNX; x+=2)
{ Lab=Comb.Grid[x+CNX*y] & 3;
if (Lab==1 || Lab==3)
{ rv=ComponLin(Comb.Grid, x, y);
if (rv<0)
{ MessageBox.Show("SearchLin, Alarm! ComponLin returned " + rv);
return -1;
}
}
}
}
// Starting the search for loops:
for (y=0; y<CNY; y+=2)
for (x=0; x<CNX; x+=2)
{ Lab=Comb.Grid[x+CNX*y] & 3;
if (Lab==2)
{ rv=ComponLin(Comb.Grid, x, y);
if (rv<0)
{ MessageBox.Show("SearchLin, Alarm! ComponLin returned " + rv);
return -1;
}
}
}
nByte++;
return nByte;
} //********************* end SearchLin ********************************
方法ComponLin跟踪并编码每个完整组件；即边的连接部分。这种方法被设计用于追踪连通图，因为边的连通分量是具有点的顶点和具有线的边的连通图。众所周知，遍历或搜索树或图数据结构可以通过广度优先算法完成，该算法使用称为队列的先进先出数据结构。
ComponLin将组件的起点放入队列。当从队列中获得一个点时，它测试与该点相关的所有四个裂纹。如果这些裂纹中的一条终止于一个端点或一个分支点，那么它被认为是一条由单个裂纹组成的短线。它将立即被编码到数据结构CListLines中，其中ComponLin是一个方法。
如果裂纹在标记为 2 的点处结束，则方法TraceLin开始。该方法跟踪一条长线，保存该线两侧的颜色索引或灰度值，并保存该线所有裂缝的方向。当TraceLin返回且处理后的图像是彩色图像时，则启动方法FraqInds，该方法计算该行的两半中最频繁的索引，并将它们分配给该行的存储代码。在灰度图像的情况下，启动AverageGrays,计算线的两半的平均灰度值。然后将平均灰度值分配给该行的代码。
下面是方法TraceLin和ComponLin的代码:
public int TraceLin(byte[] CGrid, int X, int Y, ref iVect2 Pterm, ref int dir)
/* This method traces a line in the image "Comb" with combinatorial coordinates, where the cracks and points of the edges are labeled: bits 0, 1 and 2 of a point contain the label 1 to 4 of the point. The label indicates the number of incident edge cracks. Labeled bit 6 indicates that the point already has been put into the queue; labeled bit 7 indicates that the point should not be used any more. The crack has only one label 1 in bit 0\. This method traces the edge from one end or branch point to another while changing the parameter "dir". ----------*/
{ bool atSt_P=false, BP=false, END=false;
int rv = 0;
iVect2 Crack, P=new iVect2(0,0), PixelP, PixelN, StartPO;
int iShift=-1, iCrack=0, Lab;
P.X=X; P.Y=Y;
StartPO=P;
if (nLine2==0) nByte=0;
else nByte=Line2[nLine2-1].EndByte+1;
int[] Shift={0,2,4,6};
while(true) //================================================
{ Crack=P+Step[dir];
P=Crack+Step[dir];
Lab=CGrid[P.X+CNX*P.Y]&3;
switch(Lab)
{ case 1: END=true; BP=false; rv=1; break;
case 2: BP=END=false; break;
case 3: BP=true; END=false; rv=3; break;
}
PixelP=Crack+Norm[dir];
PixelN=Crack-Norm[dir];
IndPos[iCrack]=CGrid[PixelP.X+CNX*PixelP.Y];
IndNeg[iCrack]=CGrid[PixelN.X+CNX*PixelN.Y];
if (Lab==2) CGrid[P.X+CNX*P.Y]=0;
iShift++;
iCrack++;
if (iShift==4)
{ iShift=0;
if (nByte<MaxByte-1) nByte++;
else
{ return -1;
}
}
Byte[nByte] |= (byte)(dir << Shift[iShift]);
if (P.X == StartPO.X && P.Y == StartPO.Y) atSt_P = true;
else atSt_P = false;
if (atSt_P)
{ Pterm=P;
rv=2;
break;
}
if (!BP && !END) //------------------------------------------
{ Crack=P+Step[(dir+1)%4];
if (CGrid[Crack.X+CNX*Crack.Y]==1)
{ dir=(dir+1)%4;
}
else
{ Crack=P+Step[(dir+3)%4];
if (CGrid[Crack.X+CNX*Crack.Y]==1)
{ dir=(dir+3)%4;
}
}
}
else
{ Pterm=P;
break;
} //----------------- end if (!BP && !END) --------------------
} //=========== end while =====================================
Line2[nLine2].EndByte=nByte;
Line2[nLine2].nCrack=(ushort)iCrack;
return rv;
} //*************** end TraceLin ****************************************
public int ComponLin(byte[] CGrid, int X, int Y)
/* Encodes in "CListLines" the lines of the edge component with the point (X, Y) being a branch or an end point. Puts the starting point 'Pinp' into the queue and starts the 'while' loop. It tests each labeled crack incident with the point 'P' fetched from the queue. If the next point of the crack is a branch or an end point, then a short line is saved. Otherwise the method "TraceLin" is called. "TraceLin" traces a long line, saves the color indices at the sides of the line and ends at the point 'Pterm' with the direction 'DirT'. Then the method "FreqInds" assigns the most frequent from the saved color indices to the line. If the point 'Pterm' is a branch point then it is put to the queue. "ComponLin" returns when the queue is empty. ---------------*/
{ int dir, dirT;
int LabNext, rv;
iVect2 Crack, P, Pinp, PixelN, PixelP, Pnext, Pterm=new iVect2(0, 0);
Pinp=new iVect2(X,Y); // comb. coord.
pQ.Put(Pinp);
while(pQ.Empty()==false) //======================================
{ P=pQ.Get();
for (dir=0; dir<4; dir++) //======================================
{ Crack=P+Step[dir];
if (Crack.X<0 || Crack.X>CNX-1 || Crack.Y<0 || Crack.Y>CNY-1 ) continue;
if (CGrid[Crack.X+CNX*Crack.Y]==1) //----------------------------
{ PixelP=Crack+Norm[dir]; PixelN=Crack-Norm[dir];
Pnext=Crack+Step[dir];
LabNext=CGrid[Pnext.X+CNX*Pnext.Y] & 3; //Ind0
if (LabNext==1 || LabNext==3)
{
Line1[nLine1].x = (ushort)(P.X / 2); Line1[nLine1].y = (ushort)(P.Y / 2);
Line1[nLine1].Ind0=CGrid[PixelN.X+CNX*PixelN.Y];
Line1[nLine1].Ind1=CGrid[PixelP.X+CNX*PixelP.Y];
if (nLine1>MaxLine1-1)
{
MessageBox.Show("ComponLin: Overflow in Line1; return -1");
return -1;
}
}
if (LabNext==3) pQ.Put(Pnext);
if (LabNext==2) //--------------------------------------------
{ Line2[nLine2].x=(ushort)(P.X/2);
Line2[nLine2].y=(ushort)(P.Y/2); //transf. to standard coordinates
dirT=dir;
rv=TraceLin(CGrid, P.X, P.Y, ref Pterm, ref dirT);
if (nBits3==24) FreqInds(nLine2);
else AverageGray(nLine2);
if (rv<0)
{
return -1;
}
if (nLine2>MaxLine2-1)
{ return -1;
}
else nLine2++;
if ((CGrid[Pterm.X+CNX*Pterm.Y] & 64)==0) // '64'= label for visited;
{
if (rv==3) //------------------------------------------
{ CGrid[Pterm.X+CNX*Pterm.Y] |=64;
pQ.Put(Pterm);
}
} //-------- end if ((CGrid[Pterm.X... --------------------
} // ---------- end if (LabNest==2) ---------------------------
if ((CGrid[P.X+CNX*P.Y]&3)==1) break;
} //--------------- end if (CGrid[Crack.X ...==1) ----------------
} //========== end for (dir ... ===================================
CGrid[P.X+CNX*P.Y]=0;
} //============ end while ======================================
return 1;
} //*************** end ComponLin *************************************
当SearchLin返回并且所有行的代码都保存在CListLines类的对象List中时，创建CListCode类的对象LiCod。这个类是CListLines类的简短版本。它只包含重建图像所需的信息。方法Transform将必要的信息从List复制到LiCod，并将数据转换成包含在数组ByteNew中的字节序列。为了避免使用方法Serialize在磁盘上保存代码，这个转换是必要的，因为Serialize产生的磁盘文件有时比系统方法Write产生的文件长十倍。方法Serialize因此会破坏压缩，因此我们不使用它。
下面是Transform的代码。
public int Transform(int nx, int ny, int nbits, int[] Palet, CImage Comb, CListLines L, Form1 fm1)
// Transforms the provisional list "L" to an object of the class "CListCode".
{ int i, ib, il, nCode=0;
width=nx; height=ny; nBits=nbits;
nCode+=3*4;
nLine1=L.nLine1; nLine2=L.nLine2; nByte=L.nByte;
nCode+=3*4;
nCodeAfterLine2 = nCode;
Palette=new int[256];
for ( i=0; i<256; i++) Palette[i]=Palet[i]; nCode+=256*4;
Corner[0]=Comb.Grid[1+(2*width+1)*1]; // this is a gray value or a palette index
Corner[1]=Comb.Grid[2*width-1+(2*width+1)*1];
Corner[2]=Comb.Grid[2*width-1+(2*width+1)*(2*height-1)];
Corner[3]=Comb.Grid[1+(2*width+1)*(2*height-1)];
nCode+=4;
int il1=Comb.nLoop*nLine1/Comb.denomProg;
for (il = 0; il < nLine1; il++)
{
Line1[il] = L.Line1[il];
if (((il + 1) % il1) == 0) fm1.progressBar1.PerformStep();
}
nCode += nLine1*6; // "6" is the sizeof(CCrack);
int il2 = Comb.nLoop * nLine2 / Comb.denomProg;
for (il = 0; il < nLine2; il++)
{ Line2[il]=L.Line2[il];
if (((il + 1) % il2) == 0) fm1.progressBar1.PerformStep();
}
nCode += nLine2 * 14; // sizeof(CLine);
int il3 = Comb.nLoop * nByte / Comb.denomProg;
for (ib = 0; ib < nByte; ib++)
{ Byte[ib]=L.Byte[ib];
if (((ib + 1) % il3) == 0) fm1.progressBar1.PerformStep();
}
nCode += nByte;
// The following code is necessary to avoid "Serialize":
ByteNew = new byte[nCode+4];
for (int ik = 0; ik < nCode + 4; ik++) ByteNew[ik] = 0;
int j = 0;
ByteNew[j] = (byte)(nCode & 255); j++;
ByteNew[j] = (byte)((nCode >> 8) & 255); j++;
ByteNew[j] = (byte)((nCode >> 16) & 255); j++;
ByteNew[j] = (byte)((nCode >> 24) & 255); j++;
ByteNew[j] = (byte)(nx & 255); j++;
ByteNew[j] = (byte)((nx >> 8) & 255); j++;
ByteNew[j] = (byte)((nx >> 16) & 255); j++;
ByteNew[j] = (byte)((nx >> 24) & 255); j++;
ByteNew[j] = (byte)(ny & 255); j++;
ByteNew[j] = (byte)((ny >> 8) & 255); j++;
ByteNew[j] = (byte)((ny >> 16) & 255); j++;
ByteNew[j] = (byte)((ny >> 24) & 255); j++;
ByteNew[j] = (byte)(nbits & 255); j++;
ByteNew[j] = (byte)((nbits >> 8) & 255); j++;
ByteNew[j] = (byte)((nbits >> 16) & 255); j++;
ByteNew[j] = (byte)((nbits >> 24) & 255); j++;
ByteNew[j] = (byte)(nLine1 & 255); j++;
ByteNew[j] = (byte)((nLine1 >> 8) & 255); j++;
ByteNew[j] = (byte)((nLine1 >> 16) & 255); j++;
ByteNew[j] = (byte)((nLine1 >> 24) & 255); j++;
ByteNew[j] = (byte)(nLine2 & 255); j++;
ByteNew[j] = (byte)((nLine2 >> 8) & 255); j++;
ByteNew[j] = (byte)((nLine2 >> 16) & 255); j++;
ByteNew[j] = (byte)((nLine2 >> 24) & 255); j++;
ByteNew[j] = (byte)(nByte & 255); j++;
ByteNew[j] = (byte)((nByte >> 8) & 255); j++;
ByteNew[j] = (byte)((nByte >> 16) & 255); j++;
ByteNew[j] = (byte)((nByte >> 24) & 255); j++;
for (int ii = 0; ii < 256; ii++)
{
ByteNew[j] = (byte)(Palet[ii] & 255); j++;
ByteNew[j] = (byte)((Palet[ii] >> 8) & 255); j++;
ByteNew[j] = (byte)((Palet[ii] >> 16) & 255); j++;
ByteNew[j] = (byte)((Palet[ii] >> 24) & 255); j++;
}
for (int i1 = 0; i1 < 4; i1++) ByteNew[j+i1] = Corner[i1];
j+=4;
for (int i2 = 0; i2 < nLine1; i2++)
{
ByteNew[j] = (byte)(L.Line1[i2].x & 255); j++;
ByteNew[j] = (byte)((L.Line1[i2].x >> 8) & 255); j++;
ByteNew[j] = (byte)(L.Line1[i2].y & 255); j++;
ByteNew[j] = (byte)((L.Line1[i2].y >> 8) & 255); j++;
ByteNew[j] = L.Line1[i2].Ind0; j++;
ByteNew[j] = L.Line1[i2].Ind1; j++;
}
for (int i3 = 0; i3 < nLine2; i3++)
{
ByteNew[j] = (byte)(L.Line2[i3].EndByte & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].EndByte >> 8) & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].EndByte >> 16) & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].EndByte >> 248) & 255); j++;
ByteNew[j] = (byte)(L.Line2[i3].x & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].x >> 8) & 255); j++;
ByteNew[j] = (byte)(L.Line2[i3].y & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].y >> 8) & 255); j++;
ByteNew[j] = (byte)(L.Line2[i3].nCrack & 255); j++;
ByteNew[j] = (byte)((L.Line2[i3].nCrack >> 8) & 255); j++;
ByteNew[j] = L.Line2[i3].Ind0; j++;
ByteNew[j] = L.Line2[i3].Ind1; j++;
ByteNew[j] = L.Line2[i3].Ind2; j++;
ByteNew[j] = L.Line2[i3].Ind3; j++;
}
for (int i4 = 0; i4 < nByte; i4++) ByteNew[j + i4] = L.Byte[i4];
j += nByte;
return nCode;
} //************************** end Transform ******************************
当方法Transform返回时，显示消息“图像已编码”,并显示代码长度和压缩率的值。
项目WFcompressPal还包含一个通过单击 Restore 调用的部分，这是从代码重建图像所必需的。这是必需的，因为用户必须看到重建图像的质量是否足够。否则，用户可以指定较低的阈值来获得更精细的边缘。这自然会导致压缩率变小。项目的这一部分与下一节描述的项目WFrestoreLin的相应部分相同。
单击 Save code 调用方法WriteCode，该方法将数组ByteNew和代码写入扩展名为*.dat的文件中。
如果用户没有按正确的顺序点击按钮，他或她会得到一个警告和正确顺序的指示。
项目WFrestoreLin
这个项目用于从保存的代码中重建图像，并将图像存储在磁盘上。该项目从读取一个选择的扩展名为*.dat的文件开始。系统方法Read首先读取代码的长度，这使得分配字节数组ByteNew成为可能，然后将代码读入数组ByteNew。类CListCode的构造器将ByteNew的字节转换成类CListCode的对象LiCod。该对象包含要重建的图像的参数；数组Line1、Line2和Byte的长度；调色板；图像四个角的调色板索引；这些数组。
当用户点击恢复时，用LiCod : RestoreIm中包含的参数定义两幅图像进行重建，MaskIm作为辅助图像。两幅图像的所有像素都被设置为零。然后方法Restore重建显示在pictureBox1中的图像。在项目WFcompressPal中也使用这种方法来重建编码图像，以便用户可以评估重建图像的质量。
方法Restore首先将图像角的索引或颜色放入相应的像素。如果要重建的图像是彩色图像，则Restore通过调色板将读取的索引转换成颜色。
对应于获得颜色或灰度级的像素Image的图像像素Mask总是被设置为LabMask=250。(选择这个高值是为了能够在调试时显示Mask的内容。)
然后Restore读取具有短线代码的数组Line1并将颜色或灰度值放入短线的单个裂缝两侧的像素中。如果一条线只有一条裂纹，那么这条线就是短线。其他的都是长线。Restore然后读取带有长行代码的数组Line2，并将颜色或灰度值放入每一长行开头和结尾的像素中。Restore下一步沿着长线的边执行颜色或灰度值的插值来完成它的工作。
这里是Restore的代码。
public int Restore(ref CImage Image, ref CImage Mask, Form1 fm1)
{ int dir, nbyte, x, y;
byte LabMask=250;
if (nBits==24) nbyte=3;
else nbyte=1;
fm1.progressBar1.Value = 0;
fm1.progressBar1.Visible = true;
int denomProg = fm1.progressBar1.Maximum / fm1.progressBar1.Step;
int Sum = nLine1 + nLine2;
int i1 = Sum / denomProg;
for (int i = 0; i < width * height * (nBits / 8); i++) Image.Grid[i] = 0;
for (int i=0; i<width*height; i++) Mask.Grid[i]=0;
if (nBits==24)
{ for (int c=0; c<nbyte; c++)
{ Image.Grid[c]=(byte)((Palette[Corner[0]]>>8*(2-c)) & 0XFF); // left below
Image.Grid[nbyte*(width-1)+c]=
(byte)((Palette[Corner[1]]>>8*(2-c)) & 0XFF); // right below
Image.Grid[nbyte*width*height-nbyte+c]=
(byte)((Palette[Corner[2]]>>8*(2-c)) & 0XFF); // right on top
Image.Grid[nbyte*width*(height-1)+c]=
(byte)((Palette[Corner[3]]>>8*(2-c)) & 0XFF); // left on top
}
}
else
{ Image.Grid[0]=Corner[0];
Image.Grid[width-1]=Corner[1];
Image.Grid[width*height-1]=Corner[2];
Image.Grid[0+width*(height-1)]=Corner[3];
}
Mask.Grid[0]=Mask.Grid[width-1]=Mask.Grid[width*height-1]=
Mask.Grid[width*(height-1)]=LabMask;
// Short lines:
fm1.progressBar1.Value = 0;
for (int il = 0; il < nLine1; il++) //===================================
{
if ((il % i1) == 0) fm1.progressBar1.PerformStep();
dir=((Line1[il].x>>14) & 2) | (Line1[il].y>>15);
x=Line1[il].x & 0X7FFF; y=Line1[il].y & 0X7FFF;
if (nBits==24)
{ switch(dir)
{ case 0:
if (y > 0)
{
for (int c = 0; c < nbyte; c++)
{
int Index = Line1[il].Ind1;
byte col = (byte)(Palette[Index] >> 8 * c);
Image.Grid[nbyte * (x + width * y) + 2 - c] = col;
Image.Grid[nbyte * (x + width * (y - 1)) + 2 - c] =
(byte)((Palette[Line1[il].Ind0] >> 8 * c) & 0XFF);
}
Mask.Grid[x + width * y] = Mask.Grid[x + width * (y - 1)] = LabMask;
}
break;
case 1:
for (int c = 0; c < nbyte; c++)
{ Image.Grid[nbyte*(x+width*y)+2-c]=(byte)((Palette[Line1[il].Ind0]>>8*c) & 0XFF);
Image.Grid[nbyte*(x-1+width*y)+2-c]=(byte)((Palette[Line1[il].Ind1]>>8*c) & 0XFF);
}
Mask.Grid[x+width*y]=Mask.Grid[x-1+width*y]=LabMask;
break;
case 2:
for (int c = 0; c < nbyte; c++)
{ Image.Grid[nbyte*(x-1+width*y)+2-c]=(byte)((Palette[Line1[il].Ind0]>>8*c) & 0XFF);
Image.Grid[nbyte*(x-1+width*(y-1))+2-c]=(byte)((Palette[Line1[il].Ind1]>>8*c) & 0XFF);
}
Mask.Grid[x-1+width*y]=Mask.Grid[x-1+width*(y-1)]=LabMask;
break;
case 3:
for (int c = 0; c < nbyte; c++)
{ Image.Grid[nbyte*(x+width*(y-1))+2-c]=(byte)((Palette[Line1[il].Ind1]>>8*c) & 0XFF);
Image.Grid[nbyte*(x-1+width*(y-1))+2-c]=(byte)((Palette[Line1[il].Ind0]>>8*c) & 0XFF);
}
Mask.Grid[x+width*(y-1)]=Mask.Grid[x-1+width*(y-1)]=LabMask;
break;
} //:::::::::::::::::::::: end switch ::::::::::::::::::::::::::::::::
}
else
{ switch(dir)
{ case 0: Image.Grid[x+width*y]=Line1[il].Ind1;
Image.Grid[x+width*(y-1)]=Line1[il].Ind0;
Mask.Grid[x+width*y]=Mask.Grid[x+width*(y-1)]=LabMask;
break;
case 1: Image.Grid[x+width*y]=Line1[il].Ind0;
Image.Grid[x-1+width*y]=Line1[il].Ind1;
Mask.Grid[x+width*y]=Mask.Grid[x-1+width*y]=LabMask;
break;
case 2: Image.Grid[x-1+width*y]=Line1[il].Ind0;
Image.Grid[x-1+width*(y-1)]=Line1[il].Ind1;
Mask.Grid[x-1+width*y]=Mask.Grid[x-1+width*(y-1)]=LabMask;
break;
case 3: Image.Grid[x+width*(y-1)]=Line1[il].Ind1;
Image.Grid[x-1+width*(y-1)]=Line1[il].Ind0;
Mask.Grid[x+width*(y-1)]=Mask.Grid[x-1+width*(y-1)]=LabMask;
break;
} //:::::::::::::::::::::: end switch :::::::::::::::::::::::::::::::
} //------------------- end if (nBits==24) ---------------------------
} //============ end for (il < nLine1 …============================
int first, last;
int[] Shift = new int[]{0,2,4,6};
for (int il=0; il<nLine2; il++) //======================================
{
if ((il % i1) == 0) fm1.progressBar1.PerformStep();
if (il==0) first=0;
else first=Line2[il-1].EndByte+1;
last=Line2[il].EndByte;
x=Line2[il].x;
y=Line2[il].y;
int iByte=first, iShift=0;
iVect2 P = new iVect2(), PixelP = new iVect2(), PixelN = new iVect2(); // comb. coordinates
byte[] ColN = new byte[3], ColP = new byte[3];
byte[] ColStartN = new byte[3], ColStartP = new byte[3],
ColLastN = new byte[3], ColLastP = new byte[3]; // Colors
for (int c=0; c<3; c++)
ColN[c]=ColP[c]=ColStartN[c]=ColStartP[c]=ColLastN[c]=ColLastP[c]=0;
if (nBits==24)
{ for (int c=0; c<nbyte; c++)
{ ColStartN[2-c]=(byte)((Palette[Line2[il].Ind0]>>8*c) & 255);
ColStartP[2-c]=(byte)((Palette[Line2[il].Ind1]>>8*c) & 255);
ColLastN[2-c]= (byte)((Palette[Line2[il].Ind2]>>8*c) & 255);
ColLastP[2-c]= (byte)((Palette[Line2[il].Ind3]>>8*c) & 255);
}
}
else
{ ColStartN[0]=Line2[il].Ind0;
ColStartP[0]=Line2[il].Ind1;
ColLastN[0]=Line2[il].Ind2;
ColLastP[0]=Line2[il].Ind3;
}
P.X=Line2[il].x; P.Y=Line2[il].y;
int nCrack=Line2[il].nCrack;
int xx, yy; // Interpolation:
for (int iC=0; iC<nCrack; iC++) //=====================================
{ dir=(Byte[iByte] & (3<<Shift[iShift]))>>Shift[iShift];
switch(dir) // Standard coordinates
{ case 0: PixelP=P; PixelN=P+Step[3]; break;
case 1: PixelP=P+Step[2]; PixelN=P; break;
case 2: PixelP=P+Step[2]+Step[3]; PixelN=P+Step[2]; break;
case 3: PixelP=P+Step[3]; PixelN=P+Step[2]+Step[3]; break;
}
if (PixelP.Y<0 || PixelN.Y<0 || PixelP.Y>height-1 || PixelN.Y>height-1)
{ MessageBox.Show("Restore: Bad 'PixelP' or 'PixelN'. 'Byte' is bad. iByte=" +
iByte + "; dir=" + dir + "; Byte=" + Byte[iByte]);
}
for (int c=0; c<nbyte; c++) //====================================
{ ColN[c]=(byte)((ColLastN[c]*iC+ColStartN[c]*(nCrack-iC-1))/(nCrack-1));
ColP[c]=(byte)((ColLastP[c]*iC+ColStartP[c]*(nCrack-iC-1))/(nCrack-1));
} //=============== end for (c... ===============================
// Assigning colors to intermediate pixels of a line:
xx=PixelP.X; yy=PixelP.Y;
if (xx+width*yy>width*height-1 || xx+width*yy<0)
{
MessageBox.Show("Restore: Bad 'xx,yy'="+(xx+width*yy)+"; 'Byte' is bad.");
}
if (xx+width*yy<width*height && xx+width*yy>=0)
{ for (int c=0; c<nbyte; c++) Image.Grid[c+nbyte*xx+nbyte*width*yy]=ColP[c];
Mask.Grid[xx+width*yy]=LabMask;
}
xx=PixelN.X; yy=PixelN.Y;
if (xx + width * yy > width * height - 1 || xx + width * yy < 0) return -1;
if (xx+width*yy<width*height && xx+width*yy>=0)
{ for (int c=0; c<nbyte; c++) Image.Grid[c+nbyte*xx+nbyte*width*yy]=ColN[c];
Mask.Grid[xx+width*yy]=LabMask;
}
P=P+Step[dir];
iShift++;
if (iShift==4)
{ iShift=0;
iByte++;
}
} //=============== end for (iC... ==============================
} //================ end for (il < nLine2 ===========================
Mask.Grid[0]=Mask.Grid[width-1]=Mask.Grid[width*height-1]=
Mask.Grid[width*(height-1)]=LabMask;
return 1;
} //******************* end Restore ***********************************
下一个方法Image.Smooth( )，属于类CImage，作为图像RestoreIm的方法开始。它从平滑图像的边界开始。该方法设计用于处理彩色或灰度图像。如果Image是彩色图像，变量nbyte被设置为等于 3，或者如果是灰度图像，变量nbyte被设置为 1。因此，在第一种情况下，处理三个颜色通道，而在第二种情况下，仅处理一个灰度值。
该方法将颜色或灰度值内插在角和最终存在的跨越边界的线的裂缝之间。这在该方法的四个独立部分中完成，因为四个边界具有完全不同的坐标。
Smooth然后开始平滑图像Mask的零像素区域。它首先对行进行平滑。它跟踪图像的一行，并测试图像中成对的后续像素Mask。如果实际像素( x，y )为零，前一个像素非零，那么xbeg的值被设置为x - 1，变量ColorBeg被设置为Image的像素( x - 1， y )的颜色(或灰度值)。如果实际像素和先前像素都为零，则计算通过的像素。如果实际像素( x，y )非零，并且前一个像素为零，则已经找到以零开始和结束的行的一段。变量ColorEnd被设置为Image的像素( x 、 y )的颜色(或灰度值)，所有通过的零像素被填充以ColorBeg和ColorEnd之间的插值颜色(或灰度值)。
然后沿着柱执行类似的过程。在这种情况下获得的插值值与沿行平滑期间保存在Image像素中的值进行平均。
经过这两次插值，恢复的图像已经看起来不错了。它只显示一些比周围环境稍亮或稍暗的垂直和水平线条。方法Smooth通过平滑图像Mask具有零值的区域来移除这些线。为此，它实现了拉普拉斯偏微分方程的数字解。它使用在 Press，(1990)中描述为同时过度弛豫(SOR)的方法。这种方法不是最快的，但是非常简单，很容易编程。
这里是Smooth的代码。
public int Smooth(ref CImage Mask, Form1 fm1)
/* Calculates the average colors between "gvbeg" and "gvend" and saves them in the image "this". This is a digital automation with the states S=1 at Mask>0 and S=2 at Mask==0, but S is not used. The variable "mpre" has the value of "Mask" in the previous pixel. --*/
{
int c, cnt, LabMask = 250, msk, mpre, x, xx, xbeg, xend, ybeg, yend, y, yy;
int[] Col = new int[3], ColBeg = new int[3], ColEnd = new int[3];
int nbyte;
if (N_Bits == 24) nbyte = 3;
else nbyte = 1;
// Smoothing the borders:
// Border at y=0:
y = 0; cnt = 0; xbeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[c];
for (x = 0; x < width; x++) //======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //----------------------------------------
{
cnt = 1; xbeg = x - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x - 1 + width * y) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //----------------------------------------
{
cnt++; xend = x;
for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (xx = xbeg + 1; xx < xend; xx++) //============================
{
for (c = 0; c < nbyte; c++) Grid[nbyte * (xx + width * y) + c] =
(byte)((ColBeg[c] * (xend - xx) + ColEnd[c] * (xx - xbeg)) / cnt);
Mask.Grid[xx + width * y] = (byte)LabMask;
} //============== end for (xx... ==============================
}
mpre = msk;
} //=============== end for (x=0; ... ==============================
// Border at y=height-1:
y = height - 1; cnt = 0; xbeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * width * y + c];
for (x = 0; x < width; x++) //=======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //----------------------------------------
{
cnt = 1; xbeg = x - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x - 1 + width * y) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //----------------------------------------
{
cnt++; xend = x;
for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (xx = xbeg + 1; xx < xend; xx++)
{
for (c = 0; c < nbyte; c++) Grid[nbyte * (xx + width * y) + c] =
(byte)((ColBeg[c] * (xend - xx) + ColEnd[c] * (xx - xbeg)) / cnt);
Mask.Grid[xx + width * y] = (byte)LabMask;
}
}
mpre = msk;
} //=============== end for (x=0; ... ===============================
// Border at x=0
x = 0; cnt = 0; ybeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * 0) + c];
for (y = 0; y < height; y++) //======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //-------------------------------------------
{
cnt = 1; ybeg = y - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * (y - 1)) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //----------------------------------------
{
cnt++; yend = y;
for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (yy = ybeg + 1; yy < yend; yy++)
{
for (c = 0; c < nbyte; c++)
{
Col[c] = (ColBeg[c] * (yend - yy) + ColEnd[c] * (yy - ybeg)) / cnt;
Grid[nbyte * (x + width * yy) + c] = (byte)Col[c];
}
Mask.Grid[x + width * yy] = (byte)LabMask;
}
}
mpre = msk;
} //=============== end for (y=0; ... ==============================
// Border at x=width-1
x = width - 1; cnt = 0; ybeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * 0) + c];
for (y = 0; y < height; y++) //=======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //----------------------------------------
{
cnt = 1; ybeg = y - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * (y - 1)) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //----------------------------------------
{
cnt++; yend = y;
for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (yy = ybeg + 1; yy < yend; yy++)
{
for (c = 0; c < nbyte; c++)
{
Col[c] = (ColBeg[c] * (yend - yy) + ColEnd[c] * (yy - ybeg)) / cnt;
Grid[nbyte * (x + width * yy) + c] = (byte)Col[c];
}
Mask.Grid[x + width * yy] = (byte)LabMask;
}
}
mpre = msk;
} //=============== end for (y=0; ... ==============================
// End smoothing border; Smooth on "x":
fm1.progressBar1.Visible = true;
fm1.progressBar1.Value = 0;
int Sum = height + width + 50 * 10;
int denomProg = fm1.progressBar1.Maximum / fm1.progressBar1.Step;
int i1 = Sum / denomProg;
for (y = 0; y < height; y++) //=======================================
{
if ((y % i1) == 0) fm1.progressBar1.PerformStep();
cnt = 0; xbeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * width * y + c];
for (x = 0; x < width; x++) //======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //----------------------------------------
{
cnt = 1; xbeg = x - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x - 1 + width * y) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //----------------------------------------
{
cnt++; xend = x;
for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (xx = xbeg + 1; xx < xend; xx++)
{
for (c = 0; c < nbyte; c++) Grid[nbyte * (xx + width * y) + c] =
(byte)((ColBeg[c] * (xend - xx) + ColEnd[c] * (xx - xbeg)) / cnt);
}
}
mpre = msk;
} //=============== end for (x=0; ... =============================
} //================ end for (y=0; ... ==============================
// Smooth on "y":
for (x = 0; x < width; x++) //========================================
{
if ((x % i1) == 0) fm1.progressBar1.PerformStep();
cnt = 0; ybeg = 0; mpre = 200;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * 0) + c];
for (y = 0; y < height; y++) //======================================
{
msk = Mask.Grid[x + width * y];
if (mpre > 0 && msk == 0) //----------------------------------------
{
cnt = 1; ybeg = y - 1;
for (c = 0; c < nbyte; c++) ColBeg[c] = Grid[nbyte * (x + width * (y - 1)) + c];
}
if (mpre == 0 && msk == 0) //----------------------------------------
{
cnt++;
}
if (mpre == 0 && msk > 0) //-----------------------------------------
{
cnt++; yend = y; for (c = 0; c < nbyte; c++) ColEnd[c] = Grid[nbyte * (x + width * y) + c];
for (yy = ybeg + 1; yy < yend; yy++)
{
for (c = 0; c < nbyte; c++)
{ Col[c]= (Grid[nbyte*(x+width*yy)+c]+(ColBeg[c]*(yend-yy) +
ColEnd[c]*(yy-ybeg))/cnt)/2;
Grid[nbyte * (x + width * yy) + c] = (byte)Col[c];
}
}
}
mpre = msk;
} //=============== end for (y=0; ... =============================
} //================= end for (x=0; ... =============================
// Solving the Laplace's equation:
int i;
double fgv, omega = 1.4 / 4.0, dMaxLap = 0.0, dTH = 1.0;
double[] dGrid = new double[width * height * nbyte];
double[] Lap = new double[3];
for (i = 0; i < width * height * nbyte; i++) dGrid[i] = (double)Grid[i];
fm1.progressBar1.Visible = false;
fm1.progressBar1.Value = 0;
fm1.progressBar1.Visible = true;
int it1=10;
for (int iter = 0; iter < 50; iter++) //================================
{ // Smooth Math.Abs((x-y))%2==0
if ((iter % it1) == 0) fm1.progressBar1.PerformStep();
for (y = 1; y < height - 1; y++)
for (x = 1; x < width - 1; x++)
{
if (Mask.Grid[x + width * y] == 0 && Math.Abs((x - y)) % 2 == 0)
for (c = 0; c < nbyte; c++)
{
Lap[c] = 0.0;
Lap[c] += dGrid[nbyte * (x + width * (y - 1)) + c];
Lap[c] += dGrid[nbyte * (x - 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + width * (y + 1)) + c];
Lap[c] -= 4.0 * dGrid[nbyte * (x + width * y) + c];
fgv = dGrid[nbyte * (x + width * y) + c] + omega * Lap[c];
if (fgv > 255.0) fgv = 255.0;
if (fgv < 0.0) fgv = 0;
dGrid[nbyte * (x + width * y) + c] = fgv;
}
}
// Smooth at Math.Abs((x-y))%2==1
for (y = 1; y < height - 1; y++)
for (x = 1; x < width - 1; x++)
{
if (Mask.Grid[x + width * y] == 0 && Math.Abs((x - y)) % 2 == 1)
for (c = 0; c < nbyte; c++)
{
Lap[c] = 0.0;
Lap[c] += dGrid[nbyte * (x + width * (y - 1)) + c];
Lap[c] += dGrid[nbyte * (x - 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + width * (y + 1)) + c];
Lap[c] -= 4.0 * dGrid[nbyte * (x + width * y) + c];
fgv = dGrid[nbyte * (x + width * y) + c] + omega * Lap[c];
if (fgv > 255.0) fgv = 255.0;
if (fgv < 0.0) fgv = 0;
dGrid[nbyte * (x + width * y) + c] = fgv; //(int)(fgv);
}
}
dMaxLap = 0.0; // Calculating MaxLap:
for (y = 1; y < height - 1; y++)
for (x = 1; x < width - 1; x++) //===================================
{
if (Mask.Grid[x + width * y] == 0) //------------------------------
{
for (c = 0; c < nbyte; c++) //==============================
{
Lap[c] = 0.0;
Lap[c] += dGrid[nbyte * (x + width * (y - 1)) + c];
Lap[c] += dGrid[nbyte * (x - 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + 1 + width * y) + c];
Lap[c] += dGrid[nbyte * (x + width * (y + 1)) + c];
Lap[c] -= 4.0 * dGrid[nbyte * (x + width * y) + c];
if (Math.Abs(Lap[c]) > dMaxLap) dMaxLap = Math.Abs(Lap[c]);
} //================= end for (c=0; =====================
} //------------------------------ end if (Mask... ---------------
} //================== end for (x=1; ... ==========================
int ii;
for (ii = 0; ii < width * height * nbyte; ii++) Grid[ii] = (byte)dGrid[ii];
if (dMaxLap < dTH) break;
} //==================== end for (iter... ============================
return 0;
} //********************** end Smooth **********************************
九、图像分割和连通分量
图像分割是图像分析中的一个重要步骤。以下是关于分割的已知信息(参见 Shapiro (2001 年):
在计算机视觉中*、、、图像分割、、是将一幅数字图像分割成多个片段的过程(集合像素的*，也称为超像素)。分割的目标是简化和/或改变图像的表示，使其更有意义，更易于分析。【1】【2】图像分割通常用于定位物体和边界(直线、曲线等)。)在图像中。更准确地说，图像分割是给图像中的每个像素分配标签的过程，使得具有相同标签的像素共享某些特征。**
在科学文献中有大量关于分割方法的出版物。然而，让我们记住，任何阈值图像可以被认为是一个分段的。量化图像也可以被认为是分割的，其中所有灰度值或颜色被细分为相对较少数量的组，单一颜色被分配给一个组的所有像素。分段问题的困难仅在于定义有意义的分段的可能性；即从人类感知的角度来看有意义的片段。例如，对于一个分段，我们可以说，“这个分段是房子，这个分段是街道，这个分段是汽车，等等”，这将是一个有意义的分段。据我所知，开发一个有意义的分割是一个尚未解决的问题。在这一节中，我们考虑通过将图像的颜色量化成最佳地代表图像颜色的少量组来进行分割。
通过量化颜色进行分割
量化颜色是众所周知的压缩彩色图像的方法:颜色信息不是直接由图像像素数据携带，而是存储在称为调色板的独立数据中，调色板是颜色元素的阵列。数组中的每个元素都代表一种颜色，按其在数组中的位置进行索引。图像像素不包含颜色的完整规格，而仅包含其在调色板中的索引。这是众所周知的产生索引 8 位图像的方法。
我们的方法包括产生一个对具体图像最佳的调色板。为此，我们开发了第八章中描述的方法MakePalette，并成功用于项目WFcompressPal。
连接的组件
我们假设这些片段应该是像素的连接子集。连通性是一个拓扑概念，而经典拓扑主要考虑无限点集，其中一个点的每个邻域都包含无限多个点。经典拓扑学不适用于数字图像，因为它们是有限的。
连通性的概念也用在图论中。可以将数字图像视为顶点为像素的图形，图形的边连接任意两个相邻的像素。路径是一个像素序列，其中除了起始处的一个像素和末端的一个像素之外，每个像素都恰好与两个相邻的像素相关联。
众所周知(例如，参考 Hazewinkel，(2001))，在无向图中，如果路径从 x 通向 y，则无序的顶点对{x，y}称为连通的。否则，无序的顶点对称为不连通的。
连通图是一种无向图，其中图中的每一对无序顶点都是连通的。否则，它被称为不连通图。
将图论中的连通性概念应用到图像上是可能的。我们考虑两种邻接:四邻接(图 9-1a )和八邻接(图 9-1b )。在图 9-1a 的情况下，黑色像素和白色像素的子集都是断开的；在图 9-1b 的情况下，它们都被连接。在这两种情况下，都存在一个众所周知的连通性悖论:为什么图 9-1a 中的黑色像素应该是断开的，而它们之间没有任何东西？同样，为什么图 9-1b 中的黑色像素是相连的，尽管它们之间存在从一个白色像素到另一个白色像素的路径？
图 9-1
连通性悖论:(1)四邻接，和(2)八邻接
Rosenfeld (1970)建议在二进制图像中考虑黑色像素的八个邻接和白色像素的四个邻接。这个建议解决了二进制图像的问题，但是它不适用于多色图像。据我所知，多色图像连通性悖论的唯一解决方案是将数字图像视为一个细胞复合体，如第七章所述。然后，连通性可以通过下面指定的等式规则(参见 Kovalevsky (1989))被正确地定义，并且没有悖论。
坐标为(x，y)和(x + 1，y + 1)或(x，y)和(x–1，y + 1)的两个像素组成一对对角线。在 ACCs 理论中，如果两个像素具有相同的颜色，并且位于它们之间的低维单元(一维或零维)也具有该颜色，则一对八个相邻像素被定义为连通的。自然，谈论低维细胞的颜色没有什么意义，因为它们的面积为零。然而，有必要使用这个概念来一致地定义数字图像中的连通性。为了避免使用低维单元的颜色的概念，我们建议当像素的标签等于其颜色时，用标签代替颜色。
考虑图 9-2 的例子。假设有必要以白 V 和黑 V 都连接的方式来定义连通性。这在任何邻接关系下显然都是不可能的。该目的可以通过将隶属标签(例如，对应于颜色的整数)分配给较低维度的单元的等式规则来实现。
图 9-2
一幅图像中的白色和黑色 V 形区域，由于应用了等式规则，这两个区域相互连接
二维图像的等式规则表明，一维单元总是接收与其关联的两个主要(即二维)单元中最大的标签(较浅的颜色)。
零维单元格 c 0 的标签定义如下:
- 如果 c 0 的 2 × 2 邻域恰好包含一对具有相同标号(Equ)的对角像素，则 c 0 接收该标号。如果有多个这样的对，但其中只有一个属于窄(Na)条带，则c0 接收窄条带的标签。否则c0 接收最大值；即c0 邻域内像素的较亮(Li)标号(即较大标号)。
后一种情况对应于当 c 0 的邻域不包含具有相同标号的对角对或者包含多于一个这样的对角对并且其中多于一个对角对属于窄条的情况。
为了决定在 c 0 附近的对角线对是否属于窄条纹，有必要扫描 4 × 4 = 16 像素的阵列，其中 c 0 在中间，并且计数对应于两条对角线的标签。最小的计数表示窄条纹。其他有效成员规则的例子可以在 Kovalevsky (1989)中找到。
对于分割图像的分析来说，找到并区分图像的所有连接部分是很重要的。问题是找到图像的最大连通子集，每个子集包含恒定颜色的像素。然后，有必要对组件进行编号，并为每个像素分配其所属组件的编号。让我们称分配的数字为标签。
如果对于 S 中的任意两个像素，在 S 中存在包含这两个像素的相邻像素序列，则某种颜色的像素子集 S 被称为连通的。我们认为埃奎纳利街区。
为了找到并标记组件，图像必须包含相对较少的颜色。否则，彩色图像的每个像素可能是最大连通分量。分割图像大多满足拥有少量颜色的条件。
简单的方法包括逐行扫描图像，并为在已经扫描的子集中没有相同颜色的相邻像素的每个像素分配下一个未使用的标签。如果下一个像素在已经扫描的子集中具有相同颜色的相邻像素，并且它们具有不同的标签，则将这些不同标签中最小的标签分配给该像素。具有这些不同标签之一的已扫描子集的所有像素必须获得这个最小的标签。在最坏的情况下，替换可能影响高达高度 2 /2 像素，其中高度是图像中的行数。在扫描了整个图像之后，替换的总数可以是大约高度 3 /2。这个数字可能非常大，所以这种幼稚的方法是不实际的。标记一个分量的像素的有效方法是众所周知且广泛使用的图遍历方法。
图的遍历算法及其代码
现在让我们考虑图遍历方法。图像被视为图形，像素被视为顶点。如果对应的像素具有相等的值(例如，颜色)并且相邻，则两个顶点由图的边连接。邻接关系可以是四邻接或八邻接。很容易认识到，图形的分量正好对应于图像的期望分量。
至少有两种众所周知的遍历一个图的所有顶点的方法:深度优先搜索和宽度优先搜索算法。它们彼此非常相似，这里我们只考虑广度优先搜索。这种算法被用于我们抑制脉冲噪声的项目中(第二章)。我们在这里详细描述一下。
广度优先算法采用了众所周知的数据结构，称为 *queue，*也称为先进先出:值可以放入队列，先放入的值将首先从队列中出来。该算法的过程如下。
-
将组件的编号 NC 设置为零。
-
取任意未标记的顶点 V ，增加 NC ，用值 NC 标记 V ，并将 V 放入队列。顶点 V 是具有索引 NC 的新组件的“种子”。
-
当队列不为空时，重复以下说明:
-
从队列中取出第一个顶点。
-
测试 W 的每个邻居 N 。如果 N 被贴上标签，忽略；否则用 NC 标记 N 并将其放入队列。
-
-
当队列为空时，那么组件 NC 已经被标记。
继续第 2 步，找到下一个组件的种子，直到所有顶点都被标记。
让我们描述一下这个算法的伪代码。给定一个包含N个元素的多值数组Image[]以及方法Neighb(i,m)和Adjacent(i,k)，第一个返回第i个元素的第m个邻居的索引。如果第i个元素与第k个元素和FALSE otherwise相邻，第二个函数返回TRUE。作为标记的结果，每个元素获得其所属的连接组件的标记。标签保存在一个数组Label[]中，其大小与Image[]相同。我们假设所有方法都可以访问数组Image和Label以及它们的大小N。
广度优先算法的伪代码
分配数组Label[N]。Label的所有元素现在必须用零初始化，这意味着这些元素还没有被标记。
for (i=0; i<N; i++) Label[i]=0; // setting all labels to 0
int NC=0; // number of components
for (V=0; V<N; V++) // loop through all elements of the image
{
if (Label[V]≠0) continue; // looking for unlabeled elements
NC=NC+1; // index of a new component
Label[V]=NC; // labeling the element with index V with NC
PutIntoQueue(V); // putting V into the queue
while(QueueNotEmpty()) //loop running while queue not empty
{ W=GetFromQueue();
for (n=0; n<Max_n; n++) // all neighbors of W ==========
{ N=Neighb(W,n); // the index of the nth neighbor of W
if (Label[N]==0 AND Adjacent(W,N)==TRUE
AND Image[N]==Image[W])
{ Label[N]=NC; // labeling element with index N with NC
PutIntoQueue(N);
}
} // ========= end for (n=0; ... =================
} // =========== end while ========================
} // ============= end for (V =0;... ======================
算法结束
在这种情况下，值Max_n是图像元素的所有邻居的最大可能数量，而不管它是否已经被访问过。它等于 8。需要用方法PutIntoQueue(i), GetFromQueue(),和QueueNotEmpty()提供队列数据结构。如果队列不为空，最后一个方法返回TRUE；否则它返回FALSE。
该算法标记包含起始像素的组件的所有像素。为了标记图像的所有成分，有必要在那些没有被标记为属于已经被处理的成分的像素中选择下一个开始像素。这必须重复，直到图像的所有像素都被标记。
等价类方法
还有另一种方法来解决标记连通分量的问题，即等价类(EC ),它考虑同时标记图像的所有分量。Rosenfeld 和 Pfalz (1966)是关于 EC 方法的最早出版物之一。在二维图像的情况下，该算法的思想如下进行。该算法逐行扫描图像两次。它在每个像素 P 的第一次扫描期间，根据四个或八个相邻像素并具有与 P 相同的颜色，调查与 P 相邻的所有已访问像素的集合 S 。如果没有这样的像素，那么 P 获得下一个未使用的标签。否则 P 得到最小的标签(标签集必须是有序的；整数大多用作 S 中像素的标号。同时，该算法声明 S 的像素的所有标签等同于最小的标签。为此目的，该算法采用等价表，其中每一列对应于所采用的标签之一。等效标签存储在列的元素中。
当图像被完全扫描后，一种特殊的算法处理该表并确定等效标签的类别。每个类别由最小的等价标签表示。重新排列该表，使得在对应于每个标签的列中，保留该类的标签。在第二次扫描期间，图像中的每个标签都被该类的标签所替换。
该算法的缺点之一是难以先验地估计等价表的必要大小。Kovalevsky (1990)独立于 Gallert 和 Fisher (1964)开发了 EC 方法的改进版本，称为根方法 *。*它采用了一个标签数组L，其元素数量与获得的图像中的像素数量相同。这是等价表的替代品。为了解释改进的思想，我们首先考虑二维二值图像的最简单情况，尽管该算法也适用于多值和多维图像的情况，而没有本质的改变。
建议的方法基于以下想法。图像的每个像素都有一个标签数组L的对应元素。我们姑且称之为像素的标签。数组L的类型必须能够在其元素中保存图像每个像素的索引x + width*y(类型integer对于任何大小的图像都足够)。如果一个组件包含一个像素，其标签等于该像素的索引，则该标签称为该组件的根。
在算法开始时，对应于像素( x，y )的数组L的每个元素被设置为等于索引x+width*y，其中width是图像的列数。因此，在开始时，每个像素被认为是一个组件。
该算法逐行扫描图像。如果像素在先前扫描的区域中没有相同颜色的相邻像素，则该像素保留其标签。然而，如果下一个像素 P 在已经访问过的像素中具有相同颜色的相邻像素，并且它们具有不同的标签，这意味着这些像素属于具有不同根的不同的先前发现的分量。这些不同的部件必须结合成一个单一的部件。为此，需要找到这些分量的根，并将所有这些根和标签 P 设置为最小根。重要的是要注意，只有根被改变，而不是与 P 相邻的像素的标签值被改变。像素的标签将在第二次扫描中改变。这就是这种方法如此快速的原因。在运行像素处相遇的所有分量得到相同的根，并且它们被合并成单个分量。不是根的标签保持不变。在前面提到的最坏情况下，替换可能只影响width /2 个像素。
在第一次扫描之后，每个像素的标签指向具有较小标签的相同组件的像素。因此，有一系列越来越小的标签以组件的根结束。在第二次扫描期间，所有标签都被后续的整数替换。
子集的连通性必须由邻接关系来定义。这可以是公知的( n ， m )邻接关系，其仅被定义用于二进制图像，或者是基于由方法EquNaLi定义的子集中较低维度的单元的成员的邻接关系。
可以看出，在根方法中使用的标签数组L需要太多的存储空间:为了存储索引，L的每个元素中的位数必须不小于log2(size)，其中size是图像元素的数量。然而，这个大小的L对于任何标记连通分量的方法都是必要的，如果该方法被认为适用于任何多值图像的话。事实上，图像中的每个元素都是一个组成部分。不同标签的数量等于图像元素的数量。例如图 9-3 所示的四种颜色的二维图像，在任何邻接关系下都具有这种性质。
图 9-3
具有四种颜色的图像示例，其组件数量与像素数量相同
如果已知所用标签的数量小于 size，那么每个元素具有较少比特数的数组 L 足以应用使用等价表的方法。然而，这对于现代计算机来说并不重要，因为现代计算机有充足的内存。因此，例如，在具有 32 位处理器的现代计算机上，数据类型 int 的数组可以用于包含多达 2 个 32 个元素的图像。因此，二维图像可以具有 65，536 × 65，536 像素的大小，而三维图像可以具有 2，048 × 2，048 × 1，024 个体素的大小。
让我们回到算法上来。它逐行扫描图像，并为每个像素 P 调查与 P 相邻的已经访问过的像素的集合 S 。如果在 S 中存在与 P 颜色相同的像素，并且它们具有不同的标号，那么算法找到这些像素的最小根标号，并将其分配给 P 和 S 的所有像素的根。因此，在连续像素处彼此相遇的所有分量获得相同的根，并被合并成单个分量。不是根的像素的标签保持不变。在单个组件的最坏情况下，替换只能影响width /2 个像素。
每个像素在第一次扫描后属于通向其分量的根的路径。在第二次扫描期间，根和所有其他像素的标记被随后的整数替换。
一个简单的二维例子如图 9-4 所示。给定的是一个N = 16 像素的二维二进制数组Image[N]。暗前景的两个像素相邻，当且仅当它们具有共同的入射 1 单元(四个相邻)，而亮背景的像素相邻，如果它们具有任何共同的入射单元(八个相邻)。在算法开始时，所有像素的标签等于它们的索引(图 9-4a 中从 0 到 15 的数字)。索引为 0 和 1 的像素在之前访问的像素中没有相同颜色的相邻像素。它们的标签 0 和 1 保持不变。它们是根。索引为 2 的像素具有相邻的像素 1，其根为 1。因此，标签 2 被 1 代替。在扫描第一行之后，标签 3 保持不变。根据前面的规则改变所有像素的标记，直到到达最后一个像素 15。该像素与像素 11 和 14 相邻。像素 11 的标签指向根 3；像素 14 的指向像素 0，像素 0 也是一个根。标签 0 小于 3。因此，根 3 的标签被更小的标签 0 取代(图 9-4b )。
在第二次运行期间，根的标记被随后的自然数 1 和 2 替换。所有其他像素的标签获得其根的新标签(图 9-4c )。EC 算法的代码将在后面描述。
图 9-4
标记连接组件的算法图示
根算法的伪代码
该算法可以应用于组合网格或标准网格中的图像。在组合网格中，其中两个单元的关联关系由它们的组合坐标来定义，从像素值集合中分配一个值(颜色)给任何维度的每个单元就足够了。具有相同值的两个事件单元是相邻的，并且将被分配给相同的组件。在标准网格中，必须给出一种方法来指定给定元素 e 的邻域中的哪些网格元素与 e 相邻，因此如果它们具有相同的值，则必须将其分配给相同的组件。在我们早期的简单二维例子中，我们采用了前景像素的四个相邻像素和背景像素的八个相邻像素。一般来说， n 维复合体的主单元的邻接必须由指定较低维单元的成员的规则来指定(例如，由类似于等式规则的规则来指定)，因为众所周知的( m ， n )邻接不适用于多值图像。我们在这里考虑标准网格情况下的伪代码。
给定的是一个包含N个元素的多值数组Image[]，以及方法Neighb(i,m)和Adjacent(i,k)。第一个函数返回第i个元素的第m个邻居的索引，如果第i个元素与第k个元素相邻，第二个函数返回TRUE。否则它返回FALSE。在二维情况下，这可以是实现等式规则的方法。作为标记的结果，每个元素获得它所属的连接组件的标记。标签保存在另一个数组Label[].
数组Label[N]的大小必须与Image[N].相同Label的每个元素必须至少有log2N位，其中N是Image的元素个数，并且必须用自己的索引进行初始化:
for(I = 1；I < N；i++)标签[i] = i。
为了使伪代码更简单，我们假设所有方法都可以访问数组Image and Label和它们的大小N。
在算法的第一个循环中，Image的每个元素都有一个指向其组件根的标签。值Max_n是已经被访问的元素(像素或体素)的相邻元素的最大可能数量。很容易识别出Max_n等于(3 d - 1)/2，其中 d 为图像的尺寸。例如，在二维情况下，Max_n等于 4，而在三维情况下，它等于 13 .
求根算法
for (i=0; i<N; i++)
{ for (m=0; m<Max_n; m++)
{ k=Neighb(i,m); // the index of the mth adjacent element of i
if(Adjacent(i,k) AND Image[k] == Image[i]) SetEqu(i, k);
}
} // end of the first loop
SecondRun();
子程序SetEqu(i,k)(其伪代码见下)将元素i和k的根的索引相互比较，并将具有较大索引的根的标签设置为等于较小索引。因此，两个词根属于同一个成分。方法Root(k)返回索引序列中的最后一个值，其中第一个索引k是给定元素的索引，下一个是Label[k],的值，依此类推，直到Label[k]等于k。子程序SecondRun()根据Label[k]是否等于k用元件计数器的值或用k的根代替Label[k]的值。
子程序的伪代码
subroutine SetEqu(i,k)
{ if (Root(i)<Root(k)) Label[Root(k)]=Root(i);
else Label[Root(i)]=Root(k);
} // end subroutine.
int Root(k)
{ do
{ if (Label[k]==k) return k;
k=Label[k];
} while(1); // loop runs until condition Label[k]==k is fulfilled
} // end Root
subroutine SecondRun()
{ count=1;
for (i=0; i<N; i++)
{ value=Label [i];
if (value==i )
{ Label[i]=count; count=count+1;
}
else Label[i]=Label[value];
}
} // end subroutine.
项目WFsegmentAndComp
我们已经开发了一个 Windows 窗体项目WFsegmentAndComp分割一个图像，并标记其连接组件的像素。
有可能以两种方式执行分割:或者通过将颜色聚类成相对小数量(例如，20 个)的聚类并将属于该聚类的所有像素的平均颜色分配给每个聚类，或者通过将彩色图像转换成灰度图像，将灰度值的全范围[0，255]细分成小数量的区间，并将平均灰度值分配给每个区间。这种转换不会导致信息的丢失，因为相连分量的标记是相对较大的数字，与原始颜色的颜色无关:具有某种颜色的像素集合可以由许多分量组成，每个分量具有另一个标记。
连通分量的标记是通过两种方法完成的，即广度优先搜索和根方法，以便用户可以比较这两种方法的速度。标记的结果被表示为具有 256 种颜色的人工调色板的彩色图像，并且调色板的索引是将标签除以 256 的结果。图 9-5 显示了项目的形式。
图 9-5
项目形式WFsegmentAndComp
当用户单击打开图像时，他或她可以选择要打开的图像。该图像将显示在左侧的图片框中，并将定义五个工作图像。
然后，用户应该单击下一个按钮“Segment”。原始图像将被转换为灰度图像。灰度值将除以 24，并且除法的整数结果将乘以 24。这样，新的灰度值变成 24 的倍数，灰度值的标度被量化成 256/24 + 1 = 11 个区间。分割的图像显示在右-右图片框中。
根据经验可知，即使是分割图像，分量的数量也是非常大的。在这些组件中，许多组件由单个像素或几个像素组成。图像分析对这些成分不感兴趣。为了减少小分量的数量，我们决定抑制分割图像的脉冲噪声。这将特定灰度级的许多小点与周围区域合并，并减少了小组件的数量。
用户可以选择要删除的暗点和亮点的最大尺寸，然后单击脉冲噪声。将呈现具有删除的暗点和亮点的图像，而不是分割的图像。
然后，用户可以单击广度优先搜索或根方法。带标签的图像将以错误的颜色显示；也就是说，使用包含 256 种颜色的人工调色板。调色板的索引等于标签除以 256。
第二章介绍了抑制脉冲噪声的代码。
这里，我们提出了用于标记分割图像的连接部分的宽度优先搜索的代码。如前所述，该算法处理包含起始像素的一个分量的像素。为了处理图像的所有组成部分，我们开发了方法LabelC，它扫描图像并测试每个像素是否被标记。如果不是，那么它是下一个组件的起始像素，并且调用方法Breadth_First_Search。
所有这些方法都是一个特殊类CImageComp的元素，它与类CImage的不同之处在于标签数组Lab的存在。属性nLoop和DenomProg用于控制指示方法LabelC进程的progressBar。
public class CImageComp
{
public byte[] Grid;
public int[] Lab;
public Color[] Palette;
public int width, height, N_Bits, nLoop, denomProg;
... // here are the constructors and the methods
}
下面是LabelC的代码。
public int LabelC(Form1 fm1)
{ int index, lab=0;
for (index=0; index<width*height; index++) Lab[index]=0;
int y1 = nLoop * height / denomProg;
for (int y=0; y<height; y++) //========================
{
if (y % y1 == 0) fm1.progressBar1.PerformStep();
for (int x = 0; x < width; x++) //=================
{ index=x+width*y;
if (Lab[index]==0)
{ lab++;
Breadth_First_Search(index, lab);
}
} //================= end for (int x...============
} //=================== end for (int x...==============
return lab;
} //********************* end LabelC ***************************
可以看出，该方法非常简单:它扫描图像的所有像素，并检查实际像素是否已经被标记。如果不是，则调用方法Breadth_First_Search。
这里是Breadth_First_Search的代码。
public int Breadth_First_Search(int root, int lab)
{ int light, gvP, Index, Len=1000, Nei;
bool rv;
iVect2 P = new iVect2(0, 0);
CQue Q = new CQue(Len);
light=Grid[root];
if (Lab[root]==0) Lab[root]=lab;
Q.Put(root);
while(!Q.Empty()) //=============================
{ Index=Q.Get();
gvP=Grid[Index];
P.Y=Index/width; P.X=Index-width*P.Y;
for (int n=0; n<8; n++) //====================
{ Nei=NeighbEqu(Index, n, width, height);
rv=EquNaliDir(P, n, gvP);
if (Nei<0) continue;
if (Grid[Nei]==light && Lab[Nei]==0 && rv)
{
Lab[Nei]=lab;
Q.Put(Nei);
}
} //============== end for (n... ==============
} //================ end while ====================
return 1;
} //****************** end Breadth_First_Search *************
为了显示组件标记的结果，我们使用一个简单的调色板，有 256 种颜色。调色板由方法MakePalette产生。这是它的代码。
public int MakePalette( ref int nPal)
{ int r, g, b;
byte Red, Green, Blue;
int ii=0;
for (r=1; r<=8; r++) // Colors of the palette
for (g=1; g<=8; g++)
for (b=1; b<=4; b++)
{ Red=(byte)(32*r); if (r==8) Red=255;
Green = (byte)(32 * g);
if (g == 8) Green = 255;
Blue= (byte)(64*b); if (b==4) Blue=255;
Palette[ii] = Color.FromArgb(Red, Green, Blue);
ii++;
}
nPal=4*ii;
return 1;
} //***************** end MakePalette *****************
类别CImageComp的图像BreadthFirIm包含组件和调色板的标签。为了使标签可见，我们通过方法LabToBitmap将图像BreadthFirIm转换为位图BreadthBmp，并将该位图显示在右边的图片框中。下面是LabToBitmap的代码。
public void LabToBitmap(Bitmap bmp, CImageComp Image)
{
if (Image.N_Bits != 8)
{
MessageBox.Show("Not suitable format of 'Image'; N_Bits=" + Image.N_Bits);
return;
}
switch (bmp.PixelFormat)
{
case PixelFormat.Format24bppRgb: nbyte = 3; break;
case PixelFormat.Format8bppIndexed:
MessageBox.Show("LabToBitmap: Not suitable pixel format=" + bmp.PixelFormat);
return;
default: MessageBox.Show("LabToBitmap: Not suitable pixel format=" + bmp.PixelFormat);
return;
}
Rectangle rect = new Rectangle(0, 0, bmp.Width, bmp.Height);
BitmapData bmpData = bmp.LockBits(rect, ImageLockMode.ReadWrite, bmp.PixelFormat);
IntPtr ptr = bmpData.Scan0;
int size = bmp.Width * bmp.Height;
int bytes = Math.Abs(bmpData.Stride) * bmp.Height;
byte[] rgbValues = new byte[bytes];
Color color;
int index = 0;
int y1 = nLoop * bmp.Height / 220; // denomProg;
for (int y = 0; y < bmp.Height; y++)
{
if (y % y1 == 0) progressBar1.PerformStep();
for (int x = 0; x < bmp.Width; x++)
{
color = Image.Palette[Image.Lab[x + bmp.Width * y] & 255];
index = 3 * x + Math.Abs(bmpData.Stride) * y;
rgbValues[index + 0] = color.B;
rgbValues[index + 1] = color.G;
rgbValues[index + 2] = color.R;
}
}
System.Runtime.InteropServices.Marshal.Copy(rgbValues, 0, ptr, bytes);
bmp.UnlockBits(bmpData);
} //**************************** end LabToBitmap ***********************
位图BreadthBmp用假颜色显示组件。
用户可以通过单击“保存宽度图像”将带有标记组件的图像保存为彩色图像。
如前所述，我开发了另一种标记连接组件的方法。这种方法的优点是不用一个接一个地标记组分，而是同时并行地标记。这个方法的思想在前面已经描述过了。这里我们要呈现的是重要方法ComponentE的代码，当用户点击 Root method 时，这个方法作为图片RootIm的方法被调用。
public int ComponentsE(Form1 fm1)
{ /* Labels connected components of "this" in "this->Lab" and returns
the number of components. Uses the EquNaLi connectedness. --*/
int Dir, light, i, nComp=0, rv, x1 = 0, y1 = 0;
bool adjac;
iVect2 P = new iVect2(0, 0);
for (int k=0; k<width*height; k++) Lab[k]=k;
int y2 = nLoop*height / denomProg;;
for (int y=0; y<height; y++) //==========================================
{ if ((y % y2) == 0) fm1.progressBar1.PerformStep();
for (int x = 0; x < width; x++) //=====================================
{ i=x+width*y; // "i" is the index of the actual pixel
light=Grid[i];
P.X=x; P.Y=y;
int nmax;
if (y==0) nmax=0;
else nmax=3;
for (int n=0; n<=nmax; n++) //==== the four preceding neighbors ========
{ switch(n)
{ case 0: x1=x-1; y1=y; break; // West
case 1: x1=x-1; y1=y-1; break; // North-West
case 2: x1=x; y1=y-1; break; // North
case 3: x1=x+1; y1=y-1; break; // North-East
}
if (x1<0 || x1>=width) continue;
Dir=n+4;
int indN=x1+width*y1;
int lightNeib=Grid[indN];
if (lightNeib!=light) continue;
if ((Dir & 1)==0) adjac=true;
else
adjac=EquNaliDir(P, Dir, light);
if (adjac)
rv=SetEquivalent(i, indN); // Sets Lab[i] or Lab[indN] equal to Root
} //==================== end for (int n ... ========================
} //===================== end for (int x ... ========================
} //====================== end for (int y ... ========================
nComp=SecondRun(); // Writes indexes of components from 1 to nComp to "Comp".
return nComp; // nComp+1 is the number of indexes
} //************************* end ComponentsE **************************
类别CImageComp的图像RootIm包含组件和调色板的标签。就像在Breadth_First_Search的情况中一样，方法MakePalette和LabToBitmap被调用以使标签可见。图像RootIm被转换为位图RootBmp，后者显示在右边的图片框中。然后，用户可以通过单击“根”的保存图像来保存该图像。
结论
根算法和广度优先算法的时间复杂度与图像元素的数量成线性关系。只要定义了邻接关系，这两种算法都可以应用于任何维数的图像。
作者进行的大量实验表明，根算法和广度优先算法的运行时间几乎相等。根算法的优点是它测试图像的每个元素是否已经被标记(3 d - 1)/2 次，其中 d 等于 2 或者等于 3 是图像的维度，而宽度优先算法测试它 3 d 次，这大约是两倍多。广度优先算法可以直接用于只标记包含图像的给定元素的单个分量，而根算法必须同时标记所有分量。
