Differentiable Voronoi Diagrams for Simulation of Cell-Based Mechanical Systems

一句话总结

用可微 Voronoi（幂）图把每个细胞压缩成一个 Voronoi 站点（仅 3-5 个自由度），隐式地定义界面网络的形状与拓扑，从而在细胞分裂、合并、换邻居等剧烈拓扑变化中保持连续、可用牛顿法求解，并支持与刚性/可变形边界耦合以及基于灵敏度分析的逆问题优化。

研究背景

细胞类力学系统（生物组织、肥皂泡沫）的核心特征是：细胞在运动、生长分裂或合并塌缩过程中不断改变邻接关系。这些拓扑跃迁会引起力的突变，给数值模拟带来巨大挑战。

现有方法各有短板：

• 均质化 / 抽象模型：把细胞结构平均掉，只在粗尺度上刻画力学，无法预测由局部拓扑变化主导的细胞级交互。
• 显式网络（vertex model）：显式建模细胞界面网络，能高效计算体积等几何量；但拓扑变化需要离散的网格重构和其他不连续操作，且要额外约束防止自交，对可微模拟极不友好。
• 可变形细胞模型（deformable cell model, DCM）：每个细胞用独立的表面网格表示，能精确建模接触力和拓扑变化，但每个细胞需要成百上千的自由度，且必须在每个界面显式处理碰撞，规模一旦达到几百上千个细胞就几乎不可解。

作者提出以”细胞为中心”的思路：用广义 Voronoi 图，每个细胞只用一个站点表示，隐式地定义细胞的形状与拓扑。这样自由度骤减、无需显式接触处理、并且在拓扑跃迁时保证几何连续变化。

核心方法

为什么用幂图（Power Diagram）

标准 Voronoi 图用欧氏距离定义细胞：

\[ R_i = \{p : d(p, C_i) < d(p, C_j)\ \forall j \neq i\} \]

已有理论指出，二维泡沫在平衡态几何上等价于”分段乘性 Voronoi 划分”(SMVP)：

\[ d_{\text{SMVP}}(p, C_i)^2 = \frac{1}{k_i}\left(d_{\text{Euclidean}}(p, C_i)^2 - w_i\right) \]

其中乘性权重 \(k\) 会产生弯曲界面。弯曲界面虽然能精确表达真实泡沫，但会让几何量的解析积分（3D 曲面）几乎不可能，边界裁剪也要解高次多项式。因此作者退而选择幂图：

\[ d_{\text{Power}}(p, C_i)^2 = d_{\text{Euclidean}}(p, C_i)^2 - w_i \]

它保持平面界面，同时通过权重 \(w_i\) 允许细胞大小可变，在计算可行性与表达力之间取得平衡。

闭式可微的 Voronoi 顶点

作者计算限制在包围域内的受限 Voronoi 图。在 \(m\) 维中，每个受限 Voronoi 顶点是 \(m\) 个超平面（Voronoi 平分面或边界面）的交点。3D 中存在四类顶点：无约束顶点、边界面顶点、边界棱顶点、简单边界顶点。

以无约束顶点为例，给定 Delaunay 单纯形的四个站点，平分面方程为线性方程，三个平分面的交点即求解一个 \(3\times3\) 线性系统：

\[ \begin{bmatrix} x_1-x_0 & y_1-y_0 & z_1-z_0 \\ x_2-x_0 & y_2-y_0 & z_2-z_0 \\ x_3-x_0 & y_3-y_0 & z_3-z_0 \end{bmatrix} \begin{bmatrix} x_c \\ y_c \\ z_c \end{bmatrix} = \frac{1}{2}\begin{bmatrix} x_1^2-x_0^2+y_1^2-y_0^2+z_1^2-z_0^2 \\ \cdots \\ \cdots \end{bmatrix} \]

由于 \(3\times3\) 矩阵求逆有已知闭式表达式，每个顶点坐标可写成输入自由度（站点位置和边界顶点）的闭式函数，从而解析可微。切换到幂图只需在右端加 \(\tfrac{1}{2}(w_i - w_j)\)；边界顶点只需把部分平分面方程替换成边界面方程。实现中用计算机代数系统自动生成求值和求导代码。

细胞积分与能量形式

细胞力学由各细胞的能量驱动，能量依赖于面积、体积等几何量。多面体 Voronoi 细胞的积分分片进行：把每个面三角化后逐三角求和。计算体积积分时无需显式四面体化——而是把每个表面三角形连到原点构成四面体（可正可负体积），这样单纯形积分只依赖三个顶点而非四个，减少了每个单纯形所需的一阶、二阶导数数量。

