Choosing a Plane When the Boundary Won’t #

Least-Squares Elevation Fitting for Non-Planar Platforms #

The previous article in this series showed that a closed boundary in three-dimensional space does not always define a unique surface, and worked out an exact bound on how much volume that ambiguity can cost. It stopped short of answering a more practical question, though: once you know a platform is non-planar, what do you actually do about it?

There are really only two answers worth taking seriously. Split the surface into pieces that are individually planar — triangles, most naturally, since three points always define a plane — and accept that one logical object becomes several. Or keep the object whole, and change its elevations until it fits a single plane. The first preserves the surface at the cost of the object. The second preserves the object at the cost of the surface. Neither is free, and this article is only concerned with the second.

Once you decide to fit a single plane, a question shows up immediately that sounds almost too simple to need an article: which plane? Not any old plane through the general vicinity of the points — the one that disturbs the platform’s design intent the least. That turns out to mean something quite specific, mathematically, and it’s not the plane most people would reach for first.

Conclusion first, as usual. For a set of points (x_i, y_i, z_i) with fixed plan positions, the plane

\displaystyle z = ax+by+c

that minimises the total squared elevation change

\displaystyle \sum_i \left(z_i-(ax_i+by_i+c)\right)^2

has a closed form that falls out of ordinary least squares — no iteration, no search, one pass to accumulate a handful of sums and a 2\times 2 system to solve. The interesting part isn’t the algebra. It’s the assumption baked into the objective function, what it rules out, and what happens when the input geometry pushes against the edge of what the model can represent.

---

Two Ways to Make a Boundary Behave #

Before getting into the fitting itself, it’s worth being explicit about what “resolving” a non-planar platform actually trades away, because the two available strategies fail in opposite directions.

Splitting is a choice — it doesn’t touch a single elevation the designer entered. What it costs is object identity: a platform carries a name, a mark, editing history, slope rules, reporting relationships. But the fragmentation into a pile of triangles is a problem; one has to ensure that the triangles somehow still link back to the original platform in some way. It may be fine for a severely warped platform, but overkill for one that’s a few millimetres out.

Fitting a plane is another choice. The platform stays one object, in name and in metadata, forever. What it costs is candour about the fact that the entered elevations are being changed on purpose. That’s not a defect to be minimised away quietly; it’s the whole point of choosing this method, and it should be visible to whoever clicks “apply.”

Neither method is universally correct, which is exactly why this has to be a choice and not a default. But once “fit a plane” is the chosen path, the mathematics of which plane stops being a matter of taste.

---

What “Best” Has to Mean Here #

Say the word “best-fit plane” to anyone with a numerical background and the instinct is total least squares — minimise the perpendicular distance from each point to the plane, the same idea behind principal component analysis on a 3D point cloud. It’s the right tool for fitting a plane through noisy survey data where every coordinate, x , y , and z , might be off by some measurement error.

However, it is the wrong tool here, and here’s the reason.

The operation under discussion never touches x or y . The platform’s plan outline — the shape a surveyor pegs out on site — is fixed. Only elevations move. So the quantity that should be minimised is elevation change, not 3D distance, and those are not the same thing once the plane has any slope at all.

Fitting orthogonal distance would, in effect, be allowed to nudge x and y a little to shorten that perpendicular line — which is exactly the freedom this operation doesn’t have. So the objective has to be vertical, ordinary regression of z against (x,y) , not a generic 3D plane fit. It’s a small distinction on paper and a decisive one in practice: get it backwards, and the “fitted” plane silently assumes a freedom to move the platform’s footprint that was never actually granted.

Put plainly: we’re not fitting a plane through a cloud of equally uncertain points. We’re asking, of all the planes that could sit under this fixed footprint, which one requires the smallest total disturbance to the elevations someone already typed in.

---

Setting Up the Least-Squares Problem #

