CC BY
38УДК 68i.5i8.3
МЕТОД ОПРЕДЕЛЕНИЯ ВЕКТОРА СКОРОСТИ ДВИЖЕНИЯ ПОДСТИЛАЮЩЕЙ ПОВЕРХНОСТИ
П.К. Кузнецов1, Б.В. Мартемьянов, В.И. Семавин, Е.Ю. Чекотило
Самарский государственный технический университет 443100, Самара, ул. Молодогвардейская, 244
Предлагается метод определения вектора скорости движения носителя изобразительной системы относительно подстилающей поверхности, основанный на анализе видеоинформации, содержащейся в двух последовательных кадрах изображения подстилающей поверхности. Метод не накладывает ограничений на число степеней свободы движения носителя изобразительной системы, не использует таких сложных вычислительных процедур, как поиск экстремума многомерных функций или вычисления сверток, и поэтому является перспективным для разработки «быстрых» способов и алгоритмов определения параметров движения для реализации в системах, работающих в реальном времени.
Ключевые слова: оптический поток, векторное поле скоростей, параметры движения, изображение, метод функционализации
Введение
В настоящее время активно развиваются методы обработки изображений, основанные на вычислении поля векторов скорости кажущегося движения изображения (оптического потока (optical flow) по терминологии, принятой в зарубежной литературе [1-3])), создаваемого движущимся объектом. Информация об оптическом потоке используется, например, в компьютерной графике для синтеза деформируемых и движущихся изображений. Последним обстоятельством в особенности вызван значительный интернет-интерес к методу оптического потока в профессиональной среде создателей анимационных изображений [10].
Известные способы построения оптического потока [1-3] основаны на использовании уравнения, которое устанавливает связь изменения освещенности точек изображения с параметрами вектора кажущейся скорости движения изображения. В зарубежных источниках [3] это уравнение называют brightness constancy constraint equation (BCCE). Оно имеет вид
dE (х, у) = дЕ(х1у) дЕ (х у) дЕ ^, у)
dt dt х дх у ду ’ 1 ^
где Е(х, у) - функция освещенности изображения (ФОИ);
(х, у) - координаты точки изображения на картинной плоскости;
i Кузнецов Павел Константинович, доктор технических наук, профессор.
e-mail: kurnesov@mail. ru Мартемьянов Борис Викторович, кандидат технических наук, доцент. Семавин Владимир Иванович, кандидат технических наук, доцент. Чекотило Елена Юрьевна, ассистент.
их и ьу - компоненты вектора кажущейся скорости движения изображения.
Методы определения оптического потока, основанные на использовании данного соотношения, относят классу градиентных методов.
Градиентные методы имеют существенные ограничения, связанные с необходимостью точного вычисления производных по времени и пространству от функций освещенности изображения, которые в общем случае не являются ни дифференцируемыми, ни даже непрерывными. Это приводит к появлению значительной методической погрешности в определении скорости даже при отсутствии случайных шумов в видеосигнале. В [4] показано, что градиентный метод дает оценки скорости, смещенные в сторону уменьшения модуля измеряемой скорости. Замена производных первыми разностями вводит дополнительные существенные погрешности в измерения.
Кроме того, при решении большого класса задач обработки изображений, для которых требуется знание параметров движения носителя изобразительной системы, одной лишь информации об оптическом потоке недостаточно. Примером может служить задача обеспечения требуемых измерительных свойств аэрокосмических снимков, задача обнаружения и сопровождения движущихся объектов. Дело в том, что вектор оптического потока, хотя и содержит информацию о параметрах движения наблюдаемого объекта, но в общем случае в скрытой, неявной форме, и для того, чтобы получить такую информацию, необходимо проводить дополнительные, достаточно трудоемкие вычисления.
Здесь развивается прямой метод определения параметров движения объекта, наблюдаемого изобразительной системой, который не требует предварительного вычисления оптического потока. Метод, названный в [5-9] методом функционализации оптических параметров изображения (для краткости далее будем называть его методом функционализации), является обобщением градиентного метода, но свободен от его недостатков, перечисленных выше. Метод основан на использовании уравнения функциональной связи типа (1), но такого вида, в котором коэффициентами являются величины, доступные физическому измерению (физически измеримые), именно -облученности ограниченных подобластей изображения ненулевой меры, а искомыми неизвестными - параметры полного вектора скорости относительного движения подстилающей поверхности.
1. Метод функционализации в задаче определения параметров движения
носителя оптической системы относительно подстилающей поверхности
Метод функционализации [5-9] отличается от градиентного метода тем, что в нем вместо освещенности отдельных точек изображения используют специального вида функционалы, определенные на подобластях (окнах) изображения. Это придает методу существенно новые свойства и возможности по точности, быстродействию и универсальности. Система ФС-уравнений при определенном выборе основного функционала, используемого в методе, получается линейной относительно компонент искомого вектора скорости движения носителя изобразительной системы (ИС). Ниже излагается метод функционализации в общем виде для случая, когда основной функционал, используемый в методе, имеет мультипликативный вид. Впервые идея метода изложена в [5], а в работе Шалкова (8сЬа1коу) [6] предлагается частный вариант метода, который вытекает из метода функционализации в случае, когда основной функционал имеет вид интегральной освещенности окна анализа.
1.1. Модель формирования изображения движущейся подстилающей поверхности
Будем считать, что носитель ИС совершает относительно подстилающей поверхности поступательное движение и движение вращения относительно своего центра масс и что подстилающая поверхность недеформируемая, плоская, обладает диффузной отражательной способностью и равномерно освещена сторонним источником излучения; расстояние между ИС и подстилающей поверхностью значительно больше заднего фокусного расстояния ИС, а ИС изопланатична.
На рис. 1 приведены системы координат, в которых рассматривается модель формирования изображения:
- 0ХУ2, неподвижна относительно плоскости Р подстилающей поверхности, а её оси 0Х и 0У принадлежат Р;
- система координат о'х'у'г', оси которой жестко связаны с носителем ИС, а центр о' совпадает с центром масс носителя ИС;
- охуг, оси которой жестко связаны с носителем ИС, параллельны одноименным осям системы координат о'х'у'г', центр о совпадает с передней главной точкой ИС, при этом главная ось ИС совпадает с осью ог;
- окхкук, оси которой параллельны одноименным осям системы координат оху, центр ок находится в картинной плоскости (Рк) ИС на главной оси ИС на расстоянии -/ от центра о системы координат охуг (/- заднее фокусное расстояние ИС).
Положение системы координат охуг относительно системы координат о'х'у'г' задано вектором г$=\х$,у$,г$\ (здесь и далее - знак транспонирования). Положение системы координат о'х'у'г' относительно системы координат 0ХУ2 задано вектором Кн=Кн(0=\Хя(0, Уя(0,2я(0\т и матрицей преобразования координат Л=Л(0=\аі>;(0\ (і - номер строки, у - номер столбца, і, у = 1,2,3).
Матрица Л - матрица направляющих косинусов. Ее элементы - косинусы углов между осями систем координат о'х'у'г' и 0ХУ2:
со8(о х', 0Х) со8(ох', 0У) со8(ох' , 02)
со8(о/у', 0Х) со8(о'у', 0У) со8(о'у', 02)
соБ(о'г', 0Х) со$(о'г', 0У) со8(ог', 02)
Р и с. 1. Системы координат
л = Л(0 =
(2)
— 1 т
Известно, что матрица Л ортогональна и обладает свойством Л = Л .
Каждой точке плоскости Р соответствует ее изображение в плоскости Рк (например, на рис. 1 точке М соответствует ее изображение в точке т). Положение точки М с Р в системе координат 0ХУ2 задано на рис. 1 вектором К=\Х,У,0\т, а в подвижной
системе координат охуі - вектором гР=гР(ї)[хР(ї),уР(ї),2Р(ї)] . Положение точки т с Рк в системе координат охуі на рис. 1 задано вектором г=г(ґ)=[х(ґ), у(ґ), -/\Т.
Из рис. 1 следует векторное равенство:
f = kmfP = km (A(R — RH )— fS )
(З)
где кт = кт(г) = \т\/\тр\ = -//(Азф-Яц) - 25) = - АзТ т/фи + А3Т т) |т| и |тМ - модули
векторов т и тМ;
А3Т=[а31, а32, а33] - матрица-строка.
С учетом свойства ортогональности матрицы А соотношение, задающее обратное преобразование координат точки, находящейся в картинной плоскости ИС, в
координаты соответствующей ей точки на подстилающей поверхности, имеет вид
/ \
R = AT
■ + f?
+ R
н
(4)
Освещенность Е(т,1) в точке, заданной вектором т, связана с яркостью в{я) в точке подстилающей поверхности, заданной вектором Я, известным соотношением
е( t )=N_cos4 p(R B(R, t)
(5)
где T - коэффициент пропускания ИС; 1/N - относительное отверстие ИС, равное отношению диаметра входного зрачка ИС к заднему фокусному расстоянию;
К*) - угол между оптической осью ИС и направлением на точку, заданную вектором R, на подстилающей поверхности. В дальнейшем будем считать b(*)»0, а cos4 b(i?)» 1.
1.2. Уравнение движения изображения (уравнение оптического потока)
Для получения уравнения движения изображения подстилающей поверхности ИС продифференцируем соотношение (2) по времени:
f = [u x, u y ,0] =
f rA3 aT + AAT p + km fs ) — km
—rA3 + A
f
V
н
(б)
где Г = [и х, и у ,0] - вектор скорости движения изображения в картинной плоскости
ИС;
Ун = Ян = \уш, ¥уН, У2Н ]т - вектор скорости поступательного движения носителя ИС относительно подстилающей поверхности.
Для матрицы А в силу ее ортогональности выполняется:
(7)
(8)
A A T = кососимметричная матрица вида —AA T = —fi,
' 0 — w3 «2
fi = «3 0 — w1
— «2 w1 0
где элементы м1;ю2,юз матрицы П являются проекциями вектора бО=а(0=[Ф1(0>
02(0, (0зО угловой скорости носителя ИС на оси системы координат о'х'у^' (см. рис. 1), связанной с носителем ИС .
С учетом (7-8) из (6) получим уравнение движения изображения в картинной плоскости ИС в форме, удобной для дальнейшего его применения:
Г = -
1
— гЛ3 + Л
у
н ■
(9)
Уравнение (9) позволяет определять вектор кажущейся скорости движения изображения в любой точке картинной плоскости, если известен полный вектор скорости движения носителя ИС, и тем самым вычислять оптический поток, создаваемый в картинной плоскости изображением движущейся подстилающей поверхности. При использовании уравнения движения изображения для вычисления оптического потока центр тяжести задачи смещается на вычисление параметров движения наблюдаемого объекта относительно ИС. Ниже излагается метод функционализации в приложении к задаче определения полного вектора скорости движения носителя ИС относительно подстилающей поверхности.
1.3. Уравнение функциональной связи (ФС-уравнение)
В плоскости Рк (см. рис. 1) выделим односвязную регулярную область (окно) анализа В с границей Г (В) (рис. 2).
Зададим в окне анализа на множестве функций распределения освещенности Е(т,1) функционал Ф@) вида
Ф = Ц К (г) Е(г (X), X) (10)
В
где К(г) - непрерывная, ограниченная и дифференцируемая почти всюду по всем своим аргументам функция веса;
Е(г(Х),Х) - равномерно ограниченная и дифференцируемая почти всюду по всем своим аргументам функция распределения освещенности изображения яркостного объекта в картинной плоскости Рк.
Вычислим полную производную по времени от функционала (10) в силу уравнения движения изображения (9):
Р и с. 2. Окно анализа
Ф (X) =Ф; (X )+Ц К (г) Е (г(Х), X) Г
В
где Ф'х = Л К (г )Е (г (X), X) X ds■;
(11)
В
С использованием формулы связи двойного и криволинейного интегралов (формулу Грина) из (11) получим
Ф (X) = Фг (X) +11 (X) +12 (X) + Ф1 (X) + ф 2 (X) + Ф 3 (X) + Ф 4 (X), (12)
где
) = - 1К (г ) Е(г^ ), X )Са[(гОз / / + П)(г + к«г, р;
Г( В)
I2 (t) = - f K(r)E(r (t), t)Cakm (rA3 l f + A)VHdl;
Г( D)
F,(t ) = JJ [ к ( r)]rE(r(t)>t)[(rW3 lf + W)(r + kmrs)]dS
D
Ф 2 (t )-fl [ K (r )]Г E(r (t\t)km (rA3 l f + A)VH ds;
D
F3 (t) = "ИK(r )E(r(t)’t)[- 3W3(r + kmrs )/ f + A3TWrs l(A3Trs + ZH )]dS
D
F4(t) = -IfK(r)E(r(t),t)[ГЖ l(A3Trs + Zh )-3^^ /f ]ds;
D
dv r (x)
Ca=[s/'«a; -coso; 0] (a = arc tg----------; уr (x) - локальное уравнение границы
dx
/ / /
T(D) окна анализа D); W3=[a>b ю2, W3]; K(r) r = [K(r) x,K(r) у,0], A3T=[a13, a23, a33].
В дальнейшем будем считать, что (t) = 0 , т.е. внешний источник облучения подстилающей поверхности не меняет своей интенсивности во времени.
Полученное уравнение (12) является ФС-уравнением, которое связывает компоненты вектора Vh поступательной скорости носителя ИС и вектора w угловой скорости относительно подстилающей поверхности с временными и пространственными характеристиками изображения яркостного объекта в картинной плоскости ИС.
Значения интегралов Ii(t) и I2(t) должны вычисляться на границе окна анализа T(D). Точное измерение аппаратными средствами величин, соответствующих этим интегралам, не представляется возможным, так как они определены на множестве нулевой меры, а аппаратными средствами физически может быть измерена интегральная освещенность только областей ненулевой площади. Поэтому способы определения параметров движения, вытекающие из ФС-уравнения (i2), будут корректными лишь при условии, что функция веса K(r) в функционале имеет нулевое значение на границе T(D) окна анализа.
Примером такой функций веса для случая окна анализа прямоугольной формы является «пирамидальная» функция (рис. 3). На рисунке координаты (X; у) отсчитываются от точки пересечения диагоналей окна анализа:
"(x, у): x е [a,-a], у е [- kx, kx] ® f (x, у) = u(l - |x| l a);
"(x,у): уе [b,-b], xе [-уlk, уlk]® f(x,у) = u(l-|у|lb),
где u=a/b.
у \b
Л D3\ D2/"' a / t
/D4 \Di x
4-b
б
Р и с. 3. «Пирамидальная» (а) функции веса и области D1...D4 (б) в окне анализа
Из анализа ФС-уравнения (12) следует, что оно приводится к линейному относительно вхождения компонент вектора скорости движения ИС алгебраическому уравнению вида
Ь1УХН + Ь2УУН + ЬзУ1Н + Ь4Ю1 + Ь5Ю2 + ЬбЮ3 = ^ , (13)
где УхН, УуН, У%Н, и Ю1, Ю2 и ю - компоненты вектора скорости поступательного и углового движений носителя ИС;
Н=ёФ/Л - полная производная по времени от функционала (10);
Ь - Ь6 - коэффициенты, определяемые значениями интегралов в (12). Они вычисляются из выражений
Ь1 =Ц кт
В
к X ()
X
1
— аз1 + ап
У_
1
а31 + а21
В
Ь2 = Ц кт
В
к X ()
X
1
а32 + а12
+К У ()
У -1
1
а32 + а22
В
Ь3 = Ц кт
В
-я
К X (Г)
X
1'
а33 + а13
У -1
1
а33 + а23
В
1
кта33
+К У (г)
К (г )Е (г , X )•;
у
В
11
+ К У (Г)
V
■и
В
11 К (г )Е (г, X )&;
(14)
=-л
В
к X (Г)
-И
В
3 ( I к Г ) а33XSl - а13^
— (X + kmxs) - • —
1 А3Тг1 + %н
+К У (г) К (г )Е(г, X )&;
1 1
Ьб = -Ц к (г )(У + ктУ1 ) - КУ (г )(x + kmxs )Е(г, *№ +
В
+
Л
В
а23^_+а13 К ( )Е (
А3Тг1 + %н
К (г )Е (г, X )•;
Ь
4
Ь
5
к = —
и К (г )Е (г, X
При использовании функции веса, значения которой равны нулю на границе окна анализа, все коэффициенты в ФС-уравнении (13) оказываются величинами доступными физическому измерению в реальном времени, поскольку они определены на подобластях ненулевой площади (меры). Это является существенным преимуществом метода функционализации в сравнении с градиентным методом.
Система уравнений, необходимая для отыскания компонент вектора скорости
движения Л = [Ух, Уу, У%, ®1, ю2, ю3 ]Т носителя ИС, получается параметризацией уравнения (11). Роль параметра может выполнять время, расположение, количество и конфигурация окон анализа, вид функции веса К(г) и т.п.
При такой параметризации получается система уравнений
ВЛ = Н, (15)
где В=[Ьу] - матрица коэффициентов; Н=[кг]Т - вектор-столбец; 7=1...N(Л>М), ]=М (М=1.. .6) - количество степеней подвижности носителя ИС.
С целью уменьшения погрешности получаемых оценок вектора скорости система ФС уравнений (15) переопределяется (Ы>>М), а решение ее находится методом квазиобращений:
Л = (ВТВ)-1ВТН. (16)
В общем случае метод функциональных преобразований позволяет определять не только вектор скорости Л, но и, при наличии информации о значениях элементов матрицы А в начальный момент времени X0 (А(^)=А0), пространственную ориентацию оси оптической системы. Для этого необходимо дополнить уравнение (16) кинематическим уравнением углового движения носителя ИС:
А(X) = -ПЛ (X), А(^о)=Ао. (17)
2. Влияние случайного аддитивного шума на смещение оценки скорости движения изображения
Исследуем влияние случайного аддитивного шума в изображении на смещение оценок скорости движения, получаемых методом функционализации, для случая двухмерного плоско-параллельного движения носителя ИС относительно подстилающей поверхности при направлении оптической оси ИС в надир.
При этих допущениях функционал (10) будет иметь вид
Ф(0 = Ц К (X, у,-1 ) Е( x(t), У(0,-/, №, (18)
В
где Е(x(t), у^),-1, X) = Е(x(t), у^),-1, X) + к(X, у,-1, X)- доступная измерению функция освещенности изображения (видеосигнал);
Е(x(t), у(^),—1, X)- незашумленная составляющая освещенности изображения; к( X, у,—1, X)- случайная функция (шум).
Будем полагать, что шум к(х, у,-/, X) и изображение Е(х(1), у(Х),-/, X) некорре-лированы и выполняются условия дифференцируемости реализаций случайной функции к( х, у,- /, X) по времени X.
Возьмем полную производную по времени X функционала (18):
Ф = -ихФ'х -иуф'у + я/, (19)
где и х = ^хя /1%я • и у = Ууя /1%я •
ф'х = И К'х (х, у,-/) Е (х(:), у((),-/
в
Ф у = Ц Ку (х, у,-/) Е (х(0, у(?),-/)<*;
у
в
я; =Ц к (х, у,- /) А
Уравнение (19) невозможно использовать для вычисления компонент вектора
Т /
Vх ,иу ,0] потому, что реализация случайной функции Я( недоступна измере-
нию.
Алгоритм определения оценок вектора скорости основан на использовании приближенного соотношения, получающегося из (19) отбрасыванием составляющей Я/ из ее правой части:
Ф(X)=-~хФх -~уфу, (20)
где Фх = \\ к/( х у,-/) Е( х(* X у(4),-/,
в
Ф 'у = Ц к у (х, у,-/) Е( х(0, у(0,-/, 0<*;
в
~х = Vх + АVх и ~у = Vу + АVу - оценки компонент вектора скорости движения изображения (здесь и далее Аuх и Аvу - погрешности вычисления скорости движения изображения).
Для вычисления оценок компонент Vх и ~у вектора скорости и уменьшения
влияния шума на погрешность вычислений воспользуемся методом наименьших квадратов с осреднением данных по конечному множеству окон анализа В, (■е {1...п}) и путем параметризации уравнения (20) сформируем систему уравнений вида
п п п
У (Ф1 ■ )2 + и У (Ф; ф 1 ■ + ф я1 ■ + я! ф 1 ■ + я! я! ■) = -У (Фф1 + Ф я'■ )•
х^^ч^х^ 1 ^хглл уг ллх^уг х1х± уг / хг ^1±лхг^
;=1 ;=1 ;=1 (21)
У (Ф1 ф'■ + ф' я1 + я'ф' ■ + я' я'■) + и У (Ф' ■ )2 =-У (Ф■ ф' ■ + Ф я1)
х vчtxг чtyг 1 ^упАх^ 1±х^у^ 1±х^у^^ у^ “ т)■
■=1 ■=1 ■=1
Из (21) следует, что смещения (Аих и Аиу) оценки математического ожидания
значений компонент вектора скорости ограничены, если выполняется условие отделимости от нуля значения определителя этой системы:
У [ф )2 + (я'Х1 )2] • У [(Ф'у1 )2 + (я'у1 )2] - У (фф )2).. > 5 > 0 . (22)
■=1 ■=1 ■ =1
При отсутствии случайного шума в видеосигнале последнее условие приобретает вид
(23)
У [ф■ )2 ] У [(фУ> )2 ] - У (ф>У> )2 )•• > 5 > 0.
■=1 ■=1 7=1
Это соотношение выделяет класс изображений, для которых метод функциона-лизации позволяет определять скорость движения изображений. Выделение таких классов представляет интересную задачу для дальнейших исследований.
Если изображение таково, что все функции Фу, Ф1у, яу, я у взаимно некорре-лированы, то асимптотическая оценка смещения приобретает вид
1
■и.
Аиу =-и у
1+У (ф у, )2/ У (я у ■ )2
■=1 ■=1
______________1___________
1+ У (Ф у ■ )2/ У (яу, у
(24)
■ =1
■ =1
Из (24) следует, что оценка скорости, получаемая методом функционализации, имеет смещение отрицательного знака, пропорциональное определяемой скорости движения изображения и обратное отношению мощности полезного сигнала к мощности шума. Очевидно, что это отношение увеличивается, а погрешность смещения уменьшается с увеличением полной вариации функции освещенности изображения,
которая характеризуется значениями функций У (Ф!х■ )2 и У (Ф!у■)
и количеством п
■=1
окон анализа. Это смещение будем называть «скоростной» погрешностью метода. Наличием скоростной погрешности в оценке скорости можно объяснить некоторые оптические иллюзии, в частности, оптический эффект, открытый в 1977 г. ху-дожником-графиком Оучи (Н.ОисЫ) [10]. Этот вывод аналогичен выводу, полученному в [4] при анализе погрешностей простого градиентного метода. На рис. 4 представлено изображение, вызывающее иллюзию Оучи. Малые передвижения изображения по сетчатке глаза (рисунок заимствован из [4]) приводят к кажущемуся движению внутрен -ней части рисунка относительно внешней его части. Такая ошибка в оценке скорости может объяс-
■=1
Р и с. 4. Рисунок, вызывающий иллюзию Оучи
няться различием по осям X и У полной вариации изображения внутри центрального круга и вне его.
Из (24) следует также, что погрешность смещения математического ожидания оценки скорости движения, вызванная наличием случайного шума в видеосигнале, может быть уменьшена до нуля, если измеряемую скорость V устремить к нулю, а это возможно при использовании компенсационного метода измерений. Реализовать компенсационный метод можно с применением электромеханического или электронного «подслеживания» фотоприемной структурой движущегося изображения. При этом в качестве эталона скорости движения может быть использована прецизионная система стабилизации скорости движения фотоприемной структуры, создание которой в настоящее время не встречает технических трудностей.
3. Итерационный способ определения параметров движения изображения
подстилающей поверхности
Одним из способов реализации компенсационного метода измерения скорости является способ итерационного совмещения изображений в окнах анализа, выделенных в двух последовательных кадрах изображения.
Для упрощения формулировки задачи совмещения изображений трансформируем ее так. Приемник изображения формирует последовательно, на достаточно малом интервале времени, два кадра изображения подстилающей поверхности. Кадровые изображения (рис. 5) заносятся в память вычислителя.
В каждом кадре изображения выделяют по области анализа. В первом кадре -произвольно, а во втором - в месте, положение которого на изображении может быть прогнозировано, например, по априорной информации о скорости движения изображения. Необходимо найти во втором кадре такое положение окна анализа, изображение в котором «максимально совпадает» с изображением в окне анализа, выделенном в первом кадре. Это есть задача о совмещении изображений двух последовательных во времени кадров изображения.
Расположим в первом кадре окно анализа (ОА10), положение которого относительно изображения движущейся подстилающей поверхности в момент времени X0 показано на рис. 5. Пусть в момент времени tI=t0+Аt (А?1 - интервал измерения) окно анализа во втором кадре займет некоторое положение ОА20. Необходимо переместить окно анализа во втором кадре так, чтобы изображения в окнах анализа в первом и втором кадрах совпадали наилучшим образом.
Предлагается следующая итерационная процедура совмещения изображений. На каждом ■-том шаге по информации, полученной из окон анализа ОА10 и ОА2Ъ согласно методу функционализации формируется система алгебраических уравнений, необходимая для определения оценок составляющих вектора скорости движения
изображения VX1 и . При этом вместо производной по времени функционала Ф используется разность значений этих функционалов в окнах анализа ОА10 и ОА2^ По ним определяются оценки ~х1 = Vx1Аt и ~у1 = Vy1Аt смещения окна анализа ОА2!
относительно изображения в окне ОА10. Переопределенная система ФС-уравнений получается за счет использования дополнительных окон анализа, которые получаются разбиением основного окна на подобласти. В качестве примера на рис. 5 показано разбиение окна анализа в первом кадре на четыре дополнительные подобласти. Аналогичное разбиение производится и во втором кадре. При таком разбиении переоп-
ределенная ФС-система содержит пять уравнений. Решение ее получают методом квазиобращений.
На следующем шаге окно анализа ОА21 смещается относительно изображения на величины - ~у1 и - ~х1, при этом оно (см. рис. 5) займет положение ОА21+1. По
информации, поступающей из окна анализа ОА10 и ОА21 , вновь формируется система алгебраических уравнений, из которых определяются уточненные оценки ~Х2 и ~у2 смещения окон анализа ОА10 и ОА21+1.
Р и с. 5. Процедура совмещения изображений двух последовательных кадров:
ОА; - окно анализа после і-того шага итерации; ОА4 » ОА3 ;
128x128 и 64x64 - размеры окон.
Фактическое начальное смещение изображения: (26; 48)
Данная процедура повторяется до тех пор, пока смещения ~уп и ~хп на некотором шаге не станут меньше заданной величины (пока ~уп и ~хп не войдут в «трубку достижимости» заданного радиуса р и далее ее не покинут).
В символической форме эта процедура может быть записана следующим образом:
і(к+1)= і(к)+ Лі(к), где s={sx, іу}; Лу={Д8х, Л8У}; і(0)=50; &=1,2,... - номер итерации.
На каждой итерации смещение Лу(к) определяется по информации из окон анализа ОА0 и ОА; с применением соотношения, полученного методом функционализации:
Лі(к) = -(Ф'(к))т • (Ф'(к))-1((ФХ(к))тЛФ(к)),
где ЛФ(к) = [Ф2 7 - Ф1г ]; 7=1,2,...;
Ф17 и Ф27 - значения функционалов в окнах анализа ОА0 и ОА;;
Ф Х(к) = 2[Ф1Х + Ф2Х ].
Специальным вопросом является нахождение условий сходимости предложенного итерационного процесса совмещения изображений, который здесь не рассматривается.
4. Численный анализ погрешности итерационного способа
Для получения численных значений оценок погрешностей итерационного способа было применено компьютерное моделирование. Целью моделирования являлась апробация предложенного итерационного способа и исследование влияния случайного шума на
его погрешность. Система моделирования содержала три основных блока: блок моделирования видеосигнала, снимаемого с ПЗС преобразователя изображения, блок вычисления оценок смещения изображения в двух последовательных по времени кадрах изображения и блок совмещения изображений в этих кадрах. При моделировании в сигнал, снимаемый с ПЗС-матриц, аддитивно подмешивался шум с равномерной плотностью распределения. Амплитуда сигнала шума задавалась в единицах младшего разряда (емр). Последовательные кадры, снимаемые с ПЗС, оцифровывались и пересылались в ОЗУ, в области которого формировались окна анализа. Для определения параметров движения изображения использовался функционал с «пирамидальной» функцией веса.
Эксперименты проводились с изображениями (снимками) типа «город», «горы», «море», «облачный покров», полученными при зондировании Земли из космоса. На рис. 6 представлены результаты компьютерных экспериментов по определению зависимости математического ожидания радиуса «трубки достижимости» (р) и доверительного интервала от амплитуды случайного шума Аш при динамическом диапазоне видеосигнала 127 емр для изображения типа «горы» (рис. 7). Начальное рассогласование положения изображений в последовательных кадрах изменялось от 2 до 50 пикселей, а количество итераций при совмещении изображений принималось достаточным для достоверного определения размеров трубки достижимости р.
В качестве примера на рис. 8 приведено поле направлений получаемых оценок векторов скоростей. На рисунке в пикселях заданы координаты смещения положения окна анализа во втором кадре относительно положения окна в первом кадре. Поле построено для изображения типа «город», представленного на рис. 4. Размеры окна анализа - (128x128). На рисунке показаны траектории итерационного процесса совмещения изображений для начальных смещений (26, 48) и (50, -35). В таблице приведены значения оценок вектора смещения на отдельных шагах итерационного процесса совмещения для случая начального смещения (26, 48). Положение окон анализа в кадре 2 на рис. 4 соответствует траектории процесса совмещения
Р и с.7. Изображение типа «горы»
0,8
А 0,7 пикс.
0,6 0,5 0,4 0,3
0 16 32 48 Аш 64
Р и с. 6. Зависимость математического ожидания М и доверительного интервала В0?9б (серый фон) радиуса трубки достижимости р (пикс.) от амплитуды случайного шума Аш
изображений с начальным смещением (24, 48) пикселей. Число итераций, требующихся для достижения совмещения изображений с точностью в 1 пиксель, даже при таких достаточно больших начальных смещениях не превышает трех.
д
Р и с. 8. Поле направлений вектора оценки смещения изображений, представленных на рис. 5
№ шага Сдвиг на шаге итерации Суммарный сдвиг Округленное значение суммарного сдвига
Дх Ду X (Дх) Х(Ду) ~ЦДх) »Х(Ду)
1 -40,57 -17,59 -40,57 -17,59 -41 -18
2 11,07 -35,59 -29,50 -53,19 -30 -53
3 4,26 5,06 -25,24 -48,12 -25 -48
4 -1,00 0,00 -26,24 -48,12 -26 -48
Таким образом, алгоритмы определения вектора скорости движения изображения могут быть реализованы в реальном времени и дают оценки скорости движения с высокой точностью.
Заключение
В статье развит метод функционализации параметров изображения для решения задачи определения параметров свободного движения носителя изобразительной системы относительно подстилающей поверхности. Метод основан на анализе видеоинформации, содержащейся в двух последовательных кадрах изображения, и отличается от известных сочетанием свойств универсальности и высокой точностью. Метод не накладывает ограничений на число степеней свободы движения носителя ИС, не использует таких сложных вычислительных процедур, как поиск экстремума многомерных функций или вычисления сверток, и поэтому является перспективным для разработки «быстрых» способов и алгоритмов определения параметров движения для реализации в системах, работающих в реальном времени.
1. Horn B., SchunckB. Determining Optical Flow // Artificial Intelligence 17 (1981). - P. 185-203.
2. Lucas B.D. and Kanade T. An Iterative Image Registration Technique with an Application to Stereo Vision // IJCAI '81. - P. 674-679.
3. Black M. J. and Anandan P. A framework for the robust estimation of optical flow, ICCV’93. - May, 1993. - P. 231-236.
4. Fermuller C., Shulman D., Aloimonos Y. The Statistics of Optical Flow // Computer Vision and Image Understanding 82. - 2001. - 1-32.
5. Абакумов А.М., Кузнецов П.К., Мишин В.Ю., Семавин В.И. Метод измерения параметров движения объекта // Оптические сканирующие устройства и измерительные приборы на их основе: Тез. докл. всесоюз. совещ., 4-5 июня, 1980. - Ч.1 - Барнаул, 1980. - С. 50-51.
6. Schalkov R.J. Image Motion Analysis Using Concept of Weak Solution to Distributed Parametr Systems //Proc. 1983 IEEE Computer Vision and Pattern Recognition Conf., Washington, DC. - June 1983. - P. 232-239.
7. Кузнецов П.К., Семавин В.И. Метод определения параметров движения яркостного поля // Изв. вузов. Приборостроение. - 1990. - №6. - С. 26-30.
8. Кузнецов П.К., Семавин В.И., Мишин В.Ю., Владимиров М.В. Метод функциональных преобразований в задаче определения скорости движения яркостных полей // Вестник Самар. гос. техн. ун-та.
- 1994. - №1. - С. 66-76.
9. А.с. №1742729 СССР: МКИ G 01 3 3/36. Устройство для определения составляющей вектора скорости движения объекта / В.Е. Агеев, В.Н. Войтенко, О.А. Анайкин, П.К. Кузнецов, В.И. Семавин.
- № 4774743/10. Заявл. 26.12.1989. Опубл. 23.06.1992. Бюл. №23. - 7 с.
10. Spillmann L., Tulunay-Keesey U. and Olson J. Apparent floating motiob ib bormal and stabilized vision, Investigative Ophthal. Visual Sci. Supple. 34. - 1993. - 1031.
11. http://ftp.fi.muni.cz/pub/bibliography/Graphics/flow.estimation.bib.gz.
Статья поступила в редакцию 16 сентября 2008 г.
UDC 681.518.3
METHOD FOR COMPUTING VELOCITY OF MOVING OBJECTS BY IMAGE ANALYSIS
П.К. Кузнецов1, Б.В. Мартемьянов, В.И. Семавин, Е.Ю. Чекотило
1 Samara State Technical University
244, Molodogvardeyskaya str., Samara, 443100
In this paper, we present a very accurate (of image parameters functionalisation ) and optical flow with nonuniform brightness variations. The proposed algorithm is based on a generalized dynamic image model in conjunction with a regularization framework to cope with the problem of non-uniform brightness variations. To alleviate flow constraint errors due to image aliasing and noise, we employ a least-squares method to suppress unreliable flow constraints, thus leading to robust estimation of optical flow.
Key words: optical flow, vector velocity field, motion parameters, functionalization method
1 Pavel K. Kuznesov, Doctor of Technical Sciences, Professor.
Boris V. Martemyanov, Candidate of Technical Sciences, Associate professor. Vladimir I. Semavin, Candidate of Technical Sciences, Associate professor. Elena. Yu. Chekotilo, Assistant.
