# dspa2011 5 3

Document Sample

```					                              Обработка и передача изображений
____________________________________________________________________________________________

2. Петров Е.П. Метод моделирования многомерных многозначных марковских процессов/ Е.П.Петров,
И.С.Трубин, Н.Л.Харина// Сб. трудов XII Международной конференции «Радиолокация, навигация,
связь-RNLC», Воронеж, Т. 1, 2006.- с. 122-128.

MATHEMATICAL MODELS OF VIDEOIMAGES BASED ON MULTIDIMENSIONAL MARKOV CHAINS
Medvedeva E., Harina N., Metelyov A.
Vyatka State University, Kirov

Method of synthesis of mathematical models for static and dynamic digital grayscale images based on multidi-
mensional Markov chains is offered. Digital grayscale images are g-bit binary numbers, which form the multidimen-
sional multi-valued Markov random process with2 g equally probable states and matrices of transition probabilities
from one state to another size 2  2 . Synthesis of multidimensional ( Q -dimensional) multiple-valued mathe-
g   g

matical model adequate to real digital grayscale images is divided into synthesis of mathematical model of bit binary
Q -dimensional Markov chain with two equally probable states M 1  , M 2  and
l        l
images, each of which represents a

a matrix of transition probabilities for each measurement r П  r ijl 
                   , r  1, Q , i, j  1, 2 , l  1, g .
22

To obtain the algorithm of the binary states of the elements of the                  Q -dimensional Markov chain is applied en-
tropy approach. It is as follows. Conditional entropy of the element                D  i,lj,k ,
l 
,r   with respect to the values of

   we define as the difference between the unconditional entropy of the element
l
the elements of a neighborhood           Q

D and the mutual information derived from elements of the neighborhood    :
l                                                                                                  l
Q
                                
       l    l  , ,  l   
  w D i                  r 


H  D  1  ,
  l      l      
l
   
, r   H  D   I  1  ,
   l      l              l 
, D           
  log 

                 N
  (1)
 ,

                                
  w  Dl   i l  , , r l   
                  
                                
                                
                      S         
where N  2k  1 ,     S  2m ;
for even   Q : k  0, Q  2 , m  1, Q ; for odd Q : k  0, Q  1 , m  1, Q  1 .
2              2                        2               2
The products marked with                   in (1) are computed for all possible non-matching combinations of different
subscripts                                                                             
Q -dimensional random field ( Q  2 ); w  Dl   il  , , rl  - density of the transition probabilities in


multidimensional Markov chains.
On the basis of the produced Q -dimensional mathematical model the example of obtaining three-dimensional
mathematical model adequate dynamic digital grayscale images is considered. Adequacy of the developed mathe-
matical models to real dynamic digital grayscale images is confirmed by estimations of elements of the matrix of
transition probabilities for artificial and real images.


ОЦЕНКА КАЧЕСТВА МЕТОДА СЕГМЕНТАЦИИ ИЗОБРАЖЕНИЙ НА ОСНОВЕ ДВУМЕРНЫХ
ЦЕПЕЙ МАРКОВА
Медведева Е.В., Курбатова Е.Е.
Вятский государственный университет, г. Киров

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

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              151
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

Существует достаточно большое количество алгоритмов сегментации изображений. Наиболее общим
способом нахождения контуров является обработка изображения с помощью скользящей маски (или серии
масок) с общим числом коэффициентов более четырех, вычисления характерной величины (первой, второй
производной, контрастности и т.д.), которая может принимать широкий спектр значений, и сравнения ее с
порогом.
Для уменьшения количества вычислительных затрат, приходящихся на один элемент изображения,
предлагается использовать для представления разрядных двоичных изображений (РДИ) в цифровых
полутоновых изображениях (ЦПИ) двумерную цепь Маркова с двумя равновероятными состояниями и
матрицами вероятностей перехода размером 2  2 , что позволит определить контурный элемент по двум
соседним элементам.
Другим важным показателем, влияющим на выбор алгоритма, является показатель качества его работы,
который определяется точностью определения контуров объектов интереса на изображениях. Таким
комплексным показателем может служить критерий среднеквадратичной ошибки в совокупности с
критерием минимума эмпирического расстояния между идеальным (эталонным) контурным изображением
и контурами, полученными в результате сегментации.
Поэтому разработка алгоритма сегментации, требующего небольших вычислительных ресурсов и оценка
качества его работы, позволяющая обоснованно выбрать алгоритм для конкретного приложения является
актуальной задачей.
В работе предложен алгоритм сегментации ЦПИ, на основе двумерных цепей Маркова, требующий
небольших вычислительных затрат. Проведена сравнительная оценка качества работы разработанного
алгоритма по предложенным критериям с известными методами.
Метод выделения контуров
Пусть ЦПИ представлено набором из g РДИ. Если l -е РДИ представляет собой случайное марковское
поле с разделимой автокорреляционной функций, то его можно представить как суперпозицию двух
одномерных цепей Маркова с двумя равновероятными значениями
(
M1(l ) , M 2l )           и матрицами вероятностей
переходов по горизонтали и вертикали соответственно [1]:
2                                        l 
1 l
11                12
1l                                   l
11
2       2  l
12 , l  g
1
П                                                  , 2 П l                                          (1)
1  l
 21                 l
 22
1
 21
2 l    2 l 
 22
На          рис.            1        приведен                         фрагмент                двумерного      бинарного
l 
изображения. Элемент                                             3 зависит от двух соседних элементов,
3
1
входящих в окрестность                                                                     
i ,l j   1l  , 2l  . 
Рис. 1. Фрагмент РДИ
Для выделения контуров определим величину количества информации в элементе двоичного
изображения  3 относительно элементов окрестности                                                 i , j :
l

   l 
I  3  2 , 1l     l 
   log                

w  3   1  w  3   2 
l           l
         l         l
,                             (2)
w             l 
3
l 
 ,
2
l 
1      
                      
где w  3l   1l  , w  3l   2l 


    - плотности вероятностей перехода между соседними значениями по

горизонтали и вертикали РДИ; w                         l 
3

 2l  , 1l      - плотность вероятности перехода в двумерной цепи
Маркова.
Для элементов l -го РДИ количество информации между элементом  3 и различными сочетаниями значений

   можно определить по матрице П:
l
элементов окрестности             ij
 1 2   3
1         
1
00  0
             ,
П 2            2                           01  0                                                                                       (3)
3            
3                           10  0
4        
4                           11  0

____________________________________________________________________________________________
152                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

где  i  1   i ,   i  1, 4 . Элементы матрицы П вычисляются по формулам:
1 l  2 l 

1   4     l  2 l       2   3   ii  ij ;

1
 ii l  ii ;
3
 ii

 ij 
3 l
(4)
1 l  2 l                     1 l  2 l 
 3   2 =  ij  ii ;
                         4  1 =  ij  ij ,
 ijl 
3                             3     
 ii 
l

где
3
 ii , i, j  1, 2; i  j элементы матрицы 3 П  1 П  2 П  3 iil 
                           .
22
l                 l 
В общем случае полагаем, что матрицы П                 П априорно известны и определяют
1                 2
и
корреляционные связи между двоичными элементами l -го РДИ по горизонтали и вертикали.
Количество информации в элементе РДИ будет минимально, если окрестные элементы                    1,  2   имеют
знаки одинаковые с  3 [1,2].
В случае появления на РДИ областей другой яркости, на границе области один или два окрестных
элементов будут иметь разные с  3 знаки, и количество информации в элементе  3 увеличится.
Сравнивая значения вычисленной величины количества информации в элементе изображения с порогом,
определяем, будет ли данная точка являться точкой контура.
Значение порога H вычисляют для каждого РДИ с учетом вычисленного минимального количества
информации и количества информации, когда хотя бы один из элементов окрестности будет иметь другие
I ( 4  M 1  1  M 1 , 2  M 1 )  I ( 4  M 1  1  M 1 , 2  M 2 )
знаки:      H                                                                               .         (5)
2
Толщина линии контура на одном РДИ будет составлять 1 элемент.
Для восьмиразрядного ЦПИ, представленного 256 уровнями яркости, старшему РДИ будут
соответствовать 128 уровней яркости. Поэтому по старшему РДИ можно выделить все светлые области с
яркостью от 128 до 255 на темном фоне, либо, наоборот, все темные объекты - на фоне с яркостью выше
128. Для выделения менее контрастных областей или объектов с нечетко выраженными границами
необходимо выделить контуры на следующих РДИ (7-м, 6-м или 5-м). Контурное изображение в этом
случае будет представлять сумму контурных изображений нескольких РДИ. Младшие РДИ (при
 ii  0,5 ) будут составлять фон изображения в виде двумерного шума.
Оценка качества алгоритма сегментации
В результате сегментации возникают ошибки двух типов: на сегментированном изображении точка
отмечена как контурная, а на идеальном контурном изображении она не относится к контуру; на
сегментированном изображении точка не отмечена как контурная, но она является таковой на идеальном
контурном изображении. Поэтому при оценке качества сегментации изображений выбрано два критерия:
критерий, показывающий степень сходства сегментированного и идеального контурного изображений
(FOM) и критерий, показывающий степень их отличия (RMS).
Критерий FOM (Figure of Merit) соответствует эмпирическому расстоянию между идеальным контурным
изображением, представленным в виде контуров f и контурами, полученными в результате сегментации
card ( g )
1                          1
g [3]: FOM ( f , g )                                             ,         (6)
max{card ( f ), card ( g )} i 1 1  d 2 (i )
где card ( f ) - количество пикселей в множестве f , card ( g ) - количество пикселей в множестве g ,
d (i ) - расстояние между i -м пикселем f и ближайшем к нему пикселем в g .
Критерий RMS (root mean squared error) представляет собой среднеквадратичную ошибку, определяемую
1
    1                             22
выражением [3]:            RMS ( f , g )                f ( x)  g ( x)   ,            (7)
 card ( X ) xX
                                   

где f ( x) , g ( x) - интенсивность пикселей x в f i и g i , X – множество пикселей на сегментируемом
изображении.
Для оценки качества разработанного алгоритма проводилось его сравнение с известными методами:
Робертса, Превитта, Собела и Канни. Для моделирования работы известных алгоритмов сегментации
использовалась среда MatLab. Пороговый уровень чувствительности при обработке ЦПИ известными
____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              153
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

методами устанавливался автоматически. В качестве тестовых ЦПИ и идеальных контурных изображений
(заданных экспертами) использовалось 100 ЦПИ из базы изображений университета Беркли [4]. Все
контурные изображения предлагаемым методом были получены по старшему РДИ ЦПИ.
Пример тестового ЦПИ (а), идеального контурного изображения (б) и результаты сегментации
разработанным алгоритмом (в), Робертса (г), Превитта (д), Собела (е) и Канни (ж) приведен на рис. 2.
Из приведенных на рис.2 данных следует, что разработанный алгоритм обеспечивает наиболее точное
выделение контуров, по сравнению с известными методами, что подтверждается визуальным анализом
полученных изображений и вычисленными критериями FOM и RMS.
Кроме того, разработанный алгоритм обеспечивает замкнутость контура, толщиной в один пиксел, что
является важным для следующего этапа сегментации - закрашивании объектов интереса.

а)                                     б)                     в) FOM=0,4596; RMS=0,0957

г) FOM=0,3234;              д) FOM=0,2496;           е) FOM=0,2459;                   ж) FOM=0,1091;
RMS=0,1234                  RMS=0,1361               RMS=0,1371                      RMS=0,2476
Рис.2. Пример сегментации изображения
В таблице 1 приведены усредненные оценки критериев FOM и RMS, полученные для 100 ЦПИ, для
разработанного и известных методов сегментации.
Таблица 1.
Критерий       Разработанный метод        Оператор          Оператор          Оператор            Метод
Робертса          Превитта           Собела             Канни
FOM                      0,2504               0,2380            0,2071            0,2073             0,1902
RMS                      0,2465               0,2225            0,2451            0,2459             0,3385

Анализ полученных результатов показывает, что разработанный метод сегментации на основе
вычисления количества информации по критерию FOM превышает известные методы на 5-25%, а по
критерию RMS не значительно уступает методам Робертса, Превитта и Собела.
Таким образом, разработанный алгоритм сегментации ЦПИ на основе двумерных цепей Маркова прост в
реализации, требует небольших вычислительных затрат и по точности определения контуров объектов
интереса на изображениях практически не уступает известным: Робертса, Превитта, Собела и Канни.
Литература
1. Трубин И.С. Метод моделирования цифровых полутоновых изображений / И.С. Трубин, Е.В.
Медведева, О.П. Булыгина // Инфокоммуникацинные технологии, Том 6, №1, 2008 – C.94-99.
2. Петров Е.П. Вычисление статистической избыточности статических изображений / Е. П. Петров,
Медведева Е.В. // Вопросы радиоэлектроники, сер. РЛТ, 2008, вып.3 – Москва, 2008. – С.76-83.
3. Левашкина, А.О. Сравнительный анализ супервизорных критериев оценки качества сегментации
изображений / А.О. Левашкина, С.В. Поршнев // Информационные технологии. – 2009. - № 5. – с. 52-57
4. Berkeley Segmentation Dataset – база изображений университета Беркли, – URL:
http://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench. – 17.03.2008.

EVALUATION OF QUALITY OF IMAGE SEGMENTATION METHOD BASED ON TWO-DIMENSIONAL
MARKOV CHAINS
Medvedeva E., Kurbatova E.
Vyatka state university, Kirov
The method of segmentation based on contours allocation is offered. The method is based on the representation
of digital half-tone images by two-dimensional discontinuous Markov processes. For contours allocation the amount
of information in the elements of digit binary images is calculated. A comparative evaluation of quality of algo-
rithm’s performance on the proposed criteria with known methods was held.
____________________________________________________________________________________________
154                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

Digital half-tone image is divided into g digit binary images. Given the nature of statistical relationship between
elements in the image, we assume that digital half-tone image is two-dimensional Markov process with several val-
ues. In this case, digit binary image is two-dimensional Markov process with two equiprobable values and matrixes
of transition probabilities of horizontal and vertical.
3                                                           i , j in digit binary image is
l
The amount of information in the element              concerning neighboring elements




calculated by the expression: I  3l   2l  , 1l    log

w  3   1  w  3   2 
  l         l 
     l          l
,                                 (1)
   l 
w  3  2 , 1   l    l 

where q ijl   i, j  1, 2; i  j; q  3 are elements of matrixes of transition probabilities
                                                                                                         1
 - horizontal; 2  - vertical;
   1   2   ; l  g .
3     l       l        l

The amount of information in an element of digit binary image is minimal, if neighboring elements  1 ,  2 have
the same signs with   3 . One or two neighboring elements will have different signs with  3 , if there are areas of
other brightness on digit binary image. In this case, the amount of information in the element  3 is increased. The
affiliation of a point to contour is determined by comparing calculated amount of information in the element of im-
age with the threshold.
To evaluation of quality of the developed method, it was compared with known methods (Roberts, Prewitt,
Sobel and Canny) by selected criteria FOM and RMS.
Analyses of the results showed that the developed method of segmentation based on the calculation of the
amount of information on the criterion FOM exceeds the known methods for 5-25%. It is slightly inferior to
methods of Roberts, Previtt, Sobel on the criterion RMS.


ВЫДЕЛЕНИЕ КОНТУРА КАПЛИ ЖИДКОСТИ В ЗАДАЧЕ ОПРЕДЕЛЕНИЯ ПОВЕРХНОСТНОГО
НАТЯЖЕНИЯ
Мизотин М.М.1, Крылов А.С.1, Проценко П.В.2
1
МГУ имени М.В. Ломоносова, факультет ВМК
2
МГУ имени М.В. Ломоносова, Химический факультет

Введение
Поверхностное натяжение является одним из важнейших свойств жидкости, и его точное измерение
является необходимым для изучения различных явлений и разработки технологических процессов.
Существует целый ряд способов измерения поверхностного натяжения, однако среди всех них можно
выделить метод лежащей или висящей капли. Основные достоинства метода заключаются в очень широкой
области применения – от легких текучих жидкостей до жидких металлов, и относительная простота
экспериментальной установки по сравнению с другими методами. Причем, в связи с развитием цифровой
вычислительной и фототехники стало возможным производить анализ практически мгновенно.
Суть метода состоит в следующем: капля помещается на горизонтальную подложку (метод лежащей
капли) или подвешивается на капиллярной трубке (метод висящей капли) и затем изучается ее фотография в
профиль. Измерение геометрических параметров равновесной капли, форма которой определяется
соотношением плотности и поверхностного натяжения жидкости, позволяет восстановить искомое
поверхностное натяжение. Схема установки представлена на Рис. 1.

Рис. 1. 1 – источник света (лампа или зеркало микроскопа), 2 – капля на подложке,
3 – микроскоп с цифровой камерой.
Несмотря на достаточно хорошо разработанную экспериментальную методику, до сих пор требуется
специальная дорогостоящая установка для съемки капли. В данной работе предложен алгоритм для
экспериментальной установки из широкодоступных компонентов. Недостатки установки по сравнению с
лабораторным оборудованием компенсируются предложенными методами обработки изображений.

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              155
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

Метод лежащей капли
Основное уравнение метода лежащей капли – уравнение Юнга-Лапласа, описывает поверхность капли с
симметрией вращения на горизонтальной подложке. Для решения этой задачи была предложена
эффективная методика [1], впоследствии улучшенная и дополненная [2-4].
Данная методика основана на численном дифференцировании уравнения Юнга-Лапласа. Для того, чтобы
продифференцировать уравнение Юнга-Лапласа вводится параметризация кривой {x(t ), y(t )} , где t –
длина дуги кривой от вершины капли (Рис. 2).

Рис. 2. Параметризация контура капли.
Эта параметризация удовлетворяет условию x  y  1 , и приводит к системе уравнений
2  2
 y  y0       
y    2
   x 
y                         ;
  2      x  x0 R0            (1)
 y  y0       
y     2
  y 
x                        ,
          x  x0 R0 
2

с начальными условиями y (0)  y 0 ,    x(0)  0 , x(0)  1 , y(0)  0 и дополнительным условием
          
y / x t 0  1/ R0 . В разработанном программном пакете, задача Коши (1) решается методом Рунге-Кутта

четвертого порядка точности.
Для восстановления параметров лежащей капли необходимо решить обратную задачу определения
капиллярной постоянной    2,   координат апекса капли {x 0 , y 0 } и ее радиуса кривизны R0 по функции
радиуса горизонтального сечения капли от высоты над подложкой. Эта функция измерена с ошибкой и, в
ряде случаев, доступны измерения только части контура капли. При решении данной обратной задачи
минимизируется ошибка E     d {xi , yi },{x(t ), y(t )}                               (2)
между экспериментальными точками { x i , y i } и кривой {x(t ), y(t )} , полученной в результате численного
решения задачи (2). Разность между экспериментальными точками и кривой определяется как корень из
суммы квадратов расстояний от каждой экспериментальной точки до кривой.
В связи с этим возникает следующая задача обработки изображений: автоматическое получение контура
капли { x i , y i } , что осложняется наличием пыли и мусора на снимках (что связано с применением обычной
камеры в «бытовых» условиях), а также переменными условиями освещения.
Функция ошибки
Одной из основных частей метода является вычисление функции ошибки (2). Вычисление расстояния
между точкой и кривой d {xi , y i }, {x(t ), y (t )}  min ( xi  x(t i )) 2  ( y i  y (t i )) 2 (3)
ti (  , )

в данном случае очень трудоемко, так как t i нам неизвестны, и их также необходимо находить численно
методом одномерного поиска.
Для эффективного вычисления функции ошибки предлагается следующий алгоритм. Во-первых,
необходимо все экспериментальные точки { x i , y i } отсортировать так, чтобы с возрастанием номера точки i
соответствующий ей параметр t i также увеличивался. Тогда при поиске параметра t i для каждой
следующей точки можно воспользоваться в качестве начального приближения значением параметра t i 1 , а
для первой точки начальным приближением будет t  0 . Подробнее о составлении контура капли см.
далее.
Во-вторых, вычисление функции ошибки можно проводить непосредственно в процессе интегрирования
системы (1) с помощью метода Рунге-Кутты. В самом деле, на каждой итерации нам доступны значения

____________________________________________________________________________________________
156                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

      
{x(t ), y(t ), x(t ), y(t )} , и наименьшее расстояние от точки может быть найдено с помощью решения
уравнения
d {xi , yi },{x(t ), y (t )}
 ( xi  x(t ))x  ( yi  y (t )) y  0
                       (4)
s
методом Ньютона. То есть, при численном интегрировании системы (1) нужно следить за значением
функции (4) для каждой следующей точки, и запоминать значения наименьших ошибок, при необходимости
уменьшая шаг по   t для увеличения точности результатов.
Выделение контура капли
Как было сказано выше, для эффективного расчета ошибки по формуле (4), необходимо выделить из
изображения контур капли таким образом, чтобы с возрастанием номера точки i соответствующий ей
параметр t i также увеличивался. Данная операция проводится в 2 этапа: непосредственное выделение краев
с помощью детектора Канни (Canny) [5] и выделение из полученной бинарной карты краев связанных
последовательных наборов точек.
Для отслеживания краев был разработан следующий алгоритм. Во-первых, необходимо провести
операцию утончения краев (edge thinning), поскольку детектор Канни не гарантирует, что все полученные
края будут толщиной в 1 пиксель (в основном, такая ситуация возникает в местах соединения), а такое
условие необходимо для дальнейшей обработки. Операция утончения краев может быть проведена с
использованием одного из известных методов утончения краев [6]. В данной работе использовался алгоритм
[7].
Дальнейшая обработка строится на анализе окрестности 3x3 пикселя вокруг рассматриваемого пикселя.
На рис. 3 значения пикселей в окрестности представлены переменными p i , принимающими значение 0 или
1.

Рис. 3. Окрестность 3x3 вокруг рассматриваемого пикселя p 5 , p i  0 или 1 .
Общая схема алгоритма выделения связных последовательностей точек:
1. Найти концы контуров и точки пересечения контуров. Для этого удобно оперировать понятием
«количества переходов уровня» Cx в окрестности 3x3.
Cx  p1  p2  p2  p3  p4  p1  p7  p4  p8  p7  p9  p8  p6  p8  p3  p6 Если Cx  6 и
p5  1 , то в центральном пикселе находится пересечение контуров.
Если   Cx  2 и p5  1 , то в центральном пикселе находится конец контура.
При этом, проверка этих условий может быть быстро и эффективно произведена с помощью таблиц
поиска, так как всего возможных входных значений 512 = 29.
2. Начать с одного из найденных концов контуров.
3. Добавить текущий пиксель в список пикселей контура под текущим номером и пометить на карте краев
текущий пиксель номером текущего контура.
4. Найти среди соседей текущего пикселя пиксель со значением 1.
5. Если найденный сосед не является концом контура или пересечением и не помечен на карте краев еще
ни одним номером, то передвинуть текущий пиксель в позицию найденного соседа и перейти к шагу 3.
Иначе, закончить заполнение текущего контура и перейти к следующему (шаг 2).
Заключение
Экспериментальные исследования системы парафиновое масло / декан в различных концентрациях с
помощью предложенного алгоритма показали эффективность предложенного подхода.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной
России» на 2009–2013 годы.
Литература
1. Maze C., Burnet G. A Non-linear Regression Method for Calculating the Surface Tension and Contact Angle
from the Shape of a Sessile Drop // Surf. Sci. 1969. V. 13. P. 451.
2. Krylov A. S., Vvedensky A. V., Katsnelson A. M., Tugovikov A. E. Software package for determination of sur-
face tension of liquid metals // J. Non-Cryst. Solids. 1993. V. 156-158. P. 845.

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              157
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

3. O. I. del Río and A. W. Neumann. Axisymmetric Drop Shape Analysis: Computational Methods for the Meas-
urement of Interfacial Properties from the Shape and Dimensions of Pendant and Sessile Drops // Journal of Colloid
and Interface Science, Volume 196, Issue 2, 15 December 1997, Pages 136-147.
4. M. Hoorfar and A. W. Neumann. Recent progress in Axisymmetric Drop Shape Analysis // Advances in Col-
loid and Interface Science, Volume 121, Issues 1-3, 13 September 2006, Pages 25-49.
5. Canny, J., A Computational Approach To Edge Detection // IEEE Trans. Pattern Analysis and Machine Intel-
ligence, 8(6):679–698, 1986
6. Lam L., Lee S.-W., Suen C.Y. Thinning Methodologies - A Comprehensive Survey // IEEE Transactions on
Pattern Analysis and Machine Intelligence archive, Volume 14 Issue 9, September 1992.
7. Z. Guo and R. W. Hall, «Parallel thinning with two-subiteration algorithms», Comm. ACM, vol. 32, no. 3, pp.
359-373, 1989.

DROPLET EDGE DETECTION FOR SURFACE TENSION DETERMINATION
Mizotin M.1, Krylov A.1, Protsenko P.2
1
Lomonosov Moscow State University, Faculty of Computational Mathematics and Cybernetics, Laboratory of
Mathematical Methods of Image Processing,
2
Lomonosov Moscow State University, Department of Chemistry

Surface tension is one of the key propertied of liquid, thus its measurement is crucial for studying various phe-
nomena such as wetting and development of technological processes. There sessile and pendant drop techniques are
one of the most frequently used because of their universality and simplicity of measurement process.
The method is based on studying of the axisymmetric drop profile. The balance of gravity force and surface ten-
sion forms the distinct profile shape, thus surface tension can be calculated by the solution of the inverse problem
for the Young-Laplace equation.
In this work the method of droplet contour extraction for determination of the surface tension is presented. The
key difference of the proposed method is its orientation on inexpensive experimental setup using widely available
components such as standard microscope, digital camera and substrate holder. Proposed techniques of image pro-
cessing allow to avoid most of the problems concerning inferior quality of the drop images acquired by inexpensive
setup retaining the measurement accuracy.
The work was supported by federal target program ”Scientific and scientific-pedagogical personnel of innovative
Russia in 2009-2013”.


ПРИМЕНЕНИЕ МЕТОДА МОРФОЛОГИЧЕСКИХ АМЁБ ДЛЯ ВЫДЕЛЕНИЯ
СОСУДОВ НА ИЗОБРАЖЕНИЯХ ГЛАЗНОГО ДНА
Насонов А.В.1, Черноморец А.А. 1, Крылов А.С. 1, Родин А.С.2
МГУ имени М.В. Ломоносова,
1
факультет вычислительной математики и кибернетики, лаборатория математических методов обработки
изображений http://imaging.cs.msu.ru/
2
факультет фундаментальной медицины, кафедра офтальмологии

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

1. Введение
Фотографии глазного дна используются для диагностики заболеваний сетчатки. Сегментация и
оценивание характерных величин сосудов кровеносной системы сетчатки представляют важнейший интерес
при диагностировании и лечении многих заболеваний глаз.
Выделение сосудов на изображениях сетчатки является достаточно сложной задачей обработки
изображений из-за высокого уровня шума, неравномерной освещённости, присутствия объектов, похожих
на сосуды. Среди методов обнаружения сосудов на изображения глазного дна можно выделить следующие
классы [1]:
- класс методов, использующих свёртку изображений с двумерным направленным фильтром и
последующее нахождение пиков откликов. В [2] для сегментации сосудистой сети предложен двумерный
линейный фильтр, профилем которого является гауссиан. Преимуществом данного подхода является
устойчивое нахождение прямолинейных участков сосудов и вычисление их ширины. Однако метод плохо

____________________________________________________________________________________________
158                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

детектирует тонкие и извилистые сосуды, возможны ложные срабатывания на объекты, не являющимися
сосудами, например на экссудаты.
- методы, использующие детектирование хребтов. В [3] производится нахождение примитивов —
коротких отрезков, лежащих посередине линий, затем с помощью методов машинного обучения отбираются
примитивы, соответствующие сосудам, по которым восстанавливается сосудистое дерево.
- методы, использующие трекинг сосудов, включающий в себя как соединение сосудов по паре точек, так
и продолжение сосудов [4]. К преимуществам данного подхода можно отнести высокую точность работы на
тонких сосудах и восстановление разрывных сосудов. Недостатком является сложность обработки
ветвлений и пересечений сосудов.
- попиксельная классификация, основанная на применении методов машинного обучения [5]. Здесь для
каждого пикселя строится вектор признаков, на основе которого определяется, является ли пиксель частью
сосуда или нет. Для обучения метода используются изображения глазного дна с размеченными на нём
экспертом сосудами. К недостаткам метода можно отнести большое расхождение в мнениях экспертов.
В данной работе для выделения сосудов используется метод морфологических амёб — морфологический
метод, при котором структурный элемент выбирается адаптивно для каждого пикселя.
2. Морфологические амёбы
Мы используем метод морфологических амёб, описанный в [6], с модифицированной функцией
расстояния.
Рассмотрим изображение в градациях серого I ( x, y) . Представим его в виде графа, в котором каждый
пиксель соединён с восемью соседними пикселям рёбрами с некоторым заданными весами («стоимостью»).
Тогда для каждого пикселя ( x0 , y0 ) можно найти множество всех точек ( x, y )  A(( x0 , y0 ), t ) , для
которых стоимость пути из ( x0 , y0 ) в ( x, y ) не превышает t. Полученное множество и будет являться
структурным элементом для пикселя ( x0 , y0 ) .
Мы используем следующую функцию расстояния между пикселями ( x0 , y0 ) и ( x1 , y1 ) :
d ((x0 , y0 ), ( x1 , y1 )) 
( x1  x0 ) 2  ( y1  y0 ) 2 ( I ( x0 , y0 ))2  ( I ( x1 , y1 ))2   | I ( x1 , y1 )  I ( x0 , y0 ) |,   0.
Сомножитель                   ( I ( x0 , y0 ) 2  I ( x1 , y1 ) 2 задаёт низкую стоимость перемещения по тёмным участкам и
высокую — по светлым, тем самым не давая амёбе распространяться по точкам вне сосуда, а слагаемое
 | I ( x1 , y1 )  I ( x0 , y0 ) | штрафует перемещение между пикселями с сильно различающейся интенсивностью.
Параметр  задаёт значимость штрафа при данном переходе.
Пример нахождения амёб при t  125,   0.2 приведён на рис. 1.

Рис. 1. Примеры форм морфологических амеб. Слева — исходное изображение с помеченными точками,
в которых вычисляются амёбы, справа — белым помечены найдённые структурные элементы.
3. Выделение сосудов с помощью морфологических амёб
Для прослеживания сосудов кровеносной системы на изображениях глазного дна был разработан
алгоритм, состоящий из следующих этапов:
1. Выделение зелёного канала как наиболее информативного и выравнивание освещенности на этом
канале с помощью метода [7]. Это позволяет использовать одни и те же параметры λ и t при построении
морфологических амеб для разных изображений.
2. Нахождение на полученном изображении множества точек {pn}, заведомо принадлежащих сосудам.
Мы использовали следующий метод: производилось детектирование границ сосудов с помощью метода

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              159
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

Канни, затем применялась скелетизация множества точек, не являющихся границами; связные фрагменты
скелета считались частью сосудов, если интенсивность изображения в этих точках была ниже некоторого
порогового значения.
3. Для каждой начальной точки pi  pn  строится амеба А(pi). К маске амёбы применяется ранговая
фильтрация с окном 3х3 для устранения неровностей, вызванных шумом: если у точки маски менее 3
соседей, принадлежащих маске, то точка выкидывается. Оставшиеся точки маски помечаются как точки
сосудов.
4. В случае, когда требуется не только заполнение сосудов, но и их продолжение, то шаг 3 повторяется
для всех точек, помеченных как точки сосудов, для которых ещё не была вычислена амёба.
4. Результаты
Пример работы алгоритма приведён на рис. 2.

Рис. 2. Результат выделения сосудов при помощи морфологических амеб. Слева — изображение
глазного дна (зелёный канал), по центру — точки, заведомо являющиеся точками сосудов, от которых будут
строиться амёбы, справа — результат выделения сосудов с помощью предложенного метода.
Заключение
Рассмотрено применение метода морфологических амёб для выделения сосудов на изображениях
глазного дна.
Разработанный алгоритм планируется использовать в автоматизированной системе обнаружения
заболеваний сетчатки.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной
России» на 2009–2013 годы и гранта РФФИ 10-01-00535-а.
Литература
1. R.J.Winder, P.J.Morrow, I.N.McRitchie, J.R.Bailie, P.M.Hart. Algorithms for digital image processing in diabetic retinopa-
thy // Computerized Medical Imaging and Graphics, Vol. 33, 2009, 608–622.
2. S.Chaudhuri, S.Chatterjee, N.Katz, M.Nelson, M.Goldbaum. Detection of Blood Vessels in Retinal Images Using Two-
Dimensional Matched Filters // IEEE Transactions of Medical Imaging, Vol. 8, No. 3, 1989, pp. 263–269.
3. J.Staal, M.D.Abramoff, M. Niemeijer, M.A.Viergever, B.Ginneken. Ridge-Based Vessel Segmentation in Color Images of
the Retina // IEEE Transactions on Medical Imaging, Vol. 23, No. 4, 2004, pp. 504–509.
4. M.Patasius, V.Marozas, D.Jegelevieius, A.Lukosevieius. Recursive Algorithm for Blood Vessel Detection in Eye Fundus
Images: Preliminary Results // IFMBE Proceedings, Vol. 25/11, 2009, pp. 212–215.
5. J.Soares, J.Leandro, R.Cesar Jr., H.Jelinek, M.Cree. Retinal Vessel Segmentation Using the 2-D Gabor Wavelet and Super-
vised Classification // IEEE Transactions of Medical Imaging, Vol. 25, No. 9, 2006, pp. 1214–1222.
6. M.Welk, M.Breub, O.Vogel. Differential Equations for Morphological Amoebas // Lecture Notes in Computer Science,
Vol. 5720/2009, 2009, pp. 104–114.
7. G.D.Joshi, J.Sivaswamy. Colour Retinal Image Enhancement based on Domain Knowledge // Sixth Indian Conference on
Computer Vision, Graphics and Image Processing (ICVGIP'08), 2008, pp. 591–598.

APPLICATION OF MORPHOLOGICAL AMOEAS METHOD FOR BLOOD VESSEL DETECTION IN
EYE FUNDUS IMAGES
Nasonov A.1, Chernomorets A.1, Krylov A.1, Rodin A.2
Lomonosov Moscow State University,
1
Faculty of Computational Mathematics and Cybernetics, Laboratory of Mathematical Methods of Image Pro-
cessing, http://imaging.cs.msu.ru/
2
Faculty of Fundamental Medicine, Department of Ophthalmology

An algorithm of blood vessels detection in eye fundus images has been developed. Segmentation and analysis of
blood vessels in eye fundus images provides the most important information to diagnose retinal diseases.
Blood vessels detection in eye fundus images is a challenging problem [1]. Images are corrupted by non-uniform
illumination and noise. Also some objects can be wrongly detected as blood vessels.

____________________________________________________________________________________________
160                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

The proposed algorithm is based on the method of morphological amoebas [2]. Morphological amoeba for a giv-
en pixel is a set of pixels with the minimal distance to the given pixel less than a threshold t. We use the sum of av-
erage intensity value multiplied by Euclidean distance and absolute value of difference between pixel intensity val-
ues for the distance. In this case the distance will be small for blood vessels which are usually dark and big for light
areas and edges, and the amoeba will be extended along the vessel but not through vessel walls.
The proposed algorithm of blood vessel detection consists of the following steps:
Extract the green channel as the most informative and perform illumination correction using the method [3]. It
makes possible to use unified amoebas parameters for different images.
Find the set of pixels {pn} in the obtained image which are surely the pixels of the blood vessels
Calculate the amoeba А(pi) for every pixel pi  pn  , apply rank filtering to the amoeba mask with 3x3 win-
dow: remove the pixels from the mask which have less than 3 neighbor pixels in the mask. The remaining pixels are
marked as blood vessels pixels.
If we need to extend the blood vessels, the third step is repeated for all newly added pixels to blood vessels area.
We plan to use the developed algorithm in automatic system of retinal disease detection.
The work was supported by federal target program ”Scientific and scientific-pedagogical personnel of innovative
Russia in 2009-2013” and RFBR grant 10-01-00535-а.
Literature
1. R.J.Winder, P.J.Morrow, I.N.McRitchie, J.R.Bailie, P.M.Hart. Algorithms for digital image processing in dia-
betic retinopathy // Computerized Medical Imaging and Graphics, Vol. 33, 2009, 608–622.
2. M.Welk, M.Breub, O.Vogel. Differential Equations for Morphological Amoebas // Lecture Notes in Computer
Science, Vol. 5720/2009, 2009, pp. 104–114.
3. G.D.Joshi, J.Sivaswamy. Colour Retinal Image Enhancement based on Domain Knowledge // Sixth Indian
Conference on Computer Vision, Graphics and Image Processing (ICVGIP'08), 2008, pp. 591–598.


ПРИМЕНЕНИЕ ТЕХНОЛОГИИ CUDA ДЛЯ ФРАКТАЛЬНОГО СЖАТИЯ ИЗОБРАЖЕНИЯ
Новинский Н.Б.
Федеральное государственное унитарное предприятие «Главный радиочастотный центр» (ФГУП «ГРЧЦ»)

Для получения аттракторов, близких к исходным изображениям, методом полного перебора, необходим
огромный вычислительный ресурс. Поскольку вычислительные способности, существующих на данный
момент вычислительных систем, например, на базе центральных процессоров от Intel явно недостаточно,
было принято решение попробовать многопроцессорные решения на базе технологии CUDA от nVidia.
Введение.
Фрактальное сжатие основано на получении некоторого преобразования, применение которого
несколько раз к произвольному изображению позволяет получить новое изображение, не меняющееся при
очередном применении этого преобразования. Такое преобразование называется фрактальным, а
полученное таким образом изображение называется аттрактором. Впервые такие преобразования описаны
Мандельбротом в [1]. В [2] Барнсли подробно описал один из возможных алгоритмов получения
преобразования. Он заключается в разбиении исходного изображения на множество квадратных ранговых
областей и подбора для каждой похожего домена. Домены при этом уменьшаются в два раза (до размера
ранговой области), поворачиваются при этом в восемь разных положений, приводятся по яркости и
контрасту, и сравниваются в среднеквадратичном смысле с ранговой областью. Следует заметить, что такая
методика поиска чрезвычайно затратна по вычислительным ресурсам, поскольку требует многократного
полного перебора всех возможных доменов. Причем все существующие методики ускоренных расчетов, как
правило, дают худший результат.
В проводившихся ранее исследованиях на тему фрактального сжатия, использовались изображения
небольших размеров (в основном 256х256, 512х512). Не в последнюю очередь это было связано с
ограниченными вычислительными ресурсами. Однако, в настоящее время, в связи с широким внедрением в
массовое пользование фотографической техники с большими матрицами и повсеместным внедрением новых
стандартов представления видеоинформации с большим разрешением, обработка больших изображений
приобретает особенную актуальность. С другой стороны чем больше изображение, тем выше возможность
подобрать действительно похожую пару домен – ранговая область, что должно благотворно сказываться на
качестве полученного аттрактора.
Работать предполагалось с изображениями нескольких размеров. В таблице 1 в третьей колонке
представлены размеры наиболее известных видеоформатов. В первой колонке приводится сленговое
название (оно коррелирует с числом строк), к которому будут отсылки в тексте данной статьи. В последнем
столбце размеры изображений на которых планировалось производить расчеты. Размеры выбраны исходя
____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              161
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

из кратности оных размерам ранговых областей. Формат 0,5К хоть и не совпадает ни с одним
видеоформатом, однако близок по размерам к изображениям, использовавшимся для экспериментов с
изображениями ранее, и именно на них и проводилась тестовая отработка алгоритмов.

Таб.1 Наиболее известные форматы видеоизображений.
Следует заметить, что в данной работе вообще говорить о сжатии не приходится. Поскольку результатом
расчетов является текстовый файл (для удобства анализа), из которого можно было получить аттрактор
максимально близкий к исходному изображению.
Размеры ранговых областей
4х4          8х8            16х16
CPU                                                              83 902       62 353         54 996
GPU, блоки в одну нить, использование константной памяти         11 449       8 935          7 964
GPU, блоки со множественными нитями, использование               2 518        678            393
разделяемой памяти
GPU, оптимизация вывода                                          428           180           393
Таб.2. Время получения аттрактора в сек. для изображения 0.5К.
Расчеты на CPU.
Алгоритм подбора домена для конкретной ранговой области подробно описан в [2]. Это алгоритм - один
из наиболее широко представленных в литературе и исследованиях. Когда он был реализован “в лоб” на
одном ядре на процессоре Intel Core Quard CPU Q9300 2,5 ГГц под ОС Linux , стало ясно, что это даже для
современного процессора слишком большая задача. Были произведены работы по ускорению расчетов. Во-
первых, вместо получения каждого отдельного домена, усреднением четырех соседних пикселов в один и
поворотом домена в нужное положение, были рассчитаны 32 плоскости, из которых брались исходные
значения для сравнения с ранговой областью. Это позволило сократить время расчетов примерно на 30
процентов. В табл. 2, в первой строке представлены времена расчетов для этого варианта алгоритма для
изображения 0.5К для ранговых областей разного размера. Для справки, в сутках 86400 сек. Остальные
попытки ускорить расчеты, например предварительным расчетом промежуточных сумм особого успеха не
принесли, поскольку давали не более10-15 процентов уменьшения времени расчетов, но массивы, которые
получались, категорически не влезали в оперативную память. Конечно, возможно было разделить задачу на
все 4 ядра, и получить ускорение почти в 4 раза, однако, это тоже проблемы не решало, поскольку при
росте линейного размера изображения в два раза время расчетов вырастает в 16 раз и изображение 4К
должно обсчитываться примерно 1000 суток даже на 4-х ядрах.
На этом этапе было принято решение попытаться реализовать алгоритм на графическом сопроцессоре.
Собственно реализации данного алгоритма на графическом процессоре посвящена данная статья. В
качестве тестового варианта была выбрана видеокарта от nVidia на чипсете GeForce 250, как недорогое и
достаточно производительное решение.
Расчеты c применением технологии CUDA.
Несколько слов о вычислениях на технологии CUDA.
Процессоры объединены в мультипроцессоры по 8 штук. Одна вычислительная задача в терминах
CUDA называется нить. Они могут выполняться параллельно. Нити объединяются в блоки. Число нитей в
блоке может быть, в достаточной степени, произвольным.
Плата обладает набором памяти различного назначения и свойств и грамотное использовании всего
этого многообразия несомненно обеспечит успех. Память бывает глобальной, константной, разделяемой.
Глобальная память большая (несколько гигабайт), не кешируемая и чтение одного числа занимает 400-600
тактов процессора. Что безусловно тормозит процесс расчетов. Зато все нити всех блоков имеют к ней
доступ. Ускорить чтение можно следующим образом. Если каждая нить читает последовательно числа из
глобальной памяти, причем все они принадлежат одному блоку, этих нитей 16 и чтение осуществляется с
начала банка (по 16 чисел), то тогда скорость чтения увеличивается в 16 раз.
Константная память тоже как и глобальная обеспечивает доступ с любой нити любого блока, однако она
кешируемая и маленькая. Там возможно держать небольшие объемы данных (десятки килобайт), которые
нужны во всей программе. Чтение оттуда происходит достаточно быстро.
Разделяемая память маленькая (десятки килобайт), и доступ к ней возможен только из нитей одного
блока. Зато очень быстрая.
В общем случае программа на CUDA состоит из пяти частей.
1. Инициализация программы на CPU.
2. Загрузка массивов данных на видеоплату за счет CPU.
____________________________________________________________________________________________
162                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

3. Выполнение расчетов на GPU (ядро в терминах CUDA).
4. Выгрузка массивов данных в основную память за счет CPU.
5. Окончательная обработка данных и вывод результатов на CPU.
В данном конкретном случае пункты 2-4 выполняются неоднократно.
В начале задача была написана следующим образом. Бралась очередная ранговая область, и помещалась
в константную память (пункт 2). В глобальную память подгружалась плоскость(пункт 2). В каждом блоке
была только одна нить, которая считала соответствие ранговой области к очередному домену(пункт 3).
Результатом расчета по каждой паре домен – ранговая область являлись 3 числа типа float, это
коэффициенты преобразования для яркости и контраста и среднеквадратичная разница домен – ранговая
область. Они выгружались в основную память (пункт 4). Затем из всех этих значений среднеквадратичной
разницы для данной ранговой области выбирается минимальное (пункт 5, частично). Потом перебираются
все ранговые области, а потом перебираются все плоскости. Результаты этих расчетов представлены в табл.
2 на второй строке. Как видно это уже замечательный результат. Ускорение почти на порядок.
На следующем этапе была введена множественность нитей в рамках одного блока, а количество блоков
соответственно уменьшилось. Количество нитей в блоке было выбрано как 16 штук по горизонтали и одна
высота ранговой области по вертикали. Каждая нить из блока копирует из глобальной памяти в разделяемую
четыре значения в результате получается матрица, которая копирует участок плоскости размером 32 точки
по горизонтали и два размера ранговой области по вертикали. Расчет производится уже на основе матрицы,
которая лежит в разделяемой памяти. Результат представлен в строке 3 в табл.2. Результат расчетов для
ранговых областей 16х16 получился еще луче. Время расчета для ранговых областей 4х4 оставалось все
таки слишком велико.
Были проведены замеры затрат времени на различные операции. Выяснилось, что на выполнение самого
ядра уходит примерно пятая часть времени. Еще три пятых времени уходит на копирование массивов с
видеоплаты в основную память. И еще пятая часть времени уходит на перебор массива в основной памяти
на CPU, с целью обнаружения минимального значения для среднеквадратичной разницы. Остальные
операции во всей программе занимали по времени не более одного процента.
По результатам этого исследования была произведена модернизация программы. Сортировка
результатов внутри блока была вынесена на GPU. За счет этого были уменьшены выходные массивы в
несколько десятков раз (от 64 до 256 раз в зависимости от размеров ранговой области). Это благотворно
сказалось на результатах (нижняя строчка в табл.2).
Выводы.
1. Задача фрактального сжатия, благодатна в плане обработки на GPU поскольку хорошо
распараллеливается и имеет значительные вычислительные затраты при относительно небольших входных и
еще меньших выходных массивах.
2. Применение GPU позволяет достаточно резко (200-300 раз!) ускорить решение подобного рода задач.
Расчетное время для обработки изображения 4К должно составить порядка 20 суток, что вполне приемлемо.
3. На данный момент существуют видеокарты от nVidia на новом 580-ом чипсете. Ожидается, что их
применение и может дать выигрыш минимум в 2-3 раза. Помимо этого есть профессиональные решения под
маркой TESLA, которые должны обладать еще большей производительностью. Возможно также
дополнительное распараллеливание за счет применения нескольких видеокарт в рамках одной
вычислительной системы.
Литература
1.    Mandelbrot B.B. The Fractal Geometry of Nature / W.H.Freeman, NY, 1983
2.    Barnsley M.F. Fractal Image Compression, SIGRAPH’92 Course Notes, 1992
3.    Боресков А.В., Харламов А.А. Основы работы с технологией CUDA.- M.: ДМК Пресс, 2010.
4.    CUDA Reference Manual 2.3. July 2009

USING OF CUDA TECHNOLOGY FOR FRACTAL IMAGE COMPRESSION
Novinsky N.
The biggest problem of fractal compression is tremendous time for encoding of image. In that work we used
CUDA technology for brute force attractor searching.
The device with CUDA is ideally suited for computations that can be run in parallel.
CUDA devices use several memory spaces, which have different characteristics that reflect their distinct usages
in CUDA applications. These memory spaces include global, local, shared and other.
Global memory is big but very slow. But local and shared memory is small but very fast.
Perhaps the single most important performance consideration in programming for the CUDA architecture is coa-
lescing global memory accesses. Global memory loads and stores by threads of a half warp (16 threads) are coa-
lesced by the device in as few as one transaction (or two transactions in the case of 128-bit words) when certain ac-
cess requirements are met.

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              163
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

In my investigation I worked with CUDA device in true manner, and in result was achieved acceleration in 200-
300 times.

ЦИФРОВАЯ ОБРАБОТКА ИЗОБРАЖЕНИЙ ДЛЯ ДИАГНОСТИКИ ВЫСОКОТЕМПЕРАТУРНЫХ
ПРОЦЕССОВ
Орловская С.Г., Каримова Ф.Ф., Калинчак В.В., Шкоропадо М.С., Зуй О.Н.
Одесский национальный университет им. И.И. Мечникова, Одесса, Украина

Высокотемпературные реагирующие системы, такие как пламена, характеризуются высокими
градиентами температуры и концентраций компонентов. Так как потоки тепла и компонентов зависят от
производных, то исследования процессов горения и высокотемпературного окисления требуют методов
диагностики с высоким пространственным и временным разрешением. Современные цифровые камеры
обеспечивают регистрацию процессов с высоким пространственно-временным разрешением. В настоящее
время активно разрабатываются новые методы диагностики, базирующиеся на регистрации излучения
объектов цифровой камерой и последующей компьютерной обработке изображений.
Для изучения высокотемпературного окисления тугоплавких металлов (вольфрам, молибден)
используется метод «горячей нити»: металлические проволочки нагреваются постоянным электрическим
током в контролируемой атмосфере. При этом, температуру проводника обычно измеряют либо оптическим
пирометром, либо определяют по изменению его сопротивления (электротермографический метод). Однако,
в первом случае определяется локальная яркостная температура, а во втором усредненное по длине
проводника значение температуры. Для получения более полной информации о температурном поле
проводника был разработан метод относительной яркостной пирометрии, который позволяет определять
распределение температуры по поверхности нагретого тела [1]. Метод основан на прецизионном измерении
температуры Т 0 в центре проводника x 0 эталонным пирометром с исчезающей нитью. Одновременно
излучение проводника регистрируется цифровой камерой через интерференционный фильтр в формате .tiff
либо .raw. Время экспозиции подбирается таким образом, чтобы сигнал матрицы находился в области
линейности.
Полученное изображение обрабатывается на компьютере в пакете Matlab 7.0. Соотношение между
температурой T  x  в точке x объекта и сигналом матрицы S  x  в соответствующей точке изображения
находится из соотношения, полученного из формулы Вина:         1    1     S0  ,            (1)
     ln       
T x  T0 C 2  S х  
      
где C 2  1.4388  10 2 м  К ;  – длина волны, мкм; S 0                                       Т0 .
– сигнал матрицы в точке, соответствующий
Данный метод использовался для определения температурного профиля вдоль металлических
проволочек (вольфрам, молибден, титан), нагреваемых электрическим током в воздухе [2]. На рисунке 1
приведены типичные               профили температур. Совместное использование описанного метода и
электротермографического, позволяет оценить спектральный коэффициент излучательной способности
 ( , T ) и определить значения действительной температуры в центральной (светящейся) зоне проводника.

Рис. 1. Распределение температуры по длине вольфрамового проводника диаметром d = 300 мкм при
токах нагрева: 1 – 6.8 А; 2 – 7.3 А.
Для диагностики процессов горения газообразного топлива на применен модифицированный метод
пирометрии тонкой нити (thin-filament pyrometry) [3]. Принципиальная блок-схема установки представлена
на рис.2. Для исследования температурного поля в факеле предварительно перемешанной пропан-
бутановой смеси на бунзеновской горелке (1), использовали кварцевые нити (диаметр d =150 мкм, длина l =
0,1 м), расположенных горизонтально с шагом 1 см (2), излучение от которых регистрировалось в узком

____________________________________________________________________________________________
164                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

спектральном интервале через светофильтр (   0,77 мкм) (3) с помощью цифровой камеры canon 350d (4).
Полученные изображения обрабатывались на компьютере (5).

Рис. 2. Блок- схема экспериментального стенда.
В отличие от описанных в литературе модификаций метода TFP, разработанный нами метод [3, 4]
отличается тем, что в качестве точки отсчета используется значение яркостной температуры, определенное с
помощью эталонного яркостного пирометра. В качестве излучателя мы использовали тонкие кварцевые
нити. Распределение температуры по длине нити Т(х) рассчитывали по формуле (1).
Важной характеристикой горения газообразного топлива является нормальная скорость распространения
пламени. Нормальная скорость распространения пламени рассчитывают по следующей формуле
U n  Q / S ,  где Q - расход газовоздушной смеси, см3/сек, а S - площадь поверхности факела, см2.
Как правило, фронт пламени аппроксимируется конической поверхностью и его площадь S
рассчитывается по формуле для боковой поверхности конуса. Для этого достаточно измерить только
высоту факела (диаметр среза горелки известен). Данный способ характеризуется существенной
методической ошибкой, величина которой возрастает при увеличении скорости горения и соответствующем
уменьшении высоты факела. Кроме того, заметной может быть погрешность определения положения
вершины факела, особенно при наличии пульсаций. Поэтому была разработана методика определения
границ факела (фронта горения) и точного расчета площади его поверхности.
Установлено, что изображение факела в зеленом канале (излучение радикала С 2) имеет четко
выраженные границы и максимальное отношение сигнал/шум. Поэтому для определения положения фронта
пламени использовали информацию зеленого канала. Компьютерная обработка изображения включала
медианную фильтрацию для уменьшения случайных шумов, бинаризацию с регулируемым порогом яркости
(от 0.25 до 0.35) и выделение контура бинарного изображения факела. Полученный контур
аппроксимировали полиномом шестой степени и далее вычисляли поверхность факела по формуле для
поверхности вращения. С помощью предложенной методики выполнены измерения нормальной скорости
горения предварительно перемешанной смеси пропан-бутан-воздух для разных значений коэффициента
избытка окислителя [4]. Метод обеспечивает хорошую точность определения нормальной скорости горения
и может быть использован для нестационарных пламен.
Для изучения процессов воспламенения и горения частиц твердых натуральных топлив нами разработан
комплексный метод, включающий измерения с помощью термопары и метод обработки последовательных
цифровых изображений объекта исследований. Частица угля вносилась на термопаре в нагретую до высокой
температуры печь. Сигнал с термопары поступал на аналогово-цифровой преобразователь, а затем на
персональный компьютер. Одновременно при помощи цифровой камеры производилась киносъемка
окисляющейся частицы, а так же измерение ее температуры яркостным пирометром. Так как термопары
обладают большой инерционностью, то на быстропротекающих стадиях высокотемпературного
тепломассообмена основным измерительным блоком был цифровой метод регистрации процесса с
последующей компьютерной обработкой изображений.
d, мкм
5600

4800

4000

3200

2400                                                 1
1600

800
t, c
0
0   100   200   300   400   500   600   700   800   900 1000

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              165
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

Рис. 3. Временная зависимость диаметра угольной частицы с начальным диаметром d b = 5,6 мм при
температуре газа T g = 1145 К. 1 - расчет, ● - эксперимент.

Цифровая съемка частицы в процессе высокотемпературного окисления позволила проследить динамику
изменения ее диаметра (рис. 3) и зафиксировать момент потухания. Экспериментальная зависимость d t 
хорошо согласуется с расчетами по физико-математической модели (рис.3, кривая 1), включающей
уравнения теплового баланса и химической кинетики для угольной частицы [5].
Измерение температуры частиц позволило установить стадийность процессов тепломассообмена и
химического превращения [6]: инертный прогрев частицы до температуры газа, активация химических
реакций на ее поверхности, стадия высокотемпературного окисления, потухание. На стадии
высокотемпературного окисления температура частицы равна 1250÷1300 К и медленно возрастает. Затем по
мере уменьшения диаметра окисляющейся частицы снижаются и теплопотери излучением с ее поверхности.
При помощи цифровой съемки, в конце этой стадии, обнаружено резкое увеличение интенсивности
свечения частицы, а следовательно, и ее температуры, после чего частица затухает. Этот факт
подтверждается расчетами по физико-математической модели [5] и в экспериментальном плане изучен
недостаточно [6].
Таким образом, разработанные нами методы цифровой регистрации и обработки изображений
высокотемпературных состояний излучающих объектов исследования (газообразного, металлического и
угольного топлива) позволила раскрыть закономерности и механизмы их окисления и горения.
Литература
1. Патент України № 44416. Карімова Ф.Ф., Орловська С.Г.. Спосіб визначення локальної яскравісної
температури в окремих точках нагрітого тіла та розподілу яскравісної температури по поверхні нагрітого
тіла. 12.10.2009. Бюл. №19. 2009.− 5 с.
2. S. Orlovskaya, F.Karimova, M. Shkoropado. Heat transfer during high temperature oxidation of molybdenum
filament //The seventh World Conference on Experimental Heat Transfer, Fluid Mechanics and Thermodynamics
(ExHFT-7), Krakow, Poland. - 28 June-3 July 2009. - Р.2129–2133.
3. Каримова Ф.Ф., Орловская С.Г., Шкоропадо М.С., Протас С.К., Панов В.В. Исследование температурного
поля газовоздушного пламени. Метод пирометрии тонкой нити//Физика аэродисперсных систем. – 2009,
Вып. 46, Одесса. – С. 95-100.
4. Каримова Ф.Ф., Смагленко Т.Ф. Влияние коэффициента избытка окислителя на скорость горения и спектр
излучения пламени пропан-бутан-воздух // Физика аэродисперсных систем. – 2004, Вып. 41, Одесса. – С.
165–202.
5. Орловская С.Г., Рябчук Л.И., Мирошниченко Е.В., Кысса В.Д. Влияние температурного поля по
углеродной частице на характеристики её тепломассообмена с газом // Труды пятой Российской
национальной конференции по теплообмену. В 8 томах (25—29 октября 2010 г., Москва). Т. 3. Свободная
конвекция. Тепломассообмен при химических превращениях. — М.: Издательский дом МЭИ, 2010. – С.280-
282.
6. Калинчак В.В., Орловская С.Г., Рябчук Л.И., Зуй О.Н. Влияние неоднородности поля температур внутри
углеродных частиц на их воспламенение и горение // Современная наука: исследования, идеи, результаты,
технологии. Сборник научных статей. – Вып. 2(4). – Киев: «НПВК Триакон», 2010. – С.98-103.

HIGH TEMPERATURE FAST PROCESSES DIAGNOSTICS BY DIGITAL IMAGING PROCESSING.

Orlovskaya S., Karimova F., Kalinchak V., Shkoropado M., Zuj O.

Odessa I.I. Mechnikov’s National University

To study high temperature oxidation and combustion processes new diagnostics methods are elaborated based on
digital images computer processing.
Relative brightness pyrometry method is used to define the temperature distribution along hot body surface, for
example thin filament [1]. The method is based on precise measurement of the reference temperature Т 0 in the
central point of the filament by etalon disappeared filament pyrometer. Simultaneously, the filament radiation is
registered by digital camera (tiff or raw format) through interferential filter. The exposure time is selected to provide
the linear matrix response. The image obtained is processed by use Image Processing Toolbox (MatLab 7.0). The
temperature T  x  at the point x of the filament is connected with the matrix response S  x  at the corresponding
image element by the next equation following from Wien’s formula:                          1    1    S0  , here
    ln       
T x  T0 C2  S х  
      

____________________________________________________________________________________________
166                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

 – the interferential filter wavelength, mcm; S 0 – the CCD matrix response in the
C 2  hc k  1.4388 10 2 m  К ;
reference point x 0 with temperature Т 0 .
This method was used to define the temperature distribution along refractory metals filaments (W, Mo, Ti), heat-
ed electrically in air [2].
The gas flame diagnostics method is elaborated. The temperature distribution in the propane/butane/air premixed
flame is defined by “thin filament pyrometry” method with use the formula (1) [3]. The set of parallel quartz fila-
ments is placed in the flame and the glowing filaments were recorded. Then the radial temperature profiles in the
flame cross sections are obtained by the images processing. The reference point x 0 is selected in the intersection of
one of filaments with the flame front. The temperature T0 ( x0 ) is defined by disappearing filament pyrometer. To
define the normal flame velocity the surface of the flame front is calculated with aid of image processing original
routine [4].
Ignition, burning and extinction of a single carbon particle in hot air are studied. The burning particle is regis-
tered by digital camera and its temperature is continuously measured by thermocouple. A comparison of the particle
temperature history with its diameter change allows specifying the particle combustion mechanism.

[1] Patent of Ukraine № 44416. Karimova F.F., Orlovskaya S.G.. The method of measuring the brightness tempera-
ture local values and brightness temperature distribution along hot body surface. 12.10.2009. Bull. №19. 2009 − 5
p.
[2] S. Orlovskaya, F. Karimova, M. Shkoropado. HEAT TRANSFER DURING HIGH TEMPERATURE
OXIDATION OF MOLYBDENUM FILAMENT//The seventh World Conference on Experimental Heat Transfer,
Fluid Mechanics and Thermodynamics (ExHFT-7), Krakow, Poland. - 28 June-3 July 2009. - Р.2129–2133.
[3] F. Karimova, S. Orlovskaya, M. Shkoropado, S. Protas, V. Panov. The gas flame temperature field investigation.
Thin filament pyrometry // Fizika aerodipersnyh system. V.46, Odessa.-2009. P. 95-100.
[4] F. Karimova, T. Smaglenko. Air-fuel ratio effect on burning rate and emission spectrum of propane/butane/air
flame // Fizika aerodipersnyh system. V.41, Odessa.-2004. P. 165-202.


МЕТОД ПРОЕКЦИОННОЙ ФАЗОВОЙ КОРРЕЛЯЦИИ В ЗАДАЧЕ ИДЕНТИФИКАЦИИ ЧЕЛОВЕКА
ПО РАДУЖНОЙ ОБОЛОЧКЕ ГЛАЗА
Павельева Е.А., Крылов А.С.
МГУ имени М.В. Ломоносова,
Факультет Вычислительной Математики и Кибернетики,
Лаборатория математических методов обработки изображений, http://imaging.cs.msu.ru
Предложен метод идентификации человека по радужной оболочке глаза на основе оценки фазовой
корреляции изображений. Фаза изображения считается с помощью проекционного метода Эрмита разложения
функции интенсивности изображения в ряд по функциям Эрмита - собственным функциям преобразования
Фурье.

1.Введение
Задача идентификации человека по радужной оболочке глаза является классической задачей биометрии.
Существует много методов биометрической идентификации. Большинство методов используют локальную
текстурную информацию изображений. Обычно в таких методах исследуются знаки сверток интенсивности
областей изображения с выбранным набором функций. Так же существуют методы, основывающиеся на
информации об изображениях в целом. Известно, что фаза изображения несет в себе гораздо больше
информации, чем его спектр [1]. Одним из методов, опирающихся на этот факт, является метод фазовой
корреляции изображений, использующий только информацию о фазе, называемый POC (Phase-only correla-
tion) [2]. В фазовом методе для пары изображений исследуется поведение разности фаз этих изображений и
показано, что этот параметр характеризует меру сходства изображений.
В данной работе предложен метод проекционной фазовой корреляции изображений, в котором
преобразование Фурье изображения считается с использованием разложения функции интенсивности
изображения в ряд по функциям Эрмита (собственным функциям преобразования Фурье).
2. Метод фазовой корреляции
Рассмотрим две функции                f (x )       и   g (x )      на отрезке        0, M .   Пусть   F (u )   и   G (u )   - их
i 2ux
1      M                            1   M
преобразования Фурье:         F (u )              f ( x) e          M             f ( x) (WM )ux  AF (u)eiF (u ) ,
M      x0                           M   x0

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              167
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

i 2

где WM  e          M     ,   AF -модуль преобразования Фурье, а  F                         - его фаза.
Взаимным фазовым спектром (cross-phase spectrum) двух спектральных функций                                               F (u )      и     G (u )
называется функция R (u )  F (u ) G (u )  e i (F (u )G (u )) , т.е. спектральная функция с единичным
FG
F (u ) G (u )
модулем, фаза которой равна разности фаз функций. Взяв от этой функции обратное преобразование Фурье,
1    M

получим POC-функцию фазовой корреляции:                                         POC fg ( x)            RFG (u)WMux .
M   u 0
i 2ua

Если g ( x)  f ( x  a ) , т.е. функции одинаковы с точностью до смещения, то G (u )  e M F (u ) и
POC-функция является дельта-функцией с пиком в точке x  a . Метод фазовой корреляции состоит в
следующем: если две функции                             f (x )   и   g (x )      “похожи”, то POC-функция дает четкий пик, если же
функции “не похожи”, то POC-функция не дает четкого пика. Высота пика определяет меру похожести
функций, а положение пика соответствует смещению одной функции относительно другой.
Метод фазовой корреляции применяется в различных задачах обработки изображений и биометрической
идентификации, в частности, в задаче идентификации по радужной оболочке глаза [3]. При этом показано,
что если брать обратное преобразование Фурье не от всего спектрального сигнала, а только от его части,
соответствующей низким частотам, то пик POC-функции получается более четким и устойчивым.
3. Метод проекционной фазовой корреляции
В работе предложен метод фазовой корреляции основанный на проекционном методе Эрмита [4]. Он
основан на разложении исходных функций в ряд по функциям Эрмита                                               n    – собственным функциям
преобразования Фурье, образующим полную ортонормированную систему в                                           L2 (  ,  ) :
(1) n  e  x
2
/2                          H 0 ( x)  1, H1 ( x)  2  x,
 n ( x)                                   H n ( x) , где
2 n  n!                                        H n ( x)  2  x  H n1 ( x)  2  (n  1)  H n2 ( x).
Каждая функция Эрмита                           n     является локализованной с вычислительной точки зрения на отрезке
[ 2n  1, 2n  1] .
Если функция              f (x )           задана на отрезке                A, A , то, считая, что      функции Эрмита растянуты в
A     раз (чтобы функция Эрмита с максимальным номером                                            n–   выбранным числом функций
k
2n  1
разложения ряда, была локализована на том же отрезке, что и функция                                     f (x ) ), получим
n                                        A
~                                1
f ( x)  f ( x)   ci  i ( x) , где ci   f ( x)  i ( x)dx - коэффициенты Эрмита.
i 0                    k A
~
Далее для приближенной функции                                f ( x)       применяется метод фазовой корреляции, описанный выше.
Так как функции Эрмита являются собственными функциями преобразования Фурье с собственными
значениями  1,  i , то F ( n )  i n n . Следовательно, аппроксимация преобразования Фурье имеет
n                      n                     n
следующий вид: F [ ~]  F [  c  ( x)] 
f                                                  ck F[ k ( x)]   ck i k  k ( x) ,
k k
k 0                   k 0                  k 0
т.е. вещественная часть Re  c (1) k  ( x) и мнимая часть Im  i   c
 2k                                         2 k 1 (1)  2 k 1 ( x)
k
F            2k                     F
k                                                           k
преобразования Фурье быстро считаются через коэффициенты Эрмита. Фаза                                           F   вычисляется по формуле
Im F .
 F  arctg
Re F
Далее вычисляется функция                            RFG (u )          и для подсчета POC-функции               POC fg (x)    используется
свойство     вычисления                  исходной             функции          через   дважды      примененное       преобразование           Фурье:
____________________________________________________________________________________________
168                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

f ( x)   F [ F (u )] . Таким образом, раскладывая функцию R FG (u ) в ряд по тому же количеству   n
функций Эрмита и считая преобразование Фурье от полученного разложения, получается искомая POC-
функция: POC fg ( x)  F[RFG (u)] .
Выбор количества n функций Эрмита, в ряд по которым раскладываются исследуемые функции
интенсивности радужной оболочки, обоснован ниже в таблице.
4. Метод проекционной фазовой корреляции в задаче идентификации человека по радужной
оболочке глаза
На стадии предобработки на изображении выделяется радужная оболочка глаза и переводится в
прямоугольное нормализованное изображение [5] (рис. 1 (а)). В качестве исходной функции f (x ) в
методе фазовой корреляции рассматривается функция усредненной интенсивности нормализованного
изображения по строкам верхней четверти нормализованного изображения (рис. 1 (б)). Верхняя часть
нормализованного изображения соответствует областям радужной оболочки, близким к зрачку. Эти области
содержат наиболее информативную структуру радужной оболочки и подвержены меньшим наложениям век
и ресниц, чем нижние области нормализованного изображения. Чтобы избежать возможных искажений при
определении зрачка глаза, первые 4 строки в усреднении не участвуют.

(а)

(б)
Рис. 1. (а) Исходные нормализованные изображения одной радужной оболочки; (б) Аппроксимация
функций усредненной интенсивности n  70 функциями Эрмита: тонкие линии – исходные функции f ,
толстые линии - функции f .
~

(а)

(б)
Рис. 2. (а) Исходные нормализованные изображения разных радужных оболочек; (б) Аппроксимация
функций усредненной интенсивности n  70 функциями Эрмита.
Поворот глаза соответствует циклическому сдвигу нормализованного изображения, поэтому, определив
позицию пика POC-функции для двух изображений одного глаза, автоматически находится угол поворота
одного глаза относительно второго. На рис. 3 (а) видно, что в случае изображений одной радужной
оболочки (рис. 1), пик POC-функции выражен четко, при этом позиция пика соответствует сдвигу второго
изображения относительно первого, а высота пика определяет меру похожести функций. В случае
изображений разных радужных оболочек (рис. 2) нет четкого пика POC-функции (рис. 3 (б)).

(а)                                                 (б)
Рис. 3. POC-функции изображений (а) одной (б) разных радужных оболочек.
5. Результаты
Для оценки качества предложенного метода исследовались изображения глаз базы данных CASIA-IrisV3
[6]. В данной работе использовалась часть базы данных CASIA, содержащая 600 глаз, на область
исследования которых не попадают веки и ресницы. В этой базе данных к каждой паре изображений
применялся описанный выше метод проекционной фазовой корреляции на основе проекционного метода
Эрмита. В качестве меры похожести пары изображений рассматривалась величина пика POC-функции.
Для оценки качества метода используется анализ величины EER (Equal Error Rate), где EER – величина
ошибки работы метода, при которой ошибки I и II рода равны. В таблице приведена зависимость EER от
числа n функций Эрмита, в ряд по которым раскладываются функции интенсивности изображений.
n            60      70            80           100     128
EER          0.02    0.013         0.0149       0.026   0.03
EER в %      2%      1.3%          1.5%         2.6%    3%
Таблица. Зависимость EER от числа n функций Эрмита.
Из таблицы видно, что оптимальным числом функций Эрмита для метода фазовой корреляцией является
n  70 .
6. Заключение

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              169
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

При сравнении предложенного метода проекционной фазовой корреляции с методом фазовой
корреляции POC в задаче идентификации по радужной оболочке глаза, предложенный метод показал
значительно лучшие результаты. В то же время, предложенный метод позволил значительно ускорить
расчет фазовой корреляции по сравнению с методом, вычисляющим фазовую корреляцию через
преобразование Фурье.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной
России» на 2009–2013 годы и гранта РФФИ 10-07-00433-a.

Литература
1. R. P. Millane, W. H. Hsiao The basis of phase dominance // Optics Letters v. 34, No. 17, p. 2607-2609, 2009.
2. S. Nagashima, K. Ito, T. Aoki, H. Ishii, K. Kobayashi High Accuracy Estimation of Image Rotation using 1D
Phase-Only Correlation // IEICE Trans.Fund. v. E92-A, p. 235-243, 2009.
3. K. Miyazawa, K. Ito, T. Aoki, K. Kobayashi, H. Nakajima A Phase-Based Iris Recognition Algorithm// LNCS
(ICB 2006), No. 3832, p. 356--365, 2006.
4. A. Krylov, D.Korchagin Fast Hermite Projection Method //LNCS, v.4141, p.329-338,2006.
5. Е. Павельева, А. Крылов Поиск и анализ ключевых точек радужной оболочки глаза методом
преобразования Эрмита //Информатика и ее применения, т.4, в.1, с.79-82,2010.
6. База данных CASIA-IrisV3 http://www.cbsr.ia.ac.cn/IrisDatabase.htm

PROJECTION PHASE ONLY CORRELATION METHOD FOR IRIS RECOGNITION
Pavelyeva E., Krylov A.
Lomonosov Moscow State University,
Faculty of Computational Mathematics and Cybernetics,
Laboratory of Mathematical Methods of Image Processing

An iris recognition method using phase only correlation approach has been designed. The proposed method en-
hances the POC (Phase-Only Correlation) method [1]. In our work the phase function is calculated using Hermite
Projection Method [2]. By this method the signal is expanded into series of Hermite functions being the
eigenfunctions of the continuous Fourier Transform.
The tests with CASIA-IrisV3 database show the effectiveness of the proposed method.
The work was supported by federal target program “Scientific and scientific-pedagogical personnel of innovative
Russia in 2009-2013” and RFBR grant 10-07-00433-a.
Literature
1. K. Miyazawa, K. Ito, T. Aoki, K. Kobayashi, H. Nakajima A Phase-Based Iris Recognition Algorithm// LNCS
(ICB 2006), No. 3832, p. 356--365, 2006.
2. A. Krylov, D. Korchagin Fast Hermite Projection Method//LNCS, v. 4141, p. 329-338, 2006.
3. Iris CASIA-IrisV3 database http://www.cbsr.ia.ac.cn/IrisDatabase.htm .


ФРАКТАЛЬНЫЕ МЕТОДЫ ВЫДЕЛЕНИЯ ЛОКАЛЬНЫХ ПРИЗНАКОВ САМОПОДОБИЯ НА
ЦИФРОВОМ ИЗОБРАЖЕНИИ
Привезенцев Д.Г., Жизняков А.Л.
ГОУ ВПО «Владимирский государственный университет (филиал) Муромский институт»

В настоящее время активно развиваются новые подходы к цифровой обработке изображений. Одним из
таких является обработка изображений с использованием фрактальных методов. Успешному применению
фракталов в обработке изображений способствует тот факт, что фрактал – это самоподобный объект,
обладающий свойством инвариантности к масштабу рассмотрения, кроме того установлено, что многие
изображения можно в некоторой степени считать фракталами или мультифракталами. Вычисляемые
фрактальные характеристики используются для качественного или количественного описания изображений,
и обладают очень полезным свойством – инвариантностью к изменению масштаба фрактального объекта и к
положению объекта.
Любая из процедур обработки изображений опирается на модель класса изображений –
формализованное описание, выполненное с определенной степенью абстрагирования. Роль модели
изображения в процессе извлечения информации состоит в обеспечении адекватного описания
существенных свойств класса изображений, позволяющего дать конструктивную основу, для построения
эффективных вычислительных процедур.
В частности для создания новых методов фрактальной обработки изображений необходимо создать
фрактальную модель изображений. Одним из возможных подходов к фрактальному описанию изображения

____________________________________________________________________________________________
170                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

является использование систем итерируемых функций, которые используются для построения регулярных
фрактальных объектов.
Системы итерируемых функций состоят из метрического пространства X и конечного множества
сжимающих отображений wn:XX с коэффициентами сжатия sn [1]. Если в качестве метрического
пространства X представить изображение f, а в качестве сжимающих отображений wn операторы,
переводящие одни участки изображения в другие, получим некоторое фрактальное описание изображения.
Описанное таким образом изображение разбивается на множество неперекрывающихся ранговых блоков и
множество перекрывающихся доменных блоков. Для каждого рангового блока подбирается доменный блок
и соответствующее сжимающее отображение, переводящее данный доменный блок в ранговый.
С учетом этого структуру изображения можно записать следующим образом:

                          
M
f    Bnii ,mi * Ai Bkdii,li  f   Ci
r
(1)
i 1
где f – исходное изображение, Bkn.m – оператор извлечения блока изображения, (Bkn.m)* - оператор вставки
блока, которые соответственно извлекают и вставляют блок изображения размером k×k, левый нижний угол
которого находится в точке с координатами (n,m), Ai, Ci – коэффициенты преобразования яркости
изображения.
Выражение (1) можно назвать фрактальной моделью изображения. Каждое изображение описывается
параметрами этой модели, составляющими фрактальный код изображения Ф:
       Список ранговых блоков R={Ri}; Ri={ri,ni,mi};
       Список доменных блоков D={Di}; Di={di,ki,li};
       Соответствующие преобразования T={Ti}; Ti={Ai,Ci}.
Данная модель позволяет представить изображение как множество фракталов, т.е. фактически
подчеркивает локальные свойства самоподобия.
Если для построения фрактальной модели изображений использовать один доменный блок, который
соответствует самому изображению, то состав и параметры получившихся ранговых блоков характеризуют
свойства самоподобия изображения. При таком разбиении, чем больше размер рангового блока, тем
большего размера является участок изображения подобный самому изображению. Чем меньше средний
размер ранговых блоков, тем больше изображение содержит мелких деталей, подобных целому
изображению.
Однако часто изображение не содержит участков подобных самому изображению, но содержит одни
участки, которые подобны другим участкам. Поэтому если для фрактального разложения на блоки
использовать не один доменный блок, тогда состав и параметры ранговых блоков будут отражать
проявление свойства подобия внутри изображения.
Ввиду широкого распространения цифровой техники человеку приходится работать с большим
количеством различных изображений. Поэтому, одной из важных задач в цифровой обработке изображений
является задача распознавания, включающая выделение характерных признаков и особенностей
изображений, на основе которых можно осуществлять их классификацию, упорядоченное хранение и поиск.
Предполагается применить фрактальную модель для решения этих задач. Использование методов и
средств фрактального представления изображений позволяет решать эту задачу на основе сравнения
промежуточного цифрового описания изображений, а не исходного. Здесь фрактальное представление,
выявляющее локальные признаки самоподобия внутри изображения, позволит выделить его основное
содержание.
Основной задачей здесь является разработка алгоритмов получения набора признаков для фрактальных
представлений изображений, однозначно определяющих их среди других, разработка алгоритмов
формирования фрактальных представлений классов, а также разработка алгоритмов вычисления разности
фрактальных представлений и их признаков.
Исходя из того факта, что большинство реальных изображений не полностью самоподобно, следует, что
для формирования фрактального кода используются не все участки изображения, т.е. не все доменные блоки
необходимы для восстановления изображения.
Таким образом, из множества доменных блоков D={Di} можно выделить подмножество используемых
доменных блоков Du={Dui}, отражающие самоподобные участки изображения (Рис.1).

____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              171
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

Рис.1 – используемые доменные блоки

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

Рис.2 – частота использования доменных блоков во фрактальном коде
Такое представление строится для каждого уровня доменных блоков. Объединив сведения о частоте
использования каждого участка изображения можно построить общую гистограмму использования
доменных блоков, которую целесообразно выполнить в виде трехмерной поверхности. Для лучшего
отражения непрерывного изменения исследуемой величины на трехмерных поверхностях обычно
применяются изолинии, которые отражают распределение этой величины на поверхности. Наложив
полученный рельеф на исходное изображение, можно визуально отследить степень проявления самоподобия
на изображении и выявить отдельные наиболее подобные участки, который также можно использовать в
качестве локальных признаков самоподобия изображения (рис. 3).

Рис.3 – изображение изолиний поверхности, отражающей самоподобные участки изображения
Представленная фрактальная модель цифрового изображения позволяет описать изображение с точки
зрения самоподобия. Следовательно, с помощью этой модели можно оценить различные локальные
характеристики подобия на изображении, которые, в частности, могут использоваться как информативные
признаки изображения в задачах распознавания.
Литература
1. Уэлстид С. Фракталы и вейвлеты для сжатия изображений в действии. Учебноепособ. – М.:
ИздательствоТриумф, 2003 – 320 с.: ил.

____________________________________________________________________________________________
172                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

2. Прэтт, У. Цифровая обработка изображений : в 2 кн. / У. Прэтт ; пер. – М.: Мир, 1982.
3. Новейшие методы обработки изображений. / Под ред. Потапова А.А. М.:ФИЗМАТЛИТ 2008 – 496 с.,
ил.
4. Petrou M.,Bosdogianni P. Image processing : The fundamentals, John Wiley and Sons, 1999 г. ISBN: 0-471-
99883-4.
5. Жизняков А.Л. Формирование и анализ наборов признаков многомасштабных последовательностей
цифровых изображений. / Программные продукты и системы. 2007, №4, С.24.
6. Жизняков А.Л., Фомин А.А. Классификация изображений микроструктур металлов на основе
многомасштабных моделей. / Фундаментальные проблемы современного материаловедения. 2007, Т.4, №2,
С.75-80.

FRACTAL METHODS OF SELECTION OF LOCAL FEATURES OF SELF-SIMILARITY ON THE
DIGITAL IMAGE
Privezentsev D.,Zhiznyakov A.
The Murom institute of Vladimir state university

Now new approaches to digital image processing actively develop. One of such is image processing with usage
of fractal methods. For creation of new methods of fractal image processing it is necessary to create fractal model of
images. One of possible approaches to the fractal description of the image is usage of systems of iterated functions
which are used for creation of the regular fractal objects.
The fractal model of the image can be written down as follows:

                          
M
f    Bnii ,mi * Ai Bkdii,li  f   Ci
r
(1)
i 1
where f – the source image, Bkn.m – the operator of extraction of the unit of the image, (Bkn.m)* - the operator of an
insertion of the unit which accordingly extract and insert the unit of the image in the size k×k which left lower cor-
ner is in a point with coordinates (n, m), Ai, Ci – coefficients of conversion of image brightness.
It is supposed to apply fractal model to the decision of tasks of recognition. Usage of methods and means of frac-
tal representation of images allows solving this task on the basis of comparing of the intermediate digital description
of images. Here the fractal representation revealing local signs of self-similarity in the image, will allow selecting its
The primary goal is algorithm elaboration of obtaining of a feature set for fractal representations of images, un-
ambiguously defining them among others, algorithm elaboration of formation of fractal representations of classes,
and also algorithm elaboration of calculation of a difference of fractal representations and their signs.
The presented fractal model of the digital image allows to describe the image from the point of view of self-
similarity. Hence, by means of this model it is possible to estimate various local characteristics of similarity on the
image which, in particular, can be used as the informative signs of the image in recognition tasks.


ИНТЕРПОЛЯЦИЯ КАДРОВ В ЗАДАЧАХ ПРОСТРАНСТВЕННО-ВРЕМЕННОЙ КОРРЕКЦИИ
СТЕРЕОВИДЕОПОСЛЕДОВАТЕЛЬНОСТЕЙ
Пьянков Д.И.
ГОУ ВПО «Сибирский государственный аэрокосмический университет
имени академика М.Ф. Решетнева»

В настоящее время стереосъемка в домашних условиях возможна в полуавтоматическом режиме:
необходимо синхронизировать стереопары вручную вследствие несинхронного старта видеокамер,
различного времени доступа к устройству хранения. Разработка методов пространственно-временной
коррекции стереовидеопоследовательностей позволит решить эти проблемы. Одним из методов является
коррекция видеопоследовательностей на основе преобразовании частоты кадров с помощью интерполяции.
Преобразование частоты кадров видеопоследовательностей (FRC, Frame Rate Conversion) используется
при сжатии изображений, преобразовании форматов видео, в алгоритмах повышения качества
видеоматериалов. Использование FRC делает движение объектов плавным и более приятным для зрения. С
помощью данного алгоритма возможно восстанавливать поврежденные кадры в видеопоследовательности.
Методы преобразования частоты кадров делятся на две группы: методы, не учитывающие информацию о
движении, и методы, адаптированные к вычислениям и обработке информации о движении [1]. К методам
преобразования частоты кадров на основе интерполяции без оценки движения можно отнести метод
повторения кадров (Frame Repetition) и метод линейной интерполяции (Linear Interpolation) [2]. Данные
методы не учитывают информацию о движении и являются простыми и быстрыми, но менее эффективными,
чем другая группа методов преобразования частоты кадров. Методы преобразования частоты кадров на
____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              173
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

основе интерполяции с оценкой движения основаны на интерполяции скомпенсированных кадров,
поскольку используют оценку движения (ОД, Motion Estimation) для компенсации перемещений всех
движущихся объектов на новых кадрах на нужный момент времени согласно их скорости. Исследования
проводятся по трем основным направлениям: постобработка векторов движения, маскирование артефактов
и обработка наложений. Ввиду особой важности задачи ОД в настоящее время разработано большое
количество соответствующих алгоритмов, которые можно разделить на следующие группы: методы фазовой
корреляции (Phase Correlation), методы оптического потока (Optical Flow), методы сопоставления блоков
(BMA, Block-matching algorithms).
Для наиболее точного получения интерполированных кадров в видеопоследовательности необходимо
использовать методы фазовой корреляции. Однако, поскольку это довольно медленный метод, оптимальным
выбором является алгоритм сопоставления блоков. Данный метод является быстрым, немного уступая
методу фазовой корреляции по точности интерполяции.
Получение интерполированных кадров. Пусть дана видеопоследовательность Lin(t), t{0,1,2,3,…}.
Необходимо вставить между каждыми двумя кадрами определенное число кадров, заданное параметром n.
Тогда видеопоследовательность Lout(z) после преобразования частоты кадров будет представлена
выражением:
 Lin (t ) , z  n  t , t  {0,1,2,3,. ..},
Lout ( z )                                                     (1)
 Lin (t ) , z  n  t  l , l  {1,2,3,..,n  1},
z  {0,1,2,3,.   ..}, n  {2,3,4...},
где n - параметр, определяющий, во сколько раз нужно увеличить частоту кадров; z – общее число кадров.
Общая схема получения интерполированных кадров:
1. Рассматривается текущий кадр L(t) и следующий кадр L(t+1) видеопоследовательности Lin(t),
t{0,1,2,3,…}.
2. Производится оценка движения согласно выбранному алгоритму.
3. Вычисляется поле векторов движения.
4. Производится интерполяция между текущим кадром L(t) и следующим кадром L(t+1).
5. Строятся интерполированные кадры.
Далее проводится оценка движения для нахождения поля векторов для блоков кадров (рис. 1) [3].

Рис.1. Направление векторов движения
Алгоритм интерполяции кадров. После того, как были найдены векторы движения, производится
линейная интерполяция кадров поблочно. Общая схема интерполяции кадров:
1. Рассматривается текущий кадр L(t) и следующий кадр L(t+1) видеопоследовательности Lin(t),
t{0,1,2,3,…} с найденным вектором движения move для блока B(x,y).
2. Пусть n – количество интерполированных кадров, которые необходимо вставить между кадрами L(t) и
L(t+1). Рассматривается каждый пиксель С1 блока B(x,y) на текущем кадре L(t) и каждый пиксель С2 блока
B(x,y) на следующем кадре L(t+1).

____________________________________________________________________________________________
174                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference
Обработка и передача изображений
____________________________________________________________________________________________

3.
4.   Рис.2. Линейная интерполяция блоков кадра при количестве интерполированных кадров n = 2

Согласно направлению вектора движения на каждом интерполированном кадре, данный пиксель плавно
изменяет свое значение каждой цветовой компоненты (Y,U,V) Сn. Для получения плавных
интерполированных кадров (рис.2) используется выражение, представляющее собой линейную
интерполяцию значений пикселей:      Cn  C1  (C2  C1 )  t , t  [0;1] ,
(2)
где C1 - значение цветовой компоненты Y пикселя (x,y) блока B(x,y) на текущем кадре L(t); C2 - значение
цветовой компоненты Y пикселя (x,y) блока B(x,y) на следующем кадре L(t+1); Cn - значение цветовой
компоненты Y пикселя (x,y) блока B(x,y) на интерполированном кадре n; t – параметр, задающий линейное
смещение пикселей на каждом интерполированном кадре. Параметр t вычисляется по формуле:
m                                                 (3)
t    , m {1,2,3,...,n} ,
n 1
где n – количество интерполированных кадров. Параметр t[0;1], что дает плавное смещение цветовых
компонент для каждого интерполированного кадра.
5. Шаг 2 повторяется для цветовых компонент U, V.
6. На кадрах L(t) и L(t+1) берется следующий блок. Повторяются шаги 2-3 до тех пор, пока все блоки не
будут интерполированы. Таким образом, будет получено n интерполированных кадров.
Представленный метод был реализован в среде программирования «Visual Studio 2008» на языке C++ в
виде плагина для фреймсервера «Avisynth». Плагин определяет вектора движения на каждом кадре
видеопоследовательности. Сохранив файл, получится интерполированная видеопоследовательность, частота
кадров которой возрастет в заданное количество раз. Переход между кадрами станет более плавным за счет
интерполяции.
Для оценки работы метода был разработан плагин, осуществляющий тестирование полуавтоматическом
режиме. Для тестирования использовались две видеопоследовательности: foreman, где происходит быстрое
движение; container, где происходит медленное движение. Плагин позволил собрать данные о среднем PSNR
(peak signal-to-noise ratio) интерполированных кадров каждой видеопоследовательности от текущего кадра к
следующему (PSNR1) и наоборот (PSNR2). PSNR включает в себя среднеквадратичное отклонение и
определяется следующим выражением: psnr ( z )  20 log MAX ,                                  (4)
10
MSE
где psnr(z), z{2,3,…n} – значение схожести между интерполированным кадром L(z) левой
видеопоследовательности и кадром R(t) правой видеопоследовательности;       MAX[0;255] - максимальное
значение, принимаемое пикселем изображения; MSE – среднеквадратичное отклонение для каждой цветовой
компоненты пикселя.
Исследования производились по следующим параметрам: среднее PSNR интерполированных кадров
каждой видеопоследовательности, время выполнения метода. Использовалось множество вариантов
изменения различных параметров метода: размер блока, смещение между блоками, число
интерполированных кадров.
Для быстрой видеопоследовательности foreman параметр PSNR меньше, чем для медленной
видеопоследовательности container. Это объясняется тем, что в быстрой видеопоследовательности, за счет
резкой смены кадров, будет больше шума, чем в медленной видеопоследовательности. Проведенные
эксперименты быстрой видеопоследовательности foreman показали, что лучшая точность получается при
большом блоке 1616, который попадает в области движения и интерполирует корректно. Для обеих
____________________________________________________________________________________________
Цифровая обработка сигналов и ее применение                                              175
Digital signal processing and its applications
Обработка и передача изображений
____________________________________________________________________________________________

видеопоследовательностей время выполнения алгоритма интерполяции одинаково, так как их размеры
одинаковы. При смещении между блоками в ширину блока время выполнения алгоритма уменьшается при
увеличении размера блока.
Предложенный метод позволяет получать интерполированные кадры в видеопоследовательностях и
применяться в задачах стереовидения. Однако интерполяция данного метода не устойчива к повороту,
масштабированию и сдвигу объекта при оценке движения. Это проблему планируется решить, рассмотрев
моменты изображений Зернике.
Литература
1. Горошкин, А.Н., Пьянков, Д.И. Обзор методов пространственно-временной коррекции
видеопоследовательностей в задачах стереовидения // В трудах I всероссийской научной конференции
молодых ученых "Теория и практика системного анализа", в Т. 2. - Т.1. - Рыбинск: РГАТА, 2010. - 206 с. - с.
105-111.
2. Ханова, А. А.. Интерполяция функций. // Методическое пособие для студентов Института
информационных технологий и коммуникаций. - Астрахань: Изд-во АГТУ, 2001. – с. 15-16.
3. Пьянков, Д.И., Горошкин, А.Н. Оценка движения для интерполяции кадров в задачах стереовидения.
// В материалах 16-й Международной науч.-техн. конференции «Проблемы передачи и обработки
информации в сетях и системах телекоммуникаций». - Рязань: Рязанский государственный
радиотехнический университет, 2010. - 220 с. - с. 151-152.

FRAME INTERPOLATION IN TASKS OF SPATIO-TEMPORAL CORRECTION OF STEREO-VIDEO
SEQUENCES
Pyankov D.I.
Siberian State Airspace University after academician M.F. Reshetnev (SibSAU)

Currently, stereo filming home-MADE is possible in the semi-automatic mode: it is necessary to synchronize the
stereo cameras manually due to non-synchronous launch of cameras and different access time to storage devices.
Development of methods for spatial-temporal correction of stereo-video sequences will solve these problems. One
of methods is correction of video sequences based on the frame rate conversion by interpolation.

____________________________________________________________________________________________
176                                                  Доклады 13-й Международной конференции
Proceedings of the 13-th International Conference

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 14 posted: 11/9/2012 language: Unknown pages: 26