SGI 2026

Project Members: Al Rahim Hossain, Stephanie Jung, Aleksa Milovanović, Reid Tang

Mentors: Otman Benchekroun, Ty Trusty

Introduction

Simulating deformable objects often requires optimizing over a large number of unknowns. For example, a 2D mesh with $n$ vertices has $2n$ positional variables: an $x$ and $y$ coordinate for each vertex. A 3D mesh similarly has $3n$ positional variables.

In a full order model, the deformed shape is represented by a vector $\mathbf{x}$ containing all vertex positions. It is found by minimizing the total energy of the object:

$\mathbf{x}^*=\underset{\mathbf{x}}{\text{argmin}}\left(E_{\mathrm{elastic}}(\mathbf{x})+E_{\mathrm{ext}}(\mathbf{x})\right)$

where $E_{\mathrm{elastic}}$ measures how much the object resists being deformed from its rest shape and $E_{\mathrm{ext}}$ represents external forces. The solution $\mathbf{x}^*$ is the state with the lowest total energy, where the elastic response and external forces are balanced. However, this solution may be computationally expensive as it treats every vertex position as an independently moving variable.

Reduced-Order Model

Reduced-order modelling reduces the size of this optimization by restricting the solution to a low-dimensional subspace. Instead of tracking every point independently, we choose a small set of meaningful deformation patterns. We then approximate the full configuration as $\mathbf{x} \approx \mathbf{x_0 + Bz}$ where $\mathbf{x_0}$ is the rest state of the object, the columns of matrix $\mathbf{B}$ represent the selected deformation patterns, and vector $\mathbf{z}$ contains their corresponding coefficients [1]. We can now solve for

$\mathbf{z}^* = \underset{\mathbf{z}}{\text{argmin }}E(\mathbf{x_0 + Bz})$

and the final full-space solution is given by $\mathbf{x}^* \approx \mathbf{x_0 + Bz}^*$.

The effectiveness of a reduced-order model depends heavily on the choice of the basis $\mathbf{B}$. Ideally, its columns should capture the important deformation patterns of the object using as few modes as possible. There are several ways to construct such a basis.

Proper Orthogonal Decomposition

Proper Orthogonal Decomposition constructs the basis from a dataset of observed or simulated deformations [2]. It finds the directions that best capture the variation present in the example shapes, making it a data-driven approach. It is found by solving

$\mathbf{B} = \underset{\mathbf{B}}{\text{argmin }} ||\mathbf{X – BB}^T\mathbf{X}||^2$ such that $\mathbf{B}^T\mathbf{B=I}$

Here, $\mathbf{X}$ is a matrix whose columns contain example displacement vectors from the rest configuration of the object. We use singular value decomposition, $\mathbf{X} = \mathbf{USV}^T$, to find the deformation patterns that best represent these examples. The columns of $\mathbf{U}$ are ordered from most to least important so if we want a reduced space with 20 modes, for example, we simply take the first 20 columns $\mathbf{B} = \mathbf{U}[:, :20]$.

A main drawback of POD is that it requires collecting a representative set of deformation examples beforehand. The resulting subspace can only represent deformation patterns present in, or similar to, this training data, so unseen motions may be captured poorly.

Linear Modal Analysis

Linear Modal Analysis builds the reduced basis from the object’s natural vibration modes around its rest shape [3]. The basis vectors are deformation modes obtained from the mesh’s stiffness and mass matrices, so the subspace is determined by the physics of the object rather than by example deformation data. These modes are found by solving the generalized eigenvalue problem

$\mathbf{H}\phi = \lambda\mathbf{M}\phi$.

Here, $\mathbf{H}$ is the elastic stiffness of the object and $\mathbf{M}$ accounts for its mass. Solving the eigenvalue problem gives the natural deformation patterns $\phi$, or modes, of the object. The lowest-frequency modes are usually the most important and are chosen as the columns of the reduced basis $\mathbf{B}$.

However, LMA is based on a linearization around the rest shape, so it works best for small deformations. Large rotations or strongly nonlinear deformations may not be represented well by a basis of linear vibration modes.

Actuation-Aware Subspaces

POD can capture large deformations, but it requires representative simulation data. LMA avoids this data collection, but its basis is not aware of the parameters driving the simulation. Actuation-aware subspaces build the reduced basis directly from how the object deforms as the actuation parameters change. We write

$\mathbf{x(a)} = \underset{ \mathbf{x} } { \text{argmin } } E(\mathbf{x, a})$

where $\mathbf{a}$ may represent quantities such as spring rest lengths, gravity, muscle activation, or collider parameters. The goal is to approximate how $\mathbf{x(a)}$ changes as $\mathbf{a}$ changes.

We can first approximate this relationship using a first-order Taylor expansion around a reference actuation $\mathbf{a}_0$:

$\mathbf{x}(\mathbf{a})\approx\mathbf{x}(\mathbf{a}_0)+\sum_i(a_i-a_{0,i})\frac{\partial\mathbf{x}}{\partial a_i}$.

The derivatives $\frac{\partial \mathbf{x}}{\partial a_i}$ describe the deformation caused by changing each actuation parameter. For example, if the object is driven by springs, each parameter could control the rest length of one spring. However, this is still a linear approximation and works best for small deformations near $\mathbf{a}_0$. Larger motions such as bending, rotation, and compression follow curved deformation paths that cannot be represented well using only fixed linear directions.

To capture these effects, we include second-order derivatives, $\frac{\partial^2 \mathbf{x}} {\partial a_i \partial a_j}$, for the second-order Taylor expansion:

$\mathbf{x}(\mathbf{a}) \approx \mathbf{x}(\mathbf{a}_0) + \sum_i (a_i-a_{0,i}) \frac{\partial \mathbf{x}}{\partial a_i} + \frac{1}{2} \sum_i\sum_j (a_i-a_{0,i})(a_j-a_{0,j}) \frac{\partial^2 \mathbf{x}}{\partial a_i \partial a_j}$

If the first derivatives describe the response to actuation, the second derivatives describe how that response changes as the actuation changes. These terms help capture deformation patterns that appear together during larger motions, allowing the reduced space to follow nonlinear deformation paths more closely.

The reduced subspace can therefore be formed by

