十、拉直绘画照片
我们有时对拍摄展览中展出的画作感兴趣。如果房间是纯照明，那么有必要使用闪光灯。然而，如果你把相机直接放在画的前面，闪光灯会在画的窗格中反射，照片就被破坏了。为了避免反光，有必要从侧面拍照。然后，然而，画的图像是扭曲的(图 10-1 )。
图 10-1
从侧面拍摄的绘画照片
我建议一个对这类照片进行整改的项目。如果我们知道相机与画的距离，从画的中点到相机物镜的光线方向，以及相机的焦距，那么中心投影的几何知识将是有用的。估计所有这些参数是困难的，特别是如果在摄影过程中使用变焦，因为变焦会改变相机的焦距。
因此，我建议一种基于画面框架是矩形这一假设的方法。然后，当使用照片中框架左边缘的长度与其右边缘的长度的关系时，可以估计框架的宽度与其高度的关系。上边缘的长度与下边缘的长度的关系也是有用的。
名为WFrectify的建议项目工作如下。可以使用.jpg和.bmp图像；其他格式的图像应该转换成这些格式之一。使用 IrfanView 程序很容易做到这一点，该程序是免费的，可以从互联网上下载。项目形式如图 10-2 所示。
图 10-2
项目的形式
要打开图像，用户应该单击打开图像，然后单击角。两条白线的直角标记了图像的左下四分之一。然后，用户应该决定他或她是否想要有框或无框的图片。根据这一决定，用户单击框架或框架内图像的左下角。
点选的角落会以小红点标示，白色的角会移到影像的左上角。用户应该单击相应的角，以此类推，直到所有四个角都被单击。然后用红线显示标记的四边形。如果用户单击了白线标记的角度之外的点，将显示一条消息，告诉用户该点的正确位置，用户可以关闭该消息并再次单击。
图像的中点pictureBox1必须在已处理的绘画内部。否则不可能正确点击所有四个角。在这种情况下，必须对图像进行切割，例如，使用 IrfanView 程序，以使整个图像的中点位于画作内部。四个被点击点的坐标保存在数组Point v[4]中。
用户应决定是否删除或部分保留标记方角以外的背景。相应地，他或她应选择“删除背景”或“保留背景”选项以及一个用于指示百分比的选项。
然后需要单击“拉直”。拉直的图像出现在右边的图片框中。用户可以通过单击保存结果来保存结果图像。拉直后的图像可以保存为.bmp或.jpg文件。也可以覆盖现有图像，包括刚刚打开的图像。让我们探索一下这个项目是如何运作的。
拉直的原理
为了拉直图像，我们需要拉直图像的宽度和高度。然后，我们可以通过双线性变换或中心投影将原始图像转换为拉直的图像，如下所述。
实验表明，双线性变换会产生不准确的结果。图像的不同部分被不同地放大，如图 10-3 所示。
图 10-3
通过双线性变换拉直的图像示例
如图 10-3 所示，图像左边的矩形被绘制成宽度，右边的矩形被绘制成高度。当拉直一幅画时，这个缺点并不明显，因为没有矩形，但它仍然是图像的一个不希望的失真。因此，我们决定不使用双线性变换，并开发了一种使用中心变换的方法，尽管这要困难得多。
通过众所周知的中心投影方程，可以找到变换的必要参数。然而，如前所述，我们需要知道从相机的物镜到画的距离，以及指定相机光轴方向的角度。没有办法估计这些参数。相反，我们可以使用有充分根据的假设，即这幅画的框架是一个长方形。我们在获得的图像(图 10-1 )中看到，帧的左边缘的长度不同于右边缘的长度。这些长度的关系以及上边缘和下边缘的长度的关系可以用来估计照相机与画的不同边缘的距离的关系。这些关系然后可以用于估计变换的参数。
考虑图 10-4 中的示意图，该示意图表示通过摄像机物镜 O 的平面以及在获得的图像中帧的左右边缘上的点 M 0 和 M 1 (图 10-5 )。我们将焦平面的翻转副本(M 0 ，C，M 1 )放置在物镜和被摄物体(图中未显示)之间，以便该平面中的图像不会像焦平面中的图像一样反转。
图 10-4
射线与平面的截面(M 0 ，O，M 1
为了找到变换的参数，有必要估计平行于所拍摄的绘画平面的平面的参数，因为所获得的图像在该平面上的中心投影是期望的拉直图像。该平面在图 10-4 中表示为包含该线段的直线(P 0 ，P 1 )。线段(P 0 ，P 1 )是所述平面与平面(M 0 ，O，M 1 )的截面。同样重要的是用户在所获得的图像中的框架的角上点击的四个点。这些点是非矩形四边形的顶点，代表获得的图像中的帧(图 10-5 )。
图 10-5
点 M 0 ，M 1 ，M 2 ，M 3 的说明
通过求解相应的方程可以找到变换的参数。然而，这些方程不是线性的，并且它们很难求解。因此我们用三角学的方法找到参数。
获得的图像大小为宽×高。大多数现代相机的主焦距(即没有变焦的焦距)等于光敏矩阵的宽度。因此我们假设焦距等于宽度和高度的最大值。我们用 f 表示焦距。首先我们计算角度α和β(图 10-4 ):
α = arctan(F/(宽度/2–m0)。(x)
β= arctan(f/(m1)。x–宽度/2)
其中宽度是获得的图像的宽度。
角度φ可以通过正弦定律从三角形(P 0 ，O，C)计算，其中 C 是所获得图像的中心。它的坐标是(宽/2，高/2，F):
(p【0】–o)/sin(π/2+φ)=f/sin(α–φ)
或者
(p【0】–o)=*【f】*sin(π/2+φ)/sin(α-φ)。
类似地，我们从三角形(P 1 ，O，C)得到:
(p【1】–o)=*【f】*sin(π/2–φ)/sin(β+φ)
和
(p【0】–o)= redy(p–o)
其中 RedY 为关系式(P0-O)/(P1-O)，等于得到的图像中帧的右边缘长度与左边缘长度的关系。
因此，我们得到:
f【sin(π/2+φ)/sin(α–φ)= redf【sin(π/2–φ)/sin(β+φ)。
根据等式 sin(π/2 + φ) = cos(φ)和 sin(π/2–φ)= cos(φ)，我们得到:
cos(φ)/sin(α-φ)= redⅲcos(φ)/sin(β+φ)
或减去 cos(φ)后:
sin（β + φ） = RedY⋅ sin（α – φ）。
现在我们把 sin(β + φ)换成 sin β⋅cos φ + cos β⋅sin φ，把 sin(α–φ)换成 sin α⋅cosφ–cos α⋅sinφ:
无β-cosφ+cosβ-无φ= redyⅲ(无α-cosφ–cosαⅲ无φ)。
除以 cos φ后，我们得到:
无β-1+cosβ-TGφ= redyⅲ(无α-1–cosα-TGφ)
或者
TGφⅲ(cosβ+redyⅲcosα)= redyⅲ无α–无β
和
TGφⅲ=(redyⅲ无α-无β)/(cosβ+redyⅲcosα)
或者
φ⋅ = 奥术（RedY⋅sin β – 其α）/（cos α + RedY⋅cos β）。
φ⋅的值允许我们找到包含拉直图像的平面 p 的方程的参数 CX:
(x–宽度/2)* CX+(y–高度/2)* CY+(z–F)= 0。(10.1)
因为直线(P 0 ，P 1 )是 P 与平面(M 0 ，O，M 1 )的交点，所以 CX 近似等于 tan(φ)。这是过程估计，稍后我们会给出精确的估计。
以类似的方式，我们计算平面中的角度φX(M2，O，M 3 ，其中 M 2 和 M 3 是给定图像中具有 M 2 的帧的上边缘和下边缘的点。X = M 3 。X =宽度/2(见图 10-5 )。角度φX 类似于已经计算出的位于平面内的角度φ(M0，O，M 1 )。值 RedX 是框架的上边缘长度与下边缘长度的关系。
φx = arctan(redx⋅sinβx–sinαx)/(cosαx+red⋅cosβx)。
相应的，CY 的估计值为 tan(φX)。
现在让我们计算拉直图像的宽度。我们根据正弦定律从三角形(P 0 ，O，C)获得线段的长度( P 0 ，C):
|p【0】，c | =*【f】*sin(π/2–α)/sin(α–φ)；
而从三角形(P 1 ，O，C)的长度( P 1 ，C):
| > ??【p】【1】，c | =*【f】*sin(π/2–β)/sin(β+φ)。
宽度的计算值等于:
宽度= |p【0】，【p】| = |，【c | >
**f【sin(π/2–α)/sin(α–φ)+sin(π/2–β)/sin(β+φ)】
f(cosα/sin(α-φ)+cosβ/sin(β+φ))。
以类似的方式，我们在考虑平面(M 2 ，O，M 3 时获得高度值，其中 M 2 和 M 3 是给定图像中具有 M 2 的帧的上边缘和下边缘的点。X = M 3 。X =宽度/2(见图 10-5 )。
height =f⋅(cosαx/sin(αx–φx)+cosβx/sin(βx+φx))
其中αX，βX，φX 是平面内的角度(M 2 ，O，M 3 )类似于平面内计算的α，β，φ(M0，O，M 1 )。
我们计算了平面 P 的参数 CX 和 CY 为
CX = tan(φ)；CY = tan(φX)。
实验表明，这些值是近似值。为了找到 CX 和 CY 的精确值，我们开发了方法Optimization，该方法在前面提到的近似值附近测试 CX 和 CY 的许多值。对于 CX 和 CY 的每一对值，该方法计算作为四元组 V 在平面 P 上的中心投影而获得的四元组 Q 与期望的矩形形状的三个偏差。这些偏差如下:
-
偏差“dev1”等于 Q 的水平边到其垂直边的投影。如果 Q 是一个矩形，这个投影应该等于零。
-
偏差“dev2”等于 q 的垂直边的长度差。
-
偏差“dev3”等于 q 的水平边的长度差。
这三个偏差的平方和的平方根是在所有测试对 CX 和 CY 中最小化的标准。有 11 × 11 = 121 对 CX 和 CY 的值被测试，CX 和 CY 的值改变一个步长= 0.08。在步长为 0.01 的 11 × 11 = 121 对变化值中，再次测试与标准最小值相对应的 CX 和 CY 的最佳值。实验表明，根据点击图像角落的用户指定的矢量 v 的精度，该标准的最小值达到 1.0 和 10.0 之间的值。这种优化精度是足够的。甚至通过三角计算获得的 CX 和 CY 的起始值也提供了无法从理想图像中光学区分的拉直图像。
在计算 CX 和 CY 的值之后，方法Rect_Optimal计算四个点 PC[4],这四个点 PC[4]是四个点v[4]在平面 p 上的中心投影。然后，用尺寸Width×Height定义拉直的图像ResultIm。该图像通过两个for循环扫描，坐标为像素的 X 和 Y。对于每个像素(X，Y ),通过双线性变换计算位于具有拐角PC[4]的四元组中并且对应于像素(X，Y)的点的坐标(xp，yp，zp)。坐标为(xp，yp，zp)的点被中心投影到给定图像的平面上。这个投影是一个坐标为(xf，yf，F)的点。给定图像的像素(xf，yf)的颜色被分配给拉直图像的像素(X，Y)。
一种类似的方法Rect_Retain用于产生拉直的图像，该图像包含围绕绘画框架的一部分背景。该方法与Rect_Optimal的不同之处仅在于存在一个将平面 P 中的点PC[4]从坐标为(宽度/2，高度/2，F)的中心 C 移动的过程，使得每个向量的长度(PC[i]，C)，I = 0，1，2，3；乘以参数Rel，用户可以选择值 1.1、1.15 或 1.2。作为这种倍增的结果，拉直的图像包括围绕绘画框架的原始背景的窄条。
最重要方法的代码
下面是最重要的方法Rect_Optimal的源代码。
public void Rect_Optimal(Point[] v, bool CUT, ref CImage Result)
// Calculates the corners of the straightened image and then makes a central projection.
{
bool deb = false;
int MaxSize;
double alphaX, alphaY, betaX, betaY, F, Height, phiX, phiY,
RedX, RedY, Width, M0X, M1X, M3Y, M2Y;
M0X = (double)v[1].X + (v[0].X - v[1].X) *((double)height/2 -
v[1].Y)/((double)v[0].Y-v[1].Y);
M1X = (double)v[2].X + (v[3].X - v[2].X) *((double)height/2 -
v[2].Y)/((double)v[3].Y-v[2].Y);
M3Y = (double)v[0].Y + (v[3].Y - v[0].Y) * ((double)width/2 -
v[0].X)/((double)v[3].X-v[0].X);
M2Y = (double)v[1].Y + (v[2].Y - v[1].Y) * ((double)width/2 -
v[1].X)/((double)v[2].X-v[1].X);
RedY = (double)(v[3].Y - v[2].Y) / (double)(v[0].Y - v[1].Y);
RedX = (double)(v[3].X - v[0].X) / (double)(v[2].X - v[1].X);
if (width > height) MaxSize = width;
else MaxSize = height;
F = 1.0 * (double)(MaxSize);
alphaY = Math.Atan2(F, (double)(width / 2 - M0X));
betaY = Math.Atan2(F, (double)(M1X - width / 2));
phiY=Math.Atan2(RedY*Math.Sin(betaY) -
Math.Sin(alphaY),Math.Cos(alphaY)+RedY*Math.Cos(betaY));
alphaX = Math.Atan2(F, (double)(M3Y - height / 2));
betaX = Math.Atan2(F, (double)(height / 2 - M2Y));
phiX = Math.Atan2(RedX * Math.Sin(betaX) -
Math.Sin(alphaX), Math.Cos(alphaX) + RedX * Math.Cos(betaX));
double P0X = F * Math.Cos(alphaY) / Math.Sin(alphaY - phiY);
double P1X = F * Math.Cos(betaY) / Math.Sin(betaY + phiY);
double P0Y = F * Math.Cos(alphaX) / Math.Sin(alphaX + phiX);
Width = F * (Math.Cos(alphaY) / Math.Sin(alphaY - phiY) +
Math.Cos(betaY) / Math.Sin(betaY + phiY));
Height = F * (Math.Cos(alphaX) / Math.Sin(alphaX - phiX) +
Math.Cos(betaX) / Math.Sin(betaX + phiX));
if (Width < 0.0 || Height < 0.0)
{
MessageBox.Show("The clicked area does not contain the center of the image");
return;
}
double OptCX=0.0;
double OptCY=0.0;
double CX = Math.Tan(phiY);
double CY = Math.Tan(phiX);
Optimization(F, v, CX, CY, ref OptCX, ref OptCY);
CX = OptCX;
CY = OptCY;
CImage Out;
if (CUT)
Out = new CImage((int)Width, (int)Height, N_Bits);
else
Out = new CImage(width, height, N_Bits);
Result.Copy(Out, 0);
double A = 0.0, B, C, D, Det, E, G;
double[] xc = new double[4];
double[] yc = new double[4];
double[] zc = new double[4];
for (int i = 0; i < 4; i++)
{
A = B = C = D = 0.0;
A = (F / (v[i].X - width / 2) + CX);
B = CY;
C = width / 2 * F / (v[i].X - width / 2) + CX * width / 2 + CY * height / 2 + F;
D = CX;
E = (F / (v[i].Y - height / 2) + CY);
G = height / 2 * F / (v[i].Y - height / 2) + CX * width / 2 + CY * height / 2 + F;
Det = A * E - B * D;
xc[i] = (C * E - B * G) / Det;
yc[i] = (A * G - C * D) / Det;
zc[i] = F - CX * (xc[i] - width / 2) - CY * (yc[i] - height / 2);
}
double zz;
double xp, yp, xp0, xp1, yp0, yp1, xf, yf;
for (int Y = 0; Y < Result.height; Y++) //= over the straightened image
{
xp0 = xc[1] + (xc[0] - xc[1]) * Y / (Result.height - 1);
xp1 = xc[2] + (xc[3] - xc[2]) * Y / (Result.height - 1);
for (int X = 0; X < Result.width; X++) //=============================
{
yp0 = yc[1] + (yc[2] - yc[1]) * X / (Result.width - 1);
yp1 = yc[0] + (yc[3] - yc[0]) * X / (Result.width - 1);
xp = xp0 + (xp1 - xp0) * X / (Result.width - 1);
yp = yp0 + (yp1 - yp0) * Y / (Result.height - 1);
zz = F - CX * (xp - width / 2) - CY * (yp - height / 2); // corrected
xf = width / 2 + (xp - width / 2) * F / (F - CX * (xp - width / 2) -
CY * (yp - height / 2));
yf = height / 2 + (yp - height / 2) * F / (F - CX * (xp - width / 2) -
CY * (yp - height / 2));
if ((int)xp >= 0 && (int)xp < width && (int)yp >= 0 && (int)yp < height)
if (N_Bits == 24)
{
for (int ic = 0; ic < 3; ic++)
Result.Grid[ic+3*X+3*Result.width*Y]=Grid[ic+3*(int)xf+3*width*(int)yf];
}
else
Result.Grid[X+Result.width*(Result.height-1-Y)]=Grid[(int)xf+width*(int)yf];
} //================= end for (X... ==============================
} //================== end for (Y... ==============================
} //********************* end Rect_Optimal ******************************
这里是Optimization的代码。
private void Optimization(double F, Point[] v, double CX, double CY,
ref double OptCX, ref double OptCY)
{
bool deb = false;
double A = 0.0, B, C, D, Det, E, G;
double[] xc = new double[4];
double[] yc = new double[4];
double[] zc = new double[4];
double[] xopt = new double[4];
double[] yopt = new double[4];
double[] zopt = new double[4];
double dev1, dev2, dev3, Crit;
double MinCrit = 10000000.0;
int OptIterX = -1, OptIterY = -1, IterY;
OptCX = 0.0;
OptCY = 0.0;
CX -= 0.40; CY -= 0.40;
double CX0 = CX, Step = 0.08;
for (IterY = 0; IterY < 11; IterY++)
{
CX = CX0;
for (int IterX = 0; IterX < 11; IterX++)
{
for (int i = 0; i < 4; i++)
{
A = (F / (v[i].X - width / 2) + CX);
B = CY;
C = width / 2 * F / (v[i].X - width / 2) + CX * width / 2 + CY * height / 2 + F;
D = CX;
E = (F / (v[i].Y - height / 2) + CY);
G = height / 2 * F / (v[i].Y - height / 2) + CX * width / 2 + CY * height / 2 + F;
Det = A * E - B * D;
xc[i] = (C * E - B * G) / Det;
yc[i] = (A * G - C * D) / Det;
zc[i] = F - CX * (xc[i] - width / 2) - CY * (yc[i] - height / 2);
}
dev1 = dev2 = dev3 = 0.0;
dev1 = ((xc[0] - xc[1]) * (xc[2] - xc[1]) + (yc[0] - yc[1]) * (yc[2] - yc[1]) + (zc[0] - zc[1]) * (zc[2] - zc[1])) /
Math.Sqrt(Math.Pow((xc[0] - xc[1]), 2.0) + Math.Pow((yc[0] - yc[1]), 2.0) + Math.Pow((zc[0] - zc[1]), 2.0));
dev2 =
Math.Sqrt(Math.Pow((xc[3] - xc[2]), 2.0) + Math.Pow((yc[3] - yc[2]), 2.0) + Math.Pow((zc[3] - zc[2]), 2.0)) -
Math.Sqrt(Math.Pow((xc[0] - xc[1]), 2.0) + Math.Pow((yc[0] - yc[1]), 2.0) + Math.Pow((zc[0] - zc[1]), 2.0));
dev3 =
Math.Sqrt(Math.Pow((xc[2] - xc[1]), 2.0) + Math.Pow((yc[2] - yc[1]), 2.0) +
Math.Pow((zc[2] - zc[1]), 2.0)) -
Math.Sqrt(Math.Pow((xc[3] - xc[0]), 2.0) + Math.Pow((yc[3] - yc[0]), 2.0) + Math.Pow((zc[3] - zc[0]), 2.0));
Crit=Math.Sqrt(Math.Pow(dev1, 2.0)+Math.Pow(dev2, 2.0)+Math.Pow(dev3, 2.0));
if (Crit < MinCrit)
{
MinCrit = Crit;
OptIterX = IterX;
OptIterY = IterY;
OptCX = CX;
OptCY = CY;
for (int i = 0; i < 4; i++)
{
xopt[i] = xc[i];
yopt[i] = yc[i];
zopt[i] = zc[i];
}
}
CX += Step;
} //============= end for (int IterX... ========================
CY += Step;
} //============== end for (int IterY... ========================
CX = OptCX; CY = OptCY;
CX -= 0.05; CY -= 0.05;
CX0 = CX;
double step=0.01;
for (IterY = 0; IterY < 11; IterY++)
{
CX = CX0;
for (int IterX = 0; IterX < 11; IterX++)
{
for (int i = 0; i < 4; i++)
{
A = (F / (v[i].X - width / 2) + CX);
B = CY;
C = width / 2 * F / (v[i].X - width / 2) + CX * width / 2 + CY * height / 2 + F;
D = CX;
E = (F / (v[i].Y - height / 2) + CY);
G = height / 2 * F / (v[i].Y - height / 2) + CX * width / 2 + CY * height / 2 + F;
Det = A * E - B * D;
xc[i] = (C * E - B * G) / Det;
yc[i] = (A * G - C * D) / Det;
zc[i] = F - CX * (xc[i] - width / 2) - CY * (yc[i] - height / 2);
}
dev1 = dev2 = dev3 = 0.0;
// deviation from a 90° angle:
dev1 = ((xc[0] - xc[1]) * (xc[2] - xc[1]) + (yc[0] - yc[1]) * (yc[2] - yc[1]) + (zc[0] - zc[1]) * (zc[2] - zc[1])) /
Math.Sqrt(Math.Pow((xc[0] - xc[1]), 2.0) + Math.Pow((yc[0] - yc[1]), 2.0) + Math.Pow((zc[0] - zc[1]), 2.0));
dev2 = // difference |pc[3] - pc[2]| - |pc[0] - pc[1]|:
Math.Sqrt(Math.Pow((xc[3] - xc[2]), 2.0) + Math.Pow((yc[3] - yc[2]), 2.0) + Math.Pow((zc[3] - zc[2]), 2.0)) -
Math.Sqrt(Math.Pow((xc[0] - xc[1]), 2.0) + Math.Pow((yc[0] - yc[1]), 2.0) + Math.Pow((zc[0] - zc[1]), 2.0));
dev3 = // difference |pc[2] - pc[1]| - |pc[3] - pc[0]|:
Math.Sqrt(Math.Pow((xc[2] - xc[1]), 2.0) + Math.Pow((yc[2] - yc[1]), 2.0) + Math.Pow((zc[2] - zc[1]), 2.0)) -
Math.Sqrt(Math.Pow((xc[3] - xc[0]), 2.0) + Math.Pow((yc[3] - yc[0]), 2.0) + Math.Pow((zc[3] - zc[0]), 2.0));
Crit =Math.Sqrt(Math.Pow(dev1,2.0)+Math.Pow(dev2,2.0)+Math.Pow(dev3,2.0));
if (Crit < MinCrit)
{
MinCrit = Crit;
OptIterX = IterX;
OptIterY = IterY;
OptCX = CX;
OptCY = CY;
for (int i = 0; i < 4; i++)
{
xopt[i] = xc[i];
yopt[i] = yc[i];
zopt[i] = zc[i];
}
}
CX += step;
} //============== end for (int IterX... ========================
CY += step;
} //=============== end for (int IterY... ========================
} //****************** end Optimization *******************************
结论
这个项目带来了良好的效果，并且易于使用:它为用户提供了图形支持，并在用户出错时给出提示。**
十一、区域边界和边的多边形近似
本章描述了一种将二维数字图像中的曲线表示为多边形的方法。这种曲线表示对于图像分析是有用的，因为多边形的形状可以通过简单的几何方法(例如测量长度和角度)来容易地研究。多边形近似也提出了一种估计数字曲线曲率的新方法。为此，多边形可以用圆弧和直线段的平滑序列来代替。平滑是指每条直线段都是前一条和后一条圆弧的切线。
本章包含相关出版物的简短回顾、已知多边形近似方法的简短描述、Schlesinger (2000)提出的数字曲线相似性新度量的定义以及使用该度量的新近似方法。文中还介绍了新的曲率估计方法的原理和一些实验结果。
多边形逼近问题
给出的是二维图像中的数字曲线 C 。我们希望用一个边尽可能少的多边形来近似 C ，并且与 C 相似。相似性的一个可能标准是 Hausdorff 距离:设S1 和S2 为两组点。用 d ( p 、 q )表示点 p 和 q 之间的距离。那么下列定义成立:
-
D ( p ，s)= mind(p，q)∀q∈s是从点 p 到设定 S 的距离。
-
D ( S 1 ，s2)= maxd(p，s2)∀p∈s1是距离设定值 S 的距离
-
D ( S 2 ，s1)= maxd(q，s1)∀q∈s2是距离设定值 S 的距离
-
H ( S 1 ，S2)= max(D(S1， S 2 )，D(S2， S
然而，Hausdorff 距离并不总是两条曲线相似性的合适度量。图 11-1 中红色和黑色多边形之间的 Hausdorff 距离小于 2 个像素；然而，真实的偏差要大得多:多边形并不相似。Schlesinger (2000)提出了一种更合适的曲线相似性度量方法。
图 11-1
Hausdorff 距离小于 2 的两个多边形的示例；所有箭头的长度都是 2
施莱辛格曲线相似性度量
给出了两条数字曲线C?? 1 和C2。沿C1 取 m 等间距点，沿C2 取 n 点。这样一来就得到两个点了套M1 和M2 的分。沿着曲线的后续点之间的距离越小，相似性的估计就越精确。让我们将两点 p 和 q 之间的距离记为 d ( p ， q )。
求对应点间最大距离最小值的单调映射f:m1➔m2。一个有序集 M 1 到另一个有序集 M 2 上的映射称为单调，如果对于任意两点P1<P2的 M 1 的条件 F (对应点之间的最大距离的最小值就是曲线C1 和C2 之间的施莱辛格距离 DS 。
DS ( M 1 ， M 2 ) =最大值(d【p】
Schlesinger (2000)还提出了一种计算任意两条数字曲线的距离 DS 的有效算法。在我们的例子中， DS 的值可以计算如下。设黑色曲线为 B ，红色曲线为 R 。曲线被细分为相应的线段( a 、 b )、 b 、 d )、 d 、 f 、 f 、 a )。(点 a 、 b 、 *d、*和 f 位于两个多边形上。)如图 11-2 所示，将 B 段的每一点映射到 R 对应段的最近点上。
图 11-2
施莱辛格相似性的定义
图 11-2 示例中 DS 的值由( c，d )(或( e，d )的长度定义，等于 3.6。它本质上大于等于 2 的 Hausdorff 距离，并且更精确。
近似问题的陈述
给定一条数字曲线 C 和一个公差值ε，找出边数最少的折线 P ，使得施莱辛格距离 DS ( C ， P )不大于ε。
在下一节中，我们将考虑解决这个问题的各种方法。
多边形逼近算法
我们在这里介绍一些从文献中已知的算法和一些新的算法。
拆分-合并方法
这种方法的优点是非常容易理解和编程。给定一条封闭曲线 C 和从 C 点到多边形的最大允许距离ε。在 C 上任意取一点V1。找到与V1 距离最大的点V??。现在曲线细分为两段:( V 1 ， V 2 )和( V 2 ， V 1 )。每个线段( V i ， V ， k )都有它的弦(在算法开始时两个线段有相同的弦)。对于每一段，找到与相应弦距离最大的点 V m 。如果 D > ε，将线段( V i ， V k )分割成两段( V i ， V 重复此操作，直到所有分段的 D ≤ ε。分相结束了。
现在检查每对相邻线段( V i ， V m )和( V m ， V k )的距离DV*如果是，合并段( V i ， V m )和( V m ， V k )。当所有对满足条件Dε时，停止合并阶段和算法。ε = 1.5 的例子如图 11-3 所示。*
*
图 11-3
拆分和合并方法的示例
新点的顺序是 a，b，c，d，e。线段(e，b)和(b，d)已经合并。给定曲线由粉色区域周围的黑色步进线表示。起点为 a。距离 a 最远的点为 b。距离弦(a，b)最远的上段点为 c，该点与 a 和 b 相连。距离弦(a，b)最远的下段点为 D。该点也与 a 和 b 相连。距离弦(c，b)最远的点为 e。该点与 c 和 b 相连。在此步骤之后，所有线段都满足 D ≤ ε且ε = 1.5 的条件。分相结束了。对于相邻线段对(e，b)和(b，d ),它们的公共顶点 b 到弦(e，d)的距离小于ε = 1.5。因此，这些部分应该合并。最终的近似多边形具有顶点 a，c，e，d。
扇形法
拆分-合并方法相当慢，因为每个点都必须测试多次，这取决于它所属段的拆分次数。这种方法使豪斯多夫距离最小，而不是施莱辛格距离。
威廉姆斯(1978)提出了一种有效的方法，称为扇形法。从曲线 C 的某点 V 1 (图 11-4 )开始。测试所有后续点。对于每个点 V i 用 i > 1 和D(V1、VI)>ε，画一个半径等于公差ε的圆，从 V 画两条切线
图 11-4
扇形法
切线之间的空间是一个扇形，对于位于扇形内的每条直线 L 到V1，圆心到 L 的距离小于公差ε。如果下一个点VI(V3或 V 4 在图 11-4 中)位于扇形内，则圆心在 V i 和扇形 S 的一个新圆新扇区是旧扇区与新扇区的交集。如果下一个点 V i (例如图 11-4 中的 V 5 位于新扇区之外，则没有直线 L 通过VV1VI 因此，线段( V 1 ， V i )不能作为多边形边。 V i 之前的点VI-1则是当前线段的最后一个点；[ V 1 ，VI-1]是多边形的边。下一条多边形边的构造从点VI-1开始。
扇形法的改进
扇形法很快，因为它只测试给定曲线的每个点一次。但是，它仅保证曲线和多边形之间的 Hausdorff 距离小于规定的公差ε。引入以下附加条件是为了中断应由多边形的当前边近似的曲线的当前段。当前段的点 V max 距离该段的起点最远，必须在该段的加工过程中定义。如果实际点 V i 投影到直线段上的距离 V 1 ，Vmax小于| V 1 ，Vmax|-ε这个测试保证了曲线的点在多边形上的投影顺序是单调的，因此满足了定义施莱辛格距离的条件。
在图 11-5 所示的例子中，近似多边形有顶点 a 、 b 、 c 、 *d、*和 e 。红点是算法决定中断当前段的点。公差ε = 1.5。
图 11-5
扇形法的改进
本章后面介绍的项目WFpolyArc使用了改进的扇形法。
用圆弧和直线序列替换多边形
定义和问题陈述
定义 11.1: 如果满足以下条件，一系列圆弧和直线段称为平滑:
-
每条直线段 L 跟随一条圆弧Ap并且跟随另一条圆弧 A f 。
-
线段 L 与Ap和 A f 在其公共端点相切。
-
如果一条弧后面跟着另一条弧，那么它们在公共端点处有一条公共切线。
定义 11.2: 闭合折线 P 的 ε公差带是另外两条折线P1 和P2 之间的区域，使得PI， i = 1，2，平行于的相应边闭合折线P1 内部包含 P 且 P 内部包含P2。
我们考虑以下问题:给定已知以容差ε逼近数字曲线 C 的多边形 P ，找出圆弧和直线段的平滑序列 S ，使得其以容差ε逼近的多边形等于 P ，并且圆弧曲率的绝对值之和最小。
我们假设序列 S 的弧的曲率是曲线 C 的曲率的良好估计，因为 S 和 C 之间的距离小于ε。圆弧曲率的绝对值之和必须最小，因为否则可以找到具有更高曲率值的小圆弧的另一平滑序列 S ，其中到 C 的距离小于 S 和 C 之间的距离，而 S 的圆弧曲率值可以任意高，并且与曲率无关
这个问题的确切解决方案还不知道。我们建议采用以下近似解决方案。
近似解
定义 11.3: 公差带 T 边界的一个顶点 W 若与 W 相交的 T 的边之间的钝角位于 T 之外，则称为凹入。
在下文中，我们描述了计算曲率值的过程。该程序为近似多边形的每个顶点 V (图 11-6 )计算与关联于 V 的多边形边相切的两个圆弧的曲率C1 和C2 的两个候选值。
首先考虑弧 A 穿过顶点 W 并与多边形边相切。如果圆弧 A 相切的点T1 和T2 的D1 和D2 距离 V 小于或等于半长S 然后从三角形OP1W(图 11-6 )计算圆弧半径如下:
(p【1】或=【r】【r】1】
或者
r【1】–【cosβ=ε
*或者
r1=ε/(1–cosβ)(11.1)
其中β是多边形边之间的角度γ的一半。
图 11-6
曲率值的计算
然而，如果距离 D 1 和 D 2 大于从 V 到与 V 相关的两条多边形边中较短的一条边的中点 M 的距离，则不可能使用该圆弧，因为它可能与下一个多边形顶点处的圆弧重叠。在这种情况下，有必要构建另一个不经过点 W 的圆弧，并在其中点 M 处与较短的边相切，在另一边上与 M 对称的点处与另一边相切(图 11-7 )。
图 11-7
曲率值的计算
该圆弧的半径 R 2 应从三角形 OMV 开始计算(图 11-7 ):
R2= ctgβSmin/2；(11.2)
所述计算通过接下来描述的方法Curvature进行。该方法不计算半径，而是计算作为半径倒数的曲率值:
C1= 1/R1；以及C2= 1/R2。(11.3)
该方法同时计算 C 1 和 C 2 。很容易看出，如果距离 D 1 和 D 2 比短边的半长短(如图 11-6 ，那么R1<R2， C 1 >在这种情况下，曲率的结果值 C 被设置为 C 1 并且所寻找的圆弧的两个切点到顶点 V 的距离D=R1tanβ被设置为D=R1。
在图 11-7 所示的情况下，半径R2<R1曲率的结果值 C 被设置为 C 2 并且被求圆弧的两个切点到顶点 V 的距离 D 被设置为 D
计算曲率后，该方法会计算弧的端点，这些端点等于切点。为此，在 V 之前的多边形边的长度 S 1 被细分为比例a1:(1-a1)与a1=D/S1。 V 后沿的长度 S 2 细分为比例a2:(1-a2)与a2=D/S2。
最后，该方法测试弧的长度是否为零。如果一个小的闭合多边形只包含三个顶点，就会发生这种情况。在这种情况下，程序会返回一条错误消息，并且不会保存 arc。下面是方法Curvature的代码。
public int Curvature(int x1, int y1, int x2, int y2, int x3, int y3,
double eps, ref CArc arc)
/* Calculates the curvature 'k' of an arc lying in the tolerance tube around the given two polygon edges [(x1,y1), (x2,y2)] and [(x2,y2), (x3,y3)] and having the outer boundary of the tube as its tangent. The tangency point should not be farther than the half-length of the shorter edge from (x2,y2). The radius of the arc should be as large as possible. */
{
double a1, a2, lp, // Variables for calculating the projections
dx1, dy1, dx2, dy2, len1, len2, // Increments, lengths of edges
cosgam, // cosine of the angle between the edges
sinbet, cosbet, // cosine and sine of the half-angle between the edges
strip = 0.6, // correcture of the deviation
k, cru1, cru2; // curvature for long and short edges
dx1 = (double)(x2 - x1); dy1 = (double)(y2 - y1); // first edge
len1 = Math.Sqrt(dx1 * dx1 + dy1 * dy1);
dx2 = (double)(x3 - x2); dy2 = (double)(y3 - y2); // second edge
len2 = Math.Sqrt(dx2 * dx2 + dy2 * dy2);
if ((len1 == 0.0) || (len2 == 0.0)) return (-1);
cosgam = (dx1 * dx2 + dy1 * dy2) / len1 / len2; // cosine of angle between the edges
if (Math.Abs(cosgam) <= 1.0)
sinbet = Math.Sqrt((1.0 - cosgam) * 0.5); // sine of half-angle
else sinbet = 0.0;
cosbet = Math.Sqrt((1.0 + cosgam) * 0.5);
cru1 = (1.0 - cosbet) / (eps - strip); // long edges, it is important
// that the arc goes throught the vertex
double min_len = len1;
if (len2 < min_len) min_len = len2;
if (min_len != 0.0 && cosbet != 0.0)
cru2 = 2.0 * sinbet / cosbet / min_len; // short edges, it is important
else cru2 = 0.0; // that the tangency is in the middle of the shortest edge
if ((Math.Abs(cru1) > Math.Abs(cru2)) && cru1 != 0.0)
{
if (cosbet != 0.0 && cru1 != 0.0)
lp = sinbet / cosbet / cru1; // distance of the point of tangency from (x2, y2)
else lp = 100.0;
k = cru1; // first curvature
}
else
{
lp = min_len / 2.0;
k = cru2 * 0.95; // second curvature
}
if (dx1 * dy2 - dy1 * dx2 > 0.0) k = -k; // the sign
// the first edge is divided in relation a1 : a2
a1 = lp / len1;
a2 = 1.0 - a1;
if (k != 0.0) arc.rad = 1.0F / (float)k;
else arc.rad = 0.0F;
arc.xb = (float)(x1 * a1 + x2 * a2); // first tangency (begin of the arc)
arc.yb = (float)(y1 * a1 + y2 * a2);
a1 = lp / len2; a2 = 1.0 - a1;
arc.xe = (float)(x3 * a1 + x2 * a2); // second tangency (end of the arc)
arc.ye = (float)(y3 * a1 + y2 * a2);
double AbsR, R, chord_X, chord_Y, Length, Lot_x, Lot_y;
R = arc.rad;
AbsR = R;
if (R < 0.0) AbsR = -R;
chord_X = arc.xe - arc.xb; chord_Y = arc.ye - arc.yb; // the chord
Length = Math.Sqrt(chord_X * chord_X + chord_Y * chord_Y); // length of chord
if (R < 0.0) // 'Lot' is orthogonal to chord
{
Lot_x = chord_Y; Lot_y = -chord_X;
}
else
{
Lot_x = -chord_Y; Lot_y = chord_X;
}
if (2 * AbsR < Length) return -1;
if (Length < 0.1) return -2;
double Lot = Math.Sqrt(4 * R * R - Length * Length);
arc.xm = (float)((arc.xb + arc.xe) / 2 - Lot_x * Lot / 2 / Length);
arc.ym = (float)((arc.yb + arc.ye) / 2 - Lot_y * Lot / 2 / Length);
if ((arc.xe == arc.xb) && (arc.ye == arc.yb))
{ // can be a closed line of three points ??
return -2;
}
return 0;
} // ******************* end Curvature **********************************
项目WFpolyArc
本项目的形式如图 11-8 所示。
图 11-8
项目的形式WFpolyArc
正如我们的其他项目一样，该表单包含一个打开图像按钮，使用户能够打开图像。
下一个按钮是检测边缘。点击此按钮启动第七章所述的边缘检测程序。单击多边形启动边的多边形近似过程。近似精度由近似精度设置中显示的值定义。显示值预设最大近似误差，以像素为单位。多边形显示在右边的图片框中。测试每个多边形以确定它是否平滑。这对弧的定义很重要:只为平滑多边形确定弧是有利的。平滑度的定义如下。子例程为多边形的每个顶点定义该顶点处两条边之间的角度符号是否不同于多边形的第一和第二条边之间的角度符号。角度值与平滑多边形允许的最大角度进行比较。如果角度过大或符号错误的节点的百分比超过预定阈值，则认为多边形不平滑。
单击弧启动定义平滑多边形弧的方法。绘制从弧的端点开始到曲率中心结束的半径向量。半径向量的颜色表示曲率的符号:黄色半径表示正曲率的圆弧，紫色表示负曲率的圆弧。多边形和圆弧的参数可以保存在一个文本文件中，如后面所述。用户可以通过单击保存来保存多边形和圆弧列表。当读取文本文件时，可以看到参数值(例如，使用 Microsoft Word)。
项目中使用的方法WFpolyArc
用于边缘检测的方法在章节 6 和 7 中介绍。
当在细胞复合体中使用检测到的边缘的表示时，实现边缘的多边形近似的方法可以变得更简单和更容易理解(Kovalevsky，2008)。第七章描述了细胞复合体的概念。
为了使检测到的边缘的处理更简单和更容易理解，我们使用图像CombIm，其中边缘被表示为裂纹序列。通过方法LabelCellsNew填充图像CombIm，如项目WFcompressPal(第章第节)中一样。该方法读取由方法ExtremVarColor或灰度图像的类似方法产生的图像EdgeIm。它寻找颜色差异超过预定阈值的相邻像素对。对于每一对这样的像素，位于它们之间的裂缝由1标记。此外，在这种裂纹末端的点(0-单元)被标记，并且该标记等于属于该边缘并且与该点关联的裂纹的数量。
下面是方法LabelCellsNew的代码。
public void LabelCellsNew(int th, CImage Image3)
/* Looks in "Image3" for all pairs of adjacent pixels with color differences greater than "th" and labels the corresponding cracks in "this" with "1". The points incident with the crack are also provided with labels indicating the number and locations of incident cracks. This method works both for color and gray value images. ---*/
{ int difH, difV, nComp, NXB=Image3.width, x, y;
byte Lab = 1;
if (Image3.N_Bits==24) nComp=3;
else nComp=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 the 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<nComp; c++)
{ Colorv[c]=Image3.Grid[c+nComp*((x-2)/2)+nComp*NXB*(y/2)];
Colorp[c]=Image3.Grid[c+nComp*(x/2)+nComp*NXB*(y/2)];
}
if (nComp==3) difV=ColorDif(Colorp,Colorv);
else difV=(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<nComp; c++)
{ Colorh[c]=Image3.Grid[c+nComp*(x/2)+nComp*NXB*((y-2)/2)];
Colorp[c]=Image3.Grid[c+nComp*(x/2)+nComp*NXB*(y/2)];
}
if (nComp==3) difH=ColorDif(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 LabelCellsNew ****************************
边缘的多边形近似从方法SearchPoly开始。该方法扫描图像中的点CombIm。当它找到一个标签等于 1、3 或 4 的点时，它调用方法ComponPoly。该方法在类CListLines的对象中编码包含由SearchPoly找到的起始点(X，Y)的边组件的多边形。该方法使用先进先出队列。该方法将起始点Pinp放入队列并开始一个while循环。它测试从队列中取出的与点 P 相关的每个标记裂纹。如果裂纹的下一个点是分支或终点，则裂纹被忽略。否则调用函数TraceApp。TraceApp跟踪边，直到下一个端点或分支点，并计算近似多边形的顶点。追踪在一个叫做Pterm的点结束。如果点Pterm是分支点，则将其放入队列。ComponPoly队列为空时返回。
这里是ComponPoly的代码。
public int ComponPoly(CImage Comb, int X, int Y, double eps)
{ int dir, dirT;
int LabNext, rv;
iVect2 Crack, P, Pinp, Pnext, Pterm;
Crack = new iVect2();
P = new iVect2();
Pinp = new iVect2();
Pnext = new iVect2();
Pinp.X = X;
Pinp.Y = Y; // comb. coord.
int CNX = Comb.width;
int CNY = Comb.height;
pQ.Put(Pinp);
while(!pQ.Empty()) //============================================
{
P = pQ.Get();
if ((Comb.Grid[P.X + CNX * P.Y] & 128) != 0) continue;
for (dir=0; dir<4; dir++) //======================================
{
Crack.X = P.X + Step[dir].X;
Crack.Y = P.Y + Step[dir].Y;
if (Crack.X < 0 || Crack.X > CNX - 1 || Crack.Y < 0 ||
Crack.Y > CNY - 1) continue;
if (Comb.Grid[Crack.X+CNX*Crack.Y]==1) //---- ------------ --------
{
Pnext.X = Crack.X + Step[dir].X;
Pnext.Y = Crack.Y + Step[dir].Y;
LabNext = Comb.Grid[Pnext.X + CNX * Pnext.Y] & 7;
if (LabNext==3) pQ.Put(Pnext);
if (LabNext==2) //-------------------------------------------------
{
Polygon[nPolygon].firstVert = nVert;
dirT = dir;
Pterm = new iVect2();
rv=TraceApp(Comb, P.X, P.Y, eps, ref Pterm, ref dirT);
if (rv<0)
{ MessageBox.Show("ComponPoly, Alarm! TraceApp returned " + rv);
return -1;
}
if (nPolygon>MaxPoly-1)
{ MessageBox.Show("ComponPoly: Overflow in Polygon; nPolygon=" +
nPolygon + " MaxPoly=" + MaxPoly);
return -1;
}
else nPolygon++;
if ((Comb.Grid[Pterm.X+CNX*Pterm.Y] & 128)==0 && rv>=3)
pQ.Put(Pterm);
} // ------------- end if (LabNest==2) ----------------------------
if ((Comb.Grid[P.X+CNX*P.Y] & 7)==1) break;
} //--------------- end if (Comb.Grid[Crack.X ...==1) --------------
} //========= end for (dir ... =====================================
Comb.Grid[P.X+CNX*P.Y] |=128;
} //========== end while =======================================
return 1;
} //************* end ComponPoly **************************************
当计算出近似多边形时，用户可以开始计算圆弧。方法FindArcs读取多边形列表，调查每个多边形的三个后续顶点，并通过已经描述的方法Curvature计算对应于三元组的中间顶点的弧的曲率和端点。这些计算的原理在本章前面已经描述过了。方法FindArcs产生弧的列表。列表中的每个条目都包含端点的坐标、曲率中心的坐标以及弧的曲率半径值，所有这些都以像素为单位。这里是FindArcs的代码。
public int FindArcs(PictureBox pictureBox, CImage EdgeIm, double eps, Form1 fm1)
/* The method calculates the parameters of the arcs contained in the polygons. Fills the array 'Arc []' of structures "CArc". Shows the contents of this array graphically. */
{
bool deb = false, disp = true;
int j, ip, first, last, Len, rv,
x1, y1, x2, y2, x3, y3, // three polygon vertices composing an arc
xb, yb, xe, ye, old_xe = 0, old_ye = 0; // "int" boundaries of arcs
nArc = 0;
Pen linePen = new System.Drawing.Pen(Color.LightBlue);
int marginX = fm1.marginX;
int marginY = fm1.marginY;
double Scale1 = fm1.Scale1;
for (ip = 0; ip < nPolygon; ip++) // ===========================
{
int cntArc = 0;
Polygon[ip].firstArc = -1;
Polygon[ip].lastArc = -2;
if (!Polygon[ip].smooth) continue;
first = Polygon[ip].firstVert;
last = Polygon[ip].lastVert;
old_xe = -1;
old_ye = -1;
if (last > first + 1) // ------- are there sufficient vertices?--------
{
if (disp) // here are three points of a polygon calculated but not drawn
{
x1 = (int)(Scale1 * Vert[first].X + 0.5); // starting point
y1 = (int)(Scale1 * Vert[first].Y + 0.5);
}
x2 = Vert[first].X; // will become the first point of the arc
y2 = Vert[first].Y;
x3 = Vert[first + 1].X; // second point
y3 = Vert[first + 1].Y;
Polygon[ip].firstArc = nArc;
// Points from the second one until the one before the last
Len = last - first + 1;
for (j = 2; j <= Len; j++) // =============================
{
x1 = x2; y1 = y2;
x2 = x3; y2 = y3;
x3 = Vert[first + j % Len].X;
y3 = Vert[first + j % Len].Y;
CArc arc = new CArc();
// 'Curvature' calculates and saves parameters of an arc in 'arc'.
rv = Curvature(x1, y1, x2, y2, x3, y3, eps, ref arc);
if (rv < 0) continue;
if (Math.Abs(arc.xb - arc.xe) < 1.0 &&
Math.Abs(arc.yb - arc.ye) < 1.0 || Math.Abs(arc.rad) < 2.0) continue;
// The arc is saved in the list of arcs:
Arc[nArc] = arc;
if (cntArc == 0) Polygon[ip].firstArc = nArc;
cntArc++;
nArc++;
if (disp) //-------------------------------------------------------
{
xb = marginX + (int)(Scale1 * arc.xb + 0.5);
yb = marginY + (int)(Scale1 * arc.yb + 0.5);
if (j == Len)
{
xb = marginX + (int)(Scale1 * Vert[last].X + 0.5);
yb = marginY + (int)(Scale1 * Vert[last].Y + 0.5);
}
xe = marginX + (int)(Scale1 * arc.xe + 0.5);
ye = marginY + (int)(Scale1 * arc.ye + 0.5);
rv = DrawArc(arc, true, fm1);
if (rv < 0) return -1;
Pen bluePen = new System.Drawing.Pen(Color.LightBlue);
if (j > 2 && j < Len && last - first > 2 && old_xe > 0)
fm1.g2Bmp.DrawLine(bluePen, old_xe, old_ye, xb, yb);
old_xe = xe; old_ye = ye;
} //---------------------- end if (disp) --------------------------
if (j == Len - 1) Polygon[ip].lastArc = nArc - 1;
} // ============== end for (j... ================================
} // ------------------------- end if (last > first+1 ) ---------------
} // ================ end for (ip... ================================
return 0;
} // ******************* end FindArcs ***********************************
点击保存，可将多边形和圆弧列表保存在硬盘上，生成一个扩展名为.txt的文本文件，其中包含多边形列表和圆弧列表。它包含每个多边形的参数，如firstVertex、lastVertex、firstArc和lastArc。多边形列表之后是圆弧列表。对于每个圆弧，它包含起点的坐标(xb，yb)、终点的坐标、曲率中心的坐标和半径。该文件可以通过任何文本处理程序(例如，Microsoft Word)来读取。目录TEXT和文本文件的名称根据打开的图像名称自动计算。TEXT目录应该位于包含打开图像的目录附近。在第十二章中描述的项目WFcircleReco中可以看到一个版本的保存，允许用户选择目录和文本文件的名称。第十二章中也描述了保存圆圈。
半径计算的精度
我们使用一个具有 328 和 95 像素半轴的理想椭圆的人造图像研究了所述方法的精度(图 11-9 )。
图 11-9
实验中使用的椭圆的人工图像
因为我们知道椭圆的参数，所以可以计算边界上每一点的曲率半径。计算值与FindArcs产生的值的比较表明，在真实曲率半径小于 174 像素的边界点，产生值的误差小于 15%。在更大的半径下，误差小于 20%。这种精度对于某些应用来说已经足够了。
为了测试小半径情况下的误差值，我们使用了另一个半径为 100 和 200 像素的人造图像(图 11-10 )。在两个人造图像的测试中，曲率半径的估计误差在 20%以下。这与理论研究(Kovalevsky，2001)一致，理论研究表明，数字图像中曲线导数的估计误差由数字化噪声的强度指定，使得像素坐标相当不精确。
图 11-10
实验中使用的自行车链条链节的数码照片(孔被人工填充为黑色，以使图像更简单)
这幅图像的边界包含一个膨胀点，在这里曲率改变了它的符号。该点附近的平均曲率接近于零，这对应于相当大的半径值。在其他点上，估计曲率半径的误差在 20%以下。
结论
本文报道的结果表明，该方法的精度是令人满意的。**
十二、圆形物体的识别和测量
本章介绍了一种识别和测量边界形状类似于圆形的物体参数的方法。我们使用我几年前开发的一种改进的最小二乘法来检查晶片中焊料凸点的质量。在这里，我们描述了这种方法的发展，它可以应用于任何彩色或灰度图像，以识别具有近似圆形的边界的对象。该方法计算半径和中心坐标的估计值。图 12-1 显示了根据苹果大小对其进行分类的示例。
图 12-1
(a)图像和(b)结果的示例
该方法的数学基础
该方法类似于经典的最小二乘法，其本质区别在于，我们最小化的不是给定点与所寻找圆的偏差的平方和，而是半径未知的圆 R 和中心未知的圆( x c 的面积差的平方和， y c 通过给定点(xI，y i ):
deviat =[r-x【I】-(12.1)
*根据积分的指数 i 计算总和，其中 i 得到从 0 到 N - 1 的 N 值。这个标准不同于传统的标准
crit =[r-x【I】-【c】
*我们使用半径差 2 -距离 2 来代替半径差(一个点到中点的距离)。这使得求解变得更加容易和快速，并且得到的圆与经典圆之间的差异非常小。
让我们推导出必要的方程式。我们取方程 12.1 在 x c ， y c ， R 之后的偏导数。
★deviat/x【c】=【4 (**【r】-)
★deviat/和【c】=【4 (**【r】-)
★deviat/【r =【4 ，【r】-【我】 -】***
并将它们设置为 0。我们得到三个方程:
r-x【I】-
r-x【I】-
r-x【I】-
**左边部分是包含变量xcT5、 y c 和 R 的多项式，幂为 1、2 和 3。我们将在下文中演示如何将方程 12.6 和 12.7 转换成线性方程。
我们打开方程 12.8 中的圆括号，得到
r--【I】+2 **【x】*
或者
r--【c】和【t】
或者
n (r-【c】*和(12.9)
**对从 0 到 N -1 的 i 的 N 值进行求和计算。总和∑(R2-xc2-yc2)可以替换为*N **(R2-x 既然(R2-xc2-yc2)不依赖于 i 。
现在让我们变换方程 12.6 来分离表达式R2–xc2–yc2。
r-x【I】-
r-x【I】- x 【我】-【r】-**
r-x【I】- x 【我】-(r*
第二项xc*∑(。。。)根据等式 12.8 等于 0。因此:
r-x【I】-
**打开圆括号后，我们得到:
(r--【c】和 x【c】+和【I】*【2】-2 )
*根据等式 12.9， 我们将(R2-xc2-yc2)替换为∑(xI2–2 *x yc*yI)/N 第一项与 x i 相乘第二项:
x【I】【I】【2】 n-(x【I】-2 **【I】*)*
**让我们用 S mn 来表示每一笔总和∑xIm*yIn。然后:
s10(s20–2 **x【c】 s+s*
或者
2** x*c(N * S20–S102+2 *yc(N * S11–S10* S01=(S34(12.6a)
这是一个线性方程在xc和 y c 。经过类似的变换，方程 12.7 也变成一个线性方程:
2** x*c(N * S11–S10 S01+2 *yc*(N * S02–S022=(S04(12.7a)
方程 12.6a 和 12.7a 很容易求解。这些解可以在用 S mn 符号重写的方程 12.9 中设定:
r--【c】和
*以获得 R 的值。
我们在这里展示了执行这些解决方案的方法MinAreaN2的源代码。另外两个类似的稍微简单一点的方法，MinArea2和MinAreaN，也在这个项目中使用。MinAreaN2是这两种方法的结合。
public double MinAreaN2(int ia, iVect2[] P, int Start, int np, ref double radius, ref double x0, ref double y0)
/* Calculates the estimates "x0" and "y0" of the coordinates of the center and the estimate "radius" of the radius of the optimal circle with the minimum deviation from the given set "P[np]" of points. The found values and the 13 sums used for the calculation are assigned to the arc with the index "ia". -- */
{
double SumX, SumY, SumX2, SumY2, SumXY, SumX3, SumY3,
SumX2Y, SumXY2, SumX4, SumX2Y2, SumY4;
double a1, a2, b1, b2, c1, c2, Crit, det, fx, fy, mx, my, N, R2;
int ip;
N = (double)np;
SumX = SumY = SumX2 = SumY2 = SumXY = SumX3 = SumY3 = 0.0;
SumX2Y = SumXY2 = SumX4 = SumX2Y2 = SumY4 = 0.0;
for (ip = Start; ip < Start + np; ip++) //==== over the set of points ====
{
fx = (double)P[ip].X;
fy = (double)P[ip].Y;
SumX += fx;
SumY += fy;
SumX2 += fx * fx;
SumY2 += fy * fy;
SumXY += fx * fy;
SumX3 += fx * fx * fx;
SumY3 += fy * fy * fy;
SumX2Y += fx * fx * fy;
SumXY2 += fx * fy * fy;
SumX4 += fx * fx * fx * fx;
SumX2Y2 += fx * fx * fy * fy;
SumY4 += fy * fy * fy * fy;
} //============ end for (ip...) =============================
a1 = 2 * (SumX * SumX - N * SumX2);
b1 = 2 * (SumX * SumY - N * SumXY);
a2 = 2 * (SumX * SumY - N * SumXY);
b2 = 2 * (SumY * SumY - N * SumY2);
c1 = SumX2 * SumX - N * SumX3 + SumX * SumY2 - N * SumXY2;
c2 = SumX2 * SumY - N * SumY3 + SumY * SumY2 - N * SumX2Y;
det = a1 * b2 - a2 * b1;
if (Math.Abs(det) < 0.00001) return -1.0;
mx = (c1 * b2 - c2 * b1) / det;
my = (a1 * c2 - a2 * c1) / det;
R2 = (SumX2 - 2 * SumX * mx - 2 * SumY *my + SumY2) / N + mx*mx + my*my;
if (R2 <= 0.0) return -1.0;
x0 = mx;
y0 = my;
radius = Math.Sqrt(R2);
elArc[ia].Mx = mx;
elArc[ia].My = my;
elArc[ia].R = radius;
elArc[ia].nPoints = np;
Crit = 0.0;
for (ip = Start; ip < Start + np; ip++) //======== point array =========
{
fx = (double)P[ip].X; fy = (double)P[ip].Y;
Crit += (radius - Math.Sqrt((fx - mx) * (fx - mx) + (fy - my)* (fy - my))) *
(radius - Math.Sqrt((fx - mx) * (fx - mx) + (fy - my) * (fy - my)));
} //============== end for (ip...) ==========================
elArc[ia].SUM[0] = N;
elArc[ia].SUM[1] = SumX;
elArc[ia].SUM[2] = SumY;
elArc[ia].SUM[3] = SumX2;
elArc[ia].SUM[4] = SumY2;
elArc[ia].SUM[5] = SumXY;
elArc[ia].SUM[6] = SumX3;
elArc[ia].SUM[7] = SumY3;
elArc[ia].SUM[8] = SumX2Y;
elArc[ia].SUM[9] = SumXY2;
elArc[ia].SUM[10] = SumX4;
elArc[ia].SUM[11] = SumX2Y2;
elArc[ia].SUM[12] = SumY4;
return Math.Sqrt(Crit / (double)np);
} //*************** end MinAreaN2 *********************************/
项目WFcircleReco
相应的项目WFcircleReco包含以下内容:
-
包含必要数据结构和处理方法调用的类
Form1。 -
第七章中描述的边缘检测方法。
-
第十一章中描述的多边形线近似边缘的方法。
-
一种方法
MakeArcs3，用于将多边形线溶解在称为arcs的具有近似恒定曲率的子集中。 -
为每组经过的圆弧寻找最佳匹配圆的方法
MakeCirclesEl。 -
方法
GetAngle为每个圆弧计算从匹配圆的中心到圆弧端点的半径向量之间的角度。 -
方法
MakeCircles根据圆弧点与圆的最小偏差，为每组曲率中心相近的圆弧计算最佳圆。 -
显示面、弧和最佳圆的服务方法。
我们从项目WFcircleReco的描述开始。
项目的形式WFcircleReco
项目WFcircleReco的形式如图 12-2 所示。
表单的代码由五部分组成，每一部分都在用户单击相应的按钮时开始。
图 12-2
项目的形式WFcircleReco
打开图像部分允许用户选择目录和图像。用户可以打开.bmp或.jpg图像。图像将被打开并显示在左侧的图片框中。这部分代码定义了五个工作图像，并显示原始图像。
单击边缘检测后，用户可以定义边缘检测的阈值。他或她可以看到右边图片框中的边缘，并在必要时更改阈值(标准为 20)。
边缘检测根据第七章描述的二值化梯度法调用边缘检测所需的SigmaSimpleUni和ExtremLightUni方法。图像ExtremIm中的边缘被定义为色差或灰度值差大于阈值设置的相邻像素对的序列。色差是具有亮度差符号的颜色通道的强度的绝对差之和。
为了使边缘的处理更简单、更方便，图像CombIm是用组合坐标(见第七章)和双倍的值width和height创建的。
在这个项目中，我们稍微改变了点的标记:我们现在只使用位 0、1 和 2 来指定与点相关的边缘裂纹的数量；不使用位 3、4、5 和 6。位 7 用于标记已经放入队列的点。为了指定边缘的裂缝和点，我们使用第七章中描述的方法LabelCellsSign。检测到的边缘显示在右边的图片框中。
当用户对检测到的边的质量感到满意时，他或她可以单击多边形。因为要识别的圆的种类非常多，所以我们考虑允许半径的间隔的八种变化:最小的圆的半径可以在 10 和 20 像素之间，而最大的圆的半径可以在 180 和 360 像素之间。对于这些间隔中的每一个，定义多边形近似精度的参数值epsilon被指定:它从最小圆的 1.05 像素变化到最大圆的 2.00 像素。
零件多边形包含一个贯穿所有八个变量的for循环。对于每个变量，方法SearchPoly再次开始。该方法测试包含检测到的边缘的图像CombIm，并在每个端点或分支点开始方法ComponPoly。这两个方法是类CListLines(文件CListLines.cs)的成员。当ComponPoly处理完边的所有结束点和分支点后，再次调用它，从不是结束点或分支点的点开始。然后ComponPoly对所有闭合边曲线进行近似。
方法ComponPoly与 WFcompressPal 项目(第章第八部分)中的方法几乎相同，但有以下微小差异:由方法ComponLin检测并保存的包含单个裂纹的短线现在被忽略。它们对图像压缩很重要；然而，它们对于圆识别并不重要。项目WFcompressPal中ComponLin调用的方法TraceLin在这里被方法TraceAppNew替换，它也跟踪边缘的一条线；然而，它使线的多边形近似。该方法的设计方式使其永远不会产生共线点的序列:如果三个连续的近似点似乎共线，则不保存中间点。为此，使用了一个小方法ParArea(参见下面的源代码)。它计算多边形的三个连续顶点Vert上的平行四边形的面积。
下面是TraceAppNew的源代码。省略了调试所需的指令。
public int TraceAppNew(CImage Comb, int X, int Y, double eps, 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. Bits 3 to 6 are not used. Labeled bit 7 indicates that the point should not be used any more. The crack has only one label 1 in bit 0\. This function traces the edge from one end or branch point to another one while changing the parameter "dir". It makes polygonal approximation with precision "eps" and saves STANDARD coordinates in "Vert". ----------*/
{
int br, Lab, rv = 0;
bool BP = false, END = false, deb = false;
bool atSt_P = false;
iVect2 Crack, P, P1, Pold, Pstand, StartEdge, StartLine, Vect;
Crack = new iVect2();
P = new iVect2();
P1 = new iVect2();
Pold = new iVect2();
Pstand = new iVect2();
StartEdge = new iVect2();
StartLine = new iVect2();
Vect = new iVect2();
int iCrack = 0;
P.X = X; P.Y = Y;
Pstand.X = X / 2;
Pstand.Y = Y / 2;
P1.X = Pold.X = P.X;
P1.Y = Pold.Y = P.Y;
StartEdge.X = X / 2;
StartEdge.Y = Y / 2;
StartLine.X = X / 2;
StartLine.Y = Y / 2;
int[] Shift = { 0, 2, 4, 6 };
int StartVert = nVert;
Vert[nVert].X = Pstand.X;
Vert[nVert].Y = Pstand.Y;
nVert++;
Vect = new iVect2();
int CNX = Comb.width;
int CNY = Comb.height;
CheckComb(StartEdge, Pstand, eps, ref Vect);
while (true) //==================================================
{
Crack.X = P.X + Step[dir].X;
Crack.Y = P.Y + Step[dir].Y;
if (Comb.Grid[Crack.X + CNX * Crack.Y] == 0)
{
MessageBox.Show("TraceAppNew, error: dir=" + dir + " the Crack=(" + Crack.X
+ "," + Crack.Y + ") has label 0; ");
}
P.X = P1.X = Crack.X + Step[dir].X;
P.Y = P1.Y = Crack.Y + Step[dir].Y;
Pstand.X = P.X / 2;
Pstand.Y = P.Y / 2;
br = CheckComb(StartEdge, Pstand, eps, ref Vect);
Lab = Comb.Grid[P.X + CNX * P.Y] & 7; // changed on Nov. 1
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;
case 4: BP = true; END = false; rv = 4; break;
}
if (Lab == 2) Comb.Grid[P.X + CNX * P.Y] = 0; // deleting all labels of P
iCrack++;
if (br > 0) //---------------------------------------------------------
{
if (nVert >= MaxVert - 1)
{
MessageBox.Show("Overflow in 'Vert'; X="+X+" Y="+Y + " nVert=" + nVert);
return -1;
}
if (br == 1)
{
if (nVert > (StartVert + 1) && ParArea(nVert, Pold) == 0.0)
{
Vert[nVert - 1].X = Pold.X / 2;
Vert[nVert - 1].Y = Pold.Y / 2;
}
else
{
Vert[nVert].X = Pold.X / 2;
Vert[nVert].Y = Pold.Y / 2;
}
}
else
{
if (nVert>(StartVert + 1) && ParArea(nVert, Pold)==0.0) Vert[nVert-1] = Vect;
else Vert[nVert] = Vect;
} //------------------- end if (br == 1) ----------------------------
if (nVert > (StartVert + 1) && ParArea(nVert, Pold) == 0.0)
{
StartEdge = Vert[nVert - 1];
}
else
{
StartEdge = Vert[nVert];
if (StartEdge.X > 0 || StartEdge.Y > 0) nVert++; // very important!
}
br = 0;
} //---------------- end if (br > 0) ----------------------------------
atSt_P = (Pstand == StartLine);
if (atSt_P)
{
Pterm.X = P.X; // Pterm is a parameter of TraceAppNew
Pterm.Y = P.Y;
Polygon[nPolygon].lastVert = nVert - 1;
Polygon[nPolygon].closed = true;
rv = 2;
break;
}
if (!atSt_P && (BP || END))
{
Pterm.X = P.X;
Pterm.Y = P.Y;
Vert[nVert].X = Pstand.X;
Vert[nVert].Y = Pstand.Y;
Polygon[nPolygon].lastVert = nVert;
Polygon[nPolygon].closed = false;
nVert++;
if (BP) rv = 3;
else rv = 1;
break;
}
if (!BP && !END) //---------------------------
{
Crack.X = P.X + Step[(dir + 1) % 4].X;
Crack.Y = P.Y + Step[(dir + 1) % 4].Y;
if (Comb.Grid[Crack.X + CNX * Crack.Y] == 1)
{
dir = (dir + 1) % 4;
}
else
{
Crack.X = P.X + Step[(dir + 3) % 4].X;
Crack.Y = P.Y + Step[(dir + 3) % 4].Y;
if (Comb.Grid[Crack.X + CNX * Crack.Y] == 1) dir = (dir + 3) % 4;
}
}
else break;
Pold.X = P.X;
Pold.Y = P.Y;
} //================== end while ==============================
Polygon[nPolygon].nCrack = iCrack;
return rv;
} //********************** end TraceAppNew ****************************
现在Form1调用方法CheckSmooth，检查所有多边形以确定它们是否平滑。如果多边形的边之间的大角的相对数量小于预定义的比例，则该多边形是平滑的。
然后开始方法MakeArcs3，该方法将多边形线细分成称为arcs的具有近似恒定曲率的子集。方法ThreePoints为多边形线的每三个后续顶点计算通过这三个点的圆的曲率。它的值是 1 除以这个圆的半径。该方法计算曲率中心的坐标，并返回曲率的值。在这个项目中，一个弧可以包含任意数量的顶点，这与项目WFpolyArc不同，在项目WFpolyArc中，一个弧总是正好有三个顶点。
方法MakeArcs3实现了弧的如下定义:
-
弧的所有顶点处的局部曲率的所有值具有相同的符号。曲率等于零不能出现是因为前面提到的
TraceAppNew的性质。 -
取决于弧的局部曲率和近似的精度，后续多边形的边之间的角度的绝对值具有小于阈值的近似值。
-
多边形边的长度必须小于阈值。
方法MakeArcs3为多边形线中的每个位置计算四个逻辑变量——bAngle、bCurve、bEdge和bSign——指定是否满足刚刚指示的条件。如果曲率的符号改变，变量bSign为true。
如果所有四个逻辑变量都正常，则某些值被分配给变量Start[iSign]和End[iSign]，指定弧的起点和终点位置。变量bSign指定了局部曲率的符号。圆弧的参数保存在数组elArc中，该数组是类CListLines的成员。
方法MakeCirclesEl然后计算具有相似参数的多组圆弧的最佳圆。它测试所有成对的圆弧，并检查它们的中心点彼此之间的距离是否小于预定义的距离MD以及它们的曲率半径是否近似相等:关系必须在 0.5 和 2.0 之间。满足这些条件的弧的索引保存在一个数组中，并用于计算包含所有这些弧的最佳圆。
为了计算最佳圆，该项目使用了通过MinAreaN2方法实现的最小二乘法修正程序。本章前面介绍了程序和方法。该过程非常快速和鲁棒。
项目WFcircleReco有一个重要特性，即使用多个具有相似圆心和相似半径的圆弧来计算一个圆。它们被放在一个组中，并且为该组计算最佳圆。为了找到所有具有相似中心和相似半径的弧，所有的弧对都被测试。
该项目在许多不同的图像上进行了测试，其中一些是为检查晶片中的焊料凸点而设计的原始项目的旧图像，一些新图像包含苹果、蘑菇、坚果和其他圆形物体。图 12-3 和图 12-4 是已经识别出所有圆形物体并正确估计其参数的图像示例。
图 12-4
认可苹果的项目形式
图 12-3
具有 93 个公认的晶圆凸块圆圈的表单
在左侧的pictureBox1中，识别出的圆圈用红色画在原始图像上。在pictureBox2的右侧，检测到的多边形以绿色绘制。
单击保存圆将所有已识别圆的参数列表保存在文本文件中。通过自动和定义选项，用户可以选择两种可能的方式来指定该文件的名称和位置。如果选择自动，文本文件的名称基于打开图像的名称，同时添加位于包含打开图像的目录images附近的目录的名称TEXT。打开的图像名称的扩展名.bmp或.jpg被替换为扩展名.txt。
如果选择了“已定义”,则用户可以查看目录列表，并通过单击其中一个可见文件来选择目录和文本文件的名称。如果没有合适的名称，用户可以在相应的字段中输入所需的名称。文本文件保存后，其内容通过MessageBox显示。******************
十三、交通中自行车的识别
由于所提出的圆识别方法的良好特性，我有了使用该方法来识别自行车车轮的想法，自行车车轮是理想的圆。然而，如果自行车的位置使其框架平面与观察方向成锐角，那么车轮看起来就像椭圆形而不是圆形。因此，我们也需要一种识别椭圆的方法。不幸的是，我们还没有成功地将我们的圆识别方法推广到椭圆上(见第十二章)。我们已经尝试使用众所周知的共轭梯度法，但是我们的实验表明，这种方法并不可靠:当要拟合椭圆的点不在椭圆附近时，这种方法有时会失败。
因为椭圆仅由少数几个参数定义，即五个参数，所以可以使用下一节所述的经典最小二乘法。
椭圆识别的数学基础
轴平行于笛卡尔坐标系的坐标轴且中心位于原点的椭圆具有众所周知的等式
x2/a2+y2/b2= 1。
然而，我们需要考虑移动和倾斜椭圆的一般情况。我们使用圆锥曲线的一般方程:
ax+bxy+cy**+dx+ey+f【注**
我们的目标是找到一个椭圆的参数，对于该参数，一组给定点到该椭圆的距离的平方和是最小的。因此，我们的目标函数是
)
(13.1)
括号中的表达式与椭圆上的点( x i ， y i )的距离近似成正比。它包含六个未知系数 A 、 B 、 C 、 D 、 *E、*和 F 。然而，众所周知，椭圆由五个参数唯一定义。因此，我们将方程 13.1 的所有项除以 A ，并将新系数表示如下:
B/A =2k1
C/A = k2
D/A =2k3
E/A =2k4
F/A = k5
转换后的目标函数是
)
(13.1a)
f 对k1 的偏导数为
)
将括号中的所有项乘以xI**y*I后，我们得到
)
我们将它除以 4，并用 S ( m ， n 表示每个
)，得到
)
通过将它设置为零，我们获得了未知数k1k2k3k4和 k 5 的五个方程中的第一个。
)
以类似的方式，我们得到其他四个方程:
)
)
)
)
该方程组通过使用方法Gauss_K在方法GetEllipseNew中实现的众所周知的高斯方法求解。总和S(xm， y n )通过MakeSums2方法计算。
以下是GetEllipseNew的源代码。
public int GetEllipseNew(Point[] Vert, int iv1, int nPoints1, int iv2, int nPoints2,
ref double Delta, ref double f, ref double a, ref double b, ref double c, ref double d)
{
bool deb = false;
double[,] A = new double[5, 5];
double[,] B = new double[5, 1];
int nSum = 15;
double[] Sum = new double[nSum];
MakeSums2(Vert, iv1, nPoints1, iv2, nPoints2, Sum);
A[0, 0] = 2.0 * Sum[12];
A[0, 1] = Sum[11];
A[0, 2] = 2.0 * Sum[8];
A[0, 3] = 2.0 * Sum[7];
A[0, 4] = Sum[4];
A[1, 0] = 2.0 * Sum[11];
A[1, 1] = Sum[10];
A[1, 2] = 2.0 * Sum[7];
A[1, 3] = 2.0 * Sum[6];
A[1, 4] = Sum[3];
A[2, 0] = 2.0 * Sum[8];
A[2, 1] = Sum[7];
A[2, 2] = 2.0 * Sum[5];
A[2, 3] = 2.0 * Sum[4];
A[2, 4] = Sum[2];
A[3, 0] = 2.0 * Sum[7];
A[3, 1] = Sum[6];
A[3, 2] = 2.0 * Sum[4];
A[3, 3] = 2.0 * Sum[3];
A[3, 4] = Sum[1];
A[4, 0] = 2.0 * Sum[4];
A[4, 1] = Sum[3];
A[4, 2] = 2.0 * Sum[2];
A[4, 3] = 2.0 * Sum[1];
A[4, 4] = Sum[0];
B[0, 0] = -Sum[13];
B[1, 0] = -Sum[12];
B[2, 0] = -Sum[9];
B[3, 0] = -Sum[8];
B[4, 0] = -Sum[5];
Gauss_K(A, 5, B, 1);
f = -0.5 * Math.Atan2(2.0 * B[0, 0], 1.0 - B[1, 0]);
c = (B[0, 0] * B[3, 0] - B[1, 0] * B[2, 0]) / (B[1, 0] - B[0, 0] * B[0, 0]);
d = (B[0, 0] * B[2, 0] - B[3, 0]) / (B[1, 0] - B[0, 0] * B[0, 0]);
Delta = B[1, 0] - B[0, 0] * B[0, 0];
double BigDelta = B[1, 0] * B[4, 0] + B[0, 0] * B[3, 0] * B[2, 0] +
B[0, 0] * B[3, 0] * B[2, 0] - B[2, 0] * B[1, 0] * B[2, 0] - B[3, 0] * B[3, 0] - B[4, 0] * B[0, 0] * B[0, 0];
double S = 1.0 + B[1, 0];
double a2, b2;
double aprim = (1.0 + B[1, 0] + Math.Sqrt((1.0 - B[1, 0]) * (1.0 - B[1, 0]) + 4.0 * B[0, 0] * B[0, 0])) * 0.5;
double cprim = (1.0 + B[1, 0] - Math.Sqrt((1.0 - B[1, 0]) * (1.0 - B[1, 0]) + 4.0 * B[0, 0] * B[0, 0])) * 0.5;
a2 = -BigDelta / aprim / Delta;
b2 = -BigDelta / cprim / Delta;
a = Math.Sqrt(a2);
b = Math.Sqrt(b2);
if (Delta > 0.0) return 1;
return -1;
} //***************** end GetEllipseNew *********************************
项目WFellipseBike
本项目的形式如图 13-1 所示。
图 13-1
项目的形式WFellipseBike
用户单击打开图像，然后选择文件夹和图像。图像出现在左边的图片框中。然后用户单击检测边缘。有一个数字上下工具，用于选择边缘检测的阈值。但是，预选值 20 适用于所有图像，不应更改。程序自动运行。它使用第十二章中描述的方法进行边缘检测和边缘的多边形近似。然后使用方法MakeArcsTwo将多边形细分成圆弧。该方法与第十二章中描述的MakeArcs3方法略有不同。
然后调用方法FindEllipsesMode来寻找自行车车轮的椭圆。该方法使用根据点数排序的圆弧。排序是通过方法SortingArcs将排序后的弧的索引写入数组SortArcs来执行的。点数最多的圆弧停留在SortArcs[0]。因此，圆弧可以按照点数递减的顺序来取。
椭圆的识别不像圆那么简单:如果圆弧转移到方法GetEllipseNew寻找包含该圆弧的椭圆包含的点少于十个，有时会计算出椭圆的假参数。因此，我们开发了其他方法来解决这个问题。方法QualityOfEllipseNew计算位于椭圆附近的圆弧的点数，如图 13-2 所示。
图 13-2
带有已识别椭圆的自行车图像片段
在图 13-2 中，黑线是圆弧，蓝点是圆弧包含的多边形顶点，红线是识别出的椭圆。椭圆的质量估计为椭圆周围试管中蓝色点的数量Sum。这个数字与值maxPoints进行比较，该值是作为参数传递给方法GetEllipseNew的弧ia中的点数乘以 2π并除以弧的角度。因此maxPoint是封闭曲线中最大可能的点数。如果值Sum大于0.5*maxNumber，则认为椭圆是好的。这里是QualityOfEllipseNew的代码。
public int QualityOfEllipseNew(int ia, Ellipse Ellipse, int[] SortArcs, Form1 fm1)
// Returns the sum of the numbers of points of arcs near the ellipse.
{
bool deb = false, Disp = false; //true; //
int Dif, goodDartIa, locDart, i, iv, ivm, ive, ja, Sum = 0, x, y, xm, ym, xe, ye;
double angleDart, a = Ellipse.a, b = Ellipse.b, c = Ellipse.c, d = Ellipse.d;
double maxPoints = elArc[ia].nPoints * 6.28 / elArc[ia].Angle;
int ivStart, ivMid, ivEnd, xMain, yMain;
ivStart = elArc[ia].Start;
ivMid = ivStart + elArc[ia].nPoints / 2;
ivEnd = ivStart + elArc[ia].nPoints - 1;
x = Vert[ivStart].X;
y = Vert[ivStart].Y;
double AngleStart = Math.Atan2(y - d, x - c);
xe = Vert[ivEnd].X;
ye = Vert[ivEnd].Y;
double AngleEnd = Math.Atan2(ye - d, xe - c);
xMain = Vert[ivMid].X;
yMain = Vert[ivMid].Y;
double AngleMid = Math.Atan2(yMain - d, xMain - c);
double minAngle = Math.Min(AngleStart, AngleEnd);
double maxAngle = Math.Max(AngleStart, AngleEnd), help;
bool Plus2PI = false;
if (minAngle < 0.0 && maxAngle > 0.0 && !(AngleMid >= minAngle &&
AngleMid < maxAngle))
{
Plus2PI = true;
help = maxAngle;
maxAngle = minAngle + 2 * Math.PI;
minAngle = help;
}
angleDart = 57.3 * Math.Atan2(yMain - elArc[ia].My, xMain - elArc[ia].Mx) + 15.0;
if (angleDart < 0.0) angleDart += 360.0;
goodDartIa = 6 + (int)angleDart / 30;
if (goodDartIa > 11) goodDartIa -= 12;
double AngleJa, Fx, Fxe;
for (i = 0; i < nArcs; i++) //=============================
{
ja = SortArcs[i];
if (ja == ia || elArc[ja].nPoints < 5) continue;
iv = elArc[ja].Start;
ivm = iv + elArc[ja].nPoints / 2;
ive = iv + elArc[ja].nPoints - 1;
x = Vert[iv].X;
y = Vert[iv].Y;
xm = Vert[ivm].X;
ym = Vert[ivm].Y;
xe = Vert[ive].X;
ye = Vert[ive].Y;
Fx = (x - c) * (x - c) / a / a + (y - d) * (y - d) / b / b;
Fxe = (xe - c) * (xe - c) / a / a + (ye - d) * (ye - d) / b / b;
if (Fx < 0.6 || Fx > 1.67 || Fxe < 0.6 || Fxe > 1.67) continue;
angleDart = 57.3 * Math.Atan2((ym - d) * a * a, (xm - c) * b * b);
if (angleDart < 0.0) angleDart += 360.0;
locDart = (int)angleDart / 30;
if (locDart > 11) locDart -= 12;
Dif = Math.Abs(elArc[ja].Dart - locDart);
if (Dif > 6) Dif = 12 - Dif;
if (Disp) DrawOneLongArc(ja, fm1);
if (Dif < 2)
{
if (Disp) DrawOneLongArc(ja, fm1);
for (iv = elArc[ja].Start; iv < elArc[ja].Start + elArc[ja].nPoints; iv++)
{
x = Vert[iv].X;
y = Vert[iv].Y;
AngleJa = Math.Atan2(y - d, x - c);
if (AngleJa < 0.0 && Plus2PI) AngleJa += 6.28;
if (!(AngleJa > minAngle && AngleJa < maxAngle)) Sum += elArc[ja].nPoints;
}
}
} //================== end for (i = 0; ... ===========================
return Sum;
} //********************* end QualityOfEllipseNew ************************
如果椭圆不好，则调用方法HelpArcNew。该方法获取圆弧ia作为参数，并找到位于圆弧ia的曲率圆内的不同圆弧ja。如果圆弧ja方向正确且不太靠近圆弧ia，则与圆弧ia形成一对。对每一对圆弧(ia, ja)计算一个椭圆并评估其质量。具有最佳质量的椭圆作为结果返回。这里是HelpArcNew的代码。
public int HelpArcNew(int ia, int[] SortArcs, ref Ellipse Ellipse, int SumStart, Form1 fm1)
{
bool disp = false;
int Dif, i, ivMid, ivm, ja, xMain, yMain, xm, ym;
ivMid = elArc[ia].Start + elArc[ia].nPoints / 2;
xMain = Vert[ivMid].X;
yMain = Vert[ivMid].Y;
double angleDart, a = Ellipse.a, b = Ellipse.b, c = Ellipse.c, d = Ellipse.d;
int goodDartIa;
angleDart = 57.3 * Math.Atan2(yMain - elArc[ia].My, xMain - elArc[ia].Mx) + 15.0;
if (angleDart < 0.0) angleDart += 360.0;
goodDartIa = 6 + (int)angleDart / 30;
if (goodDartIa > 11) goodDartIa -= 12;
double R = elArc[ia].R, Mx = elArc[ia].Mx, My = elArc[ia].My;
CBox Box = new CBox();
Box.minX = (int)(Mx - R) - 10;
if (c - a < Box.minX) Box.minX -= 20;
Box.maxX = (int)(Mx + R) + 10;
if (c + a > Box.maxX) Box.maxX += 20;
Box.minY = (int)(My - R) - 10;
if (d - b < Box.minY) Box.minY -= 20;
Box.maxY = (int)(My + R) + 10;
if (d + b > Box.maxY) Box.maxY += 20;
DrawRectangleSmart(Box, fm1);
int Dist2 = 0, minDist2 = (int)(1.5 * elArc[ia].R * elArc[ia].R);
int[] jBest = new int[100];
int nBest = 0;
for (i = 0; i < nArcs; i++) //===================================
{
ja = SortArcs[i];
if (!ArcInBox(ja, Box)) continue;
if (elArc[ja].nPoints < 4) continue;
if (disp) DrawOneLongArc(ja, fm1);
Dif = Math.Abs(elArc[ja].Dart - goodDartIa);
if (Dif > 6) Dif = 12 - Dif;
if (Dif > 3) continue;
ivm = elArc[ja].Start + elArc[ja].nPoints / 2;
xm = Vert[ivm].X;
ym = Vert[ivm].Y;
Dist2 = (xm - xMain) * (xm - xMain) + (ym - yMain) * (ym - yMain);
if (Dist2 > minDist2)
{
jBest[nBest] = ja;
nBest++;
}
if (nBest >= 5) break;
} //================= end for (i = 0; =============================
double Delta = 0.0, F = 0.0;
Ellipse Ellipse1 = new Ellipse();
int jbestOpt = -1, maxSum = SumStart, Sum = 0;
for (i = 0; i < nBest; i++) //=========================================
{
if (disp) DrawRedArc(jBest[i], fm1);
GetEllipseNew(Vert, elArc[ia].Start, elArc[ia].nPoints, elArc[jBest[i]].Start,
elArc[jBest[i]].nPoints, ref Delta, ref Ellipse1.f,
ref Ellipse1.a, ref Ellipse1.b, ref Ellipse1.c, ref Ellipse1.d);
Sum = QualityOfEllipseNew(ia, Ellipse1, SortArcs, fm1);
if (disp) DrawEllipse(Ellipse1, fm1);
if (!(Ellipse1.a > 5.0 && Ellipse1.b > 5.0) ||
Ellipse1.d - Ellipse1.b < fm1.height * 2 / 5) Sum = 0;
else
{
if (Sum > maxSum)
{
if (disp) DrawEllipse(Ellipse1, fm1);
if (disp) DrawRedArc(jBest[i], fm1);
maxSum = Sum;
jbestOpt = jBest[i];
Ellipse = Ellipse1;
}
}
} //================== end for (i ... i < nBest; ======================
DrawRedArc(jbestOpt, fm1);
Pen pen = new Pen(Color.Red);
DrawEllipsePen(Ellipse, pen, fm1);
return jbestOpt;
} //******************** end HelpArcNew *******************************
自行车的后轮有时会被骑自行车的人的腿遮住，因此无法检测到椭圆。如果已经识别出前轮的高质量椭圆，则将该椭圆的副本分配给后轮。为了正确完成这一分配，必须计算后轮圆弧ia1中点MP处的切线。前轮的椭圆上有一个点P，与圆弧ia1的方向相切且方向相同。放置前轮椭圆的副本，使点P位于点MP上。然后，用方法HelpArcNew指定椭圆副本的参数。在调用HelpArcNew之前和之后，使用方法CminusCel检查两个椭圆的相对位置，尤其是它们彼此之间的距离。所有这些程序都在方法FindEllipsesMode中执行。下面是这个方法的代码。
public int FindEllipsesMode(CImage SigmaIm, Ellipse[] ListEllipse, ref int nEllipse, Form1 fm1)
{
int[] SortArcs = new int[nArcs];
int maxNP = 0, k = SortingArcs(SortArcs, ref maxNP);
int i, ia, ia1, i0, i1;
nEllipse = 0;
double a = 0.0, b = 0.0, c = 0.0, d = 0.0; //, fret = 0.0;
int[,] List = new int[20, 1200];
int[] nArcList = new int[20];
SCircle[] Circle = new SCircle[20];
for (i = 0; i < 20; i++)
{
Circle[i] = new SCircle();
Circle[i].goodCirc = true;
}
Ellipse[] smalList = new Ellipse[20];
for (i = 0; i < 20; i++) smalList[i] = new Ellipse();
int Sum1 = 0;
double AnglePerPoint = 0.0, maxPoints = 0.0;
fm1.progressBar1.Visible = true;
fm1.progressBar1.Step = 1;
int jump, Len = nArcs, nStep = 20;
if (Len > 2 * nStep) jump = Len / nStep;
else
jump = 2;
double Delta = 0.0, f = 0.0, F = 0.0;
Ellipse Ellipse1 = new Ellipse();
Ellipse Ellipse2 = new Ellipse();
int[] Pattern = new int[100000];
double aa = 0.0, bb = 0.0, cc = 0.0, dd = 0.0;
for (i0 = 0; i0 < nArcs; i0++) //======================================
{
if ((i0 % jump) == jump - 1) fm1.progressBar1.PerformStep();
ia = SortArcs[i0];
DrawRedArc(ia, fm1);
if (elArc[ia].nPoints <= 5) break;
GetEllipseNew(Vert, elArc[ia].Start, elArc[ia].nPoints, 0, 0, ref Delta, ref f,
ref a, ref b, ref c, ref d);
DrawEllipse(f, a, b, c, d, fm1);
if (b < 20.0 || a < 6.0 || d + b > fm1.height || d - 4 * b < 0.0) continue;
int jbestOpt = -1;
Ellipse1.a = a;
Ellipse1.b = b;
Ellipse1.c = c;
Ellipse1.d = d;
Point P1 = new Point(0, 0);
Point P2 = new Point(0, 0);
if (a > 5.0 && b > 5.0)
{
Sum1 = QualityOfEllipseNew(ia, Ellipse1, SortArcs, fm1);
AnglePerPoint = elArc[ia].Angle / elArc[ia].nPoints;
maxPoints = 2 * Math.PI / AnglePerPoint;
if (b > fm1.height / 4 || elArc[ia].nPoints < 10) Sum1 = 0;
Pen pen = new Pen(Color.Red);
if (elArc[ia].nPoints < 10 || d + b > fm1.height * 2 / 5)
{
jbestOpt = HelpArcNew(ia, SortArcs, ref Ellipse1, Sum1, fm1);
DrawEllipse(Ellipse1, fm1);
}
}
for (i1 = i0 + 1; i1 < nArcs; i1++) //=================================
{
ia1 = SortArcs[i1];
if (!Position(ia1, Ellipse1, fm1)) continue;
if (elArc[ia1].nPoints <= 5)
{
MessageBox.Show("Finishing the search for ia1");
break;
}
CBox BoxP1 = new CBox();
CBox BoxP2 = new CBox();
int iv = elArc[ia].Start, x = Vert[iv].X, y = Vert[iv].Y;
int iv1 = elArc[ia1].Start, x1 = Vert[iv1].X, y1 = Vert[iv1].Y;
GetEllipseNew(Vert, elArc[ia1].Start, elArc[ia1].nPoints, 0, 0, ref Delta, ref f,
ref a, ref b, ref c, ref d);
DrawEllipse(f, a, b, c, d, fm1);
if (!(a > 5.0 && b > 5.0)) continue;
DrawRedArc(ia1, fm1);
double K2 = DrawTangent(ia1, ref P2, fm1);
P1 = PointWithTangent(Ellipse1, K2, elArc[ia1].Dart, fm1);
if (P1.X == 0 && P1.Y == 0) continue;
Ellipse2 = Ellipse1;
Ellipse2.c = P2.X + Ellipse1.c - P1.X;
Ellipse2.d = P2.Y + Ellipse1.d - P1.Y;
DrawEllipse(Ellipse2, fm1);
if (d + b > fm1.height || d - b < fm1.height * 2 / 5) continue;
jbestOpt = -1;
if (Ellipse2.a > 5.0 && Ellipse2.b > 5.0)
{
Sum1 = QualityOfEllipseNew(ia1, Ellipse2, SortArcs, fm1);
AnglePerPoint = elArc[ia1].Angle / elArc[ia1].nPoints;
maxPoints = 2 * Math.PI / AnglePerPoint;
Pen pen = new Pen(Color.Red);
if ((elArc[ia1].nPoints < 10 || d + b > fm1.height * 2 / 5) &&
CminusCel(Ellipse1, Ellipse2, fm1))
{
jbestOpt = HelpArcNew(ia1, SortArcs, ref Ellipse2, Sum1, fm1);
DrawEllipse(Ellipse2, fm1);
}
}
bool CMINC = CminusCel(Ellipse1, Ellipse2, fm1);
if (!CMINC) continue;
ListEllipse[nEllipse] = Ellipse1;
nEllipse++;
ListEllipse[nEllipse] = Ellipse2;
nEllipse++;
if (nEllipse >= 2)
{
fm1.progressBar1.Step = fm1.progressBar1.Maximum - Len / (100 / 6) -
fm1.progressBar1.Value;
fm1.progressBar1.PerformStep();
return 1;
}
} //============== end for (i1 = 0; ... ============================
} //=============== end for (i0 = 0; ... =============================
MessageBox.Show("FindEllipsesMode: no bike recognized");
return -1;
} //****************** end FindEllipseMode ******************************
车轮的识别是自行车识别中最重要的部分。当车轮的两个椭圆都被识别时，方法RecoFrame被调用。它能识别自行车的行驶方向。这种方法包含几种不同类型的框架的简单模型。每种车架出现两次:一次是自行车前轮的前叉在左边，另一次是前叉在右边。这些模型被一个接一个地测试。每个模型都经过变换，使得轮轴的位置与椭圆的中心相匹配。在框架的每个直线段周围形成一个狭窄的矩形。然后，经过处理的图像中的所有多边形都被通过，并且适合矩形的多边形边的长度被求和。具有最大总和的模型获胜。这样，自行车的运动方向就被识别出来了。图 13-3 中显示了框架的模型示例。
图 13-3
框架模型的示例
识别方向的另一种方法
刚刚描述的方法在框架平面与观察方向成锐角时有效，例如图 13-1 。此外，框架的某些部分经常被骑车人的腿遮挡。因此，我们开发了另一种方法黑斑，返回自行车运动的方向:0 为右，2 为左。该方法在每个车轮内制作一个框，计算亮度直方图，并对这些框中几乎水平的多边形边的长度求和，并决定自行车的绘制方向。通过比较低于阈值的直方图值的总和，即框中暗像素的数量加上两个框中的 6 *总和，来做出决定。较大的值表示后轮，也表示绘制方向:如果后轮在右手边，则自行车向左移动。
当识别出车轮的椭圆和运动方向时，车轮的椭圆和车架的型号会显示在右侧的图片框中。显示屏中间的消息框显示消息“识别到向左行驶的自行车”如果用户单击保存结果，将在磁盘上保存一个文本文件，其上下文对应于以下示例:
“向左的自行车有椭圆形的轮子。
第一:a = 113b = 173c = 864d = 796。
第二:a = 114b = 173c = 469d = 790。"
记法如下:a 是椭圆的水平半轴，b 是其垂直半轴，c 是车轮中心的 x 坐标，d 是其 y 坐标。
我们测试了大约 100 张显示自行车在城市街道上行驶的照片。除了大约 2000 × 1500 像素的大图像和大约 150 像素长的小自行车(图 13-4 )之外，自行车在所有照片中都被识别出来。
图 13-4
带有未识别自行车的图像示例
图 13-5 提供了一些已识别自行车的图像示例。
图 13-5
识别自行车的图像示例
十四、细胞分化的计算机模型
本章不属于图像处理领域。之所以把它收入本书，是因为作者的兴趣圈广泛而多才多艺。这一章的主题与生物学有关。
人们可能会问，尽管事实上不同类型的细胞都具有完全相同的脱氧核糖核酸(DNA)，但它们如何可能在有机体的生长过程中发育成不同类型的细胞。可以假定，每个细胞都必须获得一些信息，表明它在生长的有机体中的位置，也许是某种坐标。然后，DNA 必须包含不同的指令，用于具有不同坐标的细胞的发育。坐标可以在单元划分时产生:与具有坐标集( X，Y，Z )的旧单元相邻的新单元必须得到旧单元的其中一个坐标增加或减少了 1 的坐标集。改变三个坐标中的哪一个取决于向量从旧单元到新单元的方向。因此，这一章的两个主要观点是，细胞必须获得某种坐标，而指定细胞特性的 DNA 指令必须依赖于它的坐标。
我们已经开发了一个项目WFcellDivision，其中多细胞生物由 63 × 63 个细胞的阵列Org(生物)建模(奇数更好地精确定义中间的细胞)。
通过指定和保存其特征来模拟产生新单元的过程。每个细胞可以拥有坐标 X 和 Y 作为其特征，一个变量Property(这是一个颜色索引)，以及数组DNA[63×63]作为其遗传信息。单元格是CCelln类的一个对象。
int const Cwidth=63, Cheight=63;
class CCelln
{ public:
int X, Y, Property;
unsigned char DNA[Cwidth*Cheight];
}
该项目的目的是模拟一个接一个细胞的生成，并为每个细胞的变量Property赋值。该值取决于该单元格的坐标以及保存在该单元格中的DNA的内容。DNA``to Properties of all cells of Org is supposed to be impossible. It is only possible to copy the DNA from one cell to the neighboring new originating cell and to copy inside the originating cell one certain value from the DNA of the cell to the Property of this cell. The aim of the project is to demonstrate that in spite of these limitations, it is possible to assign correct values of Property to all cells of Org.
让我们描述一下这个项目的运作。方法Form1从定义和初始化拥有 63 × 63 个单元的存储器的数组Org开始。初始化Org包括将Org的每个单元的Property和数组DNA[63×63]的元素设置为−1，这意味着该单元不存在。
然后数组Org中间的一个单元格，即坐标为(31，31)的单元格被“填充”:它获得坐标 X = 31 和 Y = 31。它的数组DNA填充了一个 63 × 63 像素的小数字Image的颜色索引。该图像是任意选择的，以证明具体的内容可以分配给细胞的属性。DNA的内容成为数字图像的副本，其中每个元素是一个颜色索引，一个 0 到 255 之间的数字，可以通过颜色表Palette转换成一种颜色。在Palette中，为每个颜色索引保存一种颜色(红、绿、蓝)。数组Org的所有其他单元最初的坐标等于−1，数组DNA为空(即，用−1填充)。
接下来调用方法Grow。它模拟生物体的生长，同时从与中央细胞(31，31)相邻的细胞开始生成细胞坐标。坐标以这样的方式生成，即所有原始单元都位于围绕中心单元旋转的螺旋中。在生长过程中，越来越多的单元从相邻单元之一获得数组DNA的副本。它们还根据DNA和单元坐标的内容获得它们的坐标和某些Property值。
坐标为(31，31)的单元格的Property获取数组DNA(见图 14-1 、DNA[31, 31]的单个元素的值，即像素(31，31)的颜色索引。
图 14-1
将图像复制到Org[31, 31].DNA
生长生物体的形状由代表生长生物体的颜色区域周围的DNA的一些元素分配的黑色区域的边界指定。该数组在对应于应该存在于生物体中的细胞的每个元素中包含不同于(0，0，0)的颜色索引，并且在对应于不应该存在于生物体中的细胞的每个元素中包含指向(0，0，0)(黑色)的索引。
如果数组Org的一个单元格具有非零坐标，值Property指向不同于(0，0，0)的颜色，并且DNA的内容不同于-1，则该单元格被填充。
当生长过程开始时，数组Org中与填充单元格相邻的每个单元格都获得数组DNA的精确副本，该副本对所有单元格都是相同的。该单元的坐标获得的值与填充的相邻单元的坐标相差 1。例如，如果新单元格位于坐标为( X，Y )的已填充单元格的右侧，则它的 X 坐标获得值 X 1 = X + 1，它的 Y 坐标 Y 1 获得值 Y 。新单元获得了DNA的标准副本，其Property等于DNA[X1, Y1]，其中 X 1 和 Y 1 是新单元的坐标。
这里是Grow的代码。
int Grow(CCelln[] Org, int width, int height)
{ int cnt=0, i, k, x=(width-1)/2, y=(height-1)/2, nCount=1, X, Y;
do
{ for (i=0; i<nCount && x<width-1; i++)
{ x++;
X=Org[x+width*y].X=Org[x-1+width*y].X+1;
Y=Org[x+width*y].Y=Org[x-1+width*y].Y;
for (k=0; k<width*height; k++) Org[x+width*y].DNA[k] =
Org[x-1+width*y].DNA[k];
Org[x+width*y].Property=Org[x+width*y].DNA[X+width*Y];
}
cnt+=nCount;
if (cnt==width*height) break;
for (i=0; i<nCount && y<height-1; i++)
{ y++;
X=Org[x+width*y].X=Org[x+width*(y-1)].X;
Y=Org[x+width*y].Y=Org[x+width*(y-1)].Y+1;
for (k=0; k<width*height; k++) Org[x+width*y].DNA[k] =
Org[x+width*(y-1)].DNA[k];
Org[x+width*y].Property=Org[x+width*y].DNA[X+width*Y];
}
cnt+=nCount;
if (cnt==width*height) break;
nCount++;
for (i=0; i<nCount && x>0; i++)
{ x--;
X=Org[x+width*y].X=Org[x+1+width*y].X-1;
Y=Org[x+width*y].Y=Org[x+1+width*y].Y;
for (k=0; k<width*height; k++)
Org[x+width*y].DNA[k]=Org[x+1+width*y].DNA[k];
Org[x+width*y].Property=Org[x+width*y].DNA[X+width*Y];
}
cnt+=nCount;
if (cnt==width*height) break;
for (i=0; i<nCount && y>0; i++)
{ y--;
X=Org[x+width*y].X=Org[x+width*(y+1)].X;
Y=Org[x+width*y].Y=Org[x+width*(y+1)].Y-1;
for (k=0; k<width*height; k++)
Org[x+width*y].DNA[k]=Org[x+width*(y+1)].DNA[k];
Org[x+width*y].Property=Org[x+width*y].DNA[X+width*Y];
}
cnt+=nCount;
if (cnt==width*height) break;
nCount++;
} while(1);
return 1;
}
当数组Org的所有单元都被填充时，该过程结束。然后将所有单元格的Property值复制到一个 63 × 63 像素的新数组(图像)中，目的是使结果可见。该图像与颜色表Palette一起显示，用户可以看到该图像与复制到DNA[31, 31]的原始图像Image相同，除了数组DNA中的像素具有指向黑色的索引。这些像素显示为黑色(参见图 14-2 )。
图 14-2
从Org[*].Property复制的结果图像
因此，数组Org的所有单元都是不同的:它们包含不同的Property值，尽管它们都具有数组DNA的完全相同的副本。
结论
我们已经证明，尽管生物体的所有细胞都具有完全相同的遗传信息拷贝，并且尽管禁止将数组DNA直接复制到数组Org[x, y].Property中，但是产生具有不同细胞的生物体是可能的。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
