Author: nhu_tran

The Log Diffusion Solver

Post author By nhu_tran
Post date October 20, 2025

Teaser Figure — **Figure 1:** Comparison of distance estimation by heat methods. Standard diffuse-then-log (middle) suffers from numerical issues when the heat is small (large geodesic distances), leading to noisy and inaccurate isolines. Our proposed log-space diffusion (right) performs diffusion directly in the log-heat domain, producing smooth and accurate estimates, closely matching the ground-truth distances (left). It is worth noting that our solver outperforms others primarily for smaller time step values.

Introduction

The heat equation is one of the canonical diffusion PDEs: it models how “heat” (or any diffusive quantity) spreads out over space as time progresses. Mathematically, let $u(x,t) \in \mathbb{R}_{\geq 0}$ denote the temperature at spatial location $x \in \mathbb{R}^{d}$ and time $t \in \mathbb{R}$. The continuous heat equation is given by

\[\frac{\partial u}{\partial t} = \Delta u,\]

where $\Delta u (x, t) = \sum_{i=1}^{d} \frac{\partial^2 u}{\partial x_i^2}$ is the Laplace operator representing how local temperature gradients drive heat from hotter to colder regions. On a surface mesh, the spatial domain is typically discretized so that the PDE reduces to a system of ODEs in time, which can be integrated using stable and efficient schemes such as backward Euler. Interestingly, in many applications—including geodesic distance computation ¹ and optimal transport ²—it is often more useful to work with the logarithm of the heat rather than the heat itself. For example, Varadhan’s formula states that the geodesic distance $\phi$ between two points $x$ and $x^*$ on a curved domain can be recovered from the short-time asymptotics of the heat kernel:

\[\phi(x, x^{*}) = \lim_{t \to 0} \sqrt{-4t \,\log k_{t,x^{*}}(x)},\]

where $k_{t,x^{*}}(x) = u(x,t)$ is the heat distribution at time $t$ resulting from an initial Dirac delta at $x^*$, i.e., $u(x,0) = \delta(x-x^{*})$. In practice, this leads to a two-step workflow: (1) simulate the heat equation to obtain $u_t$, and (2) take the logarithm to recover the desired quantity. However, this approach can be numerically unstable in regions where $u$ is very small: tiny relative errors in $u$ can be greatly magnified in the log domain. For example, an absolute error of $10^{-11}$ in $u$ becomes an error of roughly $25.33$ in $\log u$.

A natural idea is therefore to evolve the logarithm of the heat directly. Let $u>0$ satisfy the heat equation and define $w = \log u$. One can show that $w$ obeys the nonlinear PDE ³

\[\frac{\partial w}{\partial t} = \Delta w + \|\nabla w\|_2^2.\]

This formulation has the appealing property that it avoids taking the logarithm after simulation, thereby reducing the amplification of small relative errors in $u$ when $u$ is tiny. However, the price of this reformulation is the additional nonlinear term $\|\nabla w\|_2^2$, which complicates time integration as each step requires iterative optimization, and introduces a discrepancy with the discretized heat equation. In practice, evaluating $\|\nabla w\|_2^2$ requires a discrete gradient operator, and the resulting ODE is generally not equivalent to the one obtained by first discretizing the heat equation for $u$ (which involves only a discrete Laplacian) up to a log transform. This non-equivalence stems from (1) there is a new spatially discretized gradient operator, and (2) the fact that the logarithm is nonlinear and does not commute with discrete differential operators.

Motivated by this discrepancy, we take a different approach that preserves the exact relationship to the standard discretization. We first discretize the heat equation in space, yielding the familiar linear ODE system for $u$. We then apply the transformation $w = \log u$ to this discrete system, producing a new ODE for $w$. Now, by construction, the new ODE is equivalent to the heat ODE up to a log transform of variable. This still allows us to integrate directly in log-space, which allows more accurate and stable computation. We show that this formulation offers improved numerical robustness in regions where $u$ is very small, and demonstrate its effectiveness in applications such as geodesic distance computation.

**Table 1:** Side-by-side procedures for computing $\log u$ in heat diffusion, note the highlighted steps are the same for standard heat diffusion and ours.
Standard (post-log)	Continuous log-PDE	Proposed discrete-log
1. Start with the heat PDE for $u$.	1. Start with the log-heat PDE for $w=\log u$.	1. Start with the heat PDE for $u$.
2. Discretize in space $\Rightarrow$ ODE for $u$.	2. Discretize in space for $\nabla$ and $\Delta$ $\Rightarrow$ ODE for $w$.	2. Discretize in space $\Rightarrow$ ODE for $u$.
3. Time-integrate in $u$-space.	3. Time-integrate in $w$-space and output $w$.	3. Apply change of variables $w=\log u$ $\Rightarrow$ ODE for $w$.
4. Compute $w=\log u$ as a post-process.		4. Time-integrate in $w$-space and output $w$.

Background

We begin by reviewing how previous approaches are carried out, focusing on the original heat method ¹ and the log-heat PDE studied in ³. In particular, we discuss two key stages of the discretization process: (1) spatial discretization, where continuous differential operators such as the Laplacian and gradient are replaced by their discrete counterparts, turning the PDE into an ODE; and (2) time discretization, where a numerical time integrator is applied to solve the resulting ODE. Along the way, we introduce notation that will be used throughout.

Heat Method

