Что такое аппроксимирующая окружность

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

n
∑	( r 2 — ( x — x i ) 2 — ( y — y i ) 2 ) 2				(1)
i = 1

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

n
r 2 =	1	∑	( ( x — x i ) 2 + ( y — y i ) 2 )				(2)
i = 1

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

n
∑	( 2 ( x i — x ) x + 2 ( y i — y ) y + xx — x i 2 + yy — y i 2 ) 2				(3)
i = 1

где

n
x =	1	∑	x i
i = 1

n
y =	1	∑	y i
i = 1

n
xx =	1	∑	x i 2
i = 1

n
yy =	1	∑	y i 2
i = 1

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

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

n
∑	( | o — p i | 2 — r ) 2				(4)
i = 1

О векторных нормах смотрите здесь. Эту формулу можно переписать так:

n
∑	( c i * ( х — x i ) + s i * ( y — y i ) — r ) 2				(5)
i = 1

где

c i = ( х — x i ) / | o — p i | 2

s i = ( y — y i ) / | o — p i | 2

Если считать в формуле (5) переменные c i и s i известными, то приходим к классической задаче МНК. Отсюда получается следующий итерационный алгоритм:
1) Получаем начальное приближение для x, y и r при помощи функции getCirclePnt22.
2) Вычисляем значения c i и s i .
3) Решая задачу (5) получаем новые значения для x, y и r.
Будем повторять шаги 2 и 3 пока алгоритм не сойдётся. В результате получаем программу:

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

n
∑	| | o — p i | 2 — r |				(6)
i = 1

Эту задачу можно решить перебором троек точек, построении по ним окружностей и выборе среди них оптимальной: Временная сложность такого алгоритма равна O ( n 4 ).

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

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

Аппроксимация дугами окружностей

Динамическое программирование может применяться и для преобразования «цепочки» в проектную линию, состоящую из круговых кривых, сопрягаемых отрезками прямых. Это частный случай более общей задачи аппроксимации плоской кривой, заданной дискретно, кривой соответствующего вида, которая может и не быть графиком однозначной функции.

Такие задачи возникают при проектировании плана трассы. Неоднозначность аппроксимируемых функций создаёт дополнительные сложности.

Рассмотрим такую задачу, в которой элементами аппроксимирующей кривой являются дуги окружностей, сопрягаемые встык или отрезками прямых.

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

• 1. Отклонения проектной кривой от исходной в заданных точках не должны превосходить заданных величин. В частности, возможно задание фиксированных точек.
• 2. Кривизна аппроксимирующей кривой во всех точках по абсолютной величине не больше заданной.
• 3. Длины элементов не должны быть меньше заданных величин.
• 4. На границах элементы имеют общую касательную.
• 5. При проектировании продольного профиля дополнительно имеется ограничение по максимальному уклону.

Считаются заданными начальная и конечная точка искомой кривой, а также начальное и конечное направления.

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

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

Для каждой нормали задаётся угол у j с осью X и разбивается сетка с шагом А относительно исходной ломаной. Декартовы координаты узлов сетки на каждой нормали вычисляются по формулам

Рис. 8.8. Первый шаг алгоритма

где п/- номер i-ой точки на j -ой нормали ( нумерация начинается с нуля).

При проектировании продольного профиля все yj =я/2 и рассматривается обычная прямоугольная сетка варьирования.

Очередной этап динамического программирования — это переход к следующей нормали, но не к следующему элементу, число которых неизвестно; и к тому же в процессе построения аппроксимирующей кривой в одну точку могут приходить варианты, состоящие из различного числа элементов. Основное понятие в динамическом программировании: “состояние системы”- в данной задаче это точка на нормали и угол касательной к окружности в этой точке с осью X. Другими словами, задача рассматривается как двухпараметрическая. Начальное состояние точка А и угол а, конечное состояние точка В и угол р (рис.8.8) Рассмотрим сначала задачу, в которой допускается сопряжение дуг окружностей без прямых вставок.

Заданную минимальную длину элемента обозначим через L_mj_n. Максимальную длину элемента L_max примем равной 2L_mj_n, так как элементы длиной больше, чем 2 L_mj_n можно рассматривать как два различных элемента.

Первый шаг алгоритма динамического программирования. Для каждого положения точки D(x_d,y_d) первой нормали по её декартовым координатам, координатам точки А и углу а строится дуга окружности (рис.8.9).

Рис.8.9. Построение первого элемента

