Аппроксимация точек окружностью

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

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

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

  n
r 2 =   1

n

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

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

n
∑ ( 2 ( xi - x ) x + 2 ( yi - y ) y + xx - xi 2 + yy - yi 2 ) 2      (3)
i = 1
где
  n
x =  1

n

 ∑xi
  i = 1
  n
y =  1

n

 ∑yi
  i = 1
  n
xx =  1

n

 ∑xi2
  i = 1
  n
yy =  1

n

 ∑yi2
  i = 1

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

Def<Circle2d> getCirclePnt22 ( CArrRef<Vector2d> p );

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

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

В третьем подходе будем искать минимум суммы абсолютных разностей между радиусом окружности r и расстоянием от точки pi до центра окружности о:

n
∑ | | o - pi |2 - r |      (6)
i = 1
Эту задачу можно решить перебором троек точек, построении по ним окружностей и выборе среди них оптимальной: 
Def<Circle2d> getCirclePnt1 ( CArrRef<Vector2d> p, nat & ix1, nat & ix2, nat & ix3 );
Def<Circle2d> getCirclePnt1 ( CArrRef<Vector2d> p );
Временная сложность такого алгоритма равна O ( n4 ).

Замечу, что окружности построенными этими тремя алгоритмами могут очень сильно отличаться.

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

Исходники находятся в approx2d.cpp

Наверх