Spatial discretization. In the heat method ¹, the heat equation is discretized on triangle mesh surfaces by representing the heat values at mesh vertices as a vector $u \in \mathbb{R}^{|V|}$, where $|V|$ is the number of vertices. A corresponding discrete Laplace operator is then required to evolve $u$ in time. Many Laplacian discretizations are possible, each with its own trade-offs ⁴. The widely used cotangent Laplacian is written in matrix form as $M^{-1}L$, where $M, L \in \mathbb{R}^{|V| \times |V|}$. The matrix $L$ is defined via

\[L_{ij} = \begin{cases} \frac{1}{2}(\cot \alpha_{ij} + \cot \beta_{ij}) \quad & i \neq j \\ -\sum_{i \neq j}L_{ij} & i = j \end{cases},\]

where $i,j$ are vertex indices and $\alpha_{ij}, \beta_{ij}$ are the two angles opposite the edge $(i,j)$ in the mesh. The mass matrix $M = \mathrm{diag}(m_1, \dots, m_{|V|})$ is a lumped (diagonal) area matrix, with $m_i$ denoting the area associated with vertex $i$. Note that $L$ has several useful properties: it is sparse, symmetric, and each row sums to zero. This spatial discretization converts the continuous PDE into an ODE for the vector-valued heat function $u(t)$:

\[\frac{\mathrm{d}u}{\mathrm{d}t} = M^{-1}L\,u.\]

Time discretization and integration. The ODE is typically advanced in time using a fully implicit backward Euler scheme, which is unconditionally stable. For a time step $\Delta t = t_{n+1} – t_n$, backward Euler yields:

\[(M – \Delta t\,L)\,u_{n+1} = M\,u_n,\]

where $u_n$ is given and we solve for $u_{n+1}$. Thus, each time stepping requires solving a sparse linear system with system matrix $M – \Delta t\,L$. Since this matrix remains constant throughout the simulation, it can be factored once (e.g., by Cholesky decomposition) and reused at every iteration, greatly improving efficiency.

Log Heat PDE

Spatial discretization. In prior work ³ on the log-heat PDE, solving this requires discretizing not only the Laplacian term but also the additional nonlinear term $\|\nabla w\|_2^{2}$. After spatial discretization, the squared gradient norm of a vertex-based vector $w \in \mathbb{R}^{|V|}$ is computed as

\[g(w) = W \left( \sum_{k=1}^{3} (G_k w) \odot (G_k w) \right) \in \mathbb{R}^{|V|}\]

where $G_{k} \in \mathbb{R}^{|F| \times |V|}$ is the discrete gradient operator on mesh faces for the $k$-th coordinate direction, and $W \in \mathbb{R}^{|V| \times |F|}$ is an averaging operator that maps face values to vertex values. With this discretization, the ODE for $w(t) \in \mathbb{R}^{|V|}$ becomes

\[\frac{\mathrm{d}w}{\mathrm{d}t} = M^{-1}L\,w + g(w)\]

Time discretization and integration. If we apply fully implicit Euler, we would find $w^{n+1}$ by solving a nonlinear equation. To avoid this large nonlinear solve, prior work ³ splits the problem into two substeps:

a linear diffusion part: $\frac{\mathrm{d}w}{\mathrm{d}t} = M^{-1}L\,w$, which is stiff but easy to solve implicitly,
a nonlinear part: the full equation, which is nonlinear but convex.

This separation lets each component be treated with the method best suited for it: a direct linear solve for the diffusion, and a convex optimization problem for the nonlinear term. Using Strang-Marchuk splitting:

Half-step diffusion: $(M – \frac{\Delta t}{2} L)\,w^{(1)} = Mw^n$
Full-step nonlinear term: Solve a convex optimization problem to find $w^{(2)}$
Half-step diffusion: $(M – \frac{\Delta t}{2} L)\,w^{(3)} = Mw^{(2)}$, and set $w^{n+1} = w^{(3)}$

A New ODE in Log-Domain

Let us revisit the spatially discretized heat (or diffusion) equation. For each vertex $i \in \{1, \dots |V|\}$ on a mesh, its evolution is described by

\[\frac{\mathrm{d} u_i}{\mathrm{d} t} = \frac{1}{m_i} (L u)_i\]

where $u$ represents our heat, and $L$ is the discrete Laplacian operator. A key insight is that the Laplacian term can be rewritten as a sum of local differences between a point and its neighbors:

\[(Lu)_i = \sum_{j \sim i}L_{ij} (u_j – u_i).\]

This shows that the change in $u_i$ depends on how different it is from its neighbors, weighted by the Laplacian matrix entries. Now, before doing time discretization, we use a change of variable $w \Leftarrow\log u$, accordingly, $u = \mathrm{e}^{w}$, to obtain the ODE

\[\mathrm{e}^{w_i}\frac{\mathrm{d} w_i}{\mathrm{d}t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w_j} – \mathrm{e}^{w_i}),\]

which can be rearranged into

\[\frac{\mathrm{d} w_i}{\mathrm{d}t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w_j \, – \, w_i} – 1).\]

In vectorized notation, we construct a matrix $\tilde{L} \Leftarrow L – \mathrm{diag}(L_{11}, \dots, L_{|V| \,|V|})$, which zeros out the diagonals of $L$. Then the log heat $w(t) \in \mathbb{R}^{|V|}$ follows the ODE

\[\frac{\mathrm{d} w}{\mathrm{d}t} = M^{-1} \left( \mathrm{diag}(\mathrm{e}^{-w})\, \tilde{L}\,\mathrm{e}^{w} – \tilde{L} 1 \right),\]

and another equivalent way to put it is

