CN103700136B - 一种利用三变量双调和B‑spline函数进行医学体数据矢量化的方法 - Google Patents

Abstract

本发明提出一种利用三变量双调和B‑spline函数进行医学体数据矢量化的方法，提出一种基于二次规划的方法构造三变量双调和B‑spline基函数，本发明采用针对专门的函数集合计算离散拉普拉斯算子，最小化误差函数。并分析了该基函数近似满足局部条件和单位分解（Partition of Unity）条件。基于双调和B‑spline函数，本发明设计了一种新的体数据矢量化方法。采用隐函数作为矢量化表示，使用线性规划方法优化该隐函数，完成三维医学体数据矢量化。

Description

一种利用三变量双调和B-spline函数进行医学体数据矢量化 的方法

技术领域

本发明涉及一种利用三变量双调和B-spline函数进行医学体数据矢量化的方法。

背景技术

以前已有很多学者提出多种离散数据插值方法。半径基函数是一类常用的基函数（Buhmann M.D.，Buhmann M.D.:Radial Basis Functions.Cambridge University Press，New York，NY，USA，2003.），如高斯函数，薄板样条等等。许多情况下，每个基函数是相互独立的，这样无法满足单位分解条件（即所有基函数的和恒等于1）。样条函数，如三次样条，多立方体样条（polycube spline）（Wang H.，He Y.，Li X.，Gu X.，Qin H.:Polycubesplines.Comput.Aided Des.40，6(June2008)，721–733.），Voronoi样条（Voronoi spline）（Mirzargar M.，Entezari A.:Voronoi splines.IEEE Transactions on SignalProcessing58，9(2010)，4572–4582.）等能够表示曲面。Feng等（Feng P.，Warren J.:Discrete bi-laplacians and biharmonic b-splines.ACM Trans.Graph.31，4(July2012)，115:1–115:11.）注意到有限差分和拉普拉斯算子之间的关联，提出了双变量调和B-spline，并证明它是局部的并满足单位分解。但是由于三维空间中复杂的邻域结构，直接将该方法推广到三维并不适用。另一方面，三维空间中的双调和B样条函数的性质还有待研究。

以前有许多工作研究离散拉普拉斯算子，Meyer（Meyer M.，Desbrun M.，Schr¨oder P.，Barr A.H.:Discrete differential-geometry operators for triangulated2-manifolds.In Proc.VisMath(2002)，pp.35–57.）等提出二维流形上的余切作为权重的离散拉普拉斯算子。Wardetzky等(Wardetzky M.，Mathur S.，K^¨alberer F.，Grinspun E.:Discrete laplace operators:no free lunch.In Proceedings of the fifthEurographics symposium on Geometry processing(2007)，SGP’07，EurographicsAssociation，pp.33–37.)综述了曲面上离散拉普拉斯算子的构造方法。对于双拉普拉斯算子，许多方法使用迭代拉普拉斯方法构造，但是，在不规则网格中，他们的精确度受到限制，Feng等通过加入三次约束，将近似精度提高到三次，但在三维空间中，这依然不够精确。本发明提出一种具体函数集相关的离散拉普拉斯算子计算方法，专门针具体的函数集合进行优化。

发明内容

本发明解决的技术问题是：提出了一种三元双调和B-spline函数，满足局部性和单位分解，能够满足离散点插值等应用。基于该函数，本发明矢量化三维医学体数据。

本发明采用的技术方案为：一种利用三变量双调和B-spline函数进行医学体数据矢量化的方法，包括以下四个步骤：

步骤（1）构造虚拟结点：在超体素外，增加虚拟结点，构造边界点的邻接结点，辅助拉普拉斯算子离散化；

步骤（2）计算离散拉普拉斯算子：针对给定的函数集合，使用二次规划方法，计算最优化的离散拉普拉斯算子；

步骤（3）计算离散双拉普拉斯算子：在步骤（2）的基础上，针对给定的函数集合，使用二次规划方法，计算最优化的离散双拉普拉斯算子；

步骤（4）矢量化医学体数据：在前3步的基础上，本发明矢量化三维医学体数据。

本发明的原理在于：

（1）在超体素外采集更多结点并用Voro++（Rycroft C.H.:Voro++:A three-dimensional voronoi cell library in c++.Chaos:An Interdisciplinary Journal ofNon-linear Science19，4(2009)，041111.）计算三维Voronoi剖分，然后计算采样结点及其一环邻域的并集得到如果包含边界点，则采样更多的虚拟结点，重新计算Voronoi图和直到不包含边界点为止。

