Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Оглавление

• Основные определения, термины и сокращения
• Введение
• Глава 1. Обзор предметной области и существующих решений
• 1.1 Обзор предметной области
• 1.1.1 Онкология
• 1.1.2 Проблема сегментации опухоли
• 1.2 Анализ существующих решения
• 1.2.1 OsiriX
• 1.2.2 Orthanc
• 1.3 Постановка задачи на разработку
• Выводы по главе
• Глава 2. Методы и алгоритмы сегментации раковых опухолей на МРТ снимках
• 2.1 Предобработка и постобработка снимков
• 2.1.1 Удаление шума
• 2.1.2 Удаление ярких краев
• 2.1.3 Морфологические операции
• 2.1.5 Метод Отцу
• 2.1.4 Удаление черепа
• 2.2 Сегментация с помощью пороговой фильтрации
• 2.3 Сегментация с помощью кластеризации
• 2.3.1 Метод k-средних
• 2.3.2 Метод нечеткой кластеризации c-средних
• 2.3.3 Метод сдвига среднего
• 2.4 Сегментация с помощью Гауссовых смесей
• 2.5 Сегментация с помощью расширяющихся областей
• 2.6 Общие вычисления на GPU
• 2.7 Оценка качества сегментации
• Выводы по главе
• Глава 3. Особенности программной реализации
• 3.1. Требования к программе
• 3.2 Инструменты разработки (сборки) и среда исполнения программы
• 3.2.1 Среда разработки
• 3.2.2 Веб-сервер
• 3.2.3 База данных
• 3.2.4 Работа с DICOM
• 3.2.5 WebGL и ThreeJS
• 3.4 Функциональные возможности программы
• 3.4.1 Загрузка МРТ снимков
• 3.4.2 Поиск МРТ снимков
• 3.4.3 Визуализация МРТ снимков
• 3.4.4 Сегментация раковых опухолей
• 3.4.5 Установка и удаление плагинов
• Выводы по главе
• Заключение
• Список использованной литературы
• Приложение а. ТЕХНИЧЕСКОЕ ЗАДАНИЕ
• Приложение Б. РУКОВОДСТВО ОПЕРАТОРА
• Приложение В. РУКОВОДСТВО ПРОГРАММИСТА
• Приложение Г. ПРОГРАММА И МЕТОДИКА ИСПЫТАНИЙ
• Приложение д. ТЕКСТ ПРОГРАММЫ


Основные определения, термины и сокращения

CUDA (Compute Unified Device Architecture) - программно-аппаратная архитектура и спецификация параллельных вычислений на графических процессорах компании Nvidia.

DICOM (Digital Imaging and Communications in Medicine) -- стандарт хранения, обработки и передачи медицинской информации и цифровых медицинских изображений.

GPGPU (General-purpose Computing for Graphics Processings Units) - техника выполнения алгоритмов общего назначения на графических процессорах.

OpenCL (Open Computing Language) - спецификация общих вычислений на графических процессорах разных производителей, разработанная консорциумом Khronos Group.

PACS (Picture Archiving and Communication System) -- сервер хранения и передачи медицинской информации и цифровых медицинских изображений.

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

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

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

Сегментация изображения -- процесс разбиения изображения на некоторые непересекающиеся области по некоторому признаку (критерию).



Введение

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

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

1. Высокая стоимость проведения процедуры;

2. Чувствительность к настройкам и особенностям аппарата;

3. Чувствительность к случайным шумам;

4. Ручной анализ занимает много времени.

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

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

- Проанализировать предметную область;

- Выявить проблематику;

- Изучить особенности существующих решений и выявить их недостатки;

- Провести анализ методов сегментации с помощью обработки изображений;

- Разработать архитектуру приложения;

- Разработать серверное приложение на Tornado;

- Разработать клиентскую часть с использованием WebGL;

- Разработать плагины для решения задачи сегментации.

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

5. 

Глава 1. Обзор предметной области и существующих решений

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

1.1 Обзор предметной области

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

1.1.1 Онкология

