CC BY
277
УДК 519.6
В.В. ВИШНЕВСКИЙ, И.К. РЫСЦОВ, М.В. ВОЛЖЕВА
ИТЕРАЦИОННЫЙ АЛГОРИТМ ПОСТРОЕНИЯ КРИВОЙ БЕЗЬЕ ПО ЗАДАННЫМ ТОЧКАМ
Abstract: In article it is considered the algorithm of aproximation of discrete date with help of parametric splines in the manner of Bezier -curves. The hypothesis is advanced about possibility of using Beziers spline in problems of pattern recognition.
Key words: Bezier splines, aproximation, iteration algorithm, pattern recognition.
Анотація: У статті розглядається алгоритм апроксимації дискретних даних за допомогою
параметричного сплайну у вигляді кривої Без'є. Висувається гіпотеза про можливість використання кривих Без'є в задачах розпізнавання образів.
Ключові слова: крива Без'є, апроксимація, ітераційний алгоритм, розпізнавання образів.
Аннотация: В статье рассматривается алгоритм аппроксимации дискретных данных параметрическим сплайном в виде кривой Безье. Выдвигается гипотеза о возможности использования кривых Безье в задачах распознавания образов.
Ключевые слова: кривая Безье, аппроксимация, итерационный алгоритм, распознавание образов.
1. Введение
Известно, что кривые Безье являются основой многих компьютерных приложений для компьютерной графики. Так, например, для CAD/CAM/CAE - систем высшего и среднего уровня, а также для коммерческих, внутренних производственных приложений и приложений для обмена данными промышленным стандартом становится геометрическое ядро Parasolid компании Unigraphics Solutions Inc., в основу которого, в свою очередь, положена теория неоднородных рациональных В-сплайнов (NURBS) и кривых Безье [1].
Интуитивно понятная геометрическая интерпретация кривых Безье позволяет предположить, что они могут быть эффективно использованы и в задачах распознавания образов. Одной из таких задач может быть, например, распознавание бинарных растровых изображений. Для этой задачи использование кривых Безье может обеспечить инвариантность алгоритмов распознавания при осуществлении афинных преобразований. Однако широкое использование кривых Безье в задачах распознавания образов сдерживается сложностью вычислительных алгоритмов аппроксимации экспериментальных данных параметрическими сплайнами, к которым относятся и кривые Безье.
В данной публикации описывается итерационный алгоритм, обладающий средней вычислительной сложностью, что позволяет аппроксимировать, например, участки контура бинарного изображения или других экспериментальных данных, которые могут быть описаны гладкой кривой с одним или двумя экстремумами.
2. Математическая постановка задачи
Обычно под кривой Безье понимается дуга плоской кривой третьего порядка Bz(t ) = (Bx(t), By(t)),
О < t < І, заданной в следующем параметрическом виде:
3
3
(1)
где БеГу (?) - базовые скалярные полиномы Бернштейна третьей степени: = (х(0;.),уОу)) -
коэффициенты кривой, а параметр ? изменяется в единичном интервале [0,1]. Полиномы
Бернштейна определяются следующим образом:
, . 3! ,
Бег (?) = С3 • (1 - 0 1 г1 =-——- • (1 - 0 и], 0 < у < 3, (2)
7 у!(3-у)!
где С3 , 0 < < 3 - биномиальные коэффициенты. Базис Бернштейна является довольно
распространенной формой аппроксимации функций, поскольку он имеет ряд преимуществ перед классическим базисом {1, ?, ?2, ?3}. В частности, он обладает большей вычислительной
устойчивостью [1].
Кривая Безье задается в параметрическом виде, поскольку в общем случае геометрические контуры не могут быть описаны в виде однозначной функции У = /( х). Например, контур может иметь вертикальные касательные или описываться обратной функцией х = /(у) и т. д. В этом состоит одна из трудностей непосредственного применения классических численных методов полиномиальной аппроксимации [2].
Заметим, что коэффициенты кривой Безье О0 и О3 определяют крайние точки, через которые проходит кривая
Бг(0) = (Бх(0), Бу(0)) = Оо , Бг(1) = (Бх(1), Бу (1)) = 03 . (3)
Коэффициенты 01 и 02 определяют величины касательных в крайних точках («усы» кривой). На
рис. 1 показан пример кривой Безье, из которого становится понятным ее геометрический смысл.
О
3
Рис. 1. Пример кривой Безье Учитывая приведенные выше определения, можно сформулировать задачу следующим образом. Пусть имеется упорядоченная последовательность различных точек на плоскости Р = (х(р),у(р)), 0 <, < N + 1, содержащая, по крайней мере, четыре точки. Требуется провести
кривую Безье таким образом, чтобы она проходила через две крайние точки Р0, PN+1 и
располагалась как можно ближе к промежуточным точкам Р,, 1 <, < N.
В качестве меры отклонения кривой от заданных точек используется квадратичная функция невязки:
1 N
/ (О, и) = - X [(Бх (г,) - х(Р, ))2+ (Бу(1,) - у(Р, ))2]. (4)
2 ,=1
Здесь суммируются квадраты евклидовых расстояний между точками Бх(г() = (Бх(г,),Бу(г,)) и Р = (х(р), у (р.)). Вектор коэффициентов кривой 0 = {00, 01, 02, 03} и параметрический вектор и = (tl,...tN) , где все принадлежат единичному интервалу [0, 1], являются аргументами этой функции. Каждую точку Бх(г,) можно рассматривать как «представителя» на кривой Безье, соответствующего исходной точке Р. Поэтому функцию невязки нужно минимизировать по 0 и
по и .
Поскольку функция невязки нелинейно зависит от параметрического вектора, то в целом эту задачу можно рассматривать как нелинейную проблему наименьших квадратов и решать ее известными численными методами [3, гл. 10]. Но реализация этих методов может потребовать в данном случае слишком большого объема вычислений, поскольку при этом придется вычислять многомерные якобианы и гессианы. Поэтому, учитывая возможную размерность задачи (от нескольких десятков до сотен точек), авторами был разработан другой метод решения.
Далее будет показано, что компоненты базовых функций
х(^ ) = Бх(г,) - х(р ) , у(^ ) = Бу(г,) - у(р,) , 1 < , < N (5)
линейно зависят от коэффициентов кривой Безье, и, очевидно, что каждая из них зависит только от одной переменной вектора и . Из (4) видно, что функция невязки, по существу, является суммой
квадратов этих компонент. Поэтому данную задачу целесообразно разбить на несколько подзадач. Сначала можно решить линейную задачу наименьших квадратов и провести минимизацию целевой функции по первому аргументу 0 при фиксированном значении второго аргумента и, а затем, наоборот, зафиксировать 0 и провести минимизацию по каждой компоненте вектора и.
Наконец, на заключительном этапе, комбинируя эти методы, можно получить итерационный алгоритм минимизации целевой функции. В соответствии с этим планом и построим наше дальнейшее изложение.
3. Линейная задача наименьших квадратов
Итак, зафиксируем параметрический вектор и = (tl,...tN ) и будем искать кривую Безье, которая минимизирует значение невязки (4) при заданных значениях промежуточных параметров.
Согласно свойству (3), коэффициенты 00, и 03 должны совпадать с крайними точками, через которые проходит кривая Безье, но по условию задачи этими крайними точками должны быть точки Р0 и Р1^+1. Поэтому для граничных значений параметра г0 = 0 и tN+1 = 1 получаем следующие краевые условия:
Бх(0) = 00 = Р0, Б1 (1) = 03 = PN+1. (6)
Остается найти коэффициенты Q1 и Q2 кривой Безье. Для этого удобно перейти к
матричной записи основных соотношений между параметрами задачи. Обозначим через Вег(и)
матрицу размера N X 4 , состоящую из значений полиномов Бернштейна [Вег^ (ї,)] ,
1 < і < N , 0 £ j < 3 , вычисленных на параметрическом векторе и . Эту матрицу назовем бернштейнианом (матрицей Бернштейна) и она для полиномов Бернштейна играет ту же роль, что и матрица Вандермонда для обычного полиномиального базиса [1]. Столбцы матрицы Бернштейна
будем обозначать через ВеГу (и ) , 0 < j < 3.
Далее пусть х(Р) = (х(Р1),... , х(PN )) - вектор абсцисс точек Рі , 1 < і < N , а
У (Р) = (У (РіХ••• , У (PN )) - вектор ординат этих точек. Аналогично пусть
х^) = (x(Q0), x(Q1), х^2), x(Q3)) - вектор абсцисс коэффициентов кривой Безье и
соответственно У^) - вектор ординат этих коэффициентов. Наконец, обозначим через х Р)
вектор абсцисс базовых функций (х(Р1),. , х(FN)) , определенных по формуле (5), и
аналогично через у(р) обозначим вектор ординат этих функций. Тогда из (1) и (5) получаем
следующие матричные равенства:
х( Р) = Вег (и) • х^) - х( Р), у(Р) = Вег (и) • y(Q) - у( Р). (7)
Функцию невязки (4) можно также представить в векторной форме:
1 2 1 2 /(Q, и ) = 2||х(Р )||2 + 2\\У (Р )||2 , (8)
где норма вектора берется в ^-мерном евклидовом пространстве [4].
Из формулы (7) видно, что вектор абсцисс базовых функций х Р) линейно зависит от
вектора абсцисс х (в) и аналогично вектор ординат базовых функций у (Р) линейно зависит от
вектора у^ ). Таким образом, получаем классическую линейную задачу о наименьших квадратах
[3]. Хорошо известен метод Гаусса для решения этой задачи, а именно, векторы х (й) и У(Я), на которых достигается минимум невязки (8), должны удовлетворять следующей системе линейных уравнений [3, § 3.6]:
ДЦ) • х(О) = Вег(Ц)Т • х(Р) , Л(Ц) • у(О) = Вег(Ц)Т • У(Р) , (9)
где Л(и) = Вег (и )Т Вег (и) , а Вег (и )Т - транспонированный бернштейниан. Чтобы убедиться в этом, достаточно продифференцировать невязку (8) по переменным х (Q) и У(Q) и приравнять все частные производные нулю, т.е. записать необходимое условие экстремума. Матрица Л (и) = [а. ] будет составлена из скалярных произведений столбцов матрицы Бернштейна:
N
а. =(Вег1 (и),Вег,.([/)) = С3 • С? • £ • (1 -1,)6-'-‘, 0 < j,к < 3 . (10)
і=1
В нашем случае система (9) будет содержать по четыре уравнения с четырьмя неизвестными. Если матрица A(U) будет обратимой, то система (9) будет иметь единственное решение, которое дается следующими формулами:
xQ) = А(и)-1 Вег(Ц)тх(Р), у@) = АЦ)-1 ВегЦ у(Р). (11)
На этом, в принципе, решение линейной задачи можно считать законченным, но возникает
вопрос об условиях, когда матрица A(U) будет обратимой. На этот вопрос ответ будет дан чуть
позднее, а сейчас обратим внимание на еще одну особенность, которая отличает наше решение от общего решения (11). Из краевых условий (6) имеем выражения для крайних коэффициентов
Q0 = P0 и Q3 = PN+1, и, следовательно, в системе (9) в нашем случае будет всего по два
неизвестных вместо четырех. Чтобы найти компактные выражения для этих неизвестных,
воспользуемся аппаратом векторной алгебры. При желании в дальнейшем читатель может рассматривать неизвестные Q1, Q2 и точки р как комплексные числа. Подставим краевые
условия (6) в уравнение (9). Тогда, раскрывая скобки и собирая все члены при Q1 и Q2, получим следующую систему уравнений относительно этих неизвестных:
а11 ' й + а12 ' Q2 = ^ , а21 " Q1 + a22 " Q2 = R2 , (12)
где свободные члены Rj вычисляются по следующим формулам:
N
^ = £ Вег}(^) •p - ао ] • p0 - aJ 3 • р+1 , 1 < ] < 2 . (13)
г=1
Здесь и в (12) коэффициенты aJ■k являются элементами матрицы A(U) и определяются по формуле (10).
Обозначим определитель системы (12) через А и предположим, что он не равен нулю. Тогда по формулам Крамера получаем решение системы (12) в следующем виде:
Й = (а22 • ^ — ^2 • ) / А , Q2 = (а11 • R2 — а12 • ^) ^ А . (14)
Из формулы (10) заключаем, что определитель А = a11a22 — a12a21 системы (12)
определяется по следующей формуле:
А = ЦВегДЦ/)||2 ЛBer2 (IГ)||2 —(Вег (U), Ber2 (IГ))2. (15)
Здесь норма и скалярное произведение по-прежнему определяются в N -мерном евклидовом пространстве. Тогда из известного неравенства Коши - Буняковского [5, § 1.3-2] заключаем, что определитель А должен быть неотрицательным. Он будет положительным тогда и только тогда, когда столбцы Ber] (Ц) и Ber2 (Ц) являются линейно-независимыми. Этот вывод является частным случаем известной теоремы [3], которая утверждает, что матрица А (и) будет обратимой, если столбцы матрицы Вег (и) являются линейно-независимыми.
Поскольку Вег1^1)/Вег2(^.) =1/ti — 1, то для положительности определителя А необходимо и достаточно, чтобы внутри открытого интервала [0,1] находились по крайней мере два различных значения t1 и t2 из параметрического вектора Ц. Вместе с граничными значениями 0 и 1 это дает четыре различных значения параметра. На практике, как правило, выполняется более сильное условие, поскольку вектор и предполагается упорядоченным:
0 = to < ^ < ••• < tN < tN + 1 = 1 . (16)
Поскольку по условию задачи N > 2, то на любом упорядоченном параметрическом векторе определитель А будет положительным.
4. Минимизация по параметрам
Рассмотрим теперь проблему минимизации невязки (4) для фиксированной кривой Безье, т.е. при фиксированном векторе Q. Тогда невязка будет зависеть только от параметрического вектора
и = (U1,...UN) и, следовательно, в этом случае имеем нелинейную задачу наименьших квадратов.
Если параметры ti, 1 < г < N, рассматривать как независимые действительные
переменные, то существование глобального минимума по всем переменным ti следует из известной теоремы анализа о том, что непрерывная функция / (и) достигает своего минимума на
компакте (многомерном кубе) [0,1^ . Свойство (16) упорядоченности вектора Ц противоречит предположению о независимости его компонент. Но нетрудно видеть, что упорядоченные параметрические векторы образуют открытое всюду плотное подмножество куба [0,1]N , поэтому к указанному минимуму можно приблизиться с любой точностью, двигаясь только по упорядоченным параметрическим векторам.
Сначала запишем необходимое условие экстремума. Для этого продифференцируем невязку (4) по всем переменным ti и в результате получим следующую систему нелинейных уравнений:
= ВХ(^ ) • (Вх (^ ) — х (р. )) + ВУ '(^ ) • (ВУ (^ ) — у (р )) = 0 , 1 < I < N . (1 7)
Для вычисления производных Вх'^1) и Ву (Х1) можно воспользоваться формулой (1) и взять обычную производную от функций Вх(Х) и Ву^) по переменной t, а затем подставить вместо нее параметр ti. Производные от полиномов Бернштейна Вег'(ti) также легко вычисляются по формуле (2) и выражаются через полиномы Бернштейна второй степени [1]. Из системы (17), раскрывая скобки, получаем следующую систему из N нелинейных уравнений с N неизвестными ti:
Вх'(Г) • Вх(Л) + Ву,(г.) • Ву^.) = Вх'^.) • х(Р) + Ву'(Г) • у (р.) , 1 < г < N . (18)
Уравнения этой системы зависят от различных переменных ti, 1 <. < N , поэтому вместо
решения многомерной нелинейной системы здесь нужно N раз решать уравнения от одной переменной. Из формул (1) и (2) видно, что левая часть каждого из уравнений (18) представляет собой, по существу, один и тот же многочлен пятой степени М(^ относительно переменной t, в который вместо этой переменной каждый раз подставляется своя переменная ti. Правая часть каждого уравнения системы (18) является многочленом второй степени М( (^.) относительно той же переменной ti. Таким образом, задача сводится к поиску вещественных корней у многочленов
пятой степени, и, значит, каждое уравнение М (^.) = М1 (^.) этой системы может иметь не более 5
и не менее одного вещественных корней, из которых можно выбрать один корень, который минимизирует невязку (4).
Но локализация и нахождение корней многочлена пятой степени является не простой вычислительной задачей. Кроме того, экстремальная точка может оказаться локальным максимумом или точкой перегиба функции невязки. К тому же, точное решение промежуточной задачи может существенно увеличить общее время минимизации невязки по всем переменным. Поэтому вместо точного решения системы (18) был использован приближенный локальный алгоритм минимизации невязки градиентного типа.
Этот алгоритм основан на геометрических соображениях и состоит в изменении значения переменной ti таким образом, чтобы уменьшалась невязка (4). С геометрической точки зрения
соотношения (17) означают перпендикулярность касательного вектора Вг' ^) = (Вх' ^), Ву (I)) в
точке t = ti вектору Вг(^.) — р = (Вх(^.) — х(Р.), Ву(^.) — у(р)) . При этом, если косинус угла
между этими плоскими векторами больше нуля, то параметр ti надо уменьшать, чтобы увеличить
угол между векторами Вг^) и Вг(^.) — р. Если же С08(ц) < 0, то параметр ti надо
увеличивать, чтобы уменьшить угол между этими векторами. Другими словами, смещение должно происходить как бы в сторону перпендикуляра, опущенного из точки ^ на кривую Безье. На рис. 2 иллюстрируются возникающие здесь ситуации. Конечно, вариацию параметра ti надо проводить
только в том случае, если невязка уменьшается, поскольку подозрительная точка может оказаться локальным максимумом. Смещение можно осуществлять методом половинного деления между соседними значениями параметра либо другими подобными методами. Метод половинного деления отрезка имеет преимущество перед смещением по равномерным узлам сетки, поскольку он позволяет сохранить упорядоченность параметрического вектора. Смещения заканчиваются, когда невязка существенно не изменяется.
Этот простой локальный алгоритм оказался достаточно быстрым. Хотя он и не гарантирует достижения глобального минимума, о котором шла речь в начале этого раздела, но позволяет быстро найти приемлемое решение.
Рис. 2. Исходные точки и кривая Безье 5. Моделирование алгоритма в пакете Мар1е
Для численного решения линейной задачи наименьших квадратов была разработана процедура Ивдиаге(и), которая по заданному параметрическому вектору и вычисляет вектор коэффициентов кривой Безье Q. Аналогично, для реализации алгоритма минимизации по параметрам была разработана процедура МШроШэ(0,и), которая по заданной кривой Безье Q и параметрическому вектору и смещает вектор и по всем параметрам ti, 1 < / < N с целью уменьшения невязки (4).
Общий алгоритм минимизации начинает работу с равномерного параметрического вектора ti = / /{Ы +1), 1 < i < N, и состоит в поочередном вызове процедур ИЭдиаге(и) и МШроМэ(0,и)
до тех пор, пока изменение невязки станет меньше заданного числа е > 0.
На рис. 3 приведены примеры работы алгоритма, характерные для задачи аппроксимации участка контура бинарного изображения, для первых 3-х итераций. Крестиками обозначены исходные значения данных, через которые должна пройти искомая кривая Безье.
Bezie Curve with points Bezie Curve with points Bezie Curve with points
n -in on "ЗЛ ЛП СП П 1П on . Tl ЛП e;n
Рис. 3. Пример моделирования итерационного алгоритма
После проведения моделирования алгоритма была разработана динамическая библиотека (dll), которая, в свою очередь, была успешно использована в эксперименте распознавания бинарных изображений нейрокомпьютером. Эксперимент подтвердил гипотезу о возможности построения алгоритмов распознавания, инвариантных к афинным преобразованиям. Описание этого эксперимента и связанной с ним информационной технологии обработки бинарных изображений является темой отдельной статьи.
6. Выводы
Задача аппроксимации участков контура кривыми Безье может быть эффективно решена на практике путем разбиения на две подзадачи: линейную процедуру аппроксимации по методу наименьших квадратов и нелинейную итерационную процедуру подбора параметров.
Итерационная процедура подбора параметров кривой Безье не обеспечивает нахождения глобального экстремума для любого типа кривой Безье. Однако для практических задач аппроксимации контуров бинарных изображений такая процедура позволяет минимизировать невязку (4) с любой вычислительной точностью.
Учитывая интуитивно понятный геометрический смысл кривой Безье, приемлемую вычислительную сложность алгоритма ее построения и проведенные модельные эксперименты, можно выдвинуть гипотезу о возможности использования кривых Безье в задачах распознавания образов.
СПИСОК ЛИТЕРАТУРЫ
1. Денискин Ю.И. Особенности аппроксимации обводов параметрическими полиномами в форме Бернштейна // Прикладная геометрия. - 1999. - Вып. 2, № 2.
2. Хемминг Р.В. Численные методы. - М., Наука, 1972. - 400 с.
3. Денис Дж., Шнабель Р. Численные методы безусловной оптимизации и решение нелинейных уравнений. -М., Мир, 1988. - 440 с.
4. Халмош П. Конечномерные векторные пространства. - М.: ФМ, 1963. - 263 с.
5. Корн Г., Корн Т. Справочник по математике. - М., Наука, 1970.- 720 с.