$\mathbf{B} = \left[ \frac{\partial \mathbf{x}}{\partial a_1} \ \cdots\ \frac{\partial \mathbf{x}}{\partial a_p} \ \frac{\partial^2 \mathbf{x}}{\partial a_1^2} \ \cdots\ \frac{\partial^2 \mathbf{x}}{\partial a_i \partial a_j} \right]$

such that this space is tailored to the parameters driving the simulation.

This idea is similar to modal derivatives [1], which describe how linear vibration modes change and interact. Here, we apply the same idea to actuation parameters, so the basis is tailored to the forces that will actually drive the simulation.

Conclusion

Reduced-order models can make deformable simulations more efficient by representing motion with a small number of deformation patterns. In this project, we explored an actuation-aware basis that is built from how the equilibrium shape changes with the parameters driving the simulation. First-order derivatives capture the main response to actuation, while second-order derivatives help represent larger and more nonlinear deformations.

There are still several directions for future work. Hyper-reduction is needed to further reduce computational cost [4], and the current method relies on a PSD-projected Hessian for stable basis construction. The number of basis terms also grows with the number of actuation parameters. Future work could address these limitations and extend the method to more complex 3D scenes and different types of actuation.

References

[1] J. Barbič and D. L. James. Real-Time Subspace Integration for St. Venant-Kirchhoff Deformable Models. ACM Transactions on Graphics, 2005.

[2] L. Sirovich. Turbulence and the Dynamics of Coherent Structures. I. Coherent Structures. Quarterly of Applied Mathematics, 45(3), 561–571, 1987.

[3] A. Pentland and J. Williams. Good Vibrations: Modal Dynamics for Graphics and Animation. SIGGRAPH, 1989.

[4] C. Brandt, E. Eisemann, and K. Hildebrandt. Hyper-Reduced Projective Dynamics. ACM Transactions on Graphics, 2018.

Introduction

Swimming with a current feels effortless; swimming across it is exhausting — even over the same distance in meters. That asymmetry isn’t a fact about water, it’s a fact about which notion of “distance” you’re using. This week’s project builds exactly that idea into a triangle mesh: a metric that makes some directions cheaper to move through than others. Along the way we needed two core tools from discrete differential geometry — intrinsic triangulations and the anisotropic Laplacian — and had to make both of them survive a real 3D shark scan, not just a hand-drawn toy mesh.

Background

A quick glossary before the derivations:

• Mesh (V, F): vertex positions and the triangles connecting them.
• Halfedge: a directed half of an edge, belonging to one triangle; its twin is the matching half from the triangle on the other side. This lets code ask “what’s across this edge?” in constant time.
• Tensor field / metric S: a small symmetric, positive-definite matrix per triangle (2×2 in 2D, 3×3 on a 3D surface) describing how expensive movement is in each direction. S = I recovers ordinary Euclidean distance.
• Intrinsic triangulation: a different triangulation of the same surface, changing only connectivity — never vertex positions.
• Edge flip: the only move allowed — swap the diagonal of the quad formed by two triangles sharing an edge.
• Delaunay triangulation: every edge satisfies the empty-circumcircle property; avoids thin, needle-like triangles.
• Cotangent weight (wij): built from the two angles opposite an edge; decides both whether to flip it and what enters the Laplacian.
• Harmonic field: a function solving $Lf=0$ with fixed boundary values — “ink diffusing smoothly” between two fixed points.

Method: Intrinsic Triangulations

An intrinsic triangulation keeps a mesh’s lengths and angles — what an ant walking on the surface could measure — while reconnecting vertices via edge flips:

  i                 i
 /|\               / \
/ | \             /   \
m |  k   ====>   m --- k
\ | /             \   /
 \|/               \ /
  j                 j

This matters because it can fix badly-shaped triangles without adding a single vertex — useful whenever an algorithm only needs lengths and angles, not raw coordinates.

The special case is the Delaunay triangulation, equivalent to:
$$w_{ij} = \tfrac12(\cot\alpha+\cot\beta) \geq 0 \iff \alpha+\beta \leq 180°$$

The proof is short: $\cot\alpha+\cot\beta = \frac{\sin(\alpha+\beta)}{\sin\alpha\sin\beta}$​, and since $\alpha,\beta\in(0°,180°)$ the denominator is always positive — so the sign of wij​ depends only on sin(α+β). Empty-circumcircle and angle-sum are the same condition. In code, this is a queue: flip any edge with negative weight, re-check its neighbors, repeat until none remain.

Method: The Anisotropic Laplacian

To make distance direction-dependent, attach a metric $M$M to each triangle: $\|v\|_M=\sqrt{v^\top Mv}$​. M must be symmetric (for real eigenvalues/perpendicular axes) and positive-definite (so no length comes out negative or imaginary).

The pleasant surprise: the Delaunay proof above never used Euclidean angles specifically — only that they’re triangle angles in (0°,180°). Measuring α,β with M instead, the same algebra holds, so the criterion is unchanged in form. On a fixed test quad, switching S from identity to diag(1,16) flips the verdict entirely (w: +0.72 → -1.37) — same vertices, different verdict on whether the edge should exist.

The same $w_{ij}$​ assembles the Laplacian ($L_{ij}=-w_{ij}$, $L_{ii}=\sum_j w_{ij}$​), which must satisfy: rows sum to zero, symmetric, constants in the kernel. When all $w_{ij}\geq0$ (i.e., Delaunay), $L$ also satisfies the maximum principle — solutions to $Lf=0$ extremize only on the boundary. That guarantee is the entire reason the flip exists before building anything on top of $L$.

Implementation

The build order, including what broke:

Synthetic 2D test bed: HalfedgeMesh with flip_edge/flip_to_delaunay, validated on a hand-built non-Delaunay quad and a randomly-scrambled disc.

Toy 2D fish: solved a harmonic field (nose=0, tail=1) with the ordinary Laplacian, took its per-triangle gradient as the “body direction” $d$, and built $S=\alpha\,dd^\top+\beta(I-dd^\top)$, $\alpha<\beta$— cheap along the body, expensive across it.

Real 3D shark scan — the actual porting work:

Cleanup: weld coincident vertices, drop degenerate (zero-area) faces, since these send cotangent weights to nan and corrupt connectivity.

Automatic nose/tail: SVD on centered vertices finds the principal axis; nose and tail are its two extremes.

