Month: July 2026

During Week 3, we had the incredible opportunity to dive into an amazing project on 3D object reconstruction using medial axis. This journey wouldn’t have been possible without the guidance of Professor Kathryn Leonard and Professor Geraldine Morin. They were truly wonderful advisors, offering us their insight and support every step of the way. We’re so excited to share our thoughts and results with you all! Let’s start with some basic background.

by Stephanie Jung, Evelyn Zhu and Anja Milutinović

• 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:

where

The guidance term is approximated by

Therefore,

The complete local Jacobian becomes

Substituting this expression into the covariance equation gives

Under a convention in which $\lambda_t\geq0$,

is negative semidefinite. It therefore contributes local contraction in directions that affect the predicted measurement. This does not imply that the full covariance must decrease, since the unconditional dynamics and coupling between state components also contribute to its evolution.

For a nonlinear $F_t$, the exact guidance Jacobian also contains second-order derivatives. Define

Since

the exact derivative is

The first term contains the Hessians of the components of $F_t$. The locally linear approximation neglects these Hessian terms and retains the first-order contribution.

9. Splitting C and D

9.1 State Decomposition

The state is now partitioned into near and far components:

The conditional mean is partitioned as

and the conditional covariance is partitioned as

Similarly, the guided velocity and its Jacobian are written as

and

9.2 Block Moment Equations

The mean equations become

and

Expanding

block by block gives the near-state covariance equation

the near-far cross covariance equation

and the far-state covariance equation

The far-state covariance therefore does not evolve independently. Its evolution depends on the far-state dynamics, the near-far cross-covariance, and the local coupling between the near and far components of the guided velocity.

9.3 Observation Influence on $D_t$

Partition the predicted clean state as

Only the near component is directly observed, so

The predicted observation is therefore

and the measurement residual becomes

Although the observation depends only on $F_t^C$, the predicted near component may depend on both $C_t$ and $D_t$. Its Jacobian can be partitioned as

The DPS correction applied to the current near component is

The correction applied to the current far component is

This expression identifies the mechanism through which a partial observation can influence the unobserved component. Although D is not measured directly, the DPS correction can modify $D_t$ whenever perturbations in $D_t$ affect the predicted near state $F_t^C(U_t).$

In the linear Gaussian model, the analogous information pathway is determined by

Both the conditional mean of $D$ and the reduction in its covariance depend on this cross-covariance term. If the observed directions of $C$ are uncorrelated with $D$, then the measurement provides no information about $D$.

Furthermore, the mean update and covariance reduction are restricted to far-state directions selected by $\Sigma_{DC}A^\top$.

In the nonlinear flow-matching setting, the local Jacobian

plays an analogous role. It identifies the directions in the current far state that influence the predicted observation and are therefore accessible to measurement-based guidance.

The Gaussian cross-covariance and nonlinear denoiser Jacobian are distinct mathematical objects, but both characterize the coupling required for information to propagate from the observed component to the hidden component.

The preceding sections characterize how partial observations influence posterior guidance and hidden-state uncertainty. The following section documents a complementary fluid-specific evaluation framework for comparing reconstructed quasi-geostrophic fields.

Within this section, $C$ in the feature-loss definitions denotes the number of feature channels, rather than the near-state component introduced earlier.

10. QG-SSL Evaluation Framework

10.1 Objective

QG-SSL is a self-supervised encoder for comparing two-dimensional quasi-geostrophic (QG) vorticity fields. It learns spatial structure and short-term dynamics from real trajectories, without quality labels or generated samples.

10.2 Data and Preprocessing

We use the original paper’s released $64$ × $64$ filtered vorticity fields. These were obtained by spectrally filtering $512$ × $512$ QG simulations, as described in the paper. We use eddy and jet flows at Reynolds numbers $10^3$ and $10^4$.

Each combination contains $500$ trajectories with $196$ saved fields. We discard the spin-up portion and use frames $101$–$195$. The random initial conditions produce an early transient; the paper reports that the energy spectrum becomes self-similar only after $t=50$. Restricting training to this later regime avoids learning initialization artifacts. Each example is a pair $(x_t,x_{t+1})$; one step is approximately $0.5$ non-dimensional time units.

