血管的三维重建

   血管的三维重建问题：

       断面可用于了解生物组织、器官等的形态。例如，将样本染色后切成厚约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个点的坐标为 (x_i,y_i)，切片轮廓线上有q个点，第j个点的坐标为(m_j,n_j),则切片内部第i点到轮廓线上第j点的距离为：

  

第i个点对应的最小距离r_i 为：

                                      

 

（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年

附录

 

 

附录一：
画出血管的轮廓点

p=zeros(512,512,100);p2=zeros(512,512,100);for i=0:99s=sprintf('D:\\program Files\\MATLAB\\R2011a\\bin\\%d.bmp',i);p(:,:,i+1)=imread(s);p2(:,:,i+1)=edge(p(:,:,i+1));imshow(p2(:,:,i+1));end

附录二：
求最大内切圆的半径

tx=zeros(101,4);for d=0:99     k=strcat('D:\Program Files\MATLAB\R2011a\bin\',int2str(d),'.bmp');    A=imread(k);    for i=1:1:512        for j=1:1:512            A(i,j)=1-A(i,j);        end    endBW=edge(A,'sobel'); BW2=bwmorph(A,'skel',inf);[x,y,z]=find(BW); [a,b,c]=find(BW2);m=length(a);n=length(x);e=zeros(m,n); f=zeros(m,2); for i=1:m    for j=1:n        p1=a(i);        q1=b(i);         p2=x(j);        q2=y(j);         e(i,j)=sqrt((p1-p2)^2+(q1-q2)^2);    end   [zr,zxh]=min(e(i,:));    f(i,1)=zr;    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 longpx=polyfit(z,x,5); x1=polyval(px,z);py=polyfit(z,y,5); y1=polyval(py,z);figure(1);  plot3(x1,y1,z)grid onxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('中轴线图');figure(2);plot(z,x1,'-r')ylabel('Z轴');xlabel('X轴');title('中轴线XOZ平面投影图');grid onfigure(3);plot(z,y1,'-b')xlabel('Z轴');ylabel('Y轴');title('中轴线YOZ平面投影图');grid onfigure(4);plot(x1,y1,'-g')xlabel('X轴');ylabel('Y轴');title('中轴线XOY平面投影图');grid on

附录四

由切片.bmp图像画出血管三维图像

a=zeros(100,512,512);for i=0:99a(i+1,:,:)=imread(strcat(num2str(i),'.bmp'));endfor i=1:100for j=106:452for k=65:476if a(i,j,k)==0plot3(j,k,i,'*')hold on;endendendendrotate3d;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;endxlabel('X轴');ylabel('Y轴');zlabel('Z轴');title('血管拟合后还原图');hold offshading interp