{"id":749,"date":"2026-08-13T18:13:37","date_gmt":"2026-08-13T18:13:37","guid":{"rendered":"https:\/\/summergeometry.org\/sgi2026\/?p=749"},"modified":"2026-08-25T16:12:01","modified_gmt":"2026-08-25T16:12:01","slug":"inverse-geodesic-design","status":"publish","type":"post","link":"https:\/\/summergeometry.org\/sgi2026\/inverse-geodesic-design\/","title":{"rendered":"Inverse Geodesic Design"},"content":{"rendered":"\n

Mentor(s):<\/strong> Dr. Ahmed Mahmoud<\/a> [MIT], Ana Dodik<\/a> [MIT]<\/p>\n\n\n\n

Fellow(s): <\/strong>Gokul Adithya Suresh<\/a> [RVCE], Yufan (Diana) Hu [Macalester College]<\/p>\n\n\n\n\n\n

Volunteer<\/strong>: Ruyu Yan<\/a> [Princeton University]<\/p>\n\n\n\n

Introduction<\/h2>\n\n\n\n

When we think about the shortest path between two points, we usually imagine a straight line (a Euclidean distance). But once those points lie on a curved surface, the notion of a straight line becomes more subtle. The shortest path must follow the geometry of the surface itself. Such paths are called geodesics<\/strong>.<\/p>\n\n\n\n

Geodesics are fundamental to understanding how geometry shapes motion and distance. Given a surface, we can ask where its geodesics go, how far apart points are, and how the geometry influences the paths between them. This is the familiar forward<\/em> perspective: given the geometry, determine the resulting geodesic structure<\/strong>.<\/p>\n\n\n\n

But what if we reverse the question?<\/p>\n\n\n\n

Instead of starting with a surface and asking what paths it produces, we can start with a desired path or geodesic behavior and ask what geometry would produce it<\/strong>. This shift, from computing geodesics on a given surface to designing the surface around desired geodesics, is the central idea of inverse geodesic design<\/strong>.<\/p>\n\n\n\n

At its core, inverse geodesic design is therefore a problem of shaping geometry through the behavior it induces. By modifying the surface, we modify its intrinsic metric, which in turn changes its geodesic distances and paths. The challenge is to find a geometry whose geodesic structure matches a desired objective.<\/p>\n\n\n\n

This simple reversal opens up an interesting computational problem: can we design a surface by optimizing it directly through its geodesics?<\/strong><\/p>\n\n\n\n

Geodesics with the Heat Method<\/h2>\n\n\n\n

To design a surface around its geodesics, we first need to measure geodesic distance on a mesh, and measure it in a way we can differentiate later on. We use the Heat Method<\/strong> for this (Crane, Weischedel and Wardetzky, 2013).<\/p>\n\n\n\n

The idea goes back to an observation by Varadhan. Put a spot of heat at a point, let it diffuse for a very short time, and the way it spreads already tells you how far away everything is. Nearby points warm up fast, far points stay cold. The Heat Method takes this and turns it into three steps, and each one is just a sparse linear solve. That last part is what we care about, since it means we can differentiate through the whole thing.<\/p>\n\n\n\n

Step 1: Diffuse heat from the source<\/strong><\/h4>\n\n\n\n

We place a unit heat source on the source set \ud835\udcae<\/mi>\\mathcal S<\/annotation><\/semantics><\/math> and let it spread for a short time t<\/mi>t<\/annotation><\/semantics><\/math>. In discrete form this is a single backward-Euler step of the heat equation:<\/p>\n\n\n\n
(<\/mo>M<\/mi>\u2212<\/mo>t<\/mi><\/mspace>L<\/mi>C<\/mi><\/msub>)<\/mo><\/mrow><\/mspace>u<\/mi>=<\/mo>\u03b4<\/mi>\ud835\udcae<\/mi><\/msub>,<\/mo><\/mrow>\\left(M – t\\,L_C\\right)\\,u = \\delta_{\\mathcal S},<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

where L<\/mi>C<\/mi><\/msub>L_C<\/annotation><\/semantics><\/math> is the cotangent Laplacian, M<\/mi>M<\/annotation><\/semantics><\/math> is the (lumped) mass matrix, and \u03b4<\/mi>\ud835\udcae<\/mi><\/msub>\\delta_{\\mathcal S}<\/annotation><\/semantics><\/math> marks the source vertices. Solving this gives u<\/mi>u<\/annotation><\/semantics><\/math>, a smooth blob of heat sitting around the source.<\/p>\n\n\n\n

Step 2: Normalize the heat gradient<\/strong><\/h4>\n\n\n\n

The gradient of the heat field points away from the source, which is the direction distance grows in. Its length is not useful to us, since the heat fades as you move away, but its direction is exactly what we want, so we normalize it to a unit vector field:<\/p>\n\n\n\n

X<\/mi>=<\/mo>\u2212<\/mo>\u2207<\/mo>u<\/mi><\/mrow>\u2016<\/mo>\u2207<\/mo>u<\/mi>\u2016<\/mo><\/mrow><\/mfrac>.<\/mi><\/mrow>X = -\\frac{\\nabla u}{\\lVert \\nabla u \\rVert}.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

The minus sign makes X<\/mi>X<\/annotation><\/semantics><\/math> point away from the source. This gives us a good guess for which way the real distance function should be increasing at each point.<\/p>\n\n\n\n

Step 3: Recover distance with a Poisson solve<\/strong><\/h4>\n\n\n\n

The last step is to find a scalar field \u03d5<\/mi>\\phi<\/annotation><\/semantics><\/math> whose gradient lines up with X<\/mi>X<\/annotation><\/semantics><\/math> as closely as it can. That is a Poisson equation:<\/p>\n\n\n\n
L<\/mi>C<\/mi><\/msub><\/mspace>\u03d5<\/mi>=<\/mo>\u2207<\/mo>\u22c5<\/mo>X<\/mi>,<\/mo><\/mrow>L_C\\,\\phi = \\nabla\\cdot X,<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

and then we shift \u03d5<\/mi>\\phi<\/annotation><\/semantics><\/math> so it reads zero at the source. What comes out is the geodesic distance, up to the usual approximation, from \ud835\udcae<\/mi>\\mathcal S<\/annotation><\/semantics><\/math> to every vertex.<\/p>\n\n\n\n

The time t<\/mi>t<\/annotation><\/semantics><\/math> sets how much the field gets smoothed. We use t<\/mi>=<\/mo>c<\/mi>\u00b7<\/mo>h<\/mi>2<\/mn><\/msup><\/mrow>t = c\\cdot h^2<\/annotation><\/semantics><\/math> with h<\/mi>h<\/annotation><\/semantics><\/math> the mean edge length and c<\/mi>=<\/mo>10<\/mn><\/mrow>c = 10<\/annotation><\/semantics><\/math> across all our grids. Since the three steps are all sparse solves, \u03d5<\/mi>\\phi<\/annotation><\/semantics><\/math> comes out as a differentiable function of the vertex positions. In iskra this is one call, heat_method_distance(V, F, src)<\/code>, and we drop it straight into the optimization loop so the gradients can flow back through it.<\/p>\n\n\n\n

This is the Heat Method on a real 3D mesh. We put one source on the bunny’s ear and colour each vertex by its geodesic distance from it, so every band is a ring of equal distance. The bands wrap around the ears and body instead of cutting straight across, which is the surface geometry bending the distances. Later we deform a surface on purpose to control exactly this.<\/p>\n\n\n\n

\"\"
Geodesic distance from a single source (pink dot) on the Stanford bunny. Colour runs from near (blue) to far (red), and each band is a level set of equal geodesic distance.<\/figcaption><\/figure>\n\n\n\n

Building Geometric Intuition<\/h2>\n\n\n\n

Formulas are easier to believe once you can poke at them and watch what happens. So we built a small interactive demo in Polyscope for the forward problem, the one we have been describing so far.<\/p>\n\n\n\n

Forward: distance follows the geometry<\/strong><\/h4>\n\n\n\n

In the forward demo you click a source vertex and the distance field and its isolines redraw live. On a flat sheet the contours come out as neat circles. Bend the surface into hills and valleys and they start to stretch and bunch up, since the shortest path now has to climb over or go around whatever is in the way. This is the idea the whole project sits on: the shape of the surface is what sets the distances.<\/p>\n\n\n\n

You drop a source and the distance field shows up right away. We start with one source, then add a second, and the field redraws. Each point just takes the distance to the nearer source, so a line forms where the two regions meet.<\/p>\n\n\n\n

\"\"
Forward demo: geodesic distance from a source on the bunny (colour + isolines). A second source is added partway through and the field re-renders; each point takes the distance to whichever source is nearer.<\/figcaption><\/figure>\n\n\n\n

Inverse Geodesics<\/h2>\n\n\n\n

So far we have only read distances off a fixed surface. The inverse problem turns that around. We decide where the source and the target points sit and what distance we want between them, and then we let the surface itself change until that is true. Nothing about the layout moves, only the heights, so the terrain rises and dips until the geodesic distance to the target lands where we asked. The rest of this section makes that concrete: how we set up the terrain, what we minimize, and how we keep the surface from turning into noise.<\/p>\n\n\n\n

Problem Formulation<\/strong><\/h4>\n\n\n\n

For our experiments, we keep the setup clean. We represent the terrain as a triangular heightfield over the x<\/mi>\u2212<\/mo>z<\/mi><\/mrow>x-z<\/annotation><\/semantics><\/math> plane. The planar coordinates and mesh connectivity remain fixed, while the vertex heights (y – coordinates) are optimized.<\/p>\n\n\n\n

Given a source set \ud835\udcae<\/mi>\\mathcal S<\/annotation><\/semantics><\/math> and target set \ud835\udcaf<\/mi>\\mathcal T<\/annotation><\/semantics><\/math>, our goal is to deform the surface until every vertex in \ud835\udcaf<\/mi>\\mathcal T<\/annotation><\/semantics><\/math> has approximately the same geodesic distance from \ud835\udcae<\/mi>\\mathcal S<\/annotation><\/semantics><\/math>.<\/p>\n\n\n\n

For the examples here, we pick a grid diagonal \ud835\udcaf<\/mi>=<\/mo>{<\/mo>(<\/mo>i<\/mi>,<\/mo>j<\/mi>)<\/mo>|<\/mo>i<\/mi>+<\/mo>j<\/mi>=<\/mo>k<\/mi>}<\/mo>,<\/mo><\/mrow>\\mathcal T=\\{(i,j)\\mid i+j=k\\},<\/annotation><\/semantics><\/math> as our target, though in practice \ud835\udcaf<\/mi>\\mathcal{T}<\/annotation><\/semantics><\/math> can be any arbitrary subsets of vertices you want to align.<\/p>\n\n\n\n

Distance Objective and Optimization Workflow<\/strong><\/h4>\n\n\n\n

We measure the geodesic distance \u03d5<\/mi>i<\/mi><\/msub>\\phi_i<\/annotation><\/semantics><\/math> at every vertex using the differentiable Heat Method. Let \ud835\udcaf<\/mi>\\mathcal T<\/annotation><\/semantics><\/math> be the target vertex set and let d<\/mi><\/mtext>target<\/mi><\/mrow><\/msub>d_{\\mathrm{target}}<\/annotation><\/semantics><\/math> be the desired geodesic distance from the source set. We define the distance loss as<\/p>\n\n\n\n
L<\/mi>distance<\/mtext><\/msub>=<\/mo>1<\/mn>|<\/mi>\ud835\udcaf<\/mi>|<\/mi><\/mrow><\/mfrac>\u2211<\/mo>i<\/mi>\u2208<\/mo>\ud835\udcaf<\/mi><\/mrow><\/munder><\/mrow>(<\/mo>\u03d5<\/mi>i<\/mi><\/msub>\u2212<\/mo>d<\/mi><\/mtext>target<\/mi><\/mrow><\/msub>)<\/mo><\/mrow>2<\/mn><\/msup><\/mrow>L_\\text{distance} = \\frac{1}{|\\mathcal T|}\\sum_{i\\in\\mathcal T}\\left(\\phi_i-d_{\\mathrm{target}}\\right)^2<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

In plain terms, this asks every target vertex to sit at the same chosen distance from the source, and it grows whenever one of them drifts off that distance.<\/p>\n\n\n\n

We optimize the mean squared distance error L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>L_{\\mathrm{distance}}<\/annotation><\/semantics><\/math> and report its square root,<\/p>\n\n\n\n
<\/mtext>RMSE<\/mi><\/mrow>=<\/mo>L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub><\/msqrt>,<\/mo><\/mrow>\\mathrm{RMSE} = \\sqrt{L_{\\mathrm{distance}}},<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

which has the same units as geodesic distance. <\/p>\n\n\n\n

We optimize the vertex heights using Adam over a fixed number of iterations. Each iteration follows this pipeline:<\/p>\n\n\n\n

Vertex heights
\u2193
Construct the heightfield mesh
\u2193
Compute Heat Method distances
\u2193
Evaluate L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>L_{\\mathrm{distance}}<\/annotation><\/semantics><\/math>
\u2193
Compute height gradients with backpropagation
\u2193
Update vertex heights with Adam
\u2193
Repeat until the final iteration<\/p>\n\n\n\n

Since the Heat Method implementation is differentiable, the gradient of L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>L_{\\mathrm{distance}}<\/annotation><\/semantics><\/math> can be computed with respect to the vertex heights. We initialize the surface with a small sinusoidal perturbation (amplitude 10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup>10^{-3}<\/annotation><\/semantics><\/math>) to break initial planarity without introducing high-frequency noise.<\/p>\n\n\n\n

Baseline: Distance Loss Only<\/h4>\n\n\n\n

We first optimized the heightfield using only the distance objective, that is, with the smoothness weight set to zero (we introduce that term in the next section):<\/p>\n\n\n\n

L<\/mi>=<\/mo>L<\/mi>distance<\/mtext><\/msub><\/mrow>L = L_\\text{distance}<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

The optimization reduces the L<\/mi>distance<\/mtext><\/msub>L_\\text{distance}<\/annotation><\/semantics><\/math> on \ud835\udcaf<\/mi>\\mathcal{T}<\/annotation><\/semantics><\/math>, making the target diagonal approximately follow a geodesic-distance contour.<\/p>\n\n\n\n
\n
\n
\"\"
(2D) Heat Method distance contours projected onto the original x -z domain.
The red points are the target vertices, and the yellow curve is the desired distance contour (\u03d5<\/mi>=<\/mo>0.4<\/mn><\/mrow>\\phi = 0.4<\/annotation><\/semantics><\/math>). <\/figcaption><\/figure>\n<\/div>\n\n\n\n
\n
\"\"
(Side view) Heightfield optimized using only L<\/mi>distance<\/mtext><\/msub>L_\\text{distance}<\/annotation><\/semantics><\/math>. <\/figcaption><\/figure>\n<\/div>\n<\/div>\n\n\n\n

While minimizing L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>L_{\\mathrm{distance}}<\/annotation><\/semantics><\/math> successfully aligns the target distance, it places no constraint on local surface geometry. Therefore, a small distance loss can still produce high-frequency surface noise and localized spiky oscillations.<\/p>\n\n\n\n

Regularizing the Surface<\/h4>\n\n\n\n

To enforce smoothness on the optimized surface, we consider two regularization schemes: height smoothness<\/strong> and normal smoothness<\/strong>.<\/p>\n\n\n\n

We first penalize height differences between neighboring vertices. Let E<\/mi>E<\/annotation><\/semantics><\/math> be the set of horizontal and vertical grid edges, and let y<\/mi>p<\/mi><\/msub>y_p<\/annotation><\/semantics><\/math> be the height of vertex p<\/mi>p<\/annotation><\/semantics><\/math>. We define the height smoothness<\/strong> as,<\/p>\n\n\n\n
L<\/mi><\/mtext>height<\/mi><\/mrow><\/msub>=<\/mo>1<\/mn>|<\/mi>E<\/mi>|<\/mi><\/mrow><\/mfrac>\u2211<\/mo>(<\/mo>p<\/mi>,<\/mo>q<\/mi>)<\/mo>\u2208<\/mo>E<\/mi><\/mrow><\/munder><\/mrow>(<\/mo>y<\/mi>p<\/mi><\/msub>\u2212<\/mo>y<\/mi>q<\/mi><\/msub>)<\/mo>2<\/mn><\/msup>.<\/mi><\/mrow>L_{\\mathrm{height}} =\\frac{1}{|E|}\\sum_{(p,q)\\in E}(y_p-y_q)^2.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

The loss is small when neighboring vertices have similar heights and large when their height difference increases. Then, the complete objective is<\/p>\n\n\n\n

L<\/mi>=<\/mo>L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>+<\/mo>w<\/mi>smoothness<\/mtext><\/msub>L<\/mi><\/mtext>height<\/mi><\/mrow><\/msub>.<\/mi><\/mrow>L =L_{\\mathrm{distance}} + w_\\text{smoothness} L_{\\mathrm{height}}.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

Height smoothness directly controls neighboring height differences. However, it does not directly measure changes in triangle orientation. A surface can therefore have relatively small height differences while still containing visible local bends.<\/p>\n\n\n\n

To fix this, we switch to Normal Smoothness<\/strong>, which penalizes orientation changes between adjacent face normals rather than absolute height differences. For each adjacent face pair (<\/mo>f<\/mi>i<\/mi><\/msub>,<\/mo>f<\/mi>j<\/mi><\/msub>)<\/mo>\u2208<\/mo>\ud835\udc9c<\/mi><\/mrow>(f_i,f_j)\\in\\mathcal A<\/annotation><\/semantics><\/math>, let \ud835\udc27<\/mi>i<\/mi><\/msub>\\mathbf n_i<\/annotation><\/semantics><\/math> and \ud835\udc27<\/mi>j<\/mi><\/msub>\\mathbf n_j<\/annotation><\/semantics><\/math> be their unit normals. <\/p>\n\n\n\n
L<\/mi><\/mtext>normal<\/mi><\/mrow><\/msub>=<\/mo>1<\/mn>|<\/mi>\ud835\udc9c<\/mi>|<\/mi><\/mrow><\/mfrac>\u2211<\/mo>(<\/mo>f<\/mi>i<\/mi><\/msub>,<\/mo>f<\/mi>j<\/mi><\/msub>)<\/mo>\u2208<\/mo>\ud835\udc9c<\/mi><\/mrow><\/munder><\/mrow>\u2016<\/mo>\ud835\udc27<\/mi>i<\/mi><\/msub>\u2212<\/mo>\ud835\udc27<\/mi>j<\/mi><\/msub>\u2016<\/mo><\/mrow>2<\/mn>2<\/mn><\/msubsup>.<\/mi><\/mrow>L_{\\mathrm{normal}}\n=\n\\frac{1}{|\\mathcal A|}\n\\sum_{(f_i,f_j)\\in\\mathcal A}\n\\left\\|\n\\mathbf n_i-\\mathbf n_j\n\\right\\|_2^2.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

The final objective is<\/p>\n\n\n\n

L<\/mi>=<\/mo>L<\/mi><\/mtext>distance<\/mi><\/mrow><\/msub>+<\/mo>w<\/mi>s<\/mi>m<\/mi>o<\/mi>o<\/mi>t<\/mi>h<\/mi>n<\/mi>e<\/mi>s<\/mi>s<\/mi><\/mrow><\/msub>L<\/mi><\/mtext>normal<\/mi><\/mrow><\/msub>.<\/mi><\/mrow>L = L_{\\mathrm{distance}} + w_{smoothness} L_{\\mathrm{normal}}.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

Pareto Analysis for Smoothness-Weight Selection<\/h4>\n\n\n\n

The regularization weight w<\/mi><\/mtext>smoothness<\/mi><\/mrow><\/msub>w_{\\mathrm{smoothness}}<\/annotation><\/semantics><\/math> dictates the trade-off between geodesic target fidelity and surface smoothness. To select more reasonable optimal hyperparameter values, we sweep w<\/mi>smoothness<\/mtext><\/msub>w_\\text{smoothness}<\/annotation><\/semantics><\/math> and construct Pareto frontiers plotting RMSE against L<\/mi><\/mtext>height<\/mi><\/mrow><\/msub>L_{\\mathrm{height}}<\/annotation><\/semantics><\/math> and L<\/mi><\/mtext>normal<\/mi><\/mrow><\/msub>L_{\\mathrm{normal}}<\/annotation><\/semantics><\/math>, respectively. <\/p>\n\n\n\n

For height regularization, each point is plotted as<\/p>\n\n\n\n

(<\/mo>L<\/mi><\/mtext>height<\/mi><\/mrow><\/msub>,<\/mo> <\/mtext><\/mtext>RMSE<\/mi><\/mrow>)<\/mo><\/mrow>,<\/mo><\/mrow>\\left(L_{\\mathrm{height}},\\ \\mathrm{RMSE}\\right),<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

and for normal regularization as<\/p>\n\n\n\n

(<\/mo>L<\/mi><\/mtext>normal<\/mi><\/mrow><\/msub>,<\/mo> <\/mtext><\/mtext>RMSE<\/mi><\/mrow>)<\/mo><\/mrow>.<\/mi><\/mrow>\\left(L_{\\mathrm{normal}},\\ \\mathrm{RMSE}\\right).<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

Both coordinates are minimized. A point is Pareto-optimal if no other run achieves both lower RMSE and lower smoothness energy.<\/p>\n\n\n\n

\"\"
Accuracy\u2013smoothness trade-off for height regularization. <\/em><\/figcaption><\/figure>\n\n\n\n
\"\"
Accuracy\u2013smoothness trade-off for normal regularization.<\/em><\/figcaption><\/figure>\n\n\n\n

These plots are used independently because the two smoothness energies have different meanings and scales. For each regularizer, we choose a weight near the knee of its Pareto curve: the point beyond which further smoothing causes a much larger increase in distance error. <\/p>\n\n\n\n

Visual Comparison<\/h5>\n\n\n\n

Based on the Pareto curves, we select one representative weight for each regularizer and rerun the optimization to visualize the resulting terrains. Both experiments use the same source, target vertices, desired distance, and initial heightfield.<\/p>\n\n\n\n

    \n
  • Set w<\/mi>height<\/mtext><\/msub>=<\/mo>0.003<\/mn><\/mrow>w_\\text{height} = 0.003<\/annotation><\/semantics><\/math> for height regularization<\/li>\n<\/ul>\n\n\n\n
    \n
    \n
    \"\"
    (2D) Heat Method distance contours projected onto the original x – z domain under height smoothness.
    <\/figcaption><\/figure>\n<\/div>\n\n\n\n
    \n
    \"\"
    (Side view) Heightfield optimized with height smoothness. <\/figcaption><\/figure>\n\n\n\n

    <\/p>\n<\/div>\n<\/div>\n\n\n\n

    It produces several sharp folds and zigzag-like oscillations near the target set. It controls neighboring height differences but does not directly constrain changes in face orientation.<\/p>\n\n\n\n

      \n
    • Set w<\/mi>normal<\/mtext><\/msub>=<\/mo>0.01<\/mn><\/mrow>w_\\text{normal} = 0.01<\/annotation><\/semantics><\/math> for normal regularization<\/li>\n<\/ul>\n\n\n\n
      \n
      \n
      \"\"
      (2D) Heat Method distance contours projected onto the original x – z domain under normal smoothness. <\/figcaption><\/figure>\n<\/div>\n\n\n\n
      \n
      \"\"
      (Side view) Heightfield optimized with normal smoothness.
      <\/figcaption><\/figure>\n<\/div>\n<\/div>\n\n\n\n

      Normal regularization produces a smoother and more coherent terrain. The deformation is distributed gradually around the target set, forming a continuous ridge without the sharp folds and zigzag-like oscillations observed under height regularization. This reflects the direct penalty on orientation changes between adjacent faces.<\/p>\n\n\n\n

      The corresponding measurements are summarized below.<\/p>\n\n\n\n

      Metric <\/td>Height regularization<\/td>Normal regularization<\/td><\/tr>
      Target-distance RMSE<\/td>6.52<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>4<\/mn><\/mrow><\/msup><\/mrow>6.52 \\times10^{-4}<\/annotation><\/semantics><\/math><\/td>3.51<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>3.51\\times10^{-3}<\/annotation><\/semantics><\/math><\/td><\/tr>
      Maximum target error<\/td>1.06<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>1.06 \\times 10^{-3}<\/annotation><\/semantics><\/math><\/td>6.35<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>6.35\\times10^{-3}<\/annotation><\/semantics><\/math><\/td><\/tr>
      Target-distance range<\/td>1.92<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/mrow>1.92 \\times 10{-3}<\/annotation><\/semantics><\/math><\/td>8.84<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>8.84\\times10^{-3}<\/annotation><\/semantics><\/math><\/td><\/tr>
      Height smoothness<\/td>1.70<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>1.70\\times10^{-3}<\/annotation><\/semantics><\/math><\/td>3.20<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>4<\/mn><\/mrow><\/msup><\/mrow>3.20\\times10^{-4}<\/annotation><\/semantics><\/math><\/td><\/tr>
      Normal smoothness<\/td>3.58<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>1<\/mn><\/mrow><\/msup><\/mrow>3.58 \\times 10^{-1}<\/annotation><\/semantics><\/math><\/td>3.07<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>2<\/mn><\/mrow><\/msup><\/mrow>3.07\\times10^{-2}<\/annotation><\/semantics><\/math><\/td><\/tr>
      Height range<\/td>1.78<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>1<\/mn><\/mrow><\/msup><\/mrow>1.78 \\times 10^{-1}<\/annotation><\/semantics><\/math><\/td>1.44<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>1<\/mn><\/mrow><\/msup><\/mrow>1.44\\times10^{-1}<\/annotation><\/semantics><\/math><\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

      Remark: Our main goal is to remove noisy folds and produce a smoother terrain. Normal regularization achieves this more effectively by directly penalizing orientation changes between adjacent faces. We therefore use normal smoothness in the following experiments.<\/p>\n\n\n\n

      Forward Solver Comparison: Convex RGD vs. Heat Method<\/h4>\n\n\n\n

      iskra’s original code computes the forward distance field with a convex-optimization solver, following the regularized geodesic distance (RGD) formulation of Edelstein et al., solved with an ADMM scheme. We compare that convex solver against the Heat Method while keeping the inverse optimization setup fixed.<\/p>\n\n\n\n

      Both methods use the same grid, source, target vertices, desired distance, stopping tolerance, and maximum number of iterations. Each experiment is repeated five times.<\/p>\n\n\n\n

      Forward solver<\/th>Converged runs<\/th>Iterations<\/th>Median time<\/th>Final L2 error<\/th>Final RMSE<\/th>Max error<\/th><\/tr>
      Convex RGD (Edelstein et al.)<\/td>0\/5<\/td>50<\/td>59.41 s<\/td>0.1307<\/td>0.0394<\/td>0.0694<\/td><\/tr>
      Heat Method<\/td>5\/5<\/td>29<\/td>3.97 s<\/td>0.00948<\/td>0.00286<\/td>0.00575<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

      Across all test runs, the Heat Method converges in all five runs and is approximately 15 times faster than the convex RGD solver. It also reaches a substantially lower distance error. We therefore use the Heat Method as the forward solver in the following experiments.<\/p>\n\n\n\n

      Before moving on, here is the whole thing running. We fix a target outline and a distance, and the flat sheet bends frame by frame until the geodesic contour traces the shape we asked for. The source ends up on top and the contours run downhill into the target.<\/p>\n\n\n\n

      \"\"
      Inverse design in action: the heightfield bends until the geodesic distance contour from the central source traces a square target. The source sits on the summit and the contours ripple downhill into a square.<\/figcaption><\/figure>\n\n\n\n

      Going Neural: Terrain as an MLP<\/h2>\n\n\n\n

      Up to this point the unknowns were the vertex heights themselves, one number per vertex. In week 6 we swapped that out. Instead of keeping a separate height at every vertex, we let a small neural network hand us the height as a continuous function of position:<\/p>\n\n\n\n

      y<\/mi>i<\/mi><\/msub>=<\/mo>f<\/mi>\u03b8<\/mi><\/msub>(<\/mo>x<\/mi>i<\/mi><\/msub>,<\/mo>z<\/mi>i<\/mi><\/msub>)<\/mo>,<\/mo><\/mrow>y_i = f_\\theta(x_i, z_i),<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      and we optimize the network weights \u03b8<\/mi>\\theta<\/annotation><\/semantics><\/math> rather than the heights. The rest of the loop does not change at all. We build the mesh from these heights, run the Heat Method, measure the distance loss, and take an optimizer step. The one thing that is different is where the heights come from.<\/p>\n\n\n\n

      There are a few reasons to do this. The network is defined everywhere on the domain, not only at the grid points, so it does not really care about the mesh resolution and we can train on a coarse grid and then read it off on a finer one. It also comes out smooth on its own, and the parameter count stays fixed no matter how many vertices we throw at it.<\/p>\n\n\n\n

      Implementation<\/strong><\/h4>\n\n\n\n

      The height field f<\/mi>\u03b8<\/mi><\/msub>f_\\theta<\/annotation><\/semantics><\/math> is a small multilayer perceptron. It takes the coordinate pair (<\/mo>x<\/mi>,<\/mo>z<\/mi>)<\/mo><\/mrow>(x,z)<\/annotation><\/semantics><\/math>, normalized to [<\/mo>\u2212<\/mo>1<\/mn>,<\/mo>1<\/mn>]<\/mo><\/mrow>[-1,1]<\/annotation><\/semantics><\/math>, runs it through two hidden layers of width 64 with tanh<\/mi><\/mrow>\\tanh<\/annotation><\/semantics><\/math> activations, and returns a single number. We squash that output so the heights cannot run off to anything extreme:<\/p>\n\n\n\n
      y<\/mi>=<\/mo>scale<\/mtext>\u22c5<\/mo>tanh<\/mi>\u2061<\/mo><\/mrow>(<\/mo>MLP<\/mtext>(<\/mo>x<\/mi>,<\/mo>z<\/mi>)<\/mo>)<\/mo>.<\/mi><\/mrow>y = \\text{scale}\\cdot\\tanh\\big(\\text{MLP}(x,z)\\big).<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      Here scale<\/mtext>\\text{scale}<\/annotation><\/semantics><\/math> is a height cap, setting how far up or down the terrain is allowed to go. The whole thing still sits on top of iskra’s differentiable Heat Method, so the gradient from the distance loss makes it all the way back to the weights \u03b8<\/mi>\\theta<\/annotation><\/semantics><\/math>. The loop now looks like this:<\/p>\n\n\n\n

      Network weights \u03b8<\/strong>
      \u2193
      Heights yi<\/sub><\/em> = f\u03b8<\/sub><\/em>(xi<\/sub><\/em>, zi<\/sub><\/em>)
      \u2193
      Construct the heightfield mesh
      \u2193
      Compute Heat Method distances
      \u2193
      Evaluate L<\/em>distance<\/sub>
      \u2193
      Update \u03b8<\/em> with Adam<\/p>\n\n\n\n

      Preliminary visuals based on MLP<\/h4>\n\n\n\n

      To validate implicit neural representations for terrain design, we fit a coordinate-based MLP on a coarse 32<\/mn>\u00d7<\/mo>32<\/mn><\/mrow>32\\times32<\/annotation><\/semantics><\/math> grid. The network converges after 378<\/mn>378 <\/annotation><\/semantics><\/math>iterations in 29.8<\/mn>29.8 <\/annotation><\/semantics><\/math> seconds. The final target-distance RMSE is 1.00<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>1.00\\times10^{-3}<\/annotation><\/semantics><\/math>, with a maximum target error of 1.68<\/mn>\u00d7<\/mo>10<\/mn>\u2212<\/mo>3<\/mn><\/mrow><\/msup><\/mrow>1.68\\times10^{-3}<\/annotation><\/semantics><\/math>.<\/p>\n\n\n\n
      \"\"
      Training resolution <\/strong>32<\/mn>\u00d7<\/mo>32<\/mn><\/mrow>32\\times32<\/annotation><\/semantics><\/math>.<\/strong> <\/figcaption><\/figure>\n\n\n\n

      After training, we evaluate the same MLP on a finer 128<\/mn>\u00d7<\/mo>128<\/mn><\/mrow>128\\times128<\/annotation><\/semantics><\/math> grid without further optimization. This coarse-to-fine evaluation produces a denser visualization of the learned continuous heightfield while preserving the height range [<\/mo>\u2212<\/mo>0.2,0.2<\/mn>]<\/mo><\/mrow>[-0.2,0.2]<\/annotation><\/semantics><\/math>.<\/p>\n\n\n\n
      \"\"
      Fine evaluation <\/strong>128<\/mn>\u00d7<\/mo>128<\/mn><\/mrow>128\\times128<\/annotation><\/semantics><\/math>. <\/figcaption><\/figure>\n\n\n\n

      The finer sampling shows that the learned terrain remains continuous across resolutions and that the target set closely follows the desired distance contour.<\/p>\n\n\n\n

      Effect of Output Scale and Activation Function<\/h4>\n\n\n\n

      The MLP output scale limits how much the terrain can deform. When the scale is too small, the network cannot produce enough height variation to match the target distance.<\/p>\n\n\n\n

      Scale<\/th>Converged<\/th>Iterations<\/th>RMSE<\/th>Normal smoothness<\/th>Height range<\/th>Runtime<\/th><\/tr>
      0.05<\/td>No<\/td>400<\/td>0.2935<\/td>0<\/td>0<\/td>30.54 s<\/td><\/tr>
      0.10<\/td>No<\/td>400<\/td>0.1432<\/td>0.00689<\/td>0.200<\/td>30.17 s<\/td><\/tr>
      0.15<\/td>No<\/td>400<\/td>0.0678<\/td>0.00866<\/td>0.300<\/td>30.91 s<\/td><\/tr>
      0.20<\/strong><\/td>Yes<\/td>27<\/td>0.0317<\/td>0.00241<\/td>0.400<\/td>2.13 s<\/td><\/tr>
      0.25<\/td>Yes<\/td>21<\/td>0.0303<\/td>0.00092<\/td>0.443<\/td>1.66 s<\/td><\/tr>
      0.30<\/td>Yes<\/td>18<\/td>0.0289<\/td>0.00079<\/td>0.505<\/td>1.40 s<\/td><\/tr>
      0.40<\/td>Yes<\/td>14<\/td>0.0315<\/td>0.00063<\/td>0.543<\/td>1.10 s<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

      The optimization fails to converge for scales below 0.2<\/mn>0.2<\/annotation><\/semantics><\/math>. At these scales, the learned height range reaches the available output range, indicating that the MLP does not have enough vertical freedom. For the following experiments, we use scale = 0.2.<\/mn>0.2.<\/annotation><\/semantics><\/math><\/p>\n\n\n\n

      We also compare several activation functions across five random seeds. The seed controls the initial values of the network parameters. Repeating each experiment with multiple seeds therefore measures the stability of each activation rather than relying on a single favorable run.<\/p>\n\n\n\n

      Activation<\/td>Convergence rate<\/td>Median iterations<\/td>Median RMSE<\/td>Median normal smoothness<\/td>Median runtime<\/td><\/tr>
      Tanh<\/td>0.80<\/td>30<\/td>0.0369<\/td>0.00255<\/td>2.39 s<\/td><\/tr>
      ReLU<\/td>1.00<\/td>30<\/td>0.0307<\/td>0.00311<\/td>2.33 s<\/td><\/tr>
      Softplus<\/td>1.00<\/td>56<\/td>0.0294<\/td>0.00350<\/td>4.37 s<\/td><\/tr>
      SiLU<\/td>1.00<\/td>39<\/td>0.0336<\/td>0.00235<\/td>3.05 s<\/td><\/tr><\/tbody><\/table><\/figure>\n\n\n\n

      No single activation function performs best across all metrics. ReLU provides the lowest median runtime, Softplus achieves the lowest median RMSE, and SiLU produces the lowest median normal-smoothness energy. Tanh is less robust, with one run failing to converge.<\/p>\n\n\n\n

      Consequently, this is a full multi-objective comparison that could use Pareto analysis to balance accuracy, smoothness, runtime, and convergence stability. We leave this analysis for future work and use Tanh<\/strong> as a fixed activation function in the following experiments<\/p>\n\n\n\n

      Fourier Analysis of the Plain MLP Heightfield<\/h4>\n\n\n\n

      Before introducing positional encoding, we first examine the frequencies represented by the plain coordinate-based MLP. Since the learned terrain is a continuous function y<\/mi>=<\/mo>f<\/mi>\u03b8<\/mi><\/msub>(<\/mo>x<\/mi>,<\/mo>z<\/mi>)<\/mo><\/mrow>y=f_\\theta(x,z)<\/annotation><\/semantics><\/math>, we can decompose its sampled heightfield into a spatial frequency components. <\/p>\n\n\n\n

      This analysis shows whether the learned surface is dominated by broad, long-wavelength deformation or also contains shorter-wavelength local detail. It also provides a baseline for evaluating how positional encoding changes the frequency content of the terrain.<\/p>\n\n\n\n

      We sample the trained MLP on an N<\/mi>\u00d7<\/mo>N<\/mi><\/mrow>N\\times N<\/annotation><\/semantics><\/math> grid and denote the resulting heights by y<\/mi>m<\/mi>,<\/mo>n<\/mi><\/mrow><\/msub>y_{m,n}<\/annotation><\/semantics><\/math>. Before computing the Fourier transform, we remove the mean height,<\/p>\n\n\n\n
      y<\/mi>~<\/mo><\/mover>m<\/mi>,<\/mo>n<\/mi><\/mrow><\/msub>=<\/mo>y<\/mi>m<\/mi>,<\/mo>n<\/mi><\/mrow><\/msub>\u2212<\/mo>y<\/mi>\u203e<\/mo><\/mover>,<\/mo><\/mspace>where<\/mtext><\/mspace>y<\/mi>\u203e<\/mo><\/mover>=<\/mo>1<\/mn>N<\/mi>2<\/mn><\/msup><\/mfrac>\u2211<\/mo>m<\/mi>=<\/mo>0<\/mn><\/mrow>N<\/mi>\u2212<\/mo>1<\/mn><\/mrow><\/munderover><\/mrow>\u2211<\/mo>n<\/mi>=<\/mo>0<\/mn><\/mrow>N<\/mi>\u2212<\/mo>1<\/mn><\/mrow><\/munderover><\/mrow>y<\/mi>m<\/mi>,<\/mo>n<\/mi><\/mrow><\/msub>.<\/mi><\/mrow>\\widetilde y_{m,n} = y_{m,n}-\\bar y, \\quad \\text{where} \\quad \\bar{y}\n=\n\\frac{1}{N^2}\n\\sum_{m=0}^{N-1}\n\\sum_{n=0}^{N-1}\ny_{m,n}.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      This removes the constant vertical offset and leaves only spatial variation. The Fourier coefficients are<\/p>\n\n\n\n

      Y<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>=<\/mo>\u2211<\/mo>m<\/mi>=<\/mo>0<\/mn><\/mrow>N<\/mi>\u2212<\/mo>1<\/mn><\/mrow><\/munderover><\/mrow>\u2211<\/mo>n<\/mi>=<\/mo>0<\/mn><\/mrow>N<\/mi>\u2212<\/mo>1<\/mn><\/mrow><\/munderover><\/mrow>y<\/mi>~<\/mo><\/mover>m<\/mi>,<\/mo>n<\/mi><\/mrow><\/msub>exp<\/mi>\u2061<\/mo><\/mrow>[<\/mo>\u2212<\/mo>2<\/mn>\u03c0<\/mi>i<\/mi>(<\/mo>p<\/mi>m<\/mi><\/mrow>N<\/mi><\/mfrac>+<\/mo>q<\/mi>n<\/mi><\/mrow>N<\/mi><\/mfrac>)<\/mo><\/mrow>]<\/mo><\/mrow>.<\/mi><\/mrow>Y_{p,q} = \\sum_{m=0}^{N-1}\\sum_{n=0}^{N-1}\\widetilde y_{m,n}\\exp\\left[-2\\pi i\\left(\\frac{pm}{N}+\\frac{qn}{N}\\right)\\right].<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      with power,<\/p>\n\n\n\n

      P<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>=<\/mo>|<\/mo>Y<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>|<\/mo><\/mrow>2<\/mn><\/msup>.<\/mi><\/mrow>P_{p,q} = \\left|Y_{p,q}\\right|^2.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      After shifting the zero-frequency coefficient to the center, coefficients near the center describe broad terrain deformation, while coefficients farther away describe finer spatial variation. We display log<\/mi>10<\/mn><\/msub>\u2061<\/mo>(<\/mo>P<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>)<\/mo><\/mrow>\\log_{10}(P_{p,q})<\/annotation><\/semantics><\/math>.<\/p>\n\n\n\n

      To obtain the radial power spectrum, we group the coefficients by normalized radial frequency,<\/p>\n\n\n\n

      \u03c1<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>=<\/mo>f<\/mi>p<\/mi>2<\/mn><\/msubsup>+<\/mo>f<\/mi>q<\/mi>2<\/mn><\/msubsup><\/mrow><\/msqrt>max<\/mi>u<\/mi>,<\/mo>v<\/mi><\/mrow><\/munder>\u2061<\/mo><\/mspace>f<\/mi>u<\/mi>2<\/mn><\/msubsup>+<\/mo>f<\/mi>v<\/mi>2<\/mn><\/msubsup><\/mrow><\/msqrt><\/mstyle><\/mfrac>,<\/mo><\/mrow>\\rho_{p,q} = \\frac{\\sqrt{f_p^2+f_q^2}}{\\displaystyle\\max_{u,v}\\sqrt{f_u^2+f_v^2}},<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      where f<\/mi>p<\/mi><\/msub>f_p<\/annotation><\/semantics><\/math> and f<\/mi>q<\/mi><\/msub>f_q<\/annotation><\/semantics><\/math> are the frequencies along the x – and z – directions. We then average the power of coefficients with similar radial frequencies. This removes directional information and shows how the spectral power changes with overall frequency magnitude.<\/p>\n\n\n\n

      To compare broad deformation with finer geometric variation, we divide the normalized frequency range at \u03c1<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>=<\/mo>0.25<\/mn><\/mrow>\\rho_{p,q}=0.25<\/annotation><\/semantics><\/math>:<\/p>\n\n\n\n
      low frequency: <\/mtext>\u03c1<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>\u2264<\/mo>0.25<\/mn>,<\/mo><\/mspace>high frequency: <\/mtext>\u03c1<\/mi>p<\/mi>,<\/mo>q<\/mi><\/mrow><\/msub>><\/mo>0.25.<\/mn><\/mrow>\\text{low frequency: } \\rho_{p,q}\\leq 0.25,\\qquad\\text{high frequency: } \\rho_{p,q}>0.25.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      This cutoff separates the region near the frequency origin from the higher-frequency portion of the spectrum.<\/p>\n\n\n\n

      \"\"
      Fourier spectrum of the MLP heightfield.<\/strong>
      The 2D log-power spectrum (left); radial power spectrum (right). The dashed line marks the normalized frequency cutoff at (0.25).<\/figcaption><\/figure>\n\n\n\n

      The radial power decreases rapidly with frequency. Nearly all spectral power lies below the 0.25<\/mn>0.25<\/annotation><\/semantics><\/math> cutoff, while the high-frequency contribution rounds to 0.00<\/mn>%<\/mi><\/mrow>0.00\\%<\/annotation><\/semantics><\/math>. This shows that the learned terrain is dominated by large-scale, low-frequency deformation.<\/p>\n\n\n\n

      So what does this tell us for the inverse problem? The plain MLP hits its distance targets almost entirely through broad, large-scale deformation rather than fine local adjustments, which is why the surface stays free of vertex-level noise even while moving large regions of the terrain. The analysis does not change the geometry itself; it pins down the limitation of this representation, that it has essentially no high-frequency content, and that is exactly what motivates positional encoding, which we turn to next to add short-wavelength detail when a target needs it.<\/p>\n\n\n\n

      The Fix: Positional Encoding<\/h2>\n\n\n\n

      The plain MLP has a built-in weakness known as spectral bias<\/strong>. Networks pick up smooth, low-frequency shapes almost immediately, but sharp, high-frequency detail comes very slowly, if it comes at all. So when the target we want has sharp features, like corners or a thin ridge, the plain MLP barely budges. There is almost no gradient pushing it toward that high-frequency shape, and it just stays smooth and flat and never gets near the target.<\/p>\n\n\n\n

      The fix is positional encoding<\/strong>. Instead of handing the network the raw coordinates and waiting for it to build high frequencies by itself (it will not), we give it those frequencies up front. Before the MLP sees anything, we expand each coordinate into a set of sinusoids at rising frequencies:<\/p>\n\n\n\n

      \u03b3<\/mi>(<\/mo>x<\/mi>,<\/mo>z<\/mi>)<\/mo>=<\/mo>[<\/mo><\/mspace>x<\/mi>,<\/mo><\/mspace>z<\/mi>,<\/mo><\/mspace><\/mspace>sin<\/mi>\u2061<\/mo><\/mrow>(<\/mo>\u03c0<\/mi>x<\/mi>)<\/mo>,<\/mo><\/mspace><\/mspace>cos<\/mi>\u2061<\/mo><\/mrow>(<\/mo>\u03c0<\/mi>x<\/mi>)<\/mo>,<\/mo><\/mspace><\/mspace>sin<\/mi>\u2061<\/mo><\/mrow>(<\/mo>2<\/mn>\u03c0<\/mi>x<\/mi>)<\/mo>,<\/mo><\/mspace><\/mspace>cos<\/mi>\u2061<\/mo><\/mrow>(<\/mo>2<\/mn>\u03c0<\/mi>x<\/mi>)<\/mo>,<\/mo><\/mspace>\u2026<\/mo>,<\/mo><\/mspace><\/mspace>sin<\/mi>\u2061<\/mo><\/mrow>(<\/mo>\u03c0<\/mi>z<\/mi>)<\/mo>,<\/mo><\/mspace><\/mspace>cos<\/mi>\u2061<\/mo><\/mrow>(<\/mo>\u03c0<\/mi>z<\/mi>)<\/mo>,<\/mo><\/mspace>\u2026<\/mo><\/mspace>]<\/mo><\/mrow>\\gamma(x,z) = \\big[\\,x,\\; z,\\; \\sin(\\pi x),\\; \\cos(\\pi x),\\; \\sin(2\\pi x),\\; \\cos(2\\pi x),\\; \\ldots,\\; \\sin(\\pi z),\\; \\cos(\\pi z),\\; \\ldots\\,\\big]<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      and feed this encoded vector to the same MLP:<\/p>\n\n\n\n

      y<\/mi>=<\/mo>scale<\/mtext>\u22c5<\/mo>tanh<\/mi>\u2061<\/mo><\/mrow>(<\/mo>MLP<\/mtext>(<\/mo>\u03b3<\/mi>(<\/mo>x<\/mi>,<\/mo>z<\/mi>)<\/mo>)<\/mo>)<\/mo>.<\/mi><\/mrow>y = \\text{scale}\\cdot\\tanh\\big(\\text{MLP}(\\gamma(x,z))\\big).<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      Now that the high frequencies are already sitting in the input, the network does not have to conjure them up, and it can actually draw sharp terrain. The number of bands L<\/mi>L<\/annotation><\/semantics><\/math> is the only real knob here, and more bands buys you finer detail. We build \u03b3<\/mi><\/mrow>\\gamma<\/annotation><\/semantics><\/math> with iskra’s HarmonicEmbedding<\/code>, so none of this adds any extra learnable parameters.<\/p>\n\n\n\n

      Result 1: sharper contours for the same budget<\/strong><\/h4>\n\n\n\n

      With the same height cap and the same training budget, the plain MLP gives us round, washed-out contours, while the positional-encoding one actually lands on the sharp target shape we asked for. The top-down view is where the gap really shows.<\/p>\n\n\n\n

      \"\"
      Plain MLP (round, blurred contours) vs positional-encoding MLP (sharp target contours), oblique and top view.<\/figcaption><\/figure>\n\n\n\n

      Result 2: the target at a smaller height budget<\/strong><\/h4>\n\n\n\n

      We also swept the height cap directly. As we shrink scale<\/mtext>\\text{scale}<\/annotation><\/semantics><\/math>, the terrain has less room to move, so reaching the target distance gets harder. The positional-encoding MLP still converges all the way down to scale 0.1, whereas the plain MLP gives up somewhere below 0.2. So positional encoding gets to the same target on roughly half the height budget.<\/p>\n\n\n\n
      \"\"
      Final RMSE to the target distance vs height cap, plain MLP vs positional-encoding MLP. PE converges down to scale 0.1; the plain MLP needs scale of at least 0.2.<\/figcaption><\/figure>\n\n\n\n

      Future work<\/h2>\n\n\n\n

      What is next: MLP-derivative smoothness<\/strong><\/h4>\n\n\n\n

      Since the terrain is now a proper continuous function, the smoothness question gets cleaner. The mesh-based terms we used earlier, penalizing height differences between neighbors or the bending between adjacent face normals, are really measuring one particular triangulation. Refine the mesh and the same surface gives a different penalty, which always felt slightly off, since what we actually care about is whether the underlying surface is smooth. With the MLP we can ask that directly. Because the height is a differentiable function of position, we hand the coordinates to autograd, get the true gradient of the surface at any point, and penalize its magnitude:<\/p>\n\n\n\n

      L<\/mi><\/mtext>deriv<\/mi><\/mrow><\/msub>=<\/mo>1<\/mn>|<\/mi>V<\/mi>|<\/mi><\/mrow><\/mfrac>\u2211<\/mo>i<\/mi><\/munder><\/mrow>\u2016<\/mo>\u2207<\/mo>f<\/mi>\u03b8<\/mi><\/msub>(<\/mo>x<\/mi>i<\/mi><\/msub>,<\/mo>z<\/mi>i<\/mi><\/msub>)<\/mo>\u2016<\/mo><\/mrow>2<\/mn><\/msup>.<\/mi><\/mrow>L_{\\mathrm{deriv}} = \\frac{1}{|V|}\\sum_{i}\\left\\lVert \\nabla f_\\theta(x_i,z_i)\\right\\rVert^2.<\/annotation><\/semantics><\/math><\/div>\n\n\n\n

      Unlike the mesh-based smoothness terms, this one does not depend on the mesh at all. It reads the smoothness of the surface itself, so training on a coarse grid and evaluating on a fine one gives a consistent notion of smooth either way. The next step is to run the same accuracy-versus-smoothness sweep we did for the mesh terms, but with this derivative penalty, and see where it lands on the trade-off curve. Our hunch is that it gives a gentler, more uniform surface for the same distance error, since it is not chasing per-triangle artifacts.<\/p>\n\n\n\n

      A couple of other directions we want to chase from here:<\/p>\n\n\n\n