血管三维重建01年A题 数学建模.doc

A题 血管的三维重建

摘要

对于血管的三维重建问题，关键是找出理论假设下血管的中轴线和血管的直径。通过这两个参数的确定就可以基本上绘制出血管图来。
首先，对给出的100张血管断面的二进制图像进行取反操作，应用软件获得图片中血管图像内部点和边缘点的坐标，根据程序求得最大内切圆的半径及其圆心坐标。具体是（1）求出内部任意一点与边缘点的距离，取距离最小的值作为以该点为圆心的内切圆半径；（2）在这些半径中找到数值最大的值即为该图像的最大内切圆半径，对应的点坐标即为最大内切圆的圆心坐标。对所有求得的半径取平均值得到。
然后，根据100个圆心坐标进行多项式曲线拟合，得到中轴线方程为：
剩余标准差为1.5518，说明此回归模型的显著性好。
绘制出曲线图，并投影到X-Y，X-Z，Y-Z坐标面上。

关键词：边界提取 最大内切圆法 多项式拟合

1 问题重述

断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约的切片，在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片， 可依次逐片观察。根据拍照并采样得到的平行切片数字图象，运用计算机可重建组织、器官等准确的三维形态。
现有某管道的相继100张平行切片图像，记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、…、 99.bmp，格式均为BMP，宽、高均为512个象素（pixel）。
试计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在X-Z、Y-Z、X-Y平面的投影图。

2 模型假设

（1）血管可视为一类特殊的管道，该管道的表面是由球心沿着某一曲线（称为中轴线）的球滚动包络而成；
（2）球半径固定；
（3）管道中轴线与每张切片有且只有一个交点；
（4）图片间距以及图像像素尺寸均为1。

3 符号说明

：切片数，从0到99。
：切片图像最大内切圆圆心范围坐标。
：切片图像边缘点的坐标。
：切片图像的最大内切圆半径。
：切片图像最大内切圆圆心的横坐标。
：切片图像最大内切圆圆心的纵坐标。
：从（-9，-9）点到（10,10）点矩形范围内所有横、纵坐标为整数的点。
：切片图像中任意一点到的最短距离

4 问题分析

本题给出了某管道的相继100张平行切片图像，记录了管道与切片的交。由于该血管是由球心沿着某一曲线的球滚动包络而成，且球的半径固定。图像格式均为BMP，宽、高均为512个象素（pixel），图像象素的尺寸为1，对图像进行取反操作，一是函数求解需要，二是得到最大内切圆半径为图像内部球体半径更符合要求。所以我们得到以下两个结论：（1）每张图片图形的最大内切圆的半径即为滚动球体的半径；（2）圆心坐标即为滚动球体的球心坐标，由于切片间距也为1，故可以得到滚动球体的三维坐标。
根据每张图片最大内切圆的圆心拟合出管道的中轴线，并绘制出来，将绘制出的中轴线分别投影到X-Y、Y-Z、Z-X平面即可得到中轴线在X-Y、Y-Z、Z-X平面的投影图。

5模型建立与求解

首先用画图板对100张图片取反色，然后导入到matlab中。
（1） 求得切片图像边缘点坐标及第一张切片图像的内部点坐标；
（2） 求中任意一点到的最短距离；
（3） 在中取最大值即为第一张切片图的最大内切圆半径，对应的内的一点的坐标即为最大内切圆的圆心坐标（，）；
（4） 求下一张图像的最大内切圆圆心范围；
（5） 求中任意一点到的最短距离；
（6） 在中取最大值即为第Z张切片图的最大内切圆半径，对应的内点的坐标即为最大内切圆的圆心坐标（，）；
（7） 检验所求圆心是否在原来预计的圆心范围内。
若不符合与
则 令 返回到（5）
（8） 利用上面得到的中心点坐标通过Matlab拟合出中轴线方程，并绘制出。
求得的最大内切圆半径和圆心坐标见附录表一。
通过对100个最大内切圆半径求平均值，得到。
利用matlab多项式拟合函数polyfit程序五进行三维曲线的拟合，得到中轴线方程为：