При этом сначала определяются координаты центра окружности х_с, у_с, а затем её радиус R.

R=AC; При ] на рис. 8.10. Тем самым “состояние системы ” то же, что и ранее (точка +угол), но углы не вычисляются, а назначаются с заданной дискретностью. Фактически изначально разбивается сетка не только по нормали с шагом А, но и по углам с шагом е.

По заданному начальному (х_с,у_с,ф|) и конечному состоянию (xa,yd, фД однозначно определяются все параметры связки: сначала вершина угла (координаты х_ь , уь точки В на рис. 8.10) и сам угол поворота со, затем длина ВС и радиус окружности R , затем длина дуги окружности L_Kp и длина отрезка прямой Ь_л.

Приведенные формулы должны быть преобразованы, если cpi или ф j равны ± я/2.

Далее проверяется выполнение всех ограничений. Допустимые соединения оцениваются по критерию оптимальности и для каждой точки и угла в конце связки запоминаются точка и угол наилучшего начала. Задача остаётся двухпараметрической, но число вариантов и сложность оценки вариантов возрастают по сравнению с аппроксимацией только круговыми элементами.

Может показаться, что изложенные алгоритмы требует полного перебора всех возможных комбинаций углов при заданных точках начала и конца элемента при сопряжении встык или при построении связки. Однако это не так. Поскольку перебор упорядочен от меньших значений к большим, то при нарушении ограничений может быть выявлена бессмысленность дальнейшего увеличения варьируемого параметра, так как при этом невязка в ограничении увеличивается. В этом случае требуется перейти к следующему состоянию. Частный случай такой ситуации показан на рис. 8.11.

На рис. 8.11. точки С и D рассматриваются как концы очередной дуги окружности (при сопряжении окружностей встык). Прямая ВС соответствует углу касательной с осью X в точке С, при котором радиус окружности R принимает минимально допустимое значение.

Рис. 8.11. Прекращение перебора при нарушении ограничения

Исследование неитерационного метода наименьших квадратов для оценивания параметров аппроксимирующей окружности Текст научной статьи по специальности « Физика»

Аннотация научной статьи по физике, автор научной работы — Михляев С. В.

Приведены аналитические выражения для смещений оценок параметров аппроксимирующей окружности, получаемых при использовании неитерационного метода наименьших квадратов. Для компенсации смещений предложена комбинированная схема оценивания по выборкам из различного количества точек.

Похожие темы научных работ по физике , автор научной работы — Михляев С. В.

Investigation of a non-iterative technique of least squares circle fitting

Analytical expressions for least squares estimation of parameters of the fitted circle and their biases produced by a non-iterative technique are given. To compensate biases, the complex estimation technique utilizing samples of various quantities of approximated points is offered.

Текст научной работы на тему «Исследование неитерационного метода наименьших квадратов для оценивания параметров аппроксимирующей окружности»

Исследование неитерационного метода наименьших квадратов для оценивания параметров аппроксимирующей окружности

Институт автоматики и электрометрии СО РАН Новосибирск, Россия

Analytical expressions for least squares estimation of parameters of the fitted circle and their biases produced by a non-iterative technique are given. To compensate biases, the complex estimation technique utilizing samples of various quantities of approximated points is offered.

Необходимость оценки параметров окружности, аппроксимирующей данные измерений, возникает в различных задачах, связанных с техническим контролем, обработкой изображений и т.д. 3. В силу нелинейности решаемой при этом оптимизационной задачи, для определения радиуса и координат центра аппроксимирующей окружности, как правило, применяются различные итерационные алгоритмы, основанные, в частности, на методе наименьших квадратов (МНК) и минимизирующие среднеквадратичное отклонение точек от аппроксимирующей окружности 7. Основной недостаток таких алгоритмов — значительные вычислительные затраты, которые не всегда допустимы в практических приложениях. Возможно также использование прямого неитерационного метода оценки параметров окружности (МНК2), являющегося по сути методом наименьших квадратов, примененным не к расстояниям, а к квадратам расстояний между точками [1, 10]. МНК2 позволяет получить решение в аналитическом виде, но дает смещенные оценки параметров [9, 10]. Смещение особенно существенно при аппроксимации фрагмента окружности в виде дуги с небольшим значением центрального угла. Для уменьшения смещения оценок параметров разработаны различные подходы, требующие, однако, дополнительных вычислительных ресурсов [10].

В работе [5] показано, что существенное уменьшение величины смещения оценок в МНК2 (МНК2-оценок) обеспечивается при использовании четырехточечной схемы аппроксимации дуги окружности. Погрешность определения параметров окружности при этом может быть допустимой для различных практических приложений. Результаты работы [5] получены путем численного моделирования в предположении независимости законов распределения погрешностей задания координат аппроксимируемых точек по координатным осям x и у.

В настоящей работе полученные в [5] результаты обобщены на случай независимых распределений погрешностей задания положений аппроксимируемых точек по радиальным и угловым координатам. Приведены аналитические выражения для смещений МНК2-оценок координат центра и радиуса окружности в зависимости от погрешностей

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

1. Анализ смещений оценок параметров аппроксимирующей окружности

Критерий оптимизации в рассматриваемом методе МНК2 записывается в следующем виде:

— Хо)2 + (уз- Уо)2 — — хо)2 + (у — Уо)2)) >> (1)

