专利名称：基于重力场模型的管道状物体拉直的图像处理方法
技术领域：
本发明涉及的是一种图像处理技术领域的方法，具体是一种基于重力场模型 的管道状物体拉直的图像处理方法。
背景技术：
虚拟可视化技术是一种观察管道状物体的三维图像可视化技术。这个技术的 优势在于无需破环或接触物体本身，通过计算机图像处理模拟切开、展平、外翻 等方式，以观察管道内壁结构和状态。虚拟拉直可以将弯的管道映射到直的柱面 上，有助于虚拟展平、外翻可视化技术的实施，既可以保持转弯处的细节信息， 又能够提供直观的整体信息。
经过对现有技术的文献检索发现，美国Iowa大学放射系王革(Ge Wang)教授 等于1995年首次提出了基于电场模型和中心路径拉直的曲截面结肠拉直展平方 法，阐述了结肠拉直在虚拟可视化技术中的重要作用(G. WangandM. W. Vannier， 〃GI tract unraveling by spiral CT, 〃 in Proc. SPIE， San Diego, CA, USA, 卯.307-315, 1995)。这个方法的不足在于对于管道状物体弯曲较大(中心路 径曲率较大)的区域，由于同种电荷的互斥作用，基于电场模型的曲截面呈现不 对称的分布，这将导致拉直变形后在该区域产生严重的失真；对于管道内部由于 凸起、凹陷等原因造成的不对称结构区域，基于中心路径拉直的方法会使管道状 物体产生不正常的扭曲导致失效，从而无法得到反映内腔的客观准确的三维拉直 图像。
发明内容
本发明针对现有技术存在的上述不足，提供一种基于重力场模型的管道状物 体拉直的图像处理方法，对于管道状物体弯曲较大的区域，以及管道内部由于凸 起、凹陷等原因造成的不对称结构区域，能够进一步获得优化的管道状物体拉直 后的图像，避免了现有技术中基于电场模型和中心路径拉直方法所造成的拉直算 法失效问题。
本发明是通过以下技术方案实现的，本发明包括以下步骤第一步、通过计算机断层扫描成像法、工业计算机断层扫描成像法、核磁共 振成像法或超声成像法获得原始三维图像数据，然后对原始三维图像数据采用滤 波去噪法、对比度调整法或插值法进行优化预处理，最后通过区域增长或水平集 方法进行图像分割，得到管道状物体的内壁三维图像数据和外壁三维图像数据。
所述的内壁三维图像数据是指通过区域增长或水平集方法获得的管道状物 体的内部空腔三维图像数据。
所述的外壁三维图像数据是指利用获取的内壁三维图像数据，通过膨胀腐 蚀或水平集方法获得的管道状物体的内部空腔三维图像数据、管道状物体的管壁 三维图像数据的总和。
第二步、将外壁三维图像数据划分为结构化网格、扫略网格或自由网格中的 一种并生成外壁三维有限元模型。
所述的划分是指通过Medial Axis算法或Advancing Front算法进行划分处理。
所述的有限元网格的网格单元为四面体单元或六面体单元。
第三步、对外壁三维有限元模型施加边界条件和重力场载荷。
所述的边界条件是指在外壁三维有限元模型的一端取一横截面，设定位于该 横截面上的所有节点的位移和速度均为0。
所述的重力场载荷是指设定垂直于施加边界约束的横截面，且指向外壁三维 有限元模型内侧的方向为重力场方向。
第四步、依据所施加的边界约束和重力载荷对外壁三维有限元模型进行非线 性有限元变形计算求解，生成拉直变形后的外壁三维有限元模型。
所述的非线性有限元变形计算求解，是指列出虚功方程，使用Newton-Raphson方法、修正的Newton-Raphson方法、准Newton方法等迭代算法中的一种 计算虚功方程的极值点。
所述的虚功方程是指:变形体中任意满足平衡的力系在任意满足协调条件的 变形状态上作的虚功等于零，即体系外力的虚功与内力的虚功之和等于零。
所述的协调条件是指在变形体内部连续可导。
第五步、使用基于特征点的非刚体配准的方法实现内壁三维图像数据的拉直 变形将外壁三维有限元模型上的各节点作为源特征点(即内壁三维图像数据的特征点)，将第四步中计算求解得到的拉直变形后的外壁三维有限元模型上的各节 点作为目标特征点(即拉直变形后的内壁三维图像数据的特征点)，根据源特征点 与目标特征点的对应关系，应用非刚体配准方法，得到内壁三维图像数据中任意 一点到拉直变形后的内壁三维图像数据中对应点的映射关系，从而得到拉直变形 后的内壁三维图像数据。
所述的非刚体配准方法是指薄盘样条法、近似薄盘样条法、体积样条法或 弹性体样条法中的一种。
第六步、采用中心路径生成算法计算外壁三维图像数据的中心路径和拉直变 形后的外壁三维图像数据的中心路径。
所述的中心路径生成算法是指根据拓扑学细化法、距离变换法、水平集方法 或手工标定法得到中心路径的方法。
所述的拉直变形后的外壁三维图像数据是指利用获取的拉直变形后的内壁 三维图像数据，通过膨胀腐蚀或非刚体配准等方法中的一种获得的拉直变形后的 管道状物体的内部空腔三维图像数据、拉直变形后的管道状物体的管壁三维图像 数据的总和。
所述的中心路径是指管道状物体内部的一根中轴轨迹，该中轴轨迹上各点都 满足最优中心点判定准则。
第七步、基于第五步中得到的特征点进行旋转失真校正处理根据第五步中 得到的源特征点与目标特征点的对应关系，计算拉直变形后的内壁三维图像数据 相邻数据层的旋转角度，对拉直变形后的内壁三维图像数据第二层以后的每一层， 进行旋转变换，得到旋转失真校正后的内壁三维图像数据。
第八步、基于第五步中得到的特征点和第七步的结果，进行长度失真校正处 理，将长度失真校正后的内壁三维图像数据重新绘制，以三维形式显示出来。观 察视角位于管道内部或外部，并可由用户任意选取。
长度失真校正的方法包括重采样方法和基于特征点的非刚体配准方法。
所述的重采样方法是指根据第五步中得到的源特征点与目标特征点的对应 关系，保持外壁三维图像数据的中心路径与长度失真校正后的外壁三维图像数据 的中心路径长度相同，对旋转失真校正后的内壁三维图像数据重采样，生成长度 失真校正后的内壁三维图像数据。所述的基于特征点的非刚体配准方法是指根据外壁三维图像数据的中心路 径与长度失真校正后的外壁三维图像数据的中心路径长度相同的条件，重新计算 第五步中得到的源特征点与目标特征点的对应关系，生成用于长度失真校正的特 征点及对应关系，并据此应用非刚体配准方法，得到旋转失真校正后的内壁三维 图像数据中任意一点到长度失真校正后的内壁三维图像数据中对应点的映射关 系，生成长度失真校正后的内壁三维图像数据。
本发明相对于现有的管道状物体拉直的图像处理方法的有益效果在于1、对 于管道状物体弯曲较大(中心路径曲率较大)的区域，避免了由于同种电荷的互 斥作用使得基于电场模型的曲截面呈现不对称的分布，导致管道状物体拉直后在 该区域产生严重失真。2、对于管道内部由于凸起、凹陷等原因造成的不对称结构 区域，避免了基于中心路径拉直的方法会使管道状物体产生不正常的扭曲导致失 效的问题。
图1为实施例原始三维图像数据。
图2为实施例采用现有技术处理后示意图。
图3为实施例采用基于重力场模型的拉直方法变形后示意图。
具体实施例方式
下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下 进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限 于下述的实施例。
本实施例具体通过以下步骤实施
第一步、对己获得的一个管道状物体CT断层图像数据进行滤波处理，图1所 示，使用三维区域生长算法从断层图像中分割出粗略的管道状物体内部空腔。再 使用三维水平集分割算法，将区域生长算法的分割结果作为初值，将原断层图像 作为速度图像，进一步分割出精确的物体内部空腔，即内壁三维图像数据。使用 方向梯度方法从原断层图像中提取出管道状物体的内外表面轮廓，此轮廓作为速 度图像，内壁三维图像数据作为初值，再次使用三维水平集分割算法提取出物体 的外壁三维图像数据。
第二步、对外壁三维图像数据划分有限元网格，生成外壁三维有限元模型。选择四面体网格单元划分结构化网格。赋值单元的密度属性、杨氏模量、泊松比。 其中泊松比可设为0，以保证拉直过程不会影响径向宽度。
第三步、对于外壁三维有限元模型的一端，取一横截面，位于截面上的所有 节点施加ENCASTRE约束，即定义其位移和速度均为0 。施加全局重力场，重力加 速度取9.8m"2，场方向为垂直于施加边界约束的截面指向模型内侧的方向。
第四步、对外壁三维有限元模型进行非线性有限元变形计算求解。设管道状 物体的占有域为/2，其边界为厂，对于参考构型中的材料点X，虚功(率)方程 为<formula>formula see original document page 10</formula>
其中，^ 为总虚功率，^mt为内部虚功率，^^'为由物体外力b(X,0和指定的
面力i(X,O所引起的外部虚功率，p(X,O为密度，"③为维数，v(X,O为材料点
的位置矢量的变化率，即速度，tr(X，/)为Cauchy应力，D(Xj)为变形率。使用
Newton-Raphson迭代算法计算虚功方程的极值点，求得拉直变形后的外壁三维有 限元模型及各节点的位置坐标。
第五步、将外壁三维有限元模型上的各节点作为源特征点(即内壁三维图像 数据的特征点)，将第四步中计算求解得到的拉直变形后的外壁三维有限元模型上 的各节点作为目标特征点(即拉直变形后的内壁三维图像数据的特征点)，根据源 特征点与目标特征点的对应关系，求得每一对特征点之间的变换向量。再对内壁 三维图像数据分段利用薄盘样条插值，得到内壁三维图像数据中任意一点到拉直 变形后的内壁三维图像数据中对应点的映射关系，使用基于特征点的非刚体配准 的方法实现内壁三维图像数据的拉直，得到拉直变形后的内壁三维图像数据。
第六步、使用距离变换对外壁三维图像数据的每一点求出该点到外壁三维图像 数据边界点的最小距离。采用Dijkstra最短路径法生成外壁三维图像数据的中心 路径即在外壁三维图像数据中选取中心路径的起点和终点，然后以外壁三维图 像数据中每一点到外壁三维图像数据边界点的距离值的倒数作为权值，进行最短 路径搜索，取所有权值最小的点连接起点和终点得到外壁三维图像数据的中心路 径。同样使用距离变换和Dijkstra最短路径法求得拉直变形后的外壁三维图像数据的中心路径。
第七步、基于第五步中得到的特征点和第六步中得到的中心路径进行旋转失 真校正处理，具体步骤包括
1 )对于拉直变形后的外壁三维图像数据的中心路径上的每一点 C,(/ = 1，2，L ,n)，在拉直变形后的外壁三维有限元模型中选择与该点距离最小的一
个节点iV,;
2) 平移拉直变形后的外壁三维图像数据的中心路径，对于平移后的中心路径 上的每一点C》、1,2,L ，")，再次在拉直变形后的外壁三维有限元模型中选择与该 点距离最小的一个节点f;
3) 在外壁三维有限元模型中，可以分别找到唯一一个与W,和《对应的节点
4) 用四组节点中序号相同的节点可以构建两组向量，5 ^和
UU1 ULI LU UUl
分别计算t/,与下一相邻向量c/,+i之间的角度《、c/,'与下一相邻向量
UULt
《+1之间的角度《；
5) 计算角度差
=《-《，z. = l,2,L ,w-1, (2) 作为中心路径点C,+i所在数据层的旋转失真的估计；
6) 对拉直变形后的内壁三维图像数据第二层以后的每一层，进行旋转变换， 得到旋转失真校正后的内壁三维图像数据对拉直变形后的内壁三维图像数据第 f(/ = l，2,L ，n)层的任意一点(x,力，与之对应的旋转失真校正后的内壁三维图像数
据点的坐标(rx,oO满足
fx, L(x — cx) cos(Zl《一)_ (少一 c_y) sin(zlD + cx」
b，O7 = S i i LL(x - cx) sin(zl《一,)一 (y - c力cos(WM ) + c_y」
其中(CJC,C力为拉直变形后的外壁三维图像数据的中心路径在第/层的坐标。
第八步、基于第五歩中得到的特征点、第六步中得到的中心路径和第七步的
if / = 1， if / > 1，
if,- = 1,
if / > 1，
(3)
(4)
11结果，进行长度失真校正处理，具体步骤包括
1) 利用第七步得到的节点集合iV—^，iV2，L ,iVj计算其中每一点
A^^ = 1,2,L ，w-l)与下一相邻节点iV^之间的直线段在拉直变形后的外壁三维图 像数据的中心路径上的投影长度《；
2) 利用第七步得到的节点集合V^7Vi,iV2，L ,《}计算其中每一点
《(A: = 1，2,L ，n-l)与下一相邻节点iVw之间的直线段在外壁三维图像数据的中心 路径上的投影长度《；
3) 计算第l)步和第2)步所得的投影长度的差
<formula>formula see original document page 12</formula>作为中心路径点C^所在数据层的长度失真的估计，使用基于特征点的非刚体配
准方法进行长度失真校正对于长度失真校正后的内壁三维图像数据的第 yt(yt = l，2，L ")层，取点(x,j，A:)作为目标特征点，取点(x，y,m)作为源特征点，其<formula>formula see original document page 12</formula>
其中X和y为拉直变形后的外壁三维图像数据的中心路径在第A层的坐标 值。根据源特征点与目标特征点的对应关系，求得每一对特征点之间的变换向量。 再对旋转失真校正后的内壁三维图像数据利用薄盘样条插值，得到旋转失真校正 后的内壁三维图像数据中任意一点到长度失真校正后的内壁三维图像数据中对应 点的映射关系，使用基于特征点的非刚体配准的方法实现长度失真校正，得到长 度失真校正后的内壁三维图像数据。
第九步、使用VTK等三维可视化工具包将长度失真校正后的内壁三维图像数 据重新绘制，以用户可控的方式显示出来，如图3所示，用户可以任意改换视角， 对管道状物体内部和外部进行观察。
与现有技术处理方式得到的图2相比，通过以上步骤，在实现管道状物体拉 直的同时，解决了高曲率区域拉直后的较大失真和不对称结构区域导致基于中心 路径拉直方法失效的问题，能够方便进行进一步可视化及其他处理操作。
权利要求
1、一种基于重力场模型的管道状物体拉直的图像处理方法，其特征在于，包括以下步骤第一步、通过计算机断层扫描成像法、工业计算机断层扫描成像法、核磁共振成像法或超声成像法获得原始三维图像数据，然后对原始三维图像数据采用滤波去噪法、对比度调整法或插值法进行优化预处理，最后通过区域增长或水平集方法进行图像分割，得到管道状物体的内壁三维图像数据和外壁三维图像数据；第二步、将外壁三维图像数据划分为结构化网格、扫略网格或自由网格中的一种并生成外壁三维有限元模型；第三步、对外壁三维有限元模型施加边界条件和重力场载荷；第四步、依据所施加的边界约束和重力载荷对外壁三维有限元模型进行非线性有限元变形计算求解，生成拉直变形后的外壁三维有限元模型；第五步、使用基于特征点的非刚体配准的方法实现内壁三维图像数据的拉直变形将外壁三维有限元模型上的各节点作为源特征点，将第四步中计算求解得到的拉直变形后的外壁三维有限元模型上的各节点作为目标特征点，根据源特征点与目标特征点的对应关系，应用非刚体配准方法，得到内壁三维图像数据中任意一点到拉直变形后的内壁三维图像数据中对应点的映射关系，从而得到拉直变形后的内壁三维图像数据；第六步、采用中心路径生成算法计算外壁三维图像数据的中心路径和拉直变形后的外壁三维图像数据的中心路径；第七步、基于第五步中得到的特征点进行旋转失真校正处理根据第五步中得到的源特征点与目标特征点的对应关系，计算拉直变形后的内壁三维图像数据相邻数据层的旋转角度，对拉直变形后的内壁三维图像数据第二层以后的每一层，进行旋转变换，得到旋转失真校正后的内壁三维图像数据；第八步、基于第五步中得到的特征点和第七步的结果，进行长度失真校正处理，将长度失真校正后的内壁三维图像数据重新绘制，以三维形式显示出来，观察视角位于管道内部或外部，并可由用户任意选取。
2、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方法，其特征是，所述的有限元网格的网格单元为四面体单元或六面体单元。
3、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方 法，其特征是，所述的边界条件是指在外壁三维有限元模型的一端取一横截面， 设定位于该横截面上的所有节点的位移和速度均为0。
4、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方 法，其特征是，所述的重力场载荷是指设定垂直于施加边界约束的横截面，且指 向外壁三维有限元模型内侧的方向为重力场方向。
5、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方法，其特征是，所述的非线性有限元变形计算求解是指列出虚功方程，使用 Newton-R邻hson方法、修正的Newton-Raphson方法或准Newton方法中的一种 计算虚功方程的极值点。
6、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方 法，其特征是，所述的虚功方程是指变形体中任意满足平衡的力系在任意满足 协调条件的变形状态上作的虚功等于零，即体系外力的虚功与内力的虚功之和等 于零。
7、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方 法，其特征是，所述的非刚体配准方法是指薄盘样条法、近似薄盘样条法、体 积样条法或弹性体样条法中的一种。
8、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方法，其特征是，所述的拉直变形后的外壁三维图像数据是指利用获取的拉直变形后的内壁三维图像数据，通过膨胀腐蚀或非刚体配准等方法中的一种获得的拉 直变形后的管道状物体的内部空腔三维图像数据、拉直变形后的管道状物体的管 壁三维图像数据的总和。
9、 根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方法，其特征是，所述的长度失真校正处理，具体步骤包括1) 利用第七步得到的节点集合iV^A^A^L ,A^计算其中每一点 A^(A; = 1，2,L ,w-1)与下一相邻节点iV^之间的直线段在拉直变形后的外壁三维 图像数据的中心路径上的投影长度《；2) 利用第七步得到的节点集合iV'^iV;,A^L ,《}计算其中每一点<formula>formula see original document page 4</formula>与下一相邻节点《+1之间的直线段在外壁三维图像数据的中 心路径上的投影长度《；3)计算第l)步和第2)步所得的投影长度的差<formula>formula see original document page 4</formula>(7)作为中心路径点(^+1所在数据层的长度失真的估计，使用基于特征点的非刚体配准方法进行长度失真校正对于长度失真校正后的内壁三维图像数据的第 A:(A: = 1,2，L ")层，取点0c,乂A:)作为目标特征点，取点(x，;;,m)作为源特征点，其中<formula>formula see original document page 4</formula>附<formula>formula see original document page 4</formula>(8)if hi.其中x和_y为拉直变形后的外壁三维图像数据的中心路径在第A:层的坐标 值，根据源特征点与目标特征点的对应关系，求得每一对特征点之间的变换向量， 再对旋转失真校正后的内壁三维图像数据利用薄盘样条插值，得到旋转失真校正 后的内壁三维图像数据中任意一点到长度失真校正后的内壁三维图像数据中对 应点的映射关系，使用基于特征点的非刚体配准的方法实现长度失真校正，得到 长度失真校正后的内壁三维图像数据。
10、根据权利要求1所述的基于重力场模型的管道状物体拉直的图像处理方 法，其特征是，所述的旋转失真校正处理，具体步骤包括1) 对于管道状物体三维数据中心路径上的每一点C,，在外壁拉直变形后的有限元模型中选择与之距离最小的节点iV,;2) 平移中心路径，对于平移后中心路径上的每一点C;，再次在外壁拉直变 形后的有限元模型中选择与之距离最小的节点《；3) 在外壁拉直变形分析前模型中，可以分别找到与iV,和《一一对应的节点4) 用四组节点中序号相同的节点可以构建两组向量，5 = $-g和uu1 uu 111 训 uul4 uui分别计算^/,与下一相邻向量[/,+1之间的角度《、《与下一相邻向5) 以角度差z^,;《-《作为中心路径点C,所在数据层的旋转失真的估计，进行旋转失真校正；6) 对拉直变形后的内壁三维图像数据第二层以后的每一层，进行旋转变换， 得到旋转失真校正后的内壁三维图像数据对拉直变形后的内壁三维图像数据第 /(/ = 1,2,L ,")层的任意一点(x，力，与之对应的旋转失真校正后的内壁三维图像 数据点的坐标(rx，o;)满足LL(x - cx) sin(z!D — (y _ cy) cos(zlD + cy」,if z' > 1,其中(cx,c力为拉直变形后的外壁三维图像数据的中心路径在第/层的坐标。if / = 1， if z' > 全文摘要
一种图像处理技术领域的基于重力场模型的管道状物体拉直的图像处理方法，包括以下步骤对管道状物体的内壁三维图像数据和外壁三维图像数据进行网格划分并生成外壁三维有限元模型；对外壁三维有限元模型施加边界条件和重力场载荷；进行非线性有限元拉直变形分析；使用基于特征点的非刚体配准的方法实现内壁三维图像数据的拉直；采用中心路径生成算法计算外壁三维图像数据的中心路径和拉直变形后的外壁三维图像数据的中心路径；进行旋转失真校正处理和长度失真校正处理。本发明能够进一步获得优化的管道状物体拉直后的图像，避免了现有技术中基于电场模型和中心路径拉直方法所造成的拉直算法失效问题。
文档编号G06T15/00GK101615299SQ20091005560
公开日2009年12月30日 申请日期2009年7月30日 优先权日2009年7月30日
发明者雷 李, 俊 赵 申请人:上海交通大学