血管的三维重建问题:
断面可用于了解生物组织、器官等的形态。例如,将样本染色后切成厚约1m m的切片,在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片, 可依次逐片观察。根据拍照并采样得到的平行切片数字图象,运用计算机可重建组织、器官等准确的三维形态。
假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。例如圆柱就是这样一种管道,其中轴线为直线,由半径固定的球滚动包络形成。
现有某管道的相继100张平行切片图象,记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、…、 99.bmp,格式均为BMP,宽、高均为512个象素(pixel)。为简化起见,假设:管道中轴线与每张切片有且只有一个交点;球半径固定;切片间距以及图象象素的尺寸均为1。
取坐标系的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)。
试计算管道的中轴线与半径,给出具体的算法,并绘制中轴线在XY、YZ、ZX平面的投影图。
第2页是100张平行切片图象中的6张,全部图象请从网上下载。(链接:http://pan.baidu.com/s/1jIazlFg)
关于BMP图象格式可参考:
1. 《Visual C++数字图象处理》第12页2.3.1节。何斌等编著,人民邮电出版社,2001年4月。
2. 《matlab图象处理》
血管的三维重建数学建模:
血管的三维重建
摘要
医学图像的三维重建涉及数字图像处理、计算机图形学以及医学领域的相关知识,是一个多学科交叉的研究领域,也是当前的一个热点研究问题。[1] 本文采用最大覆盖法、曲线拟合法对血管三维重建问题进行了分析和求解。
针对问题一:采用最大覆盖法[2]遍历每个切片的所有离散内点,得到每个内点到轮廓上所有点的最小距离,在这些最小距离中取最大值,存储最大值对应的内点即为每个切片的最大内切圆的圆心位置,最大值即为最大内切圆的半径。对100个切片的半径取平均值,得到球的半径和球心坐标,血管的半径为:
针对问题二:基于问题一已经得到100个球心坐标,利用MATLAB中的ployfit函数对这些球心坐标进行多项式拟合得到中轴线。经计算残差较小,所以选取五次多项式拟合,中轴线方程为:
针对问题三:在问题二中,已经得到了中轴线方程和中轴线的拟合曲线,利用MATLAB中的plot函数分别得到中轴线在三个坐标平面的投影图。
解决了上述三个问题后,接着根据所求的中轴线和球的半径,利用plot函数对血管进行三维重建,得到了血管的三维重建图形。让三维重建图形和由切片得到的血管三维图形对比重合度,进行模型检验。
关键词: 最大覆盖法 拟合曲线 投影图 三维重建
1 问题重述
1.1 问题背景
医学图像三维重建是目前研究的热点之一,是一个跨学科的研究领域,它将临床上通过医学影像技术得到的二维切片序列,利用计算机技术及图形图像学技术,将二维切片代表的人体组织以立体的方式展现出来,医生可以根据重建后的模型,进行旋转、缩放、移动等操作、进行多角度多层次的观察,甚至对病灶的大小、位置、形状等能得到准确的直观的认识,能获取比二维切片更多的信息,从而做出准确的诊断。[3]如要将人体全部三维信息,包含内部错综复杂的结构,完整地存储在计算机中,以现在的技术也是有一定难度的,但若改用存储人体切片信息,使用时重建再现的方法,则是利用现有技术可以解决的。假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。
1.2 所求问题
现有某管道的相继100张平行切片图象,记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、…、 99.bmp,格式均为BMP,宽、高均为512个象素(pixel)。假设:管道中轴线与每张切片有且只有一个交点;球半径固定;切片间距以及图象象素的尺寸均为1。取坐标系的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)。
问题一:计算管道的半径并给出具体的算法;
问题二:计算管道的中轴线并给出具体的算法;
问题三:绘制中轴线在X-Y、Y-Z、Z-X平面的投影图。
2 问题分析
问题一要求计算出管道的半径并给出具体的算法。因为球的任意截面都是圆,过球心的截面为最大圆,且每个切片与中轴线有且只有一个交点,所以切片必包含球的大圆。由此得到一个重要结论:切片的最大内切圆即为球的大圆,最大内切圆的半径即为球的半径,最大内切圆的圆心位置即为球心的位置。因此问题转化为求每一个切片的最大内切圆的半径和圆心。
图一:血管示意图
基本方法有两个:第一种方法为切线法,由于在最大内切圆与切片的两个切点所做切线平行,所以可以遍历切片边缘的所有平行切线,找到距离最大的一组平行切线,连结两个切点即为最大内切圆的直径。第二种方法为最大覆盖法,求任意一个内点到所有轮廓点的最小距离,再取所有内点对应的最小距离的最大值,此最大值即最大内切圆的半径,对应的内点即为最大内切圆的圆心。显然,后者在算法上更为优化,计算机在执行的过程中更为迅速。100个切片会产生100个半径的值,可以根据统计学中的平均值法,对100个数据进行平均,即可得到球的半径。
问题二要求计算出管道的中轴线并给出具体的算法。任意一条曲线都是无穷多个点的集合。当已知曲线上的有限个点时,可以采用描点作图的方法。基于问题一已经得到100个切片的最大内切圆的圆心位置,对100个圆心位置进行曲线拟合,所得曲线即为中轴线的拟合曲线。
问题三要求绘制出中轴线在X-Y、Y-Z、Z-X平面的投影图。基于问题二中已经得到中轴线的拟合曲线,在MATLAB软件中,利用plot函数即可做出其在三个坐标平面的投影图。
3 模型假设
1.假设血管无严重扭曲;
2.假设管道中轴线与每张切片有且只有一个交点;
3.假设球半径固定;
4.假设切片间距以及图象象素的尺寸均为1;
5.假设切片拍摄不存在误差,数据误差仅与切片数字图像的分辨率有关。
4 符号说明
5 模型的建立与求解
5.1.问题一
题中给出了某血管的相继100张平行切片bmp图像,记录了管道与切片的交,要求计算出管道的半径并给出具体的算法。
5.1.1模型的建立
(1)存储数据,转换格式
调用imread、imshow函数将100张切片的bmp图像录入到MATLAB中,并将其转化为512*512*100的三维0-1矩阵,其中1代表黑色像素点,0代表白色像素点。在每一个512*512像素的图像中,每一个像素都有一个确定的坐标,可以找出这100张图片的像素矩阵;
(2)求切片的轮廓
调用MATLAB中的内部函数edge可以得到所有切片的轮廓线;
(3)求出每个切片最大内切圆的半径
对切片内部任意一个点,求出它到轮廓线上所有点的距离,并取其最小值,由于所有内点都对应一个最小值,在这些最小值中取最大值,即为:最大内切圆的半径;
假设切片内部有p个点,第i个点的坐标为 (xi,yi),切片轮廓线上有q个点,第j个点的坐标为(mj,nj),则切片内部第i点到轮廓线上第j点的距离为:
第i个点对应的最小距离ri 为:
(4)将100个切片所对应的最大内切圆半径求取平均值,即为:血管的半径。
5.1.2模型求解
首先在MATLAB中利用imread函数和imshow函数画出切片的轮廓点(具体程序见附录一),接着调用内部函数edge得到所有切片的轮廓线,求出了切片内部的点到轮廓上点的最小距离,在这些最小值中取最大值,即为最大内切圆的半径。(具体程序见附录二)100张切片最大内切圆的圆心坐标及半径如表一所示:
表一:100张切片最大内切圆的圆心坐标及半径
切片序号 | 最大内切圆圆心的坐标 | 最大内切圆的半径/μm | ||||
x/μm | y/μm | z/μm | ||||
0 | -160 | 1 | 1 | 29.0688 | ||
1 | -160 | 0 | 2 | 28.2842 | ||
2 | -160 | 2 | 3 | 29.0172 | ||
3 | -160 | 2 | 4 | 29.0688 | ||
4 | -160 | 2 | 5 | 29.0688 | ||
5 | -160 | 2 | 6 | 29.0688 | ||
6 | -160 | 1 | 7 | 29.0000 | ||
7 | -160 | 4 | 8 | 29.0172 | ||
8 | -160 | 1 | 9 | 29.0000 | ||
9 | -160 | 1 | 10 | 28.8617 | ||
10 | -160 | 7 | 11 | 28.8617 | ||
11 | -160 | 8 | 12 | 28.8617 | ||
12 | -160 | 9 | 13 | 28.8617 | ||
13 | -160 | 10 | 14 | 29.0172 | ||
14 | -160 | 12 | 15 | 29.0172 | ||
15 | -160 | 13 | 16 | 29.0172 | ||
16 | -160 | 14 | 17 | 29.0172 | ||
17 | -160 | 16 | 18 | 29.0172 | ||
18 | -160 | 17 | 19 | 29.0172 | ||
19 | -160 | 18 | 20 | 29.0172 | ||
20 | -160 | 19 | 21 | 29.0172 | ||
21 | -160 | 20 | 22 | 29.0172 | ||
22 | -160 | 21 | 23 | 29.0172 | ||
23 | -160 | 22 | 24 | 29.0172 | ||
24 | -160 | 21 | 25 | 29.0688 | ||
25 | -160 | 21 | 26 | 29.0688 | ||
26 | -160 | 21 | 27 | 29.0688 | ||
27 | -159 | 30 | 28 | 29.1547 | ||
28 | -159 | 30 | 29 | 29.2745 | ||
29 | -159 | 29 | 30 | 29.2745 | ||
30 | -158 | 35 | 31 | 29.4278 | ||
31 | -157 | 40 | 32 | 29.6141 | ||
32 | -157 | 40 | 33 | 29.6141 | ||
33 | -157 | 40 | 34 | 29.6141 | ||
34 | -156 | 44 | 35 | 29.6141 | ||
35 | -153 | 55 | 36 | 29.7321 | ||
36 | -153 | 55 | 37 | 29.7321 | ||
37 | -153 | 55 | 38 | 29.7321 | ||
38 | -152 | 58 | 39 | 29.7321 | ||
39 | -152 | 58 | 40 | 29.6141 | ||
40 | -150 | 63 | 41 | 29.5465 | ||
41 | -149 | 66 | 42 | 29.5465 | ||
42 | -148 | 68 | 43 | 29.5296 | ||
43 | -148 | 68 | 44 | 29.5296 | ||
44 | -143 | 78 | 45 | 29.5296 | ||
45 | -137 | 88 | 46 | 29.4108 | ||
46 | -137 | 88 | 47 | 29.4108 | ||
47 | -116 | 115 | 48 | 29.6984 | ||
48 | -115 | 116 | 49 | 29.6984 | ||
49 | -115 | 116 | 50 | 29.6984 | ||
50 | -114 | 117 | 51 | 29.6984 | ||
51 | -114 | 117 | 52 | 29.6984 | ||
52 | -113 | 118 | 53 | 29.6984 | ||
53 | -112 | 119 | 54 | 29.6984 | ||
54 | -111 | 120 | 55 | 29.6816 | ||
55 | -111 | 120 | 56 | 29.2061 | ||
56 | -63 | 151 | 57 | 29.4108 | ||
57 | -75 | 145 | 58 | 29.5296 | ||
58 | -81 | 142 | 59 | 29.5296 | ||
59 | -51 | 156 | 60 | 29.5465 | ||
60 | -51 | 156 | 61 | 29.5465 | ||
61 | -31 | 162 | 62 | 29.6141 | ||
62 | -31 | 162 | 63 | 29.6141 | ||
63 | -31 | 162 | 64 | 29.6141 | ||
64 | -35 | 161 | 65 | 29.6141 | ||
65 | -35 | 161 | 66 | 29.6141 | ||
66 | -26 | 163 | 67 | 29.4278 | ||
67 | -35 | 161 | 68 | 29.4108 | ||
68 | -26 | 163 | 69 | 29.2745 | ||
69 | 46 | 163 | 70 | 29.4278 | ||
70 | 46 | 163 | 71 | 29.6141 | ||
71 | 46 | 163 | 72 | 29.6141 | ||
72 | 46 | 163 | 73 | 29.6141 | ||
73 | 65 | 158 | 74 | 29.6141 | ||
74 | 68 | 157 | 75 | 29.7321 | ||
75 | 65 | 158 | 76 | 29.7321 | ||
76 | 81 | 152 | 77 | 29.5465 | ||
77 | 81 | 152 | 78 | 29.5296 | ||
78 | 81 | 152 | 79 | 29.5296 | ||
79 | 135 | 118 | 80 | 29.4108 | ||
80 | 136 | 117 | 81 | 29.6984 | ||
81 | 136 | 117 | 82 | 29.6984 | ||
82 | 137 | 116 | 83 | 29.6984 | ||
83 | 138 | 115 | 84 | 29.6984 | ||
84 | 138 | 115 | 85 | 29.6984 | ||
85 | 139 | 114 | 86 | 29.6984 | ||
86 | 139 | 114 | 87 | 29.6984 | ||
87 | 139 | 114 | 88 | 29.6984 | ||
88 | 140 | 113 | 89 | 29.6984 | ||
89 | 140 | 113 | 90 | 29.6816 | ||
90 | 172 | 67 | 91 | 29.5296 | ||
91 | 172 | 67 | 92 | 29.5296 | ||
92 | 172 | 67 | 93 | 29.5296 | ||
93 | 172 | 67 | 94 | 29.5296 | ||
94 | 182 | 43 | 95 | 29.7321 | ||
95 | 187 | 24 | 96 | 29.6141 | ||
96 | 187 | 24 | 97 | 29.6141 | ||
97 | 187 | 24 | 98 | 29.6141 | ||
98 | 187 | 24 | 99 | 29.6141 | ||
99 | 188 | 18 | 100 | 29.4278 |
即血管的半径为:
5.2 问题二
题中给出了某血管的相继100张平行切片bmp图像,记录了管道与切片的交,要求计算出管道的中轴线并给出具体的算法。
5.2.1模型的建立
在问题一中,已经求出了每一个切片的最大内切圆半径,每一个最大内切圆有一个圆心,每一个切片对应一个圆心,100个切片就有100个圆心点。将100个圆心点在MATLAB中进行曲线拟合,得到的曲线即为中轴线。
具体算法如下:
(1)由表一可得:第k个切片的最大内切圆圆心的坐标为(ak,bk,ck),其中:k=1,2,3...,100;
(2)在MATLAB中,由于误差较小,可利用polyfit函数对100个圆心点做五次多项式拟合,利用plot函数得到的曲线即为中轴线的拟合曲线。
5.2.2模型求解
首先利用polyfit函数对100个圆心点做五次多项式拟合,接着利用plot函数画出中轴线的拟合曲线(具体程序见附录三),如图一所示。
5.3 问题三
题中给出了某血管的相继100张平行切片bmp图像,记录了管道与切片的交,要求绘制出中轴线在X-Y、Y-Z、Z-X平面的投影图。
5.3.1模型的建立
在问题二中,已得到了中轴线的拟合曲线和中轴线方程。可以利用plot函数分别做出其在三个坐标平面的投影图。
5.3.2模型求解
根据中轴线方程和中轴线的拟合曲线,在MATLAB中,利用plot函数分别做出其在三个坐标平面的投影图。(绘图程序见附录三)中轴线在X-Y、Z-X、Y-Z平面的投影图分别如图二、图三、图四所示。
5.4 模型检验
利用plot3函数可以将切片.bmp图像画为血管三维图像,(绘图程序见附录四),如图五所示。
对血管的三维重建图进行重新切片,与题目中所给的切片进行对比,计算出重合度,即可完成对血管三维重建模型的检验。不难发现,血管的三维重建图形较好地描述了血管的三维形态,而在一些区域,重建图和由切片得到的血管三维形态重合度较低,这表明模型存在着一定的误差。
6 模型分析
6.1模型优点
(1)解决血管三维重建问题时,应用了MATLAB数学软件对图片进行了处理,得出了大量数据;
(2)采用平均法对数据进行了科学精确地处理,保证了数据计算的精准度;
(3)本模型的算法采用遍历搜索的方法得到最大内切圆的半径的精度较高。
6.2模型缺点
(1)题目本身的像素所决定的误差不可避免;
(2)由于遍历搜索的运算次数异常庞大,导致模型的缺点是实践难度较大,整个运算耗费时间超过30分钟。
参考文献
[1]孙皓,医学图像三维重建算法研究, 上海交通大学,2005年1月
[2]丁峰平,工程数学学报,浙江工业大学,2002年
[3]董云飙,医学图像的三维重建研究,哈尔滨工业大学,2012年
附录
附录一:
画出血管的轮廓点
%构造两个零元素数组,每个数组里面有100个元素
p=zeros(512,512,100);
p2=zeros(512,512,100);
for i=0:99
%循环100次,拿到这100张图片,并读取到变量P里面
s=sprintf('D:\\program Files\\MATLAB\\R2011a\\bin\\%d.bmp',i);
p(:,:,i+1)=imread(s);
%对变量P中的100张图片进行边缘化特征处理,并把处理后的图片放在P2中
p2(:,:,i+1)=edge(p(:,:,i+1));
%显示提取边缘化特征的图片
imshow(p2(:,:,i+1));
end
附录二:
求最大内切圆的半径
%创建一个101*4的矩阵
tx=zeros(101,4);
for d=0:99
%将整形转换成字符串函数,,,并用字符串连接函数将两个字符串连接起来
k=strcat('D:\Program Files\MATLAB\R2011a\bin\',int2str(d),'.bmp');
A=imread(k);%A变量表示每一个图片
for i=1:1:512
for j=1:1:512
A(i,j)=1-A(i,j);%变为灰度图片
end
end
BW=edge(A,'sobel');%自动选择阈值用Sobel算子进行边缘检测
BW2=bwmorph(A,'skel',inf);%提取图像A的骨架,n = Inf时,移除目标边界孤立像素,保留下来的像素组合成图像的骨架
[x,y,z]=find(BW);%找到非0元素
[a,b,c]=find(BW2);
m=length(a);%BW2
n=length(x);%BW
e=zeros(m,n);
f=zeros(m,2);
for i=1:m
for j=1:n
p1=a(i);%BW2
q1=b(i);
p2=x(j);%BW
q2=y(j);
e(i,j)=sqrt((p1-p2)^2+(q1-q2)^2);%内部点到边缘轮廓线上点的距离
end
[zr,zxh]=min(e(i,:)); %内部点到边缘点的距离中取最小值
f(i,1)=zr; %重新赋值给f变量
f(i,2)=zxh;
end
[zR,zxxh]=max(f(:,1));%这些最小值中取最大值,即为最大内切圆的半径
x=a(zxxh)-256; %默认以图片中间为圆心坐标
y=b(zxxh)-256;
tx(d+1,1)=x;%圆心横坐标
tx(d+1,2)=y;%圆心纵坐标
tx(d+1,3)=d+1; %第几幅图片
tx(d+1,4)=zR; %最大内切圆的半径
end
附录三:
求中轴线的拟合曲线及其在三个坐标平面的投影
x=[-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-159;-159;-159;-158;-157;-157;-157;-156;-153;-153;-153;-152;-152;-150;-149;-148;-148;-143;-137;-137;-116;-115;-115;-114;-114;-113;-112;-111;-111;-63;-75;-81;-51;-51;-31;-31;-31;-35;-35;-26;-35;-26;46;46;46;46;65;68;65;81;81;81;135;136;136;137;138;138;139;139;139;140;140;172;172;172;172;182;187;187;187;187;188;];
y=[1;0;2;2;2;2;1;4;1;1;7;8;9;10;12;13;14;16;17;18;19;20;21;22;21;21;21;30;30;29;35;40;40;40;44;55;55;55;58;58;63;66;68;68;78;88;88;115;116;116;117;117;118;119;120;120;151;145;142;156;156;162;162;162;161;161;163;161;163;163;163;163;163;158;157;158;152;152;152;118;117;117;116;115;115;114;114;114;113;113;67;67;67;67;43;24;24;24;24;18;];
z=[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;];
format long
px=polyfit(z,x,5);
x1=polyval(px,z);
py=polyfit(z,y,5);
y1=polyval(py,z);
figure(1);
plot3(x1,y1,z)
grid on
xlabel('X轴');
ylabel('Y轴');
zlabel('Z轴');
title('中轴线图');
figure(2);
plot(z,x1,'-r')
ylabel('Z轴');
xlabel('X轴');
title('中轴线XOZ平面投影图');
grid on
figure(3);
plot(z,y1,'-b')
xlabel('Z轴');
ylabel('Y轴');
title('中轴线YOZ平面投影图');
grid on
figure(4);
plot(x1,y1,'-g')
xlabel('X轴');
ylabel('Y轴');
title('中轴线XOY平面投影图');
grid on
附录四
由切片.bmp图像画出血管三维图像
a=zeros(100,512,512);
for i=0:99
a(i+1,:,:)=imread(strcat(num2str(i),'.bmp'));
end
for i=1:100
for j=106:452
for k=65:476
if a(i,j,k)==0
plot3(j,k,i,'*')
hold on;
end
end
end
end
rotate3d;
hold off
附录五:
血管拟合后还原图
x=[-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-160;-159;-159;-159;-158;-157;-157;-157;-156;-153;-153;-153;-152;-152;-150;-149;-148;-148;-143;-137;-137;-116;-115;-115;-114;-114;-113;-112;-111;-111;-63;-75;-81;-51;-51;-31;-31;-31;-35;-35;-26;-35;-26;46;46;46;46;65;68;65;81;81;81;135;136;136;137;138;138;139;139;139;140;140;172;172;172;172;182;187;187;187;187;188;];
y=[1;0;2;2;2;2;1;4;1;1;7;8;9;10;12;13;14;16;17;18;19;20;21;22;21;21;21;30;30;29;35;40;40;40;44;55;55;55;58;58;63;66;68;68;78;88;88;115;116;116;117;117;118;119;120;120;151;145;142;156;156;162;162;162;161;161;163;161;163;163;163;163;163;158;157;158;152;152;152;118;117;117;116;115;115;114;114;114;113;113;67;67;67;67;43;24;24;24;24;18;];
z=[1;2;3;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;32;33;34;35;36;37;38;39;40;41;42;43;44;45;46;47;48;49;50;51;52;53;54;55;56;57;58;59;60;61;62;63;64;65;66;67;68;69;70;71;72;73;74;75;76;77;78;79;80;81;82;83;84;85;86;87;88;89;90;91;92;93;94;95;96;97;98;99;100;];
cenn=zeros(100,3);
cenn(:,1)=x1; cenn(:,2)=y1; cenn(:,3)=z;
t=linspace(0,pi,25);
p=linspace(0,2*pi,25);
[theta,phi]=meshgrid(t,p);
for i=1:100
x=29.4166*sin(theta).*sin(phi)+cenn(i,1);
y=29.4166*sin(theta).*cos(phi)+cenn(i,2);
z=29.4166*cos(theta)+cenn(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
shading interp