获取骨骼外轮廓点集
进入3D建模界面，主界面控制类ConstructController在工具栏上显示了5个按钮返回、网格重建、清除网格、钢板创建、保存图片，以及两个滚动条：
网格重建即是三维建模，在建模前，我们需要在左边的骨骼CT图像中框取腿骨区域，称为加速框。
加速框的作用在于缩小建模的范围，排除无关的杂质。同时只对加速框内部的图像进行轮廓计算，随着CT图像序列的计算，动态地调整加速框的大小，减少像素的计算量。
本系统的三维建模采用的是二维轮廓线表面拟合算法，因此轮廓线的提取显得尤为重要。本章重点描述单张CT图像的骨骼外轮廓点的提取算法。
用户点击三维重建，触发按钮响应函数ConstructMesh，启动伪线程ConstructCoroutine：
public void ConstructMesh()
{
StartCoroutine(ConstructCoroutine());
}
ConstructCoroutine会完成整个的三维建模过程。它主要包含两个功能，一是轮廓提取，二是轮廓三角面片生成。轮廓提取又分为5个步骤：
- 阈值边缘点提取
- 凸包检测
- 轮廓校正
- 移除过于密集的点
- 贝塞尔平滑
通过前面设置的阈值ContourValue可以初步获得边缘点。获取的方式类似于上一章所说的IsoLine.PickUpContour16方法，但是不同的是采用了加速框进行图像范围的约束。加速框会在每次计算完一张CT图像后动态的改变长宽，改变的规则是基于当前CT图像所有边缘点的最小点和最大点。当然，加速框不可能刚好包住当前CT图像的所有边缘点，这样对于下一张图像就可能显小了，因此还要稍微放大一点。基于加速框的边缘提取调用的函数是IsoLine.PickUpContour16WithAccerlateBox：
List<Vector3> listContour = IsoLine.PickUpContour16WithAccerlateBox(ref srcTexture, ref initRect, pix16, displayTunning.Lut16, (byte)contourValue, z, 1);
public static List<Vector3> PickUpContour16WithAccerlateBox(ref Texture2D texture, ref Rect rect, List<ushort> listPix, byte[] lut16, byte contourValue, float z, int step)
{
float currTime = Time.realtimeSinceStartup;
int v1, v2, v3, small, medium, large;
int[] X = new int[5];
int[] Y = new int[5];
float[] value = new float[5];
List<Vector3> listVector = new List<Vector3>();
int width = texture.width;
int height = texture.height;
Vector2 min = new Vector2(width, height);
Vector2 max = new Vector2(0, 0);
int count = 0;
for (int y = (int)rect.yMin; y < rect.yMax - 1; y++) // 遍历加速框内的像素点
{
for (int x = (int)rect.xMin; x < rect.xMax - 1; x++)
{
int index = y * width + x;
X[0] = x; Y[0] = y; value[0] = lut16[listPix[index]];
X[1] = x + 1; Y[1] = y; value[1] = lut16[listPix[index + 1]];
X[2] = x + 1; Y[2] = y + 1; value[2] = lut16[listPix[index + width + 1]];
X[3] = x; Y[3] = y + 1; value[3] = lut16[listPix[index + width]];
X[4] = (X[0] + X[1]) / 2;
Y[4] = (Y[1] + Y[2]) / 2;
value[4] = 0.25F * (value[0] + value[1] + value[2] + value[3]);
v3 = 4;
for (v1 = 0; v1 < 4; v1++)
{
v2 = v1 + 1;
if (v1 == 3) v2 = 0;
int temp;
large = v1;
medium = v2;
small = v3; // v3表示中心像素点
... // 此处省略，过程是交换small、medium、large的位置以满足大小关系
if (value[small] < contourValue && contourValue < value[large]) // 判断阈值边缘点
{
if (count++ % step == 0)
{
listVector.Add(new Vector3(x - width / 2, y - height / 2, z)); // 将点添加到数组中
min = Vector2.Min(min, new Vector2(x, y)); // 计算最小点
max = Vector2.Max(max, new Vector2(x, y)); // 计算最大点
break;
}
}
}
}
}
listTime.Add(Time.realtimeSinceStartup - currTime); // 记录时间
min = Vector2.Max(Vector2.zero, min - 10 * Vector2.one); // 加速框左上角
max = Vector2.Min(new Vector2(width, height), max + 10 * Vector2.one); // 加速框右下角
rect = new Rect(min, max - min);
return listVector;
}
通过初步的阈值边缘提取，得到的是一组边缘点List<Vector3> listContour。这样提取到的点并不是我们所期望的，这些点不是凸包的，因为骨骼中有骨髓，所以骨骼是有厚度的，而我们暂时只需要骨骼的外表，因此只需要外轮廓点。
接下来就是凸包检测，旨在消除内部轮廓点。凸包检测采用了Graham扫描算法，其核心思想是设置一个候选点的堆栈，每个点都会被压入堆栈一次，压入前，会逐一弹出栈顶点，直到栈中剩余点与当前点形成的连线是单方向折拐的（比如逆时针）。该算法需要对点集进行预处理，使得所有的点都按照逆时针或顺时针排序。
ConvexHull专门负责凸包检测：
listContour = ConvexHull.GrahamScan(listContour);
public static List<Vector3> GrahamScan(List<Vector3> listOrigion)
{
if (listOrigion.Count < 3) return listOrigion;
// GeometryPoint为点结构体，提供斜率比较
List<GeometryPoint> geometrypoints = new List<GeometryPoint>();
// 选取x坐标最小的点
int MinXPointIndex = 0;
for (int i = 0; i < listOrigion.Count; i++)
{
GeometryPoint geometryPoint = new GeometryPoint(listOrigion[i]);
geometrypoints.Add(geometryPoint);
if (geometryPoint.X < geometrypoints[MinXPointIndex].X)
{
MinXPointIndex = i;
}
else if (geometryPoint.X == geometrypoints[MinXPointIndex].X &&
geometryPoint.Y < geometrypoints[MinXPointIndex].Y)
{ // 选取x坐标最小的点，如果最小x坐标点有多个，去y最大者
MinXPointIndex = i;
}
}
Vector3 minXPoint = geometrypoints[MinXPointIndex].ToVector3();
// 计算斜率
for (int i = 0; i < geometrypoints.Count; i++)
{
if (i == MinXPointIndex)
{
geometrypoints[MinXPointIndex].SLOPE = float.MaxValue;
}
else
{
if (geometrypoints[i].X == geometrypoints[MinXPointIndex].X)
{ // 与最小x坐标的x相同的点,因为x坐标之差为零，所以取SLOPE最大值
geometrypoints[i].SLOPE = float.MaxValue;
}
else
{ // 计算斜率，注意正切函数在-0.5Pi和0.5Pi之间是单调递增的
geometrypoints[i].SLOPE = (geometrypoints[i].Y - geometrypoints[MinXPointIndex].Y) / (geometrypoints[i].X - geometrypoints[MinXPointIndex].X);
}
}
}
geometrypoints.RemoveAt(MinXPointIndex);
geometrypoints.Sort(); // 按斜率排序
List<Vector3> listOutput = new List<Vector3>();
// 前三个点先入栈
listOutput.Add(minXPoint);
listOutput.Add(geometrypoints[0].ToVector3());
listOutput.Add(geometrypoints[1].ToVector3());
// 判断与其余所有点的关系
for (int i = 2; i < geometrypoints.Count; i++)
{
Vector3 diff1 = listOutput[listOutput.Count - 1] - listOutput[listOutput.Count - 2];
Vector3 diff2 = geometrypoints[i].ToVector3() - listOutput[listOutput.Count - 1];
while (diff1.x * diff2.y - diff2.x * diff1.y < 0) // 比较斜率判断拐点方向
{
listOutput.RemoveAt(listOutput.Count - 1); // 移除拐点方向不同的点（此处为顺时针）
diff1 = listOutput[listOutput.Count - 1] - listOutput[listOutput.Count - 2];
diff2 = geometrypoints[i].ToVector3() - listOutput[listOutput.Count - 1];
}
listOutput.Add(geometrypoints[i].ToVector3());
}
return listOutput;
}
凸包检测获得了骨骼边缘的外轮廓点List<Vector3> listContour。但我们漏了一种情况，那就是外轮廓可能是有凹陷的，通过凸包检测我们可以去除内轮廓点，同时也去掉了一部分原本属于外轮廓凹陷部分的点。
凸包之后要进行轮廓校正，一是去除外轮廓点中的噪声点，二是修补外轮廓凹陷部分的点。去噪校正算法的核心思想是，骨骼是有内轮廓的，外轮廓的点必然距离内轮廓的点很近，那些偏离内轮廓的点将作为噪声点去除。修补校正算法的核心思想是，判断边缘点间距是否过长，过长表示两点间可能有凹陷部分，取两点中点，计算距离最近的内部点，加入外轮廓点，迭代此过程：
轮廓校正函数位于ConstructController中：
ContourCorrection(listContour);
public void ContourCorrection(List<Vector3> listContour)
{
if (displayTunning.ListListResult.Count > 0)
{
for (int j = 0; j < listContour.Count; j++)
{
// 移除外轮廓中的噪声点
Vector3 point1 = listContour[j];
bool legal = false;
for (int k = 0; k < displayTunning.LastListVector3.Count; k++)
{
float dist = Vector3.Distance(point1, displayTunning.LastListVector3[k]);
if (dist < 30)
{ // 存在间距较小的相邻内部点
legal = true;
break;
}
}
if (!legal)
{ // 移除与其它所有内部点较远的点
listContour.RemoveAt(j);
j--;
continue;
}
// 如果与相邻的下一个点间距较大则进行校正
Vector3 point2 = listContour[(j + 1) % listContour.Count];
if (Vector3.Distance(point1, point2) > 20)
{
Vector3 center = Vector3.Lerp(point1, point2, 0.5f); // 中点
float minDist = float.MaxValue;
Vector3 minDistPoint = center;
for (int k = 0; k < displayTunning.LastListVector3.Count; k++)
{ // 寻找离中点最近的内部点
float dist = Vector3.Distance(center, displayTunning.LastListVector3[k]);
if (dist < minDist)
{
minDistPoint = displayTunning.LastListVector3[k];
minDist = dist;
}
}
Vector2 adjustpoint = Vector3.Lerp(center, minDistPoint, 0.5f);
listContour.Insert(j + 1, minDistPoint); // 加入外轮廓点
j++;
}
}
}
}
现在我们有了正确的外轮廓点，已经可以大致地反映骨骼的外表形状了，但是过于密集的外轮廓点不利于接下来的贝塞尔平滑，所以需要做一下预处理，减少外轮廓点的数量：
IsoLine.RemoveDensePoint(listContour, meshDensityTunning.Horizontal);
public static void RemoveDensePoint(List<Vector3> listVector, float leastDistance)
{
if (listVector.Count < 3) return;
for (int i = 0; i < listVector.Count; i++)
{
int lastIndex = (i - 1 + listVector.Count) % listVector.Count;
int nextIndex = (i + 1 + listVector.Count) % listVector.Count;
if (Vector3.Distance(listVector[lastIndex], listVector[i]) < leastDistance &&
Vector3.Distance(listVector[nextIndex], listVector[i]) < leastDistance)
{ // 该点与前后相邻点距离过小则移除
listVector.RemoveAt(i);
}
}
}
贝塞尔平滑的作用是使得外轮廓点均匀平滑的描述出骨骼轮廓，避免出现尖锐的转折点，或是不均匀的点分布。这有利于三维建模时的平滑。贝塞尔平滑算法的步骤是：
首先计算出所有边的中点，然后按照顺序将这些中点依次连接成线段，这样，多边形的每个顶点就有与之对应的一个线段。依据多边形上过顶点的两条边的比例，在顶点对应的线段上找到一点，然后平移线段直到这个点与顶点齐平，线段的两个端点就是过这个点对应的两条边的控制点。通过顶点的两个相邻点和计算出来的两个控制端点，我们可以利用贝塞尔差值得到两个顶点之间的曲线表达式。依次曲线进行采样，就可以得到平滑的均匀分布的轮廓点。贝塞尔平滑方法由类Bezier负责，详细代码较多，就不再贴出了：
listContour = Bezier.CreateCurveWithDistance(listContour, 0.6f, meshDensityTunning.Horizontal);
至此，我们获得了所有CT图像的骨骼外轮廓点，数据储存在变量List<List<Vector3>> ListListResult中。接下来，就可以进行三维建模了。