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

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 4, с. 96-102 ОБРАБОТКА И АНАЛИЗ СИГНАЛОВ =

удк 621.391.837: 681.3

© Б. В. Бардин, И. В. Чубинский-Надеждин

ОБНАРУЖЕНИЕ ЛОКАЛЬНЫХ ОБЪЕКТОВ НА ЦИФРОВЫХ МИКРОСКОПИЧЕСКИХ ИЗОБРАЖЕНИЯХ

Рассмотрена задача обнаружения локальных объектов, малых по сравнению с размерами изображений. Показана необходимость двух этапов решения задачи — предварительной обработки изображения и поиска объектов. Предварительную обработку составляет обостряющая фильтрация изображения фильтром Тихонова или Винера, предназначенная для разделения слипшихся объектов и их нормализации. Кроме того, на этом этапе устраняется фоновая составляющая изображения. Описан алгоритм поиска объектов, использующий формирование пороговой поверхности обнаружения. Экспериментально показано, что эффективность обнаружения превышает 90 % при времени обработки одного изображения размером 256 х 256 менее 60 мс.

Кл. сл.: изображения живых клеток, анализ цифровых изображений, поиск объектов, преобразование Фурье, обостряющая фильтрация, пороговая поверхность

ВВЕДЕНИЕ

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

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

Первым этапом решения задачи подобных ис-

следований являются обнаружение объектов и определение их координат в пределах кадра изображения. Алгоритм автоматизированного анализа должен обеспечивать эффективность обнаружения клеток не менее 90-95 %, при относительном количестве ложных срабатываний менее 5 % во всем анализируемом диапазоне размеров и прозрачности клеток. Характерные размеры объектов (изображений клеток) составляют 3-20 пикселей, и погрешность определения центров изображений не должна превышать 1 пиксель, что необходимо для последующего адекватного определения оптических параметров. Принципиальным требованием является скорость обработки информации: результаты анализа одного кадра должны быть получены за 50-500 мс в зависимости от размера кадра.

Пример решения задачи создания алгоритма обнаружения в данной постановке является предметом настоящей работы.

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

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

ками фрагменты, количество которых значительно превышает количество клеток.

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

пазоне 5-50 %. Отношение сигнал/шум в пикселе порядка 10. Относительная неравномерность пьедестала от центра к краям кадра может составлять до 20 %.

На рис. 1 и 2 приведены два типа изображений живых клеток реальных образцов. Нетрудно видеть ряд характерных особенностей рассматриваемых задач. В пределах исследуемой пробы присутствуют объекты различных размеров. В пробе на рис. 1 видно, что часть объектов "слипается", что препятствует их разделению при обычной пороговой обработке исходного изображения, т. к. в этом изображении между подобными объектами может отсутствовать "провал". В изображении пробы на рис. 2 у объектов визуализируется только контур, да и тот может быть с разрывом.

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

Для преодоления отмеченных трудностей задача обнаружения разделяется на два этапа: 1) предварительная обработка изображения; 2) поиск объектов.

ПРЕДВАРИТЕЛЬНАЯ ОБРАБОТКА ИЗОБРАЖЕНИЯ

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

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

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

Рис. 1. Микроскопическое изображение суспензии лейкоцитов

Рис. 2. Изображение второго типа — с контрастными границами клеток

рованию параметров. Математически фильтр Тихонова (как и фильтр Винера) формулируется в терминах преобразования Фурье, в данном случае двухмерного:

2 (и,V)=

К '(и, V)

К ( и, V )| + ХР(и2 + V2)

^(и,V), (1)

где

и, V — частоты (пространственные);

S(u, V) — зарегистрированный сигнал (изображение);

£(и, V) — оценка исходного, восстанавливаемого сигнала;

К(иу) — искажающая функция, Фурье-образ объекта;

К *(и, V) — комплексное сопряжение К(и^);

К (и, V)2 = к >, V) К (и, V);

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

5 0(и, V) = £ (и, V) - Gb(u, V) х £ (и, V) = = (1 - Gb(u, V)) х £ (и, V),

(2)

Р(и2 + V2) — сглаживающий полином, обычно первого порядка, т. е. Р = (и2 + V2);

X — параметр регуляризации.