（2）在三元双调和B样条中，本发明取函数集合为φ为格林函数，函数p(x)在结点t_j离散化为其中h为待计算的系数向量，h_j表示h的第j个元素，最小化如下二次规划，

其中表示集合的中元素个数。其中其中表示顶点t_j的一环邻域。

（3）本发明采样有限函数集合并优化计算：

其中δi_j＝1如果i＝j，否则δi_j＝0，H为所有列向量h构成的矩阵，Hk_j表示矩阵H的第i行第j列的元素。

（4）使用如下线性规划优化矢量化超体素形状。

使得

-e_j≤n_j-n′_j≤e_j,

e_j≤λ₁|n′_j|,j＝1,...,dim(nⁱ),

其中α＝0.25，λ₁＝0.5λ₂＝0.05，nⁱ为待求系数向量。

本发明与现有技术相比的有点在于：

1、本发明通过增加虚拟结点，支持边界结点离散化，从而在任意超体素上计算离散拉普拉斯和双拉普拉斯算子。

2、本发明采用二次规划方法针对专门的函数集合优化离散拉普拉斯算子，从而达到很高的精度。

3、在拉普拉斯算子基础上，针对专门的函数集合，本发明采用二次规划方法，计算离散双拉普拉斯算子，实验证明可以达到比现有方法更好的精确度。

4、基于双调和拉普拉斯算子，矢量化三维医学体数据，能够更好的矢量化三维医学数据。

附图说明

图1基函数对比（a）Feng的方法，（b）我们的方法；