系统总能量写成

\[ E = F(\mathbf{c}, \mathbf{x}(\mathbf{c})) \]

其中 \(\mathbf{c}\) 是所有站点的自由度，\(\mathbf{x}\) 是所有 Voronoi 顶点坐标（由 \(\mathbf{c}\) 隐式决定）。梯度与 Hessian 通过链式法则展开：

\[ \frac{dE}{d\mathbf{c}} = \frac{\partial F}{\partial \mathbf{x}}\frac{d\mathbf{x}}{d\mathbf{c}} + \frac{\partial F}{\partial \mathbf{c}} \]

这个通用形式允许任意可微能量，同一套计算方法既能模拟泡沫也能模拟生物组织。

边界耦合

为处理外部物体与自由空间，用三角网格（顶点 \(\mathbf{v}(\mathbf{p})\)，\(\mathbf{p}\) 为边界自由度）定义包围域，令 \(\mathbf{y}=(\mathbf{c},\mathbf{p})\) 为耦合系统全部自由度，总能量为

\[ E = F(\mathbf{c}, \mathbf{x}(\mathbf{c}, \mathbf{v}(\mathbf{p}))) + F_B(\mathbf{p}) \]

其中 \(F_B\) 是边界能量（如包裹膜的弹性势）。这让模型能同时容纳弹性膜、自由表面、刚体和固定边界。

光滑性

一般构型下每个 \(m\) 维 Voronoi 顶点度为 \(m+1\)，拓扑在邻域内不变，顶点坐标是站点的局部光滑函数。拓扑跃迁时出现更高度数顶点，\(d\mathbf{x}/d\mathbf{y}\) 无定义，能量非光滑。好在细胞体积（及光滑函数的体积积分）在跃迁态跨越时是 \(C^1\) 连续的，而表面积（2D 周长）只有 \(C^0\)。若细胞能量对表面积单调且边界域为凸，梯度不连续处形成严格凹的局部极大，不妨碍最小化；作者在非凸边界的实验中均使用 \(C^1\) 连续能量以规避收敛问题。

动力学与逆问题

从静态平衡（\(F(\mathbf{y})=-dE/d\mathbf{y}=0\)）推广到牛顿第二定律并引入黏性力，运动方程

\[ M\ddot{\mathbf{y}} + \boldsymbol{\eta}\dot{\mathbf{y}} + \frac{dE}{d\mathbf{y}} = 0 \]

用有限差分离散（\(C^1\) 能量下可用 BDF2 二阶格式获得二次收敛），每一步转化为优化问题求 \(\mathbf{y}_{k+1}\)。许多细胞系统处于低雷诺数区，惯性被黏性主导，可略去动量项得到无动量（quasi-static）模型。

对于逆问题，作者用平衡约束优化框架：

\[ \min_{\mathbf{y}, \mathbf{u}} L(\mathbf{y}, \mathbf{u}) \quad \text{s.t.}\ \frac{\partial}{\partial \mathbf{y}}E(\mathbf{y}, \mathbf{u}) = 0 \]

借助隐函数定理计算灵敏度 \(d\mathbf{y}/d\mathbf{u} = -(\partial^2 E/\partial\mathbf{y}^2)^{-1}(\partial^2 E/\partial\mathbf{y}\partial\mathbf{u})\)，即可用一阶或拟牛顿方法优化。

技术细节

• 实现：C++ + Eigen；线性系统用 CHOLMOD；Voronoi 图生成借助 CGAL 与 Geogram；可视化用 Polyscope。代码开源。
• 泡沫粗化建模：每个站点五个自由度——三维坐标 \((x,y,z)\)、幂图权重 \(w\) 和体积目标 \(\bar{V}\)（代表气体质量）。通过让权重变化，细胞可以平滑收缩到零体积并自动完成局部拓扑重构，而无需 vertex model 那种删除塌缩细胞的启发式规则。
• 胚胎卵裂建模：从单个 Voronoi 细胞出发，边界为可变形膜（边长 \(\ell\) 的边能量 \(E_B=k\ell^2\)）。每 \(k\) 步把每个站点替换为一对子站点 \(c\pm\beta\bar{V}^{1/3}\mathbf{n}\)，前三次卵裂正交，之后随机。子站点继承父站点权重、静止体积取一半。