Рак занимает 7 место из 10 самых смертельных заболеваний, уступая заболеваниям сердечнососудистой системы, а также туберкулезу и другим респираторным заболеваниям [34]. Согласно исследованию Cancer Research UK, только в 2012 году было зафиксировано 14.1 миллионов случаев рака, а также около 8.2 миллиона случаем случаев со смертельным исходом [1]. Более того, ожидается рост числа раковых заболеваний в течение двух последующих десятилетий. Несмотря на то, что рак головного мозга не является самым распространённым видом рака (всего лишь около 6.4 случаев на 100000 каждый год), однако идентификация опухоли головного мозга одна из самых сложных задача, обусловленная тем, что опухоль плохо просматривается на МРТ снимках и не имеет четких краев. Более того, глиобластома и глиома, как частные виды рака мозга, одни из самых агрессивных опухолей с быстрым ростом, при котором пациент может прожить в среднем от 3 месяцев без лечения и до 2 лет с лечением [26]. Лишь меньше 5% всех пациентов могут прожить до 5 лет [26]. К тому, же лечение глиобластомы и глиомы осложняется тем, что данные типы опухоли могут образовывать вокруг себя большой отек и некроз, что делает опухоли не полностью операбельными [26]. В крайних случаях глиобластома (глиома 4 степени) не может быть удалена ни хирургически, ни с помощью химиотерапии. Опухоли могут давать рецидивы, поэтому важно все время проводить мониторинг состояния пациента и периодически фиксировать результаты диагностики, в том числе и с помощью МРТ. На МРТ снимках врач может как зафиксировать сам факт наличия опухоли и ее положение, так и процесс ее роста в течение жизни пациента.

1.1.2 Проблема сегментации опухоли

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

1.2 Анализ существующих решений

На данный момент уже существуют программные продукты, решающие выявленные проблемы. Условно можно разделить их на две большие группы: настольные программы для визуализации МРТ снимков и их обработки, и клиент-серверные приложения, предоставляющие удаленный доступ к базе данных снимков МРТ для рабочих станций в больнице. Так или иначе, все решения должны реализовывать стандарт DICOM (Digital Imaging and Communications in Medicine), который призван решить проблему по обеспечению унифицированного доступа к медицинской информации, ее передачи и обработки. Многие из решений являются либо платными, либо условно бесплатными, а также часто не имеют русской локализации. Также часто имеется проблема с переносимостью их с одной платформы на другую. В данной части мы рассмотрим два наиболее современных и продвинутых решений из числа свободно распространяемых и условно бесплатных.



1.2.1 OsiriX

OsiriX (рис. 1), как и все его производные продукты такие, как OsiriX-Lite и Horos, представляет собой настольный инструмент для просмотра, обработки и архивации DICOM файлов [22]. Также программы поддерживают 3D-реконструкцию серии изображений, поиска и загрузки снимков с PACS-серверов по протоколу DICOM. OsiriX и ее производные дают возможность врачу, в частности онкологу-диагносту, сегментировать МРТ снимки с помощью встроенного метода наращивания областей и выделять интересующие сегменты, находить их площадь или объем. Также программы позволяют расширить функциональность с помощью установки сторонних плагинов. Однако недостатками программы является отсутствие поддержки русского языка и невозможность запускать ее на операционных системах, отличных от Mac OS X. Более того, в частности, OsiriX очень требователен к техническим параметрам компьютера таким, как объем оперативной памяти. Плагины могут быть написаны только на языке Objective C/C++, что увеличивает входной порог для сторонних разработчиков.

Рисунок 1. Пример работы OsiriX



1.2.2 Orthanc