A genuinely 3D gradient: no single “perpendicular” exists on a surface, so the normal $n$ and cross product replace 2D rotation: $\nabla f = \frac{1}{2A}\sum_i f_i(n\times e_i)$, and the tensor grows to 3×3 with a transverse direction $t=n\times d$.

Hardening the flip queue: index by (face, corner) instead of hashing mutable Halfedge objects; return None instead of nan for degenerate triangles; reject flips that would duplicate an existing edge.

Results

The Polyscope UI toggles between the original mesh and the re-triangulated one, with a live counter of non-Delaunay edges before/after and total flips performed. Recomputing the direction field after flipping shows the payoff: geometry never moved, but connectivity reorganized to better respect the metric — the field still tracks the shark’s body, now on better-conditioned triangles.

By Evelyn Zhu, Mentor: Alvin Shi,

Have you ever wondered what 4D actually looks like? I do, I have always imagined what will 4D world looks like, and in this project I was able to take a small peak at this enigmatic world.

But before I show you, there’s one thing I need to introduce first. Here’s the problem: we live in three dimensions. Our eyes, our screens, our intuition are all 3D. A four-dimensional object simply doesn’t fit into the world we’re able to look at. So how could we possibly “see” one?

The trick is to translate between dimensions. There are two operations that do exactly that: one to build a dimension up, and one to peek a dimension down. They’re called extrusion and slicing, and once you have them, 4D stops being impossible to picture [1].

Extrusion: building a dimension up

Extrusion is something you already know, even if you’ve never called it that. Take a shape, make a copy, drag the copy along a brand-new direction, and fill in everything it sweeps through. A point dragged becomes a line. A square dragged becomes a cube. And a 3D object dragged along a fourth axis becomes a genuine 4D object (See figure 1).

Every rung of that ladder climbs one dimension. The secret is the direction you drag: as long as it points along a fresh axis the shape wasn’t already using, you gain a dimension. Drag a cube in a direction perpendicular to width, height, and depth all at once — a direction we can’t even point at and we’vejust built a shape in 4D.

Slicing: peeking a dimension down

Slicing is the opposite move. Instead of building up, we cut across. Pass a flat plane through a solid, and where it cuts, you get a cross-section one dimension lower. Slice an 2D triangle, you get a line segment; slicing a 3D cube, you will get a 2D square (see figure 2).

Now push that one rung higher: if we slice a 4D object with a “flat 3D plane” (a hyperplane), and the cross-section is a 3D shape, and this is finally something we can look at. That’s the whole idea. We can’t see the 4D object directly, but we can see its 3D slices.

Putting it together: watching 4D

So here’s how we can actually “see” 4D object in our 3D world. We take a 4D objects and we slice it (1 dimension down). One slice gives one 3D snapshot. But a single snapshot isn’t enough to feel a 4D shape, so we do something more: we slowly rotate the object through the fourth dimension and re-slice it at every step.

Here’s exactly what you’re looking at, step by step:

1. The object is 4D version of the stanford bunny. Every one of its points has four coordinates: the usual x, y, z, plus a fourth one we’ll call w. We can’t display it directly, because w points in a direction our world doesn’t have.
2. We fix w and take the cross-section. Holding w at a single value, we ask: which parts of the object live exactly there? That set of points forms a 3D shape — one slice. This is the slicing operation, applied to every little piece of the object and stitched together into the shape you see.
3. We rotate in 4D. Between frames, we turn the object by a small angle in a plane that involves w (here, the xw-plane). This is a true four-dimensional rotation, and it tips parts of the object that were “elsewhere along w” toward our fixed slice.
4. We re-slice, and repeat. Because the rotation carries new material through the slice, each frame’s 3D cross-section differs from the last. Play the frames in sequence and the shape appears to be growing, shrinking, splitting, and merging.

And this is the whole magic trick: we use extrusion to build a 4D object from 3D object then slicing to see it. Two simple moves and now a dimension we were never supposed to be able to view suddenly becomes something we can sit back and watch.

References

[1] Alvin Shi, Haomiao Wu, and Theodore Kim. 2025. Hyper-Dimensional Deformation Simulation. In ACM SIGGRAPH 2025 Conference Papers. ACM. https://doi.org/10.1145/3721238.3730730

By Shannon Cudworth, Mentors: Alek Fröhlich and Daniel Perazzo

Introduction

The goal of the project was to study statistical dependence from geometric perspective, where we define two random variables X, Y are defined on a surface $\mathcal{M}$, rather than in a Euclidean space. Specifically, we explore whether heat diffusion across a surface is dependent on the local curvature.

To study this relationship, we randomly sample a triangle face from a mesh and then construct one heat diffusion variable and one curvature variable for this sampled face. The former variable was defined using the Laplace-Beltrami operator, and for the latter variable we used the mean curvature value at the sampled face’s barycenter. We then employed the Hilbert–Schmidt Independence Criterion (HSIC) to investigate if faces with similar heat diffusion also have similar curvature.

To implement, we first sample faces with a probability proportional to their area, to avoid oversampling smaller triangles. After constructing the two aforementioned variables, we make curvature and heat kernel matrices to describe pairwise curvature and heat similarity, run a permutation test, and repeat for various sample sizes to estimate the power of the sample HSIC test.

Constructing the Heat Diffusion Variable

To approximate heat diffusion, we’re going to use the Laplace-Beltrami operator. For this, we need the cotangent Laplacian and mass matrix.

Cotangent Laplacian: A discrete approximation of the Laplace-Beltrami operator, used on triangle meshes. We define this matrix using the piecewise function:

Where $\alpha_{ij}$ and $\beta_{ij}$ are opposite angles for the edge (ij), and N(i) is the set of neighboring vertices to vertex i.

Mass Matrix: Represents how much of the mesh’s surface area is associated with each individual vertex, and was defined using gpy toolbox.

With both of these components, we can understand the geometric variation of the mesh’s surface (through the Laplacian), and how much the variation contributes to the overall surface area of the mesh (though the Mass Matrix).

To actually approximate the heat diffusion over the area, we need to solve the eigenvalue problem:

where L is the cotangent laplacian, M is the mass matrix, $\phi_i$ is the ith eigenvector, and $\lambda_i$ is the corresponding ith eigenvalue.

