血管的三维重建
血管的三维重建
摘 要
断面可用于了解生物组织、器官等的形态，在医学上有重要的作用。用切片机连续不断地将样本切成数十、成百的平行切片，可依次逐片观察。根据平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。
本文提出了一种基于平行切片图像的血管三维重建方法。
针对问题一求每张切片的最大内切圆半径和圆心：我们利用图像处理技术提取每张切片图像的信息，每张切片所有内点到轮廓点的最小距离中的最大值即为最大内切圆的半径，求平均值得到血管的半径，R=29.4167。
针对问题二中轴线拟合：基于问题一，得到100张切片的圆心坐标，用最小二乘法拟合中轴线方程，综合SSE与R2可知6阶多项式相对最优。中轴线方程如下：
接着在matlab中用plot函数绘出中轴线在XY,XZ,YZ,平面的投影图。
针对问题三三维血管重建及重新切片：基于问题二，对血管进行三维重建，通过遍历像素空间内所有的像素点距离中轴线的最短距离找到处于血管壁的像素点，其内部的所有像素点均为血管内部结构，转换为图像即可得到复原的血管切片。再与原切片进行重合度实验，即将像素矩阵中对应像素点进行比对，所有切片的平均重合度为80.55%。
验证了模型的有效性。实验结果表明，所提出的方法能够较为精确地重建出血管的三维形态。
关键词：最大内切圆 最小二乘法 血管切片 MATLAB图像处理 重合度
一、问题的重述
假设某些血管可视为一类特殊的由球心沿着某一曲线（称为中轴线）的球滚动包络而成的管道。
现有某血管管道的连续100张平行切片图象，记录了管道与切片的交，格式均为BMP，宽、高均为512个象素。建立坐标系，坐标系的Z轴垂直于切片，第1张切片为平面Z=0，第100张切片为平面Z=99。Z=z切片图象中象素的坐标根据它们在文件中出现的前后次序依此为：
（-256，-256，z），（-256，-255，z），…（-256，255，z），
（-255，-256，z），（-255，-255，z），…（-255，255，z），
……
（ 255，-256，z），（ 255，-255，z），…（255，255，z）。
问题一：试计算此血管管道100张切片的圆心和半径，并计算出最后的血管管道半径，需给出针对此问题的具体的算法。
问题二：根据第一问的结果，绘制中轴线以及中轴线在xy,yz,zx
问题三：画出血管的三维重建图，以可视化血管的三维形态。利用上述所有过程得到的结果重新进行模拟切片，与已有的切片模型进行重合度比对，并评价建立的模型及结果的优劣。
二、模型的假设
1.血管的表面是由半径固定、球心变化且连续的一族球运动形成的包络面。
2.管道中轴线与每张切片有且只有一个交点。
3.医学上,血管不存在过于严重的扭曲。
4.血管壁的厚度忽略不计。
三、符号说明
| 符号 | 含义 |
| R | 球的半径 |
| (xk,yk,k) | 第k张切片的最大内切圆圆心坐标 |
| rk | 第k张切片的最大内切圆的半径 |
| P | 512*512*100的三维像素矩阵 |
| A | 某张切片的内点点集 |
| B | 某张切片的轮廓点点集 |
| c | 重合度 |
四、问题的分析
4.1.问题一
问题要求我们计算出管道的半径并给出具体的算法。血管的表面是由半径固定、球心变化且连续的一族球运动形成的包络面，球心轨迹为中轴线。因为所有切片都过中轴线且与之有且只有一个交点，所以所有的切片都过球心。则有结论：每个切片都含有球心，球心到轮廓线的最短距离即为球的半径，圆心即为中轴线上的点。
对于找此最短距离，先寻找每张切片上过中轴线的圆心，则需要计算切片的内点到所有边界点的距离，寻找最小值，再找寻所有内点最小距离中的最大值，将此值作为此切片所得半径，所得此内点即为球心，做出的圆即为最大内切圆。再对这100个切片半径求平均值即为球的半径。
4.2.问题二
由问题一的结果可以得到100组球心的点，对这100组坐标运用最小二乘法进行曲线拟合，在matlab中运用plot3函数，就可以得到中轴线的轨迹，在matlab中用plot函数绘制出中轴线在
平面的投影图。
4.3.问题三
获得中轴线方程后，通过连续绘制离散的球体，模拟滚动包络过程，即可绘制出血管的近似三维模型。
以原切片方式即垂直于Z轴进行100次平行切片得到的100个512*512的像素矩阵，将这100张切片的像素矩阵与原切片的像素矩阵比对，计算重合的像素点在整张切片中的占比。最后检验模型优劣。
五、模型建立与求解过程
5.1.问题一：最大内切圆半径模型的建立与求取
5.1.1.数据处理
在matlab中引用imread和imshow函数[2]（程序详见附录一）将图象文件名依次为0.bmp，1.bmp，......，99.bmp的宽、高均为 512 个象素的100个切片图像转换为512*512*100的三维0-1矩阵P。最后会有100片512*512的0-1像素矩阵，第k片的格式如下所示，
为矩阵中的元素，取值为0或者1，0代表切片中黑色的点，1代表切片中白色的点。对于这一百张切片，我们将其编号为1，2，3，...，100，使得k=1，2，...，100。
5.1.2.最大内切圆模型
确定R以中轴线上每一点为球心,以固定常数R为半径作球,产生一球面簇,血管即为该球面簇包络，半径R由每一个切片的血管部分的最大内切圆半径决定。
以第k张切片为例建立最大内切圆模型，为此，做以下说明：
- 在切片像素矩阵中若某一点周围的数值全为0，则此点为内点，内点集即为A；若某一点周围的数值有0也有1，则此点为轮廓点，轮廓点集即为B。
- 最大内切圆为切片中以血管部分的所有内点中每一个内点到轮廓点距离的最小值中的最大值为半径，再以这个内点为圆心做出的圆即为最大内切圆。[1]
图5-1-1 以第77张切片为例的最大内切圆示意图
得到表5-2-1的数据如下所示：
表5-1-1 100张切片的最大内切圆的圆心坐标及半径
| 最大内切圆的圆心坐标 | 最大内切圆半径 | 最大内切圆的圆心坐标 | 最大内切圆半径 | ||||
| x | y | Z+1 | x | y | Z+1 | ||
| -160 | 1 | 1 | 29.06888 | -114 | 117 | 51 | 29.69848 |
| -160 | 0 | 2 | 28.28427 | -114 | 117 | 52 | 29.69848 |
| -160 | 2 | 3 | 29.01724 | -113 | 118 | 53 | 29.69848 |
| -160 | 2 | 4 | 29.06888 | -112 | 119 | 54 | 29.69848 |
| -160 | 2 | 5 | 29.06888 | -111 | 120 | 55 | 29.68164 |
| -160 | 2 | 6 | 29.06888 | -111 | 120 | 56 | 29.20616 |
| -160 | 1 | 7 | 29 | -63 | 151 | 57 | 29.41088 |
| -160 | 4 | 8 | 29.01724 | -75 | 145 | 58 | 29.52965 |
| -160 | 1 | 9 | 29 | -81 | 142 | 59 | 29.52965 |
| -160 | 1 | 10 | 28.86174 | -51 | 156 | 60 | 29.54657 |
| -160 | 7 | 11 | 28.86174 | -51 | 156 | 61 | 29.54657 |
| -160 | 8 | 12 | 28.86174 | -31 | 162 | 62 | 29.61419 |
| -160 | 9 | 13 | 28.86174 | -31 | 162 | 63 | 29.61419 |
| -160 | 10 | 14 | 29.01724 | -31 | 162 | 64 | 29.61419 |
| -160 | 12 | 15 | 29.01724 | -35 | 161 | 65 | 29.61419 |
| -160 | 13 | 16 | 29.01724 | -35 | 161 | 66 | 29.61419 |
| -160 | 14 | 17 | 29.01724 | -26 | 163 | 67 | 29.42788 |
| -160 | 16 | 18 | 29.01724 | -35 | 161 | 68 | 29.41088 |
| -160 | 17 | 19 | 29.01724 | -26 | 163 | 69 | 29.27456 |
| -160 | 18 | 20 | 29.01724 | 46 | 163 | 70 | 29.42788 |
| -160 | 19 | 21 | 29.01724 | 46 | 163 | 71 | 29.61419 |
| -160 | 20 | 22 | 29.01724 | 46 | 163 | 72 | 29.61419 |
| -160 | 21 | 23 | 29.01724 | 46 | 163 | 73 | 29.61419 |
| -160 | 22 | 24 | 29.01724 | 65 | 158 | 74 | 29.61419 |
| -160 | 21 | 25 | 29.06888 | 68 | 157 | 75 | 29.73214 |
| -160 | 21 | 26 | 29.06888 | 65 | 158 | 76 | 29.73214 |
| -160 | 21 | 27 | 29.06888 | 81 | 152 | 77 | 29.54657 |
| -159 | 30 | 28 | 29.15476 | 81 | 152 | 78 | 29.52965 |
| -159 | 30 | 29 | 29.27456 | 81 | 152 | 79 | 29.52965 |
| -159 | 29 | 30 | 29.27456 | 135 | 118 | 80 | 29.41088 |
| -158 | 35 | 31 | 29.42788 | 136 | 117 | 81 | 29.69848 |
| -157 | 40 | 32 | 29.61419 | 136 | 117 | 82 | 29.69848 |
| -157 | 40 | 33 | 29.61419 | 137 | 116 | 83 | 29.42788 |
| -157 | 40 | 34 | 29.61419 | 138 | 115 | 84 | 29.42788 |
| -156 | 44 | 35 | 29.61419 | 138 | 115 | 85 | 29.42788 |
| -153 | 55 | 36 | 29.73214 | 139 | 114 | 86 | 29.42788 |
| -153 | 55 | 37 | 29.73214 | 139 | 114 | 87 | 29.42788 |
| -153 | 55 | 38 | 29.73214 | 139 | 114 | 88 | 29.42788 |
| -152 | 58 | 39 | 29.73214 | 140 | 113 | 89 | 29.42788 |
| -152 | 58 | 40 | 29.61419 | 140 | 113 | 90 | 29.42788 |
| -150 | 63 | 41 | 29.54657 | 172 | 67 | 91 | 29.42788 |
| -149 | 66 | 42 | 29.54657 | 172 | 67 | 92 | 29.42788 |
| -148 | 68 | 43 | 29.52965 | 172 | 67 | 93 | 29.42788 |
| -148 | 68 | 44 | 29.52965 | 172 | 67 | 94 | 29.42788 |
| -143 | 78 | 45 | 29.52965 | 182 | 43 | 95 | 29.42788 |
| -137 | 88 | 46 | 29.41088 | 187 | 24 | 96 | 29.42788 |
| -137 | 88 | 47 | 29.41088 | 187 | 24 | 97 | 29.42788 |
| -116 | 115 | 48 | 29.69848 | 187 | 24 | 98 | 29.42788 |
| -115 | 116 | 49 | 29.69848 | 187 | 24 | 99 | 29.42788 |
| -115 | 116 | 50 | 29.69848 | 188 | 18 | 100 | 29.42788 |
5.1.3.半径求取
将 100 个切片所对应的最大内切圆半径
求取平均值，得血管管道半径R，R的计算公式为：
，带入表5-1-1中的第四列和第八列数据得R=29.4167。
5.2.问题二
5.2.1.中轴线投影
分别通过引用表5-1-1中的坐标
在matlab中用plot函数进行最小二乘拟合。对此，我们进行了1~9阶的多项式拟合，最后分析可得在六阶时SSE与R2相对最优（拟合1~9阶的SSE与R2结果详见附录二），得到中轴线在XY、YZ、ZX平面的投影图。
5.2.2.中轴线拟合
在matlab中运用plot3函数拟合可得拟合图。用最小二乘法拟合中轴线，由问题一得到的数据（圆心坐标）如表5-1-1所示，对XZ、YZ两两拟合可得中轴线以z为参数的参数方程，多项式系数选取详见5.2.1。
中轴线拟合后的参数方程为：
拟合图为：
图5-2-1 中轴线拟合图
5.3.问题三
5.3.1.三维血管重建
基于题目假设，该血管管道是由一个球体滚动包络而成，而该球体的半径及球心轨迹线方程（中轴线方程）已经求得，我们将该球体沿着中轴线滚动绘图，即可获得拟合血管的三维重建图。
此模型下的三维重建图如下所示：
图5-3-1 血管重建三维示意图
5.3.2.模拟切片
由上述拟合中轴线方程组结果得出预测值，在此预测值的基础上重新进行模拟切片,在重新进行血管切片时，我们考虑到重建的血管模型是由一颗半径确定的球体沿着中轴线滚动形成的包络面，那么该血管具有的一项特征就是中轴线上的点与血管壁之间的最短距离应为所求的血管半径，那么就可以通过遍历该512*512*100的像素单位空间内所有的像素点距离中轴线的最短距离以寻找到处于血管壁的像素点，寻找到血管壁以后其内部的所有像素点均为血管内部结构，使用1表示白色，0表示黑色，将对应血管壁内外使用0、1填充，转换为图像结即可得到复原的血管切片。
部分切片结果如下所示（切片结果详见支撑材料的A02bmp）
5.3.3.重合度检验
将新的100张切片在matlab中引用imread和imshow函数转换为100个512*512的01像素矩阵，0代表黑色的点，1代表白色的点。编写程序计算每组切片的重合度，重合度的计算公式为：
c=重合像素点/全部像素点
在编写程序时，将两个512*512的像素矩阵相加，将像素点对调，即0对应切片中的白点，1对应切片中的黑点。
在实际考察中，我们发现黑白像素点的数量差较大，这种情况会一定程度上影响重合度的计算，所以我们从算法上考虑将在切片中无关白色像素点进行消除。
改进方法为将两组血管切片的二值化矩阵相加，产生一组新的矩阵，这组矩阵由0、1、2组成，其中：
当像素点位值为2时，表示原血管切片与重建血管切片的对应位置均为黑色，是构成血管模型的点位，该点有效，记该点数量为cnt；
当像素点位值为1时，表示原血管切片与重建血管切片的对应位置黑色与白色各有一项，这意味着重建的血管模型与原血管模型存在误差，是模型的误差点，该点有效,记该点数量为cnt1+cnt2；
当像素点位值为0时，表示原血管切片与重建血管切片的对应位置均为白色，是不影响血管模型的无关点位，不考虑。
则重合度计算公式变为：c=cnt/(cnt1+cnt2-cnt)
最后得到整体重合度高达80.55%，其中切片编号为21.bmp~83.bmp的切片重合度高于80%，编号为0.bmp~20.bmp、84.bmp~91.bmp的切片重合度高于70%，编号为92.bmp~97.bmp的切片重合度在60%~70%之间，98.bmp和99.bmp两张图片的重合度低于60%。由此可知切片重合度有高有低，最高的重合度可达88.16%，最差的重合度只有56.02%。
表5-3-1 100张切片重合度表
| 重合度 | 切片编号 | 重合度 | 切片编号 | 重合度 | 切片编号 | 重合度 | 切片编号 |
| 72.39% | 0 | 81.94% | 25 | 86.81% | 50 | 86.90% | 75 |
| 71.39% | 1 | 82.42% | 26 | 86.90% | 51 | 86.31% | 76 |
| 71.29% | 2 | 82.58% | 27 | 87.00% | 52 | 85.82% | 77 |
| 71.65% | 3 | 82.70% | 28 | 87.01% | 53 | 85.07% | 78 |
| 71.94% | 4 | 82.49% | 29 | 87.01% | 54 | 84.30% | 79 |
| 72.33% | 5 | 82.48% | 30 | 87.09% | 55 | 83.50% | 80 |
| 72.65% | 6 | 82.68% | 31 | 87.08% | 56 | 82.50% | 81 |
| 73.23% | 7 | 82.96% | 32 | 87.13% | 57 | 81.45% | 82 |
| 73.64% | 8 | 83.44% | 33 | 87.09% | 58 | 80.48% | 83 |
| 74.10% | 9 | 83.66% | 34 | 87.10% | 59 | 79.46% | 84 |
| 74.76% | 10 | 83.87% | 35 | 87.07% | 60 | 78.26% | 85 |
| 75.41% | 11 | 84.24% | 36 | 87.11% | 61 | 77.09% | 86 |
| 75.97% | 12 | 84.68% | 37 | 87.23% | 62 | 75.84% | 87 |
| 76.26% | 13 | 84.88% | 38 | 87.27% | 63 | 74.65% | 88 |
| 76.96% | 14 | 85.20% | 39 | 87.38% | 64 | 73.33% | 89 |
| 77.59% | 15 | 85.38% | 40 | 87.40% | 65 | 71.88% | 90 |
| 77.93% | 16 | 85.83% | 41 | 87.53% | 66 | 70.42% | 91 |
| 78.43% | 17 | 85.91% | 42 | 87.65% | 67 | 68.98% | 92 |
| 78.71% | 18 | 86.10% | 43 | 87.91% | 68 | 67.37% | 93 |
| 79.31% | 19 | 86.23% | 44 | 88.06% | 69 | 65.70% | 94 |
| 79.57% | 20 | 86.35% | 45 | 88.16% | 70 | 63.91% | 95 |
| 80.04% | 21 | 86.46% | 46 | 88.16% | 71 | 62.02% | 96 |
| 80.52% | 22 | 86.61% | 47 | 88.06% | 72 | 60.13% | 97 |
| 81.20% | 23 | 86.65% | 48 | 87.81% | 73 | 58.13% | 98 |
| 81.47% | 24 | 86.76% | 49 | 87.46% | 74 | 56.02% | 99 |
显然，此模型在血管的前段有较好的拟合效果，对于血管的尾段拟合效果较差。
六、模型的评价、改进与推广
6.1.模型的优缺点
6.1.1.模型的优点
1）本模型较好地利用到题设给出的小球滚动包络形成血管这一特点，使用MATLAB中的imread函数、imshow函数、edge函数等对题目给出的血管切片图进行二值化、轮廓化，通过计算内点至轮廓的最短距离，即用最大覆盖法求解最大内切圆，极大地简化了寻找最大内切圆圆心和半径的过程。
2）在重合度计算时，应用矩阵相加后的特点，舍去无关的白色区域，增加了计算结果的精确度。
6.1.2.模型的缺点
1）模型只适应于一类特殊的管道——该管道的表面是由球心沿着某一曲线（中轴线）的球滚动包络而成。对于人体错综复杂的血管来说，适用性太低。
2）模型在重新切片后与原切片的重合度较低，对于血管模型末端的拟合较差，三维重建有失偏颇。
3）血管的重建模型虽然直观地表现出来了，但是并没能通过一组具体的数学方程进行表示，实际上血管模型仍然是离散的而非连续的，如何找到合适的方法基于中轴线和滚动球体建立曲面方程，这是有待学习的。
6.2.模型的改进
1）方法优化：
在用最小二乘法拟合时，次数越高拟合效果越好。亦或者使用稳健回归，找出那些对模型影响最大的样本点，当最小二乘法遇到此种数据样本点，去掉这些异常样本点重新进行拟合。
2）拟合方法：
在中轴线投影在XY，YZ，ZX上时，用的还是问题一得到的坐标进行拟合，而非用的是拟合中轴线后新的数据点。这些数据点中会有误差较大的数据点，从而影响结果的准确性。
七、参考文献
[1]陈凌钧 ,骆岩林. 等径管道的三维重建[J] . 高校应用数学学报. Vol. 13 Ser A Suppl. 1998 ,87 - 90
[2]《Visual C++数字图象处理》第12页2.3.1节。何斌等编著，人民邮电出版社，2001年4月。
八、附录
附录一（代码）
图像处理
clc;
clear all;
image = dir('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\');% 路径为原切片存储路径
files = dir(fullfile('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\','*.bmp'));% 处理的图片格式为bmp
SumFiles = length(files);
A=cell(SumFiles,1);% 用cell来存储每个图片所对应的矩阵
for i = 1:SumFiles
Img = imread(strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\',files(i).name));%文件所在路径
A{i,1}=Img;
disp(strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\',files(i).name)) %打印文件路径
imshow(Img)
end
for tt=1:SumFiles
if length(size(A{tt,1}))==3
A{tt,1}=A{tt}(:,:,1);
end
end
求半径
anss = zeros(101, 4); %答案矩阵
%进行灰度反转
for i = 0:99
k = strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\',int2str(i), '.bmp');
A = imread(k);
for j = 1:512
for k = 1:512
A(j,k) = 1 - A(j,k);
end
end
B = edge(A, 'sobel');% 边缘化处理
B1 = bwmorph(A, 'skel', inf);%骨架化处理
%imshow(B1);
[x, y, z] = find(B); %边缘点坐标
[a, b, c] = find(B1); % 骨架点坐标
m = length(x); %边缘点的数量
n = length(a); %骨架点的数量
temp = zeros(n, m); %距离矩阵
idx = zeros(n, 2); %存储最小值结果
for p = 1: n
for q = 1: m
x0 = a(p);
y0 = b(p);
x1 = x(q);
y1 = y(q);
temp(p,q) = sqrt((x0 - x1)^2 + (y0 - y1)^2);
end
[idx1, idx2] = min(temp(p,:));
idx(p, 1) = idx1; %最小值
idx(p, 2) = idx2; %下标
end
[maxn, idx3] = max(idx(:,1));
x = a(idx3) - 256;%相对屏幕中心
y = b(idx3) - 256;
%存储信息
anss(i+1, 1) = x;
anss(i+1, 2) = y;
anss(i+1, 3) = i + 1;
anss(i+1, 4) = maxn;
end
sum = 0;
for i = 1:100
sum = sum + anss(i,4);
end;
r = sum/100;
中轴线拟合与投影
xx = anss(:,1);
yy = anss(:,2);
zz = anss(:,3);
format long
px = [ 1.199e-08 -3.797e-06 0.0004276 -0.02049 0.443 -3.835 -151.6 ];
py = [2.729e-09 -2.8e-07 -4.688e-05 0.006523 -0.2036 2.936 -7.126 ];
x3 = polyval(px, zz);
y3 = polyval(py, zz);
figure(1);
plot3(x3, y3, zz)
grid on
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
title('中轴线图')
figure(2)
plot(zz, x3, '-');
ylabel('X轴');
xlabel('Z轴');
title('中轴线XOZ平面投影图');
grid on
figure(3)
plot(zz, y3, '-');
ylabel('Y轴');
xlabel('Z轴');
title('中轴线YOZ平面投影图');
grid on
figure(4)
plot(xx, yy, '-');
ylabel('Y轴');
xlabel('X轴');
title('中轴线XOY平面投影图');
grid on
三维重建
x = anss(:,1);
y = anss(:,2);
z = anss(:,3);
plot3(x,y,z,'+');
hold on
px = [ 1.199e-08 -3.797e-06 0.0004276 -0.02049 0.443 -3.835 -151.6 ];
py = [2.729e-09 -2.8e-07 -4.688e-05 0.006523 -0.2036 2.936 -7.126 ];
x1=polyval(px,z);%由参数方程预测出来的值
y1=polyval(py,z);
ball=zeros(100,3);%球体
ball(:,1)=x1;
ball(:,2)=y1;
ball(:,3)=z;
t=linspace(0,pi,25);%球的极角
p=linspace(0,2*pi,25);%球的方向角
[theta,phi]=meshgrid(t,p);%二维球面
for i=1:100
%球的参数方程形式
x=r*sin(theta).*sin(phi)+ball(i,1);
y=r*sin(theta).*cos(phi)+ball(i,2);
z=r*cos(theta)+ball(i,3);
hold on
surf(x,y,z);
plot3(x,y,z,'red')
axis equal; %使坐标轴按照单位长度分割开
end
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
title('血管拟合后还原图');
hold off
切片
R = 29.4167; % 血管的固定半径
z = 0:100; % Z轴坐标范围，对应100张切片
% 中轴线在x方向上的函数关系式
x = 1.199e-08 .* (z-1).^6 - 3.797e-06 .* (z-1).^5 + 0.0004276 .* (z-1).^4 - 0.02049 .* (z-1).^3 + 0.443 .* (z-1).^2 - 3.835 .* (z-1) - 151.6;
% 中轴线在y方向上的函数关系式
y = 2.729e-09 .* (z-1).^6 - 2.8e-07 .* (z-1).^5 - 4.688e-05 .* (z-1).^4 + 0.006523 .* (z-1).^3 - 0.2036 .* (z-1).^2 + 2.936 .* (z-1) - 7.126;
for m = 0:99 % m代表当前处理的切片层数，Z轴上的高度
a = ones(512, 512); % 初始化当前切片图像为全白（灰度值为1）
for i = 1:512 % 遍历图像的行
for j = 1:512 % 遍历图像的列
for p = 1:100 % 遍历血管中轴线上100个点的Z轴坐标
% 计算图像中(i, j)像素点到血管中轴线上每个点的距离
b(p) = sqrt(((i-257) - x(p))^2 + ((j-257) - y(p))^2 + (m+1 - p)^2);
end
% 如果像素点到中轴线的最近距离小于等于血管半径，则该点在血管内部，标记为黑色
if min(b) <= R
a(i, j) = 0; % 将该像素点设置为黑色（灰度值为0）
end
end
end
% 保存当前切片图像到指定路径
imwrite(a, strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A02bmp\', int2str(m), '.bmp'));
end
重合度检测
anss = zeros(100, 2);
for i = 0:99
% 构建图像路径
path = strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A01bmp\', int2str(i), '.bmp');
path1 = strcat('C:\Users\ASUS\Desktop\第一次训练\A题 血管的三维重建\A02bmp\', int2str(i), '.bmp');
% 显示正在处理的文件路径
disp(['Processing files: ', path, ' and ', path1]);
% 读取图像
A = imread(path);
B = imread(path1);
imshow(B);
% 灰度反转并转换为 double 类型
A = 1 - double(A);
B = 1 - double(B);
% 计算图像 A 和 B 的和
C = A + B;
% 初始化计数器
cnt = 0;
cnt1 = 0;
cnt2 = 0;
% 遍历图像矩阵，计算所需的计数
for j = 1:512
for k = 1:512
if C(j, k) == 2
cnt = cnt + 1;
end
if A(j, k) == 1
cnt1 = cnt1 + 1;
end
if B(j, k) == 1
cnt2 = cnt2 + 1;
end
end
end
% 计算重叠区域的比例 p
p = cnt / (cnt1 + cnt2 - cnt);
% 存储结果
anss(i+1, 1) = p;
anss(i+1, 2) = i;
% 显示计算结果
disp(['Slice ', int2str(i), ': p = ', num2str(p)]);
end
附录二（投影拟合SSE与R2数据表）
| 拟合次数 | yz | xy | xz | |||
| SSE | R2 | SSE | R2 | SSE | R2 | |
| 1 | 198700 | 0.4068 | 272500 | 0.1865 | 244500 | 0.8495 |
| 2 | 95320 | 0.7154 | 20820 | 0.9379 | 45220 | 0.9722 |
| 3 | 10720 | 0.968 | 18950 | 0.9434 | 35520 | 0.9775 |
| 4 | 10380 | 0.969 | 7093 | 0.9788 | 12780 | 0.9921 |
| 5 | 6122 | 0.9817 | 6303 | 0.9812 | 12390 | 0.9944 |
| 6 | 6055 | 0.9819 | 3962 | 9.9882 | 11110 | 0.9932 |
| 7 | 5437 | 0.9837 | 3407 | 0.9898 | 10590 | 0.9935 |
| 8 | 5451 | 0.9837 | 2819 | 0.9916 | 10560 | 0.9935 |
| 9 | 5127 | 0.9847 | 2619 | 0.9922 | 10470 | 0.9936 |