Orthanc (рис. 2) является свободно распространяемым PACS-сервером [21], полностью написанном на языке C++, что делает сервер чрезвычайно быстрым и производительным. Более того, он уже портирован на несколько операционных систем таких, как Windows, Linux и Mac OS X. Сам сервер может быть установлен путем компиляции из исходников в виде единственного бинарника или копирования готового бинарника в папку с исполняемыми файлами. Конфигурация сервера производится путем редактирования файла с параметрами в формате JSON. Функциональность Orthanc может быть расширена с помощью сторонних плагинов, написанных на C/C++. В стандартную поставку сервера идет удобный веб-интерфейсом для удаленной работой с PACS-системой. Также на официальном сайте есть плагины для добавления возможности работы с СУБД PosgtreSQL, протоколом DICOMWeb и библиотекой ImageJ. Однако минусом данной программы является то, что плагины пишутся на языке C/C++, что увеличивает порог вхождения для сторонних разработчиков. Интерфейс полностью на английском языке, не предусмотрено русской локализации. К тому же, не имеется встроенных инструментов сегментации и обработки МРТ снимков.

Рисунок 2. Просмотр тэгов DICOM в Orthanc

1.3 Постановка задачи на разработку

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

- Удаленный доступ к базе данных МРТ снимков по протоколам DICOM и DICOMWeb;

- Загрузка МРТ снимков в формате DICOM;

- Просмотр МРТ снимков с возможностью настройки визуализации (цветовая схема, поворот, увеличение, базовые фильтры);

- Поиск пациентов, обследований, серий обследований и отдельных изображений;

- Установка и удаление плагинов;

- Полуавтоматическая или автоматическая сегментация МРТ снимков для выявления опухолей и определения их размеров и положения.

Выводы по главе

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



Глава 2. Методы и алгоритмы сегментации раковых опухолей на МРТ снимках

В данной главе будут рассмотрены основные алгоритмы и методы обработки МРТ снимков и сегментации опухоли. Для каждого метода будут приведены достоинства и недостатки, а также математическое и алгоритмическое описание. В конце главы будут приведены результаты экспериментов для оценки производительности и качества работы описанных методов.

2.1 Предобработка и постобработка снимков

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

2.1.1 Удаление шума

Удаление шума может быть произведено двумя методами: гауссовым размытием (размытием по Гауссу) или медианным фильтром. Первый методы основан на нормальном распределении Гаусса и операции свертки. Свертка (ядро) представляет собой матрицу, умножая коэффициенты которой на матрицу изображения, можно получить желаемый результат. Для быстрых вычислений используется готовая матрица со значениями из выборки, полученной из нормального распределения. В результате применения размытие по Гауссу из изображения удаляются яркие зерна, т.е. шум. Достоинством данного метода является его высокая скорость, которая может быть значительно увеличена, если выполнять фильтр на графических процессорах. Однако минусом данного метода является то, что он плохо справляется с удалением ярких и острых краев. Медианный фильтр относится к нелинейным фильтрам в отличие от размытия по Гауссу. Более того, медианный фильтр решает проблему удаления острых краев и так называемого шума типа “перец и соль”. Суть фильтра заключается в том, чтобы все соседние пиксели вместе текущим пикселем заменить на медианное значение. Для вычисления медианны применяется техниках сортировки, при которой сортируются сначала строки, затем столбцы, затем главная диагональ, а после завершения сортировки в качестве медианны выбирается средний элемент на главной диагонали [28]. Как уже было сказано, медианный фильтр позволяет удалять шум "перец и соль", а также оставлять острые и яркие края, однако недостатком метода является низкая производительность по сравнению с размытием по Гауссу.

2.1.2 Удаление ярких краев

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

Края представляют собой те области на изображениях, при которых резко меняется значение яркости (интенсивности). Для определения, какие из пикселей относятся к краям, используются три основных метода нахождения краев: фильтр Лапласа, оператор Собеля и оператор Кэнни. В нашей программе мы используем оператор Собеля, так как он дает более выраженные и четкие края, чем фильтра Лапласа, и быстрее работает, чем оператор Кэнни. В операторе Собеля используется две матрицы для вычисления градиентов относительно оси абсцисс и оси ординат [35]:

-1 0 1 -2 0 2 -1 0 1	1 2 1 0 0 0 -1 -2 -1
Градиент по оси абсцисс GX	Градиент по оси ординат GY