In MiTS, every point on the platform that carries an independently entered elevation — as opposed to one interpolated from its neighbours — becomes one observation (x_i, y_i, z_i) , for i=1,\dots,n . Interpolated points don’t get a vote; they exist because of the labelled points, not alongside them as independent evidence, and letting them into the fit would mean using the platform’s own assumptions as if they were fresh data.

The model is a plane written as a height field:

\displaystyle z = ax + by + c.

For any candidate (a,b,c) , the residual at point i is how far that point’s actual elevation sits from the plane:

\displaystyle r_i = z_i – (ax_i+by_i+c).

The fitting objective is the sum of squared residuals,

\displaystyle S(a,b,c) = \sum_{i=1}^n r_i^2 = \sum_{i=1}^n \left(z_i-ax_i-by_i-c\right)^2,

and the task is to choose a , b , c to make S as small as possible.

Why squares, and not, say, the sum of absolute residuals? Mostly because the squared objective is smooth and convex — it has one minimum, reachable by ordinary calculus, with no piecewise cases to worry about. It also has the pleasant side effect of penalising a single badly-behaved point more than several mildly mismatched ones, which lines up with engineering intuition: a platform with one wild outlier vertex and three good ones is a worse candidate for “just fit a plane through it” than one where every vertex disagrees a little.

---

Deriving the Normal Equations #

Minimising S is a first-year-calculus exercise: set the partial derivatives with respect to a , b , and c to zero.

\displaystyle \frac{\partial S}{\partial c} = -2\sum_i \left(z_i-ax_i-by_i-c\right) = 0

\displaystyle \frac{\partial S}{\partial a} = -2\sum_i x_i\left(z_i-ax_i-by_i-c\right) = 0

\displaystyle \frac{\partial S}{\partial b} = -2\sum_i y_i\left(z_i-ax_i-by_i-c\right) = 0

Divide out the -2 and expand each one into three linear equations in a , b , c :

\displaystyle a\sum x_i + b\sum y_i + cn = \sum z_i

\displaystyle a\sum x_i^2 + b\sum x_iy_i + c\sum x_i = \sum x_iz_i

\displaystyle a\sum x_iy_i + b\sum y_i^2 + c\sum y_i = \sum y_iz_i

These are the normal equations — three linear equations, three unknowns, solvable directly. They always have a unique solution unless the underlying geometry is degenerate, which is a question worth its own section further down.

Solving this 3\times 3 system as written is completely valid. It’s also numerically the wrong way to do it on real project data, for a reason that has nothing to do with the mathematics and everything to do with how large the numbers get.

---

Centering: Taking the Origin Out of the Problem #

A site might sit hundreds of thousands, sometimes millions, of metres from the coordinate origin, because it’s tied to a national or state grid rather than a local one starting at zero. Plug those raw values straight into \sum x_i^2 and the numbers involved balloon — a coordinate of two and a half million squared is on the order of 6\times 10^{12} — while the actual quantity of interest, how the elevations deviate across a platform a few tens of metres wide, sits many orders of magnitude smaller. Subtracting numbers that large to recover a small, meaningful difference is precisely the situation floating-point arithmetic handles worst.

The fix doesn’t require new mathematics, only a change of variables. Work in coordinates measured from the mean of the points instead of from the world origin:

\displaystyle x_i’ = x_i-\bar x, \qquad y_i’ = y_i-\bar y, \qquad z_i’ = z_i-\bar z,

where \bar x , \bar y , \bar z are the ordinary averages over all n observations. Because a plane’s slope doesn’t depend on where you put the origin, this substitution doesn’t change a or b at all — only c shifts, and shifts in a way that’s easy to recover afterward. Substituting into the normal equations and using \sum x_i’ = \sum y_i’ = \sum z_i’ = 0 by construction, the first equation collapses entirely, and the remaining two reduce to a clean 2\times 2 system:

\displaystyle a\sum (x_i’)^2 + b\sum x_i’y_i’ = \sum x_i’z_i’

\displaystyle a\sum x_i’y_i’ + b\sum (y_i’)^2 = \sum y_i’z_i’

