CC BY
83
Вычислительные технологии
Том 18, № 5, 2013
Снижение вычислительных затрат с применением алгоритма трёхмерного рекурсивного поиска при построении векторов перемещений в оптическом методе оценки деформации*
С. В. Панин1'2, В. В. Титков1, П. С. Лювутин1 1 Институт физики прочности и материаловедения СО РАН, Томск, Россия 2Национальный исследовательский Томский политехнический университет, Россия
e-mail: svp@ispms.tsc.ru
Предложен подход к снижению вычислительных затрат и одновременному повышению помехоустойчивости построения векторов перемещений, основанный на совместном применении алгоритма трёхмерного рекурсивного поиска (3DRS) и ин-крементного подхода к определению перемещений в серии изображений. Получены оценки вычислительных затрат при комбинированном использовании указанных подходов при обработке модельных и экспериментальных оптических изображений. Проведено сравнение помехоустойчивости их функционирования при определении перемещений в сопоставлении с инкрементным алгоритмом, используемым для определения оптического потока. На основании проведённых исследований рекомендованы параметры вычислений, обеспечивающие значительное снижение вычислительных затрат и увеличение помехоустойчивости построения векторов перемещений.
Ключевые слова: вычислительные затраты, инкрементный алгоритм, вектор смещения, помехоустойчивость.
Введение
Оптический метод оценки деформации, называемый в литературе DIC (digital image correlation), включает алгоритм построения векторов перемещений и последующий расчёт компонент деформации [1]. Большинство исследований в области алгоритмов построения векторов перемещений направлены на повышение точности (увеличение помехоустойчивости) определения смещений [2, 3]. Но не менее важной является и проблема быстродействия алгоритмов. В литературе по обработке изображений и техническому зрению широкое распространение получили многомасштабные подходы, использующие
* Работа выполнена в рамках проекта исследований СО РАН № III.23.1.3, грантов РФФИ № 12-0831042 "Разработка научных основ комбинированного акустико-оптического метода диагностики состояния нагруженных материалов" и № 13-07-00009 "Развитие быстродействующих и помехоустойчивых алгоритмов обработки и анализа оптических и акустических сигналов для комбинированного метода контроля состояния нагруженных материалов", а также гранта № СП-816.2012.5 "Разработка и исследование алгоритмов обработки и анализа оптических сигналов для мониторинга состояния нагруженных деталей авиационной и космической техники" (Стипендия Президента Российской Федерации молодым учёным и аспирантам, осуществляющим перспективные научные исследования и разработки по приоритетным направлениям модернизации российской экономики).
разложение изображения в пирамиду Гаусса и позволяющие заметно сократить время вычислений. Известны также подходы, основанные на применении аппарата параллельных вычислений и многопроцессорных вычислительных систем [4 - 6]. Использование современных видеокарт для обработки данных обеспечивает повышение производительности в 2-12 раз [6]. В то же время известен ряд исследований, посвящённых непосредственно алгоритмической оптимизации вычислений при определении оптического потока. В частности, следует отметить работу [7], в которой решается проблема оптимизации вычислительных затрат для алгоритма Lucas — Kanade [8] путём вынесения повторяющихся вычислений в отдельную процедуру. Кроме того, получил развитие метод трёхмерного рекурсивного поиска (3D recursive search — 3DRS) блочной оценки перемещений, учитывающий пространственную и временную связь соседних векторов смещений [9].
Цель настоящей работы — исследовать возможности применения метода 3DRS при построении векторов перемещений в задаче оценки деформации материалов и предложить модификацию данного алгоритма путём её комбинированного использования с предложенным авторами ранее инкрементным алгоритмом [10]. Проведение сопоставительных исследований помехоустойчивости функционирования указанных алгоритмов проводилось на сериях модельных и экспериментально полученных оптических изображений.
1. Описание методов
В качестве базовых в работе применены были 3DRS (3D Recursive Search) [9] и ин-крементный алгоритмы [10]. Первый метод получил широкое распространение при обработке видеопотоков данных, и с его помощью становится возможным снизить время построения векторного поля в разы [9]. При этом построение каждого вектора смещений производится на основании выбора из совокупности векторов-кандидатов с использованием процедуры поиска минимума меры подобия, в качестве которой используется сумма абсолютных разностей SAD фрагментов изображения — текущего и предыдущего кадров. Метод имеет вероятностный характер: векторы-кандидаты указывают предположительное направление смещения фрагментов изображения. Набор векторов-кандидатов (рис. 1) содержит как временную (вектор смещения данного блока на предыдущем кадре), так и пространственную (вектора в соседних блоках) составляющие, а также векторы обновления, формирующиеся из пространственных векторов [9].
S1
S2 С
т
Рис. 1. Схема набора векторов-кандидатов. С — текущий блок. Я1, Я2 — пространственные кандидаты, Т — временной кандидат; стрелки указывают направление сканирования [9]
В инкрементном алгоритме для устранения ошибок определения перемещений, связанных с формированием деформационного рельефа и других процессов на поверхности материала, предложено оценивать смещение участков соседних изображений в серии с учётом таковых, найденных для предыдущих пар изображений. Таким образом, проводится слежение за перемещением одного и того же участка поверхности, который изменяется во времени. Положение участка оценивается с субпиксельной точностью. Поскольку изменения между соседними изображениями в серии значительно меньше, чем между начальным и текущим изображениями, то помехоустойчивость и точность определения перемещений значительно возрастают. Значения яркости пикселов искомого участка на предыдущем изображении (относительно текущего в серии) определяются с использованием B-сплайн-интерполирования. В данном методе для нахождения смещений в одной паре применяется комбинированный алгоритм [10], использующий классический корреляционный подход для определения перемещений с пиксельной точностью [11] и алгоритм Lucas — Kanade [8, 12] для уточнения перемещения до долей пиксела. Такое сочетание обусловлено прежде всего хорошей помехоустойчивостью корреляционного алгоритма при нахождении больших перемещений и точностью дифференциального алгоритма Lucas — Kanade при нахождении перемещений с субпиксельной точностью. Кроме того, в исследованиях авторов было показано [10], что использование алгоритма Lucas — Kanade в значительно меньшей степени сопровождается проявлением эффекта типа peak-locking, приводящего к заметным ошибкам в определении смещений с субпиксельной точностью. С учётом возможностей, обеспечиваемых каждым из двух описанных подходов, в настоящей работе решались следующие задачи, основанные на модификации "базовых" методов:
— использование в методе 3DRS в качестве меры подобия наряду с традиционно применяемой суммой абсолютных разностей SAD коэффициента корреляции ZNCC;
— привлечение метода 3DRS для реализации инкрементного алгоритма построения векторного поля;
— реализация инкрементного алгоритма определения векторов перемещений с применением метода 3DRS и использование в качестве меры близости коэффициента корреляции ZNCC.
Таким образом, работа основана на сопоставительном анализе вычислительных затрат и помехоустойчивости определения векторов перемещений при использовании следующих алгоритмов:
1) инкрементный алгоритм (I);
2) 3DRS с мерой подобия SAD (3DRS-SAD);
3) 3DRS с мерой подобия ZNCC (3DRS-ZNCC);
4) инкрементный алгоритм с 3DRS/SAD (I-3DRS-SAD);
5) инкрементный алгоритм с 3DRS/ZNCC (I-3DRS-ZNCC).
2. Описание исследуемых изображений
Оценка точности и помехоустойчивости построения векторов перемещений, а также вычислительных затрат исследованных алгоритмов проводилась на модельных и экспериментально полученных изображениях. Для оценки ошибки построения векторов перемещений на модельных изображениях использовались соответствующие модельные векторные поля.
2.1. Серии модельных изображений
Формирование исходного изображения поверхности. Модельное изображение получалось из заданного количества слоёв псевдослучайных отсчётов яркости; при этом каждый слой соответствовал определённой пространственной частоте [13]. Подобно описанию, приведённому в [10], имея начальный слой размером 4 х 4 пиксела, после проведения восьми итераций было получено модельное изображение 1024 х 1024 пикселов (рис. 2, а).
Двуосное растяжение (серия 1). С целью моделирования изменений, происходящих при нагружении поверхности по схеме двуосного растяжения, задавалось смещение каждой точки модельной поверхности (рис. 2, б). При этом яркость каждого пиксела изображения пересчитывалась для заданного приращения деформации. Пересчёт производился с помощью В-сплайн-интерполирования. В результате из начального изображения формируется вся серия, состоящая из 25 кадров и отражающая деформацию по схеме двуосного растяжения с конечным приращением 5 = 80 пикселов.
Зашумление изображений (серия 2). Для формирования данной серии были взяты первый и последний кадры серии 1. Далее на каждый из кадров накладывался гауссов шум с изменением амплитуды ап от 1 до 25 % с шагом 1 %. В результате получена серия
а
в
г
Рис. 2. Модельное исходное (а) и фрагментированное (в) изображения и модельное векторное поле для двуосного растяжения (б) и изображения на рисунке в (г) (серия 4)
зашумленных изображений, содержащая 25 пар кадров.
Двуосное растяжение и изменение рельефа поверхности (серия 3). Для одновременного учёта изменений, связанных с двуосным растяжением и изменением рельефа, была создана серия, в которой каждое изображение формируется по схеме, описанной в [10]. Расчёт значений пикселов изображений серии производится по формуле
Р = (1 - к)Р + кР2,
где Р\, Р2 — значения пикселов двух исходных изображений, к — весовой коэффициент, изменяющийся от 0 до 0.5 с шагом, равным обратному значению количества изображений в серии.
Изображения из четырёх фрагментов с разрывами поля перемещений (серия 4). Важным тестом быстродействия и помехоустойчивости построения векторных полей является анализ изображений, отражающих разрывы в поле смещений. С этой целью было промоделировано двуосное растяжение изображений, состоящих из четырёх фрагментов. Каждое изображение серии формировалось из четырёх фрагментов, полученных способом, описанным для серии 1. Изображения компоновались (рис. 2, в) таким образом, чтобы обеспечить максимальное количество вариантов комбинаций направлений векторов на границах фрагментов (рис. 2, г).
2.2. Серия экспериментальных изображений
Серия экспериментальных оптических изображений была получена при растяжении образцов из порошка меди, спечённого в вакууме электронным лучом (рис. 3). Спецификой анализируемых изображений было существенное изменение рельефа на поверхности при деформировании, что накладывало особые требования к алгоритмам построения векторов перемещений. Регистрацию изображений проводили по методике, описанной в [14].
а
Рис. 3. Оптическое изображение образца (а) и векторное поле (б), полученное из экспериментальных изображений
3. Методика тестирования алгоритмов
Быстродействие алгоритмов оценивалось по удельному времени расчёта t одного вектора, т. е. по отношению полного времени построения векторного поля к количеству векторов в этом поле. Для проведения расчётов использовалась ПЭВМ со следующими характеристиками: процессор Intel(R) Core(TM) i3 CPU M 350, объём оперативной памяти 2.00 Гб, операционная система Windows 7. Для исследования быстродействия алгоритмов были использованы описанные в предыдущем разделе серии изображений. Для количественного сравнения помехоустойчивости определения перемещений с применением тестируемых алгоритмов использовали коэффициент корреляции полей расстояний векторных полей [15, 16], рассчитываемый по формуле
2 wh
w h ____w h ___
E E (A* - Dx) {Dxpi. - Dxp) £ E (Dj - Dy) (D ypij Dyp) i=lj=l + i=lj=l
°DX &Dxp
aDy О Dyp
1
k
r
где Dx, Dxp, Dy,Dyp — поля расстояний соответствующих векторных полей, которые необходимы для перехода от двухкомпонентных данных векторного поля к одноком-понентным; w, h — ширина и высота поля; Dx, Dxp, Dy,Dyp — средние значения соответствующих полей расстояний; а — среднеквадратичная ошибка поля расстояний. В нашем случае каждый элемент поля расстояний находили следующим образом:
w h
Dij = - xki\ + \Vij - Vki\), г =l,...,w, j = l,...,h.
k=1 1=1
4. Результаты тестирования
4.1. Тестирование на модельных изображениях
Вычислительные затраты. Анализ времени расчёта одного вектора показал сохранение примерно на постоянном уровне значения t для всех алгоритмов вне зависимости от величины варьируемых параметров (рис. 4, а-г) для всех исследуемых модельных серий изображений. Для инкрементного I алгоритма наблюдается снижение t в конце серий 1 и 4 (рис. 4, а, г), что связано с уменьшением количества определяемых векторов в поле вследствие выхода искомых участков за пределы изображений при больших степенях деформации.
Для всех исследованных серий изображений самое малое значение t — ~ в 10 раз меньшее, чем для инкрементного метода I, зафиксировано для 3DRS-SAD и I-3DRS-SAD. В свою очередь алгоритмы, использующие ZNCC, характеризуются промежуточным значением t, поскольку они обладают как преимуществом 3DRS-алгоритма с ограниченным числом кандидатов (полный перебор не производится), так и его недостатком, связанным с большими (по сравнению с SAD) временными затратами расчёта величины ZNCC.
Время расчёта всего векторного поля (рис. 4, д) для 3DRS-SAD в 7 раз меньше, чем для метода I. Данный показатель не зависит от размера векторного поля (полного числа векторов в поле).
О 10 20 30 40 50 60 70 & пике в
Рис. 4. Зависимость времени расчёта одного вектора от величины деформации (серия 1) (а), амплитуды шума (серия 2) (б), весового коэффициента (серия 3) (в), приращения двуосной деформации фрагментированного изображения (серия 4) (г) для инкрементного I (1), ЗБИЯ-ЯЛБ (2), ЗБКЯ^СС (3), ^ЗБКЯ-ЯЛБ (4), ^ЗБИЯ^СС (5) алгоритмов; зависимость времени расчёта полного векторного поля от количества векторов в поле (д)
а
г
Помехоустойчивость определения перемещений. С целью количественного сравнения помехоустойчивости определения перемещений исследованными алгоритмами использовались векторные поля, построенные при обработке модельных серий изображений. Для этого рассчитывались значения коэффициента корреляции Кг полей расстояний векторных полей (рис. 5) для исходного и "деформированных" изображений подобно тому, как описано в разделе 3.
Для серии 1 характерна высокая "точность" построения векторных полей для всех используемых алгоритмов — ошибка не превышает 5 • 10-4, что связано с отсутствием
помех на изображениях данной серии, поэтому график Кг для данной серии не приводится.
Для серии 2 использовались только инкрементный I, ЗБКЯ-ЯАБ и 3DRS-ZNCC алгоритмы. Алгоритмы I-3DRS-SAD и I-3DRS-ZNCC в данном тесте не применяли, так как их результаты должны совпадать с таковыми, полученными для 3DRS-SAD и 3DRS-ZNCC, поскольку каждая пара кадров в серии рассматривается независимо. Видно (рис. 5, а), что зашумление не влияет на результаты функционирования 3DRS-SAD и 3DRS-ZNCC алгоритмов, а корректность определения векторов смещений при использовании инкрементного алгоритма начинает снижаться при ап = 13 %: на векторном поле появляются локальные некорректно найденные векторы. Таким образом, показано, что из всех рассмотренных алгоритмы 3DRS более устойчивы к влиянию гауссова шума.
График изменения Кг от "весового коэффициента" для серии 3 приведён на рис. 5, б. Видно, что алгоритмы 3DRS-SAD и 3DRS-ZNCC, не основанные на инкрементном подходе, не позволяют достоверно построить все векторное поле: для 3DRS-SAD величина Кг начинает снижаться при к > 0.1, в то время как для 3DRS-ZNCC — при к > 0.28, что связано с невозможностью установить соответствие областей изображений при существенном изменении рельефа поверхности. В то же время инкрементные алгоритмы показывают значение Кг, близкое к единице.
0 10 20 30 40 50 60 70 пике 0 10 20 30 40 50 60 70 5, пике
Рис. 5. Зависимость коэффициента корреляции векторных полей от амплитуды шума (серия 2) (а), весового коэффициента (серия 3) (б), приращения двуосной деформации фрагменти-рованного изображения (серия 4) (в, г) для инкрементного I (1), ЗВКЯ-ЯАБ (2), 3DRS-ZNCC (3), I-3DRS-SAD (4) и I-3DRS-ZNCC (5) алгоритмов
Снижение Kr характерно для всех алгоритмов при тестировании серии 4 (рис. 5, в). Малое значения Kr в начале серий для инкрементного алгоритма (см. кривую 1) связано с наличием ошибочно построенных векторов, которые по абсолютной величине много больше среднего значения длины корректно определённых векторов по всему полю; к середине серии количество ошибочных и корректных векторов становится примерно равным. Общее снижение Kr для всех алгоритмов связано с ростом величины смещений вдоль границ фрагментов (разрывов поля перемещений). Алгоритмы I-3DRS-SAD и I-3DRS-ZNCC показали наилучшую устойчивость к наличию разрывов в векторном поле. При этом наибольшая "точность"/помехоустойчивость при малых деформациях наблюдается при использовании SAD, в то время как в конце серии (при больших деформациях) наибольшая "точность" обеспечивается при использовании ZNCC (рис. 5, г).
Таким образом, наиболее помехоустойчивыми являются алгоритмы I-3DRS-SAD и I-3DRS-ZNCC.
4.2. Тестирование на экспериментальных изображениях
Результат исследований алгоритмов на экспериментальной серии изображений (см. рис. 3, а) показал, что для алгоритмов, использующих ЗБИБ, время расчёта £ меньше, чем для инкрементного алгоритма. Характер зависимостей параметра £ (рис. 6, а) аналогичен таковому для модельных серий (см. рис. 4). Для сравнения "точности/помехоустойчивости алгоритмов (рис. 6, б) в качестве "истинных" были взяты результаты, полученные с применением Т-ЗБКЯ-ЯАБ алгоритма (на основании данных тестирования алгоритмов на модельных изображениях, приведённых в разделе 4.1). Видно, что инкре-ментный алгоритм значительно уступает алгоритмам, реализованным с привлечением ЗБИ^-подхода, а рост значения Кг аналогичен таковому для данного алгоритма в случае модельной серии 4 (см. рис. 5, г). Графики для неинкрементных ЗБИ^-алгоритмов также имеют характер изменения, схожий с изменениями для данных алгоритмов в случае модельных серий 3 и 4 (см. рис. 5, в, г): с ростом варьируемого параметра наблюдается постепенное снижение Кг. Таким образом, результаты тестирования алгоритмов на экспериментальных и модельных изображениях согласуются.
а б
и с 0.006 0.005
0.001 0.000
Рис. 6. Зависимости времени расчёта одного вектора Ь (а) и коэффициента корреляции полей расстояний векторных полей (б) от степени деформации для инкрементного I (1), 3ВК8-8ЛБ (2), 3БК8^00 (3), I-3DRS-SAD (4) и I-3DRS-ZNCC (5) алгоритмов по результатам обработки серии экспериментальных изображений
Заключение
В работе предложены модификации алгоритма 3DRS построения полей векторов перемещений, основанные на привлечении инкрементного подхода и использовании коэффициента корреляции в качестве меры подобия. Проведённые исследования на модельных и экспериментальных изображениях позволили сформулировать следующие выводы.
1. Показана эффективность применения 3DRS с позиции уменьшения временных вычислительных затрат. В этом случае время, затрачиваемое на определение одного вектора, может быть уменьшено примерно в 10 раз при использовании I-3DRS-SAD по сравнению с инкрементным алгоритмом I.
2. Установлены зависимости коэффициента корреляции полей расстояний векторных полей Kr, позволяющие получить количественную меру качества получаемых векторных полей, от параметров расчёта. Увеличение Kr в наибольшей степени проявляется для 3DRS алгоритмов для серий изображений, отражающих приращение деформации в условиях большого зашумления (an > 13 %), при большом изменении рельефа (k > 0.1-0.28), а также наличии разрывов в векторном поле для инкрементных 3DRS алгоритмов. Характер полученных зависимостей сохранился и при обработке экспериментальных данных.
Таким образом, из исследуемых комбинаций наилучшей по вычислительным затратам и помехоустойчивости для нахождения полей смещений в оптическом методе оценки деформации является сочетание инкрементного подхода, основанного на применении алгоритма 3DRS, с мерой подобия SAD. Последняя позволяет существенно уменьшать вычислительные затраты при одновременном увеличении "точности"/помехоустойчи-вости построения векторов смещений по сравнению с применением отдельно инкре-ментного и 3DRS методов.
Список литературы
[1] Дерюгин Е.Е., Панин В.Е., Панин С.В., Сырямкин В.И. Способ неразрушающего контроля механического состояния объектов и устройство для его осуществления. Пат. РФ № 2126523 // Бюлл. изобретений № 5. 20.02.99.
[2] Xu L., Jia J., Matsushita Y. Motion detail preserving optical flow estimation // PAMI. 2012. Vol. 34, No. 9. P. 1744-1757.
[3] Sun D., Sudderth E., Black M.J. Layered segmentation and optical flow estimation over time // Computer Vision and Pattern Recognit. (CVPR). IEEE. 2012. P. 1768-1775.
[4] Tao M., Bai J., Kohli P., Paris S. Simple flow: A non-iterative, sublinear optical flow algorithm // Comput. Graph. Forum. 2012. Vol. 31, No. 2, pt 1. P. 345-353.
[5] Marzat J., Dumortier Y., Ducrot A. Real-time dense and accurate parallel optical flow using CUDA // Proc. on 17th Intern. Conf. in Central Europe on Computer Graphics, Visualization and Computer Vision (WSCG). Czech Republic, Univ. of West Bohemia, 2009.
[6] Durkovic M., Zwick M., Obermeier F., Diepold K. Performance of optical flow techniques on graphics hardware // ICME. 2006. P. 241-244.
[7] Rav-Acha A., Peleg S. Lucas —Kanade without iterative warping // Proc. of IEEE Intern. Conf. on Image Processing. 8-11 Oct. 2006. P. 1097-1100.
[8] Bruce D. Lucas Generalized Image Matching by the Method of Differences. Doctor. Diss. Tech. Report. Robotics Institute, Carnegie Mellon Univ., 1984. 167 p.
[9] Braspenning R.A., de Haan G. True-motion estimation using feature correspondence // SPIE. Proc. of VCIP. Jan. 2004. P. 396-407.
[10] Панин С.В., Титков В.В., Лювутин П.С. Инкрементный подход к определению перемещений при построении векторных полей // Автометрия. 2013. Т. 49, № 6 (в печати).
[11] Лювутин П.С., Панин С.В. Исследование точности и помехоустойчивости построения векторов перемещений при оценке деформаций оптико-телевизионным методом // Вычисл. технологии. 2006. Т. 11, № 2. С. 52-66.
[12] Lucas B.D., Kanade T. An iterative image registration technique with an application to stereo vision // Proc. of the 7th Intern. Joint Conf. on Artificial Intelligence. 1981. Vol. 2. P. 674-679.
[13] Панин С.В., Алтухов Ю.А., Лювутин П.С. и др. Применение фрактальной размерности для оценки изображений поверхности, получаемых различными датчиками // Автометрия. 2013. Т. 49, № 1. С. 42-49.
[14] Панин С.В., Бяков А.В., Лювутин П.С. и др. Многомасштабный метод изучения деформации и разрушения нагруженных твёрдых тел по данным акустической эмиссии, корреляции цифровых изображений и тензометрии // Заводская лаборатория. Диагностика материалов. 2011. Т. 77, № 9. С. 50-59.
[15] Восковойников Ю.Е., Белявцев В.Г. Нелинейные алгоритмы фильтрации векторных сигналов // Автометрия. 1999. Т. 35, № 5. С. 48-57.
[16] Панин С.В., Титков В.В., Лювутин П.С. Исследование методов фильтрации векторных полей в задаче оценки деформации материалов методом корреляции цифровых изображений // Там же. 2013. Т. 49, № 2. С. 57-67.
Поступила в 'редакцию 15 января 2013 г.
с доработки — 20 мая 2013 г.
