The Poisson solve

a side note to Painting → Bas-Relief, on how you get a surface back out of a pile of slopes
Everything here is checkable arithmetic. No calculus is assumed.

You arrived here from §6.2 of the build guide, which says the range compression “reintegrates by solving a Poisson equation with Neumann boundaries” and then moves on. This page is that sentence, unpacked.

0. The problem in one picture

The compression step in §6.2 hands you a field of slopes — for every pixel, how much the height should change going right, and how much going down. What it does not hand you is a surface. Getting from one to the other is the entire job of poisson_from_gradients.

In 1D that is trivial. In 2D it is not, for one specific reason, and the Poisson solve is the fix.

1. One dimension, where it is easy

Take five heights. The slopes between them are just the differences:

heights   0     2     3     1     4
slopes       +2    +1    -2    +3

cumulative sum from 0:
          0     2     3     1     4      <- the heights, exactly

To go back the other way you add the slopes up as you walk along — a cumulative sum. Start at 0, add 2, add 1, subtract 2, add 3, and you are back where you started from.

One thing is permanently lost. The slopes never told you how high the first point was. Start your walk at 100 instead of 0 and every height comes out 100 higher, with exactly the same slopes. This is the + C from school integration, and it never goes away — in the code it is the reason for h_hat[0, 0] = 0.0, which simply picks the version whose average is zero. For a relief it does not matter at all: norm01 renormalises straight afterwards.

2. Two dimensions, where it breaks

Now you have two slope fields: gx (how much height changes going right) and gy (going down). To find the height at some pixel, walk to it from the corner, adding slopes as you go.

But now there is more than one route. And here is the thing that makes 2D genuinely different:

If the slopes came from a real surface, every route gives the same answer. If they did not, different routes disagree — and there is no surface to find.

Take a real 3×3 surface and its true slopes, then scale gx by 0.5 and leave gy alone — which is exactly the kind of thing compress_range does:

surface H          true gx            true gy       gx scaled by 0.5
  0  1  3          1  2  0            1  1 -2          0.5  1.0  0
  1  2  1          1 -1  0            3 -2  1          0.5 -0.5  0
  4  0  2         -4  2  0            0  0  0         -2.0  1.0  0

walk from the top-left corner to the bottom-right corner:
  along the top, then down the right edge   ->  +0.500
  down the left edge, then along the bottom ->  +3.000
                                    disagreement   2.500

Two and a half units of disagreement about the height of a single corner, depending on which way round you walked. No surface in existence has those slopes. And §6.2's whole method is built on producing exactly this kind of impossible field on purpose — attenuating the steep slopes and not the shallow ones is, by construction, an independent rescaling.

No surface has these slopes, so the question becomes:

Which surface comes closest to having these slopes?

3. “Closest”, made precise

Closest in the least-squares sense: of all possible surfaces, find the one where the total squared disagreement between its actual slopes and the slopes you asked for is as small as possible. Add up, over every pixel, (actual gx − wanted gx)² plus (actual gy − wanted gy)², and make that total as small as it will go.

That is a minimisation over a third of a million unknowns, but a standard move makes it tractable: at the minimum of a smooth function, the derivative is zero. Writing “the derivative with respect to every height is zero” and simplifying gives a single equation:

∇2h = ∇ · g

That is the Poisson equation. Both symbols are simpler than they look, and both are just arithmetic on neighbouring pixels.

3.1 The right-hand side: divergence

∇ · g is the divergence of the slope field. Per pixel it is one subtraction per axis: how much gx changed since the pixel on the left, plus how much gy changed since the pixel above.

Read it as a measure of how much the slope field is spreading out at that point. Positive divergence means the arrows around this pixel are pointing away from each other, which is what happens at the bottom of a bowl. Negative means they converge — the top of a hill. It is a single number per pixel saying “this is a dip” or “this is a bump”, and it is computed directly from the slopes you were given, so the right-hand side is known before you start.

3.2 The left-hand side: the Laplacian

∇2h is the Laplacian of the surface. For a grid it is the famous five-point stencil: take a pixel's four neighbours, add them up, subtract four times the pixel itself.

            h[r-1, c]
                |
h[r, c-1] --  h[r, c]  -- h[r, c+1]        laplacian = (sum of the four neighbours)
                |                                      - 4 * h[r, c]
            h[r+1, c]

Read it as how much this pixel differs from the average of its neighbours. Zero means the pixel sits exactly on the average — locally flat, in the sense of no curvature. Positive means it sits below its surroundings, negative means it pokes above them.

So the equation says: find the surface whose local curvature everywhere matches the spreading-out of the slopes you asked for. That is all it says.

4. Boundaries, and why Neumann

The stencil needs four neighbours. On the edge of the image, one of them does not exist. On a corner, two do not. The equation is silent about this, so you have to decide — and the decision has a name.

Dirichlet boundaries fix the value at the edge: “the border is at height zero”. Neumann boundaries fix the slope at the edge: “nothing flows across the border”, which in practice means pretending the missing neighbour is a copy of the pixel itself.

A useful mental picture: a Poisson solve is a soap film settling. Dirichlet is a film on a wire loop — the rim is pinned where the wire is. Neumann is a film whose edge is free to slide up and down a frictionless wall — it can sit at any height, but it must meet the wall at right angles.