Л = — хо)2 + (у* — уо)2), (2)

где (хг, уг) — координаты исходных точек (г = 1. , N), Л2 — средний квадрат радиуса аппроксимирующей окружности, (хо,уо) — координаты центра окружности. Соотношение (1) имеет аналитическое решение:

Л2 = х2 + у2 — 2(хоХ + уоУ) + Хо + у^, (3)

где черта сверху означает усреднение по выборке из N точек. Координаты центра окружности вычисляются по формулам

Хо = (&1С2 — С1&2)/( , а угловые скобки обозначают усреднение по множеству реализаций при фиксированном значении i.

Величины смещений ДЯ2 и ДЯ для среднего квадрата Я2 и радиуса Я соответственно определяются из выражений (2)-(5) при усреднении по множеству N-точечных выборок:

ДЯ2 = -1 = a2 — 2 + , (7)

где параметр Я в (8) определяется из выражения (3).

Нелинейная зависимость соотношений (4) от Дi не позволяет выполнить усреднение непосредственно в выражениях (7), (8). Для решения этой проблемы воспользуемся разложениями xo и yo в ряды по степеням Д^

A 1 / „ „ „ Bi „ Bi „ B2 „ /Вл2

xo = 2B K 2B^(Ao + Ai — AoB0 + a2 — AiB0 — AoB0 + Ч^) + •••), (9)

У0 = 2B M 2B0(Po + pi — p0f + p2 — piB — p0Bo + Mf )’2 + . . .), (10)

где в соответствии с (4) использованы разложения

A = (biC2 — cib2)/2N2 « Ao + A^) + А2(Д2) + .

В = (ai&2 — b1)/2N2 « Bo + Bi(^) + В2(Д2) + . (11)

P = (bici — aiC2)/2N2 « Po + Pi(^) + Р2(Д2) + • • •

Из выражений (5), (9)-( 11) можно получить следующие результаты для средних значений:

/a2 = cos di/C — 2(1 — 2cos2 9i + cos 9i )(cos3 9i — cos 9i cos2 9^/sin2 9i/C2, (12)

=0, = a2/(CN), = a2/(Nsin2 9i), (13)

= a2/N + cos 9i; = a2/N, (14)

где введено обозначение

C = cos2 9i — cos 9i •

Соответствующее (10) среднее значение = 0 обусловлено выбранной симметричной схемой задания угловых координат исходных точек: sin 9i = 0.

Отметим, что при получении выражений (12)-(14) было использовано соотношение

справедливое для любых зависящих от 9 функций /1(й) и /2(9).

Подставляя (12)—(14) в (7), можно получить требуемое аналитическое выражение для величины смещения оценки среднего квадрата радиуса окружности ДЯ2:

ДЯ2 = a2[1 — 4/N + 1/(NC) + 1 /sin2 0¿/N] — 2cos 9i • (15)

Величина смещения оценки радиуса АЛ (8) может быть получена аналогичным образом, при разложении выражения для радиуса окружности Л = J1 + а в ряд по степеням А*, где

а = 2А + А2 — 2(xoX + yoy) + x0 + y2.

Можно показать, что

АЛ = а2/2[1/2 — 4/N + (1 — cos 0* (2 — cos2 0*))/sin2 0*/С] — cos 0* . (16)

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

2. Результаты моделирования

Для проверки справедливости полученных аналитических соотношений, исследования статистических характеристик оценок и их зависимостей от вариаций различных параметров проводилось численное моделирование МНК2. При моделировании для получения статистических характеристик использовалось 106 выборок исходных аппроксимируемых точек, координаты которых задавались, в соответствии с соотношениями (5),

с аддитивным нормально распределенным шумом N(0, а2) при нулевом среднем и дис-2

На рис. 1 и 2 представлены расчетные данные, полученные в результате численного моделирования и вычисленные согласно полученным выше аналитическим оценкам, для небольших значений N и Л = 20 мм, а = 0.3 мм. Отметим, что характер изменения приведенных зависимостей при вариациях Л и а близок к представленному нами в работе [5]. Абсолютные значения оцениваемых параметров, использовавшиеся при расчетах, получаются из аналитических оценок для относительных величин путем следующих замен: АЛ ^ ЛАЛ, АЛ2 ^ Л2 А Л2, x0 ^ Лх0, y0 ^ Лу0, а ^ Ла. Поведение показанных на этих рисунках зависимостей при больших значениях N не представляет интереса из-за значительных величин смещений оцениваемых параметров, присущих рассматриваемому методу аппроксимации. Для удобства сравнения оценок (7) и (8) на рис. 2 приведены зависимости для параметра АЛ1 = V — Л.

Из рис. 1 следует, что при 00 + )/2 — Л, где нижние индексы соответствуют четырех- и пятиточечной схеме аппроксимации. Из рис. 3 следует, что комбинированная схема обеспечивает получение практически несмещенных оценок для квадрата радиуса окружности при различных значениях 00. Незначительные для практических приложений отклонения расчетных данных от аналитических оценок для малых значений 00 обусловлены возрастанием дисперсии и погрешностей численных вычислений при уменьшении центрального угла аппроксимирующей дуги 00.