\[\frac{\mathrm{d} w}{\mathrm{d}t} = M^{-1} \left( \left(\tilde{L} \odot \exp(1 w^{\top} – w 1^{\top} ) \right) 1 – \tilde{L} 1 \right).\]

This ODE is equivalent to the standard heat ODE up to a log transform of the integrated results, but crucially it avoids the amplification of small relative errors in $u$ when $u$ is tiny by integrating directly in log-space.

Numerical Solvers for Our New ODE

Forward (Explicit Euler)

The most straightforward approach uses explicit Euler to update the solution at each time step. The update rule is

\[w_{n+1} = w_{n} + \Delta t\, M^{-1} \left( \mathrm{diag}(\mathrm{e}^{-w_n})\, \tilde{L}\,\mathrm{e}^{w_n} – \tilde{L} 1 \right)\]

While this is efficient and can be solved linearly, it causes instability.

Semi-Implicit

We can discretize with both explicit and implicit parts by decomposing the exponential function as a linear and nonlinear part:

\[\mathrm{e}^{x} = \underbrace{(1 + x)}_{\text{implicit}} + \underbrace{(\mathrm{e}^{x} – x – 1)}_{\text{explicit}}\]

Applying this decomposition yields a linear system for solving $w_{n+1}$:

\[\left( M + \Delta t \left(\tilde{L} – \mathrm{diag}(\tilde{L} 1) \right) \right) w_{n+1} = M w_{n} – \Delta t \left( \mathrm{diag}(\mathrm{e}^{-w_n}) \tilde{L} \mathrm{e}^{w_{n}} – \tilde{L} w_{n} + \mathrm{diag}(w_{n}) \tilde{L} 1 – \tilde{L} 1 \right)\]

This retains efficiency but does not fully address instability issues.

Backward (Convex Optimization)

Solving the backward Euler equation involves finding the vector $w^{t+1}$ that satisfies

\[\frac{w^{t+1}_i – w^{t}_i}{\Delta t} = \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w^{t+1}_j \, – \, w^{t+1}_i} – 1)\]

Taking inspiration from prior work ³, we find that relaxing this equality into an inequality leads to a convex constraint on $w^{t+1}_i$:

\[\frac{w^{t+1}_i – w^{t}_i}{\Delta t} \geq \frac{1}{m_i} \sum_{j \sim i}L_{ij} (\mathrm{e}^{w^{t+1}_j \, – \, w^{t+1}_i} – 1)\]

It can be shown that minimizing the objective function

\[f(w) = \sum_{i}M_iw_i\]

with the inequality constraint defined above yields the exact solution to the original equation. We have written a complete proof for this finding here. We implement this constrained optimization problem using the Clarabel solver ⁶ in CVXPY ⁷, a Python library for convex optimization problems.

Distance computations by different solvers — **Figure 2:** Visualization of distance computations on the mesh by different solvers. Ground Truth, HeatBackward, LogHeatBackward, NewLogHeatBackward (from left to right)

Experiments

As a common application, we demonstrate the effectiveness of log heat methods for computing fast approximations of geodesic distances. As ground truth, we used libigl’s ⁵ implementation of exact discrete geodesic distances. We used the mesh of a kitten with 10 diffusion time-steps equally spaced between $t=0$ and $t=10^{-3}$.

**Figure 3:** Correlation plot for log heat methods against an exact geodesic computation.

**Table 2:** Correlation analysis results.
Mesh	# of verts	Method	Pearson $r$ ↑	RMSE ↓	MAE ↓
Cat	9447	HeatBackward	0.986	0.285	0.241
		LogHeatBackward	-0.183	0.577	0.522
		NewLogHeatBackward (ours)	0.879	0.163	0.112

The flat region for the NewLogHeatBackward Solver probably indicates regions of numerical overflow/underflow, since the log of tiny heat values is going to be incredibly negative. Through our experimentation we observed that the solver for this our formulated PDE outperforms others mostly in cases where the diffusion time-step is extremely small.

References

[1] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. “The Heat Method for Distance Computation.” Communications of the ACM, vol. 60, no. 11, October 2017, pp. 90–99.

[2] Justin Solomon, Fernando de Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas. “Convolutional Wasserstein Distances: Efficient Optimal Transportation on Geometric Domains.” ACM Transactions on Graphics, vol. 34, no. 4, August 2015.

[3] Leticia Mattos Da Silva, Oded Stein, and Justin Solomon. “A Framework for Solving Parabolic Partial Differential Equations on Discrete Domains.” ACM Transactions on Graphics, vol. 43, no. 5, October 2024.

[4] Max Wardetzky, Saurabh Mathur, Felix Kälberer, and Eitan Grinspun. “Discrete Laplace Operators: No Free Lunch.” In Proceedings of the Fifth Eurographics Symposium on Geometry Processing, pages 33–37, Barcelona, Spain, 2007.

[5] Alec Jacobson et al. “libigl: A Simple C++ Geometry Processing Library.” Available at https://github.com/libigl/libigl, 2025.

[6] Paul J. Goulart and Yuwen Chen. “Clarabel: An Interior-Point Solver for Conic Programs with Quadratic Objectives.” arXiv preprint arXiv:2405.12762, 2024.

[7] Steven Diamond and Stephen Boyd. “CVXPY: A Python-Embedded Modeling Language for Convex Optimization.” Journal of Machine Learning Research, vol. 17, no. 83, pages 1–5, 2016.

Post author By nhu_tran
Post date October 20, 2025
Tags diffusion, pde, solver

Research

Gaussian Fluids

Post author By nhu_tran
Post date September 26, 2025

