MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]}}); function MyAutoRun() {    var topp=$(window).height()/2; if($(window).height()>450){ jQuery(".outline_switch_td").css({ position : "fixed", top:topp+"px" }); }  }    window.onload=MyAutoRun; $(window).resize(function(){ var bodyw=$win.width(); var _leftPaneInner_width = jQuery(".rich_html_content #leftPaneInner").width(); var _main_article_body = jQuery(".rich_html_content #main_article_body").width(); var rightw=bodyw-_leftPaneInner_width-_main_article_body-25;   var topp=$(window).height()/2; if(rightw<0||$(window).height()<455){ $("#nav-article-page").hide(); $(".outline_switch_td").hide(); }else{ $("#nav-article-page").show(); $(".outline_switch_td").show(); var topp=$(window).height()/2; jQuery(".outline_switch_td").css({ position : "fixed", top:topp+"px" }); } }); 基于轮廓形变的复杂表面重构
  软件学报  2016, Vol. 27 Issue (10): 2676-2690   PDF    
基于轮廓形变的复杂表面重构
张雅斓, 班晓娟, 刘旭, 刘茜     
北京科技大学 计算机与通信工程学院, 北京 100083
摘要: 提出了一种基于自由形变(free-form deformation,简称FFD)及外轴投影(external axes projection,简称EAP)的复杂表面重构算法.该算法以目标形状的切片轮廓作为输入数据,此后,轮廓被嵌入到高维空间有向距离场中,在此隐式空间中,算法主要分为以下3步:生成计算序列,计算序列由计算单元组成,每一个计算单元包含上下相邻的两个轮廓;根据相邻轮廓间的拓扑关系,进行外轴投影(EAP),以解决潜在的分支问题;在每个计算单元中,根据轮廓长度决定自由形变方向,并进行自由形变,根据自由形变结果,建立轮廓间顶点的一一对应关系,并以此进行表面重构.该方法具有以下特点:输入轮廓可具有任意拓扑结构;所生成表面与输入轮廓完全贴合,生成表面准确,无自我重叠,拓扑关系不发生改变;算法高度并行,执行效率高.实验结果表明,该算法可以解决复杂表面的重构问题.
关键词: 表面重构     自由形变     外轴投影     分支问题     对应关系    
Surface Reconstruction of Complex Shapes Based on Contour Deformation
ZHANG Ya-Lan, BAN Xiao-Juan, LIU Xu, LIU Xi     
School of Computer & Communication Engineering, University of Science & Technology Beijing, Beijing 100083, China
Foundation item: National Natural Science Foundation of China (61300074, 61272357, 61370131); Program for New Century Excellent Talents in Universities of Ministry of Education, China (NCET-10-0221)
Abstract: In this paper, a novel reconstruction method based on free-form deformation (FFD) and External Axes Projection (EAP) is presented to improve the surface smoothness effect of 3D reconstruction. The contours of the slices are implicitly embedded in a higher dimensional space of distance transforms. In this implicit embedding space, reconstruction is formulated as follows. First, an arrangement of the planar slices is computed to support the approach. The arrangement consists of cells, and each cell consists of two adjacent contours. Second, the branching problem is converted into one-to-one case by the external axes projection. Next, computing direction for each cell of the arrangement is decided by the length energy. Then, in each cell a B-spline based free-form deformation is used to establish the correspondence between the adjacent contours. Finally, the contours are stitched together based on the correspondence. The key advantage of such framework are:(1) it naturally deals with contours of arbitrary topology, and it preserves shape topology; (2) the established one-to-one correspondences between two adjacent contours can guarantee the surface is continuous and no intersection; and (3) this framework is highly parallel. Experimental results show that the proposed approach performs well and can handle complicated situations.
Key words: reconstruction     free-form deformation     external axes projection     branching problem     correspondence    

物体三维重构是计算机辅助设计、计算机图形学、计算机动画、计算机视觉、医学图像处理、科学计算和虚拟现实等领域的共性科学问题及核心技术.其在逆向工程、3D打印、无损探伤、医学图像、地理信息(地形重建)等领域中发挥着极其重要的作用.随着近年来数据采集及图像处理技术的发展, 物体表面数据的获取更加容易和准确, 表面重构技术已在更多的领域中发挥作用.表示物体表面的数据可分为点云和切片两种, 点云通常通过3D扫描仪获取(激光测距技术), 切片则通过断层扫描技术获取(超声、X光等).两种数据采集方式在不同的应用场景中各有优势, 且两种数据源可通过插值或等值面提取等方法相互转化.点云可以快速地获取物体的表面信息.但是, 由于受到采集方式的限制, 对于内部结构复杂的物体, 切片则具有更大的优势.本文算法以切片数据作为输入.基于轮廓间的几何特征对杂物体的表面重构进行研究.

基于切片的表面重构算法通常包括以下几个部分:(1)数据预处理; (2)对应关系的建立或者轮廓间插值; (3)表面生成.其中, 对应关系包括全局匹配关系和局部关系.全局匹配决定了轮廓间的对应关系, 如图 1所示, 由于全局对应关系受物体拓扑结构的影响较大, 通常情况下需要人为干预才能获取正确的全局匹配结果.局部匹配(又称缝合问题, tiling problem)决定了相邻轮廓中顶点的对应关系.如图 2所示, 错误的局部对应关系产生的表面不够光滑, 而且还会引起表面拓扑结构的变化, 如图 2(b)中的白色方形空洞.局部匹配问题的解决难度很大, 因为相邻轮廓的形状千差万别, 传统的基于推进法则和约束集的局部对应法不能作为通用方法, 而利用插值进行轮廓间表面构造常常会引起拓扑结构的变化.表面的生成及表面的光滑是计算机图形学中的另一个热点问题, 不是本文的主要内容.

