基于最优质量传输的Focus+Context可视化

doi:10.3785/j.issn.1008-9497.2023.06.003

基于最优质量传输的Focus+Context可视化

苏科华^,^,¹, 刘百略¹, 雷娜^,^,², 李可涵¹, 顾险峰³

1.武汉大学计算机学院, 湖北武汉 430072

2.大连理工大学国际信息与软件学院, 辽宁大连 116024

3.纽约州立大学石溪分校计算机系, 美国石溪 11790

Focus+Context visualization based on optimal mass transportation

SU Kehua^,^,¹, LIU Bailüe¹, LEI Na^,^,², LI Kehan¹, GU Xianfeng³

1.School of Computer Science，Wuhan University，Wuhan 430072，China

2.International School of Information Science & Engineering，Dalian University of Technology，Dalian 116024，Liaoning Province，China

3.Department of Computer Science，State University of New York at Stony Brook，Stony Brook 11790，New York，USA

通讯作者: ORCID：https：//orcid.org/0000-0003-3361-0756，E-mail：nalei@dlut.edu.cn.

收稿日期: 2023-06-12 修回日期: 2023-07-19 接受日期: 2023-07-31

基金资助:

国家自然科学基金资助项目. 62272354

Received: 2023-06-12 Revised: 2023-07-19 Accepted: 2023-07-31

作者简介 About authors

苏科华（1979—），ORCID：https：//orcid.org/0000-0002-1384-9762，男，博士，教授，主要从事计算机图形学研究，E-mail：skh@whu.edu.cn. , E-mail：skh@whu.edu.cn

摘要

在分辨率有限的显示设备上，Focus+Context技术可用于大型复杂模型的可视化。提出了一种基于最优质量传输的Focus+Context可视化方法。通过最优质量传输映射，对自身进行体积变形，将源测度（体素）转换为传输成本最小的目标测度；将求解最优质量传输问题等价于凸优化过程，转换为计算几何中经典的幂Voronoi图计算。与现有方法相比，本文方法具有坚实的理论基础，保证了解的存在性、唯一性和平滑性；允许用户精确控制目标测度，选择多个形状不规则的聚焦区域，使产生的变形是全局平滑的，并可自由翻转。用于自医学应用和科学仿真的几项体数据集，证明了所提方法是有效和高效的。

关键词： 辐射度 ; 全局光照 ; 常量时间

Abstract

In visualization field, Focus+Context techniques have been developed to visualize large, complex models on the display device with limited resolution. In this work, we propose a novel method for Focus+Context visualization based on optimal mass transportation. An optimal mass transportation map deforms a volume to itself, transforms the source measure (volumetric element) to the target measure with the minimal transportation cost. Solving the optimal mass transportation problem is equivalent to a convex optimization, and can be converted to computing power Voronoi diagrams in classical computational geometry. Comparing to existing approaches, the proposed method has solid theoretic foundation, which guarantees the existence, uniqueness and the smoothness of the solution. It allows the user to accurately control the target measure, and select multiple focus regions with irregular shapes. The resulting deformation is globally smooth and flipping free. Experiments with several volume data sets from medical applications and scientific simulations demonstrate the effectiveness and efficiency of our method.

Keywords： radiosity ; global illumination ; constant time