图2基函数之和对比(a）Feng的方法，（b）我们的方法；

图3矢量化超体素形状（a）原始图像，（b）矢量化后结果；

图4矢量化超体素亮度（a）原始图像，（b）矢量化后结果；

图5体绘制效果（a）原始图像，（b）矢量化后结果。

具体实施方式

1.双变量双调和B-spline回顾

双变量调和B-spline基于调和方程的Δφ_y(x)＝δ(x-y),其中Δ为拉普拉斯算子。δ为Dirac delta函数。在二维平面上，该方程的解为格林定理为：

其中n为边界的外法向量。如果点y∈Ω，则上述积分等于1，否则等于0。上述积分的离散形式为：

其中t_j是第j个结点，表示t_j的一环邻域。是结点t_j的维诺区域，v_ij是连接结点t_i和t_j的边，e_ij是垂直于v_ij的Voronoi图的边。

类似的，双调和方程Δ²φ_y(x)＝δ(x-y)在二维平面上的解为类似的格林公式的双调和形式为：

2.三变量B-spline和它的性质

在三维欧几里得空间中，散度定理为：

本发明用▽φ_y(x)代替F，得到：

因此

更进一步，如果用▽(Δφ_y(x))替换F，得到：

本发明称三变量双调和B-spline的特征函数为：

双调和方程Δ²φ_y(x)＝δ(x-y)在三维空间中的解是其中系数使得在特征函数在单位球上的积分为1。

3.虚拟结点

在一个闭合超体素（区域）中Ω中，采样结点本发明根据Voronoi图的邻接关系确定结点之间的邻接关系。一个结点t_i为边界点，如果它的Voronoi单元和超体素区域边界相交。离散的拉普拉斯算子一般通过结点的邻域计算，但边界附近的结点缺少邻域点，导致边界精度计算明显下降，为了解决这些问题，本发明引入虚拟结点来辅助离散微分算子计算。

本发明在超体素区域外采集更多结点并用Voro++（Rycroft C.H.:Voro++:Athree-dimensional voronoi cell library in c++.Chaos:An InterdisciplinaryJournal of Non-linear Science19，4(2009)，041111.）计算三维Voronoi剖分，然后计算和的一环邻域的并集得到如果包含边界点，则采样更多的虚拟结点，重新计算Voronoi图和直到不包含边界点为止。

4.离散拉普拉斯算子

利用虚拟结点，本发明可以离散特征函数，提出一种针对某一函数集合优化的离散拉普拉斯算子，离散拉普拉斯算子的一个重要性质是系数之和为0。函数p(x)在结点t_j离散化为其中h为待计算的系数向量，h_j表示h的第j个元素。为了求解函数集合在结点t_j处的

其中表示顶点t_j的一环邻域。

在三元双调和B样条中，本发明取函数集合为最小化如下二次规划：

其中表示集合的中元素个数。

5.离散双调和算子

本发明提出一种新的理算双调和算子计算方法，一个重要约束是公共Voronoi平面上的不同方向导数是相反数，因为只有方向不同。对于函数集合我们求解如下二次规划

其中H是矩阵，它的第j列是区域双调和基函数的系数，其它未涉及元素为0。对于无穷的三元B样条格林函数集合，本发明采样有限函数集合并优化计算：

其中δi_j＝1如果i＝j，否则δi_j＝0。将前面所计算的离散拉普拉斯算子代入得到离散双拉普拉斯算子，从而计算出三变量双调和B-spline函数。

6.三维体数据矢量化

对于输入的三维体数据（CT，MRI），所述数据由CT或MRI扫描设备获得数据，首先将它分割得到超体素，然后使用双调和B-spline将其矢量化。假设分割得到h个超体素，本发明使用隐函数来表示矢量超体素。如果当且仅当F_k(x)＞F_i(x),1≤i≤h,i≠k。对于每个F_i(x)，使用双调和B-spline将它分解为F_i(x)＝Φⁱ(x)Hⁱnⁱ，其中为个采样点。

初始化：当时，令F_i(x)近似等于1，否则近似等于。所以对于采样点{s₁,s₂,...,s_r}，解如下线性最小二乘问题：

其中如果且gⁱ(s_j)＝1，如果gⁱ(s_j)＝0.8，如果且gⁱ(s_j)＝0.2，如果且gⁱ(s_j)＝0。

初始化后，采用线性规划优化隐函数，对于每个超体素，本发明首先计算错误匹配的体素也就是说表示中的点被错误标记为属于其它超体素，表示属于其它超体素的点，被错误标记为属于本发明求解如下约束最优化，

使得

|n_j-n′_j|≤λ1|n′_j|,j＝1,...,dim(nⁱ),

其中n_j表示nⁱ的第j个元素，n_j'为n_j的初始值。以上优化问题可转化为如下线性规划求解：

使得

-e_j≤n_j-n′_j≤e_j,

e_j≤λ₁|n_j′|,j＝1,...,dim(nⁱ),

颜色矢量化：超体素的形状矢量化后，矢量化超体素的颜色。对于每个超体素i，本发明优化如下能量函数：

其中c(x_j)是点x_j的亮度。

Claims (4)

1.一种利用三变量双调和B-spline函数进行医学体数据矢量化的方法，所述医学体数据为CT或MRI图像，其特征在于包括以下四个步骤：

步骤(1)构造虚拟结点：在超体素外，增加虚拟结点，构造边界点的邻接结点，辅助拉普拉斯算子离散化；

步骤(2)计算离散拉普拉斯算子：针对给定的函数集合，使用二次规划方法，计算最优化的离散拉普拉斯算子；

步骤(3)计算离散双拉普拉斯算子：在步骤(2)的基础上，针对步骤(2)所述给定的函数集合使用二次规划方法，计算最优化的离散双拉普拉斯算子；

步骤(4)将体数据分割成超体素，使用三变量双调和B-spline，矢量化每个超体素的形状；所述三变量双调和B-spline的特征函数为：

双调和方程△²φ_y(x)＝δ(x-y)在三维空间中的解是其中系数使得特征函数在单位球上的积分为1；其中△为拉普拉斯算子，其中n为边界的外法向量，点y∈Ω；

Ω表示一个三维区域；△φ_y(x)表示格林函数φ_y(x)的拉普拉斯；S表示采样点；X表示三维点；δ(x-y)为Dirac delta函数；

步骤(5)使用三变量双调和B-spline，矢量化每个超体素颜色。

2.根据权利要求1所述的方法，其特征在于：步骤(1)中在超体素中采集采样结点得到集合在超体素外采集更多采样结点并用Voro++计算三维Voronoi图，然后计算上述所有采样结点及其一环邻域的并集得到如果包含边界点，则采样更多的虚拟结点，重新计算Voronoi图和直到不包含边界点为止。

3.根据权利要求2所述的方法，其特征在于：在三变量双调和B-spline中，取函数集合为φ为格林函数，函数p(x)在结点t_j离散化为其中h为待计算的系数向量，h_j表示h的第j个元素，最小化如下二次规划，

其中其中表示结点t_j的一环邻域；

其中表示集合的中元素个数。

4.根据权利要求3所述的方法，其特征在于：步骤(3)中采样有限函数集合并优化计算：

其中δ_ij＝1如果i＝j，否则δ_ij＝0，H为所有列向量h构成的矩阵，H_kj表示矩阵H的第k行第j列的元素。