实验结果

• 细胞内导航：刚体在类组织的弹性 Voronoi 细胞装配体中被恒力推动，出现类湍流的不规则位移；非对称形状刚体因周围细胞的力产生弯曲轨迹，验证细胞与边界自由度的动态耦合。
• 泡沫粗化：从 2000 个单分散干泡沫细胞出发，350 帧内塌缩到 2 个细胞的平衡态，隐式完成完整拓扑重构。
• 组织生长：胚胎卵裂在球形膜内模拟到 4096 个细胞，与蝾螈胚胎发育实拍图定性相似；柱形容器内的细胞增殖（带重力、自由上表面）出现频繁不规则拓扑变化。
• 从图像表征泡沫：对实验室拍摄的 2D 肥皂泡沫图像做标注，用平衡约束优化调整各细胞面积目标 \(\bar{A}\)（即压强），目标函数为标注顶点与对应 Voronoi 顶点的平方距离和。经 77 次 L-BFGS 迭代、13 秒收敛，目标值从 0.28 降到 0.034。
• 与 DCM 对比：DCM 每细胞 30 顶点（60 自由度），本方法每细胞仅 3 自由度。30 细胞准静态模拟中两者达到相同最终拓扑，但 DCM 耗时约 1000 倍、平均牛顿迭代次数近 40 倍（见运行时统计）。
• 规模扩展：随机放置 \(n\) 个站点做单次牛顿迭代，Voronoi 图生成与 Hessian 计算大致线性扩展，大规模时线性系统求解主导开销。3D 明显比 2D 昂贵，主因是平均邻居数从 2D 的 6 增至 3D 的约 15.5，导致 Hessian 更稠密、顶点更多。

运行时统计（AMD Ryzen Threadripper PRO 5995WX）：

实验	细胞数	自由度	帧数	迭代/帧 (均/最大)	每迭代耗时 (ms)	总耗时 (s)
Comparison (Ours)	30	90	100	5/25	2.32	1.23
Comparison (DCM)	30	1800	100	195/2621	65	1258
Tissue Growth (membrane)	4096	24070	700	3/12	3830	8521
Tissue Growth (cylinder)	620	6067	1101	10/47	1383	14970
Coarsening	2000	10000	350	51/522	1275	22996
Rigid Body 1	2000	4003	396	3/4	212	251
Rigid Body 2	2000	4003	304	3/5	213	194

贡献与局限

贡献：

• 提出基于可微 Voronoi（幂）图的细胞系统模拟新范式，用极其紧凑的状态（每细胞一个站点）隐式定义界面网络的形状与拓扑，天然处理分裂/合并/换邻居等拓扑跃迁且保持几何连续。
• 推导 Voronoi 顶点关于站点位置的一阶、二阶闭式导数，支持牛顿型求解器和广泛的每细胞能量函数；并统一处理细胞面与裁剪几何，实现与移动/可变形边界的耦合。
• 通过灵敏度分析获得平衡态导数，打通逆模拟问题（如从真实泡沫图像反求细胞压强）。

局限与未来工作：

• 采用各向同性距离度量且假设平面/分段线性界面，限制了可表达的细胞形状；真实上皮组织存在非凸的 scutoid 形细胞，无法表达。
• 缺乏弯曲界面导致在细胞尺寸和压强差异大的泡沫模拟中几何精度受限。
• 生物模拟要真正有信息量还需更复杂的能量模型与真实数据标定参数；推导 Voronoi 表面积的更光滑近似可支持细胞黏附等效应。
• 未来方向包括更复杂细胞形状的紧凑表示、泡沫类 3D 打印超材料的梯度优化，以及引导物体穿越细胞系统的控制算法设计。

延伸思考

这篇工作的精髓是”表示的降维”：把原本需要显式网格和碰撞处理的细胞界面，压缩成由少数站点隐式生成的 Voronoi 结构。拓扑变化不再是需要特判的离散事件，而是站点连续移动的自然结果——这与网格重构、level set 等”事后处理拓扑”的思路形成鲜明对比。用幂图（平面界面）而非乘性 Voronoi（曲面）也是典型的可微性/可解性优先于物理精确性的工程取舍，代价是泡沫几何精度受限。这种”闭式几何 + 隐式拓扑 + 灵敏度分析”的组合，对元材料设计、组织形态发生的逆向工程等需要在拓扑可变空间里做梯度优化的问题，提供了一个很有启发性的模板。