Using Scipy, we solve for the first 100 eigenvector-eigenvalue pairs, which represents the 100 smoothest solutions to the equation above. The eigenvector defines a basis function over the mesh vertices, and the eigenvalue is the rate of decay under heat diffusion.

Now, to construct our random variable X, we must recall that we have sampled a random face of the triangle mesh, and we have two arrays eigenvectors and eigenvalues, where:

eigenvectors[i]
eigenvalues[i]

returns the first 100 eigenvectors and eigenvalues at vertex i of the mesh. To utilize these arrays, we must directly evaluate the eigenvalue problem at the sampled face’s barycenter. Then for the sampled face with vertices i,j,k we calculate:

(eigenvectors[i] + eigenvectors[j] + eigenvectors[k]) / 3

which gives us a vector of the form:

where $b_i$ is the barycenter of the ith sampled face.

Then using the respective eigenvalue, we can construct a random variable $X_i$ that represents the heat diffusion from barycenter of the ith sampled face with:

Constructing the Curvature Variable

To construct the curvature variable, we calculate the mean curvature at each of the sampled face’s three vertices. Using libigl’s principal curvature value function, we return $k_1, \space k_2,$ the maximum and minimum curvature values at a vertex, respectively. We can then calculate the mean curvature value, defined as:

If the mean curvature is small or 0, we can interpret either a locally flat surface, or a saddle structure where $k_1$ and $k_2$ cancel each other out. A larger mean curvature indicates bending.

Then for the ith sampled face with vertices i,j,k, we can construct the random variable $Y_i$, which will calculate the curvature at the face’s barycenter, defined as:

where $H_i,\space H_j, \space H_k$ are the mean curvature values at vertex i, j, k, and $H(b_i)$ is the mean curvature at the barycenter of the ith sampled face.

Kernels

To check independence between our two random variables X and Y, we must construct kernel matrices, which will measure the similarity between pairs $X_i, \space X_j$ and $Y_i, \space Y_j$ (Schrab 19).

X variable: By construction of the laplacian, we can calculate the heat kernel $K_X$ by:

By definition:

Then we can say for any i,j:

Then,

which is the heat kernel formula (Mostowsky et al.).

Y variable: Unlike the heat diffusion variable X, we need to do a few more calculations to compute the curvature kernel $K_Y$, which we will do by using the Gaussian Kernel formula.

First, we define $\sigma_Y$, which represents the median pairwise distance for each $Y_i,\space Y_j$ pair. Then we can calculate the kernel $K_Y$ using the formula, for each i,j pair:

If $K_X, \space K_Y$ are large, then it follows that the ith and jth sampled faces have similar diffusion or similar curvature, respectively. Rather, if the kernels are small, then we can say the ith and jth sampled faces have different diffusion or curvature, respectively.

Hilbert-Schmidt Independence Criterion (HSIC)

The kernel matrices $K_X, K_Y$ allow us to determine if there are any similarities between pairs $X_i, \space X_j$ or $Y_i, Y_j$. Now, we want to examine that if pair $X_i, X_j$ are similar, if it is true that pair $Y_i, Y_j$ are similar as well. To acheive this, we compute the sample HSIC.

To ensure valid comparability between variables, we first make a centering matrix, defined as:

where $J_{nxn}$ is a matrix of all ones.

We can then define the centered kernels for X and Y as:

Then, by definition,

Here, a large sample HSIC value indicates dependence between heat diffusion and curvature (Schrab 25).

Permutation Testing and Power

One sample HSIC value isn’t enough to statistically determine whether or not X, Y have dependence. So, we must undergo a permutation test. To start, we define a null and alternative hypothesis:

$H_0: X, Y \space are \space independent \\ H_A: X, Y \space are \space dependent$

Our goal for this test is to simulate under the null hypothesis, and then see how likely our originally observed sample HSIC is to occur in those conditions. We permute the order of $K_Y$, recompute the sample HSIC, and repeat 200 times to create the null distribution. We then use the observed sample HSIC to calculate a p-value, which indicates whether or not we should reject or fail to reject the null hypothesis.

Taking a step further, we calculated the power of the sample HSIC test, which will give us the probability that, when X and Y are dependent, the sample HSIC will correctly detect that dependence. We do this by repeating the permutation test, and dividing the amount of times we rejected the null (p-value was less than 0.05, depending on significance level) by the total amount of retrials.

In our project, we compared the power of the sample HSIC test across sample sizes: 5, 10, 25, 50, 75, and repeated the permutation test 100 times per each sample size. We also compared the sample HSIC using the aforementioned heat diffusion kernel that was calculate with the Laplace-Beltrami operator, and another heat diffusion X variable and kernel, that was constructed using a numerical approximation and Gaussian Kernel method.

From the figure, we can see that, as the sample size increases, the power of sample HSIC test equals 1. This means that every trial rejected the null hypothesis, and the sample HSIC test is successful at detecting dependence between heat diffusion and curvature.

Possible Extensions of the Project

I think it would be fun to compare the sample HSIC across different meshes, and maybe how quickly the power of the sample HSIC test converges to one across different meshes. We could also compare different curvature formulas with heat diffusion, or vary the time that we allow heat diffusion to occur. Or, we can explore with different variables, such as letting X represent the geodesic field from a sampled face, and compare that to curvature.

Works Cited

Schrab, Antonin. “Optimal Kernel Hypothesis Testing.” University College London, 2025.

Mostowsky, Peter, et al. “The GeometricKernels Package: Heat and Matérn Kernels for Geometric Learning on Manifolds, Meshes, and Graphs.” Journal of Machine Learning Research, vol. 26, no. 276, 2025, pp. 1–14.

odedstein. sgi-introduction-course. GitHub, https://github.com/ddecatur/sgi-introduction-course.

By: Kyle Loh, Mahlet Girma, Max Dunitz

If you have ever had to flatten out a bag of chips so the checkout scanner would finally read it, you already understand the problem we worked on this week. A code printed on a surface that bends, folds, or curves is hard to read, and QR codes and barcodes carry one disadvantage here: they are obvious. In supply-chain tracking, that single visible code can be tampered with, causing a product to be moved into a market it was not meant for. Our project examines self-rectifying textures, which are patterns that look like random noise, but their autocorrelation contains a regular lattice. By measuring how that lattice bends in a photo, we can recover how the surface was deformed, without printing a single visible marker.

