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

Во введении обоснована актуальность выбранной темы диссертационной работы, охарактеризована степень ее разработанности, определены цели и задачи, осуществлен выбор предмета и объекта исследования.
В первой главе проведено исследование предметной области и дана постановка задач диссертационной работы. Выполнен аналитический обзор существующих методов выявления и анализа областей интереса на медицинских снимках, проведена систематизация и классификация существующих программных средств выделения зон интереса врача-исследователя. В частности, рассмотрены программные комплексы визуализации и обработки томографических данных, методы сегментации МРТ-изображений для выделения зон интереса, методы морфинга и сопоставления множеств точек для решения задачи трекинга скоплений стволовых клеток. Сформулированы требования по созданию программного комплекса врача-исследователя для изучения областей интереса на снимках МРТ, в которые входят зоны ишемического поражения мозга и скопления трансплантированных МСК.
Вторая глава посвящена исследованию и разработке методов обнаружения зон ишемического поражения и поиска стволовых клеток на МРТ-снимках. Для анализа влияния МСК на очаги ишемического инсульта мозга в качестве экспериментальных данных использовались серии снимков головного мозга крыс, пораженных ишемией, полученные сразу после введения стволовых клеток и через 13 дней после введения. Каждая серия содержала 30 срезов головного мозга, удобных для анализа ишемического поражения (в режиме T2), и 40 срезов – для трекинга МСК (в режиме SWI (Susceptibillity Weighted Imaging)).
В разделе 2.1 приведены разработанные методы предварительной обработки входных данных, которые включают: 1) выравнивание освещенности МРТ-снимков; 2)выделение области мозга крысы, которое необходимо для последующего вычисления информативных параметров и построения 3D- моделей; 3) нормализацию области мозга.
Выравнивание освещенности необходимо из-за неравномерной чувствительности томографа. Оно выполняется с помощью попиксельного наложения серо-белого изображения-градиента (см. Рис. 1,а) на исходное изображение (см. Рис. 1,б) по следующей формуле:
если a  127
,
иначе
где a – яркость пикселя у исходного изображения, b – яркость пикселя у изображения-градиента с теми же координатами. Значения яркостей пикселей у
 ab , f a,b 127
2a  2b  ab  255,  127
обоих изображений изменяются от 0 (черный цвет) до 255 (белый цвет). Результат выравнивания яркости приведен на рисунке 1,в. Метод был применен для последующего исследования Т2-снимков МРТ и позволил добиться повышения качества сегментации.
Алгоритм выделения области мозга основан на использовании результатов работы утилиты BET (Brain Extraction Tool). Из-за того, что BET изначально предназначена для обработки снимков человеческого мозга, потребовалась ее адаптация применительно к снимкам мозга крысы. Для решения этой задачи был разработан алгоритм улучшенного выделения области мозга.
а) изображение- б) исходное в) обработанное
градиент
изображение изображение Рис. 1. Пример выравнивания яркости
Общая схема алгоритма включает несколько шагов:
1.по результатам работы «BET» выделяется предварительный контур
области мозга, после чего вычисляется центр тяжести полученной зоны (см. Рис.2,а, контур выделен красным цветом; точка центра тяжести показана пурпурным цветом);
2. с помощью специального алгоритма строится новый уточненный контур мозга на основе текущего контура (см. Рис. 2,а, точки нового контура выделены бирюзовым цветом);
3. после вычисления координат нового контура строится предварительная невыпуклая оболочка головного мозга (см. Рис. 2,б);
4.выполняется уточнение полученной оболочки за счет применения алгоритмов эрозии и дилатации (см. Рис. 2,в).
Поскольку каждое лабораторное животное обладает различными индивидуальными особенностями строения, то при обработке снимков МРТ- исследования головного мозга возникает задача сопоставления, или регистрации, целевых изображений с некоторым эталоном.
а) б) в) Рис. 2. Выделение области мозга (режим Т2)
6

Рассмотрим алгоритм корректного сопоставления двух изображений. Пусть C  c ,c ,…,c ,…,c и Q  q ,q ,…,q ,…,q – контуры эталонного и обрабатываемого
12in 12in
изображений соответственно, где ci , qi – точки контуров, а n – количество точек.
Первый этап включает ригидную аффинную трансформацию входного изображения. Основные шаги этой трансформации отражены на Рис. 3. Серым цветом показан контур обрабатываемого изображения, черным – контур эталонного объекта.
а) начальное б) выравнивание в) поворот контура положение центров
Рис. 3. Ригидная аффинная трансформация входного изображения Алгоритм ригидной аффинной трансформации состоит из следующих
шагов:
 производится переиндексация точек контура приводимого изображения N