We train only on these $64$ × $64$ target fields, not on the paper’s $16$ × $16$ observations or on outputs from a generative model. Trajectories $0$–$399$ are used for training and $400$–$446$ for validation and for setting feature scales. All remaining trajectories are held out. For each physical configuration, we compute one scalar mean and standard deviation from the training fields and use them to standardize its inputs. Flow regime and Reynolds number are not given to the network.

The simulated square is periodic: opposite edges are connected, so a field leaving one edge re-enters from the other. We use this by cyclically rolling $x_t$ and $x_{t+1}$ by the same random multiples of 8 pixels. This changes the origin without changing their relative alignment. We then hide $50$% of $x_t$ in random 8 × 8 blocks and create a second view by rolling both the masked field and its mask again. Matching their global embeddings discourages dependence on absolute position.

The primary spatial distance nevertheless remains location-sensitive and penalizes translating only one of the two compared fields. The encoder receives two channels: the masked vorticity field and a binary visibility mask. At inference time the mask is entirely visible.

10.3 Encoder Architecture

The encoder is a compact, approximately $9$-million-parameter hierarchical transformer designed for a periodic domain. A stage is a group of transformer blocks operating at one fixed spatial resolution:

1. A $4$ × $4$, stride-4 convolution converts the input into a $16$ × $16$ grid with $96$ channels.
2. Three transformer stages operate at resolutions A $4$ × $4$, stride-$4$ convolution converts the input into a $16$ × $16$ grid with $96$ channels., 8 × 8, and $4$ × $4$, with $96$, $192$, and $384$ channels.
3. The stages contain $2$, $2$, and $4$ shifted-window attention blocks, respectively. Attention windows are $4$ × $4$, and cyclic shiftswrap across the domain boundary.
4. Each stage is projected to a $32$-channel spatial feature map. These three aligned maps form the primary representation.
5. The mean and standard deviation of every stage are concatenated and passed through an MLP to produce an optional $128$-dimensional global embedding.

Training uses a student encoder and an exponential-moving-average (EMA) teacher with the same architecture. The teacher sees complete fields; the student sees the masked current field. At each stage, separate spatial heads predict the teacher’s current features at hidden locations and its future features everywhere. Two MLP heads predict the corresponding global embeddings. The heads predict features, not pixels, and are discarded after training.

10.4 Training Losses

Setup

A training example is a pair of consecutive saved fields, $(x_t,x_{t+1})$. Let $\Omega$ be the field’s spatial domain, let $\mathcal{H}\subset\Omega$ be the randomly hidden region, and let $\mathcal{V}=\Omega\setminus\mathcal{H}$ be the visible region. The student encodes $x_t$ using only $\mathcal{V}$; the EMA teacher encodes the complete $x_t$ and $x_{t+1}$. Teacher outputs are treated as fixed targets.

At scale $s\in\{1,2,3\}$, $S_s$ is the student’s spatial feature map, $T_s^t$ and $T_s^{t+1}$ are the teacher’s current and future maps, $\Omega_s$ is the set of all spatial positions, and $\mathcal{H}_s\subset\Omega_s$ is the hidden region at that scale. The student global embedding is $g$, while the teacher embeddings are $u^t$ and $u^{t+1}$. A second student view, obtained by periodically translating both the field and its visible region by the same displacement, has embedding $\widetilde{g}$.̃

All feature errors are normalized coordinate by coordinate. For spatial maps $\widehat{F}$ and $F$ at scale $s$, each with 𝐶 channels, and any set $\mathcal{P}\subseteq\Omega_s$ of evaluated positions, define

For global vectors $\widehat{v},v\in\mathbb{R}^{D}$, define

Here $p$ indexes spatial positions, $c$ indexes spatial feature channels, $j$ and $k$ index global coordinates, $D=128$, and $\varepsilon=10^{-4}$. Each variance is taken over the training examples used to evaluate the loss; spatial variances also include all positions at the corresponding scale. Below, $\mathbb{E}$ averages over training pairs, hidden regions, and translations.

Masked Current-Feature Prediction

The head $P_s^t$ predicts the teacher’s current map from the student’s map. Only hidden positions are scored:

This forces the student to infer unobserved spatial structure from its visible context.

Future-Feature Prediction

A separate head $P_s^{t+1}$ predicts the teacher’s next-field map. The loss uses every position because the task is to predict the complete future:

Global-Feature Prediction

The heads $Q^t$ and $Q^{t+1}$ predict the teacher’s current and future global embeddings from the student’s current embedding:

