ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ
«БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ
ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ»
( Н И У «Б е л Г У »)
ИНСТИТУТ ИНЖЕНЕРНЫХ ТЕХНОЛОГИЙ И ЕСТЕСТВЕННЫХ НАУК
КАФЕДРА МАТЕМАТИЧЕСКОГО И ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ
ИНФОРМАЦИОННЫХ СИСТЕМ
РАЗРАБОТКА МЕТОДА И АЛГОРИТМА СЕГМЕНТАЦИИ
КРОВЕНОСНЫХ СОСУДОВ НА ИЗОБРАЖЕНИИ ГЛАЗНОГО ДНА
Выпускная квалификационная работа
обучающейся по направлению подготовки 02.03.03 Математическое
обеспечение и администрирование информационных систем
очной формы обучения, группы 07001402
Черноморец Дарьи Андреевны
Научный руководитель
к.т.н., доцент
Михелев В.М.
БЕЛГОРОД 2018
СОДЕРЖАНИЕ
ВВЕДЕНИЕ…………………………………………………………………..
3
1. АНАЛИЗ МЕТОДОВ СЕГМЕНТАЦИИ СОСУДОВ ГЛАЗНОГО
ДНА…………………………………………………………………………..
7
1.1 Основные задачи и проблемы при исследовании снимков
глазного дна……………………………………………………………
7
1.2 Анализ методов сегментации объектов на изображении……....
13
1.3 Особенности методов автоматической сегментации сосудов
глазного дна……………………………………………………………
19
1.4 Постановка задачи………………………………………………… 25
2. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ МЕТОДОВ, ПРИМЕНЯЕМЫХ ДЛЯ
РЕШЕНИЯ ЗАДАЧИ СЕГМЕНТАЦИИ СОСУДОВ ГЛАЗНОГО ДНА...
26
2.1 Метод повышения контрастности изображений………………..
26
2.2 Основные операции морфологической фильтрации
3.
изображения…………………………………………………………...
32
2.3 Метод кластеризации k-means…………………………………...
39
РАЗРАБОТКА
АЛГОРИТМИЧЕСКОГО
И
ПРОГРАММНОГО
ОБЕСПЕЧЕНИЯ…………………………………………………………….
46
3.1 Разработка метода и алгоритма сегментации кровеносных
сосудов на изображении глазного дна……………………………….
46
3.2 Интерфейс программной реализации…………………………....
56
3.3 Программная реализация алгоритма сегментации кровеносных
сосудов на изображении глазного дна……………………………….
58
4. ЭКСПЕРИМЕНТАЛЬНЫЕ ИССЛЕДОВАНИЯ РАЗРАБОТАННОГО
МЕТОДА СЕГМЕНТАЦИИ СОСУДОВ ГЛАЗНОГО ДНА……………… 64
4.1 Выбор критериев для экспериментальных исследований……… 64
4.2 Описание вычислительных экспериментов……………………..
65
ЗАКЛЮЧЕНИЕ……………………………………………………………… 72
СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ………………………… 74
ПРИЛОЖЕНИЕ……………………………………………………………… 79
2
ВВЕДЕНИЕ
Анализ состояния кровеносных сосудов глазного дна имеет большое
значение при диагностировании и лечении различных заболеваний.
В большинстве известных методов по исследованию глазного дна
существенное значение имеет анализ его кровеносной системы. Глазное дно
является единственной частью человеческого тела, кровеносную систему
которой
можно
наблюдать
непосредственно
без
хирургического
вмешательства [1].
При корректном выделении (сегментации) кровеносных сосудов на
изображении глазного дна становится возможным поставить более точный
диагноз, что важно при лечении пациента.
Выделение
сосудистой
системы
вручную
является
достаточно
сложным процессом, занимающим существенное количество времени и сил,
и иногда невозможным из-за слишком сложной структуры сосудистого
дерева или низкого PSNR (значения отношения сигнал-шум) на изображении
[2].
Из-за
указанных
проблем
широкое
распространение
получил
компьютерный анализ изображений кровеносных сосудов глазного дна,
который стал главным инструментом медицинских диагностических систем,
позволяющих значительно повысить качество диагностики. Современная
медицина является одной из самых высокотехнологичных отраслей,
важнейшей задачей которой является разработка новых эффективных
методик ранней диагностики различных патологий. В настоящее время
практически все методы обследования в различных областях медицины
компьютеризированы, включая диагностику по глазному дну.
Распознавание
сосудов
глазного
дна
и
определение
их
морфологических признаков – основные этапы автоматизированных методик
диагностического анализа сосудистой системы, так как от точности
3
выделения и измерения её составляющих зависят результаты диагностики
[3].
Сегментация сосудов глазного дна является основным этапом
подготовки
информации
для
медицинских
автоматизированных
диагностических систем.
Результат сегментирования сосудистой сети глазного дна очень сильно
зависит от яркости и контраста исходного изображения, а хорошее качество
изображения имеет большое значение для корректной диагностики,
проводимой вручную или автоматически.
Существующие
методы
сегментации
кровеносных
сосудов
на
изображении глазного дна имеют отдельные недостатки, в них, зачастую, не
проводится
предобработка
изображения
для
устранения
различных
артефактов, снижающих качество снимка. Также большинство методов
необходимо применять в комбинации с другими методами для достижения
наилучшего эффекта сегментации, следует отметить, что в большинстве
случае выбор комбинации методов не является тривиальным [4].
Таким образом, целью работы является совершенствование методов
сегментации кровеносных сосудов на изображении глазного дна на основе
контрастно
ограниченной
морфологической
адаптивной
фильтрации,
метода
эквализации
кластеризации
гистограммы,
k-means
и
согласованной фильтрации.
Для достижения поставленной цели необходимо решить следующие
задачи:
1. Анализ методов сегментации сосудов глазного дна;
2.
Проведение
исследования
теоретических
основ
методов,
применяемых для решения задачи сегментации сосудов глазного дна;
3. Разработка метода и алгоритма сегментации кровеносных сосудов на
изображении глазного дна;
4.
Разработка
программной
реализации
кровеносных сосудов на изображении глазного дна;
4
алгоритма
сегментации
5.
Проведение
вычислительных
экспериментов
по
проверке
работоспособности разработанного метода.
Объектом изучения данного исследования является кровеносные
сосуды глазного дна.
Предметом изучения данного исследования являются методы и
алгоритмы сегментации кровеносных сосудов на изображении глазного дна и
их программирование в среде инженерных расчетов Matlab.
Работа состоит из Введения, четырех глав, Заключения, списка
использованных источников и Приложения.
Во
Введении
показана
актуальность
выполняемой
работы,
сформулированы цель и задачи проводимого исследования, приведено
краткое содержание и структура работы.
В первой главе рассматриваются основные задачи и проблемы,
появляющиеся при исследовании снимков глазного дна, структуры глазного
дна, и диагностировании болезней, которые можно выявить по снимкам
данного
органа
зрения.
Проводится
анализ
существующих методов
сегментации объектов на изображении, а также анализ методов, которые
специализируются на сегментации непосредственно кровеносных сосудов на
изображении глазного дна.
Во второй главе рассматриваются теоретические основы методов,
применяемых для решения задачи сегментации сосудов. К таким методам
относятся контрастно ограниченная адаптивная эквализация гистограммы,
морфологическая фильтрация и метод кластеризации k-means, для которых
сформулированы основные определения, вычислительные соотношения и
алгоритмы. Во второй главе также выполнен анализ преимуществ и
недостатков данных методов.
В
третьей
главе
разработан
метод
и
алгоритм
сегментации
кровеносных сосудов на изображениях глазного дна с использованием
контрастно
ограниченной
адаптивной
эквализации
гистограммы,
морфологической фильтрации, метода кластеризации k-means, согласованной
5
фильтрации отдельно толстых и тонких сосудов. Также в данной главе
описан интерфейс программы, и приведены фрагменты листингов кода
программной реализации разработанного метода.
В
четвертой
главе
проведены
сравнительные
вычислительные
эксперименты разработанного метода и известного метода сегментации
сосудов глазного дна.
В Заключении сформулированы основные результаты проведенного
исследования.
Работа выполнена на 84 страницах и содержит 29 рисунков, 3 таблицы
и 50 литературных источников.
6
1. АНАЛИЗ МЕТОДОВ СЕГМЕНТАЦИИ СОСУДОВ ГЛАЗНОГО ДНА
В данной главе рассматриваются основные задачи и проблемы,
появляющиеся при исследовании снимков глазного дна, структуры глазного
дна, диагностировании болезней, которые можно выявить по снимкам
данного органа зрения. Также проводится анализ существующих методов
сегментации объектов на изображении, а также анализ методов, которые
специализируются на сегментации непосредственно кровеносных сосудов на
изображении глазного дна. Данные методы имеют свои преимущества, кроме
того известным методам также присущи недостатки, поэтому разработка и
совершенствование методов сегментации сосудов является достаточно
актуальной задачей, которая и определяет цель данной работы.
1.1 Основные задачи и проблемы при исследовании снимков
глазного дна
Анализ состояния глазного дна имеет большое значение при
диагностировании и лечении различных заболеваний.
Во
многих
случаях
для
диагностики
заболеваний
на
основе
исследовании глазного дна используются его снимки [1].
На снимках глазного дна обычно выделяют следующие элементы
(рис. 1.1):
- сосудистая оболочка (состоит из артерий, вен различной толщины),
- диск зрительного нерва и его углубление,
- макула (желтое пятно),
- фовеа (центральная ямка сетчатки).
7
Рис. 1.1. Изображение глазного дна
Существуют несколько методов по исследованию глазного дна. К
наиболее известным относятся следующие [2]:
1. Метод офтальмоскопии при расширенных зрачках. С помощью
офтальмоскопа исследуется состояние сетчатки последовательно от центра
до периферии по всем меридианам, тщательно осматривается диск
зрительного нерва, макулярная область, разветвления крупных сосудов.
Офтальмоскопия обязательно должна проводиться квалифицированным
специалистом.
2. Фотографирование глазного дна с помощью стандартной или
немидриатической камеры. Используя данный метод, можно получить
документальную информацию о состоянии глазного дна. Данная процедура
обследования может проводиться не только офтальмологами, но и другим
медицинским
персоналом
с
последующей
расшифровкой
снимков
специалистом.
3. Флюоресцентная ангиография (ФАГ) позволяет четко обнаруживать
циркуляцию флюоресцина в ретинальных и хореоретинальных сосудистых
системах, что стоит на первом месте для диагностики функциональных и
органических изменений в структурах глазного дна. С помощью ФАГ
обнаруживают "протекающие" капилляры.
По снимку глазного дна можно выявить и локализовать [3]:
- отслойку сетчатки;
8
- разрывы сетчатки;
- новообразования и нарушения пигментации глазного дна;
- микрокровоизлияния в сетчатку и под пигментный эпителий сетчатки;
- проявления диабетической ретинопатии (поражений сетчатки при
сахарном диабете), в том числе кровоизлияния, отеки, экссудаты;
- тромбоз сосудов глаза;
- другие патологические изменения.
Для получения снимков глазного дна применяют различные способы
[2].
Глазное дно может быть сфотографировано только через отверстие
зрачка, применяя специальную аппаратуру. Изображение глазного дна
получают двумя способами: первый способ состоит в том, что световые лучи,
отраженные глазным дном, пройдя через оптическую систему глаза,
попадают на объектив и образуют изображение на пленке (рис. 1.2а).
а
б
Рис. 1.2. Способы получения изображения глазного дна
(1 – глаз, 2 – офтальмологическая линза, 3 – объектив, 4 – пленка):
а – 1-й способ, б – 2-й способ
9
Недостатком данного способа является то, что для получения
необходимого поля зрения приходится располагать аппаратуру слишком
близко к глазу больного. Поэтому чаще применяют второй способ: лучи
света сначала проходят через офтальмоскопическую линзу и образуют
первое изображение, которое затем переносится объективом на пленку (рис.
1.2б).
В большинстве известных методов по исследованию глазного дна
существенное значение имеет анализ его кровеносной системы. Глазное дно
является единственной частью человеческого тела, кровеносную систему
которой
можно
наблюдать
непосредственно
без
хирургического
вмешательства.
При корректном выделении (сегментации) кровеносных сосудов на
изображении глазного дна становится возможным поставить более точный
диагноз, что важно при лечении пациента.
Выделение
сосудистой
системы
вручную
является
достаточно
сложным процессом, занимающим существенное количество времени и сил,
и иногда невозможным из-за слишком сложной структуры сосудистого
дерева или низкого PSNR (значения отношения сигнал-шум) на изображении
[4].
Из-за
указанных
проблем
широкое
распространение
получил
компьютерный анализ изображений глазного дна, который стал главным
инструментом
медицинских
диагностических
значительно повысить качество диагностики.
систем,
позволяющих
Современная
медицина
является одной из самых высокотехнологичных отраслей, важнейшей
задачей которой является разработка новых эффективных методик ранней
диагностики различных патологий. В настоящее время практически все
методы обследования в различных областях медицины компьютеризированы,
включая диагностику по глазному дну [5].
Распознавание
сосудов
глазного
дна
и
определение
их
морфологических признаков – основные этапы автоматизированных методик
10
диагностического анализа сосудистой системы, так как от точности
выделения и измерения её составляющих зависят результаты диагностики.
Результат сегментирования сосудистой сети глазного дна очень сильно
зависит от яркости и контраста исходного изображения, а хорошее качество
изображения имеет большое значение для корректной диагностики,
проводимой вручную или автоматически [6].
Одной из проблем, влияющих на качество снимка глазного дна,
является наличие неоднородной освещенности при фотографировании.
На рис. 1.3 приведен пример изображения глазного дна с неоднородной
освещенностью.
Рис. 1.3. Изображение глазного дна с неравномерной освещенностью.
Неоднородная освещенность на снимке может возникать вследствие
особенностей строения сетчатки. Сетчатка – это сферическая поверхность,
которая
освещена
источником
света,
расположенным
в
нескольких
сантиметрах от зрачка. Камера, которая так же находится вблизи от зрачка,
фиксирует отраженный от сетчатки свет. Из-за такого расположения
элементов съемки относительно друг друга, периферийная часть поверхности
сетчатки хуже освещается, и этот регион является более темным на
изображениях глазного дна.
Также существенную проблему представляет малая величина зрачка –
отверстия, через которое необходимо не только сфотографировать глазное
дно, но также его равномерно и хорошо осветить, устранив отражения
11
(рефлексы) от поверхности роговицы и хрусталика. Чтобы устранить
рефлексы, для введения света используют одну часть зрачка, а другая служит
для получения изображения.
Несмотря на строго контролируемые условия, в которых производится
съемка глазного дна, на качество снимков также влияют несколько факторов,
зависящих от пациента и поэтому плохо поддающихся контролю. В
результате этого на изображениях глазного дна часто наблюдается
неоднородная яркость и контраст.
Таким образом, к факторам, влияющим на качество снимков глазного
дна, относятся:
- искривленная поверхность сетчатки. В результате этого все части
сетчатки не могут быть освещены одинаково;
- съемка требует расширенного зрачка. У каждого пациента уровень
расширенности зрачка свой. Недостаточный обзор глазного дна из-за
неполного расширения зрачка ведет к снижению качества снимка;
- непредвиденные движения глаз больного. Яркий свет заставляет
человека неосознанно изменять направление взгляда;
- наличие болезней, таких как катаракта, помутнение оптических сред
глаза и многих других, предотвращающие попадание света на сетчатку.
Для достижения достаточной освещенности при фотографировании
глазного дна во всех современных приборах используют газоразрядные
импульсные лампы с энергией вспышки 100-800 Дж. Малая длительность
вспышки таких ламп исключает влияние подвижности глаза и облегчает
процедуру съемки для пациента. Постоянный свет, необходимый в процессе
фокусировки аппарата, создают низковольтные лампы накаливания [2].
Однако, зачастую использование таких методов съемок глазного дна не
достаточно,
чтобы
получить
изображение
без
каких-либо
проблем,
влияющих на качество снимка, поэтому компьютерная предобработка
изображений оказывается очень актуальной в данном случае.
12
1.2 Анализ методов сегментации объектов на изображении
Сегментация объектов на изображении – это процесс разделения
изображения на несколько сегментов по тем или иным признакам. Целью
сегментации
является
упрощение
и/или
изменение
представления
изображения, чтобы его было легче анализировать. При сегментации
происходит процесс присвоения таких меток каждому пикселю (или группе
пикселей) изображения, что пиксели с одинаковыми метками имеют
одинаковые визуальные характеристики [7].
Чаще всего сегментация изображений применяется в следующих
областях знаний:
- Медицина;
- Выделение объектов на спутниковых снимках;
- Распознавание лиц;
- Распознавание отпечатков пальцев;
- Системы управления дорожным движением;
- Машинное зрение и др.
При разработке алгоритмов сегментации объектов на изображении в
большинстве случаев используют следующие методы.
- Методы, основанные на кластеризации:
Очень
эффективными
для
рационального
использования
вычислительных ресурсов являются методы кластерного анализа. Суть
кластеризации заключается в том, что все исходные объекты (пиксели)
разбиваются на несколько не пересекающихся групп таким образом, чтобы
объекты, оказавшиеся в одной группе, имели подобные характеристики, в то
время как у объектов из разных групп эти характеристики должны сильно
отличаться. Полученные группы называются кластерами.
Например, наиболее популярным методом кластеризации является
метод k-средних (k-means). Алгоритму многие отдают предпочтение из-за его
простоты реализации и высокой скорости вычислений [8].
13
Работа алгоритма заключается в том, что он стремится минимизировать
суммарное квадратичное отклонение точек кластеров от центров этих
кластеров. Метод k-средних – это итеративный алгоритм, который делит
данное множество пикселей на k кластеров, точки которых являются
максимально приближенными к их центрам, а сама кластеризация
происходит за счет смещения этих же центров.
k
V ( x j i ) 2 ,
(1.1)
i 1 x j Si
где k – число кластеров, S i – полученные кластеры, i 1,2,3,...k и i – центры
масс векторов x j S i .
Метод k-средних чувствителен к шуму, который может существенно
исказить результаты кластеризации. Перед кластеризацией необходимо
сделать предобработку изображения с целью уменьшения шумов.
Другим известным методом кластеризации является «сдвиг среднего»
(mean shift).
Метод mean shift группирует объекты с подобными признаками.
Пиксели с похожими признаками объединяются в один сегмент, на выходе
получается изображение с однородными областями.
При применении данного метода в качестве координат в пространстве
признаков рекомендуется выбрать координаты пикселя (x, y) и компоненты
RGB пикселя. Изобразив пиксели в пространстве признаков, можно заметить
сгущения в определенных местах.
При
выборе
в
качестве
признаков
координат
пикселей
и
интенсивностей по цветам в один сегмент будут объединяться пиксели с
похожими цветами и расположенные близко друг от друга. Соответственно,
если выбрать иной вектор признаков, то объединение пикселей в сегменты
уже будет происходить на основании другого указанного признака.
14
Например, если исключить из признаков координаты, то небо и море будут
считаться одним сегментом, потому что пиксели этих объектов в
пространстве признаков попали бы в один локальный максимум.
Если объект, который необходимо выделить, состоит из областей,
значительно различающихся по цвету, то mean shift не сможет объединить
эти области в одну, и объект будет состоять из нескольких сегментов. Но
зато данный метод позволяет выделить однородный по цвету предмет на
разноцветном фоне. Также mean shift применяют при реализации алгоритма
слежения за движущимися объектами [7,8].
- Пороговая обработка:
Пороговые преобразования занимают главное место в прикладных
задачах сегментации изображений из-за интуитивно понятных свойств и
простой реализации. Обработка состоит в сопоставлении значения яркости
каждого пикселя изображения с заданным значением порога. Выбор
подходящего
значения
пороговой
величины
позволяет
выделить
на
изображении области требуемого вида.
Операция порогового разделения, в результате которой получается
бинарное
изображение,
называется
бинаризацией.
Целью
операции
бинаризации является кардинальное уменьшение количества информации,
содержащейся
полутоновое
на
изображении.
изображение,
В
имеющее
процессе
бинаризации
несколько
уровней
исходное
яркости,
преобразуется в черно-белое изображение, пиксели которого имеют только
два значения – 0 и 1. Бинарное изображение проще обрабатывать, хранить и
пересылать, так как количество информации в бинарном изображении
обычно значительно меньше, чем в совпадающим с ним по размерам
полутоновом изображении [9].
Наиболее распространенными методами бинаризации являются:
1) Бинаризация с нижним порогом:
Бинаризация с нижним порогом является самой простой операцией, в
которой используется только одно значение порога:
15
0, f (m, n) t ;
f ' (m, n)
1, f (m, n) t.
(1.2)
2) Бинаризация с верхним порогом:
Операция бинаризации с верхним порогом:
0, f (m, n) t ;
f ' (m, n)
1, f (m, n) t.
(1.3)
3) Бинаризация с двойным ограничением:
Для выделения областей, в которых значения яркости пикселей может
меняться
в заданном
диапазоне,
вводится
бинаризация
с
двойным
ограничением (t1<t2):
0, f (m, n) t1 ;
f ' (m, n) 1, t1 f (m, n) t 2 ;
0, f (m, n) t
2.
(1.4)
4) Известны методы локальной пороговой обработки. К таким методам,
например, относится метод Отцу.
Метод Отцу применяется для вычисления порога бинаризации
полутонового изображения, в задачах компьютерного распознавания образов
и обработки изображений для получения черно-белых изображений. Метод
Отцу позволяет разделить множество пикселей на два класса, рассчитывая
такой порог, чтобы внутриклассовая дисперсия была минимальной. Данный
метод заключается в следующем [10]:
Разбиение на классы производится так, чтобы достигнуть минимума
дисперсии внутри классов, которая определяется как взвешенная сумма
дисперсий двух классов:
16
2 (t ) 1 (t ) 12 (t ) 2 (t ) 22 (t ),
(1.5)
где веса 1 (t ), 2 (t ) — это вероятности двух классов, разделенных порогом t ,
12 , 22 — дисперсии этих классов.
5) При решении задач бинаризации изображений широко используется
глобальная пороговая обработка, которая основана на анализе гистограммы
изображения.
Суть
глобальной
пороговой
обработки
заключается
в
разделении гистограммы изображения на две части с помощью единого
глобального порога.
Определение величины порога с помощью гистограммы яркостей –
простой метод, который позволяет добиться «чистой» сегментации, если
гистограмма изображения имеет четко выраженный бимодальный характер.
Такая форма гистограммы означает, что на изображении можно различить
яркие и темные пиксели. При этом гистограмма легко разделяется с
помощью одного глобального порога, расположенного во впадине между
пиками гистограммы.
К методам глобальной пороговой обработки относятся метод Бернсена,
Метод Эйквила, Метод Ниблэка, Метод Яновица и Брукштейна и др.
- Сегментация методом водораздела:
Понятие водораздела базируется на представлении изображения как
трехмерной поверхности, заданной двумя пространственными координатами
и уровнем яркости в качестве высоты поверхности (рельефа). В такой
«топографической» интерпретации рассматриваются точки трех видов: точки
локального минимума, точки, находящиеся на склоне, точки, находящиеся на
гребне или пике.
Главная цель алгоритма состоит в нахождении линий водораздела.
Основная идея метода выглядит следующим образом. Предположим, что в
каждом локальном минимуме проколото отверстие, после чего весь рельеф
наполняется водой, равномерно поступающей снизу через эти отверстия, так
17
что уровень воды везде одинаков. Когда поднимающаяся вода в двух
соседних бассейнах близка к тому, чтобы слиться вместе, в этом месте
ставится
перегородка,
препятствующая
слиянию.
В
конце
концов,
заполнение достигает фазы, когда над водой остаются видны только
верхушки
перегородок.
Эти
перегородки,
соответствующие
линиям
водоразделов, и образуют непрерывные границы, выделенные с помощью
алгоритма сегментации по водоразделам [11].
- Метод выращивания областей:
При
применении
данного
метода
на
исходном
изображении
первоначально выбираются точки (центры кристаллизации), ориентировочно
принадлежащие выделяемым областям, например, это могут быть точки с
максимальным уровнем яркости. Далее из этих точек начинается рост
областей, то есть присоединение к уже имеющимся точкам области соседних,
в ходе данного процесса используется определенный критерий их близости,
например, разница в яркости, заданная некоторой пороговой величиной.
Остановка роста областей производится по какому-либо условию, например,
максимальному отклонению яркости новых точек области от уровня яркости
центра кристаллизации или максимальной площади сегментов.
Недостатком данного метода является низкая степень автоматизации,
которая проявляется в необходимости определения центров кристаллизации
вручную или задания для их обнаружения жесткого правила, применимого
для крайне узкого класса изображений. Также выращивание каждой области
происходит, как в отдельном закрытом процессе без учета ситуации на иных
участках сегментируемой сцены. Еще одним недостатком данного метода
является жесткое, весьма не универсальное правило окончания выращивания
областей [9].
- Метод слияния областей:
В основу данного метода заложена идея о том, что пиксели исходного
изображения уже являются гомогенными областями, но при этом обладают
одинаковыми
минимальными
размерами.
18
В
данном
случае
метод
сегментации должен выполнять объединение соседних областей, наиболее
схожих по какому-либо параметру (по цвету или текстуре), пользуясь
функцией определения расстояния до тех пор, пока не будет выполнено (или
нарушено) некоторое заданное условие (например, размер сегментов или их
количество). Для данного метода полностью исчезает проблема определения
центров кристаллизации, но наиболее актуальной становится проблема
определения момента завершения процесса слияний [11].
- Методы теории графов:
Суть данных методов состоит в том, что изображение представляется в
виде
взвешенного
неориентированного
графа,
где
вершины
графа
ассоциируются с пикселем или группой пикселей изображения. Вес ребра
графа
определяет похожесть (или
не
похожесть)
точек.
Разбиение
изображения моделируется с помощью разрезов графа. Каждая часть вершин
(пикселей),
получаемая
этими
алгоритмами,
считается
объектом
на
изображении. Как правило, в методах теории графов вводится функционал
«стоимости» разреза, отражающий качество полученной сегментации. Так
задача
разбиения
изображения
на
однородные
области сводится
к
оптимизационной задаче поиска разреза минимальной стоимости на графе
[8]. Наиболее часто используемыми методами теории графов для решения
задачи
сегментации
изображений
являются
жадные
алгоритмы,
нормализованные разрезы графов, алгоритм Дейкстры, метод Nested Cuts,
метод сегментации SWA, минимальный разрез и др.
1.3 Особенности методов автоматической сегментации сосудов
глазного дна
Основными направлениями разработки методов выделения сосудов на
изображениях глазного дна являются, во-первых, сегментация сосудов,
19
позволяющая выделить всё дерево сосудов в одну итерацию, и, во-вторых,
прослеживание сосудов (трассировка).
При разработке алгоритмов сегментации сосудов в большинстве
случаев применяют следующие подходы.
- Классификация отсчётов изображения (без обучения, с обучением):
Алгоритмы классификации на диагностических изображениях. Данные
алгоритмы относятся к направлению распознавания образов и основаны на
детектировании
анатомических
частей
исследуемого
органа
или
классификации сосудов и иных объектов, включая фон. Они делятся на два
типа: методы с участием пользователя, использующие заранее указанную
информацию для принятия решения о принадлежности пикселя сосуду, и
полностью независимые методы, выполняющие сегментацию без участия
пользователя.
Методы с обучением. В данных методах используются правила,
выработанные алгоритмом на основе обучающей выборки из изображений,
размеченных офтальмологом вручную с достаточно высокой точностью.
Набор таких изображений называется золотым стандартом. Тем не менее,
имеются значительные различия в разметке даже среди специалистов.
Поскольку подобные методы основываются на уже классифицированных
изображениях, их эффективность, как правило, выше по сравнению с
другими методами.
Методы без обучения. Подходы, основанные на классификации
пикселей без использования обучения, выполняют поиск шаблонов,
характерных для кровеносных сосудов на изображениях сетчатки. В данных
алгоритмах не требуется непосредственного использования обучающей
выборки.
Применение различных классификаторов для выделения сосудистой
сети глазного дна показало свою эффективность на практике. Но данный
подход требует наличия эталонного набора изображений, на котором будет
производиться обучение системы [12].
20
- Согласованная фильтрация (matched filtering):
Согласованная фильтрация предполагает свертку изображения с
двумерным фильтром. Ядро фильтра моделирует определенные детали
изображения, отклик фильтра показывает присутствие или отсутствие
искомой детали. Ядро фильтра строится на основе трех основных
предположений:
- сосуды могут быть аппроксимированы ломаными;
- диаметр сосуда уменьшается при движении от оптического диска;
- профиль поперечного сечения сосуда похож на гауссову кривую.
Ядро фильтра следует брать довольно большое, а сам фильтр
необходимо применять в нескольких направлениях, что приводит к
существенному росту сложности вычислений. Ядро фильтра должно также
соответствовать масштабу сосуда, который требуется обнаружить. Как
следствие, фильтр не даст четкого отклика на сосуды с профилем, имеющим
другой масштаб. Изменение фона изображения и наличие различных
патологических
артефактов
также
увеличивает
процент
ложных
срабатываний фильтра, поскольку патологии могут иметь параметры,
близкие к параметрам сосудов. Граница оптического диска, кровоизлияния,
повреждения тканей также могут иметь характеристики, свойственные
сосудам. Использование протяженного структурного элемента также может
затруднить
выделение
крайне
извилистых
сосудов.
Применение
согласованной фильтрации является оправданным в комбинации с другими
методами обработки изображений [13].
- Морфологическая обработка:
Морфологическая обработка подразумевает, что сосуды состоят из
соединенных между собой линейных сегментов. Морфологические операции
используются для сегментации сосудов и выделения микроаневризмов [1416].
При
сегментации
медицинских
изображений
применяются
в
большинстве случаев две морфологические операции: top-hat преобразование
21
[17] и метод водораздела [18]. Подчёркивание сосудов осуществляется с
помощью top-hat преобразования (его эффект связан с морфологическими
операциями замыкания и размыкания), которое оценивает локальный фон, а
затем вычитает его из исходного изображения, выделяя таким образом
сосуды. Заранее известно, что базовую структуру сосудистой системы можно
представить с помощью соединённых между собой линейных сегментов.
Преимуществами морфологической обработки для идентификации
определённых форм являются высокая скорость работы и устойчивость к
шуму.
Главным недостатком использования исключительно морфологических
методов является то, что они не используют информацию о форме профиля
сосудов. Кроме того, при осуществлении поиска протяженных структур, не
учитывается возможность наличия на изображении извилистых сосудов [12].
- Многомасштабный анализ:
Многочисленные исследования [19] показали, что поперечный профиль
сосуда возможно аппроксимировать функцией Гаусса, в локальном масштабе
он представляет собой прямолинейный объект с плавно уменьшающимся
диаметром. Идея, лежащая в основе использования пространственно
масштабного анализа для извлечения сосудов, заключается в определении не
зависящей от масштаба информации о сосудах [20]. В процессе разработки
фильтра для выделения сосудов имеется возможность рассматривать
многомасштабные структуры второго порядка (матрица Гессе). Алгоритм
поиска сосудов базируется на анализе собственных чисел и векторов
матрицы Гессе. Он определяет основные направления, по которым локальная
структура может быть разложена на составляющие, дающие направления
наименьшей кривизны поля вдоль русла сосуда. Результирующая оценка
сосуда вычисляется с использованием геометрических параметров сосуда, а
также собственных значений и норм матрицы Фробениуса, и проверяется на
изображениях,
полученных
магниторезонансной
с
помощью
томографии.
Подход
22
ангиографии
и
трёхмерной
многомасштабного
анализа
показал способность одновременного подавления шума и фона изображения.
На данном подходе базируются многие алгоритмы.
- Подходы, основанные на моделях:
Данные подходы используют явные модели сосудов для извлечения
структуры сосудистой системы [21]. Подходы, основанные на моделях,
делятся на два типа: (1) модели профиля сосуда; (2) деформируемые модели.
1) Модели профиля сосуда. В данном подходе поперечный профиль
сосуда аппроксимируется с помощью кривой Гаусса или комбинации
кривых, если присутствует центральный рефлекс. Также могут применяться
вторая производная функции Гаусса, кубические сплайны или полиномы
Эрмита. В более сложные модели включаются яркие или тёмные
повреждения сетчатки и характеристики фона для того, чтобы увеличить
точность сегментации на сложных изображениях. В некоторых моделях
также учитывается равномерный фон. Учет ветвления и пересечений сосудов
ведет к последующему усложнению модели.
2) Деформируемые модели. Методики сегментации сосудов на основе
деформируемых
моделей
подразделяются
на
два
типа
[22]:
(а)
параметрические модели; (б) геометрические деформируемые модели.
а)
Параметрические
модели.
Одной
из
наиболее
известных
параметрических моделей является активный контур. Метод активных
контуров, также известный как снейки [23], основан на изменении формы
контура вокруг объекта под действием внутренних и внешних сил. Внешние
и внутренние силы устроены таким образом, что под их воздействием контур
сжимается до границ объекта. Преимуществом метода по сравнению с
другими подходами является способность адаптироваться к объектам
различной формы и находить состояния с минимальной энергией, в которых
внутренние и внешние силы уравновешены. Метод активных контуров может
применяться для отслеживания объектов, как в пространстве, так и во
времени. Основным ограничением данного метода является использование
информации только о границах объекта, без учёта других параметров
23
изображения. Для корректной работы необходимо задать первоначальное
положение контура достаточно близко к объекту. Метод активных контуров
часто используется в таких областях, как трассировка, распознавание формы
объекта, сегментация и определение границ объекта.
б) Геометрические модели. Геометрическая модель активных контуров
базируется на теории эволюции геометрических форм. Данные модели, как
правило, реализуются с использованием численных алгоритмов [22]. Одна из
реализаций таких моделей, метод LSM (Level Set Method), является
численным алгоритмом для трассировки определённых поверхностей и форм
[24]. Преимуществом LSM метода является то, что численная обработка
различных кривых и поверхностей может быть проведена без указания
параметров этих объектов.
- Прослеживание сосудов (vessel tracking):
Алгоритмы прослеживания сосудов (трассировка) сегментируют сосуд
между двумя заданными точками. Данный метод работает на уровне одного
сосуда, а не на уровне сосудистой сети [25].
Преимуществом
данного
подхода
является
возможность
прослеживания тонких и разрывных сосудов, которые не обнаруживаются
при использовании иных методов. Также к преимуществу прослеживания
сосудов относится то, что алгоритмы трассировки требуют меньшего
количества вычислений. Учитывая, что все сосуды связаны между собой в
общей сосудистой системе, алгоритм трассировки может выделить все
дерево целиком без потери времени на обработку тех частей изображения,
которые не содержат сосуды. Таким образом, в отличие от других подходов
трассировка способна дать информацию о структуре сосудистого дерева: о
бифуркациях и соединениях отдельных ветвей.
Однако данный подход требует знание начальных точек, и иногда
конечных, принадлежащих детектируемому сосуду, и это затрудняет
применение данного подхода в полностью автоматических системах. Также
24
дополнительную трудность вызывают ситуации пересечения и разветвления
сосудов.
1.4 Постановка задачи
На основе проведенного анализа существующих методов сегментации
кровеносных сосудов на изображении глазного дна можно утверждать, что
данные методы имеют отдельные недостатки, и, зачастую, в них не
проводится
предобработка
изображения
для
устранения
различных
артефактов, снижающих качество снимка. Также большинство методов
необходимо применять в комбинации с другими методами для достижения
наилучшего эффекта сегментации, следует отметить, что в большинстве
случае выбор комбинации методов не является тривиальным.
Таким образом, целью работы является совершенствование методов
сегментации кровеносных сосудов на изображении глазного дна.
Для достижения поставленной цели необходимо решить следующие
задачи:
1. Анализ методов сегментации сосудов глазного дна;
2.
Проведение
исследования
теоретических
основ
методов,
применяемых для решения задачи сегментации сосудов глазного дна;
3. Разработка метода и алгоритма сегментации кровеносных сосудов на
изображении глазного дна;
4.
Разработка
программной
реализации
алгоритма
сегментации
кровеносных сосудов на изображении глазного дна;
5.
Проведение
вычислительных
работоспособности разработанного метода.
25
экспериментов
по
проверке
2. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ МЕТОДОВ, ПРИМЕНЯЕМЫХ ДЛЯ
РЕШЕНИЯ ЗАДАЧИ СЕГМЕНТАЦИИ СОСУДОВ ГЛАЗНОГО ДНА
В данной главе рассматриваются теоретические основы методов,
применяемых для решения задачи сегментации сосудов. К таким методам
относятся контрастно ограниченная адаптивная эквализация гистограммы,
морфологическая фильтрация и метод кластеризации k-means, для которых
сформулированы основные определения, вычислительные соотношения и
алгоритмы. В главе также выполнен анализ преимуществ и недостатков
данных методов. Процесс последовательного применения эквализации
гистограммы, морфологических операций и кластеризации показал свою
эффективность при решении задачи сегментации сосудов глазного дна.
2.1 Метод повышения контрастности изображений
Эквализация гистограммы – это метод обработки изображений по
улучшению контрастности с использованием гистограммы изображения [11].
Целью эквализации гистограммы является равномерное распределение
уровней серого в пределах изображения.
увеличивает
яркость
и
контрастность
Эквализация
темных
и
гистограммы
низко-контрастных
изображений. Данный подход позволяет выделять на изображении такие
детали, которые не были видны на исходном изображении. Он также
используется для нормирования яркости и контрастности изображений.
Также эквализация гистограммы используется в биологических
нейронных сетях, чтобы максимизировать скорость возбуждения нейтронов в
зависимости от исходных данных. Это было доказано, в частности, при
работе с сетчаткой мухи.
Эквализация гистограммы является отдельным случаем более общего
класса методов перераспределения гистограммы.
26
Данный метод направлен на корректировку изображения, чтобы
упростить анализ или улучшить качество изображения.
Процесс эквализации гистограммы заключается в поиске такой
функции, которая преобразует исходную гистограмму изображения в
равномерно
распределенную
выходную
гистограмму.
Эквализация
гистограммы обычно увеличивает глобальный контраст изображений.
Благодаря этой адаптации становится возможным лучше распределить
интенсивности на гистограмме. Это позволяет области с более низким
локальным
контрастом
преобразовать
в
более
высокий
контраст.
Эквализация гистограммы позволяет эффективно распределять наиболее
часто используемые значения интенсивности [26].
Данный метод полезен для изображений с одновременно яркими или
темными задним и передним планом. В частности, метод эквализации
гистограммы может привести к лучшему представлению структуры
кровеносной системы на изображениях или сильнее детализировать
фотографии с плохим освещением. Ключевым преимуществом метода
является то, что это довольно простой метод, и он является обратимым.
Теоретически, если известна функция эквализации гистограммы, то исходная
гистограмма может быть восстановлена. Объем вычислений не является
значительным с точки зрения вычислительных затрат.
Недостатком метода является то, что он является неизбирательным. Он
может увеличить яркость фонового шума и при этом уменьшить
контрастность
изображениях,
необходимого
для
которых
объекта,
что
необходимо
пространственная
учитывать
корреляция
в
важнее
интенсивности сигнала (например, изображения разделения фрагментов ДНК
дискретной длины), и для которых малые отношения сигнал/шум обычно
препятствуют визуальному анализу.
Эквализация гистограммы часто создает нереалистичные эффекты на
фотографиях, однако она очень полезна для изображений, используемых при
научных исследований, таких как тепловые, спутниковые или рентгеновские
27
изображения. Кроме того, эквализация гистограммы может создавать
нежелательные эффекты (например, градиент видимого изображения) при
применении к изображениям с низкой глубиной цвета. Например, если
эквализацию применить к 8-битовому изображению, отображаемому с 8битной палитрой серого, это еще больше уменьшит глубину цвета
(количество уникальных оттенков серого) изображения. Эквализацию
гистограммы следует применять к изображениям с гораздо более высокой
глубиной цвета, чем размер палитры, например, 16-битные изображения в
оттенках серого [27].
Существует два способа реализовать эквализацию гистограммы путем
изменения изображения или путем изменения палитры. Операция может
быть выражена как P(M(I)), где I – исходное изображение, M – операция
отображения эквализации гистограммы, а P – палитра. Если мы определим
новую палитру как P'=P(M) и оставим изображение неизменным, то
выравнивание гистограммы реализуется как изменение палитры. С другой
стороны, если палитра P остается неизменной, а изображение изменено на
I'=M(I), то реализация осуществляется путем изменения изображения. В
большинстве случаев изменение палитры лучше, так как оно сохраняет
исходные данные.
Модификации данного метода используют несколько гистограмм,
называемых субгистограммами, чтобы подчеркнуть локальный контраст, а не
общий
контраст.
Примерами
таких
методов
являются
адаптивная
эквализация гистограммы, контрастно ограниченная адаптивная эквализация
гистограммы (CLAHE), многовершинная эквализация гистограммы (MPHE)
и
многоцелевая
(MBOBHE).
Цель
бета-оптимизированная
этих
методов
состоит
эквализация
в
том,
бигистограммы
чтобы
улучшить
контрастность, не создавая артефакты сдвига среднего значения, и
детализировать утраченные признаки, изменяя алгоритм эквализации
гистограммы.
28
Процедура эквализации гистограммы основана на предположении, что
качество изображения одинаково во всех его областях, и одно уникальное
отображение
оттенков
серого
произведет
одинаковые
эффекты
по
улучшению для всех областей изображения. Однако, когда распределение
оттенков серого меняется от одной области на изображении к другой, это
предположение не является корректным. В этом случае техника адаптивной
эквализации гистограммы может значительно превзойти первоначальный
(стандартный) подход.
Адаптивная эквализация гистограммы (AHE) улучшает стандартную
эквализацию гистограммы путем преобразования каждого пикселя с
помощью функции преобразования, полученной на основе значений
пикселей из его окрестности. Данный метод был впервые разработан для
использования в дисплеях кабины самолетов [28].
Основная идея адаптивной эквализации гистограммы заключается в
нахождении отображения для каждого пикселя на основе его локального
распределения в оттенках серого. Отображение улучшений контрастности,
применяемое
интенсивности
к конкретному пикселю,
непосредственно
является
окружающих
функцией значений
пиксель.
Количество
вычислений, которые должны быть проведены, совпадает с количеством
пикселей на изображении. Это требует значительные вычислительные
затраты,
которые
даже
с
некоторой
модификацией,
не
позволяют
использовать данный подход для улучшения изображения в режиме
реального времени.
Метод,
который
является
компромиссом
между
глобальной
эквализацией гистограммы и полностью адаптивным алгоритмом – это
эквализация гистограммы отдельных областей. В этом случае изображение
делится на ограниченное число областей, и метод эквализации гистограммы
применяется к пикселям в каждой области. Этот метод, который можно
оптимизировать для обработка в реальном времени, имеет существенный
недостаток. Проблема заключается в том, что из-за различий между
29
отображениями двух соседних областей пиксели по обе стороны границы
будут отображаться по-разному и, следовательно, в результате будет получен
значительный эффект представления изображения в виде набора блоков [29].
Альтернативой
полностью
адаптивному
алгоритму
является
приближение требуемого отображения для каждого пикселя, выбирая
небольшое количество необходимых областей внутри изображения и
вычисляя
отображение
интерполяцию
для
отображений,
заданного
полученных
пикселя
из
как
билинейную
близлежащих
областей.
Другими словами, отображение для пикселя получается с использованием
взвешенной суммы отображения
его четырех ближайших областей.
Вычисленные веса основываются на близости пикселя к центрам четырех
ближайших областей.
В некоторых случаях, когда распределение оттенков серого сильно
локализовано,
нежелательно
преобразовывать
низкоконтрастные
изображения путем полной эквализации гистограммы. В этих случаях кривая
отображения может включать в себя монотонно изменяющиеся сегменты,
что означает, что два очень близких уровня серого могут быть сопоставлены
с существенно отличающимися уровнями серого. Эта проблема решается
путем ограничения контраста, изменение которого разрешено посредством
эквализации гистограммы. Сочетание данного подхода с ограничением
контрастности с результатами вышеупомянутой адаптивной эквализацией
гистограммы называется контрастно ограниченной адаптивной эквализацией
гистограммы (Contrast Limited Adaptive Histogram Equalization – CLAHE).
Contrast Limited AHE (CLAHE) отличается от обычной адаптивной
эквализации
гистограммы
ограничением
констрастности.
Данная
особенность также может быть применена к глобальной эквализации
гистограммы, которая редко используется на практике. В случае CLAHE
процедура
ограничения
контраста
должна
применяться
для
каждой
окрестности, из которой выведена функция преобразования. Метод CLAHE
30
был разработан, чтобы предотвратить избыточное усиление шума, которое
может стать результатом адаптивной эквализации гистограммы [30].
Это достигается путем ограничения
Контрастное
усиление
в
окрестности
усиления контраста AHE.
заданного
значения
пикселя
определяется наклоном функции преобразования. Оно пропорционально
наклону функции кумулятивного распределения окрестности (CDF) и,
следовательно, значению гистограммы при данном значении пикселя.
CLAHE ограничивает усиление путем отсечения гистограммы с заранее
заданным значением перед вычислением CDF. Это ограничивает наклон CDF
и,
следовательно,
функции
преобразования.
Значение,
при
котором
гистограмма обрезается, называемое предел урезания (clip limit), зависит от
нормализации гистограммы и, следовательно, от размера окрестности.
Общие значения ограничивают полученное усиление от 3 до 4 раз.
Данный метод формулируется на основе деления изображения на
несколько неперекрывающихся областей почти равных размеров. Для
изображений размером 512×512 для достижения хорошей статистической
оценки количество областей обычно выбирается равным 64, таким же как
при разделении изображения на 8 областей в каждом направлении.
Предпочтительно
не
отбрасывать
часть
гистограммы,
которая
превышает предел урезания, а распределять ее поровну между всеми
ячейками гистограммы (рис. 2.1).
Рис. 2.1. Применение контрастно ограниченной адаптивной эквализации
гистограммы
31
Перераспределение снова поместит несколько значений над пределом
урезания (на рис. 2.1 данная область выделена серым цветом), в результате
чего эффективный предел урезания больше заданного предела, и точное
значение предела урезания зависит от изображения. Если это нежелательно,
процедура перераспределения может повторяться рекурсивно, пока избыток
не будет незначительным [31].
Контрастно
ограниченная
адаптивная
эквализация
гистограммы
показала хорошие результаты при обработке медицинских изображений.
2.2 Основные операции морфологической фильтрации
изображения
Математическая морфология – теория и техника анализа и обработки
геометрических структур, основанная на теории множеств, топологии и
случайных функциях. В основном применяется в обработке цифровых
изображений, но также может применяться на графах, полигональной сетке,
стереометрии и многих других пространственных структурах [11].
Основными
базовыми
морфологическими
операциями
являются
дилатация и эрозия.
Дилатация. Операция дилатации «наращивает» или «утолщает»
объекты на двоичных изображениях. Способ и степень этого утолщения
контролируется
некоторой
формой,
которая
называется
структурообразующим элементом.
Дилатацию можно строго определить математически в терминах
операций над множествами. Дилатация множества А по множеству В (или
относительно множества В) обозначается А В и определяется по формуле:
32
^
A B {z | ( B) z A } ,
(2.1)
где обозначает пустое множество, а В – это структурообразующий
элемент.
Другими словами,
множеством,
которое
дилатация
состоит
из
А относительно В является
всех
положений
центров
структурообразующего элемента, который, будучи центрально отраженным и
сдвинутым, имеет частичное перекрытие с множеством А. Процесс сдвига
структурообразующего элемента похож на механизм построения свертки.
Операция дилатации является коммутативной, т.е. А В=В А. Однако
бывает удобно помещать на место первого операнда изображение, а на место
второго – структурообразующий элемент, который, обычно, существенно
меньше обрабатываемого изображения [32].
На рис. 2.2 показан простой пример дилатации.
а
б
Рис. 2.2. Применение дилатации: а – исходное изображение, содержащее
разорванные символы, б – изображение после применения дилатации
Эрозия. Операция эрозия «ужимает» или «утончает» объекты двоичных
изображений. Как и при дилатации, вид и размер эрозии определяется
структурообразующим элементом.
Математическое определение операции эрозии похоже на определение
дилатации. Эрозией множества А по В называется множество:
33
AB {z | ( B) z A} .
(2.2)
Другими словами, эрозия А по примитиву В – это множество все таких
точек , при сдвиге в которые множество В целиком содержится в А.
На рис. 2.3 приведен пример использования эрозии.
Рис. 2.3. Иллюстрация эрозии: а – с кругом радиуса 5, б – с кругом радиуса
20, в – с кругом радиуса 10, г – с кругом радиуса 30
В
конкретных
приложениях
обработки
изображений
операции
дилатация и эрозия часто используется совместно. Обычно изображения
подвергаются серии операций дилатации и/или эрозии с одним и тем же, а
иногда и с разными структурообразующими элементами [33].
Так, важнейшими морфологическими операциями, которые сочетают
дилатацию и эрозию, являются размыкание и замыкание.
Морфологическое размыкание А по В обозначается А◦В и определяется
как эрозия А по В, после которой выполняется дилатация результата по В:
А◦В=(A B) B.
Эту операцию можно также задать эквивалентной формулой
34
(2.3)
А◦В= {(B)z|(B)z A},
(2.4)
где {.}обозначает объединение всех множеств внутри фигурных скобок, а
формула C D означает, что множество С является подмножеством D. Такое
представление имеет простую геометрическую интерпретацию. На рис. 2.4а)
показано множество А и структурообразующий элемент В, имеющий форму
круга. На рис. 2.4б) изображены некоторые сдвиги множества В, которые
полностью содержатся в А. Объединение всех таких сдвигов выделено
темным цветом на рис. 2.4в). Эта область и является размыканием А по В. На
этом рисунке белым цветом показаны области, где структурообразующий
элемент целиком не поместился, поэтому они не принадлежат размыканию.
Морфологическое размыкание удаляет те части объектов, в которых
структурообразующий элемент полностью не помещается. Оно также
сглаживает контуры объектов, разрушает тонкие соединения и удаляет
острые выступы [34].
Морфологическое замыкание множества А по В обозначается А•В. Эта
операция представляет собой эрозию, примененную к результату дилатации:
А•В=(A B) B.
(2.5)
Геометрически множество А•В представляет собой дополнение к
объединению всех сдвигов множеств В, которые не пересекаются с А. На
рис. 2.4г) построено несколько сдвигов множества В, которые не
пересекаются с А. Взяв дополнение к объединению всех таких сдвигов, мы
получим область, затемненную на рис. 2.4г), которое есть искомое
замыкание. Подобно размыканию, морфологическое замыкание приводит к
сглаживанию контуров объектов. Однако в отличие от размыкания, оно
соединяет узкие разрывы, «заливает» длинные углубления малой ширины и
заполняет
малые
отверстия,
диаметр
структурообразующего элемента.
35
которых
меньше
размеров
Рис. 2.4. Размыкание и замыкание как объединение сдвигов
структурообразующих элементов: а – множество А и структурообразующий
элемент В, б – сдвиги В, которые полностью помещаются внутри А,
в – полное размыкание А (затемненная область), г – сдвиги В, которые лежат
целиком вне А, д – полное замыкание А (затемненная область)
Также морфологические операции замыкания и размыкания могут
применяться
комбинированно.
Такое
их
использование
приводит
к
различным эффектам на изображение, что может быть необходимо при
обработке изображений. На рис. 2.5 приведен пример последовательного
использования операций размыкания и замыкания.
Рис. 2.5. Иллюстрация размыкания и замыкания: а – исходное изображение,
б – размыкание, в – замыкание, г – замыкание изображения б
В приложениях часто бывает необходимо
уметь распознавать
определенные конфигурации пикселов, например, изолированные пикселы
переднего плана или пикселы, которые являются концами сегментов линий.
Преобразования успех/неудача помогают при решении подобных задач.
36
Преобразование успех/неудача множества А по множеству В обозначается
через А В. Здесь В обозначает структурообразующую пару элементов В =
(В1, В2), а не один элемент, как это было раньше. По определению,
преобразованием множества А по В называется множество:
A B=(A В1)∩(Ac В2)
(2.6)
Название операции «успех/ неудача» происходит от результата двух
эрозий. Например, выходное изображение на рис. 2.6 состоит из всех
позиций пикселов, которые подходят под В1 («успех») и полностью не
подходят под В2 («неудача»).
а
б
Рис. 2.6. Применение операции успех/неудача: а – исходное изображение,
б – результат применения преобразования успех/ неудача (отмеченные точки
были увеличены для лучшей заметности)
На практике широко применяется морфологическая фильтрация.
Одним из наиболее распространенных способов данной фильтрации при
обработке сосудов является top-hat преобразование [35].
В математической морфологии и цифровой обработке изображений
top-hat преобразование – это операция, которая извлекает небольшие
элементы и детали из заданных изображений. Существует два типа top-hat
преобразования:
37
- белое top-hat преобразование определяется как разница между
входным
изображением
и
его
размыканием
некоторым
структурообразующим элементом;
- черное top-hat преобразование определяется как разность между
замыканием и входным изображением.
Top-hat преобразования используются для различных задач обработки
изображений,
таких
как
извлечение
объектов,
выравнивание
фона,
улучшение изображения и другое [36].
Дадим математическое описание top-hat преобразования.
Пусть f : E R – изображение в градациях серого, где происходит
отображение точек из Евклидова пространства или дискретной сетки E в
вещественную числовую прямую. И пусть b(x ) – структурообразующий
элемент в градациях серого.
Тогда белое top-hat преобразование вычисляется по формуле:
Tw ( f ) f f b ,
(2.7)
где – морфологическая операция размыкания.
Черное top-hat преобразование вычисляется по формуле:
Tb ( f ) f b f ,
(2.8)
где – морфологическая операция замыкания.
Белое top-hat преобразование возвращает изображение, содержащее эти
«объекты» или «элементы» входного изображения, которые:
-«Меньше», чем структурообразующий элемент (то есть участки
изображения, в которых не помещается структурообразующий элемент);
- Ярче их окрестностей.
38
Черное top-hat преобразование возвращает изображение, содержащее
«объекты» или «элементы», которые:
- «Меньше», чем структурообразующий элемент;
- Темнее их окрестностей.
Размер
или
преобразование,
ширина
можно
структурообразующего
элементов,
которые
контролировать
b
элемента
.
извлекаются
top-hat
с
помощью
выбора
Чем
больше
размеры
структурообразующего элемента, тем большие по размеру извлекаются
элементы.
Оба top-hat преобразования представляют собой изображения, которые
содержат только неотрицательные значения во всех пикселях [37].
Одним из важных применений данного преобразования в сегментации
изображений
является
корректировка
неравномерного
освещения
на
изображении и обеспечение лучшего порогового значения для разделения
объектов.
2.3 Метод кластеризации k-means
Метод k-средних (англ. k-means) — наиболее популярный метод
кластеризации.
Штейнгаузом
Был
и
изобретён
почти
в
1950-х
одновременно
годах
Стюартом
математиком
Ллойдом.
Гуго
Особую
популярность данный метод получил после работы Маккуина [8].
При данном множестве объектов, цель кластеризации состоит в том,
чтобы разделить эти объекты на группы или «кластеры» такие, что объекты
внутри группы имеют тенденцию быть более похожими друг на друга по
сравнению с объектами, принадлежащими различным группам. Иными
словами, алгоритмы объединения в кластеры размещают похожие точки в тот
же кластер, в то время как непохожие точки помещают в различные
кластеры. В отличие от контролируемых задач, таких как регрессия или
39
классификация, где есть понятие целевого значения или метка класса,
объекты, которые образуют входы в процедуру кластеризации, не идут со
связанной целью. Следовательно, кластеризация часто упоминается как
бесконтрольное обучение. Поскольку нет никакой необходимости в
маркированных данных,
алгоритмы обучения
без
учителя
являются
подходящими для многих приложений, где маркированные данные трудно
получить. Бесконтрольные задачи, такие как кластеризация, так же часто
используются, чтобы исследовать и охарактеризовать набор данных перед
запуском контролируемой задачи обучения. Поскольку кластеризация
производится, не используя метки класса, некоторое представление о
сходстве должно быть определено на основании атрибутов объектов.
Описание сходства и метода, в котором точки сгруппированы, отличаются на
основании
применяемого
алгоритма
кластеризации.
Таким
образом,
различные алгоритмы кластеризации подходят для различных типов набора
данных и различных целей. Выбор «лучшего» алгоритма кластеризации
зависит от решаемой задачи [38].
Алгоритм k-means является простым повторяющимся алгоритмом
кластеризации, который разделяет определенный набор данных на заданное
пользователем число кластеров k. Алгоритм прост для реализации и запуска,
относительно быстрый, легко адаптируется и распространен на практике. Это
исторически один из самых важных алгоритмов интеллектуального анализа
данных.
Алгоритм k-means применяется к объектам, которые представляются
точками в d-мерном векторном пространстве R d . Таким образом, это
кластеры набора d-мерных векторов:
D { xi }, i 1,..., N ,
где xi R d обозначает i-ый объект или «точку данных».
40
(2.9)
k-means является алгоритмом кластеризации, который разделяет D на k
кластеров точек. То есть, алгоритм k-means объединяет все точки данных в D
так, что каждая точка xi попадает в один и только один из k разделов.
Можно отследить, какая точка находится в каком кластере, назначив каждой
точке номер кластера. Точки с таким же номером кластера находятся в одном
и том же кластере, в то время как точки с различными номерами кластера
находятся в разных кластерах. Это можно обозначить как кластерный
составной вектор m длинною N, где mi является номером кластера xi .
Значение k является основным из входных данных алгоритма. Как
правило, значение k основано на критериях, таких как предварительное
знание о том, сколько на самом деле кластеров появится в D, как много
кластеров требуется для текущего приложения, или типы кластеров,
найденные путем изучения/экспериментирования с различными значениями
k [39].
В k-means каждый из k кластеров представлен одной точкой в R d .
Обозначим
этот
C {c j }, j 1,..., k .
набор
представителей
кластера
как
множество
Эти представители k кластеров также называются
состоянием кластера или центройдами кластера.
В алгоритмах кластеризации точки группируются некоторым понятием
«близости» или «подобия». В k-means мера близости по умолчанию –
Евклидово расстояние. В частности можно с готовностью показать, что kmeans минимизирует следующую неотрицательную функцию стоимости:
N
Cost (arg min
i 1
j
xi c j
2
2
)
(2.10)
Другими словами, k-means минимизирует итоговый квадрат Евклидова
расстояния между каждой точкой xi и ее самым близким представителем
кластера c j . Уравнение (2.10) также называют целевой функцией k-means.
41
Алгоритм k-means разделяет D итеративным способом, чередующимся
между 2 шагами:
1. Переприсваивание номера кластера всем точкам в D;
2. Обновление представителей кластера, основанных на точках данных
в каждом кластере.
Алгоритм работает следующим образом. Сначала представители
кластера инициализируются, выбирая k из R d . Методы для выбора этих
начальных источников включают случайную выборку из набора данных,
устанавливая их как решение кластеризации маленького подмножества
данных, или нарушая глобальное среднее значение данных k раз. В
алгоритме происходит инициализация случайно выбранных k точек. Затем
выполняется алгоритм итерации до сходимости за 2 шага.
Шаг 1: присваивание данных. Каждой точке данных присваивается ее
самый близкий представитель притом, что связи нарушаются произвольно.
Это приводит к разделению данных.
Шаг 2: перемещение «средних». Каждый представитель кластера
перемещается к центру (т. е. среднее арифметическое) всех точек данных,
присвоенных ему. Объяснение этого шага основано на наблюдении, что
данное множество точек, единственный лучший представитель для этого
множества (в смысле минимизации суммы квадрата Евклидова расстояния
между каждой точкой и представителем) не что иное, как срединная точка
данных. Именно поэтому представитель кластера часто взаимозаменяемо
называют срединным элементом кластера или центроидом кластера, отсюда
и название этого алгоритма [40].
Алгоритм сходится, когда присвоение (следовательно, и значения с j )
больше не изменяются. Можно показать, что целевая функция k-means,
определенная в уравнении (2.10), будет уменьшаться всякий раз, когда есть
изменения в присвоении или шагах измерения, и сближение (сходимость в
одной точке) гарантировано за конечное число итераций.
42
Выбор оптимального значения k является трудной задачей. Если чтолибо известно о наборе данных, например, число разделов, которые
естественно включают набор данных, тогда эти значения могут быть
использованы при выборе k. Иначе, необходимо использовать другие
критерии выбора k, таким образом, решая проблему выбора модели. Самое
простое решение – попробовать несколько различных значений k и выбрать
кластеризацию, которая минимизирует целевую функцию k-means. К
сожалению, значение целевой функции не столь же информативно, как
можно было бы надеяться в этом случае. Для примера, стоимость
оптимального решения сокращается с увеличением k до тех пор, пока она не
достигнет нуля, когда число кластеров равняется числу отличных точек
данных. Это делает его более трудным для использования объективной
функции, непосредственно сравнивая решения с различным числом
кластеров, и нахождения оптимального значения k. Таким образом, если
требуемое k не будет известно заранее, как правило, k-means будет запущен с
различными значениями k, а затем использовать другие, более подходящие
критерии для выбора одного из результатов. В качестве альтернативы, можно
прогрессивно увеличивать число кластеров в соединении с подходящим
критерием остановки [41].
Алгоритм k-means заключается в следующем:
Этап 1. Первоначальное распределение объектов по кластерам:
- Выбор случайным образом k точек данных из D как начальное
множество представителей кластера C;
- Распределение объектов по кластерам в соответствие с формулой
(2.10).
Этап 2. Перераспределение срединных элементов:
- Вычисление центра для каждого кластера;
- Перераспределение объектов по кластерам.
Из алгоритма видно, что каждая итерация нуждается в N*k сравнений,
которые определяют временную сложность одной итерации. Число итераций,
43
требуемых
для
сходимости,
изменяется
и
может
зависеть
от
N.
Соответственно, чем больше точек во множестве N, тем дольше будет
работать алгоритм.
Сократить время работы алгоритма можно путем распараллеливания
этапа распределения точек по кластерам. Если у нас имеется Р процессоров,
то мы можем создать Р потоков. В этом случае исходное множество данных
можно разбить на Р частей и каждый поток будет работать только со своим
объемом информации, независимо от остальных данных.
На рис. 2.7 приведен пример действия алгоритма k-means в двумерном
случае. Начальные точки выбраны случайно.
а
б
в
г
Рис. 2.7. Результат работы алгоритма k-means:
а – исходные точки и случайно выбранные начальные точки, б – точки,
отнесённые к начальным центрам, разбиение на плоскости – диаграмма
Вороного относительно начальных центров, в – вычисление новых центров
кластеров (ищется центр масс); г – предыдущие шаги повторяются, пока
алгоритм не сойдётся
Потенциальная проблема алгоритма – проблема «пустых» кластеров.
При запуске k-means, особенно с большим значением k и/или когда данные
находятся в большом размерном пространстве, возможно, что в какой-то
момент исполнения, существует представитель кластера с j , такой, что все
точки xi в D ближе к некоторому другому представителю кластера, который
не является с j . Когда точки в D будут присвоены к ближайшему кластеру, j44
му кластеру будут присвоены нулевые точки. Таким образом, кластер j
является теперь пустым кластером. Стандартный алгоритм не принимает
меры
против
пустых
кластеров,
но
простые
решения
(такие
как
переинициализация представителя пустого кластера или «кража» некоторых
точек из большого кластера) возможны [42].
Другие проблемы алгоритма k-means заключаются в следующем:
1. Не гарантируется достижение глобального минимума суммарного
квадратичного отклонения V, а только одного из локальных минимумов;
2. Результат зависит от выбора исходных центров кластеров, их
оптимальный выбор неизвестен;
3. Число кластеров надо знать заранее.
Несмотря на все недостатки, k-means остается наиболее широко
используемым разделяющим алгоритмом кластеризации на практике.
Алгоритм простой, понятный и достаточно масштабируемый и может быть
легко модифицирован для решения различных задач, таких как частичное
обучение с учителем или потоковых данных. Постоянные улучшения и
обобщения основных алгоритмов обеспечили ее актуальность и постепенно
увеличивают его эффективность.
В данной главе были рассмотрены теоретические положения методов
контрастно
ограниченной
адаптивной
эквализации
гистограммы,
морфологической фильтрации и метода кластеризации k-means. Данные
методы использовались при сегментации кровеносных сосудов глазного дна.
Выявленные
преимущества
методов
показали
применения для решения поставленной задачи.
45
целесообразность
их
3. РАЗРАБОТКА АЛГОРИТМИЧЕСКОГО И ПРОГРАММНОГО
ОБЕСПЕЧЕНИЯ
В данной главе разработан метод и алгоритм сегментации кровеносных
сосудов на изображениях глазного дна с использованием контрастно
ограниченной адаптивной эквализации гистограммы, морфологической
фильтрации, метода кластеризации k-means отдельно толстых и тонких
сосудов. Для выявления тонких сосудов использована также согласованная
фильтрация, ядра которой позволяют выявить наличие отрезков заданной
длины и ориентации на плоскости. Также в данной главе описан интерфейс
программы,
и
приведены
фрагменты
листингов
кода
программной
реализации разработанного метода.
3.1 Разработка метода и алгоритма сегментации кровеносных
сосудов на изображении глазного дна
Анализ особенностей изображений глазного дна [12] показал, что
яркости пикселей, соответствующих толстым сосудам, имеют большее
отличие от яркости пикселей фона по сравнению с отличием яркости
пикселей фона от яркости пикселей, соответствующих тонким сосудам.
Учитывая данный факт, в работе предложено осуществлять выделение
на изображении отдельно тонких и толстых сосудов.
Проведенные
предварительные
вычислительные
эксперименты
показали, что при кластеризации пикселей изображения глазного дна,
используя значения их яркости, на k классов в старшем классе Sk содержится
большинство пикселей, соответствующих толстым кровеносным сосудам,
при этом значительное количество пикселей, соответствующих всем
46
кровеносным сосудам глазного дна, содержится в двух старших классах Sk 1
и Sk (что будет продемонстрировано в работе далее).
Одной из основных задач при разработке подхода к сегментации
сосудов является выделение тонких сосудов из класса Sk 1 , в котором также
присутствуют пиксели, не соответствующие сосудам. Для повышения
контрастности тонких сосудов в работе использована согласованная
фильтрация [13], ядра которой создаются с целью выявления присутствия
отрезков заданной длины и ориентации на плоскости.
Таким образом, основные шаги алгоритма сегментации сосудов
глазного дна можно сформулировать следующим образом:
1. Выполнить предварительную обработку изображения с целью
удаления шумов.
2. Осуществить кластеризацию множества пикселей обработанного
изображения I 0 на k классов.
3. Выполнить сегментацию толстых сосудов на изображении глазного
дна (сформировать бинарное изображение I B0 , используя пиксели старшего
класса Sk ).
4. Сформировать изображение I1 , содержащее пиксели класса Sk 1 .
5. Выполнить обработку, позволяющую выделить тонкие сосуды
(бинарное изображение I B1 ) на изображении I1 (выполнить повышение
контраста тонких сосудов на изображении на основе использования
согласованной фильтрации, а также выполнить повторную кластеризацию).
6. Добавить результат I B1 сегментации тонких сосудов к изображению
I B 0 , которое является результатом сегментации толстых сосудов.
Блок-схема приведенного алгоритма сегментации сосудов и блок-схема
сегментации отдельно толстых сосудов на изображении глазного дна
изображены на рис. 3.1.
47
а
б
Рис. 3.1. Блок-схемы алгоритмов:
а – блок-схема алгоритма сегментации сосудов глазного дна, б – блоксхема алгоритма сегментации отдельно толстых сосудов на изображении
глазного дна
48
При решении задач сегментации кровеносных сосудов глазного дна
зачастую возникает необходимость предобработки изображения с помощью
различных фильтров, поскольку сегментация исходного изображения обычно
не даёт ожидаемых результатов из-за неравномерной освещённости,
недостаточного контраста и др.
В работе использован алгоритм [43] предобработки изображения
глазного дна, суть которого состоит в следующем:
1. Преобразовать исходное цветное изображение, размерностью N1 N 2 ,
в изображение I c (i, j ) , i 1,2,..., N1 , j 1,2,..., N 2 , в оттенках серого (рис. 3.3а).
2. Выполнить повышение контраста полутонового изображения,
используя
метод
контрастно
ограниченной
адаптивной
эквализации
гистограммы (CLAHE). Пусть I – результат выполнения данного шага.
Гистограммы исходного и обработанного изображения приведены на
рис. 3.2.
а
б
Рис. 3.2. Гистограммы изображения:
а – гистограмма исходного изображения в оттенках серого,
б – гистограмма изображения после контрастно ограниченной адаптивной
эквализации
Как видно из рис. 3.2 значения яркости пикселей выровнялись,
распределившись по оси абсцисс в гораздо большем диапазоне. Данное
49
распределение
яркости
привело
к
повышению
контрастности,
что
продемонстрировано на рис. 3.3б.
3. Для дальнейшего выделения сосудов выполняется морфологическая
фильтрация,
которую
также
называют
модифицированным
top-hat
преобразованием (рис. 3.3в):
I 0 I ( I Se) Se ,
(3.1)
где I 0 – изображение, полученное в ходе морфологической фильтрации, I –
изображение после эквализации гистограммы, Se – структурный элемент (в
работе
применен
диск
радиусом
8
пикселей),
операции
,
–
морфологические операции размыкания и замыкания.
а
б
в
Рис. 3.3. Результаты обработки изображения:
а – исходное изображение в оттенках серого, б – результат применения
контрастно ограниченной адаптивной эквализации гистограммы,
в – результат применения морфологической фильтрации
Кластеризацию множества пикселей обработанного изображения I 0 на
k
классов S1 , S2 , …, Sk в работе предложено осуществлять на основе метода
k-means, являющимся одним из наиболее распространенных методов
кластеризации [44, 45].
50
При применении метода k-means разбиение на классы осуществляется
таким образом, чтобы минимизировать суммарное квадратичное отклонение
точек кластеров от центров этих кластеров:
k
( x j mi )2 ,
(3.2)
i 2 x j S i
где k – число кластеров, Si – полученные классы, i 1,2,..., k , mi – центр масс
точек x j Si .
Для формирования бинарного изображения I B0 , соответствующего
толстым сосудам на изображении глазного дна, используются пиксели
старшего класса Sk :
1, (i, j ) S k ,
I B 0 (i, j )
i 1,2,..., N1,
0, в противном случае,
Для
того,
соответствующие
чтобы
на
«шуму»,
изображении
из
изображения
I B0
j 1,2,..., N 2 .
удалить
(3.3)
компоненты,
I B 0 удаляются
элементы
связности, содержащие не более 2 пикселей, что способствует улучшению
выделения сосудов (рис. 3.4).
Рис. 3.4. Бинарное изображение I B0 толстых сосудов (класс S5 , k 5 )
51
Далее формируется изображение I1 (рис. 3.5), содержащее значения
пикселей, которые были включены в класс Sk 1 :
I (i, j ), (i, j ) S k 1,
I1 (i, j ) 0
i 1,2,..., N1,
0, в противном случае,
j 1,2,..., N 2 .
(3.4)
Рис. 3.5. Изображение I1 , соответствующее классу S4
Для повышения контраста тонких сосудов на изображении I1 (3.4) в
работе использована согласованная фильтрация F ( I1 , hL , ) , в основе которой
лежит свертка изображения I1 с ядром hL , имеющем размерность L L
пикселей и соответствующем отрезку, который проведен под углом к
горизонтальной оси. Параметр L принимает значения {7, 9, 11, 13, 15},
параметр принимает значения l 15(l 1) градусов, l 1,2,...,12 .
В этом случае согласованная фильтрация позволяет выделить линейно
протяженные области на изображении I1 , которые имеют различный
масштаб и различную ориентацию (на рис. 3.6 при L 15 элементы ядра,
равные 1, отображены в виде белых квадратов, равные 0 – в виде черных
квадратов).
52
а
б
в
г
д
Рис. 3.6. Примеры ядер, используемых при согласованной фильтрации,
L 15 :
а – при 0 , б – при 30 , в – при 90 , г – при 135 , д – при 150
В
работе
предложено
результат
согласованной
фильтрации
(изображение I S1 (рис. 3.7)) получать только для пикселей, входящих в класс
Sk 1 :
I (i, j ), (i, j ) S k 1 ,
I S 1 (i , j ) F 1
i 1,2,..., N1,
0
,
в
противном
случае
,
j 1,2,..., N 2 ,
(3.5)
где
7
I F1
12
F ( I1, h2 1 , l ) .
(3.6)
3 l 1
Алгоритм применения согласованной фильтрации заключается в
следующем:
1. Задать размерность ядра фильтра L=15;
2. Построить ядра фильтра размерностью L L , соответствующие
отрезкам (рис. 3.6) под углами от 0 до 165 градусов с шагом 15 градусов;
3. Применить данный фильтр к обрабатываемому изображению
(выполнить операцию свертки).
4. Изменить размерность ядра фильтра L=13.
5. Построить ядра фильтра размерностью L L , соответствующие
отрезкам под углами от 0 до 165 градусов с шагом 15 градусов;
53
6. Применить данный фильтр к обрабатываемому изображению.
Результат применения фильтра добавить к результату, полученному на шаге
3.
7. Повторить данные действия для размерностей ядра фильтра L=11, 9,
7. При этом результат фильтрации накапливается после каждого изменения
размерности ядра фильтра, то есть результат согласованной фильтрации
после применения ядра фильтра размерностью L=7 является суммой
результатов согласованной фильтрации при L=15, 13, 11, 9, 7.
Рис. 3.7. Результат I S1 выполнения согласованной фильтрации
Для выделения тонких сосудов в работе применен метод k-means для
разбиения множества пикселей изображения I S1 (3.5) на n классов R1 , …, Rn ,
при этом пиксели старшего класса Rn считаются пикселями тонких сосудов.
Тогда, формирование бинарного изображения I B1 , содержащего тонкие
сосуды (рис. 3.8), осуществляется на основе следующего выражения:
1, (i, j ) Rn ,
I B1 (i, j )
i 1,2,..., N1,
0, в противном случае,
54
j 1,2,..., N 2 .
(3.7)
Рис. 3.8. Бинарное изображение I B1 тонких сосудов (класс R3 , n 3 )
На заключительном этапе для формирования бинарного изображения
сосудов I vessel глазного дна выполняется сложение изображений I B 0 (3.3) и I B1
(3.7), соответствующих толстым и тонким сосудам (рис. 3.9):
I vessel I B 0 I B1 .
(3.8)
Рис. 3.9. Результат I vessel применения разработанного метода
В результате выполнения указанных выше действий получаем
бинарное изображение сосудов глазного дна, где черный цвет пикселей
соответствует фону, а белый – сосудам.
55
3.2 Интерфейс программной реализации
Программная
реализация
сегментации кровеносных сосудов на
изображении глазного дна реализована в виде оконного интерфейса (формы),
созданного в среде инженерных расчетов Matlab [46, 47].
Форма содержит 2 кнопки, 2 поля для вывода изображений, 2 поля для
ввода числа классов при кластеризации и 4 поля для вывода числовых
характеристик обработанного изображения (рис. 3.10).
Рис. 3.10. Первоначальный вид формы сегментации кровеносных сосудов на
изображении глазного дна
Кнопка «Загрузить изображение» служит для загрузки исходного
цветного изображения глазного дна. Данное изображение отображается слева
на форме (рис. 3.11).
Рис. 3.11. Загрузка исходного изображения глазного дна
56
Далее необходимо ввести количество классов для кластеризации
методов k-means. Верхнее поле предназначено для ввода количества классов
при кластеризации толстых сосудов. По умолчанию данное значение равно 5.
Пользователь может изменять данный параметр. Нижнее поле предназначено
для ввода количества классов при кластеризации тонких сосудов. По
умолчанию данное значение равно 3. В процессе вычислений пользователь
имеет возможность изменять и этот параметр.
После того, как исходное изображение загружено, и количество
классов для кластеризации установлено, необходимо нажать на кнопку
«Выполнить сегментацию». После чего, справа на форме отобразиться
бинарное обработанное изображение глазного дна, на котором белые пиксели
соответствуют сосудам, а черные – фону. Ниже данного изображения в
четырех полях «Точность»,
«Чувствительность»,
«Специфичность» и
«Ошибка первого рода» отобразятся соответствующие цифры, которые
характеризуют обработанное изображение (рис. 3.12) (соотношения для
определения значений данных параметров будут приведены в следующей
главе). В дальнейшем данные параметры будут использоваться для сравнения
качества сегментации сосудов на различных изображениях глазного дна.
Рис. 3.12. Отображение результата сегментации сосудов на изображении
глазного дна
Для закрытия формы следует нажать на крестик в правом верхнем углу
окна программы.
57
3.3 Программная реализация алгоритма сегментации кровеносных
сосудов на изображении глазного дна
В ходе данной работы была разработана программная реализация
приведенного ранее алгоритма сегментации кровеносных сосудов на
изображении глазного дна в среде Matlab. Исходный код программной
реализации состоит из 4 программных модулей. Схема взаимосвязи
программных модулей представлена на рис. 3.13.
Рис. 3.13. Модульная схема программной реализации
Файлом программы, содержащим точку входа, а также реализацию
интерфейса программы, является m-файл eye_fundus.m.
В данном m-файле реализованы основные операции разработанного
алгоритма, а также подключаются функция f_Param.m для вычисления
характеристик обработанного изображения (которые более подробно будут
рассмотрены в 4 главе) и функция f_getLineResponse1.m для выполнения
согласованной
фильтрации
изображения.
В
свою
очередь
функция
f_getLineResponse1.m подключает файл f_getLineMask.m для формирования
различных ядер фильтра.
В ходе программной реализации основных функций сегментации
кровеносных сосудов на изображении глазного дна разработаны следующие
фрагменты кода программы.
58
Фрагмент кода выполнения повышения контраста полутонового
изображения на основании применения метода контрастно ограниченной
адаптивной эквализации гистограммы (CLAHE) приведен в листинге 3.1.
Листинг 3.1 Применение метода CLAHE
figure('Name',['Исходная гистограмма'])
imhist(f1);
g=adapthisteq(f1,'clipLimit',0.01);
figure('Name',['g Эквализация гистограммы'])
imshow(g);
figure('Name',['g Эквализация гистограммы'])
imhist(g);
Конец листинга 3.1
здесь f1 – исходное изображение в оттенках серого. Предел урезания
(clipLimit) предлагается использовать равным 0.01, данное значение
параметра показало свою эффективность по сравнению с другими
значениями.
Выполнение морфологической фильтрации (модифицированное top-hat
преобразование) проиллюстрировано в листинге 3.2.
Листинг 3.2 Морфологическая фильтрация
se_size = 8;
se=strel('disk',se_size);
g0=imclose(g,se);
g1=imopen(g0,se);
g2=g-g1;
g3=max(g2(:))-g2;
figure('Name',['Top-Hat Transform’])
imshow(g3,[]);
Конец листинга 3.2
59
В данном фрагменте кода используется структурный элемент в виде
диска с радиусом 8 пикселей. Данное значение параметра было выбрано
путем проведения ряда вычислительных экспериментов.
Для кластеризации пикселей на классы по значениям их яркости
использовался метод кластеризации k-means. Применение данного метода
для кластеризации толстых сосудов приведено в листинге 3.3.
Листинг 3.3 Разбиение пикселей изображения на классы
X1 = g3(:);
cidx = kmeans(X1,NClass,'distance','sqeuclid');
len_cidx=length(cidx);
for k2=1:NClass
for i1=1:len_cidx
if cidx(i1)==k2
X1Class(k2)=X1(i1);
break;
end
end
end
[X1SortClass,IndX1Class]=sort(X1Class);
cidx1=cidx*0.;
for k2=1:NClass
cidx1(cidx==IndX1Class(k2))=k2;
end
Конец листинга 3.3
В данном случае разбиение изображения Х1 после морфологической
фильтрации происходило на NClass классов, значение данного параметра
указывает пользователь на форме (по умолчанию – 5). В качестве способа
вычислений
центров
каждого
класса
использовался квадрат Евклидова расстояния.
60
при
кластеризации
k-means
Формирование бинарного изображения, соответствующего толстым
сосудам на изображении глазного дна (используются пиксели старшего
класса с номером NClass), приведено в листинге 3.4.
Листинг 3.4 Формирование бинарного изображения, соответствующего
толстым сосудам
[N1,N2]=size(g3);
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
g3m(i1,i2)=cidx1(k1+i1);
if cidx1(k1+i1)== NumbClass1
g3mBI(i1,i2)=1;
end
end
end
g3mBId=g3mBI;
g3mBId=double(bwareaopen(g3mBI,2));
Конец листинга 3.4
Для формирования изображения, содержащего значения пикселей,
которые были включены в класс с номером (NClass-1), был разработан
листинг 3.5.
Листинг 3.5 Формирование изображения, содержащего значения пикселей,
которые были включены в класс с номером (NClass-1)
g3t=g3*0;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
if cidx1(k1+i1)== NumbClass2 && Mask1(i1,i2)>0
g3t(i1,i2)=g3(i1,i2);
end
end
end
Конец листинга 3.5
61
После указанных выше шагов необходимо выполнить согласованную
фильтрацию, чтобы сильнее выделить тонкие сосуды на изображении.
Фрагмент кода применения согласованной фильтрации приведен в листинге
3.6.
Листинг 3.6 Использование согласованной фильтрации
W1=15;
L1=7;
Resp1 = Resp1+f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
Resp2=Resp1;
Конец листинга 3.6
Функция
f_getLineResponse1,
выполняющая
непосредственно
согласованную фильтрацию при использовании различных ядер фильтра,
приведена в листинге 3.7.
Листинг 3.7 Согласованная фильтрация
function R = f_getLineResponse1(img,W,L)
[N1,N2]=size(img);
maxlinestrength = -Inf*ones(size(img));
for theta = 0:15:165
linemask = f_getLineMask(theta,L);
linemask = linemask/sum(linemask(:));
imglinestrength = imfilter(img,linemask);
maxlinestrength = max(maxlinestrength,imglinestrength);
end
R = maxlinestrength;
end
Конец листинга 3.7
62
Формирование бинарного изображения, соответствующего тонким
сосудам на изображении глазного дна, приведено в листинге 3.8.
Листинг 3.8 Формирования бинарного изображения, соответствующего
тонким сосудам
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
g3tm(i1,i2)=cidx1X2(k1+i1);
if cidx1X2(k1+i1)== NumbClass1X2
g3tmBI(i1,i2)=1;
end
end
end
Конец листинга 3.8
Расчет численных характеристик обработанного изображения приведен
в листинге 3.9.
Листинг
3.9
Вычисление
характерных
параметров
обработанного
изображения
[TP,TN,FP,FN]=f_Param(Etalon1,g3tmBId);
Acc_g3tm=(TP+TN)/(TP+TN+FP+FN);
Sn_g3tm=TP/(TP+FN);
Sp_g3tm=TN/(TN+FP);
AUC_g3tm=(Sn_g3tm+Sp_g3tm)/2;
FP_Error_g3tm=FP/(TP+FN);
Конец листинга 3.9
Таким образом, в данной главе был разработан метод и алгоритм
сегментации кровеносных сосудов на изображении глазного дна, а также
разработана его программная реализация в среде Matlab, включая разработку
оконного интерфейса. Полный листинг кода программы сегментации
кровеносных сосудов на изображении глазного дна приведен в Приложении.
63
4. ЭКСПЕРИМЕНТАЛЬНЫЕ ИССЛЕДОВАНИЯ
РАЗРАБОТАННОГО МЕТОДА СЕГМЕНТАЦИИ СОСУДОВ
ГЛАЗНОГО ДНА
В
данной
главе
проведены
сравнительные
вычислительные
эксперименты разработанного метода и известного метода сегментации
сосудов глазного дна. Результаты показали работоспособность и наличие
преимуществ разработанного метода по сравнению с известным методом.
4.1 Выбор критериев для экспериментальных исследований
Экспериментальные исследования проведены с целью проверки
работоспособности разработанного метода и сравнительного оценивания
результатов сегментации сосудов глазного дна, полученных с помощью
разработанного метода, а также с помощью известных методов [43].
Для оценивания работоспособности разработанного и известных
методов сегментации в работе предложено использовать критерии точность
Acc , чувствительность Sn и специфичность Sp [48]. Данные критерии
определяются следующим образом:
Acc ( Ntp Ntn ) /( Ntp N fp Ntn N fn ) ,
(4.1)
Sn Ntp /( N tp N fn ) ,
(4.2)
Sp N tn /( N tn N fp ) ,
(4.3)
где N tp – количество корректно определённых пикселей сосудов (true
positive), N tn – количество корректно определённых пикселей фона (true
64
negative), N fp – количество некорректно определённых пикселей сосудов
(false positive), N fn – количество некорректно определённых пикселей фона
(false negative).
Точность
Acc
(4.1)
показывает
общую
оценку
сегментации,
чувствительность Sn (4.2) показывает эффективность определения пикселей,
соответствующих сосудам (определяется как отношение N tp корректно
выделенных
пикселей
сосудов
к
общему
числу
N vessel
пикселей,
соответствующих сосудам),
Nvessel Ntp N fn ,
(4.4)
специфичность Sp (4.3) показывает эффективность определения пикселей,
соответствующих фону.
В качестве исходных данных для вычислительных экспериментах
использованы изображения из базы изображений глазного дна DRIVE [49].
4.2 Описание вычислительных экспериментов
Для проверки работоспособности разработанного метода сегментации
сосудов
глазного
дна
были
проведены
следующие
вычислительные
эксперименты.
Целью вычислительных экспериментов является определение степени
соответствия сосудам глазного дна пикселей, включенным в старшие классы,
которые получены при кластеризации изображения I 0 на основе метода kmeans.
65
В табл. 4.1 в качестве примера приведены для отдельных изображений
из базы DRIVE значения чувствительности Sn (4.2) и значение ошибки
первого рода EFP ,
EFP N fp /( Ntp N fn ) ,
(4.5)
EFP – отношение N fp некорректно выделенных пикселей сосудов к их
общему числу Nvessel , для двух старших классов Sk и S k 1 (3.2), полученных
при кластеризации методом k-means, для различных значений количества
классов k.
Таблица 4.1
Примеры соответствия сосудам глазного дна пикселей, включенным в
старшие классы при кластеризации изображения I 0
Изображение
04_test
08_test
15_test
16_test
18_test.
19_test
Sn
0,1897
0,3685
0,3045
0,2246
0,1501
0,3527
EFP
Sn
0,0002
0,6032
0,0228
0,8335
0,0051
0,7128
0,0017
0,6597
0,0111
0,6919
0,0042
0,8367
Sn
0,3447
0,4381
0,3666
0,4523
0,4132
0,5699
EFP
Sn
0,0031
0,8438
0,0415
0,8973
0,0096
0,7951
0,0160
0,8979
0,0396
0,9039
0,0186
0,9533
Sn
0,3856
0,0056
0,4497
0,0461
0,4906
0,0351
0,5155
0,0253
0,5120
0,0633
0,6116
0,0267
0,8886
0,9056
0,9391
0,9471
0,9516
0,9665
k
5
Sk
Sk Sk 1
4
Sk
Sk Sk 1
3
Sk
Sk Sk 1
EFP
Sn
Данные, приведенные в табл. 4.1, показывают, что при кластеризации
пикселей на k 5 классов в старшем классе Sk находится в большинстве
случаев значительное количество пикселей, соответствующих сосудам
глазного дна, при несущественной ошибке 1 рода. В двух старших классах Sk
и Sk 1 находится в большинстве случаев подавляющее количество корректно
выделенных пикселей сосудов.
66
Уменьшение
количества
классов
при
кластеризации
повышает
количество корректно выделенных пикселей сосудов, однако, при этом
существенно возрастает значение соответствующей ошибки 1 рода.
Для оценивания эффективности разработанного метода в работе
предложено использовать сравнение с одним из известных методов
сегментации сосудов, основные положения которого изложены в работе [43].
Сравнение разработанного метода с методом, предложенным в работе
[43], проводилось на основе критериев (4.1), (4.2), (4.3), используя
изображения из базы изображений DRIVE (табл. 4.2 и 4.3).
В табл. 4.2 приведены значения критериев (4.1), (4.2), (4.3)
разработанного метода.
Таблица 4.2
Примеры вычисления критериев Acc , Sn и Sp при использовании
разработанного метода
Изображение
04_test.tif
08_test.tif
15_test.tif
k
5
5
4
4
3
3
5
5
4
4
3
3
5
5
4
4
3
3
n
3
2
3
2
3
2
3
2
3
2
3
2
3
2
3
2
3
2
Acc
0,9482
0,9577
0,9031
0,9420
0,8801
0,7933
0,9066
0,9040
0,8253
0,8933
0,8083
0,7162
0,9623
0,9608
0,9588
0,9426
0,8616
0,7450
67
Sn
0,4468
0,5881
0,7717
0,7590
0,7856
0,8848
0,7181
0,7625
0,7614
0,7774
0,7648
0,8985
0,5484
0,7003
0,6674
0,7844
0,8193
0,9341
Sp
0,9990
0,9952
0,9165
0,9605
0,8896
0,7840
0,9243
0,9173
0,8313
0,9042
0,8124
0,6991
0,9942
0,9809
0,9812
0,9548
0,8648
0,7305
EFP
0,0101
0,0475
0,8247
0,3899
1,0895
2,1316
0,8039
0,8784
1,7922
1,0172
1,9926
3,1966
0,0751
0,2478
0,2438
0,5863
1,7537
3,4969
Продолжение таблицы 4.2
16_test.tif
18_test.tif
19_test.tif
5
5
4
4
3
3
5
5
4
4
3
3
5
5
4
4
3
3
3
2
3
2
3
2
3
2
3
2
3
2
3
2
3
2
3
2
0,9517
0,9594
0,9586
0,8561
0,8441
0,7345
0,9588
0,9616
0,9310
0,8562
0,8297
0,7188
0,9620
0,9690
0,9613
0,9280
0,8261
0,7165
0,4819
0,7290
0,6932
0,8924
0,8412
0,9414
0,5311
0,6410
0,8179
0,8986
0,8549
0,9470
0,5689
0,8053
0,7889
0,8919
0,8637
0,9604
0,9984
0,9822
0,9849
0,8525
0,8444
0,7139
0,9956
0,9892
0,9408
0,8526
0,8276
0,6992
0,9976
0,9839
0,9769
0,9312
0,8227
0,6944
0,0166
0,1788
0,1520
1,4857
1,5676
2,8822
0,0509
0,1255
0,6883
1,7133
2,0038
3,4959
0,0271
0,1785
0,2556
0,7601
1,9601
3,3786
Данные, приведенные в табл. 4.2, показывают, что предпочтительными
являются варианты кластеризации на 5 и 3, 5 и 2, а также 4 и 3 классов
соответственно при сегментации толстых и тонких сосудов.
При выборе кластеризации на 5 и 3 классов результат сегментации
имеет достаточно высокую точность Acc (4.1) и малую ошибку первого рода
EFP (4.5), однако в этом случае выделяются 50-60% (значение Sn (4.2))
пикселей, соответствующих сосудам.
При выборе кластеризации на 5 и 2, а также 4 и 3 классов результаты
сегментации имеют сравнимые значения критериев: при высокой точности
Acc ошибка первого рода EFP в большинстве случаев несколько меньше при
использовании варианта кластеризации на 5 и 2 классов, однако при варианте
кластеризации
на
4
и
3
классов
выделяется
больше
пикселей,
соответствующих сосудам [50].
Для сравнения эффективности применения разработанного метода в
табл. 4.3 приведены значения критериев (4.1), (4.2), (4.3) при использовании
метода, описанного в работе [43].
68
Таблица 4.3
Результаты вычисления критериев Acc , Sn и Sp при использовании
метода [43]
Acc
0,959
0,962
0,959
0,958
0,953
0,956
Изображение
04_test.tif
08_test.tif
15_test.tif
16_test.tif
18_test.tif
19_test.tif
Sp
Sn
0,704
0,782
0,821
0,726
0,808
0,760
0,984
0,975
0,970
0,977
0,964
0,977
Данные, приведенные в табл. 4.2 и 4.3, демонстрируют, что
разработанный метод в отдельных случаях имеет преимущество по
сравнению с методом [43] как в значениях точности Acc , так и в значениях
чувствительности Sn (определение пикселей сосудов) и специфичности Sp
(определение пикселей фона).
Для визуального анализа результатов применения разработанного
метода сегментации кровеносных сосудов глазного дна на тестовых
изображениях из базы DRIVE на рис. 4.1 – 4.6 приведены соответствующие
изображения.
a
б
в
Рис. 4.1. Сегментация изображения 04_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
69
a
б
в
Рис. 4.2. Сегментация изображения 08_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
a
б
в
Рис. 4.3. Сегментация изображения 15_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
a
б
в
Рис. 4.4. Сегментация изображения 16_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
70
a
б
в
Рис. 4.5. Сегментация изображения 18_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
a
б
в
Рис. 4.6. Сегментация изображения 19_test.tif:
а – исходное изображение, б – сегментация при кластеризации на 5 и 3
классов, в – сегментация при кластеризации на 4 и 3 классов
Таким
образом,
результаты
вычислительных
экспериментов,
приведенные в данной главе, продемонстрировали работоспособность
разработанного метода сегментации кровеносных сосудов глазного дна.
Точность вычислений сравнима с результатами применения известного
метода [43], и в отдельных случаях их превосходит (изображения 04_test.tif,
15_test.tif, 16_test.tif, 18_test.tif, 19_test.tif). Однако на изображении 08_test.tif
метод [43] работает лучше, что связано с наличием на данном изображении
разветвлённой сети тонких сосудов.
На полученных бинарных изображениях в результате сегментации с
помощью разработанного метода отчетливо видна сосудистая система
глазного дна.
71
ЗАКЛЮЧЕНИЕ
В процессе выполнения данной работы было выявлено, что цифровая
обработка изображений является актуальной темой в настоящее время и
находит широкое применение в различных областях медицины, включая
диагностику по глазному дну. Автоматическая сегментация кровеносных
сосудов глазного дна – важный шаг для выявления различных патологий в
глазе человека.
В ходе данной работы поставленная цель достигнута, все задачи
выполнены.
В данной работе были рассмотрены основные задачи и проблемы,
появляющиеся при исследовании снимков глазного дна, структуры глазного
дна, и диагностировании болезней, которые можно выявить по снимкам
данного
органа
зрения.
Выполнен
анализ
существующих
методов
сегментации объектов на изображении, а также анализ методов, которые
специализируются на сегментации непосредственно кровеносных сосудов на
изображении глазного дна.
Также были рассмотрены теоретические основы методов, применяемых
для
решения
ограниченная
фильтрация
задачи
сегментации
адаптивная
и
метод
сосудов.
эквализация
кластеризации
А
именно
гистограммы,
k-means,
для
контрастно
морфологическая
которых
были
сформулированы основные определения, вычислительные соотношения и
алгоритмы. Был выполнен анализ преимуществ и недостатков данных
методов.
В данной работе был разработан метод сегментации кровеносных
сосудов на изображениях глазного дна с использованием контрастно
ограниченной адаптивной эквализации гистограммы, морфологической
фильтрации, метода кластеризации k-средних отдельно толстых и тонких
сосудов. Предобработка изображения позволила более контрастно выделить
сосуды на глазном дне. Для выявления тонких сосудов также использована
72
согласованная фильтрация, ядра которой позволяют выявить наличие
отрезков заданной длины и ориентации на плоскости.
В среде инженерных расчетов Matlab разработана программная
реализация разработанного метода сегментации кровеносных сосудов,
которая продемонстрировала его работоспособность.
Метод был протестирован на изображениях из общедоступной базы
DRIVE. Вычислительные эксперименты показали наличие преимуществ
разработанного
метода
по
сравнению
с
известными
определении пикселей сосудов и фона изображения.
73
в
корректном
СПИСОК ИСПОЛЬЗУЕМОЙ ЛИТЕРАТУРЫ
1. Сетчатка / А.К. Хоу, Г.К. Браун, Д.А. МакНамара и др. – М.:
ГЭОТАР-Медиа, 2009.– 352 с.
2. Астахов, Ю. Офтальмоскопия / Ю. Астахов, Н.Даль. – СПб.:Н-Л,
2011. – 48 с.
3. Заболевания глазного дна / Дж. Дж. Кански, С. А. Милевски, Б. Э.
Дамато, В. Тэннер. – М.: МЕДпресс-информ, 2009.– 424 с.
4. Применение метода морфологических амёб для выделения сосудов
на изображениях глазного дна / А.В. Насонов, А.А. Черноморец, А.С.
Крылов, А.С. Родин. // Труды 13-й международной конференции "Цифровая
обработка сигналов и её применение" (DSPA'2011), т. 2, 2011. – с. 158–161.
5.
Батищев,
Д.С.
Инфраструктура
высокопроизводительной
компьютерной системы для реализации облачных сервисов хранения и
анализа данных персональной медицины / Д.С. Батищев, В.М. Михелев //
Научные ведомости БелГУ. Серия: Экономика. Информатика. Т. 37. № 2
(223), 2016, с. 88-92.
6. Sinthanayothin, C., Boyce, J.F., Williamson, T.H., Cook, H.L., Mensah,
E., Lal, S., Usher, D. Automated detection of diabetic retinopathy on digital fundus
images // Diabetic Medicine, Vol. 19(2), 2002, pp. 105–112.
7. Красильников, Н. Цифровая обработка 2D- и 3D- изображений / Н.
Красильников. – СПб.: БХВ-Петербург, 2011. – 608 с.
8. Кулешов, А. Методы и алгоритмы обработки изображений / А.
Кулешов, Л. Сушкова, Н. Шевченко. – Саарбрюккен: LAP Lambert Academic
Publishing, 2012. – 104 с.
9. Волосюк, В. Цифровая обработка сигналов и изображений / В.
Волосюк, О. Горячкин. – М.: ФИЗМАТЛИТ, 2007. – 552 с.
10. Otsu, N. A threshold selection method from gray-level histograms.
Automatica, 1975 Jun, 11(285–296): pp. 23–27.
74
11. Гонсалес, Р. Цифровая обработка изображений / Р. Гонсалес, Р.
Вудс. – М.: Техносфера, 2006. – 1072 с.
12. Ильясова, Н.Ю. Методы цифрового анализа сосудистой системы
человека. Обзор литературы / Н.Ю. Ильясова // Журнал «Компьютерная
оптика», т. 37 N 4, 2013, URL: http://www.computeroptics.smr.ru/.
13. Chaudhuri, S., Chatterjee, S., Katz, N., Nelson, M., Goldbaum, M.
Detection of Blood Vessels in Retinal Images Using Two-Dimensional Matched
Filters. IEEE Transactions of Medical Imaging, Vol. 8, No. 3, 1989, pp. 263–269.
14. Hoover, A., Kouznetsova, V., Goldbaum, M. Locating Blood Vessels in
Retinal Images by Piece-wise Threshold Probing of a Matched Filter Response.
IEEE Transactions on Medical Imaging, Vol. 19, No. 3, 2000, pp. 203–210.
15. Can, A., Shen, H., Turner, J.N., Tanenbaum, H.L., Roysam, B. Rapid
Automated Tracing and Feature Extraction from Retinal Fundus Images Using
Direct Exploratory Algorithms // IEEE Transactions on Information Technology in
Biomedicine, Vol. 3, No. 2, 1999, pp. 125-138.
16. Soares, J., Leandro, J., Cesar, Jr. R., Jelinek, H., Cree, M. Retinal Vessel
Segmentation Using the 2-D Gabor Wavelet and Supervised Classification. IEEE
Transactions on Medical Imaging, Vol. 25, No. 9, 2006, pp. 1214–1222.
17. Zana, F., Klein, J.C. Segmentation of vessel-like patterns using
mathematical morphology and curvature evaluation. IEEE Trans Image
Processing, Vol. 10(7), 2002. – p. 1010-1019.
18. Sun, K., Chen, Z., Jiang, S., Wang, Y. Morphological multiscale
enhancement, fuzzy filter and watershed for vascular tree extraction in angiogram.
Journal of Medical Systems, Vol. 35(5), 2011. – p. 811-824.
19. Li, Q., Jane, Y., Zhang, D. Vessel segmentation and width estimation in
retinal images using multiscale production of matched filter responses. Expert
Systems with Applications, Vol. 39(9), 2012. – p. 7600-7610.
20. Nguyen, U.T.V., Bhuiyan, A., Park, L.A.F., Ramamohanarao, K. An
effective retinal blood vessel segmentation method using multi -scale line
detection. Pattern Recognition, Vol. 46(3), 2013. – p. 703-715.
75
21. Biesdorf, A., Rohr, K., Feng, D., von Tengg-Kobligk, H., Rengier, F.,
Böckler, D., Kauczor, H.U., Wörz, S. Segmentation and quantification of the aortic
arch using joint 3D model-based segmentation and elastic image registration.
Medical Image Analysis, Vol. 16(6), 2012. – p. 1187-1201.
22. Fraz, M.M., Remagnino, P., Hoppe, A., Uyyanonvara, B., Rudnicka,
A.R., Owen, C.G., Barman, S.A. Blood vessel segmentation methodologies in
retinal images. Comput Methods Programs Biomed, Vol. 108(1), 2012. – p. 407433.
23. Osareh, A., Mirmehdi, М., Thomas, B., Markham, R. Color morphology
and snakes for optic disc localization. Pattern Recognition, Vol. 1, 2002. – p. 743746.
24. Jiang, X., Lambers, M., Bunke, H. Structural performance evaluation of
curvilinear structure detection algorithms with application to retinal vessel
segmentation. Pattern Recognition Letters, Vol. 33(15), 2012. – p. 2048-2056.
25. Patasius, M., Marozas, V., Jegelevieius, D., Lukosevieius, A. Recursive
Algorithm for Blood Vessel Detection in Eye Fundus Images: Preliminary Results
// IFMBE Proceedings, Vol. 25/11, 2009, pp. 212–215.
26. Ketcham, D. J., Lowe, R. W., Weber, J. W. Image enhancement
techniques for cockpit displays. Tech. rep., Hughes Aircraft, 1974.
27. Garg, R., Mittal, B., Garg, S. Histogram Equalization Techniques for
Image Enhancement, International Journal of Electronics and Communication
Technology, Vol.2, No, 2011. – p. 107-111.
28. Pizer, S. M., Amburn, E. P., Austin, J. D., et al. Adaptive Histogram
Equalization and Its Variations. Computer Vision, Graphics, and Image
Processing, 39, 1987. – p. 355-368.
29. Sund, T., Moystad, A. Sliding window adaptive histogram equalization
of intra-oral radiographs: effect on diagnostic quality. Dentomaxillofac
Radiol.;35(3), 2006. – p. 133-8.
76
30. Zuiderveld, K. Contrast Limited Adaptive Histogram Equalization. In: P.
Heckbert: Graphics Gems IV, Academic Press 1994, ISBN 0-12-336155-9. – p.
474-485.
31. Shome, S. K., Vadali, S. R. K. Enhancement of Diabetic Retinopathy
Imagery Using Contrast Limited Adaptive Histogram Equalization, International
Journal of Computer Science and Information Technologies, Vol. 2, No. 6, 2011. –
p.2694-2699.
32. Шапиро, Л. Компьютерное зрение / Л. Шапиро, Дж. Стокман. – М.:
БИНОМ. Лаборатория знаний, 2006. – 752 с.
33. Форсайт, Д. Компьютерное зрение. Современный подход / Д.
Форсайт, Ж. Понс. – М.: Вильямс, 2004. – 928 с.
34. Serra, J. Image Analysis and Mathematical Morphology, 1982, ISBN 012-637240-3.
35. Serra, J. Image Analysis and Mathematical Morphology, Volume 2:
Theoretical Advances, 1988, ISBN 0-12-637241-1.
36. Edward, R. D. An Introduction to Morphological Image Processing,
1992, ISBN 0-8194-0845-X.
37. Edward, R. D. Lotufo, R. Hands-on Morphological Image Processing,
2003, ISBN 0-8194-4720-X.
38. Steinhaus, H. Sur la division des corps materiels en parties. Bull. Acad.
Polon. Sci., C1. III vol IV, 1956. – p. 801—804.
39. MacQueen, J. Some methods for classification and analysis of
multivariate observations. In Proc. 5th Berkeley Symp. on Math. Statistics and
Probability, 1967. – p. 281—297.
40. Flury, B. Principal points. Biometrika, 77, 1990. – p. 33-41.
41. Gorban, A.N., Zinovyev, A.Y. Principal Graphs and Manifolds, Ch. 2 in:
Handbook of Research on Machine Learning Applications and Trends:
Algorithms, Methods, and Techniques, Emilio Soria Olivas et al. (eds), IGI Global,
Hershey, PA, USA, 2009. – p. 28-59.
77
42. Mirkes, E.M. K-means and K-medoids applet. University of Leicester,
2011.
43. BahadarKhan, K., Khaliq, A. A., Shahid, M. A Morphological Hessian
Based Approach for Retinal Blood Vessels Segmentation and Denoising Using
Region
Based
Otsu
Thresholding.
PLOS
ONE
11(7),
2016,
URL:
http://dx.doi.org/10.1371/journal.pone.0158996.
44. Steinhaus, H. Sur la division des corps materiels en parties. Bull. Acad.
Polon. Sci., C1. III vol. IV, 1956: 801—804.
45. Lloyd, S. Least square quantization in PCM’s. Bell Telephone
Laboratories Paper, 1957.
46. Гонсалес, Р. Цифровая обработка изображений в среде Matlab / Р.
Гонсалес, Р. Вудс, С. Эддинс. – М.: Техносфера, 2006. – 616 с.
47. Кривилев, А. Основы компьютерной математики с использованием
системы Matlab / А. Кривилев. – М.: Лекс-Книга, 2005. – 496 с.
48. Zhao, Y., Liu, Y., Wu, X., Harding, SP., Zheng, Y. Retinal vessel
segmentation: An efficient graph cut approach with retinex and local phase. PloS
one.
Apr
1,
2015,
10(4):e0122332.
doi:
10.1371/journal.pone.0122332.
pmid:25830353.
49. Research Section, Digital Retinal Image for Vessel Extraction (DRIVE)
Database, Utrecht, The Netherlands, University Medical Center Utrecht, Image
Sciences Institute: http://www.isi.uu.nl/Research/Databases/DRIVE/.
50. Черноморец, Д.А. О сегментации толстых и тонких сосудов / Д.А.
Черноморец, В.М. Михелев // Научные ведомости БелГУ. Сер. Экономика.
Информатика. – 2017. – № 16 (265). – Вып. 43. – С. 113-121. (из Перечня ВАК
РФ).
78
ПРИЛОЖЕНИЕ
Полный листинг кода программы сегментации кровеносных сосудов на
изображении глазного дна
global f;
global Etalon0;
global Mask1;
Etalon1=Etalon0/255;
NClass=str2double(get(handles.edit1, 'String'));
NClassX2=str2double(get(handles.edit2, 'String'));
f1 = im2double(rgb2gray(f));
[N,M]=size(f1);
figure('Name',['Исходная гистограмма'])
imhist(f1);
g=adapthisteq(f1,'clipLimit',0.01);
figure('Name',['g Эквализация гистограммы'])
imshow(g);
figure('Name',['g Эквализация гистограммы'])
imhist(g);
se_size = 8;
se=strel('disk',se_size);
g0=imclose(g,se);
figure('Name',['g0 Замыкание'])
imshow(g0);
g1=imopen(g0,se);
figure('Name',['g1 Размыкание'])
imshow(g1);
g2=g-g1;
g3=max(g2(:))-g2;
g41=g3;
figure('Name',['Top-Hat Transform без модуля g222=g3=g41'])
imshow(g3,[]);
clear g g0 g1 g2
[N1,N2]=size(g3);
X1 = g3(:);
NClass=5
NumbClass1=NClass;
79
NumbClass2=NClass-1;
cidx = kmeans(X1,NClass,'distance','sqeuclid');
len_cidx=length(cidx);
for k2=1:NClass
for i1=1:len_cidx
if cidx(i1)==k2
X1Class(k2)=X1(i1);
break;
end
end
end
[X1SortClass,IndX1Class]=sort(X1Class);
cidx1=cidx*0.;
for k2=1:NClass
cidx1(cidx==IndX1Class(k2))=k2;
end
g3m=g3*0;
g3mBI=g3*0;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
g3m(i1,i2)=cidx1(k1+i1);
if cidx1(k1+i1)== NumbClass1
g3mBI(i1,i2)=1;
end
end
end
figure('Name',['g3m Классы NClass=' int2str(NClass)])
imshow(g3m,[]);
figure('Name',['g3mBI Бинарное изображение NumbClass1=' int2str(NumbClass1)])
imshow(g3mBI,[]);
g3mBId=g3mBI;
g3mBId=double(bwareaopen(g3mBI,2));
figure('Name',['g3mBId Удаление 2 пикселей g3mBId'])
imshow(g3mBId);
[TP,TN,FP,FN]=f_Param(Etalon1,g3mBId);
Acc_g3m=(TP+TN)/(TP+TN+FP+FN);
Sn_g3m=TP/(TP+FN);
Sp_g3m=TN/(TN+FP);
AUC_g3m=(Sn_g3m+Sp_g3m)/2;
disp(['Acc_g3m=' num2str(Acc_g3m) ' ' num2str(Sn_g3m)...
' ' num2str(Sp_g3m) ' ' num2str(AUC_g3m)])
disp([' FP=' num2str(FP) ' FN=' num2str(FN)])
disp([' FP_Error=' num2str(FP/(TP+FN))...
80
' FN_Error=' num2str(FN/(TP+FN))])
FP_Error_g3m=FP/(TP+FN);
g3m2=g3m;
for i1=1:NClass-2
g3m2(g3m==i1)=0;
end
g3m2(g3m==NClass-1)=1;
g3m2(g3m==NClass)=1;
g3m2(Mask1==0)=0;
figure('Name',['g3m2 Бинарное изображение 2 старших класса'])
imshow(g3m2,[]);
[TP,TN,FP,FN]=f_Param(Etalon1,g3m2);
Acc_g3m2=(TP+TN)/(TP+TN+FP+FN);
Sn_g3m2=TP/(TP+FN);
Sp_g3m2=TN/(TN+FP);
AUC_g3m2=(Sn_g3m2+Sp_g3m2)/2;
disp(['Acc_g3m2=' num2str(Acc_g3m2) ' ' num2str(Sn_g3m2)...
' ' num2str(Sp_g3m2) ' ' num2str(AUC_g3m2)])
disp([' FP=' num2str(FP) ' FN=' num2str(FN)])
disp([' FP_Error=' num2str(FP/(TP+FN))...
' FN_Error=' num2str(FN/(TP+FN))])
Rez1=[Sn_g3m FP_Error_g3m Sn_g3m2]';
Min_g3t=Inf;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
if cidx1(k1+i1)== NumbClass1
if Min_g3t > g3(i1,i2)
Min_g3t=g3(i1,i2);
end
end
end
end
Min1_g3t=Min_g3t;
Min_g3t=Inf;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
if cidx1(k1+i1)== NumbClass2
if Min_g3t > g3(i1,i2)
Min_g3t=g3(i1,i2);
end
end
end
end
81
Min2_g3t=Min_g3t;
disp([' === Min1_g3t=' num2str(Min1_g3t) ' Min2_g3t=' num2str(Min2_g3t)])
g3t=g3*0;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
if cidx1(k1+i1)== NumbClass2 && Mask1(i1,i2)>0
g3t(i1,i2)=g3(i1,i2);
end
end
end
figure('Name',['g3t Пиксели класса NumbClass2=' int2str(NumbClass2)])
imshow(g3t,[]);
g3tt=g3t;
g3tt(Etalon1==1)=Min_g3t;
figure('Name',['g3tt Пиксели класса NumbClass2=' int2str(NumbClass2)...
' без Etalon1'])
imshow(g3tt,[]);
clear g3tt
W1=15;
L1=15;
Resp1 = f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
figure('Name',['Resp1 = f_getLineResponse L1=' int2str(L1)])
imshow(Resp1);
L1=13;
Resp1 = Resp1+f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
figure('Name',['Resp1 = f_getLineResponse L1=' int2str(L1)])
imshow(Resp1);
L1=11;
Resp1 = Resp1+f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
figure('Name',['Resp1 = f_getLineResponse L1=' int2str(L1)])
imshow(Resp1);
L1=9;
Resp1 = Resp1+f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
figure('Name',['Resp1 = f_getLineResponse L1=' int2str(L1)])
imshow(Resp1);
L1=7;
Resp1 = Resp1+f_getLineResponse1(g3t,W1,L1);
Resp1(g3t==0)=0;
82
figure('Name',['Resp1 = f_getLineResponse L1=' int2str(L1)])
imshow(Resp1);
Resp2=Resp1;
figure('Name',['Resp2+g3t '])
imshow(Resp2+g3t);
levelg3t=0.65
g3tBI = double(im2bw(Resp2,levelg3t));
figure('Name',['Пороговая обработка изображения g3tBI'])
imshow(g3tBI,[]);
g3tBId=g3tBI;
figure('Name',['g3tBId+g3mBId'])
imshow(g3tBId+g3mBId,[]);
[TP,TN,FP,FN]=f_Param(Etalon1,g3tBId+g3mBId);
Acc_g3t=(TP+TN)/(TP+TN+FP+FN);
Sn_g3t=TP/(TP+FN);
Sp_g3t=TN/(TN+FP);
AUC_g3t=(Sn_g3t+Sp_g3t)/2;
disp(['Acc_g3t=' num2str(Acc_g3t) ' ' num2str(Sn_g3t)...
' ' num2str(Sp_g3t) ' ' num2str(AUC_g3t)])
disp([' FP=' num2str(FP) ' FN=' num2str(FN)])
disp([' FP_Error=' num2str(FP/(TP+FN))...
' FN_Error=' num2str(FN/(TP+FN))])
X2 = Resp2(:);
NumbClass1X2=NClassX2;
NumbClass2X2=NClassX2-1;
cidxX2 = kmeans(X2,NClassX2,'distance','sqeuclid',...
'emptyaction','singleton');
len_cidxX2=length(cidxX2);
for k2=1:NClassX2
for i1=1:len_cidxX2
if cidxX2(i1)==k2
X2Class(k2)=X2(i1);
break;
end
end
end
[X2SortClass,IndX2Class]=sort(X2Class);
cidx1X2=cidxX2*0.;
for k2=1:NClassX2
cidx1X2(cidxX2==IndX2Class(k2))=k2;
83
end
g3tm=g3t*0+Min_g3t;
g3tmBI=g3t*0;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
g3tm(i1,i2)=cidx1X2(k1+i1);
if cidx1X2(k1+i1)== NumbClass1X2
g3tmBI(i1,i2)=1;
end
end
end
Min1_Resp2=Inf;
for i2=1:N2
k1=N1*(i2-1);
for i1=1:N1
if cidx1X2(k1+i1)== NumbClass1X2
if Min1_Resp2 > Resp2(i1,i2)
Min1_Resp2=Resp2(i1,i2);
end
end
end
end
disp([' === Min1_Resp2=' num2str(Min1_Resp2) ])
figure('Name',['g3tm Классы NClassX2=' int2str(NClassX2)])
imshow(g3tm,[]);
figure('Name',['g3tmBI Бинарное изображение NumbClass1X2=' int2str(NumbClass1X2)])
imshow(g3tmBI,[]);
g3tmBId=g3tmBI+g3mBId;
axes(handles.axes2);
set(handles.axes2);
imshow(g3tmBId,[]);title('Результат сегментации');
figure('Name',['g3tmBId'])
imshow(g3tmBId,[]);
[TP,TN,FP,FN]=f_Param(Etalon1,g3tmBId);
Acc_g3tm=(TP+TN)/(TP+TN+FP+FN);
Sn_g3tm=TP/(TP+FN);
Sp_g3tm=TN/(TN+FP);
AUC_g3tm=(Sn_g3tm+Sp_g3tm)/2;
FP_Error_g3tm=FP/(TP+FN);
set(handles.edit3,'String',num2str(Acc_g3tm));
set(handles.edit4,'String',num2str(Sn_g3tm));
set(handles.edit5,'String',num2str(Sp_g3tm));
set(handles.edit6,'String',num2str(FP_Error_g3tm));
84
Выпускная квалификационная работа выполнена мной совершенно самостоятельно.
Все использованные в работе материалы и концепции из опубликованной научной
литературы и других источников имеют ссылки на них.
«___» ________________ _____ г.
__________________________
_____________________
(подпись)
(Ф.И.О.)
85
Отзывы:
Авторизуйтесь, чтобы оставить отзыв