так, чтобы d(c ,q )  min , где d(c ,q ) – Евклидово расстояние между c и q ;
iiiiii i1
1N 1N
 пусть c  N ci , q  N qi – центры контуров; сдвиг приводимого
i1 i1
изображения выполняется таким образом, чтобы центры контуров совпадали (см.
Рис. 3,б);
 вычисляется угол поворота   N i , где i – угол между
соответствующими точками ci и qi относительно центра; для ускорения операции
выполняется целочисленный поворот эталона на угол  (см. Рис. 3,в). Следующий этап состоит из деформации входного изображения, выполняемой таким образом, чтобы контуры головного мозга совпадали.
Алгоритм деформации включает следующие шаги:
 из общего центра проводятся лучи, которые пересекают оба контура;
 у каждого луча вычисляется коэффициент растяжения (сжатия) ki , где i –
номер луча; пусть ci и qi – точки контуров, принадлежащие одному лучу, тогда ki d(qi,q)/d(ci,c),причемqc;
для каждого луча выполняется растяжение (сжатие) точек эталона, которые принадлежат данному лучу, при помощи коэффициента ki .
Раздел 2.2 посвящен описанию методов выделения зон ишемического поражения на снимках МРТ, полученных в разных режимах работы томографа (Т2 и ИКД-карты). Особенность разработанных методов – работа с входными данными, имеющими низкое разрешение, что сильно затрудняет дальнейший анализ.
1N i1
Центральную часть метода выделения области ишемии на снимках, полученных в режиме Т2, составляют алгоритмы автоматического извлечения характерных признаков Харалика и классификатор на основе обобщенного расстояния Евклида-Махаланобиса. Расстояние между классом Y (например, классом «ишемическое поражение») и точкой p (под точкой понимается вектор признаков, которые извлечены из некоторой позиции сканирующего окна)
T 1
вычисляется по формуле R p,Y p y A p y, где матрица AC E.
GY Ковариационная матрица CY для класса Y вычисляется по формуле
Y skysky iijj
Cy i, j k1 , Y 1
где sk – i -ый признак k-ой точки из класса Y , а y – мат. ожидание значения этого ii
признака в классе Y . Перед обработкой изображений необходимо задать два класса – определим область здорового мозга как «норма» (класс Y1 ), а область ишемического поражения определим как «поражение» (класс Y2 ). Каждый класс
представляется набором небольших эталонных изображений одного из срезов мозга. Из этих эталонов извлекаются признаки, которые необходимы для настройки классификатора. Алгоритм извлечения характерных признаков построен на обработке изображения размером M×N с помощью сканирующего окна и имеет сложность O(MN). Алгоритм анализа извлеченных признаков с помощью классификатора выполняется в том же цикле, что и предыдущий этап, следовательно, он также имеет сложность O(MN). Длительность однократного выполнения разработанного алгоритма оценивается следующим образом: 0,015MN. Пример работы алгоритма выделения ишемического поражения приведен на Рис. 4.
а) разметка эксперта б) результат работы в) результат фильтрации алгоритма
Рис. 4. Фильтрация артефактов
Не исключены ложные срабатывания классификатора, как отображено на
рисунке 4,б. За счет дополнительной постобработки путем фильтрации мелких областей возможно добиться улучшения точности. Суть фильтрации состоит в сравнении по яркости и размерам выделенных объектов на изображении между собой – на итоговом изображении остается самый яркий и большой объект.
Для обработки ИКД-карт используется аппарат искусственных нейронных сетей. В качестве классификатора применяется сверточная нейронная сеть (СНС), состоящая из слоев свертки и субдискретизации, которые чередуются между собой и дополнены специализированными слоями. На выходе СНС может иметь один или несколько полносвязных слоев. Последний слой нейронной сети
содержит количество нейронов, которое соответствует числу классов в решаемой задаче. В проведенных экспериментах применялась нейронная сеть на базе библиотеки nnForge. На Рис. 5 приведены результаты обработки ИКД-снимков с применением СНС.
а) исходное изображение с разметкой эксперта б) результат обработки Рис. 5. Результат обработки МРТ-данных с применением СНС
В разделе 2.3 дается описание разработанного метода выделения скоплений МСК на снимках МРТ, полученных в режиме SWI. Метод основан на многоэтапной обработке пар снимков – до введения и после введения МСК. Схема разработанного алгоритма включает следующие этапы: предварительное выделение области мозга; выделение объектов из фона за счет применения волнового алгоритма; фильтрация по яркости и удаление объектов, пересекающихся со схожими объектами на снимке до введения клеток; фильтрация обнаруженных объектов методом наложения окна специальной формы; обработка с последующим совмещением результатов; когнитивная визуализация выделенных стволовых клеток на исходном изображении.
Выделение области мозга выполняется с помощью метода, описанного в разделе 2.1. Полученное изображение головного мозга, отделенное от фона, обрабатывается с помощью волнового алгоритма. Фильтрация изображения по яркости размером M×N происходит по следующему правилу: рассматривается каждый пиксель изображения, и если яркость пикселя превышает заданный порог, то он удаляется. Удаление неинформативных объектов состоит из последовательных операций над двумя изображениями до и после введения стволовых клеток: маркировка оставшихся после предыдущих шагов объектов; подсчет площади (количества пикселей) у каждого объекта;фильтрация. На Рис. 6 представлен результат работы алгоритма выделения скоплений МСК.
Рис. 6. Результат работы алгоритма поиска МСК
Длительность однократного выполнения разработанного алгоритма оценивается следующим образом: 0,03MN + 90,78.
В третьей главе рассмотрены методы анализа и визуализации зон интереса врача-исследователя по МРТ-данным. Проведено исследование динамики изменения областей интереса с помощью разработанных методов.
В разделе 3.1 описаны разработанные методы совмещенной 2D- и 3D- визуализации головного мозга и зон ишемического поражения. Дано описание алгоритма для расчета размеров областей интереса на МРТ-снимках.
Совмещенная визуализация комбинирует результаты обработки снимков, полученных в T2- и SWI-режимах томографа (Рис. 7). Однако 2D-визуализация, несмотря на определенную информативность и выразительность, не дает полной картины о текущем состоянии головного мозга. По этой причине строятся объемные модели зон интереса, где в качестве входных данных используются результаты независимой обработки отдельных МРТ-срезов.
Рис. 7. Совмещенная 2D-визуализация
Для получения точек границы областей использовалась реализация
детектора границ Canny библиотеки OpenCV, после чего выполнялось построение на их основе невыпуклой оболочки за счет применения метода «ConcaveHull». В этом методе расчет вогнутой оболочки осуществляется на основе метода Джарвиса с использованием предложенной в работе (Moreira A., Yasmina M.S., 2007) методики выбора очередной точки вогнутого контура. В итоге точки всех полученных вогнутых оболочек сортируются таким образом, чтобы далее объединить все срезы в единую 3D-оболочку. В итоге мозг крысы визуализируется в виде сеточной модели.
Задача построения 3D-моделей областей интереса обладает следующими особенностями: на одном и том же срезе МРТ-снимка может находиться несколько областей ишемического поражения; из-за возможности наличия отверстий в областях поражения обычная триангуляция не подходит для формализованного описания их поверхностей. По этим причинам для решения поставленной задачи был разработан специальный подход, состоящий из двух этапов. На первом этапе определяются и маркируются все области поражения и отверстия в них. Второй этап – определение контуров и выполнение триангуляции Зейделя, состоящей из трех шагов: 1) разбиение полигона области поражения и полигонов-отверстий в них на трапеции; 2) разделение трапеций на монотонные полигоны с вспомогательными диагоналями; 3) триангуляция
полигонов с помощью «жадного» алгоритма, многократно отсекающего выпуклые углы. В случае построения объемной модели скоплений стволовых клеток применяются сферы, размер которых зависит от количества сгруппированных МСК в скоплениях. Совмещенная 3D-визуализация зон интереса также позволяет изучить динамику движения МСК в течение времени (Рис. 8).
В разделе 3.2 решается задача отслеживания предполагаемого движения скоплений стволовых клеток за счет применения следующих методов: 1) анализа срезов головного мозга в хронологической последовательности; 2)фиксации траекторий движения с помощью кривой Безье; 3) решение транспортной задачи.
а) сразу после введения б) через три дня в) через 13 дней МСК
Рис. 8. Отслеживание тренда скоплений МСК и динамики изменения областей ишемического поражения
Построение карт миграции МСК выполнялось за счет сочетания автоматизированных алгоритмов выделения и визуализации скоплений МСК и решения транспортной задачи. На Рис. 9 представлена когнитивная визуализация полученного решения транспортной задачи. Желтым цветом выделено скопление МСК через 7 дней после введения, которое оказалось не соединенным ни с одним другим скоплением МСК сразу после введения, поскольку данный случай является транспортной задачей открытого типа.
В разделе 3.3 решены задачи расчета информативных параметров областей интереса и выполнено 3D-моделирование трекинга скоплений МСК с применением транспортного подхода.
Рис. 9. Когнитивная визуализация решения транспортной задачи
Результаты моделирования и измерений проанализированы на предмет влияния трансплантированных клеток на очаги ишемического инсульта головного мозга. Одновременно с построением траекторий движения скоплений МСК проводится отслеживание пересечения этими траекториями зоны ишемического поражения (см. Рис. 10). Совмещенное представление ишемического поражения в разные моменты времени и визуализация решения транспортной задачи в 3D- формате приведены на Рис. 11. Зелеными сферами обозначены начальные позиции скоплений МСК, красными – конечные позиции.
Рис. 10. Траектории движения МСК
а) сразу после введения МСК б) через 13 дней после введения МСК
Рис. 11. Совмещение зон ишемического поражения и движения скоплений МСК
На Рис. 12 демонстрируется динамика изменения площадей ишемического поражения по срезам (сразу после введения и через 13 дней после введения).
3500
3000
2500
2000
1500
1000
500 0
24 25 26 27 28 29 30 31 Срез Сразу п.в. Через 13 дней п.в.
Рис. 12. Динамика изменения области ишемического поражения
Площадь очага ишемии