SGI Fellows: Haojun Qiu, Nhu Tran, Renan Bomtempo, Shree Singhi and Tiago Trindade
Mentors: Sina Nabizadeh and Ishit Mehta
SGI Volunteer: Sara Samy

Introduction

The field of Computational Fluid Dynamics (CFD) has long been a cornerstone of both stunning visual effects in movies and video games, and also a critical tool in modern engineering, for designing more efficient cars and aircrafts for example.
At the heart of any CFD application there exists a solver that essentially predicts how a fluid moves in space at discrete time intervals. The dynamic behavior of a fluid is governed by a famous set of partial differential equations, called the Navier-Stokes equations, here presented in vector form:

$$
\begin{align}
\frac{\partial \mathbf{u}}{\partial t}
+ (\mathbf{u} \cdot \nabla) \mathbf{u} &=
-\frac{1}{\rho} \nabla p
+ \nu \nabla^{2} \mathbf{u}
+ \mathbf{f},
\quad
\text{subject to}\quad \nabla \cdot \mathbf{u}
&=
0,
\end{align}
$$

where $t$ is time, $\mathbf{u}$ is the varying velocity vector field, $\rho$ is the fluids density, $p$ is the pressure field, $\nu$ is the kinematic viscosity coefficient and $\mathbf{f}$ represents any external body forces (e.g. gravity, electromagnetic, etc.). Let’s quickly break down what each term means in the first equation, which describes the conservation of momentum:

$\dfrac{\partial \mathbf{u}}{\partial t}$ is the local acceleration, describing how the velocity of the fluid changes at a fixed point in space.
$(\mathbf{u} \cdot \nabla) \mathbf{u}$ is the advection (or convection) term. It describes how the fluid’s momentum is carried along by its own flow. This non-linear term is the source of most of the complexity and chaotic behavior in fluids, like turbulence.
$\frac{1}{\rho}\nabla p$ is the pressure gradient. It’s the force that causes fluid to move from areas of high pressure to areas of low pressure.
$\nu \nabla^{2} \mathbf{u}$ is the viscosity or diffusion term. It accounts for the frictional forces within the fluid that resist flow and tend to smooth out sharp velocity changes. For the “Gaussian Fluids” paper, this term is ignored ($\nu=0$) to model an idealized inviscid fluid.

The second equation, $\nabla \cdot \mathbf{u} = 0$, is the incompressibility constraint. It mathematically states that the fluid is divergence-free, meaning it cannot be compressed or expanded.

The first step to numerically solving these equations on some spatial domain $\Omega$ using classical solvers is to discretize this domain. To do that, two main approaches have dominated the field from its very beginning:

Eulerian methods: which use fixed grids that partition the spatial domain into a (un)structured collection of cells on which the solution is approximated. However, these methods often suffer from “numerical viscosity,” an artificial damping that smooths away fine details;
Lagrangian methods: which use particles to sample the solution at certain points in space and that move with the fluid flow. However, these methods can struggle with accuracy and capturing delicate solution structures.

The core problem has always been a trade-off. Achieving high detail with traditional solvers often requires a massive increase in memory and computational power, hitting the “curse of dimensionality”. Hybrid Lagrangian-Eulerian methods that try to combine the best of both worlds can introduce their own errors when transferring data between grids and particles. The quest has been for a representation that is expressive enough to capture the rich, chaotic nature of fluids without being computationally prohibitive.

A Novel Solution: Gaussian Spatial Representation

By now you have probably come across an emerging technology in the field of 3D rendering called 3D Gaussian Splatting (3DGS). It attempts to represent a 3D scene using, not traditional polygonal meshes or volumetric voxels, but in a manner similar to an artist using paint brush strokes to approximate a colored picture, Gaussian Splatting uses 3D gaussians to “paint” a 3D scene. These gaussians are nothing more than a ellipsoid with a color value and opacity associated.

Think of each Gaussian not as a solid, hard-edged ellipsoid, but as a soft, semi-transparent cloud (Fig.2). It’s most opaque and dense at its center and it smoothly fades to become more transparent as you move away from the center. The specific shape, size, and orientation of this “cloud” are controlled by its parameters. By blending thousands or even millions of these colored, transparent splats together, Gaussian Splatting can represent incredibly detailed and complex 3D scenes.

Fig.1 Example 3DGS of a Lego build (left) and a zoomed in view of the gaussians that make it up (right).

Drawing from this incredible effectiveness of using gaussians to encode complex geometry and lighting of 3D scenes, a new paper presented at SIGGRAPH 2025 titled “Gaussian Fluids: A Grid-Free Fluid Solver based on Gaussian Spatial Representation” by Jingrui Xing et al. introduces a fascinating new approach to computational domain representation for CFD applications.

The new method, named Gaussian Spatial Representation (GSR), uses gaussians to essentially “paint” vector fields. While 3DGS encodes local color information into each gaussian, GSR encodes local vector fields. In this way each gaussian stores a vector that defines a local direction of the velocity field, where the magnitude of the vectors is greater closer the gaussian’s center, and quickly gets smaller when moving farther from it. In Fig.2 we can see an example of 2D gaussian viewed with color data (left), as well as with a vector field (right).

Fig.2 Example plot of a 2D Gaussian with scalar data (left) and vector data (right)

By splatting a lot of these gaussians and adding the vectors we can then define a continuously differentiable vector field around a given domain, which allows us to solve the complex Navier-Stokes equations not through traditional discretization, but as a physics-based optimization problem at each time step. A custom loss function ensures the simulation adheres to physical laws like incompressibility (zero-divergence) and, crucially, preserves the swirling motion (vorticity) that gives fluids their character.