Точное определение параметра регуляризации X требует применения итерационной процедуры, что требует значительного времени на вычисления. В большинстве рассматриваемых задач это крайне нежелательно. Вместе с тем многочисленные эксперименты показали, что во всех исследованных задачах хорошие результаты дает использование постоянной величины параметра регуляризации, порядка 10-5 -10-6.

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

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

В сложных случаях ширина полосы обостряющего фильтра может оказаться недостаточной для обнаружения всех объектов исследуемой пробы. В этом случае можно изображение обработать, например два раза фильтрами с разными образами К 1(и,V) и К2(и,V), и результаты сложить. Однако та же самая цель будет достигнута использованием одного фильтра с составным ядром К (и, v) = = К1(и,V) + К2(и,V).

где £(и, V) — исходное зарегистрированное изображение, Gb(u, V) — оператор усреднения, в общем случае взвешенного, например по Гауссу.

Подставив выражение для £0(и,V) в (1) вместо £(и, V), получаем окончательное выражение для предварительной обработки изображения:

2 (и, V) =

К *(и, V)

К ( и, V )| + Я(и2 + V2)

(1 - Gb(u , V)) £ (и, V). (3)

На рис. 3 приведены результаты обработки подобным фильтром изображения пробы рис. 2. Здесь (рис. 3) на нижних панелях показаны профили сигнала по строке, отмеченной на изображениях (на верхних панелях) маркером — горизонтальной линией. На левых панелях приведены исходные изображения и сигналы, на правых — после обработки фильтром (3). Размер изображений на панелях составляет 450 х 400 пикселей. В фильтре использовано составное двойное ядро, которое в геометрическом пространстве Щ(х, у) представляет сумму двух колец, сформированных, как описано выше, шириной 6 пикселей и диаметром одно кольцо 22 пикселя, что соответствует наиболее крупным объектам, и другое — 11 пикселей. Параметр регуляризации Х = 105. Оператор усреднения Gb(x, у) представлен функцией Гаусса шириной 40 пикселей. Как можно видеть на рис. 3, предварительная обработка изображения позволяет сформировать компактные усиленные образы объектов (масштабы профилей по вертикали для обеих панелей одинаковы). Побочным эффектом линейной фильтрации является появление боковых колебаний, которые в сравнении с максимумом сигнала являются не очень значительными, но которые необходимо учитывать при дальнейшей обработке.

о

о

"О

о

о Я,

о

о

1

с

О 'Л

о с

о

— Е г V

&

о

-о-

—

№

с с

о

о

о

с

о

_о

й г #

с

и» - *

% ч

'ук

А

Рис. 3. Обостряющая фильтрация с устранением фона изображения, показанного на рис. 2. а — исходное изображение; б — результат фильтрации

Рис. 4. Гистограммы сигналов изображения, показанного на рис. 2.

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

ПОИСК ОБЪЕКТОВ

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

так и низкочастотными, обусловленными фоновой составляющей сигнала.

Обычно величина порога определяется на основе гистограммы сигналов [1, 2]. Для изображения с однородными, достаточно протяженными объектами гистограмма имеет вид кривой с двумя (или более) пиками. При этом один пик, как пра-

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

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

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

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

изображения. И хотя анализ может быть очень простым — перебор всех точек окна, например для выявления максимального значения сигнала, — окно должно быть достаточно большим, порядка размера объекта, и поэтому такой анализ может потребовать много времени. Для того чтобы обойти эту трудность, обычно разбивают исходное изображение на ряд областей, в каждой из которых определяется свое значение порога [1]. Однако при существенной неравномерности фона возникает неприятная проблема состыковки областей.

Вместе с тем хорошо известно, что линейные операции значительно быстрее и эффективнее выполняются в пространстве Фурье. В данном случае в качестве линейной операции может быть использовано сглаживание (усреднение) положительной части сигнала на выходе обостряющего фильтра (рис. 3), с последующим умножением на подходящий коэффициент, чтобы результат соответствовал максимальному значению шума с некоторым запасом. И кроме того, поскольку пороговая поверхность меняется медленнее сигнала, можно использовать массивы данных меньшей размерности в сравнении с исходным изображением (в рассмотренном примере 450 х 400 пикселей), полученных посредством сжатия этого изображения, например при помощи бининга. При бининге сигнал каждого пикселя сжатого изображения представляет сумму (или среднее значение) сигналов соседних пикселей площадки размером 2 х 2, 4 х 4 и т. п. исходного изображения.

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