The usual way to undo a perspective distortion is to find distinct points (for example, the three eyes of a QR code) and use their pixel positions to solve for a homography. To stay hidden, self-rectifying textures move the landmarks into the autocorrelation of the image.

Definition of Autocorrelation: The autocorrelation measures how much an image resembles a shifted copy of itself. It is basically a dot product. For a 2D image $f: \mathbb{R}^2 \rightarrow \mathbb{R}$, the autocorrelation is defined as follows:

where $x$ is a pixel position and the lag (or shift) $\tau \in \mathbb{R}^2$. A peak is defined as the $\tau$ where $R_{f,f}(\tau)$ attains a local maximum. This operation is translation-invariant.

As autocorrelations are computationally expensive $(\mathcal{O}(N^2))$, we never compute them directly. Instead, we invoke the Wiener-Khinchin theorem, allowing us to compute the autocorrelation in terms of Fourier transforms ($\mathcal{F}$): $R = \mathcal{F}^{-1}(|\mathcal{F}(f)|^2)$ in $\mathcal{O}(N\log N)$.

Task 1 – Constructing Textures: We used a grayscale-version of Steamboat Willie with a fronto-parallel view as our base image.

Next, we superimposed three copies of the image, each offset by shift vectors $0, +u, +v \in \mathbb{R}^2$, using zero-padding. We call the superimposed image the base texture. The following autocorrelation plot would have six peaks, at $\{ \pm u, \pm v, \pm (u-v)\}$. This is known as the Fundamental Hexagon. Then, we apply a deformation $A$ to the superimposed image, causing the fundamental hexagon to become $\{ \pm Au, \pm Av, \pm A(u-v)\}.$

The Assignment Problem: Suppose you know the original shifts $\{ \pm u, \pm v, \pm (u-v)\}$ and have access to the base texture, but not the deformation. Unfortunately, the deformed hexagon $\{ \pm Au, \pm Av, \pm A(u-v)\}$ alone does not tell you which peak corresponds to $Au$ and which to $Av$, or their signs. However, if you can figure out this assignment problem, then you can immediately solve for $A$ with some linear algebra.

Given the symmetries, there are eight sign-and-order combinations. Using certain invariance properties, this can be reduced to six combinations. We find and apply six candidate inverse maps $A^{-1}$ and select the inverse map with the highest normalized cross-correlation score against the base texture. This is a computationally expensive solution, and for future work, we will explore more efficient methods for the assignment problem.

Task 2 – Full Rectification Pipeline: We introduce a method that no longer assumes a uniform linear deformation. Our goal is to recover a global inverse map $\phi^{-1}$ that maps the observed texture back to a fronto-parallel image.
For this experiment, we used blurred Gaussian white noise as our base image, constructing a texture by superimposing it with two shifted copies of itself and zero-padding. We then apply a homography deformation $A$.

Methodology: We sample small square patches from the grid. Then, we compute each patch’s autocorrelation and patches with “unreliable measurements” are discarded — criteria include: no valid fundamental hexagon detected or the inverse Jacobian’s determinant is a statistical outlier). For each patch, we detect the warped fundamental hexagon in each patch’s autocorrelation. Comparing these peaks with the known original shifts provides six candidates for the local inverse Jacobian $D\phi^{-1}(y_i)$.

Then, we construct a Delaunay triangulation mesh of the center points of the remaining patches. Two measurements are treated as neighbors when their centers share an edge in the mesh. To select on candidate Jacobian at each vertex, we assume that the physical deformation varies smoothly, so the correct matrices at neighboring vertices should be similar.

We resolve the candidate ambiguity using Minimum Spanning Tree propagation. First, we use a phase correlation on one patch with the original template to select one initial vertex and its actual local Jacobian $D\phi^{-1}(y_i)$. Suppose vertex $i$ has been already assigned the matrix $G_i^\ast$ while an adjacent vertex $j$ remains unassigned. For each of the six candidates at $j$, we compute the mismatch

and take the candidate with the minimum Frobenius norm. All edges from assigned vertices to unassigned vertices are stored in a priority queue, and we iteratively assign matrices to unassigned vertices adjacent to vertices with assignment.
Finally, we use a finite-element method that reconstructs a global inverse map whose gradient best matches the selected local Jacobians:

This determines the map up to a constant translation (the +C in integration), so we anchor the reconstructed mesh to the image boundary. We then regrid the deformed pixel values through the recovered map to produce a rectified image. For validation, we compare this result with the original template using the absolute error $|rectified – template|$.

Task 3 – Non-Planar Texture Rectification via Cylindrical Mapping: We extend our rectification pipeline to non-planar surfaces (particularly, a homography-deformed white-noise texture wrapped around a cylinder) utilizing the same MST and Finite Element Method pipeline from Task 2 to recover a global inverse deformation map.

Note: To model deformations on cylindrical geometry, spatial coordinates $(X,Y)$ in the planar texture domain are mapped to 3D surface coordinates $(X,Y, Z)$ for a cylinder of radius R=260px. To unroll or warp texture fields back into the reference domain, we map the angular coordinates $\theta = \text{atan2}(Y,X)$ back into 2D planar coordinates $(X,Y)$.

Conclusion: Our rectification pipeline has a clear strength: most internal vertices (from the Delaunay triangulation) have sub-pixel error in the rectification and are accurate intensity-wise. However, the pipeline also has two clear weaknesses: (1) there are larger intensity errors at the edges and corners, most likely from our failure to collect as many “good” observations near the edges and corners of the texture and (2) the dependency on accessing the template at least once in the process. Future work would entail resolving these two issues.

References:
[1] Bencheikh, Ismail, et al. “Autocorrelation-based Fiducial Markers for Traceability.” Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision. 2026.

• What is a Medial axis ?

The medial axis is the “skeleton” of a shape. In 2D, the medial axis of a shape is the set of all points that are equidistant from at least two points on the boundary, i.e. the centers of the largest circles you can inscribe inside it. Paired with the radius at each point, it’s a compact yet near-complete description of a shape: sweeping those inscribed disks/balls back along the axis reproduces the original subject. For a straight rod, the medial axis is a line, for a bent tube it’s a curve, and for a torus it’s a circle. The axis really captures a shape’s essential structure in far fewer numbers than the full surface. 

