【 数学建模竞赛】血管的三维重建模型g

1、 血管的三维重建模型摘要：本文对血管三维重建中，中轴线及球的半径确定问题进行了讨论。首先，根据问题及图象处理提取有效数据，给出两种可行算法，利用上述数据建立了最大最小方法和二次规划方法。搜索中心点，并给出全局和局部搜索，得到各切片中心点坐标（见表1），并通过插值方式得到中轴线图象及其各投影。最后对模型给出检验方式。一、问题的重述假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线（称为中轴线）的球（命名为包络球）滚动包络而成。现有某管道的相继100张平行切片图象,记录了管道与切片的交。假设：管道中轴线与每张图片有且只有一个交点；球半径固定；切片间距以及图象象素的尺寸均为1.取坐标

2、的Z轴垂直于切片，第1张切片为平面Z二0，第100张切片为平面Z二99.计算管道的中轴线与半径，给出具体的算法，并绘制中轴线在XY、YZ、ZX平面的投影图。二、模型假设与符号说明1、基本假设：（1）该管道的表面为一定长半径的球沿一固定的曲线运动所得曲面族包络的光滑表面。（2）该管道的中轴线连续而且光滑。（3）该管道的中轴线与每个切面有且只有一个交点。（4）图象象素的尺寸为1.（5）切片的间距尺寸为1.2、符号说明：L中轴线R包络球的半径O（x,y,z）中轴线与第i个切片的交点（定为此切片的中心）iS第i个切片切得的图形iD第i个切片的图象数据矩阵i三、问题分析及建模准备问题分析：通常血管的表面

3、可认为是连续且光滑的曲面，断面可用于了解其形态等特性。本问题给出的是一些离散的切面，要求重建出原图中轴线和求出包络球半径。因为每一个切面与中轴线L有且只有一个交点O，如果找出所有O，就可以用插值或拟合的方式作出L的ii近似图象，其在坐标平面上的投影就很容易画出。问题的关健转变为求每个平面上的O.i建模准备：1、图象的读取由于切片图象中只有黑、白两种颜色的象素，而且所给的BMP格式图象文件是512X512象素的.因此，把图象读取为一个512X512的数字矩阵；用数字1表示黑色的象素，用数字0表示白色的象素。如图1把（a）的图象读取为（b）的矩阵：图4 图1111111111(b)2、边界提取51

4、2X512的数字矩阵所占的空间较大，不便处理为节省空间,只用提取图象的边界进行处理,就可达到目的.下面给出两种提取边界的方法：(1)Roberts交叉算子方法如图2，a,a,a为相邻的象素，设黑色象素的值为1，128白色象素的值为0.对于中间象素a8当K=(a-a+(a-a0时，0426可判断其为边界点用此方法对0.bmp判别边界，的出的边界图如图3.aaa012aaa783aaa654可见，边界点重复，多达512个，对求轴心不利.(2)自行设计一种判别法:为使边界细化给出如下判别方法:对3X3膜，当中心点a为黑色象素且白象素三4时，方判其为边界点，如图4.8aoaaaaa783aaa654a

5、oaa2aaa83aaa6541)2) 而对白象素=2或3时,则要用5X5膜类似判定此时主要处理非凹图形的凹边提取.对0.bmp实施此种边界提取方法后,边界点降至82个.如图5四模型建立与求解：模型的理论基础：命题1：由球心沿着某一曲线（称为中轴线）的半径不变的球（即包络球）包络而成的曲面，若被某一平面横截，且截面与中轴线有且仅有一个交点，则此截面中包含的圆中半径最大者r即为包络球半径。证明：如图6,截面与中轴线的交点为O，球O以O为球心，以R为半径。则经过O的所有截面必包含球O的大圆，所以包络球半径Rr.若Rr，则球O与切面边界点必无交点，由中轴O线的连续性，容易知道，曲面必不是由球O包络而

6、成，产生矛盾。所以，R=r.推广：包络球半径变化时，上述命题依然成立。命题2：命题1所述截面，其中任一点到其边界点的最短的即为该点到边界点的距离，截面区域中的点到边界距离最长的即为包络球球半径，该点即为中轴线与该截面的交点。（由命题1，结论显然，证略）。模型一（最大最小法搜索）：a.全局搜索法：由上述结论，我们得到下面的数学模型寻找各切面的中心点，以及得到半径R，其算法如下：1）定义切面各点到边界点的距离。2）求出每一点到边界的距离（定义见命题2）。3）求出所有点距离的最大值，即R，并记录取得最大值的点，即中心点。用Matlab数学软件编程实现算法（程序见附录2），结果见附录3.取Rmax的平