This is the global counterpart of the two spatial prediction losses.

Periodic-Translation Consistency

The original and translated views should describe the same physical field, so their global embeddings are matched:

Variance Regularization

Let $\mathcal{G}$ be the set containing the embeddings 𝑔 and $\widetilde{g}$ from all training examples used to evaluate the loss. The variance penalty is

It prevents collapse to a constant embedding by requiring every coordinate to vary across examples.

Covariance Regularization

Let

The covariance penalty suppresses redundant correlations between distinct global coordinates:

Complete Objective

Optimization and Final Metric

We train with AdamW, batch size 128, learning rate $3\times10^{-4}$, weight decay $10^{-4}$, a cosine learning-rate schedule, bfloat16 arithmetic, and gradient clipping at 5. The teacher is updated after every batch with EMA decay 0.996. The current model is the teacher checkpoint after 10 epochs.

After training, we run the frozen encoder on the validation/calibration fields and compute the standard deviation of each feature channel across fields and spatial positions. When comparing two fields, their channel-wise feature differences are divided by these standard deviations. The primary distance is the resulting root-mean-square difference, averaged over the three spatial scales. This stops channels with naturally large numerical ranges from dominating.

For an ensemble of 𝑚 generated fields $x_1,\ldots,x_m$, a target field 𝑦, and a field distance 𝑑, we compute the energy score

The first term measures accuracy against the target. The subtracted pairwise term rewards ensemble diversity and therefore penalizes collapse. Lower is better. In later tables, “ES” means this score with 𝑑 replaced by the named distance.’

10.5 Compared Metrics and Baselines

Most metrics below define a distance $d(x,y)$ between two $64$ x $64$ fields. Pairwise experiments use that distance directly; a table entry ending in “ES” substitutes it into the energy score above. The paper’s full-cycle statistic instead compares a generated field with its given observation, while its ensemble-spread statistic has no target. Lower is better within any one metric, but absolute values cannot be compared between metrics because their scales differ.

QG-SSL Aligned

This is our primary distance. At each of the three spatial scales, corresponding feature vectors at the same location are compared after dividing each channel by its calibration standard deviation. We take the root-mean-square difference over channels and locations, then average the three scales. It therefore measures learned structure while retaining relative spatial alignment.

QG-SSL Global

This uses the optional 128-dimensional global embedding instead of the spatial maps. If $g_k(x)$ is coordinate 𝑘 and $\tau_k$ is its standard deviation on held-out real fields, then

The translation-consistency loss makes this representation approximately insensitive to the choice of spatial origin.

Pixel

The fields are flattened and corresponding grid values are compared directly. With $\sigma_p$ denoting the calibration standard deviation at grid position 𝑝,

where $\Omega$ is the set of all $64$ x $64$ positions. This baseline retains exact location and fine detail but has no learned notion of structure.

Spectrum

For each field we compute energy and enstrophy in 30 radial Fourier-wavenumber shells, take their logarithms, and concatenate them into a 60-dimensional vector. Each coordinate is standardized on held-out real fields and the distance is the root-mean-square vector difference. Because Fourier phase is discarded, this metric cannot locate structures in space.

DINOv2

DINOv2 is a self-supervised vision transformer pre-trained on natural images. We use its ViT-S/14 model without fine-tuning. A standardized vorticity field is clipped to three standard deviations, mapped to a grayscale image, resized to $224$ × $224$, and copied into three color channels. The model reduces this image to one 384-dimensional global embedding. We standardize each embedding coordinate on held-out QG fields and use the root-mean-square distance between embeddings.

LSiM

LSiM is a pretrained learned similarity metric for simulation fields. Each field is resized to $224$ × $224$, copied into three channels, and linearly mapped to $[0,225]$ using the joint minimum and maximum of the fields being compared. A five-scale convolutional network compares normalized feature maps using learned nonnegative channel weights; the square root of the summed multiscale error is the distance. The released model is used without training on our QG data.

Paper Reconstruction

This is the paper’s relative vorticity error. For a generated field 𝑥 and target 𝑦,

For an ensemble, the paper reports the mean of this error over members. Unlike an energy score, it contains no reward for ensemble diversity.

Paper Full-Cycle Consistency

