Однажды встала задача по набору измеренных точек построить аппроксимирующую плоскость в пространстве. Можно было бы конечно численными методами подогнать плоскость к точкам, но хотел получить готовую точную формулу (как для аппроксимации прямой).
Я не нашёл простой готовой формулы, либо предложенные методы были довольно сложные для понимания, например, среди решений фигурировала мнимая единица (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’ стремится к нулю.