7、均值，我们得到R=28.8974沁29根据附录3的结果，可得到中轴线的立体图及其在XY,YZ,ZX平面的投影图如下图:针对本题，该模型简单易懂，可操作性强，且结果稳定性强，误差较小。但是对数据更大或其他复杂情况效果不能保证。我们提出以下改进：b.局部搜索法：首先我们计算各切片的面积：应用Matlab中的函数bwarea直接求得到各个S（见附录i4）：其中最小者为S，由命题1,SR2,（R为包络球半径）。通过象素求和以圆的形00状，我们用手工在O.bmp图象内逼近得到面积S，则SS我们得到R的大概范围：028R30.9.现在，我们给出两种改进方法：方法1：由我们已经估计出的R值，可以肯定在x轴方

8、向从边界向内缩进27,搜索范围极大缩小。如果要搜索的范围仍然很大，我们还可以从y轴方向继续缩小。方法的选则要具体问题具题分析。方法2：由上面计算的面积图表所列面积的大小反应了切面与该管道所成角Q的大小，角度越小，则面积越大。（这里P定义为过该切片的轴心与曲线相切的直线与切平面的夹角）从而具有最小面积的S与管道所成的角最大（接近90。），故它最能反应出0管道包络球（即包络球O）大圆面，从而球半径R可通过S来估计。0由面积比确定切面与管道夹角。为方便理解，将管到理解为直圆柱，则平行切片变为不平行的切片。S=nab,S=兀a2，cos0=，0=arccos(S/S);bS将SO设为垂直与管道的切片。

9、则第I片与管道夹角0=arccos（S/S）0由附录4,可得0故可得一圆，圆心为0.（S得轴心）半径为r=ixtg0，则Si的轴心0应在圆上。实际应用中可把圆放大为圆带，搜索公共区域即可。模型二：用二次规划的方法求中心：（1）确定包络球的半径R.记切片上最大半径的内切圆的方程为：F（x,y,D,E,F）=x2+y2一2Dx-2Ey+F=0切片边界的点为（x,y）,i=1,2,n.则F（x,y）表示（x,y）与圆的位置关系，若为0iiiiii则点在圆上，大于0则点在圆外。记X=（D,E,FF（X）=x2+y22Dx-2Ey+F.iiiii二次规划的思想是：求X使工F2(X)最小，但点(x,y),

10、i=1,2,n均要在图外。iiii=1规划QP:min工F2(X)Xi=1(i)s.t.F(X)0,i=1,2,n若求得X*=(D*,E*,F为最优解，贝U该切片的轴心即为(D*,E*)，包络球半径R=；D*2+E*2-F*.下将QP变为标准形式(D)EFiiiiYF2(X)=(X；1ii其中A=(一2x,-2y,1、11，b=(、x2+y211、一2x,一2y丄nnx2+y2丿nn约束为0，f=Ab.确定X的约束：由于生成圆半径25R32.-256+30D256-30-256+30E256-30-Sb(X25R=xD2+E2-F32D2+E2一322FD2+E2一252取F的约束为1000F

11、80000.将常数项bb从目标中去掉，目标除以2,并令H=AA，而得标准的二次规划问题:min2s.t.-AXbLBXUB该规划易用Matlab求解。(程序见附录5)0.bmp的结果为：*=(160,1,24768)，而轴心坐标为(-160,1)，半径R=28.86下图给出了用QP求得的该切面图及所求生成图的示意图-220、(220、其中LB=-220,UB=220-1000/0 x0n规划为：min乂申（F（Y）化简得：miniF(Y)0i(F(Y)iii-AZ-ZZ-b+R20s.tLBZ8060、40、20、0200图13管道表面立体图五模型检验由于管道可视为球，沿中轴线运动所得的一曲面

12、族的的包络。故我们可生成某一已知中轴线的管道。然后，求出其切片，再用本文方法包络球的半径，将切片重构出其曲面的中轴线，以校验本文的方法的正确性。附介绍包络的概念：一组情形：设平面上一族曲线为fCx,y,a)二0,a为参数，包络为l。则由过l上的点P,y)必有某a使曲线fCx,y,a0)=0与之相切。可推出包络线l的方程为f(x,y,a)=0，f(x,y,a)二0空间曲面包络边有类似结果。设x=xQy=yC)z=zC)为给定的空间曲线，球心在该曲线的一族半径为r的球面方程为：(x一xQI+(y一yQI+(z一zQI=r2则其包络方程为(x一xC)2+(y一y(y)1+(z一zC)2=r2(x-xC)xr(t)+(y-y(y)yr(t)+(z-zC)z，C)=0.六模型评价本文建模思想易于理解。模型一可操作性较强，应用方便。模型二可操作性更强，所得结论可以广泛推广。本文所用的