Best fit plane for a nonplanar platform #

A platform with independently entered corner elevations may not lie on a single plane. As shown in the previous article, its boundary can then admit different surfaces, giving different earthwork volumes.

We can resolve this by dividing the platform into planar faces, usually triangles, or by adjusting its elevations to fit one plane. Triangulation retains the entered elevations, but we must choose how the points connect and keep the resulting faces associated with the platform. Fitting a plane keeps the platform whole and changes its elevations.

This article derives a plane fit that keeps every plan position fixed and minimises the sum of squared elevation changes. It gives a unique answer when the input points are not collinear in plan. Whether those changes are acceptable depends on the project’s design requirements.

Splitting also affects how the platform’s name, mark, editing history, slope rules and reports are associated with its faces. Fitting avoids that fragmentation, but the elevation changes should be shown before they are applied. The designer must choose which operation suits the platform.

What are we minimising? #

For each input point (x_i,y_i,z_i) , we keep x_i and y_i fixed and replace its elevation with a value on the plane

\displaystyle z=ax+by+c.

The residual is the original elevation minus the fitted elevation:

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

Here \hat z_i denotes the fitted elevation at point i . The coefficients a and b are the slopes in the x and y directions, and c is the plane’s elevation at (0,0) .

We choose a , b and c to minimise

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

This is ordinary least squares. Squaring the residuals prevents positive and negative differences from cancelling and gives larger differences more influence on the fit. Each entered elevation receives equal weight.

We could instead minimise the sum of absolute residuals. The squared objective is smooth and convex, which lets us solve it by setting its derivatives to zero. It has a unique minimum when the plan positions are non-collinear. Its stronger penalty on large residuals also makes the result more sensitive to an outlying elevation.

A plane fitted by perpendicular distance uses a different objective. For a sloping plane, the perpendicular distance is shorter than the vertical distance, and their ratio changes with the slope. Here, the quantity we need to control is the elevation adjustment at a fixed plan position, so we minimise the vertical residual.

Perpendicular-distance fitting is commonly called orthogonal or total least squares and is related to principal component analysis (PCA). It can suit measurements with comparable uncertainty in all three coordinates. Its closest point on the plane generally has different plan coordinates from the original point.

In MiTS, the observations are the platform points with independently entered elevations. Interpolated elevations are excluded because they are derived from those inputs; including them would give extra weight to parts of the platform simply because they contain more interpolated points.

Deriving the fitted plane #

At the minimum, the partial derivatives of S with respect to a , b and c are zero:

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

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

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

Dividing by -2 gives

\displaystyle \sum_i r_i=0, \qquad \sum_i x_i r_i=0, \qquad \sum_i y_i r_i=0.

Expanding these gives three linear equations for the three coefficients:

\displaystyle a\sum_i x_i+b\sum_i y_i+cn=\sum_i z_i,

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

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

The first equation gives

\displaystyle c=\bar z-a\bar x-b\bar y,

where \bar x , \bar y and \bar z are the means of the input coordinates. Thus the fitted plane passes through (\bar x,\bar y,\bar z) .

Centre the coordinates #

Survey coordinates can be millions of metres from the origin, while a platform may be only tens of metres wide. To reduce numerical errors when calculating the slopes, we work with differences from the means:

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

Substituting the expression for c into the residual gives

\displaystyle r_i=z_i’-ax_i’-by_i’.

Since the centred coordinates sum to zero, the remaining normal equations become

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

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

Define the five centred sums:

\displaystyle S_{xx}=\sum_i(x_i’)^2, \qquad S_{yy}=\sum_i(y_i’)^2, \qquad S_{xy}=\sum_i x_i’y_i’,

\displaystyle S_{xz}=\sum_i x_i’z_i’, \qquad S_{yz}=\sum_i y_i’z_i’.

We now have two equations for the slopes:

\displaystyle aS_{xx}+bS_{xy}=S_{xz}, \qquad aS_{xy}+bS_{yy}=S_{yz}.

In matrix form:

\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}.

Solve for the slopes #

With determinant

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

the solution, provided D>0 , 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}.

The fitted elevation at any plan position can then be evaluated directly in centred form:

\displaystyle \hat z(x,y)=\bar z+a(x-\bar x)+b(y-\bar y).

This avoids introducing the potentially large intercept c into each elevation calculation. If needed, recover it from c=\bar z-a\bar x-b\bar y .

The calculation requires the coordinate means, the five centred sums and a two-equation solve. No iterative fitting is needed. Changing the coordinate origin leaves the slopes and relative fitted elevations unchanged in exact arithmetic; centring reduces the numerical errors associated with large coordinates.

When is the solution unique? #

The determinant D is zero exactly when all the input plan positions lie on one straight line, including the case where they coincide.

To see why, let \mathbf u=(x_1′,\ldots,x_n’) and \mathbf v=(y_1′,\ldots,y_n’) . The Cauchy–Schwarz inequality gives

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

so D\geq0 . Equality holds when these vectors are linearly dependent. This means there are constants p,q , not both zero, such that px_i’+qy_i’=0 for every point: all plan positions lie on a line through their mean.

For collinear inputs, the elevations provide no information about the slope across that line. Infinitely many planes give the same minimum residual, so we cannot select a unique plane from those inputs alone. At least three non-collinear plan positions are needed.

Nearly collinear inputs also need care. Small elevation changes can cause large changes in the slope across a narrow footprint. Centring helps with large coordinates, but it cannot supply the missing spread of points across the platform.

The model z=ax+by+c represents a surface with one elevation per plan position. A vertical plane contains multiple elevations at the same plan position, so it cannot be represented by finite values of a and b . The fitted slopes can become very large for narrow footprints, but they need not do so: the elevations also determine the result.

Properties of the fit #

Elevation changes balance at the input points #

The first normal equation gives

\displaystyle \sum_i r_i=\sum_i\left(z_i-\hat z_i\right)=0.

Since the applied change is \Delta z_i=-r_i , upward and downward changes sum to zero. This is an equal-weight balance at the input points; it does not establish that earthwork volume is unchanged.

Planar inputs remain unchanged #

If all input points already lie on z=a_0x+b_0y+c_0 , that plane gives S=0 , the lowest possible value. With non-collinear plan positions the solution is unique, so the fit returns the same plane. Fitting the adjusted points again therefore makes no further change, apart from numerical rounding.

Moving the origin preserves the fit #

If we translate every input by (T_x,T_y,T_z) , the centred coordinates and slopes remain unchanged. The intercept becomes

\displaystyle c_{\mathrm{new}}=c+T_z-aT_x-bT_y.

Thus the fitted plane moves with the input points. This identity holds in exact arithmetic; centring reduces numerical errors when working with large survey coordinates.

Points can be weighted differently #

The objective used here gives each independently entered elevation equal weight. A weighted fit instead minimises

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

Positive weights allow more reliable or more important elevations to influence the fit more strongly. Using weighted means and weighted centred sums gives the same form of solution. A larger weight reduces a point’s freedom to move; holding a level exactly requires an explicit constraint.

Does the deviation have to be small? #

The derivation makes no small-deviation assumption. It minimises the full quadratic objective, with no linearisation or truncated approximation. Large elevation differences can produce large residuals, but do not destroy uniqueness when the plan positions are non-collinear.

The input elevations remain fixed during the calculation. We solve once for the plane coefficients, then calculate the residuals. Large deviations do not create a feedback loop or require repeated fitting.

Worked example: one raised corner #

Take a unit square with three corners at elevation zero and the fourth raised by d :

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

The coordinate means are

\displaystyle \bar x=\bar y=\frac12, \qquad \bar z=\frac d4,

and the centred sums are

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

Here S_{xy}=0 , so the two slope equations separate: a=S_{xz}/S_{xx}=d/2 and b=S_{yz}/S_{yy}=d/2 . With D=1 and c=\bar z-a\bar x-b\bar y=-d/4 , the fitted plane is

\displaystyle z=\frac d2x+\frac d2y-\frac d4.

Evaluating it at each corner gives:

Plan position	Original elevation	Fitted elevation	Change (fitted minus original)
(0,0)	0	-d/4	-d/4
(1,0)	0	d/4	+d/4
(1,1)	d	3d/4	-d/4
(0,1)	0	d/4	+d/4

For a 1 m × 1 m square with corner C raised by 100 mm, use d=0.1 m. With all coordinates expressed in metres, the fitted plane is

\displaystyle z=0.05x+0.05y-0.025.

Corner	Original elevation (mm)	Fitted elevation (mm)	Adjustment (mm)
A (0,0)	0	−25	−25
B (1,0)	0	25	+25
C (1,1)	100	75	−25
D (0,1)	0	25	+25

Every corner moves by 25 mm: two rise and two fall. The maximum and root-mean-square (RMS) elevation changes are both 25 mm. For example, at C the residual is r_C=100-75=+25 mm, while the applied adjustment is \Delta z_C=75-100=-25 mm.

The same formula applies for any lift d . For d=10 m, the fitted plane is z=5x+5y-2.5 (coordinates in metres), giving corner elevations -2.5,2.5,7.5,2.5 m. This deliberately large lift also checks that the calculation does not depend on small deviations.

Lifted corner d (m)	Maximum absolute adjustment (m)	RMS adjustment (m)
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

Both measures equal |d|/4 here because all four changes have that magnitude. For the 100 mm lift, RMS is \sqrt{(25^2+25^2+25^2+25^2)/4}=25 mm.

For this same 1 m × 1 m footprint, the previous article gives a maximum volume ambiguity of \Delta V_{\max}=(1\ \mathrm{m}^2)|d|/6 . A 100 mm lift therefore gives about 0.0167\ \mathrm{m}^3 of volume ambiguity, while fitting the plane changes each corner by 25 mm. The two numbers describe different effects: uncertainty in volume before choosing a surface, and the elevation edits required to make the boundary planar.