Анализ изображений непрозрачных объектов в задаче оптимизации фокусировки оптических систем Текст научной статьи по специальности «Компьютерные и информационные науки»

Анализ изображений непрозрачных объектов в задаче оптимизации

фокусировки оптических систем

А. И. Чуличков1а, Д. А. Попов2-6

1 Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра компьютерных методов физики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

2 Московский институт электроники и математики НИУ ВШЭ, факультет прикладной мате-

матики, кафедра «Кибернетика». Россия, 109028, Москва, Б. Трехсвятительский пер., д. 3.

E-mail: а achulichkov@gmail.com, 6 me@d-popov.ru

Статья поступила 26.09.2014, подписана в печать 29.10.2014.

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

Ключевые слова: оптическая микроскопия, автоматическая фокусировка, оптимизация, операторы Собела, гауссов пучок.

УДК: 519.6. PACS: 02.70.-c, 06.90.+V, 42.30.-d.

Введение

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

Экспериментальный процесс можно разделить на три основных этапа, определяющих общий технологический цикл [1]:

1) загрузка кассеты с набором образцов в механизированный оптический микроскоп;

2) получение (регистрация) подробных снимков каждого образца;

3) автоматический анализ полученных снимков с использованием специализированных вычислительных комплексов.

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

Обычно в современных установках, например IMSTAR PATHFINDER, процесс первичной фокусировки оптической системы занимает приблизительно одну минуту. Более того, для каждого образца необходимо сделать определенную серию снимков.

При стократном увеличении и в предположении, что размер исходного образца имеет размер порядка 1 см2 , для обработки всего одного образца требуется порядка 10 000 снимков. Кроме того, при перемещении предметного столика микроскопа происходит нарушение фокусировки оптической системы как из-за естественной неравномерности движения самого предметного столика, так и из-за неидеально гладкой поверхности образца. Поэтому для каждого снимка необхо-

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

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

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

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

В работах [2-4] описаны, протестированы и ранжированы по специально разработанной методологии восемнадцать алгоритмов программной автофокусировки, наиболее часто используемых в цифровой микроскопии при исследовании различных медико-биологических проб [5-14]. Анализ совокупности указанных работ позволил прийти к следующим выводам.

Во-первых, ни один из восемнадцати применяемых на практике алгоритмов фокусировки не является уни-

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

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

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

1. Описание предлагаемого подхода

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

Для построения такой модели ранее в статье [15] в рамках положений волновой оптики описана математическая модель оптической системы, формирующей изображение плоского объекта в ситуации, когда положение фокуса не совпадает с плоскостью объекта [16]. Эта модель характеризуется распределением яркости изображения точечного источника, расположенного вне плоскости фокусировки оптической системы. Такое изображение назовем пятном размытия. Как показано

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

Рис. 1. Пример исследуемого изображения

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

Рис. 2. Изображения / и ¡в (вверху), (внизу)

и £о,

едшу

2. Класс исследуемых объектов и класс их изображений

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

I(x, r) =

Io

1 Fx +( А + Jpi)

x exp

Ро 1 - 2x + U + '

(1)

g(i, j; t) = с

f (i', j') J_ e-(B/t2)[(i-i')2+(J-j')2] di' dj' t2

(i, j) G X, t G [tb t2],

где в соответствии с (1)

t(x) = 1 - Ix + (F2 + k2p0)

(2)

множество Л£ есть подмножество поля зрения X, занимаемое £ -м объектом. На рис. 1 объекты изображены черным цветом, им в модели приписывается яркость, равная 1, а фон, изображенный белым, имеет яркость, равную нулю. Тогда регистрируемое изображение имеет вид

go(i, j; t) = C

" i

Y,Xk(i', j') 1 e-Bt2)[(-i )2+(-j')2] di' dj' k=1

где х — продольная координата, определяющая точку пересечения оси пучка с плоскостью изображения и зависящая от положения оптической системы вдоль оптической оси; 10 — интенсивность светового пучка; г — радиус пучка; Р — фокусное расстояние оптической системы; £ — волновое число; р0 — радиус (ширина) исходного гауссова пучка.

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

Пусть распределение яркости плоского объекта на поле зрения X = {(г,]) & [а, Ь] х [с, d]} (геометрический формат снимка) задается функцией [(•, •): ^ Ль тогда результат регистрации изображения объекта имеет вид

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

Пусть рассматриваемые объекты являются плоскими, расположенными в предметной плоскости, перпендикулярной оптической оси видеосистемы (оптической системы микроскопа). Математическую модель такого объекта зададим в виде

п

I(г, = Х£(г,;), (г,¡) & х.

и=1

Здесь х£(•, •) — индикатор множества Л£ С Л2, £ = 1,..., п, причем эти множества не пересекаются,

(i, j) G X, t G [t1, t2]. (3) Пусть расстояние между множествами Ak с R2 и их характерный размер много больше характерного радиуса пятна размытия для каждой локальной композиции множеств Ak, а их границей являются гладкие замкнутые кривые без каких-либо наложений. Такие предположения оправдываются имеющимися массивами данных (см. рис. 1). Тогда вдоль направления внешней нормали к касательной к границе фиксированной области Ako в точке (i, j) яркость размытого изображения goO, ■; t) можно аппроксимировать функцией

—z(ij)

q(z; i, j) = C0 e-u2/(2t2) du.

—^o

Здесь точка z(i, j) = 0 соответсвует координатам (i0, j0) границы области Ako. Эта аппроксимация возникает при замене области интегрирования в (3) на полуплоскость с границей, совпадающей с касательной к границе области Ako в точке (i0, j0).

Вычислим распределение яркости на поле зрения размытого изображения, т. е. меры (площади) подмножеств B(I, I+S) поля зрения, в точках которых яркости размытого изображения лежат в пределах от I до I + S. Мера множества B(I, I+S) при указанных предположениях может быть оценена путем изменения упорядочения точек поля зрения: выберем прямоугольник Y в R2, одна сторона которого (вдоль оси OX) равна суммарной длине границ всех множеств Ak с R2, k = 1,..., n, а вторая (вдоль оси OY ) такова, что площади множеств X и Y совпадают. На этом прямоугольнике выделим подобласть D прямоугольной формы, одна из сторон которой (вдоль оси OX ) совпадает со стороной прямоугольника Y , а площадь прямоугольника D равна площади всех объектов, и построим изображение

fD(i, j) = Xd(i, j), (i, j) G X, где xD — индикатор множества D. Изображение fD(i, j) эквивалентно изображению f(i, j) в том смысле, что меры множеств Лебега Mf (c) = {(i, j) : f (i, j) < c}) и MfD (c) = {(i, j) : fD (i, j) < c}) изображений f (i, j) и fD(i, j) совпадают:

mes Mf (c) = mes MfD (c), c G Ri.

Пример эквивалентных изображений f (i, j) и fD (i, j) показан на рис. 2, вверху.

Тогда при отмеченных условиях это же свойство эквивалентности с некоторой точностью будет выполнено для изображения g0(i, j; t), определенного в (3):

go, ецшуО\ j; t) = C

Xd (i'

j,) t2e

-wm- ' )2+(j—j ' )2] di 'dj ',

(i, j) G X, t G [t1, t2].

X

2

x

2

1

r

2

x

2

x

На рис. 2, снизу также показан пример эквивалентных изображений go и go, ецшу .

Последующее вычисление яркости изображения g0, I; 0 приводит к следующему результату:

то 0

g0, ецшу(г',/; 0 = С

е-(в//2т-')2+^-п2\ м ^ =

—оо —оо

2= С0Ф ("О

(4)

где функция Лапласа Ф(г) определяется стандартным

г 2

образом: Ф(г) = \ ^ле—и /2 йи, -то < г < то. Следовательно, искомая мера подмножеств В(1,1+6) может быть оценена, если в качестве размытого изображения выбирается изображение g0> еЦшУ(г,|; ^. На рис. 3 показана зависимость функции goI ечи1у(-,|; ^ от аргумента г при фиксированном I. Так как g0(г,|; ^ не зависит от |, показана лишь зависимость от одной переменной г.

'1

1-

0.5 \

I

-4 -2 0 2 4

Рис. 3. Вид функции goIеЧи1у(г,О в зависимости от г

3. Количественная характеристика распределения яркости дефокусированного изображения

Рассмотрим в качестве характеристики распределения яркости дефокусированного изображения меру множества Q(t, /), в точках которого квадрат нормы градиента яркости этого изображения не меньше, чем задаваемое пороговое значение /:

Q(t, /) = {('• I) е У: ||Vgo(г,|; Щ2 > 1}.

(5)

Поскольку при указанных предположениях производная изображения [(г, |) по направлению нормали к границе объекта практически совпадает с производной изображения [о(г, |) по направлению оси ОУ, а по направлениям, перпендикулярным к ним, равна нулю, то, заменяя норму градиента изображения g0(•, •; t) на норму градиента введенной функции goI ечи1у0, •; t) и учитывая (4), получим

Vgo, еяи!у(г,I; 0 = (С0е—2/^,0),

На рис. 4 показаны графики первой координаты градиентов Vg0(i, |; t) как функций аргумента г для нескольких фиксированных значений параметра t. Чем меньше значение t, тем меньше мера множества Q(t,/).

Рис. 4. Графики первой координаты градиентов Vgo(í, 0 как функций аргумента г для нескольких фиксированных значений параметра t

Для дальнейшего удобно выбрать значение порога в виде / = (5С0)2, 0 < 6 < 1, тогда множество (5) запишется в виде

Q(t> 1)/=^ = Q6(t) = {('"• I) е У: е-1'2/2t2 ^ 6}, 0< 6< 1.

Мера множества Q6(t) равна 2й6(t)h, где к — размер множества У по оси ОУ, а с!6(^ определяется из соотношения 6 = е-'й(^ , эквивалентного при t ^ 0, й ^ 0 равенству

= (6)

Как видно из (6) и из рис. 5, минимум параметрической функции й^)6 по t (а значит, и мера множества Q6(t)) достигается при минимальном значении радиуса пятна размытия t = ^¡п, а t определена в (2) как функция

Рис. 5. Графики функций +й(^х)) и -й^(х))

е

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

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

На практике в качестве меры фокусировки удобнее пользоваться отношением p(x) меры дополнения множества Qs(t(x)) до всего поля зрения к мере Qs(t(x)) при фиксированном S. Действительно, согласно (2), (6), с учетом равенства mes Qs(t) = 2ds(t)h это отношение можно записать в виде

p(x) =

Sy

Sy - 2t(x)h^/\n(S) =

2t(x)hy/\n(ô) = 2t(x)hy/\n($j

- 1, (7)

где Sy — площадь всего поля зрения X. Далее, полагая, что функция t(x) в окрестности точки минимума локально квадратична и заменяя ее на полином второй степени вида t(x) = k\(x — х^)2 + k2, получим из (7) аппроксимацию для p(x), пропорциональную функции, задающей распределение Коши

p(x) '

a

1 +

(.X XminïL b2

(8)

с параметрами а, Ь и хт^п, причем значение параметра хтп дает искомое положение фокуса.

Применим описанные здесь свойства для нахождения положения фокуса по набору изображений Л(хп), полученного с разными значениями Х[,...,хN координаты пересечения плоскости изображения с оптической осью, п = 1,..., N.

Для практического нахождения положения фокуса системы при сделанных предположениях для вычисления градиента Vgo(i,¡; 0) целесообразно использовать фильтр Собела [17]. Результат фильтрации осуществляется применением двух операторов Собела:

Sx =

y

-1 0 1

-2 0 2

V-1 0 VJ -1 -2 -1 000

* Л(Хп),

1

2

1

(9)

* Л(хп).

Символ * обозначает операцию свертки изображения Л(хп) с соответствующим оператором Собела. Построим изображение

Sn(i, ]) = у/&хп а, ]))2+ №, ]))2.

Значение этого изображения в точке ]) поля зрения дает оценку градиента исходного изображения Л(хп) в этой же точке.

В получаемом изображении Sn наибольшую яркость будут иметь точки, лежащие на границе объектов на изображениях Л(хп). Далее преобразуем изображение Sn, нормируя его так, чтобы максимальная яркость результата равнялась единице, и приводя его к бинарному изображению Бтп, в котором «белыми»

становятся пиксели, в которых яркость нормированного изображения Sn превышает порог 5, и черными — остальные. В результате каждое изображение Binn будет содержать в себе bn «черных» пикселей и wn «белых» пикселей. Каждому изображению Л(хп) сопоставим значение

bn

p(xn) = —.

Wn

Согласно приведенным в начале этого параграфа соображениям, зависимость p(xn) описывается распределением Коши (8).

Процесс поиска сфокусированного изображения (фокуса) из имеющейся дискретной серии цифровых изображений {Л(хп) , n = 1,..., m} по набору значений p(xn) сводится к выбору параметров распределения Коши (8), наилучшим образом приближающего зависимость p(xn), n = 1,..., n.

На рис. 6 точками показаны значения p(xn), полученные из экспериментальных данных с оптического микроскопе IMSTAR PATHFINDER. Сплошная кривая является графиком функции p(x) в (8) с параметрами, подобранными методом наименьших квадратов. По оси абсцисс расположены номера снимков: первый снимок (на рис. 6 ему соответствует точка с номером —25) — начальное положение микроскопа вдоль оптической оси, пятидесятый снимок — конечное положение. Различия между экспериментальными данными и графиком показывают расхождения реальных данных с предложенной моделью.

о p(xJ

0.025 <

/ " 7 ♦ о/ У 0.015 о/ А

У 0.010

У* 0.005 О-1-1-1-1- X, X >Vl "

-20

-10

0

10

20

Рис. 6. Значения p(xn), полученные из экспериментальных данных с оптического микроскопа IMSTAR PATHFINDER (точки) и их наилучшее приближение мнк по формуле (8) (сплошная кривая)

Заключение

В настоящей работе был проведен анализ структуры яркости изображения непрозрачных объектов, формируемых расфокусируемой оптической системой, в зависимости от удаления положения оптической системы от точки ее фокуса. На основе функции рассеяния точки (1), полученной в [15], была введена функция, характеризующая близость оптической системы к сфокусированному состоянию. Функция использует суперпозицию операторов Собела (9) и операцию поро-

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

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

Работа выполнена при финансовой поддержке РФФИ (грант 14-07-00409а).

Список литературы

1. SchlagJ., Sanderson A., Neuman C., Wimberly F. // Robotics Institute. 1983. Technical Report CMU-RI-TR-83-14. P. 510.

2. Sun Y., Duthaler S, Nelson B. // Microscopy Research and Technique. 2004. 65. P. 139.

3. Santos A., Solorzano C., Vaquero J., Pena J. // J. of Microscopy. 1997. 188. P. 264.

4. Liu X., Wang W., Sun Y. // J. of Microscopy. 2007. 227. P. 15.

5. Беззубик В.В., Белашенков Н.Р., Устинов С.Н. // Опт. журн. 2009. 76, № 10. С. 16.

6. Brenner J., Dew B., Horton J. et al. // J. of Histochemistry and Cytochemistry. 1971. 24. P. 100.

7. Krotkov E. // Intern. J. of Computer Vision. 1987. 1. P. 223.

8. Yeo T., Jayasooriah S., Sinniah R. // Image and Vision Computing. 1993. 11. P. 629.

9. Nayar S., Nakayawa Y. // IEEE Tras. on Pattern Analysis and Machine Intelligence. 1994. 16. P. 824.

10. Subbaro M., Choi T., Niksad A. // Optical Engineering. 1993. 32. P. 2824.

11. Yang G., Nelson B. // Proc. of IEEE/RSJ. Intern. Conf.on Intelligent Robots and Systems. 2003. P. 2143.

12. Yang G, Nelson B. // Proc. of IEEE/RSJ. Intern. Conf. on Intelligent Robots and Systems. 2003. P. 3200.

13. GroenF., Young I., Lighart G. //Cytometry. 1985. 6. P. 81.

14. Vollath D. // Journal of Microscopy. 1987. 147. P. 279.

15. Попов Д.А. // Качество, инновации, образование. 2014. № 7. С. 62.

16. Захарченко А.А., Чуличков А.И. // Измерительная техника. 2007. № 1. С. 14.

17. Гонсалес Р., Вудс Р. Цифровая обработка изображений. М., 2005.

An analysis of opaque objects in the problem of the optimization of the focusing of optical systems A.I. Chulichkov1a, D.A. Popov2b

1 Department of Computer Methods in Physics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

2 Department of Cybernetics, Faculty of Applied Mathematics, HSE Moscow Institute of Electronics and Mathematics, Moscow 109028, Russia.

E-mail: a achulichkov@gmail.com, b me@d-popov.ru.

An urgent problem of the optimization of automated optical microscopy is considered. An approach is suggested during which variations in the brightness of a Gaussian beam in microscope optics are analyzed. A discrete (digital) model, which correlates with the continuous one, is introduced. An analytical approach to finding the maximum of a discrete function that estimates the measure of optical system focusing is suggested; it is based on a discrete model of Gaussian beam behavior and requires a much smaller number of optical microscope images, which speeds up the process of the automated focusing of an optical microscope.

Keywords: optical microscopy, automated focusing, optimization, Sobel operators, Gaussian beam. PACS: 02.70.-c, 06.90.+v, 42.30.-d. Received 26 September 2014.

English version: Moscow University Physics Bulletin 1(2015). Сведения об авторах

1. Чуличков Алексей Иванович — доктор физ.-мат. наук, профессор, профессор; тел.: (495) 939-41-78, e-mail: achulichkov@gmail.com.

2. Попов Дмитрий Андреевич — аспирант; e-mail: me@d-popov.ru.