Выделение центров полос на интерферограмме Текст научной статьи по специальности «Математика»

ОБРАБОТКА ИЗОБРАЖЕНИЙ

Н. U. Крайнюков, А.Г. Храмов ВЫДЕЛЕНИЕ ЦЕНТРОВ ПОЛОС НА ИНТЕРФЕРОГРАММЕ 1. ВВЕДЕНИЕ

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

Если интерферограмма получена методом двух экспозиций, то при отсутствии шумов интенсивность интерференционной картины выражается формулой [1-5]:

1(х,у) = а(х, у) + Ь(х, у) • (1 + cos(5 (х, у))), где 1(х,у) - интенсивность в точке с координатами (х, у), а(х. У) ~ интенсивность фона, Ь(х, у) - контрастность (Ь(х, у) > 0), 5(х, у) - фазовая функция.

Считают, что функции а(х, у) и Ь(х, у) медленно меняющиеся. При этом, если значение фазовой функции 5(х, у) - л + 2-л-к, то точка Р с координатами (х, у) является центром темной полосы. Если представить функцию 1(х, у) в виде графика поверхности, то точка Р лежит на дне оврага. Значение фазовой функции S(х, у) = 2-л-к соответствует точке Q, лежащей в центре светлой полосы (точка Q лежит на хребте).

Таким образом, выделение центров полос сводится к задаче определения хребта или оврага функции двух переменных. Ниже описывается алгоритм нахождения центров полос.

- ISO -

2. ХРЕБТЫ И ОВРАГИ ФУНКЦИИ ДВУХ ПЕРЕМЕННЫХ

Прежде чем описывать алгоритм нахождения центров полос, дадим определение хребта (оврага) функции двух переменных. Заметим, что это понятие возникает только для функций нескольких переменных.

3 12

Пусть х = f (X -х ) - дважды непрерывно дифференцируемая функция. График этой функции в евклидовом пространстве К3 с декартовыми координатами (х ,х , х ) представляет собой гладкую поверхность.

Разложим функцию f(x\ х2) в ряд Тейлора в точке Р « (х1, х2),

о о

тогда

где h

f(x\x2) = f(xj.x2) + Df. Й + |-Йт. D2f-Ä + о(|ЙЦ2).

(2)

hi 1

1 112 2 2 , h = х -xq, h = x -xq - вектор-столбец,

йт

[„• ь'] -

вектор-строка,

Df

df af

дХ дх

-- - градиент функции f (х1, х2),

D2£ =

rVx' fxV

fxV fxV

- второй дифференциал функции fix1, х2), o(J|h||2) - остаточный член ряда Тейлора функции f(x\x2).

Предположим, что точка Р является невырожденной

критической точкой * 0.

о' о' #

в окрестности точки Р

функции ^х\х2), т.е. О^х*. х2) = (0 0) и <до[о2£(

В этом случае поведение определяется вторым дифференциалом х2).

Матрица х2) симметричная, поэтому ее собственные значения

действительны. Пусть Л, Л - собственные числа и в .е- соответствующие

12 1 «

собственные векторы матрицы х2) единичной длины.

Предположим, что в точке Р максимум функции, тогда точка Р является критической точкой эллиптического типа и собственные числа А^ Л2 отрицательны. Пусть 0>Л1>Х2, тогда в направлении собственного вектора е1 Функция убывает медленнее, чем в любом другом направлении. Поэтому направление вектора е является направлением хребта функции £(хо> хо) в точке максимума Р. В направлении вектора е2 функция убывает быстрее, чем в любом другом.

Рассмотрим теперь случай, когда в точке Р собственные значения Л^ Л2 Разных знаков. Пусть \>0>\, в этом случае точка Р является критической точкой гиперболического типа. Направление хребта определяется собственным

вектором е^ В направлении вектора ефункция будет возрастать, а в направлении вектора убывать.

Если в точке Р минимум функции f(x, у), тогда точка Р лежит на дне оврага функции f(x, у).