Рис. 5. Функциональная схема алгоритма поиска объектов

На рис. 5 приведена функциональная схема алгоритма поиска объектов. Здесь использованы следующие обозначения:

F и F 1 — соответственно прямое и обратное преобразования Фурье;

х, y — координаты геометрического пространства (изображения);

u, v — координаты пространства Фурье (частоты);

MN и mn — размерности массивов данных, причем m < M и n < N;

MN^ mn — операция уменьшения размерности данных, сжатие данных;

mn^-MN — операция увеличения размерности данных, восстановление данных.

Данные с выхода обостряющего фильтра Z(x, y) подаются на узел Dd (diode), который выделяет положительную полярность сигнала. После усреднения полученного сигнала фильтром с ядром Gb(u, v) и ограничения результата по минимуму величиной hmin формируется пороговая поверхность h(x, y).

По результатам сравнения сигнала Z(x, y) с пороговой поверхностью h(x, y) (операция Cmp) формируется битовая карта В(х, y), в которой значение каждого бита определяет больше (1) или меньше (0) пороговой поверхности величина сигнала для соответствующего пикселя, т. е. принадлежит ли данный пиксель изображению какого-либо объекта или нет.

По этой карте производится определение областей связности (FLR — find of link region) для пикселей с В(х, y) =1. Каждая из областей представляет отдельный объект. Здесь Obj [k] — список объектов. Алгоритм определения областей связности описан в литературе по морфологической обработке видеоинформации [1].

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

Таким образом, в результате выполнения операций поиска формируется список объектов с координатами их центров Obj[k](x, y).

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

На рис. 6 и 7 приведены результаты выявления областей изображения, в которых локализованы объекты, для проб, показанных на рис. 1 и 2 соответственно. Выявленные области закрашены белым цветом. Эти данные приведены здесь для наглядной иллюстрации результатов работы описан-

Рис. 6. Выявление областей изображения первого типа (рис. 1), в которых локализованы объекты

ного алгоритма. Анализ достаточно большого массива результатов анализа реальных образцов показал, что погрешность определения центров объектов не превышает 1 пиксель (когда не вызывает сомнений собственно правильность интерпретации истинного положения центра). Время анализа составляло менее 60 мс при размере изображения 256 х 256. Эффективность обнаружения достигала 90 % при анализе изображений обоих типов, показанных на рис. 1 и 2.

С нашей точки зрения результаты вполне удовлетворительны, и основные элементы описанного алгоритма были использованы при создании нового прибора — детектора субпопуляций клеток. Тем не менее надо отметить, что для таких сложных совокупностей объектов, как например, в пробе на рис. 1, и при наличии большого количества мелких фрагментов вряд ли возможно полностью избежать ошибок. Однако для обсуждаемых здесь задач исследования это в определенной мере допустимо, т. к. окончательными результатами анализа являются различные статистические распределения свойств объектов и возможна дополнительная коррекция и/или фильтрация ложных данных по совокупности признаков объектов.

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

СПИСОК ЛИТЕРАТУРЫ

1. Гонсалес Р., Вудс Р. Цифровая обработка изображений. Пер. с англ. М.: Техносфера, 2005. 1070 с.

2. Ярославский Л.П. Цифровая обработка сигналов в оптике и голографии. М.: Радио и связь, 1987. 296 с.

Институт аналитического приборостроения РАН, Санкт-Петербург

Материал поступил в редакцию 29.07.2009.

DETECTION OF LOCAL OBJECTS ON DIGITAL MICROSCOPIC IMAGES

B. V. Bardin, I. V. Chubinskiy-Nadezhdin

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The problem of detection of local objects, small in comparison with the image size is considered. It is shown that two-stage solution of the problem is necessary — preliminary processing of the image and the search for objects. Preliminary processing includes peaking filtration of the image by Tikhonov or Wiener filter, intended for discerning of conglutinated objects and for their normalization. Besides at this stage the background component of the image is eliminated. A search algorithm using formation of a threshold surface of detection is described. It is experimentally shown, that the detection efficiency exceeds 90 %, and the processing time is less than 60 ms for the image size 256 x 256.

Keywords: : live cell imaging, analysis of digital image, search of objects, Fourier transform, peaking filtering, threshold surface

Рис. 7. Выявление областей изображения второго типа (рис. 2), в которых локализованы объекты