血管的三维重建

血管的三维重建问题 摘要 越来越多的在医学领域及生物领域需要大量的切片工作 来研 究切片表面的组织形态 因而如何能将这些切片在利用计算机进行 复原 显得很重要 在对管道的半径的求解中 我们采用的是枚举法 即逐个将 100 切片的最大内切圆的半径都求出来 最后取其平均值 即为我 们所求的管道的半径 在半径求出的同时把最大内切圆的圆心也求 出来 要求得切片的最大内切圆 首先对切片图像 0 1 矩阵转化 然后利用 edge 函数找到图形的轮廓 利用 bwmorph 函数求得图形 的轮廓 然后利用 find 函数求出轮廓与骨架上的坐标 利用两点间 的距离公式 我们就能求出骨架上每个点到轮廓的距离 在所有点 到轮廓的最短距离中的那个最大值 就是我们要求的最大内切圆半 径 此圆心即为中轴线上的点 在中轴线与中轴线在 XY YZ ZX 平面的投影的计算中 我们主要 采用多项式拟合 polyfit 的方法 从得出拟合曲线可以看出与实际图 形很接近 最后我们拟合出了血管的三维形态 关键词 最大内切圆 半径 中轴线 多项式拟合 一 问题重述 如今在医学领域 生物学领域等都需要了解生物组织 器官 等的断面形态 以便于能够准确的研究该样品 得出具有说服力的 研究成果 例如 将样本染色后切成厚约 1 m 的切片 在显微镜下 观察该横断面的组织形态结构 如果用切片机连续不断地将样本切 成数十 成百的平行切片 可依次逐片观察 这时候我们如果我们 需要知道血管的三维形态 根据拍照并采样得到的平行切片数字图 象 运用计算机就可重建组织 器官等的三维形态 假设某些血管 可视为一类特殊的管道 该管道的表面是由球心沿着某一曲线 称 为中轴线 的球滚动包络而成 本文给出了某管道的相继 100 张平行切片图象 记录了管道与 切片的交 图象文件名依次为 0 bmp 1 bmp 99 bmp 格式 均为 BMP 宽 高均为 512 个象素 pixel 我们要做的就是根据 这些切片计算管道的中轴线与半径 给出具体的算法 并绘制中轴 线在 XY YZ ZX 平面的投影图 二 问题分析 本题是一个求最优解的问题 通过具体的算法得出最接近实际 血管的中轴线与半径 并能较准确的绘制出中轴线在 XY YZ ZX 平 面的投影图 得出直观的图形 需要处理大量的图片 得出大量的 数据信息 所以为处理该问题的程序中需要大量的循环 所以程序 的完成与运行皆很复杂 本题管道为由球心沿着某一曲线的球滚动包络而成 且管道中 轴线与每张切片被假设有且只有一个交点 则每张切片一定包含球 的最大截面 即过球心的圆面 则每张切片的最大内切圆就是过球 心的圆面 由于本题中取坐标系的 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 则每个切片的最大内切圆的圆心的 z 坐标已给 再加上 x y 则中心 轴就是每个切片的最大内切圆的圆心 x y z 连成的光滑曲线 而管 道的半径就是这个球的半径 即最大内切圆的半径 进而再求出中 轴线在 XY YZ ZX 平面的投影图 那么这道题的问题最后就转化为 求每张切片的最大内切圆的圆心与半径问题 三 模型假设 1 管道的表面是由球心沿着某一曲线 称为中轴线 的球滚动包络 而成 2 管道中轴线与每张切片有且只有一个交点 3 球半径固定 4 切片间距以及图象象素的尺寸均为 1 四 符号说明 一张切片内骨架点 i 到边缘点 的距离 第 i 张图的最大内切圆半径 半径的平均值 五 模型的建立与求解 问题一 求管道的半径与每张切片的最大内切圆的圆心 要求得管道的半径 首先需要求出每个切片的最大内切圆 为了得出比较精确的结果 我们采用枚举法 对与 100 张切片 考 虑到求得的最大内切圆半径由于误差的存在肯定会有所不同 我们 将所有切片的最大内切圆半径求出后 然后取其平均值 为得到每张图片的最大内切圆半径 首先先将二值图片转换 成 0 1 矩阵 矩阵横纵坐标对应原图像的直角坐标系位置 1 代表 黑像素 0 代表白像素 为了后面 find 函数寻找矩阵中为 1 的坐 标并记录 通过相关函数找出轮廓 edge 和骨架 bwmorph 记录 轮廓 q i j 和骨架 p i j 的坐标 利用两点间的公式去求骨架上的点到边界上所有点的最小距离即骨 架上此点的内切圆的半径 1 骨架上所有点的最大内切圆即为该切片的最大内切圆 2 这里 i 代表第 i 张切片 程序见附录 1 其半径即为最大内切圆半径 通过循环找出最大半径并记录该 半径所对的坐标点 基于此我们能求出所有切片的最大内切圆的半 径与圆心 数据见附录 5 通过求得的 100 张切片的最大内切圆的半径我们取其平均值得 出更加精确的管道的半径 问题二 求中心轴与中心轴在 XY YZ ZX 平面的投影 x y z 为切片的最大内切圆的圆心坐标 1 利用问题一中求得数据 x y z 然后根据多项式拟合 polyfit 得出 x y z 的拟合方程为 中轴线在 Z Y 平面的投影拟合多项式 中轴线在 X Y 平面的投影拟合多项式 中轴线在 Z X 平面的投影拟合多项式 程序见附录 2 得到的拟合曲线如下 图 1 为中轴线在 XY 平面的投影 曲线代表多项式拟合的结果 折线代表真实数据的图形 图 2 为中轴线在 ZY 平面的投影 曲线代表多项式拟合的结果 折线代表真实数据的图形 图 为中轴线在 ZX 平面的投影图 曲线代表多项式拟合的结果 折线代表真实数据的图 形 从以上拟合图形中我们可以看出 通过多项式拟合的投影图与 通过坐标点得出的投影图基本相同 所以可以近似认为拟合的投影 图即为我们所要求得的中轴线在 XY YZ ZX 平面的投影 2 同样利用多项式拟合 polyfit 得出管道的中轴线在三维坐标中的 拟合曲线如下 程序见附录 3 图 为中轴线在三维坐标的立体图 3 利用球坐标以及以上求出的最大半径 以及拟合时在 z 时的坐标 x1 y1 利用 plot3 程序如附录四 图 血管的三维重建 六 模型评价及总结 解决本题的首要前提是要理解像素的概念 在解此题的最初由 于对图片的像素没有明确的理解而浪费了很长时间 在我们的模型 中还有尚有不完美的地方 在求切片的最大内切圆半径与圆心的时 候 没有编写出能够一下就将所有切片的数据都求出的程序 而是 仅仅编写出一个图片的程序 然后逐一求出 在这个问题上花费的 时间长 且过程繁琐 但是在进行曲线拟合时 由于运用 知道了相关函数的用法 运用 polyfit 进行拟合 很大程度上节省 了时间 拟合出的相对来说比较满意 立体切片的三维重建在医学 与生物学等领域的作用很大 因为这些领域往往会研究一些物体切 片的组织结构 所以此问题的实现能够解决很多现实问题 为一些 研究领域提供很大的帮助 七 参考文献 1 董辰辉 MATLAB2008 全程指南 电子工业出版社 2 张红云等 图像骨架的提取的应用 2010 年第四期 3 4 附录 1 jieguo zeros 100 4 J0 imread E 99 BMP for i 1 512 for j 1 512 j0 i j 1 J0 i j 转化为黑色为 1 白色为 0 为了后面 find 函数寻找矩阵中 为 1 的坐标并记录 end end lk edge j0 sobel gj bwmorph j0 skel inf x0 y0 v0 find lk 找到轮廓灰色区域 a0 b0 c0 find gj 找到骨架灰色区域 m length a0 骨架灰色区域个数 n length x0 轮廓灰色区域个数 jl zeros m n 建立 0 矩阵为求内切圆半径 cf zeros m 2 建立两列 0 矩阵为存放中心点坐标 for i 1 m for j 1 n p1 a0 i q1 b0 i p2 x0 j q2 y0 j 骨架 轮廓坐标赋值 jl i j sqrt p1 p2 2 q1 q2 2 end zx zxxh min jl i 骨架上一点到轮廓的最短距离 即骨架上各个点的内切圆的半径 cf i 1 zx cf i 2 zxxh end zd zdxh max cf 1 找到其中最大的半径 并把半径和圆心坐标存储 x a0 zdxh 256 与题目所给坐标轴对应 y b0 zdxh 256 与题目所给坐标轴对应 jieguo k 1 1 x x 轴坐标 jieguo k 1 2 y jieguo k 1 3 k 1 jieguo k 1 4 zd 2 polyfit Z X 7 plot Z X 0 10 99 polyval ans 0 10 99 p s polyfit Z X 7 polyfit Z Y 7 plot Z Y 0 10 99 polyval ans 0 10 99 p s polyfit Z Y 7 plot X Y 160 10 188 polyval ans 160 10 188 3 format long px polyfit z x 7 x1 polyval px z 取拟合后再 z 点时的点 py polyfit z y 7 y1 polyval py z figure 1 plot3 x1 y1 z r t linspace 0 pi 25 存储平均分 pi 为 25 份存储 p linspace 0 2 pi 25 存储平均分 2pi 为 25 份存储 theta phi meshgrid t p 存储角度 for i 1 100 ceen 1 x1 ceen 2 y1 ceen 3 z x 29 49288 sin theta sin phi ceen i 1 y 29 49288 sin theta cos phi ceen i 2 z1 29 49288 cos theta ceen i 3 hold on 使图像不被覆盖 plot3 x y z1 axis equal hold off end 5 切片号半径x 坐标y 坐标z 坐标切片号半径x 坐标y 坐标z 坐标 129 069 160005130 414 11411750 228 284 160215230 414 11411751 329 160125330 414 11311852 429 069 160135430 414 11211953 529 069 160245530 11112054 629 069 160255629 732 11112055 729 160165729 698 11112056 829 017 160475829 698 11112057 929 160185929 53 8114258 1028 862 160196029 547 5115659 1128 862 1607106129 547 5115660 1228 862 1608116229 614 3116261 1328 862 1609126329 614 3116262 1429 017 16010136429 614 3116263 1529 017 16012146529 614 3516164 1629 017 16013156629 614 3516165 1729 017 16014166729 428 2616366 1829 017 16016176829 411 3516167 1929 017 16017186929 275 2616368 2029 017 16018197029 4284616369 2129 017 16019207129 6144616370 2229 017 16020217229 6144616371 2329 017 16021227329 6144616372 2429 017 16022237429 6146515873 2529 069 16021247529 7326815774 2629 069 16021257629 7326515875 2729 069 16021267729 5478115276 2829 155538115277 2929 275538115278 3029 27569813511779 3129 42869813611680 3229 61469813611681 3329 61469813711582 3429 61469813811483 3529 61469813811484 3629 73269813911385 3729 73269813911386 3829 73269813911387 3929 73269814011288 4029 61469814011289 4129 5475317267