Предположим, что точка Р является вырожденной критической точкой. Пусть ранг матрицы D2f(xJ, х^) равен единице, тогда точка Р - критическая точка параболического типа. В окрестности точки Р функция f(x, у) представляет собой параболический цилиндр. Пусть собственные значения

тогда если Х^ 0, направление хребта определяется вектором t . В случае,

если Л2= О, в точке Р - овраг. Если ранг матрицы D2f(x*, х2) равен

нулю, то для определения поведения функции в окрестности точки Р

необходимо исследование производных более высокого порядка. Например,

функция z - х2у2 в точке Р = (0 0) имеет два оврага вдоль оси х и вдоль оси у.

Теперь мы можем дать определение точки, принадлежащей хребту функции f(x\x2).

Точка Р = (xj, х2) является точкой хребта функции f(x\x2), если производная функции f(x\x2) по направлению собственного вектора квадратичной формы D2f(x\ х2), соответствующего минимальному собственному значению, равна нулю, причем это минимальное собственное значение меньше нуля.

Аналогично можно определить точку, принадлежащую оврагу функции f(x\x2).

Точка Р= (х1, х ) является точкой оврага функции f(x\x2), если

12

производная функции f (х , X ) по направлению собственного вектора квадратичной формы D2f(x*, х2), соответствующего максимальному собственному значению, равна нулю, причем это максимальное собственное значение больше нуля.

Рассмотрим примеры.

Пример 1. Пусть г « f(x, y) = а-х2 + у2/а, а>1. в этом случае функция f(х, у) представляет собой параболоид с главными осями вдоль осей х и у. Овраг образуют точки, лежащие на оси у.

Градиент функции f(x, у):

D£ - ("К" w} - (2а'х 2-У/а).

Второй дифференциал функции f(x, у):

D2f = i 2;а 0 1

I 0 2/а J*

Собственные значения квадратичной формы D2f(x, у) равны Х^ - 2-а и = 2/а(Л1>Л2). Собственные векторы, соответствующие этим собственным

значениям, = ( о )' = 1 )' В точках> принадлежащих оврагу функции

f(x, у), производная по направлению вектора ej должна равняться нулю:

= (Df, е ) = 2• а-х = 0.

de

1

Таким образом, точки, лежащие на оси у, являются точками оврага функции f(x, у).

Пример 2. Пусть z = f(x, у) = а-х2 + b-y2 + с-х-у + d-x + е-у + f. Найдем точки, принадлежащие хребту функции f(x, у). Градиент функции f( х, у) :

Г df df ) - i

[ дх dy } i

Второй дифференциал функции f(x, у):

я. f 2-а с ) D f = 1 с 2 • b J'

Собственные значения квадратичной формы D2f(x, у) равны;

= а + Ь + /(а - Ь)2 + с2 , Л2 = а + b - /(а - Ь)2 + с2 (X!>X2)'

Для того чтобы в точке (х, у) был хребет функции f(x, у), собственное

значение Л2 должно быть меньше нуля, т.е. а + b<v(a - b) 2+ с2. Если а + Ь<0, то это условие выполняется. Если а + Ь>0, то должно выполняться условие с2>4-а-Ь. Собственный вектор, соответствующий собственным значени-

ям Л2 , равен

е =

2

-с 2-а-А,

. В точках, принадлежащих хребту функции

/(2-а-Л2)2+с2

f(x,y), производная по направлению вектора е2 должна равняться нулю:

М. = (Df, е) = 0.

de2

Подставляя соответствующие значения, получим уравнение прямой: -с-Л2'Х + (4-а-Ь - с2 + Л2) -у - d-c + е-(2-а - Л2) = 0.

Таким образом, точки, лежащие на прямой, определяемой этим уравнением, являются точками хребта функции f(x, у).