剩余标准差为1.5518，数值小，说明此回归模型的显著性好。
中轴线三维图为：

图一 中轴线三维图
中轴线X-Y面上的投影：

图二 中轴线X-Y面上的投影
中轴线X-Z面上的投影：

图三 中轴线X-Z面上的投影

中轴线Y-Z面上的投影：

图四 中轴线Y-Z面上的投影

6 模型检验

画出由图片信息绘制三维点图程序，根据求得中轴线和半径绘制三维血管图。
实际血管切片图像边缘包络出的血管图在各个坐标面上的投影：

图五 实际血管X-Y平面投影

图六 实际血管X-Z平面投影

图七 实际血管Y-Z平面投影

通过球体沿着中轴线滚出的血管图在各个坐标面上的投影：

图八 理论血管X-Y平面投影

图九 理论血管X-Z平面投影

图十 理论血管Y-Z平面投影

将模拟出来的实际血管投影图和理论血管投影图进行重叠，重合的部分不显示，进行投影比较：
实际血管和理论血管X-Y平面投影图比较：

图十一 X-Y平面投影图比较
实际血管和理论血管X-Z平面投影图比较：

图十二 X-Z平面投影图比较
实际血管和理论血管Y-Z平面投影图比较：

图十三 Y-Z平面投影图比较
选取Y-Z平面投影，用matlab求出图形面积。

图十四 实际Y-Z平面投影 图十五 X-Z平面投影图比较（经处理）
算的实际血管投影面积为：
实际与理论面积偏差：
偏差比
误差较小，说明计算得到的理论血管比较符合实际血管。

7 模型评价

此模型优点在于对原始图片进行取反操作能较为准确的得到血管的三维图，多使用matlab函数进行求解，方便快捷。
缺点在于模型检验是通过整体面积比较。没有做出算出的血管的切片和原始图片进行比较，准确性不高。

8 模型推广

此模型可以由实际的图片切片得到三维图的关键数据信息已达到重建三维模型的目的。可广泛应用于有数学特征的模型的重建上。比如，可以去点等半径球体假设，球体半径按照规律可变，或者血管的法向切面为其他形体的实体的重建。

参考文献

[1] 姜启源，谢金星，叶俊，《数学模型》（第三版），北京：高等教育出版社，2003.8
[2] 阳明盛、熊西文、林建华《MATLAB基础及数学软件》，大连理工大学出版社，2003年8月第一版
参考网站：
[1] 中国数学建模网， 

附录：

表一 最大内切圆半径和圆心坐标
序号（Z坐标）
半径
X坐标
Y坐标
0
29
257
416
1
29
257
416
2
29
257
416
3
29
257
416
4
29
257
416
5
29
257
416
6
28.8617
257
416
7
29
258
416
8
29
258
416
9
29
258
416
10
29
258
416
11
29
259
416
12
29
259
416
13
29
259
416
14
29
260
416
15
29
261
416
16
29
261
416
17
29
262
416
18
29
263
416
19
29
264
416
20
29
265
416
21
29
266
416
22
29
267
416
23
29
268
416
24
29.0172
277
416
25
29.0172
276
416
26
29.0172
286
415
27
29.0689
286
415
28
29.0689
292
414
29
29.0689
292
414
30
29.1204
297
413
31
29.1548
292
414
32
29.1204
305
411
33
29.1204
305
411
34
29.1204
308
410
35
29.1204
314
408
36
29.1204
317
407
37
29.1204
317
407
38
29.1548
322
405
39
29.1548
322
405
40
29.1548
322
405
41
29.1548
345
393
42
29.1548
348
391
43
29.1548
351
389
44
29.4109
348
391
45
29.4109
348
391
46
29.4109
348
391
47
29.6985
371
372
48
29.6985
373
370
49
29.6985
374
369
50
29.6985
375
368
51
29.6985
375
368
52
29.6985
375
368
53
29.6985
375
368
54
29.4109
376
367
55
29.2062
389
351
56
29.2062
389
351
57
29.4109
403
328
58
29.4109
404
326
59
29.1548
409
315
60
29.1548
409
315
61
29.1204
418
286
62
29.1204
419
281
63
29.1548
420
275
64
29.1548
420
275
65
29.1548
420
275
66
29.1204
419
281
67
29.0689
421
268
68
29.0172
421
268
69
29.0172
422
236
70
29.1204
419
211
71
29.1204
416
198
72
29.1204
413
188
73
29.1204
410
180
74
29.1204
410
180
75
29.1548
408
175
76
29.1548
405
168
77
29.1548
388
139
78
29.1548
388
139
79
29.2062
382
131
80
29.2062
376
124
81
29.4109
376
124
82
29.6816
376
124
83
29.6816
375
123
84
29.6816
363
111
85
29.6816
362
110
86
29.6816
362
110
87
29.4109
362
110
88
29.2062
347
98
89
29.2062
347
98
90
29.1548
333
89
91
29.1548
310
78
92
29.1548
310
78
93
29.1548
310
78
94
29.1548
310
78
95
29.1204
283
70
96
29.1204
278
69
97
29.1204
278
69
98
29.1204
278
69
99
29.1204
288
71