Every sum in this system is now a sum of centred quantities — differences from the mean — which stay comfortably small regardless of how far the site sits from the origin. The intercept is recovered afterward from the means:

\displaystyle c = \bar z – a\bar x – b\bar y.

This isn’t a numerical nicety tacked on for polish. It’s the difference between a fit that behaves identically whether the site sits at the origin or two and a half million metres away, and one that quietly degrades as coordinates grow — the exact translation-invariance property the previous article proved for the volume ambiguity bound turns up again here, for the fitting problem, for the same underlying reason: real project coordinates are large, and any formula that cares where zero happens to be is a formula waiting to misbehave on a real project.

---

The Closed-Form Solution #

Write the centred sums with shorthand, the way a statistics text would:

\displaystyle S_{xx}=\sum (x_i’)^2, \quad S_{yy}=\sum (y_i’)^2, \quad S_{xy}=\sum x_i’y_i’, \displaystyle S_{xz}=\sum x_i’z_i’, \quad S_{yz}=\sum y_i’z_i’.

The 2\times 2 system from the previous section is

\displaystyle \begin{pmatrix} S_{xx} & S_{xy} \ S_{xy} & S_{yy} \end{pmatrix} \begin{pmatrix} a \ b \end{pmatrix} = \begin{pmatrix} S_{xz} \ S_{yz} \end{pmatrix},

and inverting a 2\times 2 matrix is Cramer’s rule, nothing more exotic. With determinant

\displaystyle D = S_{xx}S_{yy}-S_{xy}^2,

the solution is

\displaystyle a = \frac{S_{xz}S_{yy}-S_{yz}S_{xy}}{D}, \qquad b = \frac{S_{yz}S_{xx}-S_{xz}S_{xy}}{D},

with c=\bar z-a\bar x-b\bar y as before.

That’s the entire algorithm. One pass over the points to accumulate five sums and three means, a 2\times 2 solve, and a final substitution. No iteration, no initial guess to seed, no convergence tolerance to fuss over — because the objective function is a quadratic bowl with exactly one minimum, calculus hands you the answer directly rather than making you search for it.

---

When the System Refuses to Invert #

The closed form above has one visible failure point: D=0 , where Cramer’s rule divides by zero. Let’s examine what it means.

Lemma #

D = S_{xx}S_{yy}-S_{xy}^2 = 0 if and only if every point’s plan position (x_i,y_i) lies on a single straight line.

Proof #

D is the Gram determinant of the two centred vectors \mathbf{u}=(x_1’,\dots,x_n’) and \mathbf{v}=(y_1’,\dots,y_n’) in \mathbb{R}^n , since S_{xx}=\mathbf u\cdot\mathbf u , S_{yy}=\mathbf v\cdot\mathbf v , and S_{xy}=\mathbf u\cdot\mathbf v . The Cauchy–Schwarz inequality gives

\displaystyle (\mathbf u\cdot\mathbf v)^2 \le (\mathbf u\cdot\mathbf u)(\mathbf v\cdot\mathbf v),

with equality exactly when \mathbf u and \mathbf v are parallel — that is, when there’s a constant \lambda with y_i’=\lambda x_i’ for every i (or \mathbf u=\mathbf 0 , all points sharing one x ). Either condition means every centred plan position is a scalar multiple of a single direction, which is exactly the statement that all the original (x_i,y_i) points lie on one line.

∎

So the fit doesn’t fail because of a numerical accident. It fails exactly when the plan footprint degenerates into a line — a platform with zero effective width. That’s a sensible place for the method to give up: with every point’s plan position falling on one line, there are infinitely many planes achieving the same minimal residual (tilt the plane any way you like across the direction perpendicular to that line, and nothing in the objective function notices), so “the” best-fit plane simply isn’t unique anymore. Rather than pick one arbitrarily and pretend it was forced, the right answer is to say the fit is undefined and let the caller decide what to do about a degenerate footprint. Refusing to guess at that point is not a shortcoming of the method — it is the method correctly reporting that the question has no single answer left to give.