Each point represents a Gaussian splat, with its color mapped to the logarithm of its anisotropy ratio. This visualization highlights how individual Gaussians are stretched or compressed in different regions of the domain.

Enforcing Physics Through Optimization: The Loss Functions

The core of the “Gaussian Fluids” method is its departure from traditional time-stepping schemes. Instead of discretizing the Navier-Stokes equations directly, the authors reframe the temporal evolution of the fluid as an optimization problem. At each time step, the solver seeks to find the set of Gaussian parameters that best satisfies the governing physical laws. This is achieved by minimizing a composite loss function, which is a weighted sum of several terms, each designed to enforce a specific physical principle or ensure numerical stability.

Two of the main loss terms are:

Vorticity Loss ($L_{vor}$): For a two dimensional ideal (inviscid) fluid, vorticity, defined as the curl of the velocity field ($\omega=\nabla \times u$), is a conserved quantity that is transported with the flow. This loss function is designed to uphold this principle. It quantifies the error between the vorticity of the current velocity field and the target vorticity field, which has been advected from the previous time step. By minimizing this term, the solver actively preserves the rotational and swirling structures within the fluid. This is critical for preventing the numerical dissipation (an artificial smoothing of details) that often plagues grid-based methods and for capturing the fine, filament-like features of turbulent flows.
Divergence Loss ($L_{div}$): This term enforces the incompressibility constraint, $\nabla \cdot u=0$. From a physical standpoint, this condition ensures that the fluid’s density remains constant; it cannot be locally created, destroyed, or compressed. The loss function achieves this by evaluating the divergence of the Gaussian-represented velocity field at numerous sample points and penalizing any non-zero values. Minimizing this term is essential for ensuring the conservation of mass and producing realistic fluid motion.

Beyond these two, to function as a practical solver, the system must also handle interactions with boundaries and maintain stability over long simulations. The following loss terms address these requirements:

Boundary Losses ($L_{b1}$, $L_{b2}$): The behavior of a fluid is critically dependent on its interaction with solid boundaries. These loss terms enforce such conditions. By sampling points on the surfaces of geometries within the scene, the loss function penalizes any fluid velocity that violates the specified boundary conditions. This can include “no-slip” conditions (where the fluid velocity matches the surface velocity, typically zero) or “free-slip” conditions (where the fluid can move tangentially to the surface but not penetrate it). This sampling-based strategy provides an elegant way to handle complex geometries without the need for intricate grid-cutting or meshing procedures.
Regularization Losses ($L_{pos}$, $L_{aniso}$, $L_{vol}$): These terms are included to ensure the numerical stability and well-posedness of the optimization problem.
- The Position Penalty ($L_{pos}$) constrains the centers of the Gaussians, preventing them from deviating excessively from the positions predicted by the initial advection step. This regularizes the optimization process, improves temporal coherence, and helps prevent particle clustering.
- The Anisotropy ($L_{aniso}$) and Volume ($L_{vol}$) losses act as geometric constraints on the Gaussian basis functions themselves. They penalize Gaussians that become overly elongated, stretched, or distorted. This is crucial for maintaining a high-quality spatial representation and preventing the numerical instabilities that could arise from ill-conditioned Gaussian kernels.

However, when formulating a fluid simulation as a Physics-Based optimization problem, one challenge that soon becomes apparent is that choosing good loss functions can be really tricky. The reliance on “soft constraints”, such as the divergence loss function, in the optimization process means that small errors can accumulate over time. Also, the handling of complex domains with holes (non-simply connected domains) is also a big challenge.
When running the Leapfrog simulation we can see that the divergence of the solution, although relatively small, isn’t as close to zero as it should be (ideally it should be exactly zero).

This “Leapfrog” simulation highlights the challenge of “soft constraints.” The “Divergence” plot (bottom-left) shows the small, non-zero errors that can accumulate, while the “Vorticity” plot (bottom-right) shows the swirling structures the solver aims to preserve.

Divergence-free Gaussian Representation

Now we want to extend on this paper a bit. Let’s focus on the 2D case for now. The problem of regularizing divergence is that the velocity vector field won’t be exactly divergence free. I can perhaps help this with an extra step every iteration. Define $\mathbf{J} \in \mathbb{R}^{2 \times 2}$ as a rotation matrix

$$
\mathbf{J} = \begin{bmatrix}
0 & -1\\
1 & 0
\end{bmatrix}.
$$

It is known that for any potential function $\phi: \mathbb{R}^{2} \to \mathbb{R}$, we have $\mathbf{J} \ \nabla\phi(\mathbf{x})$ divergence free:

$$
\nabla \cdot (\mathbf{J} \ \nabla \phi) = –
\frac{\partial }{\partial x} \frac{\partial \phi}{ \partial y} + \frac{\partial}{\partial y} \frac{\partial \phi}{\partial x} = 0.
$$

This potential function can be constructed with the same gaussian mixture by carrying a $\phi_{i} \in \mathbb{R}$ at every gaussian, and thus a continuous function having full support

$$
\phi(\mathbf{x}) = \sum_{i=1}^{N} \phi_i G_i(\mathbf{x}).
$$

similarly we can replace $G_i(\mathbf{x})$ with the clamped gaussians $\tilde{G}_i(\mathbf{x})$ for efficiency. Now we can construct a vector field that is unconditionally divergence-free

$$
\mathbf{u}_{\phi}(\mathbf{x}) = \mathbf{J} \ \nabla \phi(\mathbf{x})
$$

