ЧИСЛЕННЫЕ МЕТОДЫ И АЛГОРИТМЫ ПОСТРОЕНИЯ ПОЛЕЙ НАПРАВЛЕНИЙ КВАЗИПЕРИОДИЧЕСКИХ СТРУКТУР
Н.Ю. Ильясова, A.B. Устинов. А.Г. Храмов* Институт систем обработки изображений РАН, *) Самарский государственный аэрокосмический университет
Аннотация
Описываются новые численные методы построения поля направления, приводится полная классификационная схема методов и алгоритмов построения поля направлений, описание ее классов. Дается сравнительный анализ алгоритмов по точности и быстродействию в отдельных классах, а также между" классами указанной схемы. Приводятся примеры практических задач анализа изображений, решаемых с использованием поля направлений (анализ дактилограмм, кристаллограмм биологических сред).
1. Введение кровяной плазмы и многие другие изображения ее-
Представителями* изображений, характеризуе- тественного и искусственного происхождения
мых наличием квазипериодических структур, яв- <Рис1)- Квазипериодические структуры определя-
ляются интерферограммы, дактилограммы,' крис- ются наличием многокотурнои упорядоченной
таллограммы слезной жидкости, крисгаллограммы текстуры с вьраженной ориентацией.
а) б) в) г)
Рис.1. Иллюстрация различных квазипериодических структур: а) дактилограмма, б) кристаллограмма слезной жидкости, в) кристаллограмма плазмы крови, г) интерферограмма.
В каждой малой области квазипериодичсской структуры функция яркости является периодической вдоль определенного направления и визуально выражается в системе параллельных полос. Основными параметрами таких изображений являются преимущественное направление полос и их густота в каждой точке изображения. Кроме того, на изображениях выделяются особые точки и линии (сингулярности), в которых нарушается периодическая структура. Поле направлений является удобным носителем информации для решения задач интерпретации, распознавания и диагностики изображений с квазипериодической структурой. Для формирования признаков данного класса изображений предлагается использовать подход, основанный на вычислении поля направлений.
2. Метод поля направлений
Впервые идея метода поля направлений для интерпретации и анализа диагностических изображений была предложена в работе [1] применительно к анализу дактилограмм. В работах [2, 3] метод использовался для диагностики кристаллограмм слезной жидкости [4]. Несколько позднее в работе [5] были впервые представлены теоретические аспекты метода поля направлений и определены некоторые
подходы к построению алгоритмов его расчета (оптический метод и др.).
Под полем направлений будем понимать поле углов преимущественного направления полос в локальной окрестности точки изображения. Угол направления полосы Ц1(х,у) в данной точке равен по определению углу касательной к линии уровня функции яркости (рис.2).
Традиционные методы обработки изображений (линейная обработка окном, нелинейная фильтрация, и т.п.) не могут бьпъ непосредственно применены к полям направлений. Это связано с особенностями арифметики направлений, которая является периодической по значениям, причем период равен я, а не 2п, как в векторной геометрии. В [5] показано, как можно решить некоторые проблемы арифметики направлений, позволяющие применять традиционные методы обработки для полей направлений.
Комплексное поле направлений определяется
Рис.2. Описание поля направлений.
как:
v(x,y) = w(x,y)exp(j2y(x.y)).
(1)
отгппгиг»«»« и
Рис.3. Стандартные конфигураций полос и соответствующие им поля направлений.
методы расчета поля направлений
параметрическая
аппроксимация
дифференциальные
ЧЕТЫРЕХ'-ГРАДАЦИОННЫЙ ПЕРЕБОРНЫЙ
обработка маской 3x3
члтырех-градациоиныйс интерполяцией
параболическая аппроксимация
обработка маской 5x5
ЦИЛИНДРИЧЕСКАЯ АППРОКСИМАЦИЯ
¡ые методы I
V I
МЕТОД основной ЛОКАЛЬНОЙ ЧА! готъ:
обработка токальными масками (бпф)
квадратная маска
РАЗНОСТНЫЕ СХЕМЫ
КРЕСТОВИДНА» масха
СИНУСОИДАЛЬНЫЕ АППРОКСИМАЦИИ
УСРЕДНЕНИЕ ПО УГЛАМ
УСРЕДНЕНИЕ ПО ЧАСТОТАМ
Рис.4. Классификационная схема методов построения поля направлений.
Весовая функция м>(х,у) в (1) имеет физический смысл достоверности (выраженности, надежности определения) поля направлений в данной точке. Очевидно, она должна иметь верхнее граничное значение, равное 1, на участках изображения с отчет ливо выраженным направлением периодической структуры. Такими участками являются области с резкими перепадами яркости (контуры). Минимального значения, равного нулю, весовая функция должна достигать на участках с постоянной яркостью. На рис.3 приведены примеры нескольких видов стандартных конфигураций полей направлений (шаблонов), используемых в задачах диагностики. Поле направлений изображено градациями серого, причем нулевым значениям соответствует черный цвет, а значению я - белый цвет.
Непосредственное использование численного дифференцирования для цифрового построения поля направлений возможно лишь для "гладких" изображений при условии полного отсутствия шумов наблюдения. На практике необходимо использовать помехоустойчивые методы построения поля направлений, основанные на аппроксимации и усреднении. В статье представлены различные подходы к построению алгоритмов расчета поля направления, проведена их систематизация и сравнительное исследование их точности и быстродействия.
3. Классификация методов построения поля поправлений
Используя различные модели представления поля направлений, выделим пять обширных классов методов его построения: методы параметрической аппроксимации и локальных градиентов, про-скционно-дисперсионные, дифференциальные и спектральные методы. На рис.4 представлена классификационная схема методов построения полей направлений, разработанная в ходе экспериментальных исследований. Целью классификации является выбор наилучших методов для конкретных приложений.
3.1. Методы параметрической аппроксимации
Класс методов параметрической аппроксимации основан на аналитическом представлении функции яркости в симметричном прямоугольном скользящем окне. Пусть 1(х,у) - функция яркости изображения в точке (х,у) е G, где G некоторая область изображения. Аппроксимация изображения 1(х,у) в симметричном прямоугольном окне W={(u,v): -\i < и <М, -М < v <М} представляет собой двумерный полином степени Р:
¡(x-u,y-v)= У Oj Ax,y)u'vJ,(x,y)eG,(u,v)fzW,
ЬЛ Q
где 0= {(/, j)0 <. i й Р,0 <, jüP,Ozi + jüP) множество показателей степеней. Коэффициенты полинома а^{х,у) являются функциями положения центра окна (х,у).
Рассматриваемая локальная аппроксимация функции яркости позволяет сгладить изображение для вычисления производных и получить некий геометрический инвариант в виде поверхности первого или второго порядка. В качестве аппроксимирующей функции используем следующие поверхности:
плоскость: l(u,v)-au+ bv + c \ квадратичную поверхность:
l(u,v) = au2 + bv2 + cuv + du+ev+/; цилиндрическую поверхность:
l(u,v)= а{у + kti)2 + b(v + ku) + c.
Коэффициенты аппроксимирующих полиномов определим по методу наименьших квадратов из условия минимума функционала, представляющего собой среднеквадратичную ошибку аппроксимации;
J = jj[/(x - и,у - v)- ¡(u.v^dudv ->■ min .
(u.v)EiT at>-c-
В дискретном виде функционал имеет следующее выражение:
J = Xx[jr(x-M>V-v)-/(w,v)}'arM£/v-> min .
(,U.vajbfi...
Приравнивая к нулю частые производные функционала по полиномиальным коэффициентам, получим систему линейных уравнений относительно неизвестных коэффициентов: аппроксимация плоскостью:
¿Sv- с^1 = \Уоо (2)
где ц^Дх.у)- локальные степенные моменты изображения в окне обработки:
VpqM= Tupv4{x-u,y-v), (3)
(a.vjeW
быстрые алгоритмы рекурсивного вычисления моментов приведены в [6];
аппроксимация квадратичной поверхностью;
<с^и2\2 = =Ию,еТч2 = И-о// (4)
а^и2 +/•£! = »оо'>
аппроксимация цилиндрической поверхностью:
а^У +бвЛ£и V +ак4£и* +с£у2 ^ск2У,и2 =
= ц02+2кци+к211201
ЬЪ»2+Ьк2Ъи2=Ио1+к»10> (5)
а^У +ак21"2 +с^1 = [100, 13)
ба2к^и2у2 + 2а2к3£и4 + 2аскУи2 +Ь2к£и2 = = 2а[1,, +2ак\х20+Ьц10.
Аппроксимация плоскостью. В системе координат, связанной с положением скользящего окна И', зададим одномерный полином / (х,у). аппроксимирующий изображение 1(х,у) в виде плоскости:
¡(х, у) = ах + Ьу + с{х, у) 6IV , коэффициенты которого имеют следующий вид:
а= 6У-ю ь= 6У-01
м(м+фм+¡)' м{м+фм +7)'
Используя определение поля направлений и систему (2) получим для его определения следующее вьфажение:
Квадратичная аппроксимация. Метод локальной квадратичной аппроксимации основан на аппроксимации функции 1(х,у) в пределах окна IV квадратичной поверхностью:
= ах2 + Ьу2 + сху + сЬс + еу + у)е IV (6)
Рис.5. Линии уровня квадратичной поверхности
За значение поля направления в точке (х,у) принимается направление главной оси параболоида
(6) (рис.5), когда центр окна IV находится в точке (х,у). Численное решение данной задачи приведено в [6), где предложен эффективный алгоритм вычисления коэффициентов с использованием рекурсив-ного метода вычисления моментов (3). Дня определения направления главной оси параболоида перейдем к каноническому виду уравнения поверхности
[7]:
/(У,У)=ОУ2+/>У2+/\
Для этого используем следующие аффинные преобразования координат:
\х = х соха- у зта+ х0, [у = х' ш а+ у' а>5 а+ у 0,
где
Хп =■
се - 2М
уу
4аЬ - с' сс1 - 2ае 4аЬ-с2
0)
являются координатами вершины аппроксимирующего параболоида, а угол а определяется из следующего уравнения: tg2ct = с/(а-Ь).
В результате решения данного уравнения получим два взаимно перпендикулярных направления аI и а2, соответствующие главным осям квадратичной поверхности (6). Выбор направления, соответствующего большей полуоси, осуществляется следующим образом: а = ш*£ тт Ы, где
а' = асой2а + Ьят2 а ч-с^шасо^а .
Данное выражение получено при переходе к каноническому виду уравнения поверхности с помощью аффинного преобразования координат (7). Вектор параметров аппроксимирующего полинома найдем из системы (4):
_1^20_
м(м + фм + 1У[4М2 + 4М-3)
_Ну-оо_
(2М + 1)2(4М2 + 4М-3}'
а =
Ъ =
_4^02_
М(М + фм + ¡У [4М 2 + 4М - 3)
15\Уоо_
4М-3)'
с =
(2М + ])2{
4М2 + ■
М2^ + 1)2(2М + 1)2 ''
е -
м(м + фм + ¡У '
м{м+фм + ])2 '
\xJl4M2 + 14М-з)-15{[120+^02) (2М + 1У(4М2 +4М -З)
/
Цилиндрическая аппроксимация. Общая идея метода предложена в [8]. Мы рассмотрим построение системы нелинейных уравнений в явном виде. В системе координат скользящего окна рассмотрим полином второго порядка одной переменной, который повернут на плоскости изображения на некоторый угол у, аппроксимирующий изображение 1(х,у) в виде цилиндра 1{х,у) =о{у +Ь)2+Ь(у , где
у. За направление в центре окна (х,у) принимаем направление, вдоль которого повернут цилиндр, а именно, угол Параметры а, Ь, с, к аппроксимирующего полинома определим из системы
(5):
М(М + фМ+1У(1+к2)' (8)
с = Ш-¿А-¿, 6а2к*22+2а2к3з4 +
+ 2аск$2 +Ь2кяг = +2ак\120 +Ь\х10,
где =
[х.у^УУ
•У у*,
¿2 ~ 5 20 ~ 302 ~
ъ22
_м{м+фм + 1)2
3
м2(м +1)2(2М +¡У
9
_ м(м + фм + 1)2(зм2 + зм-])
5 4 - х40 ^04--у^-■
В отличие от аппроксимации плоскостью и квадратичной поверхностью система (5) является нелинейной, что требует использования для ее решения соответствующих численных методов.
Метод комбинированной аппроксимации. При проведении исследований рассмотренных методов было замечено, что при аппроксимации изображения плоскостью возникает большая погрешность в точках экстремума функции яркости, а при аппроксимации квадратичной поверхностью - в точках перегиба. Из рис.6 видно, что в точках экстремума функции (участок 1) погрешность можно уменьшить, используя параболическую аппроксимацию, в точках перегиба (участок 2) - аппроксимацию прямой (для двумерного случая - плоскостью). Таким образом, воспользуемся методом комбинированной аппроксимации. Определим две величины, характеризующие соответственно параметр крутизны и кривизны аппроксимируемой функции:
2 (дЦх.у))2 , (д!(х,у)Х
дх
ду
'(?Н*у)} 2
1 аг2 -л- эу2)
(9)
И . (Ю)
Рис.6. Метод комбинированной аппроксимации
Величина (9) достигает высоких значений в точках перегиба, а величина (10) в точках экстремума. Согласно (6) данные характеристики будут
иметь следующие выражения: - с12 + е2,
Д/(х-, у) = 4(а2 + Ь 2 )(2М +1)2 . Использование
конкретного способа аппроксимации определяется следующими условиями: аппроксимация
плоскостью: параболическая:
а) б) в) г)
Рис. 7. Иллюстрация работы методов параметрической аппроксимации: а) изображение интерферограммы; поле направлений, полученное методом: б) аппроксимации плоскостью, в) параболической аппроксимации, г) комбинированным.
Экспериментальные исследования. Класс методов локальной параметрической аппроксимации обладает высокой устойчивостью к шумам, так как метод наименьших квадратов имеет шумоподав-ляющее свойство. Недостатком же является необходимость адаптивного выбора размеров окна обработки IV (размер окна должен быть меньше расстояния между соседними полосами).
Проведем исследование рассмотренных выше методов на тестовых изображениях. Из-за специфики значений отсчетов поля направлений невозможно использовать понятие среднеквадратической ошибки в обычном смысле. Воспользуемся представлением поля направления в комплексной форме (1) с единичной весовой функцией. В данном случае под среднеквадратичной ошибкой будем понимать следующую величину:
где £> - область изображения. \у(х,у) - истинное
значение поля направления, (\>(х,у) - его оценка. Используя (1), получим:
И (x,ILd
(11)
\и\ (Х.УНО
Значение среднеквадратичной погрешности находится в диапазоне от нуля до четырех. При этом нулевое значение соответствует совпадению полей направлений, а максимальное - расхождению на 90°. В качестве значения погрешности будем использовать угловую среднеквадратическую погрешность: = агсвш^Е2 / 4 . Исследование точности проводилось на серии тестовых изображений, отличающиеся наличием шума, равномерностью и неравномерностью пространственной частоты и видом квазипериодической структуры (рис.8). Рассматривалась зависимость точности оценивания
от отношения шум/сигнал h' дисперсия функции яркости:
a2 =r^U(l(x,yyi)2dxdy. По
а„ / а.
где СУ2'
б) в) г)
Рис.8. Тестовые изображения квазипериодических структур: а) с константным полем направления; б)-г) различные виды полос; д) с шумам 50%.
Полученная ошибка зависит от уровня шума, от пространственной частоты (рад/пиксел) (или от периода анализируемой структуры г/ = 2п/со (пиксел/период)) и от метода построения поля направлений. Для полноты экспериментальных исследований данного класса методов проведем также тести-
рование переборного метода цилиндрической полиномиальной аппроксимации [8]. Сравнивая между собой методы параболической аппроксимации и аппроксимации плоскостью (рис.9, 10). можно сделать вывод, что метод аппроксимации плоскостью.
обладая наименьшей вычислительной сложностью, обеспечивает наилучшую оценку поля направления.
0,2% 0,5% 1% 2% 5% 10% ................ПЛОСКОСТЬ
—о— параболоид
--комбинированный
------переборный цилиндрический
Рис.9. Зависимость точности оценивания поля направлений аппроксимационными методами от интенсивности шума (маска 11x11, со=0.07тф.
Это объясняется тем, что на тестовом изображении площадь участков, соответствующих перепаду функции яркости (рис.6) больше, чем участков, соответствующих точкам экстремума. Можно заметить также существенное уменьшение погрешности оценки поля направлений при использовании комбинированной аппроксимации.
0.02л 0.03я 0.04я 0.07л Пространственная частота
ОЛтг 0.13л 0.17я мI рад / пиксел)
...............плоскость
—о-параболоид
-комбинированный
-------переборный цилиндрический
Рис. 10. Зависимость погрешности оценивания поля направления аппроксимационными методами от пространственной частоты квазипериодической структуры (маска 9x9).
3.2. Методы локальных градиентов
Рассматриваемый класс методов основан на том факте, что градиент функции в любой точке перпендикулярен касательной к линии уровня в этой точке. Методы локальных градиентов основаны на вычислении градиента функции яркости при различных положениях локальной маски внутри сканируемого по изображению симметричного прямоугольного внешнего окна IV размером М х N (локальный
градиент): (/хи где 1 <к ¿К 1<1 <М.
Использовались два типа локальной маски (рис.11). За направление полос в точке сканирования принимаем угол у, вычисленный на основе найденных локальных градиентов:
Ъ Ч^-У) = ~/х//у,0 < у(х,у) < тг, где =
ох оу Рассматриваемый класс методов нахождения поля направления Рис.11. можно разбить на два подкласса, характеризуемые различным способом использования локальных градиентов в общей схеме расчета поля направлений: методы усреднения проекций градиента и методы усреднения локальных углов направления. Метод усреднения проекций градиентов основан на использовании локальных градиентов
{/*■ соответствующих положению (к,Г) ло-
кальной маски , в расчете градиента функции яркости в центре внешнего окна
/ N М
мм
Метод усреднения локальных углов направлении использует локальные градиенты (/хи,/^1) для расчета локальных углов:
(12)
кг!
Тогда направление полосы в центре внешнего окна можно вычислить усредняя поле локальных углов:
к,1
. 1 N М .
(13)
Значение весовой функции поля направлений соответственно будет иметь следующее выражение:
1 N м ( \ V = — £• ММ к-Н-1
Фактически метод состоит из двух этапов: 1) определение градиента, либо направления (на основе градиента) в каждой точке локальной маски (2x2, крест); 2) сглаживание поля градиентов или поля направления в окне. Процесс классификации методов градиентов можно продолжить, базируясь на различных способах нахождения локальных градиентов: использование разностных схем и локальной аппроксимации гармоническими функциями. Общая классификационная схема для класса методов локальных градиентов представлена на рис. 12.
МЕТОДЫ ЛОКАЛЬНЫХ ГРАДИЕНТОВ |
СИНУСОИДАЛЬНАЯ АППРОКСИМАЦИЯ 1 [ РАЗНОСТНЫЕ СХЕМЫ
УСРЕДНЕНИЕ ПО ПРОЕКЦИЯМ | УСРЕДНЕНИЕ ПО ПРОЕКЦИЯМ
- КРЕСТОВИДНАЯ МАСКА) - КРЕСТОВИДНАЯ МАСКА]
" КВАДРАТНАЯ МАСКА 1 | КВАДРАТНАЯ МАСКА | |
УСРЕДНЕНИЕ ПО УГЛАМ | УСРЕДНЕНИЕ ПО УГЛАМ ^
Рис.12. Классификация методов локальных градиентов.
Градиентный метод с использованием разностных схем. В этом методе для нахождения градиента используется его определение. При вычислении частных производных применяется разностный метод со сканированием локальной маской двух видов: квадратная размером 2x2 и крестообразная. Эти методы, как и все градиентные, позволяют получить непрерывные значения угла направлений, что обеспечивает повышение точности оценивания по сравнению с переборными алгоритмами. Быстродействие методов существенно зависит от размера внешнего окна сканирования V/. Недостатком метода является повышенная чувствительность к шумам.
Методы локальной синусоидальной аппроксимации. Методы описываемого подкласса градиентных методов основаны на аппроксимации в пределах локального окна функции яркости синусоидальной функцией двух переменных: 77(х,у) = = А$т[сохх + (оуу + <р). Значение локальных проекций градиента в центре окна (при х=у=0): с™
и.
/"-А
"а," со.
И-
Параметры А, <ох , ооу, ср, можно определить методом наименьших квадратов, который сводится к системе нелинейных уравнений, которую аналитически решить невозможно. Численные методы очень сложны в реализации, а в данном случае ненадежны в связи с неоднозначностью обратных тригонометрических функций. Поэтому данная задача решается без использования метода наименьших квадратов и формулируется таким образом, что допускает аналитическое решение. Величина А является константой, равной наибольшему значению функции яркости. Необходимо определить параметры сох, соу, ф. Применим метод сведения переопределенной задачи к такой, где число уравнений равно числу неизвестных. Для этого перейдем от используемого нами вида локальной маски к маскам нескольких видов. Маска первого типа соответствует наклону полосы в центре окна на 45° и используется при определенной конфигурации значений яркости, соответствующей условию (рис. 13):
[('/ «зН<, <0)М(', >'зН'/ >Ь)1 (14)
Выберем стандартную декартову систему координат, обеспечивающую наименьшее значение фазы (что определяется неоднозначностью арксинуса при выполнении условия (14)). Для данной задачи важен лишь знак частного сох / юу. Поэтому частоту ©х считаем всегда положительной, а юу с данной маской при выполнении условия (14) - отрицательной. Это всегда можно сделать выбором начала отсчета и начальной фазы. Таким образом, если условие (14) выполняется, то приходим к системе следующего вида:
Аятср Аып(-а)у -нр) =/2
А51П((йх+<р) =(3
решение которой^
Ж
Ф = агс51П — А
агс$т—~ ф Г А 41
Рис.13. Миска первого типа > 4
Если условие (14) не выполняется, воспользуемся маской второго типа (рис.14), соответствующей наклону полосы в центре окна на 135° (14). Частоты сох ,юу являются положительными и система уравнений будет иметь вид:
Абш ф -I! \Ахт(- со„ +ср) =/2 Ахт(-ах + ф,> =/3 решение которой:
И
Рис.14. Маска второго типа.
Ф = агс.чт—,
У А
. ь • ь
Ф - агс.ип — А , ®у = Ф -агсзт— А
Маска третьего типа соответствует наклону полосы в центре окна на 90° (рис.15), что соответствует выполнению следующего условия:
[(/, </>(/, <'з)Н('2 > > <з)}- (15)
В этом случае частоты: о~>х>0, <эу =0, и система уравнений будет иметь следующий вид:
I /1л'Шф =
Рис.15 Маска типа три
<яг = агсят'——--ф
Аят( ®х + у) = —,
А
решение которой:
■ Ч +*з Ф = агат-— Т 2А
2А
Маска четвертого типа соответствует горизонтальному расположению полос (рис.16), что соответствует выполнению следующего условия:
[('/ «3Н2 <ОМ('; ><зН2 >'<)]- (16)
Аналогично частоты:
сох=0, ©у>0, и система уравнений имеет вид:
{и+Ц)/2
Рис.16. Маска четвертого типа.
л • {1+*2 ЛШПф = —--
Азт(- со +(р) =
решение которой:
Ф = агс5т ■
2А
-> со,
Ф - агс5т
2А
, со =0
Во внешней маске W укладывается (М-1)х(Кг-1) локальных масок размером 2x2. Таким образом для определения частоты функции яркости, соответствующей центру глобального окна используются усредненные по окну значения модулей локальных
k I к I
частот со 'со'
М IN-11
к J\
J М 1N-1,
(17)
к,1\
При этом частотным компонентам сох и со^
присваивается тот знак, который преобладает во множестве локальных частот. На рисунках 17 - 19 приведено исследование точности оценивания поля направления градиентными методами в зависимости от интенсивности шума и пространственной частоты квазипериодической структуры.
Рис. 17. Зависимость погрешности оценивания поля направления от интенсивности шума (маска 7x7) для класса градиентных методов, использующих стандартные разностные схемы.
5«°
45 -------------
40 35 + 30 25
20
1
15 f 10 5 О
синусоидальная аппроксимация
< углы)
— раэност.сх. (проекции )
0 0,2% 0,5% 1% 2% 5% 10%
----------синус.аппрок. (проекции)
■- ........синус.аппрою(углы)
-разност.схемы, (крест- проекции)
---разност ,схемы.(крест- углы)
----разност .схемы.(2x2- проекции)
—о— разност.схемы.(2х2- углы)
Рис.18. Зависимость погрешности оценивания поля направления от интенсивности шума (маска 7x7) для различньис градиентных методов.
Синусоидальный метод является устойчивым к шуму при увеличении пространственной частоты квазипериодической структуры.
Методы, основанные на разностных схемах, использующие усреднение проекций обеспечивают более высокую помехоустойчивость при достаточно большой интенсивности шума, однако методы, использующие усреднение углов обладают тем преимуществом, что являются более точными при высоких пространственных частотах.
0.04ц 0.01п 0.1 Он 0.1 Зя 0.18л 03* 0.5*
а)
25
0.04* 0.07л 0.10л 0.13я 0.18л 03л б)
0.04л 0.07л 0.10л 0.13* 0.18* 0.3* 0.5л в)
г)
16 14 12 10 3 6
h о Í-
Ч-.:
0.04* 0.07* 0.10* 0.13* 0.18* 0.3« 0.5* 0.04* 0.07* 0.10* 0.13п 0.18* 0,3* 0.5*
д) е)
Рис.19. Зависимость погрешности по- -строения полей направлений от частоты 5*5 полос. Разностные схемы: а) крест-проек- 9x9 ции, б) крест-углы, в) 2х2-проекции, г) 2х2-углы. Синусоидальные методы: д) усреднение углов, е) усреднение проекций.
3.3. Дифференциальные методы
Дифференциальных методы основаны на том факте, что производная функции яркости по направлению, совпадающем)' с направлением полосы, имеет наименьшее по модулю значение среди производных по направлению в текущей точке: у = argmin\f'( фу|. При использовании дифференциального метода с маской 3x3 по отсчетам внутри окна вычисляем производные вдоль направлений, определяемых углами 0. 45, 90, 135° (1). При этом используются симметричные разностные схемы:
ft-+0 2h
При использовании маски размеров 5x5 вычисляются производных по восьми направлениям, определяемых углами: 0, 26=arctg(l/2), 45, 63=arctg2, 90, 116, 135, 153° |1|
Основным преимуществом дифференциальных методов является максимально возможное быстродействие среди всех представленных классов методов построения поля направления.
20
15 усреднение проекций
усреднение углов
—I--1---\---1—
0,5% 0.01 0.02 0.05 2x2 с усреднением углов крест с усреднением углов 2x2 с усреднением проекций крест с усреднением проекций
Весовую функцию определим соотношением: И' = (тах\Г(у)\-пип\Г(ч^)' тах Существенным недостатком дифференциальных методов является сильная чувствительность к шумам, что даст более низкое качество оценки по сравнению с дру гими методами. На рис.20 показана зависимость погрешности оценки поля направлений для ннтерферограмм от интенсивности шума.
45 ,0" 40 • Э5 I 30 . 25 • 20 •*' 15 ; 10 ♦ 5 Г
О V
3*3
1-1
5*5 4
— » —i— »--
0.24 0.5% 1% 2%
54 104
Рис.20. Зависимость погрешности построения поля направлении дифференциальны.ми методами от интенсивности шума для двух типов квашпериодических структур.
Основным ограничением методов является принципиальная невозможность уменьшения максимальной ошибки оценивания до значений меньших чем 90"п. где п • количество производных по направлению (п 4 для маски 3x3, п Л - 5x5).
3 4 11роекционно-<)исперсионные методы
Проскцнонно-дисперсионнын метод построения поля направлений основан №1 применении преобразовании Радона |9| хи вычисления "томографических" проекций функции яркости изображения 1<х.у> внутри скользящего окна вдоль направления, перпендикулярного направлению, задаваемому углом ц' (ряс 21)
Г 2
R М= i
Г 3
+ z сол( ц) + arctg — ),
V/' * : л/п(ц1 t arclg — ))ik
(18)
где предполагается, что центр окна рал мера 7x7 находятся в точке (0.0) Рассматривая функцию (18) при фиксированном шаченин параметра у и значениях аргумента 1*(-Т 2. Т 2/, можно измерить разброс шаченнй относительно среднего (дисперсию)
г :
г ;
la оценку значения поля отравлений принимается значение угла лоставляюшес максимум функции />{ф) ф = org max/>(ц/) Вссовую функцию комплексного поля направления определим из следующего соотношения
и = max Div)- mm lAv) max ¿K^)
, у у J у
Рис.21. Проекционно-дисперсионный метод.
Непосредственное использование преобразование Радона (18) для дискретного случая вызывает трудности, связанные с невозможностью вращения окна обработки Практическая реализация дисперсионного метода построения поля направлений сводится к вычислению функции D(vj/) для частного случая четырех (вертикального, горизонтального н диагональных) направлений на дискретном множестве отсчетов функции яркости 1(х,у) внутри окна
Произведем сканирование изображения 1(х,у) квадратной А/хА/ маской и ддя каждого из ее положений сформируем четыре последовательности средних значений для отсчетов, расположенных вдоль линий, параллельных одному из четырех направлений: 0, 45. 90. 135°
) = ~ J + U)M = U....2M - /.
k j,t
— = 1.2......
А/ ,-/
= 1,2,... M; (19)
А/
k ,,,
k = I.2....JM - i.
В результате получаем дискретный вырнант
к
(20)
к = 1.2......V/. vy = 0 или ф = 90'
I А/
JL /_ м{
.и-к
А/
, к = 12,...-1,
>
Ф = 45' или \|/ = 135*.
Основным недостатком такой реализации 111 проекционно-дисперсионного метода является кван-тованность результирующих значений поля направлений (0е: 4.5° 90°, 135°) Можно увеличить число фадаций до восьми, обобщив выше приведенные соотношения за счет рассмотрения четырех дополнительных направлений 26" агс!%(1 2). 63° /т(ц2, 116е, 153". Однако такой метод увеличивает трудоемкость и громоздкость алгоритмов Рассмотрим подход к увеличению точности за счет интерполяции (рис.22).
Diy)
IX : Ч 45
i :\ Щ.
о 2<ii аь \ :
Рис.22. Иллюстрации дисперсионного метода с интерполяцией.
Дисперсионный метод с интерполяцией. Допустим, что функция 1)(\\/) имеет вид синусоиды /j(v}/) = А ту(2(\|/ - »jy))-г Н Зная значения функции D0,D4S,D90.O,)s в точках \j/=0°, 45°.90о. 135°. можно оценить неизвестные значения параметров А,\\>,Н :
A cos 2 у + В - I) 0
Aco.s\--2yj + B = D4} D0+D4}+D„0+Dns=B.
откуда полу чаем
lg2y =( 04S - И)'( D„ -ti)
(21)
На отрезке [0°, 180°| данное уравнение имес! два решения, соответствующих минимальному и максимальному значению функции Щу)- А'1Я ного метода весовую функцию поля направления определим следующим образом и 2А (А Н>
I
переборной
ИиТерПОГЯ ц ИОННЫЙ
0 2% 0 5% 14
2%
Рис. 23. Зависимость погрешности построении полей направ.1ений проекционно-дисперсионными методами от интенсивности шума д.ля Овух типов квазипериодических структур
Исследование точности проекпионно-дисперсионных методов (рис.23-25) показало, что интерполяционный метод в 2 - 12 раз увеличивает точность
оценивания поля направления по сравнению с четы-рс.хградационным переборным методом (рис 23) Интерполяционный метод проитрывает переборно-м> только для больших значений пространа венной частоты и большой интенсивности шума
3x3
с
(10ре<м>рныи
J в
О
ИИТОр'ЮЛИЦНОНКМЙ
0.2% 0,54 1% перЫ>ор 3*3 поре<*>р 13*13 Интерпол 91(9
2% 5% nepebop 9г9 Интерпол 3x3 имгврпог. 13« 13
h
Рис.24. Зависимость погрешности оценки полей
н оправлении проекц ионно-ди сп ере и он н и методами от интенсивности шуми От рахшч-ных мисок.
16 ■ 14 5
4 5 [ 14 •
•I3.5 • 13 t f2 Ь J "
12 i , 115 1 • 1 •
11 I..... 0......
' МГц U.ЛI* ai«» 0.1» 0Л* ч/W* МГц i Um 0 I V« 1H* 0.1« U.Sn
<V
(ч uimitрптнционный
S«', mi
о/ переборный
Рис. 2 5. Зависимость погрешности "уценки поля направлений проекционными методами от пространственной частоты квашпериодической структуры
3 5 ( пектрсиьные метопы
Спектральный метод основан на модели изображения с локально-периодической функцией яркости П соответствии с згой моделью производится Спектральный анализ функции яркости изображения
''*.}') » пределах скользящего окна 1У(х(„у^ (*„. О ¡(.Г, у)-\х-х0\<7/2.ц - <. /'/2}, где
* <. }„/ - координаты центра окна "/ч/-размеры окна I Рассмотрим коэффициенты ряда <1>\рьс функции М зкости в пределах окна
Г ^я/
С(т.п)- \\1{х. V )едг/а (тх + пу) |<1х<Ь
Т
' определим координаты максимального по модулю о оффнциента
{т,,) - игу тих\( (т. п)
I Л», л I'<0 01
Индексы (т0 п^ (рис 26) определяю! превали-Уюшис в данном фрагменте изображения простран-п кгнные частоты м/хсуо/, о)/ХаУо> с точностью до н; «ка и позволяют оценить значения поля нагтравле-И1 й и поля пространственных частот в точке (х^уы
5
г
3
- 1 1 2 'у
— 1 __
-7 - 6- 5 - 4 - . 1 3 -2- 1 1 Ч-! \ ] < н
т
С(ш,п)|
Рис.2б. Иллюстрация спектрального метода. ГПг
по
(22)
¿>(хО'Уо) = у^'по +
В отличие от метода локальной квадратичной аппроксимации здесь размеры окна ограничены снизу: при каждом положении на плоскости изображения оно должно пересекаться несколькими "интерференционными" полосами. Чем больше полос попадает в окно и чем меньше искривленность полос внутри окна, тем более точными будут оценки (22). Практическая реализация спектрального метода расчета поля направлений основана использовании рекурсивных алгоритмов быстрого преобразования Фурье. Но, несмотря на это, данный метод чрезвычайно трудоемок с вычислительной точки зрения. Рассмотрим метод основной локальной частоты, основанный на непосредственном "геометрическом" способе оценивания пространственных частот.
Метод основной локальной частоты. Этот метод также как и представленный выше спектральный метод, основан на поиске главной гармоники. Отличие заключается в процедуре поиска. Компоненты пространственной частоты для малой области изображения квазипериодической структуры равны о)х=2л7(1 а а)у-2л/с1>, где с1х, (1У - периоды функции яркости 1(х, у) в вертикальной и горизонтальной развертке, задача сводится к определению периода сечений функции яркости по горизонтали и вертикали._
На рис.27 отмечены по вертикали и горизонтали периоды функции яркости, равные удвоенному расстоянию между соседними экстремумами (минимума-Рис.27. Иллю- ми и максимумами). В силу страция метода наличия шумов на рассмат-основной локальной риваемых сечениях функ-частоты. дни яркости присутствуют
побочные минимумы и максимумы. Поэтому необходима процедура их отсечения. Отбросим минимумы, лежащие выше порогового значения Птпт и
максимумы - ниже порогового значения Птах (на рис.28 они соответственно обозначены - 0 и *). В результате на ярких участках группируются максимумы. не разделенные минимумами, ни темных -минимумы.
X
Рис. 28. Удаление ложных экстремумов.
Поскольку экстремумы разных типов должны чередоваться, то воспользуемся процедурой прореживания. При этом оставим один первый максимум (минимум), остальные до ближайшего минимума (максимума) отсечем (рис.30).
Для определения пространственной частоты просканируем исходное изображение. Размер маски должен составлять приблизительно два - три периода функции яркости. При увеличении размера маски нарушается условие однородности квазипериодической структуры, при уменьшении - возрастает влияние краевых эффектов, х
Рис.29. Удаление группы экстремумов.
Для каждой текущей точки, соответствующей центру окна, определим значения пространственных частот сох и од. В соответствии с размером маски сформируем периодически обновляемый буфер из строк (столбцов) изображения. Выберем в маске первую, среднюю и последнюю строки (столбцы). В каждой строке (столбце) определим количество и номера минимумов и максимумов (отсчетов с наименьшими и наибольшими значениями яркости).
X
Рис.30. Удаление соседних экстремумов.
Определим среднее расстояние между пиками (максимумами) в строке (причем расстояние равно нулю, если пиков в строке (столбце) нет, и размеру маски, если пик один), и определим соответствующее значение частоты. Результирующую частоту
определим как среднее значение найденных частот по выбранным строкам (столбцам) сканируемой маски. Оценка значения поля направления равна отношению найденных компоненту и ц,. В отличие от спектрального или градиентных методов в данном случае возникает неоднозначность, связанная с невозможностью определения знаков компонент пространственной частоты. В данном случае мы не сможем отличить направления а и (180 - а). Для определения знака частоты предлагается анализировать сдвиг экстремумов функции яркости в соседних линиях развертки изображениях.
На рис.31 в приведено исследование точности оценивания поля направления в зависимости от пространственной частоты квазипериодической структуры для различных размеров окна сканирования. На рис.31а.б приведена точность оценивания поля направления при различной интенсивности шума для двух сечений данного графика, соответствующих низкой и высокой частоте структуры. При увеличении пространственной частоты повышается точность оценивания поля направления. При наличии шу мов погрешность достаточно высока, но не сильно увеличивается при увеличении интенсивности шума (это объясняется природой накладываемого шума (высокочастотной) и тем, что спектральные методы хорошо работают на высокочастотных квазипериодических структурах (размер окна обработки ограничен снизу и должен содержать как минимум несколько "интерференционных" полос).
50 40
зо 20 10 О
Г
О 0.2% 0,5% 0,01 0,02 0.06 0,1
О 0,2% 0,5% 1% 2% 5% 10%
0.04л 0.07я 0.10л 0.13л 0.18л 0.3 л 0.5 л Пространственная частота ы [ рад / пиксел]
50 30 ,20 , 15 , II . 6 ,
пиксел / период d
маска 5x5---маска 7x7
-маска 9x9 ----маска 11x11
—о-маска 13x13
4
Рис.31. Зависимость погрешности построения полей направлений методом "основной локальной частоты" а),б) от интенсивности шума для различных масок и различных частот: а) о)=О.Зп. б) а>=0.08л; в) от пространственной частоты.
4. Экспериментальные исследования методов построения поля направлений На рис.32 приведено исследование точности оценивания поля направления различными методами в зависимости от пространственной частоты, а так же от периода анализируемой квазипериодической структуры при обработке окном 9x9 (дифференциальные методы окном 5x5) и в зависимости от интенсивности шу ма для высокочастотной интерфе-рограммы рис.32в и низкочастотной рис.326, соответствующих сечениям графика рисунка 32а.
0.04я
50
I—
0.07л 0. Юл 0.13л 0.18я I О.Зя
|
Пространственная частота (рад I пиксел]
О.Зя
30
20
15
11
пиксел I период
-4~
а)
Г"
60 50 40
30
20 <Г 10
0
'-Ь
h'
0.5%
2%
10% 30%
5С%
0,5% 2% 10% дифференциальные градиентные (2x2) спектральные
30%
50% проекционные аппроксимационные
Рис.32. Зависимость погрешности оценивания поля направления разными методами (маска 9x9): а) от пространственной частоты; б), в) от интенсивности шума: б) при со~0.07л рад/пиксел, изображение 256x256, в) при (о-О.Зк, размер изображения 128x128 пикселей.
Максимальную теоретическую точность (при отсутствии шумов) обеспечивают градиентные и аппроксимационные методы построения полей на-
правлений. При уровне шу ма 0.2% и выше преимуществом по точности оценивания поля направления (в диапазоне низких частот) обладают проекционно-дисперсионные методы. Для высокочастотных квазипериодических структур преимущество данных методов начинает проявляться при уровне шума 0.5% и выше, при этом в случае наличия большого шума методы являются достаточно помехоустойчивыми. и ошибка оценивания практически не зависит от уровня шума. Для высокочастотных квазипериодических структу р наилучшую точность оценивания поля направления обеспечивают спектральные методы. Влияние шума начинает сказываться на работе спектральных методов в диапазоне высоких частот при отношении шум/сигнал 5%. Для дифференциальных методов наблюдается резкое возрастание ошибки оценивания при пространственных частотах квазипериодичсской структуры выше О.Зя рад/пиксель, для данных методов молено также отметить сильную зависимость от шумов. Аппроксимационные методы обеспечивают хорошую точность оценивания поля направления, если окно обработки захватывает не более одной полосы квазипериодической структуры. Наибольшим быстродействием обладают аппроксимационные, дифференциальные методы, метод основной локальной частоты.
10 11 12 13
I-дисперсионный переборный. 2-дисперсион. интерпоч.; 3-дифференциальный; 4-градиентный (синусоидальн. Ан-прокс.) с усреднением углов: 5-градиентный (синусоидальн. аппрокс.) с усреднением проекций;6-градиешпный (разностные схемы) с усреднением проекций (2x2); 7-градиентный (разностные схемы с усреднением углов (2x2); ^-градиентный (разностные схемы с усреднением проекций (крест); 9-градиентный (разностные схемы) с усреднением углов (крест); 10 - цилиндрическая аппроксимаций; Л-аппроксимация плоскостью; /2-аппроксимация квадратичной поверхностью; !3-комбинированная аппроксимация; 14-спектральный (основной локальной частоты<).
Рис.33. Диаграмма затрат времени на обработку квашпериодической структуры размером 256x256 различными методами построения полей направлений дня различных масок.
5. Примеры практических задач, решаемых с использованием поля направлений
5.1. Анализ дактилограмм.
Дактилоскопические изображения широко используются в криминалистике для идентификации
личности правонарушителя, а также в некоторых типах устройств для ограничения доступа персонала. Представленная в работах [1, 10] методика кодирования изображений отпечатков пальцев основана на обнаружении глобальных и локальных особенностей и расчете их геометрических характеристик (расстояния и углы между особыми точками). К глобальных особенностям относятся тип узора ("петля", "дельта", "спираль'" и т.п.). Локальные особенности - это точки разрыва, ветвления и слияния линий папиллярного узора. Метод поля направлений может быть применен для обнаружения и определения координат глобальных особенностей, для оценки их геометрических характеристик и для построения процедур поиска локальных признаков дактилограмм. Для вычисления поля направлений в данном случае возможно использование спектрального или дисперсионного методов (метод локальной квадратичной аппроксимации не применен из-за высокой густоты полос, сравнимой с частотой дискретизации дактилографического изображения). В данных исследованиях был использован проекци-онно-дисперсионный метод, дающий четыре градации направления, а также градиентный метод с использованием синусоидальной аппроксимации. На рис.34 приведены исходное изображение и результаты расчета и фильтрации поля направлений.
а) б) в)
Рис.34. Поле направлений дактилокарты
а) функция яркости исходного изображения;
б) четырехградационпое поле направлений; в) результаты линейного сглаживания. Распознавание типа папиллярного узора основано на анализе конфигураций поля направлении на границе скользящего окна. При этом используется четырехгра-дационное поле направления. Для приведения непрерывных значений к четырем градациям воспользуемся номерами секторов (0, 1, 2, 3), определяемые следующим диапазоном углов: [0°-22°, 158°-180°]; [23°-67°]; [68°-1120]; [113°-157°]. Центры приведенных диапазонов соответствуют углам 0°, 45°, 90°, 135°. Окрестность каждой точки дактило-граммы характеризуется
е.>ьты
т ■ш
■ ^я
V Я шРнВёр
. -
Спираль Рис.35. Глобальные особенности
определенным порядком смены значений поля направлений. Для глобальных особенностей этот порядок отличается от порядка произвольной точки отпечатка (рис.35,36).
Оценивание геометрических параметров глобальных особенностей
производится на основе геометрии секторов поля направлений. В качестве геометрических параметров глобальных особенностей используются; угловые характеристики дельт, наклон оси спирали, характеристики закручивания спирали: степень отклонения от вертикали, размах секторов поля направлений слева и справа от центра, признак несимметричности расположения дельт (рис.37).
Применение поля направлений позволяет существенно упростить алгоритм распознавания и анализа глобальных и локальных особенностей.
5.2. Анализ кристаллограмм
Конфигурация Тип точки
1-0-3-2-1-0-3-2 0-3-2-1-0-3-2-1 3-2-1-0-3-2-1-0 2-1-0-3-2-1-0-3 Центр спирали
0-1-2-3/1-2-3-0 2-3-0-1/3-0-1-2 Дельта
1-0-3-2/-3-2-1 3-2-1-0/-1-0-3 Петля
2-1-0-1 дв,левая петля
3-2-3-0 дв. прав, петля
1-0-3-0 Шатер
Рис.36
Рис. 3 7. Геометри ча ние характерасшина глобальных особенностей.
Для поиска локальных особенностей (точек разрыва, ветвления и слияния папиллярных линий) используются традиционные методы получения бинарных препаратов изображений и дальнейшая их логическая обработка с использованием поля направлений (рис.38). Классификации подвергались точки, которые лежат на качественных областях, определяемых критериями качества, полученных на основе анализа изображения скелетизированного папиллярного рисунка и весовой функции поля направления. __
Рис.38. А пали >• проклассифицированных . ¡окольных особенностей папиллярного piuynur
В настоящее время в офтальмологии существуют методы диагностики глазных болезней по состоянию кристаллограммы слезной жидкости пациента. Основными признаками патологичности являются большое количество центров кристаллизации, большой разброс направлений и высокая плотность кристаллов [2, 3] (рис.39).
Криаиа i икрим.чм
Лрасти сш.'рчммы пита и/.ч¡ей Рис.ЗУ. Характерные к юбра женин криапаио. рччм слезной жидкости <■ порче и /¡/-о пат.'логичегких излнчи пнях ор.чгнах '.рения.
Данными признаками также обладают кристаллограммы плазмы крови [11], используемые аналогично для диагностики ряда заболеваний. На рис.40 приведены образцы кристаллограмм плазмы крови и соответствующие поля направлений с основными характеристиками. Для анализа степени разброса направлений кристаллизации используется графический препарат поля направлений, на основе которого определяется коэффициент однонаправленности:
Код*Що\
дцt(x,y) 1 \ i
дх 1 Oy У
В качестве признаков классификации были выбраны также следующие характеристики кристаллограмм, построенные с использованием поля направлений: относительная площадь областей с выраженным направлением Кт /Я. где .9,,- площадь участков, выделенных на основе анализа весовой функции поля направлений. $ - площадь кристаллограммы; густота лучей, определяемая на областях четких линий Д.: Кг = т-^-г || («~ + г л ^ \lxdy,
|Д| о,
прозрачность кристалла. Кпр = ¿7У ) /, , где
1=1тп \lHx.yjdxdy, К ={1тах +/тт)/2. \ич\ (х.у/ьР,
норма
пито. югин
■ Г\ ' " . ;
Исходное Сглаженное Весовая функция
изображение поле направлении поля направлении
Рис.40. Анализ кристаллов плазмы крови
К опт) ■рн ы ù преп арат поля направлений
5. Заключение Эффективным методом анализа изображений, характеризуемых наличием квазипериодических структур, является метод поля направлений. Он позволяет сократить структурную избыточность таких изображений, упрощает обработку и анализ, повышает качество распознавания. Метод поля направлений дает возможность выявить существенные признаки изображений, позволяющие провести эффективный анализ. В представленной работе предложены численные методы оценивания поля направления и сравнительное исследование их точности и быстродействия. Все методы разбиты на пять классов в зависимости от используемой модели представления поля направлений и подхода к его вычислению: методы параметрической аппроксимации. методы локальных градиентов, дифференциальные методы, проекционно-дисперсионные методы, спектральные методы. Для оценивания точности методов в качестве критерия использовалась приведенная угловая среднеквадратичная погрешность, рассчитанная с использованием комплексного представления поля направлений. Разработанная классификационная схема методов построения поля направлений позволяет осуществить выбор наилучших методов для конкретных задач обработки изображений.
Список литературы
1. Ильясова Н. Ю., Устинов А. В , Храмов А. Г. Методы анализа дактилоскопических изображений на основе поля направлений// Научное приборостроение, т.З, с.89-101. 1993, Санкт-Петербург.
2. Ильясова Н. Ю., Устинов А. В. Компьютерный анализ изображения кристаллов слезы// Тезисы доклада на 2-ой международной конференции "Распознавание^", Курск, 1995, с.248-250.
3. ДворяноваТ.П., Ильясова Н.Ю, Устинов А.В., Храмов А.Г. Компьютерная система анализа диагностических кристаллограмм // Компьютерная оптика, вып. 16, с.90-96, 1996.
4. Ченцова О. Б., Прокофьева Г. Л. Кристаллографический метод обследования при некоторых заболеваниях глаз // Методические рекомендации. М.: 1988г.
5. Soifer V. A, Kotlyar V V., Khonina S. N„ and Kliramov A. G., The Method of the Directional Field in the Interpretation and Recognition of Images with Structure Redundancy, Pattern Recognition and Image Analysis, v.6, No.4, pp.710-724 (19%)
6. Glumov N.I. Krainukov N.I. Sergeyev V.V. Kliramov A G. The Fast Algorithm of Image Approximation in a Sliding Window// Pattern Recognition and Image Analysis, 1991, N 4, p.424-426
7. Бронштейн И. H., Семендяев К.А. Справочник по математике, 1956 г.
8. Сергеев В. В.. Фролова Л. Г. Разработка и применение алгоритма цилиндрической полиномиальной аппроксимации изображения в скользящем окне// Автометрия, 1, 1996, 22-30с.
9. Хермен Г. Восстановление изображений по проекциям: Основы реконструктивной томографии. - М.: Мир, 1983 - с352.
10. "Руководство по системе "Фрагмент"; УВД Куйбышевского горисполкома, 1976
11. Hozman J., Kubinec R., Traka J., Varenka J. Automatic computer evaluation of blood serum crysiallograms//' "Biomedical Image Processing Applications", in Biomedical Engineering & Biotechnology, (A.Strejc, ed.), (Praha), Publishing House of the Czech Technical University, 1995.