One of the nice properties of the medial axis is that it preserves topology and geometry of the original shape. Given the medial axis of a 2D shape and the corresponding radius of each medial axis point, the original shape can be reconstructed.

• How do we Compute Medial axis using Voronoi diagrams?

The medal axis is defined with respect to the continuous boundary of a shape and computing exact medial axis is then difficult. It’s then easier to approximate it using Voronoi diagrams of sample points on the boundary.

If we have a set of sample points P = {p1, p2, … , pn}, the Voronoi cell of the sample point pi is a set of all points that are closer to pi, then any other pj:

The Voronoi diagram is a collection of all Voronoi cells.

The medial axis and Voronoi diagram’s edges both have in common that they are defined by points that are equidistant from multiple boundary locations. The algorithm to compute medal axis then becomes:

1. Sample the boundary of the given object
2. Construct the Voronoi diagram of the sample points
3. Keep only the Voronoi vertices that are inside the object
4. Connect the neighboring vertices with existing Voronoi edges

With this foundational understanding of medial axes and Voronoi diagrams under our belts, we were ready to tackle the core of our Week 3 project. We wanted to see if we could take these 2D ‘skeletons’ and use them to solve a much more complex problem: bringing flat images to life by reconstructing a full 3D object from them.

Reconstructing 3d object from its 2d images using medial axis

Method 1:

• We reconstruct a 3D object from a set of binary silhouettes taken from viewpoints spread over a sphere. For each view we compute an orthographic projection of the mesh to get a silhouette, then extract its 2D medial axis and radius via the Voronoi diagram of the boundary. We then fuse the views in 3D by carving a voxel grid (keeping only voxels whose projection lands inside every silhouette) and thin the result with a 3D skeletonization to recover the object’s 3D medial axis. Finally, we rebuild the surface by sweeping balls of the recovered radius along that axis, which reproduces the original shape. 

Pipeline Visualization:

Method 2:

• Each point on a 2D medial axis represents the center of a maximal inscribed circle. All of the possible centers lie along the camera’s viewing direction. So then, every circle naturally defines a cylinder in 3D space whose axis is the same direction as the camera viewing direction and whose radius equals the radius of the medial circle. Each orthographic image generates an entire set of cylinders, or as we called them “swept cylinders”.

• We position the “swept cylinders” by setting the medial axis from each camera view in a world space and using the camera parameters (orthographic scale, location with respect to the object, viewing direction) to position each 2D medial axis about the world origin such that it reflected the camera view about the original 3D object. 

Intersecting the “sweep cylinders” from multiple views in 3D space, gives us a feasible reconstruction region:

Method 2a:

To compute the radii of each sphere we use a signed distance function for a cylinder. Since the center point can belong to multiple cylinders, in each “swept cylinder”, we take the one where the point has a maximum SDF. The sphere must be contained in the intersection of each chosen “swept cylinder”, so the final radius is taken as minimum.

Method 2b:

In an alternative method, to compute the radii of each sphere we use the signed distance field of the binary image for each view, which evaluates a point’s distance from the nearest boundary of the shape. We sample points in 3D space at regular intervals and project them down into each view’s 2D image plane, then compute the signed distance for the point in that view. This method is functionally equivalent to using a 3D signed distance field in world space; it evaluates a point’s distance from the nearest boundary of the swept cylinder for that camera view. Then, the radius of a sphere centered at a point is the minimum signed distance across all camera views, since the sphere must be contained in the intersection of all swept cylinders.

SGI Mentor: Akhil Sadam

SGI Fellows: Santoshi Yadagiri, Pietro Palombini

1. Introduction

Many geophysical inverse problems require reconstruction of a high-dimensional physical state from observations that are incomplete, noisy, or available only over part of the spatial domain. In oceanic and atmospheric applications, observations may be coarse, sparse, or separated by large unmeasured regions. The objective is not only to produce a plausible reconstruction, but also to characterize how much information the measurements provide about the unobserved portion of the state.

Let the clean physical state be partitioned as

where $C$ denotes the observed or near component and $D$ denotes the unobserved or far component. Measurements are assumed to depend directly only on $C$

Equivalently, defining

the observation model becomes

The operator $A$ determines which spatial directions or scales of the near state are measured, while $R$ describes the measurement-noise covariance. No component of $D$ is observed directly. Information about the far region can therefore be recovered only through statistical or dynamical coupling between $C$ and $D$.

This problem is motivated by the reconstruction of quasi-geostrophic flow fields from coarse, sparse, and gappy observations. Diffusion-based generative models can provide probabilistic reconstructions without explicitly solving for a single deterministic inverse. However, recent numerical results show that guided unconditional methods such as diffusion posterior sampling may have difficulty propagating observational information into unobserved regions [$1$]. This raises the question of how information from a partial observation influences uncertainty and reconstruction quality outside the observed region.

First, a linear Gaussian model is used to derive the conditional distribution, uncertainty reduction, and information transfer exactly. Second, the same structure is extended to a nonlinear flow-matching model through local linearization of the denoising map and the guided velocity field.

2. Background

2.1 Geophysical Inverse Problems

An inverse problem seeks to infer an unknown physical state from indirect measurements. In the present setting, the forward observation process maps the clean state $U_0$ to data $Y$:

For a linear observation operator, this reduces to

The inverse problem is generally ill posed because multiple clean states may produce similar observations, particularly when the data are low resolution, spatially incomplete, or noisy. Rather than selecting a single reconstruction, a probabilistic method seeks the posterior distribution

This posterior represents both the states that are compatible with the measurement and the remaining uncertainty after conditioning.

In the near-far decomposition, the corresponding far-state posterior is

A principal objective is to determine when this distribution differs meaningfully from the prior distribution of $D$. If it does not, then the observation provides no usable information about the unobserved region.

2.2 Diffusion-Based Posterior Sampling

Diffusion models introduce a family of noisy states connecting clean data to an approximately Gaussian terminal distribution. A standard linear forward-noising model is

The reverse process is governed by a score function. The unconditional score is

while the conditional score is

Bayes’ rule gives the score decomposition

The first term is given by an unconditional generative model. The second term incorporates the measurement and directs sampling toward states that are compatible with the observation.