Мы проходим по все пикселям с заранее заданным окном и для каждого из них находим всех соседей. Все найденные соседи и текущий пиксель образуют активный регион. Каждый пиксель из региона умножаем на значение из матрицы Собеля. Вычисленные произведения суммируем. Полученные суммы для оси абсцисс и ординат являются градиентами для текущего пикселя. После того, как были посчитаны градиенты для отдельных осей, мы должны найти величину градиента, которая может быть вычислена как [17]:

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



Рисунок 3. (1) - исходное изображение, (2) - обнаруженные края оператором Собеля, (3) - удаленные края из исходного изображения

2.1.3 Морфологические операции

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

2.1.5 Метод Отцу

Метод Отцу успешно применяется в задаче нахождения бинарного порога и широко используется в компьютерном зрении и обработки изображении. Метод был разработан японским ученым Nobyuki Otsu в 1979 [12]. Метод основан на предположении, что лучшее значение бинарного порога достигается только тогда, когда достигается минимальная внутриклассовая дисперсия. Nobyuki Otsu показал, что задача минимизации внутриклассовой дисперсии равносильна максимизации межклассовой дисперсии [12]. Рассмотрим алгоритм по шагам:

1. Пусть имеется изображение с уровнями интенсивности и размером ;

2. Мы хотим разбить изображение на два класса относительно бинарного порога : фон и перед, тогда фон и перед равны и ;

3. Вычислим гистограммы для изображения ;

4. Найдем вероятности попадания пикселей и в классы и и средние интенсивности и этих классов соответственно:

5. Найдем дисперсии и для классов и :

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

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

2.1.4 Удаление черепа

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

2.2 Сегментация с помощью пороговой фильтрации

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

При выполнении пороговой сегментации в нашей программе используется два предложения: раковая опухоль наиболее яркая область на снимке и головной мозг имеет линию симметрии. Исходное изображение головного мозга подвергается предобработки, при которой удаляется шум, яркие края и череп. Затем изображение разбивается на две части относительно линии симметрии головного мозга. В каждой из половин считается среднее значение яркости. Затем вычисляется первое опорное значение порога. Затем пиксели делятся на два класса: в первый попадают те, которые ниже порога, а во второй - выше порога. Затем снова считаются средние значения яркости для каждого из классов. Порог считается как среднее значение средних интенсивностей двух классов. Затем новое полученное значение порога сравнивается с пороговым значением на предыдущем шаге. Если разница между порогами лежит в пределах ошибки, то алгоритм останавливается работу, в противном случае - итерация продолжается. Достоинством данного метода является его простота реализации и высокая скорость выполнения. Недостатком данного метода является то, что он чувствителен к расположению опухоли (невозможно выявить опухоль, которая находится вблизи центрального канала, а также несколько очагов опухолей на разных участках мозга), случайному шуму и артефактам. Более того, метод разбивает изображение только на два сегмента, но в некоторых случаях для выявления опухоли нужно изображение разбивать больше, чем на два сегмента.

2.3 Сегментация с помощью кластеризации

Сегментация изображений МРТ может быть успешна решена с помощью алгоритмов кластеризации. При кластеризации мы пытаемся разбить изображение на несколько различимых по некоторому критерию кластеров, причем так, чтобы пиксели с схожими свойствами были как можно ближе друг к другу, а пиксели с разными характеристиками - как можно дальше. Условно методы кластеризации можно разделить на две большие группы: параметрические и непараметрические. В параметрических методах кластеризации указывается количество кластеров, а при непараметрических - количество кластеров вычисляется на основании анализа исходной выборки. В нашей программе реализованы три метода сегментации, основанных на кластеризации: кластеризация k-средних (K-Means), нечеткая кластеризация c-средних (Fuzzy C-Means) и кластеризация сдвигом по средней (Mean Shift). В данной части мы рассмотрим их математическую базу, алгоритм реализации, достоинства и недостатки.



2.3.1 Метод k-средних