Let 𝐴 be the paper’s observation operator, which filters a $64$ × $64$ field to the resolved $16$ × $16$ information, and let 𝑜 be the given observation. The cycle error of a generated field is

“Full-cycle” uses observations covering the full periodic domain; the partial version additionally restricts the comparison to observed regions. We retain this diagnostic only when 𝑜 is the actual observation used to condition the generated sample; it is not used as a distance between arbitrary pairs of fields. Ensemble results average the member-wise errors.

Paper Log-Energy

Let 𝐸(𝑥) be the 30-shell kinetic-energy spectrum of field 𝑥. The paper’s metric is

The paper averages this member-wise error for an ensemble.

Paper Log-Enstrophy

With 𝑍(𝑥) denoting the corresponding 30-shell enstrophy spectrum, this metric is

It is also averaged member-wise and, like log-energy, ignores Fourier phase.

Paper Ensemble Standard Deviation

For ensemble $X=\{x_1,\ldots,x_m\}$, the paper reports the mean pointwise spread

This has no target and is not a quality score by itself: either too little or too much spread can be wrong. In the controlled diversity benchmark, “ensemble-std discrepancy” is the relative difference between 𝑠(𝑋) and the spread of a clean reference ensemble.

10.5 Initial Metric Benchmarks

Temporal Neighborhood Self-Consistency

Each metric ranks seven real candidate trajectories relative to a query, then ranks the same trajectories h frames later. Spearman correlation measures preservation of the complete ranking. Each early, middle, or late origin panel contains 192 rankings across the four physical cases; the table averages available panels. Shuffling future identities gives correlations near zero.

h	Time	QG-SSL aligned	Paper reconstruction	Pixel	LSiM	DINOv2
1	0.5	0.968	0.970	0.971	0.906	0.605
4	2	0.899	0.884	0.875	0.715	0.332
8	4	0.847	0.782	0.764	0.595	0.285
16	8	0.783	0.598	0.554	0.387	0.156
32	16	0.681	0.461	0.376	0.267	0.182
80*	40	0.428	0.209	0.087	-0.015	0.096

Controlled Phase and Diversity Failures

We randomize Fourier phases while preserving magnitudes, destroying spatial structure without changing the spectrum. We divide this distance by that caused by a periodic translation. A high ratio means phase sensitivity with translation insensitivity; absolute values cannot be compared between metrics.

Metric	Phase-randomized	Translated	Ratio
QG-SSL aligned	1.161	1.294	0.90
QG-SSL global	0.877	0.035	24.78
DINOv2	3.249	0.639	5.09
LSiM	0.623	0.644	0.97
Pixel	1.597	1.668	0.96
Spectrum	2.56 × 10-6	7.24 × 10-7	3.54
Paper reconstruction	1.381	1.432	0.96
Paper log-energy	2.24 × 10-7	6.22 × 10-8	3.60
Paper log-enstrophy	3.40 × 10-7	9.03 × 10-8	3.76

M1-M4 Selection by Future Physical Utility

For each of the four physical cases, we take one 16-member ensemble from each of M1–M4. We evaluate every metric on each ensemble at the current time. Lower is better, so these four scores produce a current-time ranking of $M1-M4$

Independently, we lift every generated field to $512$ × $512$, evolve it with the matching QG solver, filter it back, and compute the ensemble’s future pixel-space ES against the ground-truth future field from the same QG trajectory that supplied its conditioning observation. This produces a reference ranking of $M1-M4$ by future physical utility.

There are six unordered pairs among four methods. “Pairwise order agreement” is the fraction of the $4$ x $6$ = $24$ method pairs for which the current-time ranking and the future ranking choose the same method as better; ties are excluded. “Correct winner” counts the physical cases in which the method with the lowest current-time score is also the method with the lowest future score.

Solver error on held-out real states is $0.003-0.010$ at $h=1$ and $0.017-0.064$ at h = 8. The future $M1-M4$ order is identical at $h=1, 4, 8$, so the table applies to each horizon.

Metric evaluated before rollout	Pairwise order agreement	Correct winner
QG-SSL aligned ES	1.000	4/4
Pixel ES	1.000	4/4
Paper reconstruction	1.000	4/4
LSiM ES	0.917	2/4
QG-SSL global ES	0.875	1/4
Paper cycle-consistency	0.875	1/4
Spectrum ES	0.750	1/4
Paper log-energy	0.708	0/4
Paper log-enstrophy	0.708	0/4
DINOv2 ES	0.625	0/4