Рис. 1. Зависимости смещения АЯ оценки радиуса окружности от количества аппроксимируемых точек N при различных значениях центрального угла дуги 0о: сплошные линии — аналитические оценки, символы — результаты моделирования

Рис. 2. Зависимости параметра АЯ1 = у/ -Я от количества аппроксимируемых точек N при различных значениях центрального угла дуги сплошные линии — аналитические оценки, символы — результаты моделирования

Рис. 3. Зависимости смещения AR45 при комбинированной схеме оценки радиуса окружности от центрального угла дуги во: сплошная линия — аналитические оценки, символы — результаты моделирования

Расчеты показывают, что введение дополнительного шума в угловое распределение точек со среднеквадратичным отклонением (СКО) в несколько градусов слабо сказывается на полученных результатах.

При выводе выражений для оценок параметров предполагалось симметричное относительно оси x расположение аппроксимируемых точек, при этом не накладывалось никаких ограничений на угловые координаты точек. Полученные соотношения поэтому могут быть использованы для исследования зависимости оцениваемых параметров от вариаций углового положения точек. Наиболее интересными представляются случаи с малым количеством аппроксимируемых точек N = 4 и N = 5, обеспечивающими достаточно малое смещение параметров при значениях в0 20°. Для N = 5 смещение оценки радиуса при ^ > 20° прак-

Рис. 4. Зависимости смещения АЯ и СКО-оценки радиуса окружности от углового отклонения р внутренних аппроксимируемых точек от крайних при N = 4, N = 5 и во = 90°: символы — результаты моделирования; сплошные и штриховые линии для АЯ — аналитические оценки

Рис. 5. Зависимости смещения АЯ45 при комбинированной схеме оценки радиуса окружности от углового отклонения р внутренних аппроксимируемых точек от крайних при одинаковых значениях р для N = 4 и N = 5: сплошная линия — аналитические оценки, символы — результаты моделирования

тически постоянно и по абсолютной величине существенно превосходит аналогичный показатель для N = 4.

Зависимости смещения ДЯ45 от р показаны на рис. 5. Отметим, что на этом рисунке, в отличие от рис. 3, угловые положения внутренних точек для N = 4 и N = 5 одинаковы, и аппроксимируемые выборки точек для различных N отличаются лишь центральной точкой при в = 0. Для во = 90° минимальное смещение достигается в окрестности значения р = 30°, соответствующего равномерному угловому распределению точек при N = 4.

Расчеты показывают, что наличие шумовой составляющей в несколько градусов в угловом положении точек слабо влияет на характер представленных на рис. 4 и 5 зависимостей. Значительные изменения имеют место лишь для N = 4 при близком расположении внутренних и крайних точек, соответствующем малым значениям

