不规则曲面上两点距离求取

 背景

在CT中求皮肤上两点间的弧长。由于人体表面并不是规则的曲面，不可能用圆的弧长求取方法来计算出两点间的弧长。

而在不规则的曲面上求两点的距离，都可以用类似测地线距离求取的方式来求取（积分），而转化为搜索路径问题。

弧长计算思路

以下在轴位CT的皮肤表面的两点间弧长计算思路。

CT图像是以一个个像素组成的一个规则的图像，通常接触到的都是512*512分辨率的CT，像素间距（一般像素的x和y轴间距相同，若不同需要特殊考虑，本文仅考虑x和y轴像素间距相同的情况）为零点几mm到5mm左右。

那么也就是说，其实CT图像就是一个512*512的网格，每个网格的大小即像素间距d毫米。那么计算两个点在不规则皮肤上的距离，即可找到点A到点B的前进方向，然后沿着前进方向在皮肤与空气的交界处一个像素一个像素的前进到B点，然后按若相邻点是斜相邻则距离为√2个像素间距，竖直或水平相邻则为1个像素间距，最后将累加起来的像素距离转化为真实距离。

A在左下，B在右上

具体则如下，若A点为在左下角，B点在右上角，那以A点为起点。

1.判断A点的左上角像素点是否皮肤，是，则将距离加上√2，然后A前进到左上像素点；

2.若不是，则判断上方像素点是否是皮肤，是，则将距离加上1；

3.若不是，则判断右上方像素点是否是皮肤，是，则将距离加上√2；

4.若不是，则判断右方像素点是否是皮肤，是，则将距离加上1，否则A、B点放置是有问题的。在A、B放置正确的情况下，若前述像素点都不是皮肤，那么最后一个点肯定是皮肤（A、B都是皮肤上的一个点，那那一定能从A点到B点）。

A前进到下一个点后，然后再重复上述4个步骤，一直到到最后A点的x与y坐标等于B点的x与y坐标（由于是在轴位上，故不考虑Z轴坐标）。

注意：最后比较时x、y的值需要转为整数，因为位置移动是按一个一个像素来移动的。

如下示意图所示，A点为棕色点，B点为蓝色点；黑色为空气，黄色为皮肤。

 A在左上，B在右下

若A点在左上角，B点在右下角，那么A到B，那么A点移动时的逻辑，则是应该判断A右方的像素点，A右下角像素，A点下方像素，A点左下方像素。

具体实现

代码如下

internal double ObliqueSearchSkinPoint(Vector3[] directions,Vector3 startPoint, Vector3 endPoint)
{var stack = new Stack<(Vector3 Point, double Distance)>();var visited = new HashSet<Vector3>();stack.Push((startPoint, 0));visited.Add(startPoint);while (stack.Count > 0){var (current, distance) = stack.Pop();if ((int)current.x == (int)endPoint.x && (int)current.y == (int)endPoint.y){return distance; // 找到终点，返回累计距离}// 按方向顺序检查foreach (var dir in directions){Vector3 next = current + dir;if (!visited.Contains(next) && CheckPointOnSkin((int)next.z, (int)next.x, (int)next.y)){visited.Add(next);double cost = (dir.x != 0 && dir.y != 0) ? Math.Sqrt(2) : 1; // 斜向距离√2，直线距离1stack.Push((next, distance + cost));break; // 找到第一个可行点就跳出循环，不再检查其他方向}}}return double.PositiveInfinity; // 无可行路径,点标记错误
}

上述代码中添加了一个HashSet来记录访问过的点，这个后续发现在当前的场景下，其实是可以不需要的。 

以上CheckPointOnSkin方法的用途是用于判断像素点是否在皮肤上，主要是思路是皮肤上的像素点的CT值与空气的CT值不同，皮肤的CT是100左右的，而空气是-1000左右。

此仅实现了轴位CT图像皮肤上两点在同一张CT上的弧长计算；

若两点不在相同一张CT上的话，思路基本与此一致，不过需要考虑Z轴的变化，搜索前进方向则需要根据AB的方向向量来决定，搜索方向的数组就相对更复杂一些，实现当然也会复杂一些。以后再做讨论。

搜索方向

A点在左下，B点在右上时

   // 定义搜索方向（按优先级顺序）Vector3[] directions = {new Vector3(-1, -1, 0),  // 左上（斜向，距离√2）new Vector3(0, -1, 0),   // 上（直线，距离1）new Vector3(1, -1, 0),   // 右上（斜向，距离√2）new Vector3(1, 0, 0)    // 右（直线，距离1）};

注意：左上为什么y轴是-1，不是1；是因为在Canvas坐标系中，原点在左上角，x轴为从左到事，y轴为从上到下。

A点在左上，B点右下时

   // 定义搜索方向（按优先级顺序）Vector3[] directions = {new Vector3(1, 0, 0),  // 右（斜向，距离1）new Vector3(1, 1, 0),   // 右下（直线，距离√2）new Vector3(0, 1, 0),   // 下（斜向，距离1）new Vector3(-1, 1, 0)    // 右下（直线，距离√2）};

当然上述ObliqueSearchSkinPoint方法也可以递归来实现，不过递归时一定要处理好结束条件，否则极易造成溢出问题，

同时递归一定要注意可能的性能问题，主要是重复计算要尽量避免。

优化

同时上述是单方向的从A搜索向B，若需要加速搜索，是可以同时A向B与B向A双向搜索的，直到最后新的A与B重合即可。