Another interesting question: what about a plane that’s exactly vertical?

The plain answer is that the model z=ax+by+c simply cannot represent one, and that’s not an oversight to patch — it’s a structural fact about writing z as a function of (x,y) . A function assigns exactly one output to each input; a vertical plane, by contrast, contains an entire line’s worth of elevations sitting above a single (x,y) position, which is precisely a function of (x,y) failing to exist. No value of a or b , however large, reaches it — you’d need an infinite slope, which isn’t a number this system can produce.

But this is not really relevant as far as we are concerned. A platform, by its very definition, is a plan footprint with one elevation recorded per plan position — there is no way to hand the fitting routine two different elevations at the same (x,y) in the first place, so “what if the true best-fit plane is vertical” can’t arise from valid platform data at all. The closest a real, degenerate input can get is the collinear case proved above: as a platform’s footprint narrows toward a line, the fitted plane’s slope grows without bound, tilting ever closer to vertical — and exactly at the point where the footprint truly collapses to a line, the determinant hits zero and the method declines to answer, rather than returning an implausibly steep plane and calling it a fit.

---

What the Solution Means #

A few properties fall out almost for free once the closed form is in hand, and each has a direct engineering reading.

Corollary: Residuals Sum to Zero #

Rearranging the first normal equation, a\sum x_i+b\sum y_i+cn=\sum z_i , directly gives

\displaystyle \sum_i \left(z_i-(ax_i+by_i+c)\right) = \sum_i r_i = 0.

The fitted plane doesn’t systematically sit above or below the data — positive and negative elevation changes exactly cancel. Plainly: some points get nudged up, some get nudged down, and on balance the plane runs straight through the middle of them, not off to one side.

Corollary: An Already-Planar Platform Is a Fixed Point #

If the input points already lie exactly on some plane z=a_0x+b_0y+c_0 , every residual is zero for (a,b,c)=(a_0,b_0,c_0) , which drives the objective S to its global minimum of zero — and since S is a strictly convex quadratic with one minimiser, that minimiser has to be (a_0,b_0,c_0) itself. So a planar platform run through this fitting routine comes back unchanged, to floating-point precision, and running the fit a second time on an already-fitted platform changes nothing further. That idempotence matters in practice — it means the operation is safe to reapply, or to leave in place as a standing correction, without the result drifting further with each use.

Corollary: Large Coordinates Don’t Distort the Fit #

This was really proved already, in the centering section, but it’s worth stating as a standalone guarantee: shifting every point by a constant vector T=(T_x,T_y,T_z) leaves a and b completely unchanged, and shifts c by exactly T_z – aT_x – bT_y — the fitted plane simply rides along with the translation. A platform sitting near a national grid origin and the same platform re-expressed relative to a local site datum fit to the identical relative plane. The mathematics doesn’t know or care where the surveyor decided to put zero.

Corollary: Every Point Counts Equally — By Choice, Not by Necessity #

The objective function \sum r_i^2 gives every labelled point the same say in the outcome, regardless of how the points happen to be distributed across the footprint. That’s a deliberate modelling choice, not a mathematical inevitability — a weighted variant, minimising \sum w_i r_i^2 for some chosen weights w_i , is exactly as solvable in closed form, and would let a designer, say, discount a point known to be less reliable, or hold a corner more firmly in place than the rest. Equal weighting is simply the least assumption available: absent a stated reason to trust one entered elevation less than another, treating them identically is the right default, but it’s worth naming as a choice rather than letting it pass as the only option.

---

Worked Example: The Same Warped Face, a Different Question #

The previous article’s closing worked example lifted one corner of a unit square by a height d , to show how much volume ambiguity that single warped face could introduce. The same four points make an equally good worked example here, because the question this time is different: not how much does the volume disagree between interpretations, but what single plane best absorbs this platform’s non-planarity, and how much does each corner have to move to reach it?

Take the footprint

\displaystyle (0,0,0),\ (1,0,0),\ (1,1,d),\ (0,1,0),

three corners flat at zero, one corner lifted by d . The means are \bar x=\bar y=\tfrac12 , \bar z=\tfrac{d}{4} . Working through the centred sums:

\displaystyle S_{xx}=S_{yy}=1, \qquad S_{xy}=0, \qquad S_{xz}=S_{yz}=\frac d2.

S_{xy}=0 is worth noticing on its own — the plan footprint is a perfect square, so x and y carry no linear relationship to each other, and the 2\times2 system decouples into two independent one-variable fits, one for a and one for b . With D=S_{xx}S_{yy}-S_{xy}^2=1 ,

\displaystyle a=\frac{S_{xz}S_{yy}-S_{yz}S_{xy}}{D}=\frac d2, \qquad b=\frac{S_{yz}S_{xx}-S_{xz}S_{xy}}{D}=\frac d2, \qquad c=\bar z-a\bar x-b\bar y=-\frac d4.

So the best-fit plane is

\displaystyle z=\frac d2 x+\frac d2 y-\frac d4,

and evaluating it back at the four corners gives fitted elevations

\displaystyle -\frac d4,\quad \frac d4,\quad \frac{3d}{4},\quad \frac d4.

Set d=10 and this is z=5x+5y-2.5 , with fitted corners -2.5,\ 2.5,\ 7.5,\ 2.5 — clean numbers, which is why this particular footprint makes a convenient closed-form check for an implementation of the method.

The elevation changes at each corner — fitted minus original — come out to +\tfrac d4,\ -\tfrac d4,\ +\tfrac d4,\ -\tfrac d4 . They cancel exactly, as the residual-sum corollary promised, and by symmetry every corner moves the same distance, d/4 , in one direction or the other. For d=10 , that’s a maximum elevation change of \pm2.5 at every corner, and consequently an RMS elevation change of exactly d/4=2.5 as well — same figure, because every residual happens to have the same magnitude here.

Lifted corner d	Max elevation change	RMS elevation change
0.0	0.000	0.000
0.1	0.025	0.025
0.3	0.075	0.075
1.0	0.250	0.250
10.0	2.500	2.500

Compare this against the previous article’s result for the identical face: the maximum volume ambiguity there was \Delta V_{\max}=d/6 . Here, the maximum elevation change needed to remove that ambiguity entirely is d/4 . Different quantities, different units, and — worth sitting with for a moment — the elevation fix is proportionally larger relative to d than the volume ambiguity it resolves. Flattening a platform isn’t a free rounding operation; it’s a real, quantifiable edit to the design, and now there’s a number attached to exactly how large that edit is, for any platform, before anyone has to accept it.

---

Discussion #

But won’t large deviation change the definition of “best fit” plane? #

What if the elevations are deviating a lot from the fitted plane? Big deviations change the fitted plane, the changed plane changes what counts as a deviation, and round it goes — plane depends on deviation, deviation depends on plane, nothing ever pins down. If that were actually how the method worked, “undefined” would be the right word for it, and the closed form from earlier would be some kind of lie.

It isn’t how the method works, and the fix is to look at what a , b , and c are actually functions of. The objective function is

\displaystyle S(a,b,c)=\sum_i\left(z_i-ax_i-by_i-c\right)^2,

and every z_i in that sum is a fixed, already-known number — the elevation the designer typed in, before any fitting has happened. Nothing about minimising S changes those z_i . The minimisation solves for a , b , c once, treating the input data as constant throughout, and the residuals only get computed afterward, by plugging the now-fixed a , b , c back into the same fixed z_i . Data goes in, coefficients come out, residuals get read off the result. One direction, one pass, no step where the output is fed back in as new input.

Real feedback loops in curve and surface fitting do exist — iteratively reweighted least squares and RANSAC both refit against their own previous output, and both need a stopping rule and a convergence argument because of it. Neither is what’s happening here. This method is ordinary least squares, run exactly once, against elevations that never move during the computation. So “the deviation is too big” doesn’t send the calculation into an unresolved spiral — it just produces one specific, fully determined plane, whose residuals happen to be large. Large is not the same as undefined.

That’s really the whole correction: a big deviation changes the outcome of the fit, in a perfectly ordinary cause-and-effect way — it doesn’t change whether the fit has an outcome. The single failure mode that actually leaves a , b , c undetermined was already covered above, and it has nothing to do with deviation size: it’s the plan footprint collapsing onto a line. Short of that, however violently the input elevations disagree with any plane, the one-shot solve still lands on exactly one answer.

---

Does This Require the Deviation to Be Small? #

The derivation can look like it’s quietly leaning on that assumption somewhere. It isn’t, and here’s exactly where “small” does and doesn’t enter.

Nowhere in the calculus does a smallness assumption appear. Minimising S(a,b,c)=\sum(z_i-ax_i-by_i-c)^2 is an exact global minimisation for whatever values the z_i happen to take — a platform whose worst corner is out by a millimetre and one whose worst corner is out by ten metres are solved by the identical closed form, with no linearisation, no tangent-plane approximation, no truncated series anywhere in sight. The worked example even said as much without dwelling on it: the formulas hold for any d , not just small ones, and the table ran the same closed form out to d=10 without it needing a caveat.

What isn’t answered anywhere above is a different question: given that this method will always dutifully return a plane and a residual, at what residual does returning one stop being the sensible thing to do? That threshold doesn’t exist in this article, and it can’t be manufactured out of the mathematics, because it isn’t a mathematical question. A plane fit to a platform where every corner needs to move a millimetre is obviously fine. A plane fit to a platform where the worst corner needs to move a metre is a different design, wearing the original platform’s name. Somewhere between those two, a line gets crossed, and the fitting algorithm has no opinion about where — it will compute a , b , c with equal confidence on either side of it.

That’s not a hole in the theory so much as a boundary the theory was never going to fill. The earlier article in this series landed on the same conclusion from the opposite direction, when it pointed out that a 0.5 mm deviation might be everything on a machine part and noise on a hectare-sized earthwork platform — there’s no universal number, only a number relative to what the project actually needs. What this article contributes is the other half of that judgement call: an exact max and RMS elevation change, computed before a single point moves, to weigh against whatever tolerance the project has decided matters.

The mathematics supplies the number. Deciding how large a number is still acceptable stays an engineering call, not something a closed-form solution gets to hand down. But the previous article can answer a related one: from that same elevation deviation, it lets us compute the potential volume deviation directly.

---

Appendix: Summary of Key Equations #

Objective function:

\displaystyle S(a,b,c)=\sum_i\left(z_i-ax_i-by_i-c\right)^2.

Normal equations:

\displaystyle a\sum x_i+b\sum y_i+cn=\sum z_i, \displaystyle a\sum x_i^2+b\sum x_iy_i+c\sum x_i=\sum x_iz_i, \displaystyle a\sum x_iy_i+b\sum y_i^2+c\sum y_i=\sum y_iz_i.

Centred sums:

\displaystyle S_{xx}=\sum(x_i-\bar x)^2,\ S_{yy}=\sum(y_i-\bar y)^2,\ S_{xy}=\sum(x_i-\bar x)(y_i-\bar y), \displaystyle S_{xz}=\sum(x_i-\bar x)(z_i-\bar z),\ S_{yz}=\sum(y_i-\bar y)(z_i-\bar z).

Closed-form solution:

\displaystyle a=\frac{S_{xz}S_{yy}-S_{yz}S_{xy}}{S_{xx}S_{yy}-S_{xy}^2}, \qquad b=\frac{S_{yz}S_{xx}-S_{xz}S_{xy}}{S_{xx}S_{yy}-S_{xy}^2}, \qquad c=\bar z-a\bar x-b\bar y.

Degeneracy condition:

\displaystyle S_{xx}S_{yy}-S_{xy}^2=0 \iff \text{all } (x_i,y_i) \text{ collinear.}