Если с = 0, то собственные значения квадратичной формы D f(х, у) равны Х= 2-а и Х2 = 2Ь. Точки хребта функции f(x, у) в этом случае лежат на прямой х = - d/(2-a), если а>Ь, и у = -е/(2-Ь), если а<Ь. Если а = Ь, то хребтом можно считать любую прямую, проходящую через точку с координатами х = - d/(2-a) и у = - е/(2•а).

3. АЛГОРИТМ ВЫДЕЛЕНИЯ ЦЕНТРОВ ПОЛОС

Выделение центров полос производится в несколько этапов:

- аппроксимация функции яркости i(x, у) в пределах окна W некоторой функцией 1(х, у)?

- определение точек, принадлежащих хребту (оврагу) функции 1(х, у)?

- вычисление расстояния р от центра окна до хребта (оврага) функции 1(х,у) ;

- определение того, принадлежит ли центр окна хребту (оврагу) функции яркости i(x,у) или не принадлежит. Если расстояние р меньше или равно пороговому R, то считается, что центр окна принадлежит хребту (оврагу). Если расстояние р больше порогового R, то центр окна не принадлежит хребту (оврагу) функции яркости 1(х, у).

Нетрудно заметить, что ширина линии центров полос в этом случае будет равна R/2.

Пусть I(m, п) - дискретные значения функции яркости 1(х, у) в точках целочисленной решетки (т, п). Методом наименьших квадратов локально аппроксимируем функцию яркости I(т, п) в окне W функцией 1(т, п):

I (т, п) = а-ш2 + b-n2 + с-го-п + d-m + e-n + f . (3)

Коэффициенты а,Ь,с, d, е, f определяются из условия минимума средней квадратической погрешности:

с2 =1 (1(ш, п) - 1(ш, п))2 (4)

( m,n)€W

Предположим, это координаты центра окна m = 0, п = 0 и окно симметрично относительно центра, тогда для определения коэффициентов получим систему уравнений:

го4 + ъ-1 m2-n2+ ш2 = У ш2-1(ш, п),

Н W WW

а.у го2.п2 + Ь.У»4+ f.y п2 = У п2-1(ш, п), * Www

O'l Го2-п2 = у m- n-1 (m, n), H w

m2 = £ m» I (m, n), и w

n2 = £ n-1 (m, n), w w

a-Y m2 + b-£ n2 + f 1 = £ I(m, n). w w w w

Для окон размерами 3x3 и 5x5 маски для вычисления коэффициентов а, Ь,

с, d, е, f приведены в приложении 1.

Размеры окна следует выбирать так, чтобы расстояние между полосами было приблизительно равно размерам окна.

Зная коэффициенты а, Ь,с, d, е, f, можно определить прямую, на которой лежат точки хребта аппроксимирующей функции 1(ю, п) (см. пример 2). Найдем расстояние р от центра окна W до этой прямой. Если расстояние р меньше или равно пороговому значению R, то центр окна принадлежит хребту функции.

Вычисление коэффициентов при аппроксимации функции яркости I(m, п) отнимает большую часть времени обработки.

Априорная информация о виде фазовой функции 6(х, у) дает возможность выбрать форму окна W и аппроксимирующую функцию i(m. п) специальным образом, тем самым сократить время обработки.

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

4. АЛГОРИТМ СКЕЛЕТИЗАЦИИ БИНАРНЫХ ИЗОБРАЖЕНИЙ С СОХРАНЕНИЕМ ВОСЬМИТОЧЕЧНОЙ СВЯЗНОСТИ

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

В зависимости от локального окружения центрального элемента окна по таблице выбирается следующее значение центрального элемента. Алгоритм

обеспечивает симметричное утончение линий, для этого используется набор переключающихся масок.

В приложении 2 представлен набор масок, которые изменяют значение Центрального элемента окна. Маски представлены в порядке возрастания

количества точек, принадлежащих линии. Маски указываются с точностью до поворота вокруг центрального элемента.

Значение 1 соответствует точке, принадлежащей линии, а значение О

соответствует фону.