The likelihood term is usually intractable because the observation is defined on the clean state $U_0$, while the reverse process evolves through the noisy state $U_t$. Diffusion posterior sampling approximates this term by applying the observation operator to a denoised estimate of the clean state and differentiating the resulting data-fidelity loss [$2$].

2.3 QG Sampler and Observation Operator

The quasi-geostrophic (QG) implementation provides a concrete realization of the likelihood-guided reverse process. The sampler assumes a variance-preserving stochastic differential equation and supports unconditional sampling, conditional sampling, classifier-free guidance, SDEdit, and diffusion posterior sampling [$3$].

At diffusion time $t$, the unconditional network predicts a noise field

Using the variance-preserving parameterization, the corresponding estimate of the clean field is

The implementation then applies a low-resolution observation operator to this clean-state estimate. In Fourier space, the field I first multiplied by a Gaussian filter of the form

where $s$ is the coarsening scale, $\Delta x$ is the grid spacing, and $k_r$ is the radial wavenumber. Additional spectral cutoffs remove modes above the coarse-grid Nyquist limit. The filtered field is transformed back to physical space, average pooled on an $s\times s$ grid, and repeated to the original resolution.

The resulting operator may be represented abstractly as

The implementation also includes an optional gappy-observation mask that sets selected spatial swaths of the coarsened field to zero [$3$].

This construction makes the spatial bandwidth of the observation operator explicit. The Gaussian spectral filer, hard spectral cutoff, pooling scale, and spatial mask collectively determine which directions of the high-resolution state are visible to the likelihood.

2.4 Observation Support and Information Propagation

For a partial observation,

so the likelihood depends directly only on the near component. The spatial support of bandwidth of the observation kernel is encoded by $A$. A wider kernel may observe more spatial directions, while a restricted kernel may leave large subspaces unmeasured.

Direct observation is not the only mechanism through which information can propagate. If the near and far components are statistically correlated, observing $C$ can reduce uncertainty in $D$. In a nonlinear generative model, an analogous effect occurs when the predicted observed region depends on the hidden coordinates of the current state. The linear Gaussian model isolates this mechanism in a form that can be derived exactly.

3. Linear Gaussian Model

3.1 State and Observation Model

Let

with

Assume a centered Gaussian prior

with block covariance

The diagonal blocks are the marginal covariances

and the off-diagonal blocks are the cross-covariances

The observation is

where

and $N$ is independent of $U_0$. With

the observation equation is

3.2 Gaussian Conditioning

The conditioning calculation uses the following standard result. Let $X$ and $Z$ be jointly Gaussian with zero means,

Then [$4$]

Because $Y$ is a linear function of independent Gaussian variables, the pair $(U_0,Y)$ is jointly Gaussian. Its observation covariance is

and the state-observation cross-covariance is

Applying the Gaussian conditioning formula [$4$] gives

where

and

The posterior mean depends linearly on the realized observation $y$. The posterior covariance does not depend on the particular observed value because the model is linear Gaussian.

4. Information Transfer from C to D

4.1 Far-State Posterior

The cross-covariance between the far state and the observation is

Applying Gaussian conditioning directly to $(D,Y)$ gives

with posterior mean

and posterior covariance

These expressions separate the two mechanism that determine information transfer. The observation operator $A$ selects directions of the near state, while the cross-covariance $\Sigma_{DC}$ determines which of those observed directions are correlated with the far state.

If

then

In this case, conditioning on the observation does not change the distribution of $D$.

4.2 Reduction in Far-State Uncertainty

Define

Then

For any $v\in\mathbb{R}^{n_D}$,

Therefore,

and

Thus, conditioning cannot increase posterior uncertainty in any linear direction of the far state.

Let

Since

the rank of the covariance reduction satisfies

Consequently, $m$ measurements can reduce far-state uncertainty in at most $m$ independent directions. In addition,

Both the posterior-mean update and covariance reduction are therefore restricted to directions selected by $\Sigma_{DC}A^\top$.

4.3 Mutual Information

For a Gaussian random variable $X\in\mathbb{R}^k$ with covariance $S\succ0$, the differential entropy is

Applying

causes the constant terms to cancel and gives the determinant ratio below.

The reduction in uncertainty can also be expressed through mutual information. For Gaussian variables.

This quantity measures the information about the far state contained in the observation. It vanishes exactly when

Therefore, the observation carries information about $D$ only through the near-far covariance directions that are also visible to $A$.

Define the normalized covariance reduction

Then

and

The eigenvalues of $K_D$ quantify the fractional uncertainty reduction along informative far-state directions.

4.4 Information Under Forward Noising

For the forward-noised far state

the unconditional and conditional covariances are

and

Hence,

For $\sigma_t>0$, define the signal-to-noise ratio

Scaling by a nonzero constant preserves mutual information, so

To compare two signal-to-noise ratios, suppose

and define

Let

where $\varepsilon_1$ and $\varepsilon’$ are independent standard Gaussian variables independent of $(D,Y)$. Expanding $Z_2$ gives

The Gaussian noise in parentheses is standard, so $Z_2$ has the lower-signal-to-noise law and is obtained from $Z_1$ by adding noise. Therefore,

is a Markov chain. By the data-processing inequality [5],

Hence $I(D_t;Y)$ is nondecreasing in $\lambda_t=\alpha_t^2/\sigma_t^2$. Equivalently, forward noising cannot increase the information about $Y$ available in the far state. This gives an information-theoretic interpretation of the reverse process. As the signal-t0-noise ratio increases, the process can progressively recover the information present in the conditional clean-state distribution.

5. Gaussian Score and DPS

5.1 Forward-Noised Posterior

The forward-noised state is

with independent of $(U_0,Y)$. Since $U_0\mid Y=y$ is Gaussian,

The far marginal satisfies

5.2 Exact Conditional Score

For a Gaussian random variable

the score is

Applying this identity to the conditional distribution of $U_t$ gives

Unconditionally,

where

The unconditional score is therefore

The score of the far marginal is

This far-marginal score is not generally equal to the far coordinates of the full-state score. The full-state score may depend jointly on the current near and far coordinates, while the marginal score depends only on $d$.

5.3 Exact Denoiser and Noisy Likelihood