and the gradient of the potential can be computed

$$
\begin{align}
\nabla \phi(\mathbf{x}) &= \sum_{i=1}^{N} \phi_i \nabla G_i(\mathbf{x}) \\
&= – \sum_{i=1}^{N} \phi_i G_i(\mathbf{x})
\mathbf{\Sigma}_{i}^{-1}
(\mathbf{x} – \mathbf{\mu}_i).
\end{align}
$$

so the equation can be written as

$$
\mathbf{u}_{\phi}(\mathbf{x})
= \sum_{i=1}^{N}
\underbrace{\left(
– \phi_i \mathbf{J} \mathbf{\Sigma}_{i}^{-1}(\mathbf{x} – \mathbf{\mu}_i)
\right)}_{\text{per-gaussian vector}}
G_i(\mathbf{x})
$$

We want to somehow fit this $\mathbf{u}_{\phi}(\mathbf{x})$ to the vector field $\mathbf{u}(x)$ directly from the gaussian velocities

$$
\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{N} \mathbf{u}_i G_i(\mathbf{x}).
$$

This can be a simple loss

$$
\underset{\{\phi_i\}_{i=1}^{N}}{\arg\min} \frac{1}{d|\mathcal{D}|}
\int_{\mathcal{D}}
\|\mathbf{u}(\mathbf{x}) – \mathbf{J} \ \nabla \phi(\mathbf{x}) \|_{2}^{2}
\ \mathrm{d}\mathbf{x}
$$

and we use the similar approach of Monte Carlo way of sampling $\mathbf{x}$ uniformly in space.

Doing this, we managed to achieve a divergence much closer to zero, as shown in the image below.

References

Jingrui Xing, Bin Wang, Mengyu Chu, and Baoquan Chen. 2025. Gaussian Fluids: A Grid-Free Fluid Solver based on Gaussian Spatial Representation. In Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques Conference Conference Papers (SIGGRAPH Conference Papers ’25). Association for Computing Machinery, New York, NY, USA, Article 9, 1–11. https://doi.org/10.1145/3721238.3730620

Post author By nhu_tran
Post date September 26, 2025
Tags fluid-simulation, gaussian-splatting

Research

Star-Shaped Mesh Segmentation

Post author By nhu_tran
Post date August 15, 2025

SGI Fellows: Lydia Madamopoulou, Nhu Tran, Renan Bomtempo and Tiago Trindade
Mentor. Yusuf Sahillioglu
SGI Volunteer. Eric Chen

Introduction

In the world of 3D modeling and computational geometry, tackling geometrically complex shapes is a significant challenge. A common and powerful strategy to manage this complexity is to use a “divide-and-conquer” approach, which partitions a complex object into a set of simpler, more manageable sub-shapes. This process, known as shape decomposition, is a classical problem, with convex decomposition and star-shaped decomposition being of particular interest.

The core idea is that many advanced algorithms for tasks like creating volumetric maps or parameterizing a shape are difficult or unreliable on arbitrary models, but work perfectly on simpler shapes like cubes, balls, or star-shaped regions.

The paper by Hinderink et al. (2024) demonstrates this perfectly. It presents a method to create a guaranteed one-to-one (bijective) map between two complex 3D shapes. They achieve this by first decomposing the target shape into a collection of non-overlapping star-shaped pieces. By breaking the hard problem down into a set of simpler mapping problems, they can apply existing, reliable algorithms to each star-shaped part and then stitch the results together. Similarly, the work by Yu & Li (2011) uses star decomposition as a prerequisite for applications like shape matching and morphing, highlighting that star-shaped regions have beneficial properties for many graphics tasks.

Ultimately, these decomposition techniques are a foundational tool, allowing researchers and engineers to extend the reach of powerful geometric algorithms to handle virtually any shape, no matter how complex its structure. Some applications include:

Guarding and Visibility Problems. Star decomposition is closely related to the 3D guarding problem: how can you place the fewest number of sensors or cameras inside a space so that every point is visible from at least one of them? Each star-shaped region corresponds to the visibility range of a single guard. This is widely used in security systems, robot navigation and coverage planning for drones.
Texture mapping and material Design. Star decomposition makes surface parameterization much easier. Since each region is guaranteed to be visible from a point, its much simpler to unwrap textures (like UV mapping), apply seamless materials and paint and label parts of a model.
Fabrication and 3D printing. When preparing a shape for fabrication, star decomposition helps with designing parts for assembly, splitting objects into printable segments to avoid overhangs. It’s particularly useful for automated slicing, CNC machining, and even origami-based folding of materials.
Shape Matching and Retrieval. As mentioned in Hinderink et al. (2024), matching or searching for shapes based on parts is a common task. Star decomposition provides a way to extract consistent, descriptive parts, enabling sketch-based search, feature extraction for machine learning, and comparison of similar objects in 3D data sets.

In this project, we tackle the problem of segmenting a mesh into star-shaped regions. To do this, we explore two different approaches to try to solve this problem:

An interior based approach, where we use interior points on the mesh to create the mesh segments;
A surface based approach, where we try to build the mesh segments by using only surface information.

We then discuss results and challenges of each approach.

What is a star-shape?

You are probably familiar with the concept of a convex shape, that is a shape such that any two points on inside or on the boundary can be connected by straight line without leaving the shape. In simpler terms we can say that a shape is convex if any point inside it has a direct line of sight to any other point on its surface.

