血管的三维重建摘要本文探讨血管的三维重建,由血管的相继100张平行切片图像计算血管的中轴线与半径,并绘制血管在三个坐标平面上的投影。由于血管的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成,由此我们得出结论:每个切片一定包含滚动球的大圆,并且它一定为切片的最大内切圆,而最大圆所对应的半径即为血管半径,所以求血管半径就转化为求每一个切片内部的点到切片外部轮廓的最大半径。首先,读取100张血管切面图,把它们转换成Logical矩阵,从中提取切片截面轮廓点构成一个新的矩阵。然后找到原图片矩阵中像素点的内点(切片图片中轮廓线中的点),从而得到内点到切片轮廓点的最小距离矩阵和最小距离中的最大值矩阵,最大值即为血管半径。最后计算所有切片的血管半径,并对这些半径求平均值,得到平均血管半径为:29.6799μm。由100张切片的最大内切圆圆心坐标拟合得出中轴线方程以及其在三个坐标平面内的投影曲线方程。由中轴线得到血管的三维立体重建图,用平面)74,49,24,0(,iiZ去截血管的三维立体重建图,得到新的4张截面图。把它们分别与题设中的对应截面进行内点个数对比。我们定义两张切片所共同拥有的内点个数与原切片内点个数的比值为重合度。计算得到平均重合度为:98.19%。关键词:血管半径中轴线切片重建(来自作者:欢迎各界人士批评指正,学术交流邮箱nibz@qq.com。文章作于2011年8月10日,陕西科技大学理学院实验室)1问题的重述断面可用于了解生物组织、器官等的形态。例如,将样本染色后切成厚约1μm的切片,在显微镜下观察该横断面的组织形态结构。如果用切片机连续不断地将样本切成数十、成百的平行切片,可依次逐片观察。根据拍照并采样得到的平行切片数字图象,运用计算机可重建组织、器官等准确的三维形态。假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。例如圆柱就是这样一种管道,其中轴线为直线,由半径固定的球滚动包络形成。现有某管道的相继100张平行切片图像,记录了管道与切片的交。图像文件名依次为0.bmp、1.bmp、…、99.bmp,格式均为BMP,宽、高均为512个象素(pixel)。为简化起见,假设:管道中轴线与每张切片有且只有一个交点;球半径固定;切片间距以及图像象素的尺寸均为1。计算管道的中轴线与半径,给出具体的算法,并绘制中轴线在XY、YZ、ZX平面的投影图。2模型假设(1)假设血管中轴线与每张切片有且只有一个交点。(2)假设血管的半径固定不变。(3)假设在对切片拍照的过程中不存在误差。(4)假设血管不存在严重扭曲。3符号说明kC第k张切片中轮廓线的最大内切圆圆心kR第k张切片中轮廓线的最大内切圆半径ijkd第k张切片中第i个内点与第j个轮廓点的距离100512512PA100张切片图片转换以后的三维0-1矩阵100512512PB100张切片图片的轮廓线生成的矩阵3100cen100张切片轮廓的最大内切圆圆心坐标R平均血管半径k第k张切片图片中所有内点的集合k第k张切片图片中轮廓点的集合S原切片图片的上内点及边界点的集合T重新切片得到的内点及边界点的集合R-square多项式拟合的指标数4问题的分析根据题目整个管道是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。基于几何原理可以得出以下两条定理。(1)球的任意截面都是圆。(2)经过球心的球截面是所有的截面圆当中半径最大的圆。基于上述两个定理可以得出:每张切片的最大内切圆的圆心位于血管的中轴线上,该圆的半径等于血管半径。分别求出100张切片轮廓线的最大内切圆圆心kC与半径kR。对所有最大内切圆的圆心kC做拟合,即可求出血管的中轴线。为减少误差,我们用求均值的方法1001001)(kkRR来得到平均血管半径R。根据拟合出的中轴线方程,即可求出中轴线在XY、YZ、ZX平面上的投影。5模型的建立与求解5.1模型建立5.1.1求血管的半径(1)导入数据,转换存储方式为了计算方便,利用MATLAB[1]软件将100张BMP文件转换存储为一个三维0-1矩阵100512512PA(0代表黑色像素点,1代表白色像素点)。(2)求切片轮廓线上各点坐标利用MATLAB软件的内部函数edge,求得所有切片图片的轮廓线生成的矩阵100512512PB(3)求切片轮廓线的最大内切圆的半径首先求出每个内点距轮廓线的距离,取其中的最小值,即为以这些内点为圆心的轮廓线的内切圆半径;其次找出所有内点确定的最小内切圆中半径最大的一个,即为轮廓线的最大内切圆半径。具体步骤:Step1在三维0-1矩阵100512512PA中遍历搜索内点,将其存储于集合)(k中;Step2在轮廓线矩阵100512512PB中遍历搜索轮廓点,将其存储于集合)(k中;Step3计算集合)(k中第i个内点yxA,到集合)(k中第j个轮廓点qpB,的距离22qypxdijk,首先取每个内点到轮廓点的最小距离,即为以这些内点为圆心的内切圆半径;其次找出所有内点确定的最小内切圆中半径最大的一个,即为轮廓线的最大内切圆半径}{minmaxijkjikdRStep4计算平均血管半径1001001)(kkRR(4)求切片轮廓线的最大内切圆心坐标将100张轮廓线的最大内切圆半径所对应的圆心坐标),(wzCk记录于矩阵)3,100(cen中。)100,...,2,1(,,,)3,()2,()1,(kkcenwcenzcenkkk5.1.2由圆心坐标求中轴线方程及曲线投影(1)利用MATLAB画图工具由100组圆心坐标画出中轴线以及中轴线在XY、YZ、ZX面的投影,由于图像文件以象素为单位存储数据,这使得图像的原始数据是离散化的。对于离散的位图来说,将变得不可分辨,可能切片轮廓线的最大圆不止一个。为了减少误差,将中轴线在XY、YZ、ZX面的曲线方程进行多项式拟合,建立拟合方程:)(1XFY)(2YFZ)(3ZFX(2)由建立的投影曲面得到中轴线的方程);();();(213TFZTFYTFX5.2模型求解根据所建模型求出100张切片的最大内切圆的圆心坐标、半径以及血管半径表1最大内切圆半径及圆心坐标切片序号最大内切圆半径(R/μm)最大内切圆圆心坐标横坐标(X/μm)纵坐标(Y/μm)竖坐标(Z/μm)029.069-16010129.069-16011229-16022329.069-16023429.069-16024529.069-16025629.155-16026729.155-16027829.275-16028929.275-160391029.428-1603101129.428-1604111229.614-1604121329.614-1605131429.833-1606141529.833-1607151629.833-1608161730-1609171830-16010181930-16011192030-16012202130-16013212230-16014222330-16015232430-16015242529.833-16016252629.614-15927262729.614-15927272829.614-15928282929.614-15834293029.614-15834303129.614-15740313229.614-15548323329.614-15548333429.614-15548343529.732-15355353629.732-15258363729.732-15258373829.732-15258383929.614-15258394029.547-14868404129.547-14966414229.53-14084424329.53-14084434429.682-13592444529.682-13592454629.682-13592464729.698-115116474829.698-113118484929.698-112119495029.698-111120505129.698-111120515229.698-111120525329.698-111120535429.682-111120545529.53-66150555629.53-59153565729.732-70148575829.967-70148585929.732-54155596029.732-54155606129.614-31162616229.614-31162626329.614-6166636429.833-6166646529.83381676566301016766673011167676830121676869301216769703012167707129.83326166717229.61426166727329.61446163737429.73265158747529.73265158757629.73271156767729.68296145777829.68296145787929.698127124798029.698127124808129.732122128818229.732122128828329.732122128838429.698127124848529.732150101858629.732150101868729.73215496878829.698137115888929.698139113899029.68216088909129.68216285919229.54717659929329.61418340939429.73218049949529.61418340959629.61418340969729.614185339798301900989930190199(1)血管半径为R=29.6799μm。(2)中轴线在XY平面投影曲线的拟合方程(R-square:0.9939)25566-1041.7X18.059X-0.17693X0.00108X-X104.26X101.09-X101.74X101.58-X106.232345-66-87-118-149-18Y(3)中轴线在YZ平面投影曲线的拟合方程(R-square:1)86253456-57-88-119-15103.60Y109.75-Y101.17813.88Y-3.6293Y0.01075Y-Y102.11Y102.66-Y101.95Y10-6.32Z(4)中轴线在ZX平面投影曲线的拟合方程(R-square:1)92.9882.4609Z0.55402Z-0.053686Z0.00269Z-Z107.56Z101.23-Z101.16Z105.96-Z101.292345-56-67-88-119-13X(5)中轴线方程:,92.9882.4609Z0.55402Z-0.053686Z0.00269Z-Z107.56Z101.23-Z101.16Z105.96-Z101.29,103.60Y109.75-Y101.17813.88Y-3.6293Y0.01075Y-Y102.11Y102.66-Y101.95Y10-6.3225566,-1041.7X18.059X-0.17693X0.00108X-X104.26X101.09-X101.74X101.58-X106.2323455-66-78-811-913-862534565-78-811-915-2345-66-87-118-149-18XZY0100200300400500250300350400450020406080100x中轴线曲线yz010020030040050