Метод кластеризации k-средних впервые был предложен математиками Гуго Штейнгаузом и Стюартом Ллойдом в 1950-х прошлого века [30]. Метод k-средних считается одним из самых быстрых и популярных в задаче кластеризации, в том числе и в области сегментации изображений. Метод вычисляет евклидово расстояние от точек из выборки до ближайшего кластера, пытаясь минимизировать среднеквадратичную отклонение точек кластеров от центров этих кластеров. Другими словами, алгоритм находит минимум целевой функции:

Для понимания работы метода в задаче сегментации приведем алгоритм:

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

2. Зададим количество кластеров . Инициализируем центроиды кластеров случайными значениями из нормального распределения;

3. Вычислим матрицу расстояний между пикселями и центроидами . Так как мы имеем одномерное пространство, то для увеличения производительности будем считать не квадрат расстояния, а модуль разности между интенсивностью пикселя и центроида. Матрица расстояний может быть вычислена как:

,

где ;

4. Затем для каждого пикселя найдем тот кластер, расстояния до которого наименьшее, то есть найти вектор меток :

5. Вычисляем значение целевой функции как сумма модулей разностей между центроидами и точками из выборки:

;

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

7. Выбираем кластер с наибольшим значением яркости. Этот кластер и будет регионом опухоли:

8. Применяем бинарную пороговую фильтрацию: всем пикселя, относящимся к кластеру, присваиваем максимально возможную яркость, а остальным - нулевую:

Результат работы алгоритма может быть продемонстрирован на следующем примере сегментации головного мозга с выделенной опухолью (рис. 4):

Рисунок 4. (1) - исходное изображение, (2) - после предобработки, (3) - результат кластеризации, (4) - конечный результат сегментации (опухоль красного цвета)

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

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

2. Зададим количество кластеров . Инициализируем центроиды кластеров случайными значениями из нормального распределения;

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

4. Вычислим матрицу расстояний между пикселями и центрами кластеров . Матрица расстояний может быть вычислена как:

раковая опухоль снимок сегментация

,

где ;

5. Затем для каждого пикселя найдем тот кластер, расстояния до которого наименьшее, то есть найти вектор меток :

6. Вычисляем значение целевой функции как квадратичную сумма разностей между центрами и точками из выборки:

;

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

8. Выбираем кластер с наибольшим значением яркости. Этот кластер и будет регионом с опухолью:

9. Применяем бинарную пороговую фильтрацию: всем пикселя, относящимся к кластеру, присваиваем максимально возможную яркость, а остальным - нулевую:

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

Рисунок 5. (1) - исходное изображение, (2) - 2 кластера, (3) - оптимальный вариант при 3 кластерах, (4) - 7 кластеров



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

2.3.2 Метод нечеткой кластеризации c-средних

Метод нечеткой кластеризации c-средних очень похож на метод кластеризации k-средних за тем исключением, что при нечеткой кластеризации точки из выборки (пиксели на изображении) относятся к разным кластерам с разной степенью принадлежности [6]. Данный метод является более робастным на шумных изображениях, имеющих ярко выраженные отклонения, чем алгоритм четкой кластеризации k-средних. Так же, как и метод k-средних, нечеткая кластеризация является параметрическим методом, требующим указания числа кластеров [6]. Однако, в отличие от k-средних, нечеткая кластеризация зависит также и параметра "нечеткости" (параметра m). Параметр m влияет на "размытость" результата и устойчивость к отклонениям. При m, близкому к 1, результат работы метода почти ничем не отличается от результата работы четкой кластеризации. Нечеткая кластеризация пытается оценить коэффициенты принадлежности точек к кластерам таким образом, чтобы минимизировать квадратичную сумму отклонений этих точек от кластеров. Другими словами, метод минимизирует следующее выражение [6]:



Опишем работу алгоритма кластеризации c-средних в задаче сегментации изображений:

1. У нас имеется чёрно-белое изображение МРТ снимка. Создадим множество однокомпонентных векторов со значением яркости пикселей длиной, равной произведению ширины изображения на высоту:

.

2. Зададим число кластеров и некоторый коэффициент .

3. Инициализируем матрицу принадлежности случайными значениями из нормального распределения так, что сумма по строкам была равна 1:

4. Вычислим центры кластеров следующим образом:

5. Вычислим целевую функцию:

6. 

7. Если модуль разницы между значением целевой функции на текущем шаге и значением целевой функции на предыдущем шаге меньше заданной ошибки, то останавливаем алгоритм кластеризации и переходим к шагу 8. В противном случае мы переходим к шагу 7. Критерий останова может быть вычислен как:

8. Пересчитываем матрицу принадлежности и переходим обратно к шагу 4:

9. После завершения кластеризации производим дефаззификацию индексов кластеров для пикселей с помощью метода FOM (First of Maximum):

10. Выбираем кластер с наибольшим значением яркости. Этот кластер и будет регионом с опухолью:

10. Применяем бинарную пороговую фильтрацию: всем пикселя, относящимся к кластеру, присваиваем максимально возможную яркость, а остальным - нулевую:

Как уже было сказано выше, нечеткая кластеризация позволяет получать более точную сегментацию по сравнению с методом k-средних, особенно на шумных изображениях с большим количеством отклонений. Но за точность приходится платить производительностью - для нечеткой кластеризации требуется умножение матриц, а также дополнительные действия для определения индекса кластера для пикселя, что увеличивает время подсчета. Так или иначе, оба метода требуют указание количества кластеров. Более того, методы не способны выявлять опухоль на МРТ снимках типа T1, на которых область опухоли имеет невысокие значения интенсивности (опухоль преимущественно темнее других тканей на МРТ снимке).

2.3.3 Метод сдвига среднего

Для решения проблемы выбора числа кластеров существует метод сдвига среднего (Mean Shift). Данный метод строит для исходной выборки функцию плотности и находит максимумы этой функций для областей с высокой плотностью точек. Впервые метод был предложен учеными K. Fukunaga и L.D.Hostetler в 1975, а затем развит Y. Cheng в 1995 [5]. В основе метода лежит так ядра (kernel), которые позволяют построить непараметрическую оценку плотности распределения (kernel density estimation). В качестве ядер могут быть использованы простые функции, так и нелинейные функции (функция Гаусса, ядро Епанечникова). Примеры ядер представлены внизу (первая формула относится к плоскому ядру, вторая - к ядру Гаусса):



Чтобы функция была ядром, она должна удовлетворять следующим условиям [5]:

1. Ядро должно быть неотрицательным, то есть ;

2. Ядро должно быть невозрастающим, то есть ;

3. Ядро должно быть кусочно-непрерывной функцией, причем

Распишем работу алгоритма в задачи сегментации изображения:

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

2. Введем оценочную функцию распределения:

;

3. Найдем средневзвешенную среднюю функции плотности в точки с заданным окном :

4. Введем понятие сдвига по средней (сдвига среднего) ;

5. Вектор сдвига среднего совпадает с максимальным ростом плотности в точке ;

6. Кластеры соответствуют локальным максимумам (модам) функции плотности;

7. Вычисляем градиент для оценочной функции плотности;

8. Итеративно присваиваем каждому пикселю значение двигаемся по градиенту до тех пор, пока не будет удовлетворено условие останова [19];

9. После завершения нахождения максимумов объединим пришедшие в одну окрестность локального максимума пиксели в один сегмент (кластер) [37].

10. Завершаем алгоритм. Мы получили кластеры, которые и являются результатом сегментации. Закрашиваем каждый кластер в свой цвет и производим смешивание цветов сегментации с цветами исходного изображения. Результат сегментации представлен на рис. 6:

Рисунок 6. Результат сегментации с помощью Mean Shift (самые яркие кластеры относятся к опухоли)

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

2.4 Сегментация с помощью Гауссовых смесей