Если значение центрального элемента окна равно 1 (линия), то из таблицы выбирается соответствующее следующее значение центрального элемента окна. Если значение центрального элемента окна равно 0 (фон), то никаких действий не производится.

Таким образом, точки,принадлежащие линии, не добавляются. Если при скелетизации не удалено ни одной точки после двух проходов, то этот процесс прекращается.

Маски выбираются таким образом, чтобы сохранялась восьмиточечная связность непрерывных линий единичной толщины. Рассмотрим маски, содержащие одну, две, три точки линии (см. маски 1, 2, 3 приложения 2). Одиночные точки удаляются маской 1. Концевые точки линий удаляются масками 2 и 3.

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

— 0 I 0 0

0 1 А 1 в 0

0 0 1 0

Рис. 1. Применение переключающейся маски

5. РЕЗУЛЬТАТЫ ОБРАБОТКИ ИНТЕРФЕРОГРАММ

Ниже приводятся результаты обработки интерферограммы лопатки компрессора газотурбинного двигателя.

Ввод изображения интерферограммы осуществлялся специализированным модулем ввода изображений в ПЭВМ типа 1БМ РС /АТ /ХТ. Размеры изображения

256 х 256 отсчетов, 256 градаций яркости (8 бит на отсчет). Результаты обработки интерферограммы лопатки газотурбинного двигателя представлены на рисунках 2,3,4.

Исходное изображение интерферограммы лопатки компрессора газотурбинного двигателя представлено на рис.2. Интерферограмма получена методом усреднения во времени. Результат выделения центров полос окном 7x7 с пороговым расстоянием И = 0, 5 представлен на рис. 3.

Рис. 2. Исходное изображение интерферограммы

Рис. 3. Результат применения алгоритма выделения центров полос

На рис.3 видны темные линии - центры темных полос, а также большое количество ложных точек. Ложные точки удаляются алгоритмом скелетизации. Результат представлен на рис.4.

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

Рис. 4. Результат применения алгоритма скелетизации

Литература

1. Островский V.U., Щепнов В. П., Яковлев В. В. Топографические методы измерения деформаций. М. : Наука, 1988. - 248 с..

2. Ярославский Л. П., $аянс A.M. Исследование возможностей обработки и анализа интерферограмм на ЦВМ. // Иконика. Цифровая голография. Обработка изображений. И. : Наука, 1975, с. 24-49.

3. Храмов А.Г. Методы восстановления интерферограмм на ЭВМ. // Оптическая запись и обработка изображений: Сб. науч. тр. Куйбышев: КуАИ, 1988, с. 123-128.

4. Виттих В.А., Сергеев В.В., Сойфер В. А. Обработка изображений в автоматизированных системах научных исследований. М. : Наука, 1982. - 214 с.

5. Вест Ч. Топографическая интерферометрия. М. : Мир, 1982.-504 с.

ПРиЛОЖЕШЕ 1

Маски для вычисления коэффициентов a,b,c,d,e,f для окна 3x3 и 5x5

111 ■2 -2 -2 111

ь-4-

1 -2 1 -2 1 -2

10-1 ООО -10 1

'-2 -1 0 1 2 '-13 2 7 2 -13

1 -2 -1 0 1 2 •1 2 17 22 17 2

50 -2 -2 -1 -1 0 0 1 1 2 2 Ь: 175 • 7 22 2 17 27 22 22 17 7 2

-2 -1 0 1 2 -13 2 7 2 13

ПРиЛОЖЕШЕ 2

Набор масок логического фильтра. Приводятся маски, которые изменяют значение центрального

элемента окна

ООО 0 10 ООО

2.

10 0 0 10 0 0 0

3.

4.

110 110 0 0 0

Переключающиеся маски

0 10 0 11 0 0 1

10 0 110 0 10

5.

111 0 11 0 0 0

6.

111 0 1 1 0 0 1

Переключающиеся маски

111 1 1 1 0 0 0

0 0 0 1 1 1 1 1 1

7. Переключающиеся маски

111 1 1 1 10 0