Fig. 1 Global correspondence 图 1 全局对应关系

Fig. 2 Local correspondence 图 2 局部对应关系

本文提出了一种全新的切片重构方法, 该方法基于自由形变技术(free-form deformation, 简称FFD)可以很好地解决局部匹配问题, 并且以任意拓扑结构的轮廓作为输入, 具有很强的鲁棒性, 通过外轴投影(external axes projection, 简称EAP)技术可以很好地解决分支问题.该方法生成的表面完全贴和于轮廓, 因此生成的表面光滑、准确.

本文第1节介绍相关工作.第2节详细介绍本文算法、数据整理及序列的计算、基于自由形变的顶点匹配、帽形区域.第3节介绍分支问题的处理-外轴投影.最后第4节给出实验结果及讨论.

1 相关工作

研究者们提出了很多基于切片的表面重构算法, 其中一些方法可归类于栅格插值法, 即水平集方法, 该类方法的主要思路是对相邻切片间的空隙进行插值计算, 即计算切片间体素的水平值, 然后通过提取等值面的方法构造表面, 例如步进三角形法.Cline[1, 2]等人运用步进三角形法根据体素的水平集构造表面, 在输入的数据质量很高并且所选择的步进正方体尺度很小时, 该方法可以得到很好的效果, 并且可以自然应对分支情况.段黎明等人[3, 4]利用工业CT图像直接重建得到STL格式的三角网格模型, 保存了三角形面片的外法向量和顶点的位置信息, 能够实现对具有复杂内腔结构零件的三维模型重建.Masala等人[5]通过蒙特卡罗算法重建使得到的体元大小一致, 模型中顶点的分布近似均匀.Zagorehev[6]利用隐式曲面重构方法, 针对空间采样中不规则点进行全局优化.然而, 该方法的缺点也十分明显, 由于最终表面的形成完全依赖于体素的水平值, 因此无法控制最终表面的拓扑结构.而且步进三角形法往往会产生额外的错误表面或洞.此外, 该方法的计算开销巨大.

另外一些解决方案是完全依赖于几何特性的, 包括本文所述的方法.这类方法首先将轮廓从切片中提取出来, 并通过建立相邻轮廓上顶点的匹配关系来生成表面.最初, 基于几何特性的重构方法只能应对“一对一”的现象, 即每一层切片只有一个轮廓, 在这种简单的情况下, 研究者普遍利用最优化、推进规则、限制集等方法来获取轮廓间顶点与顶点的对应关系.Keppel[7]用环形图中的最短路径来缝合相邻轮廓.Fuchs[8]等人提出了Divide-and-Conquer算法来寻找切片排列中的瓶颈, 并通过最小面积的最优化方法来插值表面, 与Fuchs寻找重构瓶颈类似, Sloan[9]提出了一种度量重构成本的方法, 并对Fuchs的方法进行了改进.但是, 以上方法缺点十分明显, 只有当临近轮廓具有很高的相似度时才适用.Christiansen和Sederberg[10]使用局部推进法则应对临近轮廓不相似的情况, 此外, 他们还尝试使用添加桥的办法来解决分支问题, 如图 3(a)所示, 通过在上层两轮廓间添加桥, 使上层的两个轮廓融合, 使分支问题变成一对一问题.但是, 桥的添加经常会引起相邻轮廓间拓扑结构的冲突, 如图 3(b)所示.Barequet通过应用cleft技术来解决分支问题, 并通过最小面积的方式进行表面剖分[11], 除了最小面积以外, 他还尝试了一些其他的成本度量方式, 进行表面的三角化处理[12, 13], Barequet的主体思想是用cleft技术消除轮廓间相似的部分, 用动态规划技术填充轮廓间的非相似部分.Chandrajit建立了一套十分复杂的约束条件以建立顶点间的匹配关系, 并根据匹配关系进行表面剖分.Hong[14]通过弗莱纳坐标和边界融合技术成功地重建了血管系统.与Hong类似, 解立志[15]等人利用B样条对人脑部血管进行了建模, 张舜德[16]等人用分段分面表面重构方法(DS-P)对多嵌套、多分支的情况进行了研究, 很好地解决了轮廓嵌套、轮廓严重偏置等问题, 魏亦文[17]等人使用平板能结合3次样条差值的方法进行地表重构.刘等人[18]以及Barequet[19]和Bermano[20]近期将他们的工作拓展到了非平行切片序列和多区域的重构问题, 但是目前他们使用的方法主要还是体素插值法.

Fig. 3 The application of bridge in the branch 图 3 桥在分支情况中的应用

针对分支问题, Chandrajit等人将分支问题的解决方法分为4大类, 如图 4所示.4类方法的共同点均是将分支问题转化为“一对一”问题.图 4(b)表示了一种利用组合轮廓代替分支轮廓的方法[18], 用一个凸多边形的包围盒代替c1和c2, 从而使分支简化为“一对一”问题, 该方法的优点是实现简单, 缺点则是细节损失过于严重.图 4(c)所示正如我们上面所述, 通过引入桥的方法来简化问题, 缺点是容易引入拓扑冲突.近期, 如图 4(d)图 4(e)所示的方法成为研究热点, 两种方法均是通过插入一条曲线来解决分支问题, 区别是图 4(d)在底层轮廓上插入曲线, 将底层轮廓分开, 实现与上层轮廓的一一对应.图 4(e)是在两层轮廓之间插入一条曲线, 该曲线同上层所示的两个轮廓一起组合为一个整体, 从而将分支转化为“一对一”问题.实际上, 由于现有扫描技术的发展, 切片的层间距离已经和层片间的分辨率处于同一数量级上, 因此, 两种方法的实现效果区别不大.在此问题上, Barequet[19]、Chandrajit[21]和Geiger[22]等人选用straight skeleton的方法进行直线的插入来消除分支, Jeong[23]则使用颜色场在轮廓间插值引入新的轮廓来简化分支问题.

