Метод итеративной компенсации проективных искажений изображений
ИКОНИКА – НАУКА ОБ ИЗОБРАЖЕНИИ
УДК 621.397.3: 536.3
МЕТОД ИТЕРАТИВНОЙ КОМПЕНСАЦИИ ПРОЕКТИВНЫХ ИСКАЖЕНИЙ ИЗОБРАЖЕНИЙ
© 2009 г. В. Р. Луцив, канд. техн. наук
НПК “Государственный оптический институт им. С.И. Вавилова”, Санкт-Петербург
E-mail: vluciv@list.ru
Представлен метод нормализации изображения относительно полной группы проективных преобразований на плоскости, основанный на итеративной оптимизации решения в пространстве двух параметров проективного преобразования, отличающих его от аффинного. Эта пара параметров представляется вектором. Каждой итерации нормализации предшествует измерение параметров и компенсация аффинного преобразования согласно формулам, выведенным и опубликованным автором ранее. Для частичной компенсации проективного преобразования на каждой итерации используется такой по направлению вектор нормализующего преобразования, что изменению его направления на противоположное соответствует наибольшее по модулю изменение скорости смещения центра тяжести преобразуемого изображения. Итеративный процесс сходится к состоянию изображения, эталонному относительно проективного преобразования, при котором вектор нормализующего преобразования любого направления приводит к одинаковому по модулю смещению центра тяжести.
Ключевые слова: афинные преобразования, проективные преобразования, смещение центра тяжести, вектор нормализующего преобразования
Коды OCIS: 100.2980, 3000.3020
Поступила в редакцию 27.02.2009
Введение
Измерение параметров и компенсация аффинных преобразований (АП) и проективных преобразований (ПП) и поиск инвариантов этих преобразований являются важными направлениями автоматического анализа изображений. Это связано с тем, что в большинстве систем компьютерного зрения изображение на матрице фоточувствительных элементов формируется путем проецирования через сферическую линзу. Вследствие этого наблюдаемые изображения двумерных проекций объектов претерпевают на плоской матрице фотоприемников, а затем и в памяти компьютера геометрические трансформации, описываемые проективной группой преобразований на плоскости. Подгруппой ПП является группа АП, соответствующая случаю, когда объектив системы компьютерного зрения находится достаточно далеко от наблюдаемой сцены. Известны методы [1, 2] вычисления инвариантов группы АП на плоскости. Известен также альтернативный подход [3], в соответствии с которым сначала могут быть в аналитической фор-
ме определены и измерены параметры АП, а затем это преобразование может быть компенсировано, причем при измерении параметров преобразования информация о типе анализируемого объекта (об эталоне) не используется. ПП в отличие от АП нелинейно. Это создает препятствия как для вычисления инвариантов ПП, так и для построения аналитических выражений, пригодных для определения параметров ПП по содержанию изображения. В статье предлагается метод итеративной нормализации изображений относительно их ПП. Применяемый подход является расширением подхода к компенсации АП, предложенного автором в работах [3, 4]: в настоящей работе в результате нормализации изображение приводится к состоянию, эталонному относительно ПП, и измеряются параметры этого преобразования.
1. Измерение параметров и компенсация АП
В работах [3, 4] показано, что АП |A| может быть разложено на простые геометрические трансформации следующим образом:
“Оптический журнал”, 76, 7, 2009
53
a1 a2 a3 A = a4 a5 a6 =
001
= A1 ⋅ A2 ⋅ A3 ⋅ A4 ⋅ A5 ⋅ A6 ,
(1)
где a1–a6 – параметры АП, представленного в матричной форме, а матрицы | Ai| описывают в том порядке, в котором они пронумерованы в
выражении (1), смещение центра тяжести (ЦТ)
по двум декартовым координатным осям, вытя-
гивание или сжатие в каком-либо направлении,
зеркальное отражение относительно прямой,
проходящей через начало координат, масшта-
бирование и вращение относительно начала
координат.
Представленный в выражении (1) порядок
простых преобразований может быть изменен,
однако, именно такой порядок позволяет произ-
водить поэтапное измерение их параметров и их
компенсацию, начиная с | A1| и заканчивая | A6|, поскольку именно при компенсации простых
преобразований | Ai| в порядке от | A1| к | A6| последующие операции компенсации не нарушают
результатов предыдущих. При этом преобразова-
ния | A1| и | A2| взаимно ортогональны и могут быть переставлены местами в выражении (1). Анало-
гично могут быть переставлены местами преоб-
разования | A5| и | A6|. Преобразование | A4| может быть переставлено местами с преобразованием
| A5|, и оно может измеряться и компенсироваться либо одновременно с преобразованием | A3|, либо одновременно с преобразованием | A6|.
Итак, в качестве первого шага нормализации
декартовы координаты xc и yc ЦТ объекта вычисляются по следующим формулам:
xc = ∫∫ xB(x,y)dxdy ∫∫ B(x,y)dxdy,
ΩΩ
yc = ∫∫ yB(x,y)dxdy ∫∫ B(x,y)dxdy,
ΩΩ
(2)
где Ω – двумерная область определения изображения объекта, B(х, у) – яркость пиксела с координатами {х, у}. Это дает возможность компенсировать линейные смещения | A1| и | A2| изображения объекта.
Следующим преобразованием, подлежащим
измерению и компенсации, является масштабирование | A3| вдоль некоторого направления. В работах [3, 4] получены следующие аналити-
ческие выражения для вычисления параметров преобразования | A3| – направления Θ и величины μ изменения масштаба изображения объекта:
( ) ( )μ = ± D + C2+ E2 D − C2+ E2 , (3)
Θ =0,5 arctg(C/E)+πn,
(4)
E = ∫∫ B(r,ϕ)r2cos(2ϕ)drdϕ,
Ω
C = ∫∫ B(r,ϕ)r2sin(2ϕ)drdϕ,
Ω
(5) (6)
D = ∫∫ B(r,ϕ)r2drdϕ,
Ω
(7)
где B(r,ϕ) – яркость пиксела центрированного изображения объекта, представленного в поляр-
ной системе координат, Ω – область определения изображения объекта в полярной системе коор-
динат, {r, ϕ} – соответственно радиус и азимут пиксела в полярной системе координат.
После компенсации преобразования | A3| остаются некомпенсированными преобразования | A4|, | A5| и | A6| зеркального отражения, равномерного масштабирования и вращения. Преобразование
| A5| не зависит от углового положения пиксела и не влияет на | A6| и | A4|, поэтому параметр М
преобразования
A5
=
M 0
0 M можно измерить
по проекции изображения на радиус полярной системы координат. Коэффициент М может быть измерен согласно следующему выражению по смещению ЦТ проекции изображения на радиус относительно эталонного положения А:
M = ∫∫ rB(r,ϕ)drdϕ ∫∫ B(r,ϕ)drdϕ/A. (8)
ΩΩ
Вращение центрированного изображения относительно начала координат и его зеркальное отражение относительно прямой, проходящей через начало координат, не изменяют радиальные координаты пикселов, поэтому параметры этих преобразований можно измерить по проекции изображения на окружность, т. е. на угловую полярную координату. Параметры Ψ вращения и ε зеркального отражения могут [4] быть вычислены по формулам (9) и (10)
Ψ = ϕk − ϕm − ε(ψk − ψm ),
(9)
ε = (2ϕk − ϕn − ϕm )/(2ψk − ψn − ψm ), (10)
где ϕm, ϕk, ϕn – фазы трех соседних Фурьегармоник (m = k – 1, n = k + 1) проекции анализируемого изображения на окружность, ψm, ψk,
54 “Оптический журнал”, 76, 7, 2009
ψn – фазы аналогичных гармоник для соответствующего эталонного изображения. Возможны и другие пути вычисления параметров Ψ и ε [4], которые могут быть использованы, если в спектре проекции изображения на окружность подходящих соседних гармоник не обнаружено.
2. Итеративная компенсация и измерение параметров ПП
ПП на плоскости описывается следующими выражениями, связывающими координаты исходного X = |x, y|Tи преобразованного X′ = |x′, y′|T изображений
x′ = (a1x + a2y + a3 )/(a7x + a8y +1)= = (a1x + a2y + a3 )/(X⋅ P +1),
(11)
y′ = (a4x + a5y + a6 )/(a7x + a8y +1)= = (a4x + a5y + a6 )/(X⋅ P +1),
(12)
где a1–a6 – параметры АП, упомянутые в формуле (1), а параметры a7 и a8 характеризуют отличие ПП от АП. Из выражений (11) и (12)
видно, что знаменатель формулы ПП можно выразить через скалярное произведение X·P
вектора декартовых координат точки и вектора параметров P = |a7, a8|T. ПП – нелинейно, поэтому не удается определить вектор P его параметров в
аналитическом виде, как это было сделано для
параметров АП. Проблема решается путем ите-
ративной оптимизации решения в пространстве
параметров a7 и a8. Разработанный метод основан на том, что ПП в отличие от АП приводит к
перераспределению точек изображения объекта относительно его центра тяжести, что несложно
показать аналитически. Представим формулу (11) в следующем виде:
x′ = α/(1+ sign(β)⋅ | β |),
(13)
где α соответствует числителю выражения (11), | . | обозначает модуль, β = X·P. Пусть для определенности вектор P антипараллелен оси абсцисс, т. е. его скалярное произведение с ортом оси абсцисс отрицательно. Пусть ЦТ изображения находится в начале координат. При β ≠ 0 x′ = = α/(1 – |β|) справа от оси ординат, а слева от оси ординат x′ = α/(1 + |β|). Путем несложных преобразований легко показать, что α/(1 + |β|) – α/1 < 0, α/(1 – |β|) – α/1 > 0 и |α/(1 + |β|) – α/1| < < |α/(1 – |β|) – α/1|. Из этого следует, что абсциссы точек, находящихся левее оси ординат, в результате приложенного ПП уменьшаются по модулю,
а абсциссы точек, находящихся правее оси ординат, увеличиваются по модулю, что приводит к смещению ЦТ в область положительных абсцисс. При этом правее оси ординат приращение абсцисс происходит быстрее, а по мере удаления от начала координат увеличивается степень растяжения (разрежения) изображения вдоль оси абсцисс в области положительных абсцисс, и увеличивается степень сжатия вдоль оси абсцисс в области отрицательных абсцисс. В связи с этим после компенсации смещения ЦТ слева от оси ординат будет сосредоточено больше точек изображения (эта часть изображения будет более массивна), но в среднем они будут иметь меньшие по модулю абсциссы, а справа от оси ординат изображение будет более разреженным: оно будет содержать меньшее число точек, имеющих в то же время большие абсциссы. Указанные построения справедливы и для произвольного направления вектора P, поскольку такие случаи можно свести к рассмотренному заменой системы координат.
Рассмотрим произвольное изображение, ЦТ которого имеет нулевые координаты. Разделим изображение на два полу-изображения произвольной прямой, проходящей через начало координат. Полученные полу-изображения можно представить в виде половинок диполя, массы M1 и M2 которых сконцентрированы в ЦТ полуизображений, отстоящих от начала координат соответственно на расстояния L1 и L2. Диполь уравновешен относительно ЦТ полного изображения, но в общем случае не симметричен, т. е. M1L1 = M2L2, но M1 ≠ M2 и L1 ≠L2, причем в каком-то направлении Θ отношение L2 к L1 максимально. Тогда изображение можно представить диполем, ориентированным в направлении Θ, и описать его вектором L, длина которого равна отношению L2 к L1, а направление соответствует направлению плеча L2 диполя. Эталонное состояние изображения такое, при котором
M1 = M2 , L1 = L2,
(14)
т. е. при котором оно представляет собой диполь, симметричный и уравновешенный относительно ЦТ изображения при любом направлении прямой, делящей его на два полу-изображения. В эталонном состоянии ||L|| = 1при любом направлении Θ плеча диполя. Примером эталонного изображения может служить круг с центром в начале координат. В соответствии с вышесказанным любое ПП с параметрами a7 ≠ 0 и/или a8 ≠0, приложенное к изображению круга, вызывает перераспределение его точек относительно ЦТ, что приводит к нарушению условия (14) и к сме-
“Оптический журнал”, 76, 7, 2009
55
щению ЦТ. В случае круга при вращении вектора Р возникает изменение направления смещения ЦТ, но модуль его смещения остается неизменным. Если изображение круга подверглось ранее ПП с параметрами P с последующим перемещением его ЦТ в начало координат, его изображение будет соответствовать несимметричному диполю, для которого ||L|| > 1. Тогда вращение вектора Pt параметров дополнительно приложенного к нему “тестового” ПП приводит к различным по модулю смещениям ЦТ, причем скорость смещения максимальна при некотором направлении Θ = Θmax вектора Pt, которое для изображения круга соответствует направлению вектора P ранее приложенного к нему ПП. Чтобы понять это, заметим, что круг симметричен, поэтому достаточно рассматривать его изображение, к которому было применено ПП, как диполь. Пусть вектор диполя L совпадает по направлению с осью абсцисс, а вектор P ранее приложенного к изображению ПП имеет противоположное направление. Докажем для диполя следующую лемму:
Лемма 1: модуль смещения ЦТ несимметричного диполя под действием ПП, вектор Pt параметров которого противоположен направлению диполя больше, чем модуль смещения ЦТ под действием ПП, для которого вектор Pt параллелен диполю.
Доказательство: Пусть М1 – масса (яркость) левой части диполя, L1 – модуль абсциссы его левой части, М2 – масса его правой части, L2 – модуль ее абсциссы. M1 = NM2, L2 = NL1, N > 1. Тогда в результате дополнительного применения к такому диполю тестового ПП с вектором параметров Pt (пусть Pt = || Pt || ), направленным также как и вектор P противоположно направлению оси абсцисс, преобразованные длины L1′ и L2′ плеч диполя в соответствии с выражениями (11) и (12) будут иметь вид (15)
L1′ = L1/(1+ L1Pt ), L2′ = L2/(1− L2Pt ). (15)
Приращения ΔL1 и ΔL2 координат половинок диполя вычисляются тогда как
ΔL1 = L1 − L1/(1+ L1Pt ) = L21Pt /(1+ L1Pt ), (16)
ΔL2 = L2/(1− L2Pt ) − L2 = L22Pt /(1− L2Pt ). (17)
На основе равенств (16) и (17) модуль ΔC1 приращения абсциссы ЦТ диполя можно описать формулой (18), при выводе которой использовались подстановки M1=NM2, L2=NL1.
ΔC1
=
ΔL1M1 + ΔL2M2 M1 + M2
=
=
(1
+
NM2L21Pt (1+ N) L1Pt )(1− NL1Pt )(M1
+
M2
)
.
(18)
Аналогично выводится формула (19), пока-
зывающая модуль приращения ΔC2 абсциссы ЦТ, если направление вектора Pt изменится на противоположное.
ΔC2
=
(1−
NM2L21Pt (1+ N) L1Pt )(1+ NL1Pt )(M1
+
M2 )
.
(19)
Сравним модули ΔC1 и ΔC2 смещений ЦТ, соответствующие противоположным направлениям
Pt (в выражении (20) знак ~ означает неизвестный знак неравенства)
(1
+
NM2L21Pt (1+ N) L1Pt )(1− NL1Pt )(M1
+
M2
)
~
~
(1−
NM2L21Pt (1+ N) L1Pt )(1+ NL1Pt )(M1
+
M2
)
.
(20)
1 – NL1Pt > 0 и 1 – L1Pt > 0, что необходимо, чтобы точки диполя не перемещались в бесконечность.
Поэтому, не нарушая неравенства, можно свести
(20) к виду (21)
(1 – L1Pt)(1 + NL1Pt) ~ ~ (1 + L1Pt)(1 – NL1Pt),
(21)
откуда после выполнения умножения и приведения подобных членов получаем N ~ 1. Выше предполагалось, что N > 1, поэтому, когда вектор параметров Pt прикладываемого тестового ПП совпадает по направлению с вектором параметров ПП, имевшего место ранее, скорость смещения ЦТ больше, чем при противоположно направленном векторе параметров тестового ПП (тогда вектор Pt параллелен диполю, а тестовое ПП частично компенсирует исходное). Итак, лемма 1 доказана.
Из доказательства леммы 1 видно, что чем больше асимметрия N диполя, тем больше разница скоростей смещения ЦТ при противоположных направлениях Pt.
Полученный результат корректен для изображения круга, однако, реальные изображения характеризуются более сложными законами пространственного распределения яркости пикселов, в общем случае не симметричными. Чтобы показать справедливость утверждения,
56 “Оптический журнал”, 76, 7, 2009
доказанного в лемме 1, для реальных образов
допустим, что произвольное изображение имело
эталонный вид, соответствующий условию (14),
и что в результате некоторого проективного пре-
образования и последующего переноса его ЦТ в
начало координат это изображение в целом при-
обрело вид несимметричного уравновешенного
диполя, ориентированного вдоль оси абсцисс, у
которого правое плечо L2 в N раз длиннее левого (M1 = NM2, L2 = NL1, N > 1). Покажем, что такое произвольное изображение может быть пред-
ставлено в виде совокупности уравновешенных
несимметричных диполей, соответствующих
отдельным парам пикселов или их фрагментов,
для каждой из которых M1 = NM2, L2 = NL1 при неизменном N > 1. Для этого выполним приве-
денный ниже алгоритм.
Алгоритм 1:
1. Выбираем любой оставшийся пиксел пра-
вее оси ординат, имеющий абсциссу l и яркость
(“массу”) ml. Если такой пиксел не найден, переходим к п. 6.
2. Выбираем любой пиксел левее оси орди-
нат, имеющий абсциссу L и яркость ML, для которого
LML ≤ lml.
(22)
Если такого пиксела не найдено, переходим к
п. 4, в противном случае “уравновешиваем” от-
носительно оси ординат найденный левее оси ординат пиксел с параметрами {L, ML} “фрагментом” ранее выбранного пиксела, имевшего параметры {l, ml} и находящегося правее оси ординат. Исключаем из дальнейшего рассмотрения
полностью уравновешенный пиксел слева от оси
ординат и “уравновешивающий” его фрагмент яркостью LML/l пиксела справа.
3. Оставшаяся неуравновешенной часть
пиксела справа от оси абсцисс имеет яркость ml′ = ml – LML/l. Если новая яркость этого пиксела ml′ > 0, переходим к п. 2, в противном случае переходим к п.1.
4. Выбираем любой пиксел слева с параметрами {L, ML}, для которого LML ≥ lml. Такой пиксел найдется, поскольку неравенство (22)
в п. 2 выполнено не было. Уравновешиваем пиксел справа, имеющий параметры {l, ml}, фрагментом найденного слева пиксела. Исклю-
чаем из дальнейшего рассмотрения полностью
уравновешенный пиксел справа от оси ординат
и “уравновешивающий” его фрагмент яркостью lml/L пиксела, найденного слева.
5. Оставшаяся неуравновешенной часть пиксела слева имеет яркость M′ = M – lml/L. Выби-
раем любой пиксел справа с параметрами {l, ml}, для которого LML′ ≥ lml. Если такой пиксел не найден, переходим к п. 1, в противном случае переходим к п. 4.
6. Конец выполнения алгоритма. Каждый раз при выполнении п. 2 и п. 4 из дальнейшего рассмотрения исключаются уравновешенные относительно оси ординат диполи, соответствующие парам пикселов или их частей. Исключаемые диполи имеют соотношение плеч N в направлении исходного диполя, который до начала выполнения настоящего алгоритма соответствовал всему изображению. Не трудно убедиться, что часть исходного диполя, не исключенная еще из рассмотрения после каждого выполнения п. 2 или п. 4, также уравновешена относительно оси ординат и имеет аналогичное соотношение плеч N в направлении исходного диполя. Поэтому, если при выполнении п. 1 оказалось, что правее оси ординат не найдено ни одного пиксела или фрагмента пиксела, еще не исключенного из рассмотрения, можно утверждать, что и слева от оси ординат не осталось ни одного подлежащего рассмотрению пиксела или фрагмента пиксела. Таким образом, выполнением п. 6 действительно успешно завершается выполнение алгоритма 1, представляющего изображения в виде набора “точечных” диполей, “масса” каждой половинки каждого из которых сосредоточена в одной точке, а соотношение плеч которых постоянно и равно N в направлении диполя, соответствующего всему изображению. В лемме 1 было доказано, что смещение ЦТ, вызванное тестовым ПП, направленным “навстречу” преобразованию, ранее исказившему анализируемое изображение, меньше по модулю, чем смещение, вызванное параллельно направленным тестовым ПП, если изображение представлено диполем, “масса” каждой из половинок которого сосредоточена в одной точке. Под направлением ПП здесь и далее подразумевается направление соответствующего ему вектора параметров P. Смещение диполя, перпендикулярное направлению ПП, согласно формулам (11) и (12) на полученные результаты влияния не оказывает. Важно только, чтобы ЦТ диполя не был смещен относительно начала координат в направлении, параллельном или антипараллельном направлению ПП. Представим произвольное двумерное изображение, ЦТ которого расположен в начале координат, в виде совокупности уравновешенных несимметричных “точечных” диполей, имеющих одинаковое соотношение длины плеч в направлении вектора Рt параметров тестового
“Оптический журнал”, 76, 7, 2009
57
ПП, и не имеющих смещения ЦТ относительно начала координат в параллельном или антипараллельном Pt направлениях, как это было показано в алгоритме 1. Тогда для любого изображения с ограниченной областью определения доказательство, приведенное в лемме 1, должно быть справедливо для каждого составляющего его “точечного” диполя и, следовательно, для всего изображения в целом, если ЦТ изображения имеет нулевую абсциссу. В общем случае для произвольного двумерного изображения направление соответствующего ему уравновешенного несимметричного диполя (направление максимальной асимметрии диполя) может не совпадать по направлению с осью абсцисс, однако, как следует из леммы 1, чем больше асимметрия диполя в направлении тестового ПП, тем больше будет разница скоростей смещения его ЦТ в направлениях параллельном и антипарллельном тестовому ПП. Это справедливо для любого точечного диполя, составляющего изображение, и для всего изображения в целом. Таким образом, вращая вектор Pt параметров тестового ПП, можно найти направление Pt, соответствующее максимальной разности скоростей смещения ЦТ, и это направление будет совпадать с направлением диполя, имеющего максимальный коэффициент N асимметрии.
В начале этого раздела на основании формул (11) и (12) было продемонстрировано, что под действием ПП происходит перераспределение точек изображения относительно его ЦТ. Если исходное изображение соответствовало симметричному уравновешенному диполю, т. е. для него выполнялось условие (14), то в результате действия некоторого ПП условие (14) нарушается, и преобразованное изображение будет соответствовать несимметричному диполю, направление которого (направление, соответствующее максимальному значению коэффициента N асимметрии) в общем случае не совпадает по направлению с вектором параметров P имевшего место ПП. Однако из формул (11) и (12) следует (достаточно для наглядности направить вектор P вдоль оси абсцисс), что направление смещения точек под действием ПП, направление максимального разрежения изображения под действием ПП и направление возникающего в результате действия ПП несимметричного диполя (направление большего плеча диполя) всегда имеет неположительную проекцию на вектор P параметров этого ПП. Следовательно, небольшое по величине тестовое ПП, вектор Pt параметров которого совпадает по направлению с более длинным плечом диполя,
будет иметь неположительную проекцию на вектор P параметров исходного ПП. Такое тестовое ПП (назовем его компенсирующим ПП) будет частично компенсировать исходное ПП. Назовем это утверждение леммой 2 и докажем его.
Доказательство леммы 2. Рассмотрим сначала, как строится вектор параметров итогового ПП, являющегося результатом нескольких последовательно выполненных ПП. Пусть P1, P2 и P3 – векторы параметров ПП, | A| и | B| – матрицы АП. Тогда смысл формул (11) и (12) можно выразить единственной формулой (23), а при отсутствии “аффинной составляющей” преобразование выразится формулой (24). Будем называть преобразование вида (24) непосредственно проективным преобразованием (НПП).
X′= ( A ⋅ X)/(X⋅ P +1),
(23)
X′ = X / (X⋅ P +1).
(24)
Тогда последовательное выполнение первого
НПП с вектором параметров P1, АП, описываемого матрицей | A|, и второго НПП с вектором па-
раметров P2 можно записать выражением
X′= ⎡⎢⎣( A ⋅X)/(X⋅P1+1)⎤⎥⎦/⎣⎢⎡( A ⋅X)/(X⋅P1+1)⋅P2 +1⎥⎦⎤ = = (B ⋅X)/(X⋅P3 +1),
из которого после выполнения несложных преобразований следует, что матрица | B| аффинной составляющей результирующего ПП равна матрице | A| исходного, а вектор P3 параметров НПП, соответствующего результирующему преобразованию, может быть вычислен по формуле (25). Формула (25), если пренебречь матрицей | A| аффинного преобразования, показывает, что вектор параметров результирующего НПП равен сумме векторов параметров последовательно примененных НПП.
P3 = A ⋅ P2 + P1.
(25)
Рисунок показывает результат действия компенсирующего ПП, вектор Pt параметров которого имеет неположительную проекцию на вектор P исходного ПП. Видно, что вектор параметров результирующего ПП становится короче. Как следует из формул (11) и (12), теоретически возможно, что направление диполя почти ортогонально вектору P. Это может произойти, если размер изображения в направлении вектора P был на много меньше, чем в ортогональном вектору P направлении. Однако вектор Pt никогда не ортогонален P, поэтому и в этом случае с помощью выбора достаточно короткого (но не бес-
58 “Оптический журнал”, 76, 7, 2009
gf
d
a
b c
z
Различные варианты суммирования векторов параметров основного P и компенсирующего Pt ПП: a – вектор параметров основного ПП, z – окружность радиуса ||a|| , описываемая при вращении этого вектора, b, d, f – различные положения вектора компенсирующего ПП (бесконечно короткий вектор b ортогонален вектору a и лежит на окружности z, длинный вектор d ортогонален вектору a, но не лежит на окружности z), c, e, g – векторы результирующего ПП, являющиеся соответственно результатами сложения вектора a с векторами b, d, f компенсирующего ПП.
конечно короткого) вектора Pt можно добиться чтобы ||P + Pt|| < ||P||. В то же время, если Pt почти ортогонален P, скорость уменьшения ||P + Pt|| и сходимость всего итеративного алгоритма компенсации ПП будет очень низкой. Во избежание таких ситуаций на каждой итерации частичной компенсации ПП перед оцениванием направления компенсирующего вектора Pt необходимо устранить существенное удлинение изображения объекта в каком-либо направлении. Это не входит в противоречие с предлагаемой процедурой компенсации ПП, поскольку является обязательным этапом компенсации АП, описанной в разделе 1. Кроме того, оценивание направления диполя и соответствующего ему направления вектора Pt производится после выполнения такого компенсирующего АП, поэтому выполненным на каждой итерации АП можно пренебречь, как это было предложено выше при выводе формулы (25). Таким образом, применение компенсирующего ПП с вектором параметров Pt, выбранным описанным выше способом, обеспечивает укорачивание вектора параметров результирующего ПП по сравнению с вектором параметров P исходного ПП (частично компенсирует исходное ПП), поэтому лемма 2 доказана.
Частичная компенсация исходного ПП на каждой итерации обеспечивает сходимость предлагаемой итеративной процедуры компенсации ПП изображения. Соответствующий этой процедуре алгоритм можно описать следующим образом.
Алгоритм 2: 1. Компенсировать сдвиг ЦТ изображения объекта относительно начала координат и компенсировать преобразование |A3| неравномерного удлинения изображения объекта, как это было описано в разделе 1. 2. Оценить направление вектора Pt параметров небольшого компенсирующего ПП (это такое направление, при котором модуль ΔC1 смещения ЦТ под действием ПП с вектором параметров Pt максимально отличается от модуля ΔC2 смещения ЦТ под действием ПП с противоположно направленным вектором Pt). Оценить величины ΔC1 и ΔC2, соответствующие выбранному вектору Pt. 3. Закончить выполнение алгоритма, если |ΔC1 – ΔC2| много меньше ΔC1 и |ΔC1 – ΔC2| много меньше ΔC2. В противном случае выполнить ПП, соответствующее выбранному Рt и перейти к п. 1 (здесь |.| обозначает вычисление абсолютной величины). Параметры P имевшего место ПП изображения могут быть вычислены по окончании выполнения описанного итеративного алгоритма как результат последовательного сложения согласно формуле (25) векторов Pt компенсирующих ПП.
Заключение
Предложенный алгоритм (алгоритм 2) реализован в виде программной модели и успешно испытан при компенсации ПП разнообразных изображений. Описание модели и проведенных экспериментов можно найти в более ранней работе автора [5]. В то время не удалось обосновать сходимость предлагаемого алгоритма при нормализации произвольных изображений, в связи с чем и написана настоящая статья. В заключение интересно отметить, что схожие “алгоритмы” используются в человеческом зрении, в котором изменение наклона поверхности и связанное с ним ПП определяются по скорости смещения текселов на изображении [6].
ЛИТЕРАТУРА
1. Qin-Sheng Chen, Defrise M., Deconinck F. Symmetric phase-only matched filtering of Fourier-
“Оптический журнал”, 76, 7, 2009
59
Mellin transforms for image registration and recognition // IEEE Transactions on Pattern Analysis and Machine Intelligence. 1994. V. 16. № 12. Р. 1156–1167.
2. Maitra S. Moment invariants // Proc. IEEE. 1979. V. 67. № 4. P. 697–699.
3. Луцив В.Р. Визуальные системы роботов. Глава VIII / Программное управление оборудованием. Мясников В.А., Игнатьев М.Б., Покровский А.М. 2-е изд. Л.: Машиностроение, 1984. C. 354–381.
4. Луцив В.Р. Методы и средства управления промышленным оборудованием на основе визуаль-
ной информации // Диссертация на соискание ученой степени кандидата технических наук. Л.: ЛИАП, 1985. 183 с.
5. Луцив В.Р. Разработка методов оптимизационного системного проектирования механических конструкций, приводов, систем управления и адаптации роботизированных комплексов // Отчет по Г.Б. НИР–Г-3-4-14-86 “Роботы”. Промежуточный. Инв. № 089857, № Г.Р. 01860016170 / Л.: ЛИАП, 1987. С. 34–50.
6. Марр Д. Зрение. Информационный подход к изучению представления и обработки зрительных образов. Пер. с англ. / М.: Радио и связь, 1987. 400 с.
60 “Оптический журнал”, 76, 7, 2009
УДК 621.397.3: 536.3
МЕТОД ИТЕРАТИВНОЙ КОМПЕНСАЦИИ ПРОЕКТИВНЫХ ИСКАЖЕНИЙ ИЗОБРАЖЕНИЙ
© 2009 г. В. Р. Луцив, канд. техн. наук
НПК “Государственный оптический институт им. С.И. Вавилова”, Санкт-Петербург
E-mail: vluciv@list.ru
Представлен метод нормализации изображения относительно полной группы проективных преобразований на плоскости, основанный на итеративной оптимизации решения в пространстве двух параметров проективного преобразования, отличающих его от аффинного. Эта пара параметров представляется вектором. Каждой итерации нормализации предшествует измерение параметров и компенсация аффинного преобразования согласно формулам, выведенным и опубликованным автором ранее. Для частичной компенсации проективного преобразования на каждой итерации используется такой по направлению вектор нормализующего преобразования, что изменению его направления на противоположное соответствует наибольшее по модулю изменение скорости смещения центра тяжести преобразуемого изображения. Итеративный процесс сходится к состоянию изображения, эталонному относительно проективного преобразования, при котором вектор нормализующего преобразования любого направления приводит к одинаковому по модулю смещению центра тяжести.
Ключевые слова: афинные преобразования, проективные преобразования, смещение центра тяжести, вектор нормализующего преобразования
Коды OCIS: 100.2980, 3000.3020
Поступила в редакцию 27.02.2009
Введение
Измерение параметров и компенсация аффинных преобразований (АП) и проективных преобразований (ПП) и поиск инвариантов этих преобразований являются важными направлениями автоматического анализа изображений. Это связано с тем, что в большинстве систем компьютерного зрения изображение на матрице фоточувствительных элементов формируется путем проецирования через сферическую линзу. Вследствие этого наблюдаемые изображения двумерных проекций объектов претерпевают на плоской матрице фотоприемников, а затем и в памяти компьютера геометрические трансформации, описываемые проективной группой преобразований на плоскости. Подгруппой ПП является группа АП, соответствующая случаю, когда объектив системы компьютерного зрения находится достаточно далеко от наблюдаемой сцены. Известны методы [1, 2] вычисления инвариантов группы АП на плоскости. Известен также альтернативный подход [3], в соответствии с которым сначала могут быть в аналитической фор-
ме определены и измерены параметры АП, а затем это преобразование может быть компенсировано, причем при измерении параметров преобразования информация о типе анализируемого объекта (об эталоне) не используется. ПП в отличие от АП нелинейно. Это создает препятствия как для вычисления инвариантов ПП, так и для построения аналитических выражений, пригодных для определения параметров ПП по содержанию изображения. В статье предлагается метод итеративной нормализации изображений относительно их ПП. Применяемый подход является расширением подхода к компенсации АП, предложенного автором в работах [3, 4]: в настоящей работе в результате нормализации изображение приводится к состоянию, эталонному относительно ПП, и измеряются параметры этого преобразования.
1. Измерение параметров и компенсация АП
В работах [3, 4] показано, что АП |A| может быть разложено на простые геометрические трансформации следующим образом:
“Оптический журнал”, 76, 7, 2009
53
a1 a2 a3 A = a4 a5 a6 =
001
= A1 ⋅ A2 ⋅ A3 ⋅ A4 ⋅ A5 ⋅ A6 ,
(1)
где a1–a6 – параметры АП, представленного в матричной форме, а матрицы | Ai| описывают в том порядке, в котором они пронумерованы в
выражении (1), смещение центра тяжести (ЦТ)
по двум декартовым координатным осям, вытя-
гивание или сжатие в каком-либо направлении,
зеркальное отражение относительно прямой,
проходящей через начало координат, масшта-
бирование и вращение относительно начала
координат.
Представленный в выражении (1) порядок
простых преобразований может быть изменен,
однако, именно такой порядок позволяет произ-
водить поэтапное измерение их параметров и их
компенсацию, начиная с | A1| и заканчивая | A6|, поскольку именно при компенсации простых
преобразований | Ai| в порядке от | A1| к | A6| последующие операции компенсации не нарушают
результатов предыдущих. При этом преобразова-
ния | A1| и | A2| взаимно ортогональны и могут быть переставлены местами в выражении (1). Анало-
гично могут быть переставлены местами преоб-
разования | A5| и | A6|. Преобразование | A4| может быть переставлено местами с преобразованием
| A5|, и оно может измеряться и компенсироваться либо одновременно с преобразованием | A3|, либо одновременно с преобразованием | A6|.
Итак, в качестве первого шага нормализации
декартовы координаты xc и yc ЦТ объекта вычисляются по следующим формулам:
xc = ∫∫ xB(x,y)dxdy ∫∫ B(x,y)dxdy,
ΩΩ
yc = ∫∫ yB(x,y)dxdy ∫∫ B(x,y)dxdy,
ΩΩ
(2)
где Ω – двумерная область определения изображения объекта, B(х, у) – яркость пиксела с координатами {х, у}. Это дает возможность компенсировать линейные смещения | A1| и | A2| изображения объекта.
Следующим преобразованием, подлежащим
измерению и компенсации, является масштабирование | A3| вдоль некоторого направления. В работах [3, 4] получены следующие аналити-
ческие выражения для вычисления параметров преобразования | A3| – направления Θ и величины μ изменения масштаба изображения объекта:
( ) ( )μ = ± D + C2+ E2 D − C2+ E2 , (3)
Θ =0,5 arctg(C/E)+πn,
(4)
E = ∫∫ B(r,ϕ)r2cos(2ϕ)drdϕ,
Ω
C = ∫∫ B(r,ϕ)r2sin(2ϕ)drdϕ,
Ω
(5) (6)
D = ∫∫ B(r,ϕ)r2drdϕ,
Ω
(7)
где B(r,ϕ) – яркость пиксела центрированного изображения объекта, представленного в поляр-
ной системе координат, Ω – область определения изображения объекта в полярной системе коор-
динат, {r, ϕ} – соответственно радиус и азимут пиксела в полярной системе координат.
После компенсации преобразования | A3| остаются некомпенсированными преобразования | A4|, | A5| и | A6| зеркального отражения, равномерного масштабирования и вращения. Преобразование
| A5| не зависит от углового положения пиксела и не влияет на | A6| и | A4|, поэтому параметр М
преобразования
A5
=
M 0
0 M можно измерить
по проекции изображения на радиус полярной системы координат. Коэффициент М может быть измерен согласно следующему выражению по смещению ЦТ проекции изображения на радиус относительно эталонного положения А:
M = ∫∫ rB(r,ϕ)drdϕ ∫∫ B(r,ϕ)drdϕ/A. (8)
ΩΩ
Вращение центрированного изображения относительно начала координат и его зеркальное отражение относительно прямой, проходящей через начало координат, не изменяют радиальные координаты пикселов, поэтому параметры этих преобразований можно измерить по проекции изображения на окружность, т. е. на угловую полярную координату. Параметры Ψ вращения и ε зеркального отражения могут [4] быть вычислены по формулам (9) и (10)
Ψ = ϕk − ϕm − ε(ψk − ψm ),
(9)
ε = (2ϕk − ϕn − ϕm )/(2ψk − ψn − ψm ), (10)
где ϕm, ϕk, ϕn – фазы трех соседних Фурьегармоник (m = k – 1, n = k + 1) проекции анализируемого изображения на окружность, ψm, ψk,
54 “Оптический журнал”, 76, 7, 2009
ψn – фазы аналогичных гармоник для соответствующего эталонного изображения. Возможны и другие пути вычисления параметров Ψ и ε [4], которые могут быть использованы, если в спектре проекции изображения на окружность подходящих соседних гармоник не обнаружено.
2. Итеративная компенсация и измерение параметров ПП
ПП на плоскости описывается следующими выражениями, связывающими координаты исходного X = |x, y|Tи преобразованного X′ = |x′, y′|T изображений
x′ = (a1x + a2y + a3 )/(a7x + a8y +1)= = (a1x + a2y + a3 )/(X⋅ P +1),
(11)
y′ = (a4x + a5y + a6 )/(a7x + a8y +1)= = (a4x + a5y + a6 )/(X⋅ P +1),
(12)
где a1–a6 – параметры АП, упомянутые в формуле (1), а параметры a7 и a8 характеризуют отличие ПП от АП. Из выражений (11) и (12)
видно, что знаменатель формулы ПП можно выразить через скалярное произведение X·P
вектора декартовых координат точки и вектора параметров P = |a7, a8|T. ПП – нелинейно, поэтому не удается определить вектор P его параметров в
аналитическом виде, как это было сделано для
параметров АП. Проблема решается путем ите-
ративной оптимизации решения в пространстве
параметров a7 и a8. Разработанный метод основан на том, что ПП в отличие от АП приводит к
перераспределению точек изображения объекта относительно его центра тяжести, что несложно
показать аналитически. Представим формулу (11) в следующем виде:
x′ = α/(1+ sign(β)⋅ | β |),
(13)
где α соответствует числителю выражения (11), | . | обозначает модуль, β = X·P. Пусть для определенности вектор P антипараллелен оси абсцисс, т. е. его скалярное произведение с ортом оси абсцисс отрицательно. Пусть ЦТ изображения находится в начале координат. При β ≠ 0 x′ = = α/(1 – |β|) справа от оси ординат, а слева от оси ординат x′ = α/(1 + |β|). Путем несложных преобразований легко показать, что α/(1 + |β|) – α/1 < 0, α/(1 – |β|) – α/1 > 0 и |α/(1 + |β|) – α/1| < < |α/(1 – |β|) – α/1|. Из этого следует, что абсциссы точек, находящихся левее оси ординат, в результате приложенного ПП уменьшаются по модулю,
а абсциссы точек, находящихся правее оси ординат, увеличиваются по модулю, что приводит к смещению ЦТ в область положительных абсцисс. При этом правее оси ординат приращение абсцисс происходит быстрее, а по мере удаления от начала координат увеличивается степень растяжения (разрежения) изображения вдоль оси абсцисс в области положительных абсцисс, и увеличивается степень сжатия вдоль оси абсцисс в области отрицательных абсцисс. В связи с этим после компенсации смещения ЦТ слева от оси ординат будет сосредоточено больше точек изображения (эта часть изображения будет более массивна), но в среднем они будут иметь меньшие по модулю абсциссы, а справа от оси ординат изображение будет более разреженным: оно будет содержать меньшее число точек, имеющих в то же время большие абсциссы. Указанные построения справедливы и для произвольного направления вектора P, поскольку такие случаи можно свести к рассмотренному заменой системы координат.
Рассмотрим произвольное изображение, ЦТ которого имеет нулевые координаты. Разделим изображение на два полу-изображения произвольной прямой, проходящей через начало координат. Полученные полу-изображения можно представить в виде половинок диполя, массы M1 и M2 которых сконцентрированы в ЦТ полуизображений, отстоящих от начала координат соответственно на расстояния L1 и L2. Диполь уравновешен относительно ЦТ полного изображения, но в общем случае не симметричен, т. е. M1L1 = M2L2, но M1 ≠ M2 и L1 ≠L2, причем в каком-то направлении Θ отношение L2 к L1 максимально. Тогда изображение можно представить диполем, ориентированным в направлении Θ, и описать его вектором L, длина которого равна отношению L2 к L1, а направление соответствует направлению плеча L2 диполя. Эталонное состояние изображения такое, при котором
M1 = M2 , L1 = L2,
(14)
т. е. при котором оно представляет собой диполь, симметричный и уравновешенный относительно ЦТ изображения при любом направлении прямой, делящей его на два полу-изображения. В эталонном состоянии ||L|| = 1при любом направлении Θ плеча диполя. Примером эталонного изображения может служить круг с центром в начале координат. В соответствии с вышесказанным любое ПП с параметрами a7 ≠ 0 и/или a8 ≠0, приложенное к изображению круга, вызывает перераспределение его точек относительно ЦТ, что приводит к нарушению условия (14) и к сме-
“Оптический журнал”, 76, 7, 2009
55
щению ЦТ. В случае круга при вращении вектора Р возникает изменение направления смещения ЦТ, но модуль его смещения остается неизменным. Если изображение круга подверглось ранее ПП с параметрами P с последующим перемещением его ЦТ в начало координат, его изображение будет соответствовать несимметричному диполю, для которого ||L|| > 1. Тогда вращение вектора Pt параметров дополнительно приложенного к нему “тестового” ПП приводит к различным по модулю смещениям ЦТ, причем скорость смещения максимальна при некотором направлении Θ = Θmax вектора Pt, которое для изображения круга соответствует направлению вектора P ранее приложенного к нему ПП. Чтобы понять это, заметим, что круг симметричен, поэтому достаточно рассматривать его изображение, к которому было применено ПП, как диполь. Пусть вектор диполя L совпадает по направлению с осью абсцисс, а вектор P ранее приложенного к изображению ПП имеет противоположное направление. Докажем для диполя следующую лемму:
Лемма 1: модуль смещения ЦТ несимметричного диполя под действием ПП, вектор Pt параметров которого противоположен направлению диполя больше, чем модуль смещения ЦТ под действием ПП, для которого вектор Pt параллелен диполю.
Доказательство: Пусть М1 – масса (яркость) левой части диполя, L1 – модуль абсциссы его левой части, М2 – масса его правой части, L2 – модуль ее абсциссы. M1 = NM2, L2 = NL1, N > 1. Тогда в результате дополнительного применения к такому диполю тестового ПП с вектором параметров Pt (пусть Pt = || Pt || ), направленным также как и вектор P противоположно направлению оси абсцисс, преобразованные длины L1′ и L2′ плеч диполя в соответствии с выражениями (11) и (12) будут иметь вид (15)
L1′ = L1/(1+ L1Pt ), L2′ = L2/(1− L2Pt ). (15)
Приращения ΔL1 и ΔL2 координат половинок диполя вычисляются тогда как
ΔL1 = L1 − L1/(1+ L1Pt ) = L21Pt /(1+ L1Pt ), (16)
ΔL2 = L2/(1− L2Pt ) − L2 = L22Pt /(1− L2Pt ). (17)
На основе равенств (16) и (17) модуль ΔC1 приращения абсциссы ЦТ диполя можно описать формулой (18), при выводе которой использовались подстановки M1=NM2, L2=NL1.
ΔC1
=
ΔL1M1 + ΔL2M2 M1 + M2
=
=
(1
+
NM2L21Pt (1+ N) L1Pt )(1− NL1Pt )(M1
+
M2
)
.
(18)
Аналогично выводится формула (19), пока-
зывающая модуль приращения ΔC2 абсциссы ЦТ, если направление вектора Pt изменится на противоположное.
ΔC2
=
(1−
NM2L21Pt (1+ N) L1Pt )(1+ NL1Pt )(M1
+
M2 )
.
(19)
Сравним модули ΔC1 и ΔC2 смещений ЦТ, соответствующие противоположным направлениям
Pt (в выражении (20) знак ~ означает неизвестный знак неравенства)
(1
+
NM2L21Pt (1+ N) L1Pt )(1− NL1Pt )(M1
+
M2
)
~
~
(1−
NM2L21Pt (1+ N) L1Pt )(1+ NL1Pt )(M1
+
M2
)
.
(20)
1 – NL1Pt > 0 и 1 – L1Pt > 0, что необходимо, чтобы точки диполя не перемещались в бесконечность.
Поэтому, не нарушая неравенства, можно свести
(20) к виду (21)
(1 – L1Pt)(1 + NL1Pt) ~ ~ (1 + L1Pt)(1 – NL1Pt),
(21)
откуда после выполнения умножения и приведения подобных членов получаем N ~ 1. Выше предполагалось, что N > 1, поэтому, когда вектор параметров Pt прикладываемого тестового ПП совпадает по направлению с вектором параметров ПП, имевшего место ранее, скорость смещения ЦТ больше, чем при противоположно направленном векторе параметров тестового ПП (тогда вектор Pt параллелен диполю, а тестовое ПП частично компенсирует исходное). Итак, лемма 1 доказана.
Из доказательства леммы 1 видно, что чем больше асимметрия N диполя, тем больше разница скоростей смещения ЦТ при противоположных направлениях Pt.
Полученный результат корректен для изображения круга, однако, реальные изображения характеризуются более сложными законами пространственного распределения яркости пикселов, в общем случае не симметричными. Чтобы показать справедливость утверждения,
56 “Оптический журнал”, 76, 7, 2009
доказанного в лемме 1, для реальных образов
допустим, что произвольное изображение имело
эталонный вид, соответствующий условию (14),
и что в результате некоторого проективного пре-
образования и последующего переноса его ЦТ в
начало координат это изображение в целом при-
обрело вид несимметричного уравновешенного
диполя, ориентированного вдоль оси абсцисс, у
которого правое плечо L2 в N раз длиннее левого (M1 = NM2, L2 = NL1, N > 1). Покажем, что такое произвольное изображение может быть пред-
ставлено в виде совокупности уравновешенных
несимметричных диполей, соответствующих
отдельным парам пикселов или их фрагментов,
для каждой из которых M1 = NM2, L2 = NL1 при неизменном N > 1. Для этого выполним приве-
денный ниже алгоритм.
Алгоритм 1:
1. Выбираем любой оставшийся пиксел пра-
вее оси ординат, имеющий абсциссу l и яркость
(“массу”) ml. Если такой пиксел не найден, переходим к п. 6.
2. Выбираем любой пиксел левее оси орди-
нат, имеющий абсциссу L и яркость ML, для которого
LML ≤ lml.
(22)
Если такого пиксела не найдено, переходим к
п. 4, в противном случае “уравновешиваем” от-
носительно оси ординат найденный левее оси ординат пиксел с параметрами {L, ML} “фрагментом” ранее выбранного пиксела, имевшего параметры {l, ml} и находящегося правее оси ординат. Исключаем из дальнейшего рассмотрения
полностью уравновешенный пиксел слева от оси
ординат и “уравновешивающий” его фрагмент яркостью LML/l пиксела справа.
3. Оставшаяся неуравновешенной часть
пиксела справа от оси абсцисс имеет яркость ml′ = ml – LML/l. Если новая яркость этого пиксела ml′ > 0, переходим к п. 2, в противном случае переходим к п.1.
4. Выбираем любой пиксел слева с параметрами {L, ML}, для которого LML ≥ lml. Такой пиксел найдется, поскольку неравенство (22)
в п. 2 выполнено не было. Уравновешиваем пиксел справа, имеющий параметры {l, ml}, фрагментом найденного слева пиксела. Исклю-
чаем из дальнейшего рассмотрения полностью
уравновешенный пиксел справа от оси ординат
и “уравновешивающий” его фрагмент яркостью lml/L пиксела, найденного слева.
5. Оставшаяся неуравновешенной часть пиксела слева имеет яркость M′ = M – lml/L. Выби-
раем любой пиксел справа с параметрами {l, ml}, для которого LML′ ≥ lml. Если такой пиксел не найден, переходим к п. 1, в противном случае переходим к п. 4.
6. Конец выполнения алгоритма. Каждый раз при выполнении п. 2 и п. 4 из дальнейшего рассмотрения исключаются уравновешенные относительно оси ординат диполи, соответствующие парам пикселов или их частей. Исключаемые диполи имеют соотношение плеч N в направлении исходного диполя, который до начала выполнения настоящего алгоритма соответствовал всему изображению. Не трудно убедиться, что часть исходного диполя, не исключенная еще из рассмотрения после каждого выполнения п. 2 или п. 4, также уравновешена относительно оси ординат и имеет аналогичное соотношение плеч N в направлении исходного диполя. Поэтому, если при выполнении п. 1 оказалось, что правее оси ординат не найдено ни одного пиксела или фрагмента пиксела, еще не исключенного из рассмотрения, можно утверждать, что и слева от оси ординат не осталось ни одного подлежащего рассмотрению пиксела или фрагмента пиксела. Таким образом, выполнением п. 6 действительно успешно завершается выполнение алгоритма 1, представляющего изображения в виде набора “точечных” диполей, “масса” каждой половинки каждого из которых сосредоточена в одной точке, а соотношение плеч которых постоянно и равно N в направлении диполя, соответствующего всему изображению. В лемме 1 было доказано, что смещение ЦТ, вызванное тестовым ПП, направленным “навстречу” преобразованию, ранее исказившему анализируемое изображение, меньше по модулю, чем смещение, вызванное параллельно направленным тестовым ПП, если изображение представлено диполем, “масса” каждой из половинок которого сосредоточена в одной точке. Под направлением ПП здесь и далее подразумевается направление соответствующего ему вектора параметров P. Смещение диполя, перпендикулярное направлению ПП, согласно формулам (11) и (12) на полученные результаты влияния не оказывает. Важно только, чтобы ЦТ диполя не был смещен относительно начала координат в направлении, параллельном или антипараллельном направлению ПП. Представим произвольное двумерное изображение, ЦТ которого расположен в начале координат, в виде совокупности уравновешенных несимметричных “точечных” диполей, имеющих одинаковое соотношение длины плеч в направлении вектора Рt параметров тестового
“Оптический журнал”, 76, 7, 2009
57
ПП, и не имеющих смещения ЦТ относительно начала координат в параллельном или антипараллельном Pt направлениях, как это было показано в алгоритме 1. Тогда для любого изображения с ограниченной областью определения доказательство, приведенное в лемме 1, должно быть справедливо для каждого составляющего его “точечного” диполя и, следовательно, для всего изображения в целом, если ЦТ изображения имеет нулевую абсциссу. В общем случае для произвольного двумерного изображения направление соответствующего ему уравновешенного несимметричного диполя (направление максимальной асимметрии диполя) может не совпадать по направлению с осью абсцисс, однако, как следует из леммы 1, чем больше асимметрия диполя в направлении тестового ПП, тем больше будет разница скоростей смещения его ЦТ в направлениях параллельном и антипарллельном тестовому ПП. Это справедливо для любого точечного диполя, составляющего изображение, и для всего изображения в целом. Таким образом, вращая вектор Pt параметров тестового ПП, можно найти направление Pt, соответствующее максимальной разности скоростей смещения ЦТ, и это направление будет совпадать с направлением диполя, имеющего максимальный коэффициент N асимметрии.
В начале этого раздела на основании формул (11) и (12) было продемонстрировано, что под действием ПП происходит перераспределение точек изображения относительно его ЦТ. Если исходное изображение соответствовало симметричному уравновешенному диполю, т. е. для него выполнялось условие (14), то в результате действия некоторого ПП условие (14) нарушается, и преобразованное изображение будет соответствовать несимметричному диполю, направление которого (направление, соответствующее максимальному значению коэффициента N асимметрии) в общем случае не совпадает по направлению с вектором параметров P имевшего место ПП. Однако из формул (11) и (12) следует (достаточно для наглядности направить вектор P вдоль оси абсцисс), что направление смещения точек под действием ПП, направление максимального разрежения изображения под действием ПП и направление возникающего в результате действия ПП несимметричного диполя (направление большего плеча диполя) всегда имеет неположительную проекцию на вектор P параметров этого ПП. Следовательно, небольшое по величине тестовое ПП, вектор Pt параметров которого совпадает по направлению с более длинным плечом диполя,
будет иметь неположительную проекцию на вектор P параметров исходного ПП. Такое тестовое ПП (назовем его компенсирующим ПП) будет частично компенсировать исходное ПП. Назовем это утверждение леммой 2 и докажем его.
Доказательство леммы 2. Рассмотрим сначала, как строится вектор параметров итогового ПП, являющегося результатом нескольких последовательно выполненных ПП. Пусть P1, P2 и P3 – векторы параметров ПП, | A| и | B| – матрицы АП. Тогда смысл формул (11) и (12) можно выразить единственной формулой (23), а при отсутствии “аффинной составляющей” преобразование выразится формулой (24). Будем называть преобразование вида (24) непосредственно проективным преобразованием (НПП).
X′= ( A ⋅ X)/(X⋅ P +1),
(23)
X′ = X / (X⋅ P +1).
(24)
Тогда последовательное выполнение первого
НПП с вектором параметров P1, АП, описываемого матрицей | A|, и второго НПП с вектором па-
раметров P2 можно записать выражением
X′= ⎡⎢⎣( A ⋅X)/(X⋅P1+1)⎤⎥⎦/⎣⎢⎡( A ⋅X)/(X⋅P1+1)⋅P2 +1⎥⎦⎤ = = (B ⋅X)/(X⋅P3 +1),
из которого после выполнения несложных преобразований следует, что матрица | B| аффинной составляющей результирующего ПП равна матрице | A| исходного, а вектор P3 параметров НПП, соответствующего результирующему преобразованию, может быть вычислен по формуле (25). Формула (25), если пренебречь матрицей | A| аффинного преобразования, показывает, что вектор параметров результирующего НПП равен сумме векторов параметров последовательно примененных НПП.
P3 = A ⋅ P2 + P1.
(25)
Рисунок показывает результат действия компенсирующего ПП, вектор Pt параметров которого имеет неположительную проекцию на вектор P исходного ПП. Видно, что вектор параметров результирующего ПП становится короче. Как следует из формул (11) и (12), теоретически возможно, что направление диполя почти ортогонально вектору P. Это может произойти, если размер изображения в направлении вектора P был на много меньше, чем в ортогональном вектору P направлении. Однако вектор Pt никогда не ортогонален P, поэтому и в этом случае с помощью выбора достаточно короткого (но не бес-
58 “Оптический журнал”, 76, 7, 2009
gf
d
a
b c
z
Различные варианты суммирования векторов параметров основного P и компенсирующего Pt ПП: a – вектор параметров основного ПП, z – окружность радиуса ||a|| , описываемая при вращении этого вектора, b, d, f – различные положения вектора компенсирующего ПП (бесконечно короткий вектор b ортогонален вектору a и лежит на окружности z, длинный вектор d ортогонален вектору a, но не лежит на окружности z), c, e, g – векторы результирующего ПП, являющиеся соответственно результатами сложения вектора a с векторами b, d, f компенсирующего ПП.
конечно короткого) вектора Pt можно добиться чтобы ||P + Pt|| < ||P||. В то же время, если Pt почти ортогонален P, скорость уменьшения ||P + Pt|| и сходимость всего итеративного алгоритма компенсации ПП будет очень низкой. Во избежание таких ситуаций на каждой итерации частичной компенсации ПП перед оцениванием направления компенсирующего вектора Pt необходимо устранить существенное удлинение изображения объекта в каком-либо направлении. Это не входит в противоречие с предлагаемой процедурой компенсации ПП, поскольку является обязательным этапом компенсации АП, описанной в разделе 1. Кроме того, оценивание направления диполя и соответствующего ему направления вектора Pt производится после выполнения такого компенсирующего АП, поэтому выполненным на каждой итерации АП можно пренебречь, как это было предложено выше при выводе формулы (25). Таким образом, применение компенсирующего ПП с вектором параметров Pt, выбранным описанным выше способом, обеспечивает укорачивание вектора параметров результирующего ПП по сравнению с вектором параметров P исходного ПП (частично компенсирует исходное ПП), поэтому лемма 2 доказана.
Частичная компенсация исходного ПП на каждой итерации обеспечивает сходимость предлагаемой итеративной процедуры компенсации ПП изображения. Соответствующий этой процедуре алгоритм можно описать следующим образом.
Алгоритм 2: 1. Компенсировать сдвиг ЦТ изображения объекта относительно начала координат и компенсировать преобразование |A3| неравномерного удлинения изображения объекта, как это было описано в разделе 1. 2. Оценить направление вектора Pt параметров небольшого компенсирующего ПП (это такое направление, при котором модуль ΔC1 смещения ЦТ под действием ПП с вектором параметров Pt максимально отличается от модуля ΔC2 смещения ЦТ под действием ПП с противоположно направленным вектором Pt). Оценить величины ΔC1 и ΔC2, соответствующие выбранному вектору Pt. 3. Закончить выполнение алгоритма, если |ΔC1 – ΔC2| много меньше ΔC1 и |ΔC1 – ΔC2| много меньше ΔC2. В противном случае выполнить ПП, соответствующее выбранному Рt и перейти к п. 1 (здесь |.| обозначает вычисление абсолютной величины). Параметры P имевшего место ПП изображения могут быть вычислены по окончании выполнения описанного итеративного алгоритма как результат последовательного сложения согласно формуле (25) векторов Pt компенсирующих ПП.
Заключение
Предложенный алгоритм (алгоритм 2) реализован в виде программной модели и успешно испытан при компенсации ПП разнообразных изображений. Описание модели и проведенных экспериментов можно найти в более ранней работе автора [5]. В то время не удалось обосновать сходимость предлагаемого алгоритма при нормализации произвольных изображений, в связи с чем и написана настоящая статья. В заключение интересно отметить, что схожие “алгоритмы” используются в человеческом зрении, в котором изменение наклона поверхности и связанное с ним ПП определяются по скорости смещения текселов на изображении [6].
ЛИТЕРАТУРА
1. Qin-Sheng Chen, Defrise M., Deconinck F. Symmetric phase-only matched filtering of Fourier-
“Оптический журнал”, 76, 7, 2009
59
Mellin transforms for image registration and recognition // IEEE Transactions on Pattern Analysis and Machine Intelligence. 1994. V. 16. № 12. Р. 1156–1167.
2. Maitra S. Moment invariants // Proc. IEEE. 1979. V. 67. № 4. P. 697–699.
3. Луцив В.Р. Визуальные системы роботов. Глава VIII / Программное управление оборудованием. Мясников В.А., Игнатьев М.Б., Покровский А.М. 2-е изд. Л.: Машиностроение, 1984. C. 354–381.
4. Луцив В.Р. Методы и средства управления промышленным оборудованием на основе визуаль-
ной информации // Диссертация на соискание ученой степени кандидата технических наук. Л.: ЛИАП, 1985. 183 с.
5. Луцив В.Р. Разработка методов оптимизационного системного проектирования механических конструкций, приводов, систем управления и адаптации роботизированных комплексов // Отчет по Г.Б. НИР–Г-3-4-14-86 “Роботы”. Промежуточный. Инв. № 089857, № Г.Р. 01860016170 / Л.: ЛИАП, 1987. С. 34–50.
6. Марр Д. Зрение. Информационный подход к изучению представления и обработки зрительных образов. Пер. с англ. / М.: Радио и связь, 1987. 400 с.
60 “Оптический журнал”, 76, 7, 2009