The pair $(U_0,U_t)$ is jointly Gaussian, with

and

Conditioning gives

where

and

The exact unconditional denoiser is therefore

Given $U_t=u$, the observation distribution is

where

The matrix $\Gamma_t$ is the effective observation covariance at time $t$. It combines the original measurement noise with the remaining uncertainty about the clean state after conditioning on $U_t=u$.

Differentiating the exact noisy likelihood gives

5.4 Exact Gaussian DPS Identity

Combining the unconditional score with the noisy likelihood gradient yields

In the linear Gaussian model, this identity is exact. It has the same structure as diffusion posterior sampling:

1. estimate the clean state from the current noisy state
2. apply the observation operator to the estimate
3. compute the measurement residual
4. propagate the residual back to the current state.

A standard point-estimate DPS approximation replaces the effective covariance $\Gamma_t$ with a simpler measurement weighting and replaces the exact linear denoiser $B_tu$ with a learned denoising estimate. The Gaussian analysis therefore provides both an exact benchmark and a direct motivation for the nonlinear flow-matching extension.

5.5 Exact Reverse-Process Distribution

Let $\widehat{U}$ denote the output of the exact continous reverse process initialized from the exact terminal distribution, driven by the exact conditional score, and simulated without numerical error. Because this is the exact reverse process,

In the linear Gaussian setting, this gives

Let $P_D$ denote projection onto the far-state coordinates and define

Then

Thus, under exact sampling assumptions, the reverse process reproduces both the complete conditional distribution and its far-state marginal.

5.6 Relation to the QG DPS Implementation

The exact Gaussian identity can be compared directly with the quasi-geostrophic sampler. The exact linear denoiser

is replaced by the learned clean-state estimate

The theoretical observation operator $H$ is replaced by the filtered and coarsened operator

The implementation defines a mean-squared measurement loss

Automatic differentiation is used to compute

The implemented DPS correction is

with time-dependent scale

The correction is added to the model score during both the reverse updated and the Langevin correction steps [3].

This implementation has the same computational structure as the exact Gaussian likelihood correction:

1. estimate the clean state from the current noisy state
2. apply a measurement operator to the estimate
3. compare the predicted measurement with the observed field
4. differentiate the discrepancy with respect to the current state
5. add the resulting correction to the unconditional score

The principal difference is the likelihood weighting. In the exact Gaussian model, the residual is weighted by the full effective covariance

The implementation instead uses the scalar factor

Thus, the Gaussian calculation identifies the matrix-valued uncertainty correction that is approximated in practice by a time-dependent scalar weighting. It also clarifies how the observation operator’s spectral bandwidth and spatial mask determine which state-space directions receive guidance.

6. Flow-Matching Extension

The linear Gaussian model provides an exact setting in which the posterior distribution and conditional likelihood gradient can be computed analytically. The clean-state estimator is linear, and the noisy likelihood gradient takes the form

where $B_tu$ is the conditional estimate of the clean state and $\Gamma_t$ accounts for measurement noise and the remaining uncertainty in that estimate.

The same structure can be extended to a nonlinear flow-matching model by replacing the linear clean-state estimator with a time dependent denoising map.

6.1 Denoising Map

Let $U_t\in\mathbb{R}^n$ denote the state at time $t$. The unconditional flow is controlled by

where $v_t:\mathbb{R}^n\rightarrow\mathbb{R}^n$ is the unconditional velocity field.

Assume that a time-dependent map $F_t$ predicts the clean state $U_0$ from the current state

The map $F_t$ is treated as given. No separate differential equation for $F_t$ is introduced.

The observation model is

where $H$ is the observation operator and $R\succ0$ is the measurement-noise covariance.

Applying the observation operator to the predicted clean state gives the predicted measurement

The corresponding measurement residual is

This residual measures the disagreement between the observed data and the observation predicted from the current state.

6.2 DPS Guidance

Diffusion posterior sampling introduces an observation-dependent correction through the gradient of the likelihood. Using $F_t(U_t)$ as a point estimate of the clean state, the Gaussian likelihood is

The DPS guidance term is defined as

where $\lambda_t$ controls the time-dependent strength and sign convention of the guidance.

Define the Jacobian of $F_t$ by

Applying the chain rule gives

Equivalently,

where

The residual is evaluated in observation space. Multiplication by $H^\top$maps it to the clean-state space, while $J_{F_t}(U_t)^\top$propogates the correction from the predicted clean state back to the current state.

Although the measurement is defined on $U_0$, the likelihood gradient can guide intermediate states throughout the flow.

6.3 Guided Flow

After adding the DPS correction, the state evolves according to

Define the full guided velocity as

The guided dynamics are therefore

The complete observation-guided update follows the sequence

The current state is first mapped to a clean-state prediction. This prediction is passed through the observation operator, and the resulting residual is propagated back through the denoising map.

7. Moment Evolution

The DPS equation describes the evolution of an individual trajectory. To characterize the conditional distribution of trajectories given $Y=y$, consider the conditional mean and covariance.

7.1 Conditional Mean

Define the conditional mean as

Since each trajectory satisfies

differentiating the conditional expectation gives

Substituting the definition of the guided velocity yields

When either $v_t$ or $F_t$ is nonlinear, the expectation depends on the full conditional distribution of $U_t$, rather than only on $m_t$. The mean equation is there not closed in general.

7.2 Conditional Covariance

Define the conditional covariance by

Equivalently,

Let

Then

Differentiating

gives

This covariance equation is also exact, but it is not closed for a nonlinear guided velocity field.

The flow-matching dynamics considered here are deterministic. Consequently, there is no separate process-noise covariance term. The conditional covariance instead comes from the distribution of the initial state and from conditioning on the observation.

A local approximation is therefore required to obtain closed evolution equations for $m_t$ and $\Sigma_t$.

8. Local Linearization

8.1 Local Approximation

Locally linearize the full guided velocity around the conditional mean:

where

Because

the Jacobian can be decomposed as

where

and

Since

the locally closed mean equation becomes

Expanding the guidance term gives

The covariance equation becomes

These equations have the same structure as the moment equations for a locally linear deterministic system. The effective Jacobian additionally includes an observation-dependent DPS contribution.

8.2 DPS Jacobian

To identify this contribution, locally linearize the denoising map: