Однажды встала задача по набору измеренных точек построить аппроксимирующую плоскость в пространстве. Можно было бы конечно численными методами подогнать плоскость к точкам, но хотел получить готовую точную формулу (как для аппроксимации прямой).

Я не нашёл простой готовой формулы, либо предложенные методы были довольно сложные для понимания, например, среди решений фигурировала мнимая единица (i² = -1), которую предлагалось потом каким-то образом обходить.

Поэтому я решил вывести формулу самостоятельно. Оказалось всё довольно просто.

(Не пугайтесь, в конце статьи будет готовый к применению алгоритм расчета)

Итак:

Задача такая – есть некоторое количество точек M(x,y,z) в трехмерном пространстве, нужно найти плоскость такую, чтобы сумма квадратов расстояний от точек M до этой плоскости была минимальной.

1) Уравнение плоскости:

Ax+By+Cz+D = 0

2) Расстояние от произвольной точки M(x,y,z) до плоскости:

d = Ax+By+Cz+D

3) Потребуем выполнения условия:

S(A,B,C,D) = Σd² = Σ(Ax+By+Cz+D)² = min

4) Раскроем скобки для d² (разными способами, комбинируя разные слагаемые)

d² = (Ax+(By+Cz+D))² = A²x²+2Ax(By+Cz+D)+(By+Cz+D)²

d² = (By+(Ax+Cz+D))² = B²y²+2By(Ax+Cz+D)+(Ax+Cz+D)²

d² = (Cz+(Ax+By+D))² = C²z²+2Cz(Ax+By+D)+(Ax+By+D)²

d² = (D+(Ax+By+Cz))² = D²+2D(Ax+By+Cz)+(Ax+By+Cz)²

4) Функция S(A,B,C,D) должна иметь локальный минимум. Найдём такие A,B,C,D, чтобы производные функции S по её аргументам A,B,C и D равнялись нулю.

dS/dA = Σ(2Ax²+2x(By+Cz+D)) = 0

dS/dB = Σ(2By²+2y(Ax+Cz+D)) = 0

dS/dC = Σ(2Cz²+2z(Ax+By+D)) = 0

dS/dD = Σ(2D+2(Ax+By+Cz)) = 0

5) Получилась система из четырёх уравнений с четырьмя неизвестными A,B,C,D

AΣx²+BΣxy+CΣzx+DΣx = 0

AΣxy+BΣy²+CΣyz+DΣy = 0

AΣzx+BΣyz+CΣz²+DΣz = 0

AΣx+BΣy+CΣz+DN = 0, (N - число точек)

6) Из уравнения 4. выразим D через A,B,C и подставим его в уравнения 1.2.3.

D = −(AΣx+BΣy+CΣz)/N

AΣx²+BΣxy+CΣzx−Σx(AΣx+BΣy+CΣz)/N = 0

AΣxy+BΣy²+CΣyz−Σy(AΣx+BΣy+CΣz)/N = 0

AΣzx+BΣyz+CΣz²−Σz(AΣx+BΣy+CΣz)/N = 0

Или

A(Σx²−(Σx)²/N)+B(Σxy−ΣxΣy/N)+C(Σzx−ΣzΣx/N) = 0

A(Σxy−ΣxΣy/N)+B(Σy²−(Σy)²/N)+C(Σyz−ΣyΣz/N) = 0

A(Σzx−ΣzΣx/N)+B(Σyz−ΣyΣz/N)+C(Σz²−(Σz)²/N) = 0

7) Введём обозначения для известных чисел (их можно просто вычислить из массива точек, просуммировав соответствующие выражения):

"XX" = Σx²−(Σx)²/N

"YY" = Σy²−(Σy)²/N

"ZZ" = Σz²−(Σz)²/N

"XY" = Σxy−ΣxΣy/N

"YZ" = Σyz−ΣyΣz/N

"ZX" = Σzx−ΣzΣx/N

8) Система уравнения принимает вид:

A∙XX+B∙XY+C∙ZX = 0

A∙XY+B∙YY+C∙YZ = 0

A∙ZX+B∙YZ+C∙ZZ = 0

9) Дальше я руками ввёл некоторую δ, чтобы избежать тривиального решения A=B=C=0, а дополнительное условие нормировки позволит определить эту δ. Это - тонкое место, я не могу объяснить, в чем здесь фокус, но результат получается правильный.

A∙XX+B∙XY+C∙ZX = δ

A∙XY+B∙YY+C∙YZ = δ

A∙ZX+B∙YZ+C∙ZZ = δ

A²+B²+C² = 1 дополнительное условие нормировки

10) Введём обозначения для определителей матриц:

11) тогда решение системы уравнений можно записать так:

δ = 1/√(DetA(1)²+DetB(1)²+DetB(1)²)

A = DetA(δ)

B = DetB(δ)

C = DetC(δ)

мы также помним: D = −(AΣx+BΣy+CΣz)/N

На самом деле не нужно вычислять величину δ, вместо этого следует действовать по следующему алгоритму:

АЛГОРИТМ:

Шаг 1. Вычисляем числа "XX", "YY", "ZZ", "XY", "YZ", "ZX" в соответствии с пунктом 7

Шаг 2 Вычисляем определители, составленные из этих чисел и единиц:

Шаг 3 Вычисляем ненормированные коэффициенты:

A1=DetA(1)

B1=DetB(1)

C1=DetC(1)

Шаг 4 Вычисляем длину вектора (A1,B1,C1)

L = √(A1²+B1²+C1²)

Шаг 2’ Вычисляем определители другие (в первых строчках минусы перед единицами):

Шаг 3’ Вычисляем ненормированные коэффициенты со штрихами:

A1’=DetA’(1)

B1’=DetB’(1)

C1’=DetC’(1)

Шаг 4’ Вычисляем длину вектора (A1’,B1’,C1’)

L’ = √(A1’²+B1’²+C1’²)

Шаг 5 Сравниваем, что больше L, или L’

Шаг 6 Если L > L’ тогда:

A = A1/L

B = B1/L

C = C1/L

Шаг 6’ Если L < L’ тогда:

A = A1’/L’

B = B1’/L’

C = C1’/L’

Шаг 7 Вычисляем D

D = −(AΣx+BΣy+CΣz)/N

Искомое уравнение плоскости: Ax+By+Cz+D = 0

Здесь я алгоритм разбил на две ветки со штрихом и без штриха, потому что L может оказаться близко к нулю для некоторых углов (если очень не повезет, то L может быть вообще равна 0). Тогда при делении на величину, близкую к нулю мы потеряем точность.

Такое случается, например, при повороте плоскости на 45 градусов вокруг одной из осей,

Например, плоскости A=−√2, B=0, С=√2 и при малом случайном разбросе точек около плоскости величина L стремится к нулю, в то время, как L’ как раз максимальна. И наоборот для плоскости A=√2, B=√2, С=0 величина L’ стремится к нулю.