Для сегментации опухоли может быть применена техника, основанная на так называемых Гауссовых смесях или моделях Гауссовых смесей. Гауссова смесь представляет собой взвешенную сумму двух или нескольких Гауссовых распределений [27]. Каждое отдельное распределение представляет собой кластер. Гауссовы смеси могут справляться лучше с нелинейными структурами в данных, чем k-средних. Более того, метод, основанный на модели гауссовой смеси, учитывает разброс данных и пытается максимально точно описать "форму" кластеров. Кластеризация с помощью Гауссовых смесей представляет собой EM-алгоритм по аналогии с кластеризацией k-средних и состоит из двух шагов [27]: вычисление ожидаемого значения функции правдоподобия (E-шаг) и оценка ожидаемого значения функции правдоподобия (M-шаг). Но вместо того, чтобы вычислять центры кластеров как среднее и присваивать метки отдельным пикселям, мы будем определять параметры каждого из распределений, входящих в Гауссову смесь. Рассмотрим алгоритм по шагам:

1. Пусть имеется изображение , где ;

2. Пусть есть Гауссова смесь , где ;

3. Инициализируем параметры ;

4. Пусть имеется логарифмическая функция правдоподобия:



5. Вычислим новые значения весов компонент для пикселя (E-шаг);

6. Оценим параметры (вес, среднее и дисперсию) функции правдоподобия (M-шаг):

7.

8. Вычислим логарифмическую функцию правдоподобия:

9. Если модуль разности между значением функции правдободия на текущем шаге и значением функции правдоподобия || на предыдущем шаге меньше некоторого заданного , то переходим к шагу 9, в противном случае - к шагу 5;

10. Вычислим метки для каждого пикселя:

11. Найдем индекс региона опухоли:

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

2.5 Сегментация с помощью расширяющихся областей

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

1. Выберем опорный пиксель и зададим значение порога ;

2. Создадим пустой стек и положим в него опорную точку;

3. На этапе итерации вытащим элемент из стека и найдем всех его соседей ;

4. Найдем модуль разности яркостей текущей точки и его соседа. Если мы еще не просматривали соседа и модуль разности меньше, чем порог, то помечаем соседа как пройденного и кладем его в стек.

5. Возвращаемся к шагу 3 и продолжаем итерацию до тех пор, пока стек не опустеет.

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

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

2.6 Общие вычисления на GPU

Для выполнения методов обработки изображения и сегментации требуется много времени. Современные графические процессоры позволяют сократить время работы и распараллелить алгоритмы. На данный момент существует две основные технологии программирования на графических картах: CUDA и OpenCL. CUDA предоставляет удобный API для работы с графическими картами компании Nvidia, архитектуру параллельных вычислений и набор готовых библиотек. Достоинством данной технологии является то, что она хорошо оптимизирована для работы на GPU Nvidia, а также богатый набор библиотек для решения задач линейной алгебры, генерации псевдослучайных чисел и Фурье-преобразования. Однако данная технология может быть использована только на видеокартах Nvidia. Другая технология OpenCL дает возможность программисту создавать программы для видеокарт разных производителей, в том числе и Nvidia и AMD. Более того, OpenCL может задействовать вычислительные мощности центрального процесса, что позволяет создавать гетерогенные программные комплексы для проведения параллельных вычислений общего назначения. Данная технология полностью открытая и имеет хорошую спецификацию. Но не смотря на все плюсы, OpenCL хуже оптимизирован для видеокарт Nvidia, поэтому в бенчмарках программы, написанные на OpenCL, работают медленнее, чем на CUDA [11].

Наша программа не зависит от типа карты, так как реализует механизм позднего связывания через систему плагинов. Каждый из плагинов может быть реализован как на CPU, так и на GPU. Единственное ограничение заключается в том, на компьютере пользователя или сервер должны быть установлены современные графические карты и соответствующий набор библиотек и исполняемая среда. В рамках выпускной квалификационной работы было решено использовать CUDA 9.0. Тестирование алгоритмов производилось на видеокарте Nvidia GeForce GT 650M с размером видеопамяти 512 Мб.

