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

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

Методы поиска и выделения зон ишемии на ангиографических изображениях глазного дна в литературе найти не удалось. Были изучены подходы к решению задачи выделения сосудов на изображениях глазного дна человека. В работе J. Soares и соавт. [9] проводится классификация пикселей на «сосуды» и «не сосуды». В качестве признаков авторы данной статьи используют двумерные вейвлеты Морле, а в качестве классификатора — байесовские классификаторы. В работе K. Vermeer и соавт. [10] используется метод «слежения по линии»: от начальной точки, принадлежащей сосудам, идет движение по соседним пикселям, если они также принадлежат сосудам. Другие авторы [11] используют метод «согласованных фильтров» для выделения сосудов на изображении. Расширенный согласованный фильтр с первой производной от гауссова фильтра позволяет выделить границы сосудов. Во всех перечисленных статьях в качестве изображений глазного дна использовались изображения из открытого набора данных «The STARE» (STructured Analysis of the REtina; https://cecas.clemson.edu/~ahoover/stare/), сделанные на ретинальную (фундус) камеру. Изображения (флюоресцентные ангиограммы глазного дна) получены с помощью фундус-камеры после внутривенного введения контрастного вещества. Применение перечисленных методов на данных изображениях не приводит к хорошим результатам. Анализ данных литературы показал, что методы, используемые для идентификации сосудов на двумерных изображениях, базируются на локальных признаках изображения, учитывающих специфические свойства васкулярных сегментов. Выделяют два подхода к идентификации таких областей, основанные на двух различных методах классификации: методы, основанные на обработке пикселей, и методы «отслеживания». Первый класс методов разделяется на методы, «основанные на ядрах» и «основанные на классификаторах» [10], второй класс — на методы, «основанные на правилах», и «обучаемые», упомянутые в работе J. Staal и соавт. [12]. S. Chaudhuri и соавт. [13] была предложена концепция обнаружения «согласованного фильтра». Методы отслеживания начинаются с локализации (вручную или автоматически) точек сосуда, используемых для выделения сосудистой сети, путем измерения некоторых локальных свойств изображения. Методы, оценивающие только пиксели, близкие к исходным позициям, также называются поисковыми алгоритмами. Большинство методов используют для характеризации контуров сосуда гауссовы функции. Некоторые авторы [14] предлагают алгоритм, который инициализируется определением начальной и конечной точек, и за ним автоматически следует применение фильтра для определения границ сосуда, отслеживания средней линии и извлечения клинически важных параметров.

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

Материал и методы

Исходные данные для исследования (ангиограммы с различными заболеваниями глаз; n=117) предоставлены отделом современных методов лечения в офтальмологии ФГБНУ «Научно-исследовательский институт глазных болезней». Экспериментальная проверка нового метода и интерпретация результатов проведены совместно с группой ученых из ФГУ «Федеральный исследовательский центр «Информатика и управление» РАН». Разработка математического метода велась совместно с учеными из Белорусского государственного университета.

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

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

A. Предобработка изображения — необходима для повышения качества изображений таким образом, чтобы повысить эффективность их дальнейшего анализа и распознавания [4].

Шаг A1. Яркостная и цветовая коррекция изображений. Исходное изображение (рис. 1, а) дано в монохроматическом виде, но с тремя каналами цветов. Поэтому исходное изображение переводится в одноканальное изображение в оттенках серого с помощью уравнения, где R, G, B — красный, зеленый и синий каналы изображения соответственно.

Шаг A2. Нормализация изображения. Для дальнейшего удобства работы с изображением оно нормализуется на отрезок [1].

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

,

где gain — фактор масштабирования, Input — пиксели входного изображения, Output — пиксели выходного изображения.

Рис. 1. Предобработка изображений.

а — исходное изображение (флюоресцеиновая ангиограмма глазного дна); б — предобработанное изображение; в — ядра с откликом на линии в определенной ориентации; г — найденные края; д — обработанное изображение; е — бинаризация Саувола.

Шаг A4. Фильтрация шума на изображении. Исходное изображение достаточно шумное, и с целью уменьшения шума применяется фильтр Гаусса [15] (рис. 1, в). Фильтр Гаусса также хорошо способствует выделению границ, поскольку уменьшает чувствительность детектора границ к шуму.

Шаг A5. Повышение контрастности изображения. Мелкие сосуды представляют собой тонкие структуры и являются низкоконтрастными, поэтому их выделение — трудоемкая задача. Для улучшения их отделимости от фона и лучшего, более точного выделения на изображениях делали свертку ядрами для поиска краев [16], соответствующих четырем ориентациям (горизонтальный, вертикальный, 45° и 135°). После свертки данными ядрами для изображения находили края (рис. 1, г) и суммировали с обработанным на предыдущем шаге изображением (рис. 1, д).

Шаг A6. Бинаризация изображения. К полученному изображению применяется бинаризация Саувола [17]. Так как изображение неоднородное по освещенности, хорошо подходит адаптивная бинаризация. Для данного изображения применялась бинаризация Саувола (рис. 1, е).

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

Шаг B1. Удаление шумовых объектов по площади. Вместе с сосудами также выделяются шумы (мелкие компоненты) из-за неоднородности освещенности. Эти компоненты удаляют из изображения по площади (рис. 2, а).

Рис. 2. Анализ изображений.

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

Шаг B2. Удаление шумовых объектов по удлиненности. В получившемся после очистки на предыдущем шаге изображении все еще остаются компоненты, которые не являются сосудами. Все эти элементы удаляются по порогу удлиненности области:

.

Поскольку эти объекты менее вытянуты и больше похожи на круг, подобрать такой порог можно без риска потери сосудов. Результат в сравнении с исходным изображением показан на рис. 2, б, в.

Шаг B3. Построение остова области. Строим скелет для выделенных сосудов для дальнейшего измерения их ширины [18]. На рис. 2, г показан скелет некоторой области изображения, а на рис. 2, д — скелет, наложенный на изображение с сосудами.

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

Шаг C1. Подсчет плотности расположения сосудов. Отношение выделенных на шаге B2 пикселей, относящихся к сосудистому руслу, к общей площади исходного изображения является величиной, по которой можно судить о плотности расположения сосудов. Автоматический расчет показателей плотности сосудов определяется по формулам:

,

где P_vessel — сосудистая плотность, P_vessel — топологическая плотность сосудов, площади скелета сосудов), Area_common — общая площадь исследуемого региона. Эти характеристики играют важную роль при описании глобальных свойств сосудистой сети и являются производными для ряда других характеристик. Например, средний диаметр сосуда (D_vessel) определяется как отношение этих плотностей:

.

Шаг C2. Подсчет диаметра сосудов. По построенному остову считается ширина сосудов для каждой точки скелета. Определение реальной ширины сосудов — сложная процедура, что связано, в первую очередь, с непостоянным диаметром сечения кровеносных сосудов. На первом шаге по бинарному изображению сосудов, полученному ранее, строится дистанционная карта, в которой каждому пикселю ставится в соответствие расстояние до края сосуда в соответствии с метрикой Чамфера [19] (оптимальная для задачи метрика в соотношении время/качество). На втором шаге выполняется утоньшение [20] бинарного изображения сосудов, полученного ранее. На третьем шаге остовы сосудов и дистанционная карта позволяют получить оптимальный по точности набор пикселей со значениями, соответствующими значениями ширины сосудов. Анализ дистанционной карты позволяет определить характер изменения диаметра сосуда по его длине (рис. 3, а, б).

Рис. 3. Подсчет морфометрических характеристик сосудов.

а — исходное изображение; б — карта расстояний.

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

Рис. 4. Распределение ширины сосудов.

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

По известным диаметрам сосудов в каждой точке скелета вычисляются статистические характеристики: средний диаметр и среднеквадратичное отклонение.

Шаг C3. Подсчет длины сосудистого русла. При известном масштабе исходного изображения длина построенного остова в пикселях переводится в длину сосудистого русла в микрометрах. Длина обычного объекта определяется как максимальное расстояние между двумя точками, в то же время основной физический смысл имеют длина и толщина развернутого, вытянутого сосуда. Такая длина соответствует размерам средней линии — скелету, а толщина — минимальному расстоянию между двумя граничными точками внутри объекта. В самом грубом случае длина определяется как площадь скелета:

.

Однако изгибы сосуда вносят изрядную долю погрешности, поэтому более точным считается следующее определение длины:

,

где perimeter — периметр сосуда.

На основе разницы между этими понятиями длины можно определить ряд новых параметров, например, скрученность (curl) волокна — отношение длины объекта к длине волокна:

,

где Length — простая длина сосуда как бинарного объекта.

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

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

Шаг A0. Выделение области интереса. Для поиска зон ишемии рассматривается область картинки с обрезанными краями в 17% (высчитано эмпирически). Такое урезание картинки связано с тем, что при получении ангиографических изображений глазного дна края изображения из-за округленности глазного яблока освещаются слабо и мешают определению зон ишемии.

Шаги A1—A4 аналогичны описанным выше в методе выделения сосудов.

B. Анализ изображений (выделение зон ишемии).

Шаг B1. Построение гистограмм яркостей изображения. Изображение на предыдущих шагах нормализовано на отрезок [1]. Будем делить отрезок на 256 бинов (bins) и строить гистограмму для этих подотрезков (рис. 5, а).

Рис. 5. Анализ изображений.

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

Шаг B2. Анализ гистограмм яркостей изображения. Теоретически, учитывая, что зоны ишемии — самые темные области на изображении, на гистограмме первый локальный максимум будет отделять эти области от остальных освещенных областей. Валидность данного деления была проверена на выборке из нескольких изображений. Экспериментальным путем было доказано, что действительно данный экстремум будет порогом для нахождения зон ишемии.

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

Шаг B3. Наложение сосудов на область зоны ишемии. Зоны ишемии находятся в другой плоскости относительно сосудов, поэтому необходимо учитывать (добавлять) зоны ишемии, пересекающиеся с сосудами. Выполняем заливку области ишемии, пересекающейся с сосудами. На рис. 5, в показано наложение найденных сосудов на выделенные зоны ишемии.

Шаг B4. Расширение зоны ишемии до пересечения с сосудами. Для заполнения данных участков применим операцию дилатации для расширения зоны ишемии до пересечения с сосудами (рис. 5, г).

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

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

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

Результаты

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

Рис. 6. Алгоритмическая схема анализа ангиографических изображений глазного дна человека.

Программная реализация выполнена на языке программирования Python версии 3 с использованием библиотеки scikit-image. Интерфейс написан на языке программирования Python версии 3 с использованием библиотеки Qt. Основными шагами работы с интерфейсом являются открытие файла с исходным изображением, выделение области интереса на изображении и анализ изображения с выводом графиков и статистических данных. Реализованная программа является кроссплатформенной и запускается на любой стандартной ПЭВМ с поддержкой языка программирования Python 3.

Заключение

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

Данная работа была выполнена при частичной поддержке Российского фонда фундаментальных исследований (грант №20-57-00025).

Участие авторов:

Концепция и дизайн исследования: И.Г., В.Я., С.А.

Разработка и математический анализ, создание программного обеспечения: А.О., А.Н., А.Т.

Сбор материала, подготовка текста статьи: Н.Ф.

Авторы заявляют об отсутствии конфликта интересов.

The authors declare no conflicts of interest.