PDF (3910KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

苏科华, 刘百略, 雷娜, 李可涵, 顾险峰. 基于最优质量传输的Focus+Context可视化. 浙江大学学报(理学版)[J], 2023, 50(6): 681-691 doi:10.3785/j.issn.1008-9497.2023.06.003

SU Kehua, LIU Bailüe, LEI Na, LI Kehan, GU Xianfeng. Focus+Context visualization based on optimal mass transportation. Journal of Zhejiang University(Science Edition)[J], 2023, 50(6): 681-691 doi:10.3785/j.issn.1008-9497.2023.06.003

0　引言

Focus+Context（F+C）框架可在有限分辨率屏幕上对大型复杂模型进行可视化，通过模拟光学镜头或变形视体放大感兴趣区域，并允许用户在不丢失模型形状和拓扑结构的全局视图下查看感兴趣区域的细节，在放大感兴趣区域的同时不会剪掉其他区域。这种可视化技术可将模型的整体感知传达给用户，并引导用户将注意力集中在局部区域。

0.1　空间变形法

基于网格的空间变形是F+C可视化最常用的技术之一^［1］。体表示一个规则的网格，并被划分为具有不同重要性的区域，赋予体素对应的重要性值。通过重新定位顶点对个体进行变形，使体素与其重要性成比例扩大，同时保持其局部形状，变形是全局平滑的。此外，变形没有翻转或自交叉。

用 $Ω \subset R^{3}$ 表示体积，其欧几里得体积元为 $d x \land d y \land d z$ 。指定函数 $μ : Ω \to R^{+}$ ，将其视为体积元的放大因子，目标体积元由 $μ (x, y, z) d x \land d y \land d z$ 给出。由于变形不超过体积，故可得自映射 $φ : Ω \to Ω$ 。显而易见，映射 $φ$ 是光滑的、一对一的、映上的，因此 $φ$ 为微分同胚。更重要的是，映射将初始的欧几里得体积元变形为目标体积元：

D φ : d x \land d y \land d z \to μ (x, y, z) d x \land d y \land d z 。

这就需要 $φ$ 的雅可比矩阵的行列式等于 $μ$ ，即

d e t (D φ) (x, y, z) = μ (x, y, z) 。

（1）

在现有方法中，式（1）中的雅可比关系可用最小二乘能量近似，首先在能量中加入惩罚项，实现微分同胚条件，然后对能量进行优化得到变形，但这些方法有以下缺点：（1）不能准确地给出放大因子。没有理论保证可以得到一个完全满足式（1）的解；（2）通过优化惩罚能量可能无法满足微分同胚条件，得到的映射仍有可能存在翻转或自交；（3）优化可能会停止在非线性能量的局部最小值，而不是全局最小值。因此，优化依赖于初始条件的选择，结果的可重复性较差。

0.2　最优质量传输法

为了克服这些缺点，本文提出了一种基于最优质量传输（optimal mass transportation，OMT）的F+C可视化方法。先假设 $Ω$ 上的测度的定义，将其视为所需的体积元，进一步假设总测度等于初始体积。满足式（1）的微分同胚体有无限多个，但仅有一个可使得传输成本最小化：

C （ φ ） : = {\int_{Ω} |p - φ (p)|}^{2} d x \land d y \land d z,

其中， $φ$ 为最优质量传输映射。此外，凸函数 $f : Ω \to R$ ，其梯度映射 $p \mapsto \nabla f (p)$ 给出了最优质量传输映射。此时，式（1）变为Monge-Amper´e方程：

d e t H e s s f (x, y, z) = μ (x, y, z),

其中，Hess f是函数f的Hessian矩阵。

在本文方法中，Monge-Amper´e方程被离散化，通过凸优化进行求解。优化是迭代的，每步的算法可归结为计算一组超平面的上包络，并投影包络以获得幂Voronoi图和对偶幂Delaunay三角剖分，通过计算几何中的经典算法实现^［2］。

与现有方法相比，OMT方法具有以下优点：（1）可以直接找到式（1）的精确解，且其存在性、唯一性有理论保证。这使得用户可以精确控制放大因子；（2） OMT映射是微分同胚的，因其雅可比矩阵处处为正；（3）寻找OMT映射的方法等价于凸优化过程。由于能量的凸性，该方法存在唯一全局最小值；结果与初始条件无关，且易复现。

第1节简要回顾F+C可视化的相关工作；第2节介绍理论背景；第3节详细介绍最优质量传输映射算法；第4节介绍F+C可视化；第5节呈现实验结果；第6节讨论算法的优点和缺点；第7节为总结和展望。

1　相关工作

F+C可视化已被用于图形、城市和地图、嵌套网络、3D模型等领域进行交互式数据探索。

F+C可视化方法之一是基于光学透镜和放大镜的仿真，例如由MANOJIT等^［3］提出的用于文本、图像和图形可视化的fisheye方法。fisheye镜头提供了模型的整体视图，虽保留了空间关系，但在边缘引入了明显的失真，无法形式化地控制焦点区域，也不能保留上下文区域的局部特征。KUMAR等^［4］对fisheye相机模型进行了统一分类处理，提供了一个自成体系的fisheye感知参考。KUMAR等^［5］开发了一种用户界面Tool glass and magic lenses，增强了感兴趣区域的焦点特征，压缩了感兴趣程度低的区域。AGLEDAHL等^［6］受文献［5］的启发，提出了一种新型的用于虚拟环境的放大工具，但是性能并没有显著提高。CARPENDALE等^［7-8］提出了几种畸变模式，如拉伸正交和非线性径向和多尺度方法，为焦点区域分配了更多空间，实现了3D视点独立畸变，通过统一展示空间，高效地利用了可用的显示空间^［9］。VIOLA等^［10］提出了基于重要性驱动的自动显示物体的F+C体绘制。利用预先确定的物体重要性，对体素可见性优先级进行编码，用于指导渲染，避免重要区域被不重要区域遮挡；另外还提出了一种自动聚焦特征的技术^［11］，从一组预定义的特征中挑选焦点，自动确定特征中最具表现力的视图。PIETRIGA等^［12］提供了一种新型的Sigma透镜，利用时间和其半透明特点实现更高效的过渡；之后，又提出了一种就地放大和独立展示的系统^［13］。

另一种方法是基于体积变形的仿真。LAMAR等^［26］使用硬件加速对渲染的2D或3D图像体积的变形过程，动态计算网格顶点的纹理坐标，并将坐标投影到同质空间上渲染纹理，得到理想的结果。MCGUFFIN等^［14］将变形技术应用于体数据浏览，允许用户打开、散开或剥开外层撕裂的小牛肉隐藏结构。CORREA等^［15］使用物理和光学图解运算符操纵数据对象的几何形状。WANG等^［16］提出了一种自由形式的体透镜函数方法高亮、曝光和非线性放大目标，同时提出通过变形突出3D体数据集中特征的解决方案。对于大型表面模型，WANG等^［17］提出了一种能量优化模型，在变形上下文区域的同时放大ROI；进一步将该模型应用于体数据集，以保留感兴趣的特征，作为F+C可视化结果^［1］。ZHAO等^［18］提出的共形放大镜方法适用于体数据集和表面网格模型，能消除局部角度失真，并以一致和全局的方式保持视觉连续性。TAO等^［19］将F+C可视化技术推广至基于网格空间变形的流场可视化领域。然而，大多数技术无法精确控制放大因子。有些方法是启发式的，没有坚实的理论基础，解的存在性、唯一性和平滑性也没有理论保证。相比之下，本文方法具有严格的数学基础，并可在体素级别上实现控制放大。

ZHAO等^［20］通过最优质量传输获得了保面积映射，并将其应用于医学图像和图形的可视化，通过生成平滑、连续的过渡体，控制放大的精确测度，且不会丢失视觉一致性上的任何上下文信息，但方法的焦点在2D表面。而本文方法主要用于体可视化，具有推广至更高维度的潜力。

2　理论背景

2.1　最优质量传输

Monge问题　设X和Y是分别具有概率测度μ和ν的测度空间。假设X和Y具有相等的总测度，为

\int_{X} μ = \int_{Y} ν ，

如果映射 $T : X \to Y$ 满足：对任何可测度的集合 $B \subset Y$ ，有

\int_{T^{- 1} (B)} μ = \int_{B} ν ，

则T为保测度映射，说明由T引出的μ的正向推量等于ν，表示为 $ν = T_{#} μ$ 。

设 $c (x, y)$ 为从 $x \in X$ 到 $y \in Y$ 的传输成本，则T的总传输成本为

E (T) : = \int_{X} c (x, T (x)) d μ (x) 。

（2）

Monge问题引出了最优质量传输问题：如何找到一个保测度映射 T， $T_{#} μ = ν$ ，使得式（2）中的传输成本最小？^［21］

KANTOROVICH^［22］引入了Monge问题的松弛，并使用线性规划解决了该问题。BRENIER^［23］证明了以下定理。

定理1　假设X，Y是 $R^{n}$ 中的子集，源X是凸域，传输成本是二次欧氏距离，即

c （ x, y ） = {|x - y|}^{2} 。

分别给定在X和Y上的概率测度μ和ν，则存在唯一的最优传输映射 $T : (X, μ) \to (Y, ν)$ ，此外还存在凸函数 $f : X \to R$ ，其可根据一个常数确定唯一性，最优质量传输映射由梯度映射 $T : X \mapsto \nabla f (X)$ 给出。

假设μ和ν是光滑的，f是二阶光滑的， $f \in C^{2} (X, R)$ ，如果f是保测度的，则其满足Monge-Amper´e方程：

d e t (\begin{matrix} \frac{\partial f}{\partial x_{1}^{2}} & \frac{\partial f}{\partial x_{1} x_{2}} & \dots & \frac{\partial f}{\partial x_{1} x_{n}} \\ \frac{\partial f}{\partial x_{2} x_{1}} & \frac{\partial f}{\partial x_{2}^{2}} & \dots & \frac{\partial f}{\partial x_{2} x_{n}} \\ ⋮ & ⋮ & ⋮ \\ \frac{\partial f}{\partial x_{n} x_{1}} & \frac{\partial f}{\partial x_{n} x_{2}} & \dots & \frac{\partial f}{\partial x_{n}^{2}} \end{matrix}) = \frac{μ}{ν \circ \nabla f} 。

（3）

一般来说，Monge-Amper´e方程是高度非线性的，采用常规的有限元方法无法求解此类型的偏微分方程。然而，基于其几何性质，可采用变分方法，通过凸优化求解。

2.2　离散最优质量传输

在实践中，采用离散点集确定目标域，以表述离散设置下的最优传输问题，当采样密度趋于无穷时，离散解收敛至光滑解。假设μ在X上有一个紧支撑，定义为

Ω = S u p p μ = \{x \in X | μ (x) > 0\},

并假设Ω是X中的凸域。空间Y通过狄拉克测度被离散为 $Y = \{y_{1}, y_{2}, \dots, y_{k}\}$ ，狄拉克测度表示为

ν = \sum_{i = 1}^{k} ν_{i} δ (y - y_{i}) 。

定义一个由k个实数组成的高度向量 $h = (h_{1}, h_{2}, \dots, h_{k}) \in R^{k}$ 。对于每个 $y_{i} \in y$ ，构造定义在X上的超平面：

π_{i} (h) : 〈x, y_{i}〉 + h_{i} = 0,

（4）

其中， $〈,〉$ 是 $R^{n}$ 中的内积。定义分段线性函数

u_{h} (x) = \underset{1 \leq i \leq k}{m a x} \{〈x, y_{i}〉 + h_{i}\},

（5）

则 $u_{h}$ 为凸函数。用 $G (h)$ 表示凸函数 $u_{h}$ 的图， $G (h)$ 是无限凸多面体，可支撑平面 $π_{i} (h)$ ，所以 $G (h)$ 是平面 $\{π_{i} (h)\}$ 的上包络线。由 $G (h)$ 的投影引出多面体划分Ω：

\begin{matrix} Ω = \cup_{i = 1}^{k} W_{i} (h), \\ W_{i} (h) : = \{x \in X | u_{h} (x) = 〈x, y_{i}〉 + h_{i}\} \cap Ω 。 \end{matrix}

（6）

单元 $W_{i} (h)$ 是凸多面体 $G (h)$ 的一个面片在Ω上的投影。单元 $W_{i} (h)$ 上的凸函数 $u_{h}$ 是线性函数，可描述为 $π_{i} (h)$ ，因此，梯度映射为

\nabla u_{h} : W_{i} (h) \mapsto y_{i}, i = 1,2, \dots, k,

（7）

可将单元 $W_{i} (h)$ 映射至点 $y_{i}$ 。

在当前离散设置约束下，定理1可以表述为

定理2（离散OMT）　对于具有凸支撑 $Ω \subset X$ 的X上的任何给定测度μ和Y上的狄拉克测度ν，有

\int_{Ω} d μ = \sum_{i = 1}^{k} ν_{i}, ν_{i} > 0,

必须存在唯一的高度向量 h，直到添加一个常数向量 $(c, c, \dots, c)$ ，由式（5）中的凸函数引出如式（6）所示的单元划分Ω，使得所有单元都满足保测度约束：

\int_{W_{i} (h)} d μ = ν_{i}, i = 1,2, \dots, k 。

（8）

此外，梯度映射 $\nabla u_{h}$ 优化了传输成本：

E (T) : = {\int_{Ω} |x - T (x)|}^{2} μ (x) d x 。

（9）

ALEXANDROV^［24］用拓扑学方法证明了解的存在性和唯一性。之后存在性也被AURENHAMMER^［25］证明。

GU等^［2］基于变分原理也给出了存在性和唯一性证明。首先，定义高度向量的可容许空间：

H : = \{h | \int_{W_{i} (h)} d μ > 0\} ⋂ \{\sum_{i = 1}^{k} h_{i} = 0\} 。

（10）

然后，定义能量 $E (h)$ 为以图 $G (h)$ 为界的凸多面体和通过Ω边界的圆柱体的体积减去一个线性项：

E (h) = \int_{Ω} u_{h} (x) μ (x) d x - \sum_{i = 1}^{k} ν_{i} h_{i} 。

（11）

能量梯度为

\nabla E (h) = (\int_{W_{i} (h)} μ - ν_{i}) 。

（12）

假设单元 $W_{i} (h)$ 和 $W_{j} (h)$ 在某表面处相交，即

f_{i j} (h) = W_{i} (h) ⋂ W_{j} (h) ⋂ Ω,

那么， $E (h)$ 的Hessian矩阵各元素（i为行数，j为列数）为

\frac{\partial^{2} E (h)}{\partial h_{i} \partial h_{j}} = \{\begin{array}{l} \frac{f_{i j} (h) μ}{|y_{j} - y_{i}|} ， W_{i} (h) ⋂ W_{j} (h) ⋂ Ω \neq \emptyset ， \\ 0 ， 其他 。 \end{array}

（13）

可容许空间H是凸的和Hessian矩阵在H上是正定的已得到证明，因此在式（11）中能量 $E (h)$ 是凸的。此外，全局最小值 $h^{*}$ ^{是H的一个内部点。在最小值处， $\nabla E (h^{*}) = 0$ ，表示梯度映射 $\nabla u_{h}$ 满足式（8）的保测度约束，而且梯度映射是最优质量传输映射。}

利用牛顿法可求得能量（式（11））的全局最小值。与文献［22］中 $k^{2}$ 个未知数相比，该方法只有k个未知数。

3　算法

源域Ω是 $R^{3}$ 中的标准单位立方体，目标是一组离散点 $Y = \{q_{1}, q_{2}, \dots, q_{k}\}$ ，这些离散点密集均匀地对单位立方体进行采样。立方体上的源域测度处处为μ=1。Y上的目标测度由用户指定， $ν = \{ν_{1}, ν_{2}, \dots ， ν_{k}\}$ 。对于每个目标点 $q_{i} \in Y$ ，先在 $R^{4}$ 上构造超平面：

π_{i} (h, p) : = 〈q_{i}, p〉 + h_{i}, i = 1,2, \dots, k,

再计算超平面的上包络线 $G (h)$ 。

3.1　幂Voronoi图和幂Delaunay三角剖分

对于每个超平面 $π_{i} (h)$ ，首先构造对偶点 $π_{i}^{*} (h) \in R^{4}$ ，假设 $q_{i} \in R^{3}$ 的坐标为 $(x_{i}, y_{i}, z_{i})$ ，则对偶点为

π_{i}^{*} (h) = (x_{i}, y_{i}, z_{i}, - h_{i}), i = 1,2, \dots, k 。

然后利用增量凸包算法^［26］，计算 ${π_{1}^{*} (h), π_{2}^{*} (h), \dots,$ $π_{k}^{*} (h)}$ 的凸包，并将得到的凸包表示为 $𝒞 (h)$ 。最后检查 $R^{4}$ 中凸包 $𝒞 (h)$ 边界上的四面体 $t_{i j k l}$ ：

t_{i j k l} (h) \in \partial 𝒞 (h) 。

该四面体的超平面方程为

|\begin{matrix} x_{i} & y_{i} & z_{i} & - h_{i} & 1 \\ x_{j} & y_{j} & z_{j} & - h_{j} & 1 \\ x_{k} & y_{k} & z_{k} & - h_{k} & 1 \\ x_{l} & y_{l} & z_{l} & - h_{l} & 1 \\ x & y & z & w & 1 \end{matrix}| = 0,

可得超平面的法向量，记为 $n_{i j k l}$ 。如果法向量 $n_{i j k l}$ 的w分量为负，则将四面体 $t_{i j k l}$ 投影至 $(x, y, z)$ -超平面。凸包 $𝒞 (h)$ 的投影过程产生了点集Y的幂次Delaunay三角剖分，记为 $𝒯 (h)$ 。

超平面 $\{π_{i} (h)\}$ 的上包络表示为 $ℰ (h)$ ，是凸包 $𝒞 (h)$ 下半部分的对偶，对于凸包 $\partial 𝒞 (h)$ 边界上的每个四面体 $t_{i j k l}$ ，其法向量是负w分量，其对偶是上包络 $ℰ (h)$ 上的一个顶点，这是4个超平面 $\{π_{i} (h), π_{j} (h), π_{k} (h), π_{l} (h)\}$ 的交集。四面体 $t_{i j k l}$ 中的每个三角形 $△_{i j k}$ 对偶于 $ℰ (h)$ 的一条边， $ℰ (h)$ 是3个超平面 $\{π_{i} (h), π_{j} (h), π_{k} (h)\}$ 的交集。四面体 $t_{i j k l}$ 中的每条边 $e_{i j}$ 对应 $ℰ (h)$ 中的一个面，是2个超平面 $\{π_{i} (h), π_{j} (h)\}$ 的交集。四面体 $t_{i j k l}$ 中的每个顶点 $ν_{i}$ 对应 $ℰ (h)$ 的一个单元，被超平面 $π_{i} (h)$ 所支撑。通过计算凸包 $𝒞 (h)$ 的对偶，得到上包络 $ℰ (h)$ ，将上包络投影至 $(x, y, z)$ -超平面，得到该超平面的幂Voronoi图，每个幂Voronoi单元与Ω相交得到Ω的幂Voronoi单元分解，记为 $𝒱 (h)$ 。

实际上，上包络 $ℰ (h)$ 是凸函数 $G (h)$ 的图，幂Voronoi图 $𝒱 (h)$ 是在式（6）中通过投影 $G (h)$ 得到的多面体划分Ω。上包络 $ℰ (h)$ 、凸包 $𝒞 (h)$ 、幂Voronoi单元分解 $𝒱 (h)$ 、幂Delaunay三角剖分 $𝒯 (h)$ 及其关系的二维类比如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 上包络 $ℰ (h)$ 、凸包 $𝒞 (h)$ 、Voronoi单元分解 $𝒱 (h)$ 、幂Delaunay三角剖分 $𝒯 (h)$ 及其关系的二维类比

Fig.1 The upper envelope $ℰ (h)$ ， the convex hull $𝒞 (h)$ ， the power voronoi cell decomposition $𝒱 (h)$ ， the power Delaunay triangulation $𝒯 (h)$ and two-dimensional analogies of their relationship

3.2　最优质量传输映射

设离散点集Y包含在单位立方体Ω中，初始高度向量设为

h_{i} = \frac{1}{2} 〈q_{i}, q_{i}〉, i = 1,2, \dots, k 。

初始幂Delaunay三角剖分 $𝒯 (h)$ 是传统的Delaunay三角剖分，幂Voronoi单元分解 $𝒱 (h)$ 是传统的Voronoi单元分解。在每一步中，首先，计算幂 Delaunay三角剖分 $𝒯 (h)$ 和幂Voronoi单元分解 $𝒱 (h)$ 。式（11）中的体积能量梯度由式（12）给出，体积能量的Hessian由式（13）给出。然后，求解线性方程

\nabla E (h) = [\frac{\partial^{2} E (h)}{\partial h_{i} \partial h_{j}}] δ h

。（14）

当线性约束 $\sum_{i = 1}^{k} h_{i} = 0$ 时，解存在且唯一。最后，用牛顿法更新高度向量：

h \leftarrow h + λ (δ h),

其中， $λ$ 为步长参数。理论上，步长参数的选择应使高度向量保持在可容许空间H（式（10））内，即在幂Voronoi单元分解 $𝒱 (h)$ 中，每个单元 $W_{i} (h)$ 是非空的。在优化过程中，允许 h 超过可容许空间H。体能的凸性自动引导高度向量返回可容许空间。

4　F+C可视化

体原始图像数据的分辨率为512×512×512，每个体素的强度范围为0~255。源域Ω是以原点为中心、边长为2的立方体。

4.1　焦点选择和目标测度约束

使用文献［1］中的方法，通过显著性人工选择或自动计算得到焦点区域。人工选择的一种简单直接的方式是：指定2个同心球 $S (c, r)$ 和 $S (c, R)$ ， $r < R$ ，以及参数 $λ \geq 1$ 。首先，计算目标测度函数 $μ : Ω \to R$ ，较大的球体 $S (c, R)$ 外部等于1，较小的球体 $S (c, r)$ 内部等于 $λ$ 。此外，μ是两个球体之间区域的调和函数，满足具有狄利克雷边界条件的拉普拉斯方程：

\{\begin{array}{l} Δ μ (p) = 0, & r < |p - c| < R ， \\ μ (p) = 1, & |p - c| = r ， \\ μ (p) = λ, & |p - c| = R 。 \end{array}

（15）

然后，对目标测度归一化，计算总体积：

V = \int_{Ω} μ （ p ） d p = \sum_{t_{i j k l}} \int_{t_{i j k l}} μ (p) d p 。

其中， $t_{i j k l}$ 表示由 $v_{i}$ ， $v_{j}$ ， $v_{k}$ 和 $v_{l}$ 组成的四面体。最后，将目标测度按4/V缩放， $μ \leftarrow 4 μ / V$ 。

4.2　离散化

本文算法将目标域离散为一组点 $Y = \{q_{1}, q_{2}, \dots, q_{k}\}$ 。用Delaunay细化算法^［2］对Ω进行三角剖分，通过指定每个四面体的最大体积，获得均匀的采样；用TetGen^［27］计算Ω的四面体网格。四面体网格仍记为 $Ω = (V, E, F, T)$ ，其中V，E，F，T分别表示顶点、边、面和四面体的集合。

用户可以选择聚焦区域，选择两个近似于球体的三角形网格。内部和外部的球体仍表示为 $S (c, r)$ 和 $S (c, R)$ 。离散Laplace-Beltrami算子采用有限元方法进行建模^［28］。设 $e_{i j} \in E$ 是一条边，则余切边权重定义为

w_{i j} = \frac{1}{12} \sum_{k l} l_{k l} c o t θ_{i j}^{k l},

（16）

其中， $l_{k l}$ 为边 $e_{k l}$ 的长度， $θ_{i j}^{k l}$ 为四面体 $t_{i j k l}$ 中边 $e_{k l}$ 上的二面角，对所有与边 $e_{i j}$ 相邻的四面体求和。

目标测度μ近似于一个分段线性函数，表示为定义在顶点集 $μ : V \to R$ 上的函数，并且该测度对内部点进行线性插值，即假设点p在四面体 $t_{i j k l}$ 内，则

μ (p) = λ_{i} μ (v_{i}) + λ_{j} μ (v_{j}) + λ_{k} μ (v_{k}) + λ_{l} μ (v_{l}),

（17）

其中， $（ λ_{i}, λ_{j}, λ_{k}, λ_{l} ）$ 分别为点p关于 $v_{i}, v_{j}, v_{k}, v_{l}$ 的质心坐标，是下列线性方程的唯一解：

\{\begin{array}{l} λ_{i} v_{i} + λ_{j} v_{j} + λ_{k} v_{k} + λ_{l} v_{l} = p ， \\ λ_{i} + λ_{j} + λ_{k} + λ_{l} = 1 。 \end{array}

所以，离散Laplace-Beltrami算子为

Δ μ (v_{i}) = \sum_{e_{i j} \in E} w_{i j} (μ (v_{j}) - μ (v_{i})) 。

将式（15）离散为线性方程组。可以证明，离散Laplace-Beltrami算子是线性空间 $\sum_{i} x_{i} = 0$ 上的正定矩阵，因此，离散拉普拉斯方程有唯一的常数解。在获得初始测度 $μ : V \to R$ 后，首先计算总测度：

V = \int_{Ω} μ (p) d p = \sum_{t_{i j k l}} \int_{t_{i j k l}} μ （ p ） d p 。

因为μ是分段线性的，故每个四面体上μ的积分为

\begin{array}{l} \int_{t_{i j k l}} μ (p) d p = \frac{1}{4} [μ (v_{i}) + μ (v_{j}) + μ (v_{k}) + \\ μ (v_{l})] V o l (t_{i j k l}) 。 \end{array}

（18）

然后，通过 $μ \leftarrow 1 / V μ$ 对目标曲率归一化。以顶点为中心计算Ω的Voronoi单元分解：

Ω = \cup_{i = 1}^{k} W_{i}, W_{i} = \{p \in Ω | |p - v_{i}| \leq |p - v_{j}|, j = 1,2, \dots, k\} 。

最后，计算狄拉克测度：

μ_{j} = \int_{W_{j}} μ （ p ） d p 。

先计算Voronoi单元 $W_{i}$ 与Delaunay三角剖分中每个四面体 $t_{i j k l}$ 的交点，该交点为凸多面体；再将凸多面体分解为一对四面体，用式（18）计算每个四面体上μ的积分，得到样本的狄拉克测度 $\{(v_{1}, μ_{1}), (v_{2}, μ_{2} ）, \dots ，（ v_{k}, μ_{k} ）\}$ 。

4.3　离散OMT映射算法

首先，计算Ω的一个幂Voronoi单元划分，将单元 $W_{i} (h)$ 映射至样本点 $v_{i}$ 。然后，将此表示转换为从四面体网格Ω到自身的分段线性映射。将每个Voronoi单元 $W_{i} (h)$ ，分解为一对四面体 $t_{i}^{j}$ ，四面体的质心是其4个顶点的均值，记为 $c_{i}^{j}$ ， $W_{i} (h)$ 的质心为

c_{i} (h) = \frac{\sum_{j} c_{i}^{j} V o l (t_{i}^{j})}{\sum_{j} V o l (t_{i}^{j})} 。

如果 $W_{i} (h)$ 与Ω的边界相交，则将其质心投影至边界曲面 $\partial Ω$ ，并用投影图像替换质心 $c_{i} (h)$ 。

立方体Ω通过幂Delaunay三角剖分 $𝒯 (h)$ 进行三角剖分，顶点位置在源域中为 $c_{i} (h)$ ，在目标域中为 $v_{i}$ 。通过分段线性映射 $φ : （ Ω, 𝒯 (h) ） \to （ Ω, 𝒯 (h) ）$ ，将顶点从 $c_{i} (h)$ 映射至 $v_{i}$ ，有 $φ （ c_{i} (h) ） = v_{i}$ 。

OMT映射算法：

输入　一个凸域 $Ω \subset R^{3}$ 和一组离散点 $Y = \{q_{1}, q_{2}, \dots, q_{k}\}$ ，离散目标测度 $ν = \{ν_{1}, ν_{2}, \dots, ν_{k}\}$ ，使得 $\sum_{i} ν_{i} = V o l (Ω)$ 。

输出　 $Ω$ 的一个划分， $Ω = \underset{i}{\cup} W_{i}$ ，使得 $W_{i} \mapsto q_{i}$ 是OMT映射。

while true do

for $i \leftarrow 1$ tokdo

建立超平面 $π_{i} (h) : 〈q_{i}, p〉 + h_{i}$ 。

计算超平面的对偶点 $π_{i}^{*} (h)$ 。

建立对偶点集 $\{π_{i}^{*} (h)\}$ 的凸包 $𝒞 (h)$ 。

计算凸包的对偶，得到超平面集 $\{π_{i} (h)\}$ 的上包络 $ℰ (h)$ 。

投影 $𝒞 (h)$ 获得Y的幂Delaunay三角剖分 $𝒯 (h)$ 。

投影 $ℰ (h)$ 获得 $Ω$ 的幂Voronoi单元分解 $𝒱 (h)$ 。

for $i \leftarrow 1$ tokdo

计算 $W_{i} (h)$ 的体积 $w_{i} (h)$ 。

建立梯度 $\nabla E (h) = （ ν_{i} - w_{i} (h) ）^{T}$ 。

for 每条边 $e_{i j} \in 𝒯 (h)$ do

计算对偶面 $f_{i j}$ 的面积 $A_{i j}$ 。

建立Hessian矩阵 $\frac{\partial^{2} E (h)}{\partial h_{i} \partial h_{j}} = \frac{A_{i j}}{|q_{j} - q_{i}|}$ 。

求解线性方程 $H e s s (h) δ h = \nabla E (h)$ 。

λ \leftarrow 1

。

计算 $Ω$ 的幂Voronoi图 $𝒱 (h + λ δ h)$ 。

while 存在 $w_{i} (h + λ δ h)$ 为空 do

λ \leftarrow 1 / 2 λ

。

计算 $Ω$ 的幂Voronoi图 $𝒱 (h + λ δ h)$ 。

h \leftarrow h + λ δ h

。

if $\forall |w_{i} (h) - ν_{i}| < ε$ then break。

return 映射 $\{W_{i} \mapsto q_{i}, 1,2, \dots, k\}$ 。

4.4　三线性体素重采样

假设给定一个体图像数据S，使用最优质量传输映射 $φ : Ω \to Ω$ 生成另一个具有不同分辨率和体积测度的体图像数据T。将体图像数据表示为三维数组，对于T中的每个体素 $T (i, j, k)$ ，在Ω中定位其位置p，并找到包含p的四面体 $t_{p q r s}$ ，用式（17）计算重心坐标 $(λ_{p}, λ_{q}, λ_{r}, λ_{s})$ ，p的原像可表示为

φ^{- 1} (p) = λ_{p} c_{p} (h) + λ_{q} c_{q} (h) + λ_{r} c_{r} (h) + λ_{s} c_{s} (h) 。

然后，找到S中与 $φ^{- 1} (p)$ 相邻的体素，用三线性插值得到 $φ^{- 1} (p)$ 的强度，并将其赋值于 $T (i, j, k)$ 。

5　实验结果

实验设备为一台具有双处理器2.9 GHz CPU和8 GB内存的笔记本电脑。在Windows平台采用通用C++开发系统。基于包含计算幂Voronoi图和Delaunay三角剖分模块的CGAL库实现OMT。在GPU上使用基于raycasting算法的Voreen完成体绘制。

5.1　大变形和全局光滑

体动脉瘤数据的分辨率为256×256×256。焦点区域为实心球，中心c=（0.06，0.32，0.04），半径r=0.1。放大因子为1，27，125，343，例如放大因子为343表示将感兴趣区域（region of interest，ROI）的体积放大343倍，如图2所示，可知，OMT能够设计大变形并保持映射的平滑性。因为OMT映射是凸函数的梯度映射，所以OMT映射的雅可比矩阵就是函数的Hessian矩阵。式（3）中的Monge-Amer'e方程表明雅可比矩阵处处为正，因此OMT映射是全局微分同胚的。即便放大因子很大，整个变形仍处处光滑，更进一步，映射没有折叠，因此空间和拓扑关系都得到了很好的保留。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 用OMT技术放大的动脉瘤

Fig.2 The volumetric aneurism magnified by large scales using the OMT technique

5.2　精确控制

图3展示了体元变形准确性的验证结果。盆栽模型的分辨率为512×512×512，强度等级为255。焦点区域是一个实心球，圆心c=（0，0，0），半径为0.2，如图3（b）中红色实心球所示。采用OMT映射算法，在放大因子为27的条件下进行体积变形，焦点区域如图3（e）中的实心球所示，对于每个幂Voronoi单元 $W_{i} (h)$ ，计算其最终值v与期望测度 $v_{i}$ 的比值，取对数，并作如图3（g）所示的直方图。可知，直方图高度集中在原点，证明本文算法能在体素水平上精确控制变形。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 对盆栽模型测度的精确控制

Fig. 3 Measure accurate control for Bonsai model

5.3　变形的光滑性

体积变形由目标测度的连续性控制，如图4所示。在非调和方式下，焦点区域是以c=（-0.42，0.13，0.26）为中心，半径r=0.26的实心球。目标测度μ在球内为64，在球外为1，因此μ是不连续的。如图4（b）~（d）所示，这种测度会引起较大的形变，特别是焦点区域的边界处。在调和方式下，选择2个同心球，中心为c=（-0.42，0.13，0.26），半径分别为r=0.15和r=0.26。目标测度μ在球内为64，在球外为1，是两个球体之间的调和函数。因此，μ是一个光滑函数，其引起的形变明显较小，如图4（e）~（g）所示。表明体的全局变形受目标测度的平滑性影响较大。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 目标测度μ的平滑程度对脚模型变形平滑程度的影响

Fig. 4 The effectiveness of the smoothness of the target measure μ on the smoothness of the deformation

5.4　不规则感兴趣区域

图5为NCAT phantom模型的F+C可视化结果，可知，本文算法允许用户定义高度不规则的感兴趣区域。焦点区域不是规则的实心球，而是包裹了脊柱的凹形不规则形状，严格拟合了体数据中的感兴趣内容。OMT映射仍是微分同胚的，且可以完全控制目标测度。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 NCAT phantom模型的F+C可视化结果

Fig. 5 F+C visualization for NCAT phantom model

5.5　多焦点区域

采用基于OMT的F+C方法，选择多个焦点区域，可为不同的兴趣区域分配不同的放大因子。在图6所示的CT膝关节模型中，选择两个兴趣区域：一个是左髌骨区域的球，中心 $c_{1}$ =（-0.4，-0.36，-0.50），半径为0.19，放大因子 $λ_{1}$ =8的球；另一个是右胫骨的实心球，中心 $c_{2}$ =（0.46，-0.06，-0.71），半径为0.15，放大因子 $λ_{2}$ =16。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 多焦点多视角放大的CT膝关节模型的F+C可视化结果

Fig. 6 F+C visualization for CT knee model with multiple focus magnification in different views

5.6　时间复杂度

表1列出了本文方法在各模型上的运行时间。Ω是边长为2的标准立方体，离散样本数约为10 000。手动选择感兴趣的区域，OMT映射的计算是离线执行的。

表1 运行时间

Table 1 Running time

模型	数据来源	运行时间/s	分辨率
动脉瘤	Philips Research，Hamburg，Germany	118	512×512×512
盆栽	Rosttger S，VIS，University of Stuttgart	156	512×512×512
脚	Philips Research，Hamburg，Germany	182	256×256×256
NCAT phantom	Segars WP，Tsui BMW	156	512×512×512
CT膝关节	Department of Radiology University of Iowa	134	440×440×440

新窗口打开| 下载CSV

6　优点和缺点

提出了一种基于最优质量传输的F+C可视化方法。与现有方法相比，本文方法：（1）具有坚实的理论基础，保证了解的存在性、唯一性和正则性；（2）允许用户精确控制目标体积元（测度）；（3）保证全局微分同胚，变形是全局平滑的；（4）为形状不规则的多个焦点区域提供了灵活性；（5）可推广至任意维度。

另外，该方法基于非线性凸优化，需幂 Delaunay三角剖分和幂Voronoi图的复杂几何数据结构，且组合结构是动态变化的。由于优化主要在CPU上进行，故无法实时计算。但Voronoi图可以仅使用GPU进行计算^［29］，意味着未来有可能实现实时计算。

7　结语

提出了一种基于最优质量传输理论的F+C可视化方法，保证了解的存在性、唯一性和平滑性。允许用户精确控制目标体积元，可选择多个形状不规则的焦点区域，且可通过计算几何中的幂Voronoi图和Delaunay三角剖分实现优化。此外，还可将该方法推广至更高维度。

未来，将在GPU上实现最优质量变换映射算法，并将可视化技术推广至更高维度。

http://dx.doi.org/10.3785/j.issn.1008-9497.2023.06.003

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

WANG

Y S

， WANG

， LEE

T Y

， et al.

Feature-preserving volume data reduction and Focus+ Context visualization

［J］. IEEE Transactions on Visualization and Computer Graphics， 2010， 17（2）： 171-181. DOI：10.1109/tvcg.2010.34

[本文引用: 3]

[2]

， LUO

， SUN

， et al.

Variational principles for Minkowski type problems， discrete optimal transport， and discrete Monge-Ampere equations

［J］. Asian Journal of Mathematics， 2016， 20（2）： 383-398. doi:10.4310/ajm.2016.v20.n2.a7

[本文引用: 3]

[3]

SARKAR

， BROWN

M H

Graphical fisheye views of graphs

［C］// SIGCHI Conference on Human Factors in Computing Systems. Monterey： SIGCHI， 1992： 83-91. DOI：10.1145/142750.142763