Для графических карт пишутся специальные программы, называемые вычислительными ядрами или просто ядрами [14,36], на языке программирования, являющегося расширением C++. Особенность данных программ заключается в том, что один и тот же код выполняется на нескольких ядрах или параллельных узлах внутри графического процессора (принцип SIMT - Single-Instruction, Multiple-Thread) [14,36]. Язык программирования C++ CUDA помимо стандартных средств языка C++ предоставляет специальные ключевые слова и функции. Каждое ядро представляется собой обычную C++ функцию, помеченную ключевым словом “__global__” [36]. Ядро может вызывать и другие ядра или локально определенные функции на GPU, помеченные “__device__” [36]. Основное устройство или хост может запускать только ядра с CPU. Каждый вызов ядра асинхронен и записывается в очередь вызовов на GPU. Сами ядра не возвращают результат, т.е. являются “void”-функциями. На вход они могут принимать как примитивные типы (int, float, double), так и указатели, массивы и структуры. Ниже представлено определение функции-ядра для вычисления суммы двух векторов (рис. 7):

Рис. 7 - Определение функции-ядра для вычисления суммы векторов

В данной функции (рис. 7) для каждого отдельного потока вычисляется сумма двух векторов, которая присваивается третьему вектору. Позиция в векторе вычисляется как сумма номера потока (thread) и позиция блока потоков (blocks) в сетке (grid).

CUDA скрывает внутреннюю реализацию графического процессора и предоставляет трехуровневую абстракцию “поток-блок-сетка” [36]. Каждый поток представляет собой единицу параллельного выполнения. Максимальное количество потоков варьируется от модели к модели GPU, но стандартное значение равно 1024. Потоки объединяются в блоки потоков или просто блоки. С точки зрения внутренней реализации каждый блок представляет собой отдельный мультипроцессор. Максимальное количество блоков в каждой карте равно 65536. Все блоки объединяются в сетку. С точки зрения реализации алгоритмов, поток представляет собой столбец в векторе или матрице, а блок - строку. Все потоки могут представлять собой одномерные, двумерные и трехмерные вектора. Блоки могут быть как одномерными, так и двумерными. Для каждого потока в блоке и блока в сетке можно узнать их индекс. Также имеется возможность получить длину блока (число потоков) и длину сетки (число блоков). Трехуровневая модель представлена на рис. 8:



Рисунок 8. Архитектура CUDA

Помимо трехуровневой модели потоков CUDA предоставляет различные типы памяти [18]: регистровая, локальная (приватная), разделяемая, глобальная, константная и текстурная. Также GPU предоставляет два типа кеша: L1-кеш на уровне одного мультипроцессора и L2-кеш на уровне глобальной памяти. Среднестатистические карты, стоящие на настольных компьютерах, не могут обращаться напрямую к DRAM, однако последние модели Nvidia располагают возможностью работы с основной памятью компьютера. Самая быстрая память считается регистровая, но она недоступна напрямую программисту. Локальная память выделяется на каждом потоке и обычно равна одному или нескольким килобайтам. Разделяемая память доступна потокам из блока и обычно равна от 4 до 48 килобайт [16]. Глобальная память доступна всем блокам из вычислительной сетки и имеет объем, равный объему видеопамяти. Также имеется возможность оптимизировать частые обращения к памяти с помощью константной памяти, которая имеет размер, ограниченный от несколько десятков до несколько сотен килобайт. Текстурная память применяется в задачах обработки изображений из-за частого обращения к одним и тем же участкам памяти. Использование текстурной памяти особенно полезно, когда поток обращается к близким участкам памяти по отношению к той памяти, к которой привязан поток. Текстурная память, как и глобальная ограничена размером видеопамяти, но за счет кеширования обращение к ней происходит быстрее. Однако текстурная память является только читаемой и недоступна для чтения. В последних версиях CUDA имеется возможность для графических процессоров обращаться к данным непосредственно через оперативную память. Данная память называется закрепленной или общей. Подобный тип памяти позволяет обмениваться данными между CPU и GPU напрямую, минуя копирования в и из видеопамяти. Для каждого из типов памяти для переменных предусмотрены различные модификаторы, которые извещают компилятор NVCC, как переменные должны быть расположены в памяти. Ниже представлен список типов памяти их характеристики (табл. 1):