Now a shape is called star-shaped if there’s at least one point inside it from which you can see the entire object. Think of it like a single guard standing inside a room; if there’s a spot where the guard can see every single point on every wall without their view being blocked by another part of the room, then the room is star-shaped. This point is called a kernel point, and the set of all possible points where the guard could stand is called the shape’s visibility kernel. If this kernel is not empty, the shape is star-shaped.

The image above shows an example in 2D.

On the left we see an irregular non-convex closed polygon where a kernel point is shown along with the lines that connect it to all other vertices.
On the right we now show its visibility kernel in red, that is the set of all kernel points. Also, it is important to note that the kernel may be computed by simply taking the intersection of all inner half-planes defined by each face, i.e. the half-planes of each face whose normal vector points to the inside of the shape.

As a 3D example, we shown in the image bellow a 3D mesh of the tip of a thumb, where we can see a green point inside the mesh (a kernel point) from which red lines are drawn to every vertex of the mesh.

Point Sampling

Before we explore the 2 approaches (interior and surface), we first go about finding sample points on the mesh from which to grow the regions. These can be both interior and surface points.

Surface points

To sample points on a 3D surface, several strategies exist, each with distinct advantages. The most common approaches are:

Uniform Sampling: This is the most straightforward method, aiming to distribute points evenly across the surface area of the mesh. It typically works by selecting points randomly, with the probability of selection being proportional to the area of the mesh’s faces. While simple and unbiased, this approach is “blind” to the actual geometry, meaning it can easily miss small but important features like sharp corners or fine details.

Curvature Sampling: This is a more intelligent, adaptive approach that focuses on capturing a shape’s most defining characteristics. The core idea is to sample more densely in regions of high geometric complexity. These areas, often called “interest points,” are identified by their high curvature values, such as Gaussian curvature ($\kappa_G$) or mean curvature ($\kappa_H$). By prioritizing these salient features (corners, edges, and sharp curves), this method produces a highly descriptive set of sample points that is ideal for creating robust shape descriptors for matching and analysis.

Farthest-Point Sampling (FPS): This iterative algorithm is designed to produce a set of sample points that are maximally spread out across the entire surface. The process is as follows:

A single starting point is chosen randomly.
The next point selected is the one on the mesh that has the greatest geodesic distance (the shortest path along the surface) to the first point.
Each subsequent point is chosen to be the one farthest from the set of all previously selected points.

This process is repeated until the desired number of samples is reached.

In this project we implemented the curvature and FPS methods. However, none of these approaches gave us the desired distribution of points. So we came up with a hybrid approach: Farthest-Point Curvature-Informed sampling (FPCIS) . It aims to achieve a balance between two competing goals:

Coverage: The sampled points should be spread out evenly across the entire surface, avoiding large unsampled gaps. This is the goal of the standard Farthest Point Sampling (FPS) algorithm.
Feature-Awareness: The sampled points should be located at or near geometrically significant features, such as sharp edges, corners, or deep indentations. This is where curvature comes into play.

We start by selecting the highest curvature point in the mesh. Then, similar to FPS we compute the geodesic distances from this selected point to all other points using the Dijkstra algorithm over the mesh. However, instead of just taking the farthest point, we first create a pool of the top 30% farthest points, and from this pool we then select the one with highest curvature. The subsequent points are chosen in a similar manner to FPS.

The image above shows the sampling of 15 points using pure FPS (left) and using the FPCIS (right). we see that although FPS gives a pretty well distributed set of sample points we found that some regions that we deemed important where being left out. For example, when sampling 15 points using FPS on the cat model one of the ears was not being sampled. When using the curvature-informed approach we can see that the important regions are correctly sampled.

Interior points

To achieve a list of points that lie well inside the mesh and have good visibility out to the surface, we cast rays inward from every face and recording where they first hit the opposite side of the mesh. Concretely:

For each triangular face, we compute its centroid.
We flip the standard outward face normal to get an inward direction $d$.
From each centroid, we cast a ray along $d$, offset by a tiny $\epsilon$ so it doesn’t immediately self‐intersect.
We use the Möller–Trumbore algorithm to find the first triangle hit by each ray.
The midpoint between the ray origin and its first intersection point becomes one of our skeletal nodes—a candidate interior sample.

Casting one ray per face can produce thousands of skeletal nodes. We apply two complementary down‐sampling techniques to reduce this to a few dozen well‐distributed guards:

Euclidean Farthest Point Sampling (FPS)
- Start from an arbitrary seed node.
- Iteratively pick the next node that maximizes the minimum Euclidean distance to all previously selected samples.
- This spreads points evenly over the interior “skeleton” of the mesh.

k-Means Clustering
- Cluster all the midpoints into k groups via standard k-means in $\mathbb R^3$.
- Use each cluster centroid as a down‐sampled guard point.
- This tends to place more guards in densely populated skeleton regions while still covering sparser areas.

Interior Approach

Building upon the interior points sampled earlier, Yu & Li (2011) proposed 2 methods to decompose the mesh into star-shaped pieces: Greedy and Optimal.

1. Greedy Strategy

To achieve this, we solve a fast set-cover heuristic:

Repeatedly pick the guard whose visible-face set covers the most currently uncovered faces, remove those faces, and repeat until every face is covered.

For each guard $g$, we shoot one ray per face-centroid to decide which faces it can “see.” We store these as sets
$$ C_j=\{i: \text{face i is visible from } g_j\}$$
Let $U = { 0,1,…, |F|−1}$ be the set of all uncovered faces.
Repeat the following steps until $U$ is empty:
- For each guard j, compute the number of new faces it would cover:
  $$ s_j = | C_j \cap U | $$