Matlab:

程序一：
求中心点范围程序：
function a=g(xx,yy)
for j=1:20
for i=1:20
a(j*20+i-20,1)=xx-10+i;
a(j*20+i-20,2)=yy-10+j;
end
End

程序二：
求最大内切圆半径和圆心坐标程序：
BW10=imread('0.bmp');
BW20=bwmorph(BW10,'remove');
[i0,j0]=find(BW10==1);
b(:,1)=j0;
b(:,2)=512-i0;
[r xx yy]=ff(a,b)

程序三：
计算最大半径程序：
function [r xx yy]=ff(a,b)
m=size(a,1);
n=size(b,1);
w=zeros(m,1)+40;
h=0;
for i=1:m
for j=1:n
D(i,j)=sqrt((a(i,1)-b(j,1))^2+(a(i,2)-b(j,2))^2);
if D(i,j)<w(i,1)
w(i,1)=D(i,j);
x(i,1)=a(i,1);
y(i,2)=a(i,2);
end
end
end
[r index]=max(w);
xx=a(index,1);
yy=a(index,2);

程序四：
求每张图片中个点的坐标程序：
BW10=imread('0.bmp');
BW20=bwmorph(BW10,'remove');
[p0,q0]=find(BW20==1);
[i0,j0]=find(BW10==1)

程序五：
曲线拟合程序：
p=polyfit(z,x,6);q=polyfit(z,y,6);
xx=polyval(p,z);yy=polyval(q,z)
plot3(x,y,z,'.r',xx,yy,z)

程序六：
连续读图程序：
I=zeros(99,0);
for i=0:99
url=strcat(num2str(i));
imageurl=strcat(url,'.bmp');
I=imread(imageurl);
I=imresize(I,[512 512],'bilinear');
figure;imshow(I);
Imgstore(:,:,:,i)=I;
end

程序七：
轮廓提取程序：
BW1=imread('0.bmp');
imshow(BW1);
BW2=bwmorph(BW1,'remove');
BW3=bwmorph(BW1,'skel',inf);
imshow(BW2)
figure,imshow(BW3)
[p,q]=find(BW2==1);
plot(q,-p,'.r')

程序八：
根据图片信息绘制三维点图程序：
clear;
for b=0:99
m1=imread([int2str(b),'.bmp']);
m(:,:,b+1)=edge(m1,'sobel');
end
for k=0:99
for i=1:512
for j=1:512
if(m(i,j,k+1)==1)
plot3(j,512-i,k+1,'b-.');hold on
end
end
end
end
程序九：
根据求得中轴线和半径绘制三维血管图程序：
[X0,Y0,Z0]=sphere(30);
surf(29*X0+257,29*Y0+416,29*Z0+0);
axis equal;
hold on;
for i=1:99;
[Xi,Yi,Zi]=sphere(30);
surf(R(i)*Xi+O(i),R(i)*Yi+P(i),R(i)*Zi+Q(i));
axis equal;
end
14