Несмещенные оценки ДЯ45 в комбинированной схеме могут быть получены для различных сочетаний значений углов р = и р = соответствующих N = 4 и N = 5, что поясняется рис. 6. На этом рисунке показаны зависимости и от р для N = 4 и N = 5, из которых следует, что компенсация смещения оценки достигается для любой пары значений углов соответствующих одинаковым значениям ДЯ; = = из заштрихованной области. В частности, значению ДЯ; = 0.8 соответствуют два угла = 34.7° и = 17.4°. Заметим, что при одина-

Рис. 6. Зависимости смещений и при комбинированной схеме оценки радиуса окружности от углового отклонения р внутренних аппроксимируемых точек от крайних для N = 4 и N = 5: сплошные линии — аналитические оценки, символы — результаты моделирования

ковом угловом расположении внутренних точек р = р5 = 17.4°, как следует из рис. 5, смещение ДЯ45 существенно превосходит нулевое значение и составляет около 0.02 мм.

Приведенные выше результаты соответствуют независимому распределению шума координат аппроксимируемых точек по Я и 0, что подтверждает общие закономерности, полученные нами ранее для независимого закона распределения шума по х- и у-координатам [5]. Выводы относительно величин смещений оценок параметров окружности, обеспечиваемые для малых значений N при МНК2-аппроксимации, остаются справедливыми и в случае зависимости дисперсии шума по координатам х и у, отличной от (6). Интерес представляют, в частности, задачи, связанные с обработкой изображений и измерением координат точек границ объектов (контуров) вдоль скан-линий, соответствующих анализируемым сечениям изображения х = хс или у = ус. Расчеты показывают, что, в предположении пропорциональной зависимости дисперсии измеряемой координаты границы от ширины контура изображения в анализируемом сечении, величины смещений и дисперсии МНК2-оценок для малых N существенно возрастают лишь с приближением скан-линий к положениям касательных к границе.

Полученные аналитические оценки хорошо совпадают с результатами численного моделирования и подтверждают сделанные нами ранее выводы об эффективности применения неитерационного МНК-оценивания и четырехточечной схемы аппроксимации для оценки радиуса и координат центра аппроксимирующей окружности [5]. Метод аппроксимации МНК2 при N = 4 обеспечивает минимальные смещения оценок, близкие к нулевым значениям, и может применяться в условиях ограниченных вычислительных ресурсов для аппроксимации дуги окружности. Справедливость данного положения подтверждается результатами численного моделирования различных практически важных координатных зависимостей дисперсии шума координат аппроксимируемых точек.

Приведенные в работе аналитические оценки могут использоваться для исследования зависимости оцениваемых параметров от углового положения точек и оптимизации схемы аппроксимации.

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

Автор выражает признательность В.П. Косых за полезные дискуссии.

[1] Баранов В.Г., Ильясов Н.Ю., Котляр Р.В. и др. Бесконтактное измерение радиуса кривизны сферических поверхностей // Распознавание образов и анализ изображений: новые информационные технологии (РОАИ-5-2000). Тр. Междунар. конф. Самара, 2000. Т. 3. С. 458-462.

[2] Гривов М.Г., ХАЧУМОВ В.М. Определение геометрических параметров объектов по растровым изображениям // Автометрия. 2001. № 1. С. 40-49.

[3] Demeyere M., Dereine E., Eugene C., Naydenoy V. Measurement of cylindrical objects through laser telemetry: application to a new forest caliper // IEEE Trans. Instrument. Measur. 2002. Vol. 51, N. 4. P. 645-649.

[4] Mikhlyaey S.V. Method for measuring the diameter of a growing crystal // Pattern Recogn. Image Analysis. 2005. Vol. 15, N. 4. P. 690-693.

[5] МихляЕВ С.В. Аппроксимация окружности при измерении диаметра кристалла // Вы-числ. технологии. 2007. Т. 12, № 1. С. 61-71.

[6] Zhang Z. Parameter estimation techniques: a tutorial with application to conic fitting // Image Vision Comput. 1997. Vol. 15, N. 1. P. 59-76.

[7] Marquardt D. An Algorithm for least squares estimation of nonlinear parameters // SIAM J. Appl. Math. 1963. Vol. 11. P. 431-441.

[8] Landau U.M. Estimation of a circular arc center and its radius // Comput. Vision, Graphics Image Proces. 1987. Vol. 38. P. 317-326.

[9] Spath H. Least-squares fitting by circles // Comput. 1996. Vol. 57, N. 2. P. 179-185.

[10] Chernov N., Lesort C. Least squares fitting of circles //J. Math. Imaging Vision. 2005. Vol. 23. P. 239-251.

🎥 Видео