- Pick the guard j with the largest $s_j$.
- Record one star-patch of size $s_j$.
- Remove those faces: $U = U \backslash C_j$.
Output: Get a list of guard indices (= number of star-shaped pieces) that covers the entire mesh.

2. Optimal Strategy

To get the fewest possible star-pieces from the pre-determined finite set of candidate guards, we cast the problem exactly as a minimum set-cover integer program:

Decision variables: Introduce binary $x_j \in \{0,1\}$ for each guard $g_j$.
Every face i must be covered by at least one selected guard (constraint) and we aim to minimize the total number of guards (thus pieces), hence together we solve:

The overall pipeline of our implementation is as follows:

We build the $|F| \times L$ visibility matrix in sparse COO format.
We feed it to SciPy’s milp, wrapping the “≥1” rows in a LinearConstraint and all variables as binary in a single Bounds object.
The solver returns the provably optimal set of $x_j$. The optimization problem can be stated as $\min \sum_{i=1}^m x_i$ subject to $x_i=0,1$ and $\sum_{j\in J(i)} x_j\ge 1, \forall i\in\{1,\dots,n\}$.

Experimental Results

We ran both methods on a human hand mesh with 1515 vertices and 3026 faces. The results are shown below:

Results from Greedy algorithm (left) and Optimal ILP solver (right)

Down-sampling method	Strategy	# of candidate guards	# of star-shaped pieces	Piece connectivity
FPS sampling	Greedy	25	13	No
	Greedy	50	9	Yes
	Optimal	25	9	No
	Optimal	50	6	No
K-means clustering	Greedy	25	8	No
	Greedy	50	8	No
	Optimal	25	7	No
	Optimal	50	6	No

Some key observations are:

Both interior approaches have no constraints to ensure connectivity among faces of the same guard, hence causing difficulties in closing the mesh of each piece or inaccurate calculation of final connected piece number.
Down-sampling methods: k‐means clustering tends to slightly better distribute guards in high‐density skeleton regions, hence having a smaller gap between greedy and optimal strategy.

Surface Approach: Surface-Guided Growth with Mesh Closure

One surface-based strategy we explored is a region-growing method guided by local surface connectivity and kernel validation, with explicit mesh closure at each step. The goal is to grow star-shaped regions from sampled points, ensuring at every expansion that the region remains visible from at least one interior point.

This method operates in parallel for each of the N sampled seed points. Each seed initiates its own region and grows iteratively by considering adjacent faces.

Expansion Process

Initialization:
- Begin with the seed vertex and its immediate 1-ring neighborhood—i.e., the three faces incident to the sampled point.

Iterative Expansion:
- Identify all candidate faces adjacent to the current region (i.e., faces sharing an edge with the region boundary).
- Sort candidates by distance to the original seed point.
- For each candidate face:
  - Temporarily add the face to the region.
  - Close the region to form a watertight patch using one of the following mesh-closing strategies:
    - Centroid Fan – Connect boundary vertices to their centroid to create triangles.
    - 2D Projection + Delaunay – Project the boundary to a plane and apply a Delaunay triangulation.
  - Run a kernel visibility check: determine if there exists a point (the kernel) from which all faces in the current region + patch are visible.
  - If the check passes, accept the face and continue.
  - If not, discard the face and try the next closest candidate.
- If no valid candidates remain, finalize the region and stop its expansion.

Results and Conclusions

The kernel check was not robust enough—some regions were allowed to grow into areas that clearly violated visibility constraints.
Region boundaries could expand too aggressively, reaching parts of the mesh that shouldn’t logically belong to the same star-shaped segment.
Small errors in patch construction led to instability in the visibility test, particularly for thin or highly curved regions.

As such, the final segmentations were often overextended or incoherent, failing to produce clean, usable star regions. Despite these issues, this approach remains a useful prototype for purely surface-based decomposition, offering a baseline for future improvements in kernel validation and mesh closure strategies.

References

Steffen Hinderink, Hendrik Brückler, and Marcel Campen. (2024). Bijective Volumetric Mapping via Star Decomposition. ACM Trans. Graph. 43, 6, Article 168 (December 2024), 11 pages. https://doi.org/10.1145/3687950

Yu, W. and Li, X. (2011), Computing 3D Shape Guarding and Star Decomposition. Computer Graphics Forum, 30: 2087-2096. https://doi.org/10.1111/j.1467-8659.2011.02056.x

Yusuf Sahillioğlu. “digital geometry processing – mesh comparison.” YouTube Video, 1:11:26. April 30, 2020. https://www.youtube.com/watch?v=ySr54PMAlGo.

Post author By nhu_tran
Post date August 15, 2025
Tags mesh segmentation, star-shape

Standard (post-log)	Continuous log-PDE	Proposed discrete-log
1. Start with the heat PDE for \(u\).	1. Start with the log-heat PDE for \(w=\log u\).	1. Start with the heat PDE for \(u\).
2. Discretize in space \(\Rightarrow\) ODE for \(u\).	2. Discretize in space for \(\nabla\) and \(\Delta\) \(\Rightarrow\) ODE for \(w\).	2. Discretize in space \(\Rightarrow\) ODE for \(u\).
3. Time-integrate in \(u\)-space.	3. Time-integrate in \(w\)-space and output \(w\).	3. Apply change of variables \(w=\log u\) \(\Rightarrow\) ODE for \(w\).
4. Compute \(w=\log u\) as a post-process.		4. Time-integrate in \(w\)-space and output \(w\).