Аппроксимация точек сферой

Пусть в пространстве дано множество из n точек pi = ( xi, yi, zi ) и нужно аппроксимировать его сферой. Рассмотрим несколько подходов к решению этой задачи. В первом подходе координаты центра o = ( x, y, z ) и радиус r будем искать, как аргументы при которых функция

n
( r 2 - ( x - xi ) 2 - ( y - yi ) 2 - ( z - zi ) 2 ) 2      (1)
i = 1

принимает минимальное значение. Сведём эту задачу к задаче наименьших квадратов. Для этого продифференцируем (1) по r и приравняем производную нулю. Отсюда получим выражение r 2 через x, y и z:

  n
r 2 =   1

n

( ( x - xi ) 2 + ( y - yi ) 2 + ( z - zi ) 2 )     (2)
  i = 1

Подставим это выражение в (1) и приведём подобные члены. Получится задача наименьших квадратов:

n
( 2 ( xi - x ) x + 2 ( yi - y ) y + 2 ( zi - z ) z + xx - xi 2 + yy - yi 2 + zz - zi 2 ) 2      (3)
i = 1

где

  n
x 1

n

xi
  i = 1
  n
y 1

n

yi
  i = 1
  n
z 1

n

zi
  i = 1
  n
xx 1

n

xi2
  i = 1
  n
yy 1

n

yi2
  i = 1
  n
zz 1

n

zi2
  i = 1

Следующая программа решает эту задачу:

Def<Sphere3d> getSpherePnt22 ( CArrRef<Vector3d> p );

Основное преимущество рассмотренной задачи - это простота решения, но более интересной для практики является минимизация другой функции отклонений:

n
( | o - pi |2 - r ) 2      (4)
i = 1
О векторных нормах смотрите здесь. Эту формулу можно переписать так:
n
( cxi * ( х - xi ) + cyi * ( y - yi ) + czi * ( z - zi ) - r ) 2      (5)
i = 1
где
cxi = ( х - xi ) / | o - pi |2
cyi = ( y - yi ) / | o - pi |2
czi = ( z - zi ) / | o - pi |2
Если считать в формуле (5) переменные cxi, cyi и czi известными, то приходим к классической задаче МНК.
Отсюда получается следующий итерационный алгоритм:
1) Получаем начальное приближение для o и r при помощи функции getSpherePnt22.
2) Вычисляем значения cxi, cyi и czi.
3) Решая задачу (5) получаем новые значения для o и r.
Будем повторять шаги 2 и 3 пока алгоритм не сойдётся. В результате получаем программу:
Def<Sphere3d> getSpherePnt2 ( CArrRef<Vector3d> p );

Описание класса Vector3d находится здесь.
Описание класса Sphere3d находится здесь.
Описание шаблона классов Def находится здесь.
Описание шаблона классов CArrRef находится здесь.

Исходники находятся в файле approx3d.cpp

Наверх