For a relief panel, Neumann is obviously right. You do not know in advance how high the edges of your relief should be, and you certainly do not want them forced to zero. Pinning the border would warp the entire field to reach it — on the same test problem, the Neumann solve recovers the true surface to 4 × 10−15 while the Dirichlet solve is off by 4.03 near the edge. The recessed border you do want comes later and separately, from apply_rim, as a deliberate cosmetic choice rather than a constraint baked into the physics.

Neumann has one quirk that explains a line of the code. Since nothing flows across the boundary, the total of the whole right-hand side must come to zero — whatever spreads out somewhere has to converge somewhere else:

f              3.  1.  4.  1.  5.
padded     3.  3.  1.  4.  1.  5.  5.      <- edge value repeated
laplacian     -2.  5. -6.  7. -4.
                                sum = 0.0  <- always, under Neumann

That also means the constant surface satisfies the equation with zero on the right. Which is the same + C from §1, arriving from a different direction: the solution is unique only up to an overall height.

5. Solving it without a solver

Written out, the Poisson equation is one linear equation per pixel — 329,120 of them at the build guide's grid, each coupling a pixel to its four neighbours. As a matrix that is 329,120 × 329,120. Storing it densely would take about 0.9 TB, and solving it the naive way is on the order of 1016 operations.

Nobody does that. The usual approach is an iterative solver — guess, refine, repeat until it stops changing — which works but introduces a tolerance to tune and a convergence to worry about.

Under Neumann boundaries there is something much better available, and it rests on one fact:

Cosine waves are eigenvectors of the Laplacian. Apply the five-point stencil to a cosine and you get the same cosine back, just scaled.

“Eigenvector” only means: this operation does not change my shape, it only changes my size. Check it on a 5×5 grid — the predicted scale factor, and the one measured by actually applying the stencil:

basis fn     predicted lambda     measured        spread
cos(0,0)          0.0000          0.0000       0.0e+00
cos(0,1)         -0.3820         -0.3820       5.0e-16
cos(1,0)         -0.3820         -0.3820       1.7e-16
cos(2,3)         -4.0000         -4.0000       1.3e-15
cos(4,4)         -7.2361         -7.2361       6.2e-15

The scaling factor for the cosine with i wiggles down and j across is

λij = 2(cos(πi/M) − 1) + 2(cos(πj/N) − 1)

and it turns the whole problem inside out. In ordinary space the Laplacian is a huge tangle of coupled equations. In the cosine basis every equation is independent — each cosine component of the answer is just that component of the right-hand side divided by one number. A matrix inversion becomes a division.

That gives the recipe, which is the entire solver:

1.  DCT the right-hand side          -> how much of each cosine it contains
2.  divide each one by its lambda    -> how much of each cosine the answer contains
3.  inverse DCT                      -> the answer, back in pixel space

The DCT is the transform that decomposes an image into exactly these cosines, and it runs in O(n log n) — about 6 × 106 operations at the real grid rather than 1016, and the result is exact to floating point.

5.1 The two fiddly lines

Look at the eigenvalue grid for a 5×5 and one entry stands out:

 0.000  -0.382  -1.382  -2.618  -3.618
-0.382  -0.764  -1.764  -3.000  -4.000
-1.382  -1.764  -2.764  -4.000  -5.000
-2.618  -3.000  -4.000  -5.236  -6.236
-3.618  -4.000  -5.000  -6.236  -7.236

The corner is exactly zero, and it is the only one. That is the constant surface — the + C — showing up as an eigenvalue of zero, because a flat surface has no curvature anywhere. Dividing by it would be dividing by zero.

Hence the two lines in the code that otherwise look like superstition:

    lam[0, 0] = 1.0                                   # make the division legal
    h_hat = dctn(div, type=2, norm="ortho") / lam
    h_hat[0, 0] = 0.0                                 # discard the meaningless result

Set the offending eigenvalue to anything non-zero so the division is legal, then throw away whatever landed in that slot afterwards. The value 1.0 is arbitrary; nothing survives it. What you are choosing is the zero-mean member of the family of valid answers.

6. Does it actually work?

Two checks, both cheap.

When the slopes are consistent, it is exact. Take a real surface, compute its slopes, reintegrate, and you get the surface back to floating-point noise — 1.1 × 10−15 on the 3×3 example above, and 2.7 × 10−13 at 256×256 in §6.5's self-test.

When they are not, it finds the genuine minimum. Feed it the impossible field from §2 and it returns a surface with a squared error of 2.4479 against what was asked for. Nudge that surface 20,000 times at random and see how many nudges do better:

squared error of the Poisson answer:  2.4479
of 20,000 random nudges to it, how many did better:  0

None — it is the provably closest surface to a request that could not be satisfied.

7. The code, mapped

Everything above is fourteen lines in heightfield.py. Each line now has a name:

Try it once

Open a Python prompt with heightfield.py importable and paste the 3×3 example from §2. Print the surface it gives back, take its differences, and compare them with what you asked for. Ten seconds, and the idea stops being abstract — you can see it splitting the difference between two irreconcilable demands.

8. Where else this turns up

Once you recognise the shape of it, gradient-domain editing is everywhere. The pattern is always the same: it is easier to say what the slopes should be than what the values should be, so you specify slopes and solve for values.

References

Back to §6.2, or to the build guide from the top.