В таблице 1 отражена общая информация о результатах эксперимента. Таблица 1 – Общая информация об изменении состояния областей мозга и ишемического поражения
1810740 141314 7,80 % 1810740 60569 3,34 %
Из полученных данных, приведенных в таблице выше, следует, что очаг ишемии существенно уменьшился за 13 дней. Полученные результаты эксперимента согласуются с врачебными исследованиями по изучению движения и влияния стволовых клеток на очаги ишемического инсульта (Кониева А.А. Влияние экзогенных мезенхимальных стволовых клеток плаценты человека на динамику некоторых патологических процессов ЦНС в эксперименте // Дисс. … канд. медиц. наук. – Москва, 2010. – 117 с.).
В четвертой главе рассмотрена структура разработанного программного комплекса врача-исследователя, в основу которого взяты методы, описанные ранее в главах 2 и 3. Описаны формы графического интерфейса программного комплекса, функциональные возможности, особенности реализации. Проведено сравнение разработанного комплекса с аналогами. На Рис. 13 представлена одна из основных форм интерфейса.
Рис. 13. Основное окно интерфейса
При обработке серии снимков можно задать диапазон обрабатываемых
срезов мозга, после чего его можно уточнить, указав первый и последний кадры. После того, как все задачи выделения областей интереса выполнены, появляется возможность просмотра статистических данных. Эта форма содержит результаты автоматического расчета важных характеристик зон интереса.
Состояние сразу после введения МСК
Состояние через 13 дней после введения МСК
Объем мозга, воксели
Объем ишемии, воксели
Доля повреждения, %
Объем мозга, воксели
Объем ишемии, воксели
Доля повреждения, %
13

Сравнение c аналогами показывает, что разработанные средства автоматического выделения зон интереса по функциональным возможностям не уступают имеющимся за рубежом средствам. Наиболее близкие по возможностям программные решения содержатся в пакетах FMRIB Software Library («FSL»), 3D Slicer и NITRC. Сравнение с ними показало преимущество разработанного комплекса в части возможностей моделирования, анализа и визуализации трекинга стволовых клеток.
В заключении приведены основные результаты, полученные в диссертационной работе.
Основные результаты работы
1. Разработаны методы предварительной обработки МРТ-снимков, позволяющие решить задачи выравнивания освещенности изображения, выделения области мозга из фона и приведения выделенной области мозга к эталону.
2. Разработаны методы выделения зоны ишемического поражения головного мозга по данным МРТ, полученных в разных режимах работы томографа (Т2 и ИКД-карты), за счет применения спектрографического анализа и аппарата искусственных нейронных сетей.
3. Приближенно решена задача трекинга скоплений стволовых клеток с помощью разработанного метода выделения скоплений МСК, основанного на эвристических критериях выделения стволовых клеток на сложном фоне, воссоздания траектории движения МСК с помощью кривой Безье и решения транспортной задачи.
4. На основе разработанных методов создан графический интерфейс врача- исследователя, функциональные возможности которого позволяют отслеживать динамику движения стволовых клеток на фоне ишемического инсульта за счет комбинированной визуализации и построения 2D- и 3D-моделей зон интереса. Одновременно с построением моделей вычисляются их информативные параметры.
5.Проведенный анализ динамики изменения областей интереса с применением разработанного инструментария показал, что стволовые клетки, попадая в очаги ишемического инсульта, уменьшают их размеры.