基于NURBS曲面的人体胸腰部骨组织三维重构

基于NURBS曲面的人体胸腰部骨组织三维重构 系 统 仿 真 学 报? Vol. 20 No. 14 第 20 卷第 14 期 2008 年 7 月 Jul., 2008 Journal of System Simulation 基于 NURBS 曲面的人体胸腰部骨组织三维重构 聂文忠,叶 铭,王成焘 (上海交通大学机械与动力工程学院，上海 200240) 摘 要: 提出从中国数字化可视人体数据集中获取骨组织曲线、曲面模型的完整技术路线。首先采 用多边形逼近分割图像,轮廓跟踪获得目标组织轮廓的数据点作为型值点,采用周期B样条曲线逼 近拟合成封闭曲线,利用B样条曲线控制点影响的局部性性质,调整曲线的控制多边形,获得与目 标骨组织边界吻合的轮廓曲线,再将完整的骨组织曲线模型通过直纹面法建立NURBS曲面模型。 该技术路线对人体骨组织及其它组织器官的建模提供了一种有效实用的方法。 关键词:骨组织,轮 廓跟踪,多边形逼近,周期 B 样条,曲面重构 中图分类号: TP391.9 文献标识码:A 文章编号:1004-731X (2008) 14-3738-04 3D Reconstruction of Human Thoracolumbar Skeleton Tissue Based on NURBS Surface NIE Wen-zhong, YE Ming, WANG Cheng-tao (School of Mechanical Engineering, Shanghai Jiaotong University, Shanghai 200240, China) Abstract: A whole technology way of reconstructing curve and surface model of skeletal tissue based on the Chinese digital human data was proposed. Polygonal approximation was first applied to image segmentation and the points required from contour tracing as data result. Periodic B-spline was used to approximate the curve as closed curve. Secondly the boundary curve fitted with the objective skeletal tissue can be acquired by adjusting the controlling polygon according to the local influence property of B-spline. The NURBS surface model of total skeletal tissue was established by ruled method. This method can offer an effective and validate reconstruction method for human skeletal & other tissue and organs. Key words: skeletal tissue; contour tracing; polygon approximation; periodic B-spline; surface reconstruction 为此本文提出了从中国数字化可视人体数据集中获取 引 言骨组织曲线、曲面模型的完整技术路线。依次对每张冷冻切 人体骨骼组织的曲面重建是骨组织生物力学分析、医学 片图像进行分割处理，处理每张冷冻切片时，在分割对象的 内植物设计和计算机辅助手术等应用的基础。目前，人体数 边界上交互选择少量的特征点，采用三次有理 B 样条将系 据多数是通过 CT 或 MRI 扫描获取，尽管 CT 和 MRI 对人 列特征点拟合成封闭曲线，利用样条曲线控制点影响的局部 性性质，调整曲线的控制多边形，使其与目标骨组织边界吻 体组织的检查都有自身的特点和优势，但在成像中又都有一 合，再将完整的骨组织曲线模型通过直纹面法建立 NURBS 定的局限性。CT 对硬组织检查的分辨率很高，但对于密度 曲面模型。建立后的曲线、曲面模型可以通过 IGES 标准格 非常接近的软组织的分辨有困难，同时容易产生骨性伪影; 式文件输出到 CAD 和 CAE 软件中，为人体骨组织的有限元 而 MRI 检查对软组织有很强的显示能力，但却对骨质病变 分析和运动学、动力学分析提供模型。 及钙化病灶显示差。因此，如果想从 CT 和 MRI 来重建人 体的曲面模型，就必须将同一部位的 CT 和 MRI 两种成像 融合在一起，从而弥补其中任何一种单项检查成像的不足， 1 目标组织轮廓的获取 全面地反映正常的组织结构和异常改变。但 CT 和 MRI 的 融合是一个仍在克服的技术难关。因此，自从 1989 年美国 1.1 目标组织分割 国立医学图书馆建立了采集人体横断面 CT、MRI 和组织学 数据项目以来，世界上很多国家在开展本国的数字人研究， 本文采用第三军医大学中国数字化可视人体数据集中都是从人尸体的冷冻切片开始，这样做的优点在于，它可以 的首例男性数据，胸腰部的冷冻切片层间距为 1mm，共有 包括人体的所有信息，如:骨骼、韧带、肌肉、血管等，因 [2]570 层。对冷冻切片图像中骨组织轮廓的提取虽然可以在 此直接从人体冷冻切片重建人体骨骼及组织的三维模型是 每个断层上进行，但目标组织必须作为整体通过 3D 图像分 [1]项极有意义的工作。 割得到，因为如果对每个断层进行独立图像分割处理,容易 错误分割或遗漏目标对象位于当前断层的区域，同时也会破 坏目标组织 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 面的自然连续性，所以必须根据断层图像的切 片位置，将所有断层在 3 维空间内排列起来，用插值的方法 在切片间填充体素，构造 3 维数据场，对几何建模关键的位 [3,4]收稿日期:2007-03-24 修回日期:2007-05-29 置进行重新切片处理。 基金项目:国家自然科学基金重点资助(30530230) 作者简介:聂文忠(1971-), 男, 江西省高安市人, 博士生, 研究方向为医 学 CAD/CAE。 • 3738 • 1.2 轮廓跟踪和多边形逼近 在分割出目标组织之后, 需要对每个断层进行轮廓跟 踪、多边形逼近和 B 样条曲线逼近。对轮廓操作的图像处 理方法很多，其中 Canny 边缘检测利用非最大化抑制技术可 以得到单像素宽的边缘，但很难将得到的轮廓像素按位置顺 序连通为一条封闭轮廓;此外，轮廓线存在间断和分支,难 以完成自动连通。采用轮廓跟踪算法容易得到目标的封闭轮 廓，但存在描述轮廓的像素点数据量大的问题，如果直接进 [5]行样条逼近，则反算控制点的运算量太大;因此本文采用 b) 多边形逼近后的轮廓线 如下方法:先进行多边形逼近，压缩边缘数据量并去除边缘 图 1 冷冻切片的轮廓跟踪和多边形逼近 本身的毛刺。 2.1.2 节点矢量构成 设 p 次周期 B 样条闭曲线 C (u) 的控制多边形由 n + 1 个 2 周期 B 样条曲线逼近 互异控制顶点{P, P,…, P} 组成，那么对应的节点矢量为: 0 1 n 经多边形逼近后，对轮廓曲线有插值和逼近两种处理方 = 0, u ,…, u , u ,…, u ,{u ,…, u 0 p p +1 n n +1 n + p 法可以选择。采用 B 样条曲线逼近有如下优点:?利用 B = 1,…, uu }; n + p +1n + 2 p +1 样条曲线的局部影响性质,抑制轮廓上噪声点对整条轮廓线 曲线方程:的影响,使曲线趋于光滑，使轮廓曲线的表示更加紧凑，减 n + p 小了数据量;?能够将最后结果写成 IGES 格式文件,容易与 (1) (0 ? u ? C (u ) = N(u )Pi , p iMOD ( n +1)? i = 0 各种 CAD/CAM/CAE 软件进行数据交换。 1) 如图 2 所示，现在的任务是如何将 m + 2 个参数值合理 由于断层上目标组织的轮廓是封闭的，所以需要采用周 期 B 样条曲线对其进行逼近。设多边形逼近后得到 m +1 个 分配到节点矢量的 n + 1 个内节点区间内，由于是逼近问题， 轮廓特征数据点: {Q, Q,…, Q} ;逼近目标:在控制误差 0 1 m 只考虑 m > n ? p 的情况。参照文献[7]按以下方法确定节点 范围内，得到表示轮廓曲线的最简周期 B 样条曲线所对应 矢量中的全部内部节点: 的控制点列{P} 和节点矢量{u} 。 i i m + 2 设:浮点数 f = 2.1 参数化和节点矢量构成 n + 1 2.1.1 参数化 数据点参数化的方法有等距参数化、弦长参数对于整数 j ,取整数: i = int( j ? f ) 化和向心 取余数: α = j ? f ? i 参数化等，考虑人体骨组织的特点，这里采用向心参数法对内部节点: u = (1 ? α)u+ αu; j = 1, 2,…, n p + j i ?1 i 按此种方法构成的节点矢量可以保证每个内节点区间 u至少包含一个参数值 ;两端的未知节点确定方法如下: k u = u? 1[6] p ?t n + p +1?t 轮廓特征数据点进行参数化。设与轮廓数据点 Q 相对应的i t = 1,…, p u= 1 + u n + p +1+ t p + t参数值为 u，共有 m + 1 个参数值;考虑封闭曲线首末端点 i 重合( Q= Q)，我们增加一个参数值 u= 1 。其它各点m +1 0 m +1 uuuuu u u u m-2 m-1 m m+1 0 1 2 3 0 1 的参数值确定方法如下: m +1 uuuuun-1 uun n+1 n+p n+p+1 uuu pp+1n+2p+10 1 令: 设: d = Q? Qu= 0 0 kMOD ( m +1) k ?1? 图 2 节点分布 k =1 Q? Q k k ?1 2.2 周期 B 样条闭曲线最小二乘逼近 , k = 1, 2,…, m= u+u k k ?1 d 取逼近曲线 C (u) 上具有参数值 u( k = 0,1,…, m )的点 k C (u) 与轮廓上对应的数据点 Q之间的距离平方和构造标 k k 量目标函数: 2 2 m m n + p f = Q? C = Q? N(u)P? ? ? k k i , p k iMOD ( n = 0k = 0i =0k (u)k +1) 为了推导方便，将上式改写为: 2 m n n f = Q? N(u)P? N ? ? ? k i , p k i n +1+i , p =0i =0i =0k (u)Pk i 令: M (u) = N(u) + N(u) ; i k i , p k n +1+i , p k a) 从切片图像中提取的边界 其中当 j > n + p 时， N (u) ?0 ; j , p k 第 20 卷第 14 期 Vol. 20 No. 14 2008 年 7 月系统仿真 学 报Jul., 2008 2 m n 如果上一次的控制点列为空，逼近失败;否则，N ? N ， low now f = Q? M ? ?k i k =0i = 0 N ? ( N + N，返回(2)。 ) / 2 now low upper (u)Pk i m n n = (Q? M (u)P) ? (Q? M ?? ?k i k i k i 3 胸肋骨曲面重构 k =0 i =0i =0 (u)P) k i m m n 3.1 轮廓线的转正 = Q? Q? 2[M (u)P] ??? k k i k i 由于本文方法是单方向上的闭曲面构造,因此如果轮廓 k =0k =0 i =0 ? Q+m n n k 线起始点位置相差较大的话,曲面就会出现扭转。为此,在作 [(M (u) P) ? (M (u) ?? ? i k i i k 直纹面重构之前,需要将轮廓线起始点作对齐处理,使曲面上 k =0 =io i =0 P)]不同轮廓线的起始点位置最为接近,以消除扭转。采用最短 i 将目标函数对控制点 P求偏导数: l 距离的方法来寻找最接近的轮廓起始点，过程如下: m n m ?f (A) 将所有轮廓线按照 z 坐标排序得到轮廓线序列 C; i = ?2M (u) ? Q+ 2M (u)P? (2) l k k i k i ? ?? ?P lk =0i = ok =0 M (u)l k 为了使目标函数最小，令:(B) 从第一条轮廓线开始,对于其中的每一个数据点 Q j ?f = 0l = 0,1, 2,…, n 建立一个列表 L 和距离变量 D ,寻找下一条轮廓线所有数 j j ?P l据点中距离 Q 最近的点,并加入列表 L ,将距离累加给变量 j j 从而可将式(2)改写为: D ，依此类推,按照轮廓线顺序,从上至下, 直到最后一条轮 j n m m 廓线; [M (u)M (u)]P= M (3) l k i k i l ???i =o k =ok =0 (C) 将 D 中的最小值 D所对应的列表 L中的数据 j (u)Qmin min k k 每个 l 对应一个式(3)，所以得到反算控制点矩阵方程: 点作为各自轮廓线的起始点,并将 L中所有点组成的线作 min T (4) (M M )P = Q 为基准线。 … (u)(u) M M P? ?? ? 0 0 n 0 0 3.2 胸、肋骨直纹面构造 ?? ? ? # # fi其中 M = ; P = # ?? ? ? 给定两条一般的 NURBS 曲线 ????… M (u)M (u)P0 m n m n ? ? ?? m1 1M (u)Q+ … + M (u)Q?? 0 0 0 0 m m c (u) =R(u), i , K1 1 ?d i?? i =0 Q = # 0 ?u ? ?? m 2 2 1 ?M (u)Q+ … + M (u)Q?c(u) = R(u),?? n 0 0 n m m 2 i , K ? 2 d ii =0 由于节点矢量的每个内节点区间至少含有一个参数值上。欲在该两曲线之间生分别定义在节点矢量 u 与 u 2 1 T u，所以矩阵 (M M ) 是正定的，用高斯消去法可以求解式 k 成一张直纹面: m l(4)的方程。 p(u, v) = dR(u, v)i , j i , k ; j ,l ?? 2.3 逼近过程 i = o j = o ωN(u) N (v) 2 i , j i , k j ,l 三次 B 样条曲线( p = 3 )在节点处是 C 连续， 足以满 其中: R (u, v) = i , k ; j ,1 m lωN(u) Nr , s r , k s ,l 足描述轮廓曲线的要求，所以本文采用三次周期 B 样条曲 ??r = 0 s = 0 [8](v) 线进行逼近。为了达到轮廓曲线的最简表达，采用的算法 这里 v 向节点矢量 v = 0, 0,1,1。现在必须决定 u 参数 () 如下: 的次数 k，u 向节点矢量，控制顶点 d 与权因子 ω ，以致 i , j i , j 设最少容许控制点数起始值为 L ，最多容许控制点数 j m m p(u, 0) = dR(u) = c(u) , p(u,1) = dR(u) = c(u) 。i ,0 i , k 1 i ,1i , k 2 ?? 起始值 N= m ? 1 ; N 为上一次逼近的控制点数， upper previous i =0i = 0 U 为上一次的节点矢量， P为上一次的控制点列。 previous previous 因此，必须找到 c (u) 和 c (u) 的 k 次有理基函数表示，两者 1 2 逼近过程如下: 都有公共的节点矢量 u 和相同数量 m + 1 个控制点。 (1) 从 N = ( N + N) / 2 开始; now low upper 设 K < K ，则可使 K = K 。然后将 c (u) 的次数从 K 升 1 2 2 1 1 (2) N ? N ，U ? U，参数化并构造节阶到 K 。若 K> K，则可使 K = K。然后将 c(u) 的次数 previous now previous now 1 2 1 2 从 K升阶到 K。升阶的结果，部分老控制顶点及其权因子; 点矢量U now 2 (3) 最小二乘逼近得到 P; 被新控制顶点及其权因子所替代。现在，两条曲线具有了相now 同的次数 K 升阶过程也改变老节点矢量为新节点矢量。若 (4) 检查 C (u) 各数据点位置的误差是否满足要求; 新节点矢量与另一未升阶曲线的节点矢量不一样，可取两者 (5) 如果满足要求，执行(6);否则，执行(7); 的并集作为公共的节点矢量 u。通过插入节点，使两条曲线 (6) 如果 N? N?1 ，成功退出;否则，N? now low upper 有统一的节点适量 u。同时也使两条曲线有相同数量 m + 1 个 N， now N? ( N+ N) ，返回(2);，ω ，i = 0,1,…m 用于 c(u) ，新控制顶点与新权因子:d now low upper i ,0 i ,0 1 (7) 如果 ( N?N ) ? 1, 以上一次的逼近曲线退出， upper now ?3740? [9] 及 d， ω ， i = 0,1,…m 用于 c(u) ，至此得到了建立直以较快的速度为人体运动学、动力学和有限元分析提供所需 i ,1 i ,1 2 的几何模型，为人体的生物力学分析提供了极有价值的建模 纹面所要求的全部数据，人体胸腰部骨组织的曲线、曲面重 方法。 构模型如图 3 所示。 参考文献: [1] Ackerman M J. The visible human project: A resource for education [J]. Acad Med (S1040-2446), 1999, 74(6): 667-670. 张绍祥, 刘正津, 谭立文, 等. 首例中国数字化可视人体完成[J]. [2] 第三军医大学学报, 2002, 24(10): 1231-1232. 叶铭, 杨庆铭, 于力牛, 等. 断层医学图像处理技术及其在骨科手 [3] 术中的应用[J]. 上海交通大学学报, 2003, 37(1): 30-33. 聂生东, 郑斌, 李雯. CT 图像肺结节计算机辅助检测与分类系统 [4] 设计[J]. 系统仿真学报, 2007, 19(5): 935-944. 刘斌, 关亮, 关履泰. 人脸造型中的细分网格样条技术[J]. 系统 仿真学报. 2003, 15(4): 497-498. [5] Piegl L, Tiller W. The Nurbs Book [M]. 2nd ed New York: Springer-Verlag, 1997: 181-185. [6] (a) (b) Les A. Piegl, Wayne Tiller. The Nurbs Book [M]. 2nd Edition, 图 3 胸腰部骨组织的曲线曲面模型 Springer, 1996: 207-212. [7] 高鹰, 谢胜利. 复杂目标的X-样条曲面拟合[J]. 系统仿真学报, 4 结论 2001, 13(2): 163-165. [8] 本文给出了从人体冷冻切片图像提取目标组织及几何 施法中. 计算机辅助几何设计与非均匀有理 B 样条(CAGD— NURBS) [M]. 北京航空航天大学出版社, 1994: 452-453. 建模过程中的关键技术,通过轮廓跟踪、多边形逼近和非均 [9] 匀周期 B 样条逼近得到目标组织轮廓的 B 样条曲线，再通 过直纹面重构人体胸肋骨的曲面模型，并按照 IGES 标准, M Tarokh, G McDermott, S Hayati, et al. Kinematics Modeling of a [3] (上接第 3737 页) High Mobility Mars Rover. [C]// Proceedings of the IEEE 4 结论 International Conference on Robotics and Automation. Detroit, MI, USA, 1999. USA: IEEE, 1999: 本文指出，在噪声丰富的 3D 地形上，对于月球机器人 P F Muir, C P Neumann. Kinematics Modeling of Wheeled Mobile [4] 的虚拟现实仿真，按照通常的单纯基于运动学的方法难以克 Robots [J]. Journal of Robotic Systems (S0741-2223), 1987, 4(2): 服误差的事实，并分析了这种误差不可避免的根源。为了消 282-340. 除误差，作者提出绕开机理层面直接从数值层面消除误差的 Mahmoud Tarokh, Gregory J McDermott. Kinematics Modeling and [5] 思想，并利用优化的方法建立一种误差补偿机制:在运动学 Analyses of Articulated Rovers [J]. IEEE Transactions on Robotics 反解的基础上，对位姿参数进行优化，补偿各种原因产生的 (S1042-296X), 2005, 21(4): 539-552. 误差，实现车体与地面的良好交互效果。最后，选取虚拟地 Yong Chang, Dalong Tan, Hongguang Wang, et al. Kinematics [6] 形和重建的真实地形对两种算法进行对比验证，体现出优化 Analysis of a Six-Wheeled Mobile Robot. [C]// Proceedings of the 方法的优越性。在基于 VC 平台的虚拟现实中，基于运动学 2006 IEEE/RSJ International Conference on Intelligent Robots and Systems. Beijing China, 2006. USA: IEEE, 2006: 的优化方法获得了良好的仿真效果，消除了车体下陷和悬浮 S Hayati, Richard Volpe, Paul Backes, et al. The Rocky 7 rover: A [7] 于地面的现象，充分证明了该算法在月球机器人仿真中的正 Mars Sciencecraft Prototype. [C]// Proceedings of the IEEE 确性和适用性。 International Conference on Robotics and Automation, Albuquerque, 参考文献: NM, 1997. USA: IEEE, 1997: 2458-2464. [1] Burdea G. Invited Review: The Synergy between Virtual Reality and Hubbard P M. Approximating Polyhedral with Sphere for [8] Robotics [J]. IEEE Transaction on Robotics and Automation Time-critical Collision Detection [J]. ACM Transaction on Graphics (S1042-296X), 1999, 15(3): 400-410. (S0730-0301), 1996, 15(3): 179-210. [2] J C Alexander, J H Maddocks. On the Kinematics of Wheeled Mobile Martinez Salvador B, del Pobil A P, Perez-Francisco M. Very Fast [9] Robots [J]. Robotic Research (S0278-3649), 1989, 8(5): 15-26. Collision Detection for Practical Motion Planning [C]// IEEE International Conference on Robotics and Automation, 1998, 1(5): 624-629.