CC BY
34
РЕКУРСИВНЫЙ МЕТОД ПОСТРОЕНИЯ ПОЛЯ ЧАСТОТ КВАЗИПЕРИОДИЧЕСКИХ ИЗОБРАЖЕНИЙ
А.В. Устинов Институт систем обработки изображений РАН
1. Введение
Среди изображений, возникающих в реальных задачах обработки, существенную долю составляют квазипериодические изображения. Примерами квазипериодических изображений являются интерфе-рограмма, кристаллограмма, дактилоскопические изображения. Важными параметрами таких изображений являются преимущественное направление полос и их частота в каждой точке изображения. Полем частот будем называть поле мгновенных значений частоты полос в данной точке. Ниже описывается эффективный алгоритм его построения, основанный на предварительном нахождении точек экстремума функции яркости.
2. Метод основной локальной частоты
Во многих задачах обработки квазипериодических изображений (например, восстановление фазовой функции по интерферограмме) можно считать, что функция яркости в локальной области имеет вид ^,у) = A + Bsm(®xx + ®yy + р). При этом возникает необходимость оценки частот ах и а у. Наиболее логически простым, но вычислительно очень неэффективным методом является нахождение преобразования Фурье от функции яркости в пределах скользящего окна и определение номеров максимального по модулю отсчёта спектра (кроме нулевого). Предлагается другой подход к поиску главной гармоники - метод основной локальной частоты. Его преимущество именно в том и состоит, что нам не надо искать все гармоники, чтобы выделить главную. Оба метода применительно к построению поля направлений были кратко описаны в [1], хотя точность алгоритмов вычисления частоты в [1] не исследовалась.
Предлагаемый метод использует следующую процедуру поиска главной гармоники. Компоненты пространственной частоты для малой области изображения квазипериодической структуры равны ах=2п/й х, ау=2п/й у , где й х , й - периоды функции яркости 1(х, у) в вертикальном и горизонтальном направлении. Поэтому задача сводится к определению периода сечений функции яркости по горизонтали и вертикали.
На рис.1 по вертикали и горизонтали отмечены периоды функции яркости, равные удвоенному расстоянию между соседними экстремумами (минимумами и максимумами). В силу наличия шумов на рассматриваемых сечениях функции яркости присутствуют побочные минимумы и максимумы. Поэтому необходима процедура их отсеивания.
Отбросим минимумы, лежащие выше порогового значения П^ и максимумы - ниже порогового значения П^ (на рис.2 они соответственно обозначены
- 0 и *). Пороговые значения определяются для каждой строки и столбца матрицы отсчётов функции яркости и вычисляются по формулам Пmm.=A-B/2, П^ =A+B/2. В результате на светлых участках группируются максимумы, не разделенные минимумами, на темных - минимумы (рис.3). Если значения A и B в формуле для 1(х,у) существенно изменяются по полю изображения, то при программной реализации предусмотрена возможность адаптивного оценивания значений A и B.
Рис.2. Удаление ложных экстремумов
П„
ГГ.
Рис.3. Результат удаления ложных экстремумов
Рис.1. Метод основной локальной частоты
Рис.4. Окончательное расположение экстремумов
Поскольку экстремумы разных типов должны чередоваться, то воспользуемся процедурой прореживания. При этом оставим один первый максимум (минимум), остальные до ближайшего минимума (максимума) отсечем (рис.4). Такая процедура поиска экстремумов производится для каждой строки (горизонтальное сечение) и столбца (вертикальное сечение). В результате формируются два вспомогательных массива с размерами исходного изображения, каждое значение которых равно 0, если в данной точке функция яркости не имеет экстремума, 1, если есть минимум, и 2, если максимум. В первом массиве находятся экстремумы горизонтальных сечений, во втором - вертикальных.
Отметим близость описанной процедуры алгоритму выделения центров полос Рамеша - Прамода [2,3]. А именно, сформированные два вспомогательных массива, по сути эквивалентны частичным скелетам в направлении 0° и 90° в алгоритме Рамеша -Прамода. Отсюда вытекает следующий вывод. Если при обработке квазипериодического изображения уже проведено выделение центров полос, причём выделены и минимумы, и максимумы яркости, а применённый алгоритм не использует аппроксимации функции яркости в пределах окна, то процедуру поиска экстремумов можно пропустить, а в качестве вспомогательных массивов использовать изображение центров полос.
Для определения пространственной частоты про-сканируем исходное изображение. Размер маски желательно использовать не больше двух -трёх периодов функции яркости, так как при увеличении размера маски нарушается условие однородности квазипериодической структуры. Для каждой текущей точки, соответствующей центру окна, определим значения пространственных частот сох и соу. Для этого во вспомогательных массивах анализируем строки, соответствующие первой, средней и последней строке (столбцу) маски. Таким образом, определим количество и положение минимумов и максимумов.
Определим среднее расстояние между экстремумами в пределах маски и вычислим соответствующее значение частоты. Расстояние можно найти, если в пределах маски имеется хотя бы два пика, что, как правило, не выполняется на низкочастотном изображении. В таких случаях поступаем следующим образом: ищем ближайшие пики за пределами маски слева(сверху) и справа(снизу). Если внутри маски пиков нет, то используем расстояние между найденными пиками, если внутри маски пик один, то используется расстояние между пиком в маске и ближайшим из двух внешних пиков. Рассмотрим вырожденные случаи: при отсутствии пиков на линии частоту полагаем равной нулю; при одном пике на линии расстояние между пиками принимаем равным 3М4, где N - длина линии (такой выбор можно оправдать). Результирующую частоту определим как среднее значение найденных частот по выбранным трём строкам (столбцам) сканируемой маски. Наличие вспомогательных массивов обеспечивает рекурсивное определение положения пиков при
сдвиге маски на одну позицию, что существенно уменьшает число арифметических операций. Оценка
значения поля частот равна с =
-с.
3. Исследование точности алгоритма
Точность алгоритма исследовалась на тестовых изображениях размером 256х256 с постоянной частотой: с прямыми полосами (рис.5) и круговыми полосами (рис.6). В качестве меры точности использовалась относительная среднеквадратичная погреш-
ность в = -
1
М N
МИ
-с)
1=1 }=1
с
где М^ - размеры
изображения.
Рис.5. Угол наклона полос 15° т=0,6
Рис.6. с=0,6 Поля частот (фильтрованные) для изображений с круглыми полосами и постоянной частотой приведены на рисунках 7 и 8, на которых яркость пропорциональна частоте, а чёрная кайма означает, что на краях изображения вычисления не производятся.
Рис.7. а=0,6; размер маски 21
Рис.8. а=1,26; размер маски 21
Таблица 1. Зависимость погрешности от размера маски. Частота а=1.26, угол наклона полос а=15°
Размер маски 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39
в,% 21,18 11,15 8,39 4,45 2,50 2,64 1,47 1,26 1,64 1,18 0,95 0,94 0,65 0,56 0,96 0,44 0,52 0,55 0,77
Таблица 2. Зависимость погрешности от частоты при четырёх размерах маски. Угол наклона полос а=15
Частота 0,03 0,09 0,15 0,30 0,45 0,60 0,75 0,90 1,20 1,50 2,25 3,00
е,%; м. 11 16,87 19,99 11,87 2,21 2,56 4,17 3,50 2,96 2,67 2,59 8,26 15,80
е,%; м. 15 17,59 17,05 7,82 3,04 2,95 2,25 2,29 2,10 2,58 1,31 7,17 18,61
е,%; м. 21 15,26 12,97 2,12 1,68 1,48 1,66 1,22 1,09 1,00 1,14 7,66 25,14
е,%; м. 29 15,55 6,05 1,94 1,40 0,79 0,95 0,64 1,13 0,76 0,39 5,96 26,54
Таблица 3. Зависимость погрешности от размера маски. Частота а=1.26
Размер маски 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37
е,% 20,65 14,41 10,39 7,70 6,93 6,47 6,15 6,00 5,88 5,81 5,72 5,74 5,60 5,59 5,53 5,48 5,54 5,46
Таблица 4. Зависимость погрешности от частоты при четырёх размерах маски
Частота 0,03 0,09 0,15 0,30 0,45 0,60 0,75 0,90 1,20 1,50 2,25 3,00
е,%; м. 11 35,17 34,49 14,48 5,35 3,94 4,71 4,51 3,29 3,01 3,02 6,10 15,00
е,%; м. 15 34,97 31,25 10,83 4,66 4,07 3,58 3,17 2,80 2,22 2,20 5,98 17,01
е,%; м. 21 34,69 25,67 7,68 4,83 3,22 2,83 2,55 2,37 2,15 2,12 6,10 18,82
е,%; м. 29 34,45 18,30 6,37 4,20 3,31 2,73 2,62 2,65 2,50 2,49 6,04 20,55
Таблица 5. Зависимость погрешности от радиуса кривизны линии для двух масок. Частота а=1.26.
Радиус кривизны 3 5 7 9 12 15 20 25 30 40 50 60 80 100
е,%; м. 15 21,05 10,36 8,08 8,04 9,62 8,72 8,19 7,64 6,82 6,33 6,31 6,22 5,92 5,84
е,%; м. 25 20,26 21,85 19,84 13,68 6,22 6,13 7,62 6,79 6,26 6,25 5,96 5,84 5,75 5,60
Результаты, показывающие зависимость точности от размера маски и частоты приведены в таблицах 1,2 (прямые полосы) и 3,4 (круглые полосы). Результаты позволяют заключить, что наибольшая точность достигается, когда период равен 4-7 пикселам, а в маске укладывается 2-4 периода. Сравнение первой и второй пары таблиц показывает, что при прочих равных условиях погрешность на изображениях с круглыми полосами выше, чем с прямыми. Причина этого в следующем: хотя частота у круглых полос всюду одинакова, в центре изображения велика кривизна линий, что приводит к противоречию с предположением о приблизительной локальной прямолинейности полос. Зависимость погрешности от радиуса кривизны приведена в таблице 5 (погрешность считалась в точках, лежащих на окружности данного радиуса). Видно, что удовлетво-
рительные результаты получаются только, если радиус кривизны больше 1/2 -г- 2/3 размера маски. Коэффициент несколько колеблется в зависимости от соотношения частоты и размера маски в силу эффекта биения, влияние которого особенно чётко видно при величине периода, очень близкой к целому числу (см. таблицы 3 и 4). Из-за этого эффекта погрешность при частоте ю=1,26 более чем вдвое больше погрешности при частоте ю=1,20. Эффект биения возникает в связи с тем, что количество пиков при движении маски изменяется даже при постоянной частоте. Отметим, что значительное увеличение погрешности при малом радиусе кривизны возникает и во всех методах построения поля направлений, описанных в [1]. Свободен от этого недостатка оптический метод построения поля частот [4].
Угол наклона 0° 15° 30° 45° 60° 75° 90° 105° 120° 135° 150° 165°
в,% 2,33 1,19 1,58 1,23 1,58 1,19 2,33 1,19 1,58 1,22 1,58 1,19
45°
135°
225°
315°
'А,. „Л, ■А
0,04 О
i --0,04
0° ' 90° ' 180° ' 270° ' 360°
Рис.9. Профиль частоты вдоль окружности R=60,
точная частота со=0,6
0°
90°
270°
1360°
-г* 4 j- 4-
V J У—
0,04 О
0,04
45° ' 135° ' 225° ' 315° ' Рис.10. Профиль частоты вдоль окружности R=25,
точная частота со=1,2 45° 135° 225° 315°
- 0,04
- О
i ~-0,04
0° ' 90° ' 180° ' 270° ' 360°
Рис.11. Профиль частоты вдоль окружности R=60,
точная частота со=1,2
При исследовании обнаружился также несколько неожиданный момент - зависимость погрешности от угла наклона линии. В таблице 6 даны значения погрешности для прямых линий, а на рис. 9-11 показан профиль частоты для круговых линий, за нуль взято истинное значение частоты. Причём в случае кругов эта зависимость различна в области высокой кривизны и низкой кривизны, где она похожа на зависимость для прямых линий. В случае прямых линий или низкой кривизны круговых линий можно дать следующее объяснение. Можно было ожидать, что точность будет выше (но не настолько, как оказалось на самом деле), когда линии параллельны сторонам маски, но при таком расположении сильнее проявляются биения, так как появление/исчезновение пиков происходит сразу во всех трёх линиях. Поэтому наибольшая точность достигается при небольшом отклонении направления линий от сторон маски. Кроме того, измерения показали (см. также рис. 9,11), что измеренная частота в среднем больше истинной.
Таблица 7. Зависимость погрешности от уровня шума. Угол наклона полос 15 ° размер маски 15, а=1,26 (Первая строка - при гауссовском шуме, вторая - при равномерном шуме).
h 0 0,001 0,002 0,005 0,01 0,02 0,05 0,10 0,15 0,20 0,27 0,35 0,50 0,65 0,80
e,% 1,47 1,47 1,47 1,47 1,47 1,52 1,79 2,51 3,30 4,62 7,13 9,67 14,93 22,57 33,14
e,% 1,47 1,47 1,47 1,47 1,47 1,52 1,80 2,55 3,09 3,88 7,65 10,39 14,17 21,68 35,89
Таблица 8. Зависимость погрешности от уровня шума. Угол наклона полос 15 ° размер маски 19, а>=0,63. Шум гауссовский..
h 0 0,005 0,02 0,05 0,10 0,15 0,20 0,27 0,35 0,55 0,80
e,% 1,03 1,11 1,51 2,61 3,71 4,87 6,08 11,27 24,99 81,11 129,4
Зависимость погрешности от уровня шума исследовалась на изображениях с прямыми линиями, на которые накладывался аддитивный шум. Результаты приведены в таблицах 7 и 8, в первой строке которых приведено отношение шум/сигнал (по амплитуде). Видно, что качество вычисления частоты практически не зависит от распределения шума. Сравнивая таблицы 7 и 8, можно увидеть, что при более низкой частоте с увеличением дисперсии шума увеличение погрешности происходит быстрее. Причина в том, что при низкой частоте (малой крутизне функции яркости) больше вероятность того, что шумовой выброс приведёт к появлению ложного экстремума функции яркости, который не будет отсеян.
Сравним теперь наш метод и метод преобразования Фурье по скорости. Пусть изображение имеет размеры МхЫ, а маска КхЬ. В нашем методе при любых условиях число операций примерно равно ММ(3Ь+К+18), из них 12ММ с плавающей точкой. На низкочастотном изображении добавляются ещё 8
целочисленных операций для каждой точки, где требуется поиск пиков за пределами маски. Если такой поиск производится в каждой точке, то число добавочных операций равно 8МК Таким образом, имеем оценку сверху для общего числа операций: МЫ(3Ь+К+26). В методе преобразования Фурье (при нерекурсивной реализации) число операций (все с плавающей точкой) примерно равно ММ[р1(КЬ)1оя(КЬ)+(Р2+4)КЬ], где О: и 02 - константы, определяющие число операций для выполнения одного двумерного преобразования Фурье при наличии быстрого алгоритма: 01(КЬ)^(КЬ)+ 02(КЬ). Эти константы зависят от размеров К и Ь. При рекурсивной реализации число операций примерно равно ММ(16КЬ). Очевидно, что применение нашего метода требует значительно меньших вычислительных затрат. Рассмотрим благоприятный для метода Фурье случай, когда 01=0, 02«6 (то есть нерекурсивная реализация лучше, но такие и 02 бывают лишь при строго определённых К и Ь), а изображение низкочастотное. Пусть
К=Ь, тогда даже при К=3 имеем выигрыш в 2,37 раза, а при К=15 выигрыш в 26,16 раза. Правда метод преобразования Фурье может дать более точный результат при высоком уровне шума в силу его фильтрующего свойства.
4. Заключение Поле частот является одним из удобных средств описания квазипериодических изображений. В представленной работе предложен метод построения поля частот, проведено исследование его точности. Показано, что по быстродействию данный метод существенно эффективнее использования преобразования Фурье. Алгоритм нахождения точек экстремумов функции яркости применим также и для выделения центров полос. Поэтому можно рассчитывать на широкое использование данного метода.
Литература
1. Ильясова Н.Ю., Устинов А.В., Храмов А.Г. Численные методы и алгоритмы построения полей направлений квазипериодических структур // Компьютерная оптика, вып. 18, с.150-164, 1998.
2. K.Ramesh and B.R.Pramod. Digital image processing of fringe patterns in photomechanics // Optical Engineering, vol. 31(7), p.1487-1498, 1992.
3. K.Ramesh, R.K.Singh. Comparative performance evaluation of various fringe thinning algorithms in photomechanics // Electronic Imaging, vol. 4(1), p.71-83, 1995.
4. V.A.Soifer, V.V.Kotlyar, S.N.Khonina, A.G.Khra-mov, and R.V.Skidanov. Image Recognition Using a Directions Field Technique // Proceedings of SPIE, vol.3346, p.238-258, 1998.