Fig. 4 Sketch map of branch reconstruction 图 4 各类分支重构示意图

2 基于自由形变和外轴投影的复杂表面重构

综上所述, 基于几何特性的表面重构方法, 重点是如何建立相邻轮廓上顶点的匹配关系.因为相邻轮廓的拓扑关系千差万别, 很难对轮廓间的差异进行分类, 以建立完善的推进法则或约束集, 所以无论是通过推进法则还是约束集来建立顶点间的匹配关系, 都很难作为通用算法, 但是, 如果相邻轮廓具有很高的相似度甚至完全相同的话, 局部间匹配关系的建立将十分简单, 且匹配也更加准确.因此, 本文提出了一种全新的重构方法, 该方法以使相同轮廓具有相同外形为目标, 以自由形变技术为核心, 很好地解决了局部匹配问题.

2.1 问题描述

目标物体为${\mathrm O}\in {{\Re }^{3}}, $n个平行切片SO相交, 则得到m(mn)个轮廓Ci, j=SjO, i=1, ..., m, j=1, ..., n.算法目标是对于给出的${{C}_{1}}...{{C}_{m}}$重构出一个物体表面$R\in {{\Re }^{3}}$, 并且R满足RSj=Ci, j, i=1, ..., m, j=1, ..., n.所重构的应属于流形, 没有自交叉等情况.

2.2 数据整理及序列的计算

在本文所述算法中, 切片集S为基本输入, 每一个切片Sj中含有至少一个轮廓Ci, j, 即切片集S已做过边界提取预处理, 并且轮廓Ci, j的顶点P是按顺序排列的(顺时针或逆时针排列), 这样的组织形式有助于进行基于顶点间对应关系的表面三角剖分.

A.轮廓嵌入有向距离场

在本文中, 设Ω为一个矩形像素区域$\Omega \subset {{\Re }^{2}}$, 一个轮廓Ci, j⊂Ω, 则(Ω, d)构成了一个度量空间, 欧几里德距离d(x, Ci, j)为空间中点x到轮廓Ci, j的距离, 有向距离场函数Φ的定义如下式:

$ {{\Phi }_{i, j}}(x)=\left\{ \begin{align} & d(x, {{C}_{i, j}}), \text{ }\{x|x=(x, y), x\in {{w}_{i, j}}\} \\ & 0, \text{ }\{x|x=(x, y), x\in \partial {{w}_{i, j}}\} \\ &-d(x, {{C}_{i, j}}), \text{ }\{x|x=(x, y), x\in \Omega \backslash {{{\bar{w}}}_{i, j}}\} \\ \end{align} \right. $ (1)

其中, wi, j⊂Ω, Ci, j=∂wi, j, 则有:

$ d(x, {C_{i, j}}) = \mathop {\inf }\limits_{y \in \partial {w_{i, j}}} d(x, y)\{ x|x = (x, y), x \in \Omega \} $ (2)

Ci, j的有向距离场记作Φi, j.

B.计算轮廓间面积覆盖率

首先定义单位阶跃函数H:

$ H(z) = \left\{ \begin{array}{l} {\rm{1}}, {\rm{if}} z \ge {\rm{0}}\\ {\rm{0}}, {\rm{if}} z < {\rm{0}} \end{array} \right.. $

当轮廓Cm, j, Cn, j+1分别属于两切片Sj, Sj+1时, 轮廓间面积覆盖率R(Cm, j, Cn, j+1)定义如式(3):

$ R({{C}_{m, j}}, {{C}_{n, j+1}})=\frac{\int\limits_{\Omega }{H(\phi ({{\Phi }_{m, j}}(x, y), {{\Phi }_{n, j+1}}(x, y)))\text{d}x\text{d}y}}{\int\limits_{\Omega }{H(\phi ({{\Phi }_{m, j}}(x, y), 1))\text{d}x\text{d}y\centerdot \int\limits_{\Omega }{H(\phi (1, {{\Phi }_{n, j+1}}(x, y)))\text{d}x\text{d}y}}} $ (3)

其中,

$ \phi (x, y) = \left\{ \begin{array}{l} 1, {\rm{if}} x > 0, y > 0\\ 0, {\rm{else}} \end{array} \right.. $

C.计算序列

本文中将所有相邻切片中轮廓间面积覆盖率达到一定阈值的轮廓组Cm, j, Cn, j+1, 认为其可以配对, 将其加入一个序列.序列由k个计算单元组成, 每一个计算单元中包括相邻切片Sj, Sj+1中的两个轮廓Cm, j, Cn, j+1, 计算单元kiCm, j, Cn, j+1满足以下条件, 见式(4):

$ R({C_{m, j}}, {C_{n, j + 1}}) = \mathop {Sup}\limits_{{{C'}_{n, j + 1}} \subset {S_{j + 1}}} R({C_{m, j}}, {C'_{n, j + 1}}){\rm{ and }}R({C_{m, j}}, {C_{n, j + 1}}) > \delta $ (4)

其中, Cm, jSj, Cn, j+1Sj+1, δ为面积覆盖率阈值, 当两轮廓覆盖面积过小时我们可以认为两轮廓没有相互覆盖, 从而不将其配对, 当两轮廓覆盖面积满足阈值δ时, 将其加入序列.

由于最后一层切片Sp不用在配对计算上, 即只有p-1层切片需要进行计算, 轮廓Cm, j可能会没有配对轮廓, 因此, 序列k的单元个数小于所有轮廓的个数q, 此处算法最坏情况为$o\left({\frac{{{p^2}}}{q}} \right).$

图 5所示, 一个Y形物体被切片S1S2所截产生3个轮廓C1C2C3, 这是一个典型的分支问题.根据上述步骤计算得到计算序列k含有两个计算单元, 分别为k1(C2, C1)和k2(C3, C1).从此可以看出, 我们通过计算单元将原本的分支问题变为“一对一”的重构问题, 关于分支的顶点间对应关系问题, 我们将在外轴投影部分详细加以介绍.

Fig. 5   图 5  

2.3 基于自由形变的顶点匹配

本节首先介绍基于自由形变的顶点匹配, 外轴投影等技术则基于自由形变的工作机理而提出, 将放到第3节介绍.本节只讨论“一对一”的重构问题, 如图 6所示.在第2.2节中已获得了计算序列, 本节将对序列中每一个计算单元ki进行自由形变, 以建立轮廓间的顶点对应关系.

Fig. 6 The matching process based on free form deformation 图 6 基于自由变形的匹配过程

自由形变技术可以使一个图形的形状发生局部变化但不影响拓扑结构, 其通过改变源形状S的有向距离场Φs使其与目标形状T的有向距离场Φt相同, 从而建立起形状ST之间的对应关系.有若干可对局部进行非刚性变换的算法, 如光流跟踪[24, 25]、平板能(thin plate splines, 简称TPS)[26, 27].此外, Snake方法等光流跟踪和自由形变两种方法可以满足我们的需求, 均可将一个轮廓变形成另一轮廓, 并且两种方法均支持由粗到细的变尺度计算策略.自由形变技术本身保证了平滑、连续、抗噪等特点, 并且还可以保证被变形轮廓的拓扑结构不发生改变, 这些特性从本质上确保了轮廓匹配时顶点间的一一对应.与自由形变技术不同, 光流跟踪技术采用预测流动趋势方式进行形变, 即使找到合适的规则和平滑函数, 光流跟踪技术也不是很适合表面重构, 因为技术本身不能确保拓扑不变性, 因此它也不能保证所建立的顶点间的匹配关系是一一对应的.

本文应用一种基于3次B样条曲线的自由形变模型, 通过最小化Sum-of-Squared-Differences(SSD)能量来驱动轮廓形变.

由于轮廓均被嵌入到有向距离场中, 有向距离场具有更高的维度, 不仅可以隐式地表达出轮廓的形状, 还可以给出任意像素距轮廓的最近距离, 因此轮廓的形变十分适合于在有向距离场中进行.自由形变的本质是将网格P嵌入到有向距离场中, 网格的节点即为变形的控制点, 形变的过程就是操控网格P中节点位移的过程.网格P的形变由网格中所有节点{Pm, n}的位移Θ构成, 如式(5)所示.因此, 将每一个节点{Pm, n}看作是B样条曲线的控制点, 将Θ位移作用在控制点上, 此时将轮廓上的顶点网格化到网格P中, 与各顶点相关的网格顶点{Pm, n}视为3次B样条的控制点.通过3次B样条插值, 轮廓将由于节点{Pm, n}的牵动发生变化, 如图 6(c)所示.

$ \Theta = \{ \delta {P_{m, n}}\} = \{ \partial p_{m, n}^x, \partial p_{m, n}^y\}, (m, n) \in [1, M] \times [1, N] $ (5)

其中, δPm, n代表Pm, n的变化量.

根据B样条的性质, 原始图像上的任何像素X在嵌入空间中可表示为公式

$ L(X) = X = \sum\limits_{k = 0}^3 {\sum\limits_{l = 0}^3 {{B_k}(u){B_l}(v)\partial {P_{i + k, j + l}}} } $ (6)

其中, 根据3次B样条曲线定义:

$ i=\left\lfloor \frac{x-\min x}{\max x-\min x}\centerdot (M-1) \right\rfloor, j=\left\lfloor \frac{y-\min y}{\max y-\min y}\centerdot (N-1) \right\rfloor $ (7)

●  $P_{i + k, j + l}^0, X=\{ x, y\} $是网格P中的控制点, P0表示原始图像.

● Bk(u)(Bl(v))是3次B样条曲线的kth(lth)基.

像素X变换后的位置可以表示为

$ L(X:\Theta ) = X + \sum\limits_{k = 0}^3 {\sum\limits_{l = 0}^3 {{B_k}(u){B_l}(v)\partial {P_{i + k, j + l}}} } . $

根据自由形变的定义, 如何使一个轮廓变形成另一个轮廓就相当于找到每个控制点Pi+k, j+1的位移Θ.例如, 欲将图 6(a)中上层轮廓(绿色)变形成下层轮廓(红色), 应首先将两个轮廓的有向距离场分别嵌入到网格P中, 如图 6(b)所示.随后, 通过梯度下降法求得目标方程的最小值(见公式(8)), 从而获取网格控制点Pi+k, j+1的位移Θ.能量方程由内部能量、式(9)和外部能量、式(10)构成, 其中, Φu代表上层轮廓的有向距离场, Φl代表下层轮廓的有向距离场.自由形变问题即可表示为$\mathop {inf}\limits_\Theta J{\rm{(}}\Theta {\rm{)}}$:

$ J(\Theta ) = {E_{data}}(\Theta ) + {E_{smoothness}}(\Theta ) $ (8)
$ {{E}_{data}}(\Theta )=\iint\limits_{\partial \Omega }{({{\Phi }_{u}}(X)-{{\Phi }_{l}}{{(L(\Theta :X))}^{2}}dX} $ (9)
$ {{E}_{smoothness}}(\Theta )=\iint\limits_{\partial \Omega }{\left( \frac{\partial \delta L{{(\Theta :X)}^{2}}}{\partial x}+\frac{\partial \delta L{{(\Theta :X)}^{2}}}{\partial y} \right)dX} $ (10)
$ \left. \begin{align} & \frac{\partial }{\partial \Theta }J(\Theta )=\frac{\partial }{\partial \Theta }{{E}_{data}}+\frac{\partial }{\partial \Theta }{{E}_{smoothness}}= \\ & \text{ }-2\iint\limits_{\partial \Omega }{({{\Phi }_{u}}(X)-{{\Phi }_{l}}(L(\Theta :X)))(\nabla {{\Phi }_{l}}(L(\Theta :X))\centerdot \frac{\partial }{\partial {{\theta }_{i}}}\delta L(\Theta :X))dX} \\ & \begin{matrix} {} & \begin{matrix} {} & {} \\ \end{matrix} \\ \end{matrix}\text{ }+2\alpha \iint\limits_{\Omega }{\frac{\partial }{\partial x}\delta L(\Theta :X)\centerdot \frac{\partial }{\partial {{\theta }_{i}}}(\delta L(\Theta :X))}+\frac{\partial }{\partial y}\delta L(\Theta :X)\centerdot \frac{\partial }{\partial {{\theta }_{i}}}(\delta L(\Theta :X))dX \\ \end{align} \right\} $ (11)

图 6(c)所示, 经过自由形变后, 上下两层轮廓外形变得完全相同, 图中网格的变化代表了网格控制点Pi+k, j+1的位移Θ.从图 6(d)中可以看出, 经过形变后可以简单地获得顶点间的一一对应关系.图 6(e)展示了最终表面.

2.4 帽形区域

帽形区域是指相邻轮廓的相对覆盖面积比较小, 两个轮廓组成了棒球帽的形状, 较长轮廓的一部分在另一轮廓中没有对应关系.传统的重构方法在面对帽形区域时会发生形状退化的现象.Barequet[11, 12]和Jeong[23]等人详细讨论了这种情况, 一些方法是将无法对应的部分填充上平面, 但是合理的方法是构造出合适的斜面来连接两个轮廓.如果直接对计算单元中的轮廓进行自由形变也无法很好地处理帽形问题, 如图 7单数行图所示, 上层轮廓经过自由形变后收缩得十分严重, 从图 7(b)图 7(c)的单数行图可以看出, 此时根据自由形变所建立的顶点间匹配关系是错误的, 因此, 最终所重构的表面也是错误的.实际上, 图 7所示的是一种分支问题, 图 7的(5, 6)表示了整个分支, 而图 7的(1, 2)和图 7的(3, 4)分别表示了分支中的一个单支(Jeong[23]将帽形区域叫做单支问题).

Fig. 7 The degradation problem of hat shaped area 图 7 帽形区域退化问题

本方法在解决帽形问题时十分简单, 只需在对轮廓执行自由形变前测量轮廓的长度, 并确保将较短的轮廓向较长的轮廓变形即可.图 7所示双数行图可以证明该方法十分有效, 如图 7(b)图 7(c)的双数行图所示, 轮廓的缩进问题被解决了, 并且除了图 7(b)中红框所示区域外, 其他部分的匹配关系也完全正确.图 7(b)中红框所示的缺口是因为在该区域中红色轮廓在形变时(因为红色轮廓较短)在绿色轮廓内部找不到相对应的顶点, 因为轮廓内部不存在顶点, 所以该区域内红色顶点均匹配到了绿色边界上, 我们将在外轴投影中解决该问题.其中, 图 7所示有两层切片, 底层切片含有两个轮廓(红色), 上层切片含有一个轮廓(绿色).图中单数行(1、3、5)所示为不进行长度检测的自由形变结果, 双数行为进行长度检测的自由形变结果.(1)、(2)为底层轮廓中较小的轮廓参与运算.(3)、(4)为底层轮廓中较大的轮廓参与运算.(5)、(6)为两个轮廓同时参与运算.图 7(b)中所标出的红块展示了对应上的缺口, 将在外轴投影中解决.图 7(c)中标出的红块, 展示出错误轮廓自由形变后溃缩到一起.图 7(c)所示的两幅图片为在没有长度检测的情况下, 上层轮廓(绿色)完全形变后与底层轮廓重合的结果, 由于形状溃缩严重因此分开展示.

3 外轴投影(external axes projection, 简称EAP)解决分支问题

在通过自由形变进行顶点匹配工作之前, 需要对每个计算单元进行分支检测, 若存在分支问题, 则需要进行外轴投影.上一节介绍了自由形变的工作机理, 从图 6可以看出, 将自由形变技术应用到表面重构中的优点突出——可以准确地找到相邻轮廓中顶点的匹配关系.但是, 由于在一个计算单元中只包含两个轮廓, 在没有其他信息的情况下, 如果直接进行自由形变, 上层轮廓在能量方程的驱使下将逐渐变形为下层轮廓, 此时分支将引起表面重叠等问题.为了更好地说明问题, 图 8表示了一个简单的分支问题(其中, 对比图 8(c)中的红色区域, 可以证明, 由于外轴投影的引入, 对应关系的缺口已得到弥补.从图 8(d)可以看出, 由于分支问题引起的重叠问题也已经解决), 上层切片包含两个简单的圆形轮廓, 下层轮廓包含一个椭圆.在不采取任何措施的情况下, 经过自由形变后, 将会产生拓扑上的重叠, 这种现象在绝大多数情况下是不可取的.此外, 即使是单分支情况, 所构表面也将会留下缺口, 如图 7(b)所示.

Fig. 8 Simple branching problem 图 8 简单的分支问题

对上述问题分析如下:自由形变的形变能量主要来自于源形状S上一点vs(x, y)在有向距离场中的值Φs(vs)与其变形后其在目标形状T的有向距离场的值Φt(L(Θ:vs))的差, 即减少式(6)的值.自由形变后, 原形状变形后与目标形状相似, Φs(vs)≈Φt(L(Θ:vs)).然而, 在分支问题中, 原形状的形变目标是与目标形状中的一部分相似, 而不是与整个目标形状相似.此时需要控制原形状的变形结果, 即需要更改目标轮廓的有向距离场Φt, 使原形状变形到与目标形状一部分相似时已使得式(6)的能量值达到最小.

3.1 外轴投影

本文利用外轴投影技术来解决上述问题, 图像的外轴推广自Vronoi图, 与Vronoi图不同的是, 图像的外轴以图像中的轮廓作为种子进行空间分割, 其定义如下.

在度量空间(Ω, d)中, Ω是一个非空的空间, d是距离函数, 轮廓(Ck)kK作为有序元组, Ck⊂Ω, K为索引.对于轮廓集C的外轴E, 满足:

$ E=\{x\in \Omega, x\notin {{w}_{i}}|d(x, {{C}_{i}})=d(x, {{C}_{j}})\text{ for all }i\ne j\}, 其中, {{w}_{i}}\subset \Omega, {{C}_{i}}=\partial {{w}_{i}}. $

图像的外轴可以代表同一层轮廓间的拓扑关系, 将相邻切片的外轴投影到目标形状上, 即对目标形状进行合理的分割, 此时重新计算目标形状的有向距离场得到$\Phi _{t}^{e}, $将其替换$\Phi _{t}^{{}} $进行自由形变计算, 即可实现原形状形变后只与目标形状的一部分相似的功能.

设原形状$C_{1}^{1}$在切片S1的轮廓集$C_{1}^{*}$中, 且$C_{1}^{*}$中包含n个轮廓$(C_{1}^{1}, C_{1}^{2}, ..., C_{1}^{n}), $目标形状$C_{2}^{1}$${{S}_{2}}$的轮廓集$C_{2}^{*}$中, $C_{2}^{*}$中含有m个轮廓$(C_{2}^{1}, C_{2}^{2}, ..., C_{2}^{n}), $在进行外轴投影时, 只计算集合$ C_{1}^{related}$的外轴, 其满足:

$ C_{1}^{related}=\{C\subset C_{1}^{*}|R(C_{_{1}}^{i}, C_{2}^{1})>\delta, i<n\}, $

其中, $R(C_{_{1}}^{i}, C_{2}^{1})>\delta $与式(1)中的定义相同, 因为将上层切片的外轴投影到目标形状后, 将对目标轮廓产生影响, 因此只需计算与$ C_{2}^{1}$相交的$ C_{1}^{related}$的外轴即可, 如果参与计算外轴的轮廓集不够准确, 则将直接影响对目标轮廓的分割结果, 严重影响自由形变结果.

图 9所示, 图 9(2)说明了外轴投影的过程.图 9(b)中的(2)所示为经过外轴投影后的自由形变结果, 图中红色点线为外轴投影结果, 对比图 9(c)中的红色区域内的对应关系, 可以证明, 由于外轴投影的引入, 对应关系的缺口已得到弥补.从图 9(d)可以看出, 由于分支问题引起的重叠问题也已经解决.

Fig. 9 Real branch data. The red frame section indicates the function area of the outer projection 图 9 真实的分支数据.图中红框部分标明了外轴投影的作用区域

3.2 算法总结

本文所述算法以切片作为基本输入, 且假设切片已进行过预处理操作, 即已完成轮廓提取及顶点排序的工作.对物体表面进行重构, 首先需要计算出执行过程中的计算序列, 计算序列由计算单元组成(见第3.2节); 随后, 对每一个计算单元进行帽形区域检测(见第3.4节); 对于存在分支的计算单元需要进行外轴投影(见第3.5节), 以消除分支带来的重叠及顶点对应缺口等问题; 最后对计算单元进行基于自由形变的顶点匹配计算, 得出轮廓间顶点的对应关系.根据顶点间的对应关系可以准确地构造出物体表面.此外, 从图 10可以看出, 图 10(b)~图 10(f)均在各计算单元中完成, 其各个计算单元不相互影响, 因此, 本算法具有极高的并行度.

Fig. 10 Algorithmic process 图 10 算法流程

4 实验结果

这里, 展示算法结果的程序由C++编写, 直接按算法中点与点之间的对应关系进行三角剖分, 生成最终表面.对比实验使用VTK中的Marching Cube算法, 由于Marching Cube算法只针对于体数据, 因此使用vtkVoxelModeller进行图像的体素化.本实验的输入数据为ply文件, 对比实验中双方使用的数据相同.

图 11展示了兔子和龙的图像.图 11(1)的兔子的原始文件有362 272个点.图 11(b)中的(1)为本算法实验结果, 输入数据是原始文件中抽取的200层数据(实际使用198层数据, 个别抽取轮廓过短而舍去)共105 038个点.图 11(c)中的(1)为Stanford大学官网展示的重构结果, 本算法重构所用数据量不到其三分一.图 11中的(2)展示了龙的重构结果, 龙的原数据为2 748 318个顶点.图 11(b)中的(2)为本算法实验结果, 输入数据是原始文件中抽取的400层数据(实际使用数据389层), 共572 655个顶点.图 11(c)中的(2)为Stanford大学官网展示的重构结果, 本算法重构所用数据量为其20.8%.

Fig. 11 The reconstruction results of Stanford bunny and Stanford dragon 图 11 Stanford bunny及Stanford dragon展示

图 12展示了在相同数据的情况下, 本算法与Marching Cube算法进行的实验对比, 为了保证实验的准确性, 此处两种算法输入的文件相同, 且Marching Cube算法完全基于VTK来实现.

Fig. 12 The comparison between our method and Marching Cub method 图 12 本算法与Marching Cube算法的实验对比

图 12(a)中的(1)展示了泰式雕像的重构对比, 其原始数据点为19 400 000个, 本次实验数据从中抽取200层, 共847 414个顶点, 包含轮廓358个.本文算法实际使用198层.Marching Cube算法使用的vtkVoxelModeller的采样维度为(100, 100, 100), 使用的水平值为0.1.

图 12(b)展示了弥勒佛的重构结果, 该原文件共4 586 124个顶点, 本实验从中抽取200层, 共280 801个顶点, 包含轮廓369个.本算法实际使用198个点, Marching Cube算法的体数据采样维度为(100, 100, 100), 使用的水平值为0.1.

图 12(c)展示了Kiss的重构对比, 其扫描数据点数不详, 本次实验数据来自对其三维轮廓的离散化提取, 共400层, 125 528个顶点, 包含轮廓526个.本文算法实际使用398层.Marching Cube算法使用的vtkVoxelModeller的采样维度为(100, 100, 100), 使用的水平值为0.1.

图 12(d)展示了Lucy的重构对比, 其原始数据为58 241 932个顶点, 本次实验数据从中抽取200层, 1 422 873个顶点, 包含轮廓354个.本文算法实际使用198层.Marching Cube算法使用的vtkVoxelModeller的采样维度为(100, 100, 100), 使用的水平值为0.1.

图 13展示了在相同数据的情况下, 本算法与Crust表面重构方法进行的实验对比.Crust算法采用Delaunay的方法进行表面三角剖分.为了保证实验的准确性, 此处两种算法输入的文件相同, 图 13(a)展示了心脏的重构结果, 图 13(b)展示了矿石的重构结果, 从其结果不难看出, 本文的算法可以很好地解决分支问题.

Fig. 13 The comparison between our method and point-cloud method 图 13 本算法与点云的实验对比

此外, 本算法的计算集中在计算单元中, 且各个计算单元的数据不相互依赖, 因此本算法具有很高的并行度.我们对图 11所示的兔子进行了并行性验证.兔子的原始文件有362 272个点, 输入数据是原始文件中抽取的200层数据(实际使用198层数据, 个别抽取轮廓过短而舍去)共105 038个顶点.该实验在本算法下, 运行需要115.905s, 非并行情况下运行需要469.707s.

5 总结

本文描述了一种基于切片的表面重构算法, 该算法创新性地将自由形变技术应用到复杂表面的重构技术中.本算法使用计算序列的方式, 将轮廓组(上下有重叠面积的轮廓)插入到计算单元以备计算, 但是, 由于计算单元中只有两个轮廓, 无分支信息, 因此引入外轴投影技术来控制轮廓的变形范围, 以应对分支问题.本算法可以简单地解决帽形区域问题.实验结果表明, 本算法对输入数据的数据密度要求较低, 使用稀松的数据源仍可以获得较好的物体表面.此外, 对算法的并行度进行了验证.

本算法可以构造出光滑、准确的表面, 所用的自由形变技术符合由粗到细的算法策略.此外, 本算法具有一定的局限性, 其准确程度取决于预处理中边缘提取的准确程度, 在未来我们将会作进一步改进.

参考文献
[1] Cline HE, Lorensen WE, Ludke S, Crawford CR, Teeter BC. Two algorithms for the three-dimensional reconstruction of tomograms. Medical Physics, 1988, 15 (3) :320–7. [doi:10.1118/1.596225]
[2] Lorensen WE, Cline HE. Marching cubes:A high resolution 3D surface construction algorithm. Computer Graphics, 1987, 21 (4) :163–169. [doi:10.1145/37402.37422]
[3] Duan LM, Yan ZM, Chen Z, Gu MH. Arbitrary planes rendering algorithm of industrial CT serial image. High Power Laser and Particle Beams, 2013, 25 (11) :3040–3044(in Chinese with English abstract). [doi:10.3788/HPLPB20132511.3040]
[4] Duan LM, Bai Y, Wang WL, Wang Y, Wang HY, Liu FL. Simplification with feature preserving for mesh model reconstructed by industrial CT serial images. High Power Laser and Particle Beams, 2014, 26 (11) :183–188(in Chinese with English abstract). http://www.cnki.com.cn/Article/CJFDTOTAL-QJGY201411038.htm
[5] Masala GL, Golosio B, Oliva P. An improved marching cube algorithm for 3D data segmentation. Computer Physics Communications, 2013, 184 (3) :777–782. [doi:10.1016/j.cpc.2012.09.030]
[6] Zagorchev LG, Ardeshir G. A curvature-adaptive implicit surface reconstruction for irregularly spaced points. IEEE Trans.on Visualization&Computer Graphics, 2012, 18 (9) :1460–1473. [doi:10.1109/TVCG.2011.276]
[7] Keppel E. Approximating complex surfaces by triangulation of contour lines. IBM Journal of Research & Development, 1975, 19 (1) :2–11. [doi:10.1147/rd.191.0002]
[8] Fuchs H, Kedem ZM, Uselton SP. Optimal surface reconstruction from planar contours. Communications of the ACM, 1977, 20 (10) :693–702. [doi:10.1145/359842.359846]
[9] Sloan KR, Painter J. Pessimal Guesses may be optimal:A counterintuitive search result. IEEE Trans.on Pattern Analysis & Machine Intelligence, 1988, 10 (6) :949–955. [doi:10.1109/34.9117]
[10] Christiansen HN. Conversion of complex contour line definitions into polygonal element mosaics. ACM Siggraph Computer Graphics, 1978, 12 (3) :187–192. [doi:10.1145/965139.807388]
[11] Barequet G, Sharir M. Piecewise-Linear interpolation between polygonal slices. In:Proc.of the 10th Symp.on Computational Geometry, 1994 :251–272. [doi:10.1145/177424.177562]
[12] Barequet G, Shapiro D, Tal A. History consideration in reconstructing polyhedral surfaces from parallel slices. In:Proc.of the Conf.on Visualization.IEEE Computer Society Press, 1996 :149–156. [doi:10.1109/VISUAL.1996.567804]
[13] Barequet G, Shapiro D, Tal A. Multilevel sensitive reconstruction of polyhedral surfaces from parallel slices. The Visual Computer, 2000, 16 (2) :116–133. [doi:10.1007/s003710050201]
[14] Hong QQ, Li QD, Tian J. Implicit reconstruction of vasculatures using bivariate piecewise algebraic splines. IEEE Trans.on Medical Imaging, 2012, 31 (3) :543–553. [doi:10.1109/TMI.2011.2172455]
[15] Xie LZ, Zhou MQ, Wu ZK, Tian Y. Cerebrovascular expression model based on ball B-spline curves and tree structure. Journal of Computer-Aided Design & Computer Graphics, 2013, 25 (4) :541–549(in Chinese with English abstract). http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201304015.htm
[16] Zhang SD, Lu BH. 3D Triangulation surface reconstruction from complex slicing contours by DS-P method. Journal of Computer-Aided Design & Computer Graphics, 2011, 23 (2) :339–349(in Chinese with English abstract). http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201102018.htm
[17] Wei YW, Wang YC. A hybrid interpolation for reconstructing near-surface model. Journal of Computer-Aided Design & Computer Graphics, 2012, 24 (4) :466–470(in Chinese with English abstract). http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201204008.htm
[18] Liu L, Bajaj C, Deasy JO, Low DA, Ju T. Surface reconstruction from non-parallel curve networks. Computer Graphics Forum, 2008, 27 (2) :155–163. [doi:10.1111/j.1467-8659.2008.01112.x]
[19] Barequet G, Vaxman A. Reconstruction of multi-label domains from partial planar cross-sections. In:Proc.of the Computer Graphics Forum.Blackwell Publishing Ltd., 2010 :1327–1337. http://www.researchgate.net/publication/220505984_Reconstruction_of_Multi-Label_Domains_from_Partial_Planar_Cross-Sections
[20] Bermano A, Vaxman A, Gotsman C. Online reconstruction of 3D objects from arbitrary cross. ACM Trans.on Graphics, 2011, 30 (5) :1317–1329. [doi:10.1145/2019627.2019632]
[21] Bajaj CL, Coyle EJ, Lin KN. Arbitrary topology shape reconstruction from planar cross sections. Graphical Models & Image Processing, 1997, 58 (6) :524–543. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.53.2267
[22] Geiger B. Three-Dimensional modeling of human organs and its application to diagnosis and surgical planning. Technical Report, Institute National de Recherche Informatique et Automatique, 1993 :2105. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.56.9652
[23] Jeong J, Kim K, Park H, Cho H, Jung M. B-Spline surface approximation to cross-sections using distance maps. The Int'l Journal of Advanced Manufacturing Technology, 1999, 15 (12) :876–885. [doi:10.1007/s001700050145]
[24] Paragios N, Rousson M, Ramesh V. Non-Rigid registration using distance functions. Computer Vision & Image Understanding, 2003, 89 (2-3) :142–165. [doi:10.1016/S1077-3142(03)00010-9]
[25] Xu C, Pilla JJ, Isaac G, Gorman Iii JH, Blom AS, Gorman RC, Ling Z, Dougherty L. Deformation analysis of 3D tagged cardiac images using an optical flow method. Journal of Cardiovascular Magnetic Resonance, 2010, 12 (4) :508. [doi:10.1186/1532-429X-12-19]
[26] Belongie S, Malik J, Puzicha J.Matching shapes.In:Proc.of the 8th IEEE Int'l Conf.on Computer Vision (ICCV 2001), Vol.1.2001.IEEE, 2001.454-461.[doi:10.1109/ICCV.2001.937552]
[27] Roth HR, McClelland JR, Hampshire TE, Boone DJ, Hu YP, MarcModat, Zhang H, Ourselin S, Halligan S, Hawkes DJ, Roth H.Registration of prone and supine CT colonography datasets with differing endoluminal distension.In:Abdominal Imaging, Computation and Clinical Applications.Berlin, Heidelberg:Springer-Verlag, 2013.29-38.[doi:10.1007/978-3-642-41083-3_4]
[3] 段黎明, 颜志明, 陈中, 谷明辉. 工业CT切片序列任意方向剖面的绘制方法. 强激光与粒子束, 2013,25 (11) :3040–3044. [doi:10.3788/HPLPB20132511.3040]
[4] 段黎明, 白洋, 王武礼, 王勇, 王浩宇, 刘丰林. 基于工业CT图像重构的网格模型的保特征简化方法. 强激光与粒子束, 2014,26 (11) :183–188. http://www.cnki.com.cn/Article/CJFDTOTAL-QJGY201411038.htm
[15] 解立志, 周明全, 武仲科, 田沄. 球B样条曲线结合树状结构的脑血管表示模型. 计算机辅助设计与图形学学报, 2013,25 (4) :541–549. http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201304015.htm
[16] 张舜德, 卢秉恒. 复杂断层轮廓集分段分面三角化表面重构. 计算机辅助设计与图形学学报, 2011,23 (2) :339–349. http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201102018.htm
[17] 魏亦文, 王彦春. 混合插值法重构近地表模型. 计算机辅助设计与图形学学报, 2012,24 (4) :466–470. http://www.cnki.com.cn/Article/CJFDTOTAL-JSJF201204008.htm