11. Conclusion

This analysis examined how partial observations influence both observed and unobserved components of a high-dimensional physical state. The clean state was decomposed into a near component $C$, which is directly measured, and a far component $D$, which is not directly observed. The observation model was

The linear Gaussian model provides an exact description of this inverse problem. The posterior distribution of the far state is

with

and

These expressions show that information reaches the hidden component only through directions that are both visible to the observation operator $A$ and correlated with the far state through $\Sigma_{DC}$. If

then the observation does not change either the posterior mean or covariance of $D$. More generally, the number of far-state directions in which uncertainty can be reduced is limited by the rank of the observation and cross-covariance operators.

The Gaussian setting also yields an exact conditional score. The likelihood correction is

where $B_tu$ is the exact clean-state estimate and $\Gamma_t$ combines measurement noise with the remaining uncertainty in that estimate. This result provides an exact benchmark for diffusion posterior sampling.

The flow-matching analysis extends the same structure to a nonlinear denoising map $F_t$. The DPS-guided dynamics are

After local linearization around the conditional mean, the first two moments evolve approximately as

and

The local DPS contribution to the Jacobian is

where

After splitting the flow state into near and far components, the correction acting on the far coordinates is

Thus, a measurement applied only to the near component can influence the current far component whenever the predicted near state depends on the far coordinates. In the linear Gaussian model, this coupling is represented by $\Sigma_{DC}A^\top$. In the nonlinear flow-matching model, it is represented locally by

These are distinct mathematical objects, but they express the same central principle: information can propagate from an observed region to an unobserved region only through statistical or dynamical coupling between the two.

The analysis therefore provides a framework for interpreting the spatial reach of posterior guidance. It separates direct observational support from indirect information transfer and helps explain why reconstruction quality can deteriorate outside the observed region when the relevant coupling is weak.

12. Future Work

Future work should focus on testing how observation design, nonlinear coupling, and physical dynamics affect information transfer into unobserved regions.

12.1 Observation Bandwidth

Different observation operators $A$ measure different spatial scales and regions. Their effect can be studied through the Gaussian coupling

and the nonlinear sensitivity

Comparing these quantities with far-region reconstruction error could show how observation bandwidth and placement determine the spatial reach of guidance.

12.2 Numerical Validation

The locally linear moment equations should be compared with ensemble estimates from a trained model:

This would determine when the first-order approximation is accurate and when nonlinear Hessian terms become important.

12.3 Irregular Domains and Stability

The analysis should also be extended to irregular masks, such as land-ocean boundaries, while avoiding numerical artifacts near mask edges. Reconstructed states should then be evolved with the governing fluid solver to test long-time stability, forecast skill, energy and enstrophy behavior, and ensemble calibration.

12.4 Improved Posterior Guidance

The exact Gaussian model weights the residual using

whereas practical DPS often uses a simpler scalar weighting. Future work could develop ensemble-based or low-rank approximations to this effective covariance and compare them with standard DPS, conditional generative models, and PDE-constrained approaches.

13. References

[$1$] A. N. Suresh Babu, A. Sadam, and P. F. J. Lermusiaux, “Guided Unconditional and Conditional Generative Models for Super-Resolution and Inference of Quasi-Geostrophic Turbulence,” Journal of Advances in Modeling Earth Systems, vol. 18, no. 3, $e2025MS005324$, 2026. DOI: $10.1029/2025MS005324.$

[$2$] H. Chung, J. Kim, M. T. McCann, M. L. Klasky, and J. C. Ye, “Diffusion Posterior Sampling for General Noisy Inverse Problems,” International Conference on Learning Representations, 2023.

[$3$] A. N. Suresh Babu et al., “quasi-geostrophic-beta-plane-super-resolution,” GitHub repository, Models/samplers.py. This implementation contains the VP-SDE reverse sampler, DPS correction, Fourier/coarsening observation operator, and gappy-observation mask.

[$4$] R. A. Johnson and D. W. Wichern, Applied Multivariate Statistical Analysis, 6th ed., Pearson, Result 4.6, p. 160.

[$5$] Y. Polyanskiy and Y. Wu, Information Theory: From Coding to Learning, Cambridge University Press, data-processing inequality, Theorem 3.7(c).