Таблица 1

Тип памяти	Модификатор	Размер	Назначение
Локальная или приватная	Нет	Несколько килобайт	Доступна отдельному потоку
Разделяемая	__shared__	От 4 до 48 килобайт	Доступна всем потокам из блока
Глобальная	__device__ или нет	Максимальный размер видеопамяти VRAM	Доступна всем блокам
Константная	__constant__	От нескольких десятков до нескольких сотен килобайт	Полностью читаемая кэшируемая недоступная для записи память
Текстурная	Нет	Максимальный размер видеопамяти VRAM	Полностью читаемая кэшируемая недоступная для записи память, оптимизирована для задач обработки изображений
Общая или закрепленная	Нет	Размер RAM	Доступна как CPU, так и GPU



Из-за параллельного обращения к глобальной и разделяемой памяти выполнение задач на GPU затруднено проблемой состояния гонки (race condition), при котором каждый из потоков одновременно работает с одним участком памяти и конкурирует за право записи в него или чтение из него. Подобная неопределенность ведет к непредсказуемым результатам. Для синхронизации потоков внутри блока используется специальная функция __syncthreads() [36], которая приостанавливает поток в текущем блоке, если все остальные не закончили работу. Данная функция служит своего рода барьером для уже закончивших работу потоков и не дает им продолжить работу, пока остальные потоки не закончили задачу. Данный механизм синхронизации не работает на уровне всех блоков внутри сетки, так как реализация синхронизации на уровне всего GPU очень затратно и неэффективно [36]. Данная проблема особо остро встает при реализации алгоритмов редукции. Однако она может быть решена с помощью замены одной редукции последовательностью частичных редукций. На каждом этапе мы выбираем кусочек данных и производим редукцию на нем в пределах одного блока. Затем для всех частичных редукций вычисляем общий результат путем глобальной редукции на уровне всего GPU. Глобальная редукция производится только в том случае, если размер массива результатов частичных редукций меньше, чем размер одного блока. Таким образом, мы добиваемся высокой скорости выполнения задачи без риска испортить данные непредсказуемостью порядка параллельного обращения к памяти. Рассмотрим данный подход на примере реализации метода сегментации k-средних. В ходе анализа алгоритма, было выяснено, что метод может быть разделен на три параллельных блока: нахождение индекса кластера для каждого пикселя, подсчета суммы и количества элементов кластера в отдельном блоке и вычисление центров. Для каждого из этих блоков были разработаны функции-ядра. Рассмотрим все этапы отдельно. Первый этап состоит из параллельного поиска для каждого пикселя ближайшего кластера. Приведем работу алгоритма и код реализации (рис. 9):

1. Пусть мы имеем номер потока, номер блока и размер блока , и соответственно. Вычислим позицию в исходном изображении ;

2. Если позиция больше, чем размер изображения, то останавливаем работу ядра, в противном случае - переходим к шагу 3;

3. Пусть минимальное расстояние равно 1. Для каждого пикселя ищем ближайший кластер. Для этого мы пробегаемся по всем кластерам, вычисляем Эвклидово расстояние от кластера до пикселя. Если текущее расстояние меньше минимального, то считаем его за минимальное и обновляем номер найденного кластера. Номер найденного кластера записываем в , где ;

4. Завершаем алгоритм. Мы получили номера кластеров для всех пикселей.

Рис. 9 - Вычисление ближайшего кластера

На следующем шаге мы частично вычисляем сумму всех элементов, входящих в кластер, и их число. Приведем работу алгоритма и код реализации (рис. 10):

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

2. Выделим разделяемую память размером байт, где , ;

3. Инициализируем разделяемую память нулями. Если текущий пиксель относится к заданному кластеру, то в разделяемую память и записываем его яркость и единицу соответственно;

4. Приостановим работу ядра, пока все потоки не закончат этап 3;

5. Пусть смещение равно .

6. Последовательно для каждого потока считаем сумму элемента из разделяемой памяти с индексом, равного номеру потока, и элемента, имеющего индекс, равный номера потока плюс смещение, то есть:

7. Уменьшаем смещение в два раза и возвращаемся к этапу 6, пока смещение не станет нулевым, в противном случае - переходим к шагу 8.

8. Если номер потока равен нулю, то записываем в результат частичных редукций , где , , то есть:



9. Инкрементируем номер кластера. Если номер кластер меньше числа кластеров, переходим к шагу 1, иначе заканчиваем работу.

Рисунок 10. Вычисление центров (частичная редукция)

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

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

2. Выделим разделяемую память размером байт, где , ;

3. Запишем в разделяемую память результат частичной редукции, то есть:

4. Приостановим работу ядра, пока все потоки не закончат шаг 3.

5. Пусть смещение равно .

6. Последовательно для каждого потока считаем сумму элемента из разделяемой памяти с индексом, равного номеру потока, и элемента, имеющего индекс, равный номера потока плюс смещение, то есть:

7. Уменьшаем смещение в два раза и возвращаемся к этапу 6, пока смещение не станет нулевым, в противном случае - переходим к шагу 8.

8. Если номер потока равен нулю, то вычисляем центр кластера, то есть

9. Инкрементируем номер кластера. Если индекса кластер меньше числа кластеров, переходим к шагу 1, иначе заканчиваем работу.



Рисунок 11. Вычисление центров (полная редукция)

Когда завершены все итерации, мы получаем оптимальную кластеризацию пикселей в исходном изображении. Таким образом, за счет особенностей алгоритма кластеризации k-средних и возможностей графического процессора задачу сегментации можно распараллелить и повысить скорость вычислений в несколько раз по сравнению с вычислениями на центральном процессоре. Однако при выполнении алгоритма на GPU нужно избегать частых обменов между ним и CPU, так как процесс копирования занимает много времени и снижает эффективность работы методов. Чтобы не производить частые копирования, нужно за раз загрузить все данные, выполнить действия над ними, и только после полного завершения работы производить обратный обмен с CPU. В случае с методом k-средних мы пропустили этап вычисления целевой функции и преждевременного останова метода. Вместо того, чтобы все время вычислять целевую функцию и ошибку, а затем копировать на CPU и принимать решение о продолжении работы, мы делаем столько итераций, сколько задается пользователем. Поэтому мы делаем только два обмена. При первом копируем в VRAM все изображения, изначальные центры и параметры, а при втором - метки и центры обратно на CPU.

Так как управление памятью GPU напрямую недоступно CPU, то для ее выделения и инициализации применяются специальные функции cudaMalloc и cudaMemcpy. Первая функция выполняет ту же функцию, что и malloc на CPU. Она принимает два аргумента: указатель на переменную и размер переменной в байтах. Вторая функция производит копирования из одного участка памяти в другой. Она поддерживает три режима копирования: от центрального процессора на графическую карту, от графической карты на центральный процессор и между блоками в графическом процессоре. Функция принимает четыре аргумента: указатель на получателя, указатель на отправителя, размер в байтах и режим копирования. Помимо выделения памяти, как и в CPU, на GPU нужно очищать выделенную память для предотвращения утечек. Для этого предусмотрена функция cudaFree, которая принимает единственный аргумент - указатель на очищаемую память.

Важно заметить, что выполнение ядер происходит асинхронно, но последовательно в пределах одного потока (stream). Соответственно CPU не ожидает завершения ядра, однако при копировании данных с GPU центральный процессор блокируется до тех пор, пока все ядра не выполнят свою работу и не вернут данные обратно CPU. В некоторых случаях требуется ручная блокировка центрального процессора, поэтому для этого предусмотрена функция cudaDeviceSynchronize. В силу асинхронности ядер нельзя измерить время их выполнения с помощью стандартных средств измерения времени на CPU [15,36]. Для этого предусмотрен механизм так называемых событий, при котором в определенные промежутки времени GPU записываются события в журнал [15,36]. Затем после завершения всех ядер вычисляется промежуток времени между первым и последним событиями. Для создания и удаления событий предусмотрены функции cudaEventCreate и cudaEventDestroy [15,36]. Для вычисления времени между двумя событиями используется функция cudaEventElapsedTime [15,36], возвращающая время в миллисекундах. В следующей табл. 2 представлены все определения функций, описанных выше:



Таблица 2

Функция	Параметры	Возвращаемое значение	Назначение
cudaMalloc()	void** src - указатель на переменную, size_t - размер в байтах	cudaError_t - код ошибки	Выделение памяти
cudaMemcpy()	void* dst - указатель на получателя, const void* src - указатель на отправителя size_t - размер в байтах, cudaMemcpyKind - режим обмена	cudaError_t - код ошибки	Обмен данных
cudaFree()	void* src - указатель на переменную	cudaError_t - код ошибки	Очищение памяти
cudaDeviceSynchronize()	Нет	cudaError_t - код ошибки	Синхронизация центрального процессора
cudaEventCreate()	cudaEvent_t* - указатель на хендлер события	cudaError_t - код ошибки	Создание события
cudaEventDestroy()	cudaEvent_t - хендлер события	cudaError_t - код ошибки	Удаление события
cudaEventElapsedTime()	float* msc - указатель на переменную, которая хранит время, cudaEvent_t start - хендлер события, cudaEvent_t end - хендлер события	cudaError_t - код ошибки	Вычисление времени между двумя событиями

2.7 Оценка качества сегментации

Качество и производительность методов сегментации варьируется в зависимости от алгоритма, реализации, вычислительной системы и параметров. В большинстве случаях мы хотим узнать, насколько алгоритм сегментации справился со своей задачей, насколько он хорошо предсказывает опухоль и как долго выполняется. Обычно за истинный регион опухоли принимают ручную сегментацию врача (Ground Truth). Большая проблема в оценке качества и интерпретации результата заключается в выборе метрики оценивания. В задаче оценки качества сегментации принято использовать следующие метрики: коэффициент Dice-Sшrensen, индекс Jaccard Index, F1-Score, ROC-кривая, кривая Precision-Recall. В нашей программе мы используем коэффициент Dice-Sшrensen. Данная метрика показывает, насколько хорошо наш алгоритм угадывает область опухоли. Метрика вычисляется, как доля мощности удвоенного пересечения двух классов к сумме мощностей самих классов. Чем ближе данный коэфициент к 1, тем лучше наш алгоритм угадывает область опухоли. Коэффициент может быть вычислен как:

Также при оценке качества работы алгоритмов необходимо производить измерение времени работы. Каждый алгоритм прогоняется по несколько раз на каждом снимке из исходной выборки, а затем вычисляется среднее время работы алгоритма как:

.

Для вычисления 95% доверительного интервала для среднего с неизвестной дисперсией должно быть использовано распределение Стьюдента. Для величины со средним и вычисленной дисперсией 95%-интервал равен:

где .

Для проведения оценки качества работы методов были выбраны 5 алгоритмов: кластеризация k-средних (K Means), многомерная кластеризация k-средних (Multivariate K Means), нечеткая кластеризация c-средних (Fuzzy C Means), Гауссова смесь (Gaussian Mixture) и кластеризация k-средних на GPU (CUDA K Means). В качестве выборки были использованы снимки 44 пациентов с глиомой 2 степени размером 256 на 256 пикселей, взятые с TCIA [15]. Все замеры производились на Interl Core i7 с тактовой частотой 2,3 ГГц (CPU) и Nvidia GeForce GT 650M (GPU) с видеопамятью объемом 512 Мб и расширяемой памятью объемом 48 Кб. Оценки производились на разном числе числе кластеров: от 3 до 9. Для подсчета среднего значения времени и точности для каждого снимка производилось 10 итераций. Таким образом, производилась оценка времени и точности работы одного алгоритма в зависимости от числа кластеров. В итоге были получены следующие результаты для времени в миллисекундах (табл. 3) и точности в процентах (табл. 4):

Таблица 3

Алгоритм	Число кластеров
	3	4	5	6	7	8	9
K Means	32,53	48,81	60,54	84,75	98,05	133	120,53
Fuzzy C Means	672,44	1214,11	1646,82	2100,9	2505,93	2878,68	3215,06
Gaussian Mixture	253,99	273,08	327	363,89	427,29	489,68	509,64
Multivariate K Means	67,82	99,01	171,89	184,54	273,94	365,63	1189,01
CUDA K Means	21,13	25,19	30,78	35,78	41,2	46,52	52,01

Таблица 4

Алгоритм	Число кластеров
	3	4	5	6	7	8	9
K Means	66,72	79,7	78,17	75,17	70,2	63,55	59,55
Fuzzy C Means	63,96	73,91	76,71	75,52	73,94	71,2	68,65
Gaussian Mixture	58,03	79,04	77,38	74,35	69,42	64,56	59,95
Multivariate K Means	72,5	78,36	77,18	74,2	68,53	63,7	60,63
CUDA K Means	69,82	80,09	77,49	72,11	66,73	60,67	55,67

После получения первичных результатов был вычислен 95% доверительный интервал для времени (табл. 5) и точности (табл. 6):

Таблица 5

Алгоритм	Число кластеров
	3	4	5	6	7	8	9
K Means	[24.74, 40.33]	[38.90, 58.72]	[47.51, 73.57]	[70.98, 98.52]	[81.36, 114.74	[89.00, 177.00	[96.06, 144.99
Fuzzy C Means	[555.01, 789.8	[965.10, 1463.	[1383.32, 1910	[1712.07, 2489	[2034.09, 2977	[2333.74, 3423	[2640.47, 3789
Gaussian Mixture	[209.20, 298.7	[215.17, 330.9	[273.07, 380.9	[303.00, 424.7	[345.87, 508.7	[402.42, 576.9	[420.45, 598.8
Multivariate K Means	[51.34, 84.30]	[69.77, 128.26	[93.44, 250.33	[145.33, 223.7	[193.48, 354.3	[250.90, 480.3	[0.00, 2878.78
CUDA K Means	[17.91, 24.36]	[21.40, 28.98]	[26.18, 35.37]	[30.39, 41.18]	[35.01, 47.38]	[39.52, 53.53]	[44.16, 59.85]



Таблица 6

Алгоритм	Число кластеров
	3	4	5	6	7	8	9
K Means	[57.26, 76.18]	[74.76, 84.65]	[73.71, 82.63]	[70.71, 79.63]	[65.95, 74.44]	[58.90, 68.19]	[55.14, 63.95]
Fuzzy C Means	[54.04, 73.89]	[66.18, 81.65]	[71.39, 82.03]	[71.02, 80.02]	[69.66, 78.23]	[67.14, 75.26]	[64.68, 72.62]
Gaussian Mixture	[48.93, 67.14]	[74.91, 83.16]	[72.87, 81.89]	[69.71, 78.99]	[64.79, 74.05]	[60.02, 69.11]	[55.54, 64.36]
Multivariate K Means	[64.75, 80.25]	[73.42, 83.30]	[72.45, 81.91]	[69.08, 79.32]	[63.16, 73.91]	[58.70, 68.69]	[55.52, 65.73]
CUDA K Means	[61.66, 77.99]	[75.25, 84.94]	[72.85, 82.13]	[67.01, 77.22]	[61.81, 71.65]	[55.82, 65.52]	[51.17, 60.17]

Также были построены графики зависимости времени работы (рис. 12) и точности алгоритма от числа кластеров (рис. 13) :

Рисунок 12. Зависимость времени от числа кластеров



Рисунок 13. Зависимость точности от числа кластеров

Для сравнения производительности реализаций метода кластеризации k-средних на CPU и GPU были построены отдельные графики зависимости времени (рис. 14) и точности от числа кластеров (рис. 15)

Рисунок 14. Графики зависимости времени от числа кластеров для k-средних на CPU и GPU

Рисунок 15. Графики зависимости точности от числа кластеров для k-средних на CPU и GPU

Исходя из полученных результатов можно убедиться, что использование графических карт позволяет ускорить выполнение алгоритма кластеризации k-средних почти в 3 раза, что подтверждается в табл. 3, табл. 5 и на рис. 14. Маленькая разница между времени исполнения метода кластеризации на CPU и GPU при числе кластеров 3 можно объяснить большими задержками в обмене данных между VRAM и DRAM, которые занимают большую часть времени работы алгоритма на GPU. Однако, как можно убедиться исходя рис. 14, при увеличении числа кластеров использование GPU позволяет получить прирост в скорости, причем при самом увеличении числа кластеров показатель времени работы алгоритма на GPU растет медленнее, чем на CPU. Худшие показатели производительности показывает метод сегментации, основанный на нечеткой кластеризации. В среднем данный метод работает от 672.44 (число кластеров 3) до 3215.06 миллисекунд (число кластеров 9). Это может быть объяснено дорогостоящими операциями умножения матриц и дополнительным этапом по вычислению центров кластеров через деффазификацию. Использование Гауссовых смесей дает прирост по скорости по сравнению с нечеткой кластеризацией в среднем в 2-7 раз. Использование k-средних дает прирост по скорости по сравнению с нечеткой кластеризацией в среднем в 21-27 раз. Несмотря на различия в скорости работы методов, показатель точности у всех методов держится на уровне 55%-81%. Максимальный показатель точности достигаются при числе кластеров 4. При 4 кластерах мы получили в среднем 78.22% по метрике Dice-Sшrensen. Также, можно заметить исходя из графика рис. 13, что при увеличении числа кластеров точность имеет тенденцию падать. Графики всех методов, кроме нечеткой кластеризации c-средних, имеют схожий характер падения при увеличении числа кластеров. Также можно заметить, что график метода нечеткой кластеризации c-средних имеет более пологую правую часть, чем у всех остальных графиков. Чтобы объяснить особенность поведения метода нечеткой кластеризации и его отличия от других методов сегментации, были принято решение провести дополнительные эксперименты.

В отличие от остальных методов, нечеткая кластеризация c-средних зависит не только от числа кластеров, но и от параметра "нечеткости" или "размытия" (параметра m). При m, близкому к 1, результат становится более "четким" и неотличимым от результата работы метода k-средних. При m, близким к бесконечности, результат становится более "размытым" и "нечетким". Более того, при увеличении параметра m алгоритм нечеткой кластеризации становится более устойчивым к ошибкам и выбросам, что делает график зависимости точности от числа кластеров более пологим при большим числе кластеров. Для определения влияния параметра m на результат сегментации, были произведены дополнительные эксперименты, в ходе которых были получены данные влияния числа кластеров и параметра m на точность выходной сегментации. Все замеры производились на том же оборудовании и той же выборке. После получения первичных результатов были построены графики зависимости точности от числа кластеров для каждого значения параметра m (рис. 16).

Рисунок 16. Графики зависимости точности от параметра m

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

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

Таблица 7.

	Пространственная сложность	Временная сложность
Кластеризация k-средних
Нечеткая кластеризация c-средних
Бинарный порог
Гауссовы смеси

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

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

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

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



Глава 3. Особенности программной реализации

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

3.1 Требования к программе

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

- Программа должна уметь загружать и хранить МРТ снимки в формате DICOM;

- Программа должна уметь показывать МРТ снимки с возможность настраивать параметры визуализации;

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

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

- Программа должна уметь сегментировать опухоль на снимке МРТ;

- Программа должна загружать истинную сегментацию для оценки качества работы алгоритмов сегментации.

3.2 Инструменты разработки (сборки) и среда исполнения программы

В данной части будут рассмотрены инструменты для разработки, отладки и сборки программы, а также программное окружение для обеспечения ее исполнения.

3.2.1 Среда разработки

Для разработки программы было решено использовать следующий набор технологий:

- Python 3 [25] - мультипарадигменный криптовый язык;

- Jetbrains PyCharm [9] - среда разработки под Python;

- Jetbrains Webstorm [10] - среда разработки веб-приложений;

- Nvidia Nsight [18] - среда разработки GPU-программ для видеокарты Nvidia с использование CUDA;

- Clang G++ [2] - компилятор языка программирования C/C++;

- NVCC [19] - специализированный компилятор программ на языке CUDA C/C++ для графических процессор Nvidia;

- Make [8] - система сборки программ на C/C++;

- CMake [3] - система сборки программ на C/C++;

- Django [7] - веб-фреймворк для создания веб-приложений на Python;

- Numpy [13] - библиотека для работы с многомерными массивами;

- OpenCV [20] - библиотека для обработки изображений;

- CUDA [17] - библиотека для создания программ для графических карт Nvidia;

- Cython [4] - компилятор расширений на C/C++ для программ на Python;

- Tornado [33] - асинхронный неблокирующий веб-сервер;

- pydicom [24] - библиотека для чтения данных из DICOM-файлов;

- Scikit-Learn [29] - библиотека, содержащая реализации алгоритмов сегментации;

- PyCUDA [32] - обертка библиотеки CUDA для программ на языке Python;

- PostgreSQL [23] - реляционная база данных для хранения информации о снимках.

Язык программирования Python 3 был выбран, потому что является наиболее удобным вариантом для написания сложной системой, а также то, что для этого языка предусмотрена большое количество библиотек для решения поставленных задач (сегментация, обработки изображений, сетевое взаимодействие). Однако в будущем, когда программа пройдет апробацию, основную часть программы следует переписать на более быстрый язык программирования C++. Также язык Python 3 подходит для разработки большего числа сторонних плагинов, которые могут быть легко внедрены в программу.

3.2.2 Веб-сервер

Так как программа является клиент-серверной, то было решено использовать легкую в обслуживании и в то же время быструю и асинхронную библиотеку Tornado. В нашей программе мы решили произвести скрещивание Django, Django RestFramework и Tornado. Django и Django RestFramework предоставляет следующий набор возможностей:

- Объектно-реляционное отображение;

- Легкое и быстрое построение REST API;

- Панель администрирования;

- Механизм встроенных команд и возможность написания своих;

- Сериализация данных в JSON;

- Автоматическая конфигурация окружения.

Tornado представляет собой сетевую библиотеку для создания быстрых и асинхронных серверных приложений. При запуске приложения Tornado создает несколько процессов, каждый из которых обрабатывает свой набор клиентских запросов. Также библиотека позволяет работать с низкоуровневым протоколом TCP/IP, что полезно при работе с протоколом DICOMNet. Запросы в библиотеки могут быть как синхронными, так и асинхронными. Последние особенно полезны, когда процесс обработки данных занимает больше времени, чем обычные запросы. Также достоинством Tornado является поддержка веб-сокетов.

3.2.3 База данных

Наша программа не зависит от конкретного типа СУБД, так как в Django имеется удобный механизм объектно-реляционного отображения моделей в таблицы БД. По умолчанию используется СУБД SQLite, главным достоинством которой является то, что она представляет собой автономную базу данных, т.е. не требуется развертывать целый сервер и настраивать его. Однако для более высокой производительности и дополнительных возможностей, которые не предоставляет SQLite, такие, как работа с JSON, процедуры, эффективная фильтрация, рекомендуется использовать СУБД PostgreSQL. Для того, чтобы программа могла работать с другим видом СУБД, нужно установить соответствующий драйвер СУБД и указать в главном конфигурационном файле название базы данных, хост, имя пользователя и его пароль.

3.2.4 Работа с DICOM

Как уже было сказано ранее, для хранения, обработки и передачи медицинской информации используется стандарт DICOM (Digital Imaging and Communications in Medicine). Данный стандарт включает в себя описание файла хранения медицинских изображений, а также сетевой протокол взаимодействия. На данный момент наша программа успешно работает с DICOM файлами, содержащими 8-, 16- и 32-битные чёрно-белые несжатые изображения, а также реализует сетевой протокол DICOMWeb поверх HTTP/HTTPS.

Каждый DICOM-файл имеет две основные части: информационные тэги и изображение или серию изображений. DICOM-файл реализует информационную модель типа “один-ко-многим” и хранит информацию о четырех сущностях таких, как пациент, обследование, серия обследований и экземпляр изображения. Сущность “Пациент” может как хранить данные о пациенте такие, как пол, ФИО, возраст, год рождения, а может быть полностью анонимизирована. Сущность “Обследование” включает в себя сведения о дате обследования, типе обследования, враче, который проводил обследование, и назначение обследование. Сущность “Серия обследований” содержит сведения о типе оборудования, на котором проводили обследование, дата проведения, положении пациента и т.д. Сущность “Экземпляр изображения” содержит данные о формате изображения, ее размеры, разрешения и цветовая гамма. По стандарту DICOM-файл может содержать как одно изображение, так и несколько. Изображение может быть как черно-белым, так и цветным. Также предусматривается возможность хранить сжатые данные в формате JPEG.

Для работы с DICOM мы используем библиотеку PyDICOM, которая позволяет открывать файлы, просматривать тэги и выгружать исходное изображение. Также библиотека позволяет редактировать файл и сохранять его обратно.

3.2.5 WebGL и ThreeJS

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

3.4 Функциональные возможности программы

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

3.4.1 Загрузка МРТ снимков

Программа предусматривает загрузку снимков в формате DICOM. На данной этапе не имеет значение, какой тип снимка будет загружен в базу данных. Важно, чтобы снимки имели корректный формат в соответствии со стандартом. Для загрузки снимка пользователю нужно в главном меню найти кнопку “Upload DICOM” и нажать на нее (рис. 17). В окне откроется форма для загрузки снимков (рис. 18). В данной форме предусматривается возможность загрузки снимка путем перетаскивания или открытия через файловый диалог. Чтобы перетащить файлы в форму, пользователь должен выбрать файлы на компьютере и медленно переместить иконки файлов в форму для перетаскивания. После перетаскивания под формой появится список загруженных файлов с их названиями. Другой способ загрузки заключается в том, чтобы нажать на область для перетаскивания, после чего откроется файловый диалог. В открывшемся диалоге нужно найти папку с файлами и выбрать файлы в ней. В результате выполнения этого действия будет сформирован список загруженных файлов. Последний шаг заключается в том, чтобы подтвердить отправку файлов на сервер. Для этого нужно нажать на кнопку “Upload”. Если процесс загрузки завершится удачно, то будет показано сообщение об успешном выполнении (рис. 19).

Рисунок 17. Загрузка DICOM

Рисунок 18. Загрузка DICOM

Рисунок 19. Загрузка DICOM



3.4.2 Поиск МРТ снимков

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

Рисунок 20. Просмотр пациентов



Рисунок 21. Просмотр обследований

Рисунок 22. Просмотр серий обследований

Рисунок 23. Просмотр изображений



3.4.3 Визуализация МРТ снимков

В окне визуализации (рис. 24) имеется панель инструментов и область со снимком. На верхней панели пользователь настраивать цветовую гамму и фильтры, прокручивать вперед и назад снимки, а также выбирать плагины для работы со снимком. Для прокручивания назад снимка нужно нажать на кнопку (рис. 25), а для прокручивания вперед кнопку (рис. 26). Врач может изменить цветовую гамму в выпадающем меню (рис. 27), а также выбрать интересующий фильтр. Для применения плагина нужно в выпадающем меню (рис. 28) выбрать плагин и нажать на строку с ним. После этого откроется окно для работы с плагином.

Рисунок 24. Просмотр изображений

Рисунок 25. Просмотр изображений

Рисунок 26. Просмотр изображений



Рисунок 27. Просмотр изображений

Рисунок 28. Просмотр изображений

3.4.4 Сегментация раковых опухолей

В открытом окне для работы с плагином (рис. 29) имеется область просмотра изображения и панель инструментов. Вначале пользователь должен настроить процесс сегментации. Для этого он должен нажать на кнопку (рис. 30) и в открытом диалоговом окне в форме задать нужные значения параметров (рис. 31). После завершения заполнения нужно нажать на кнопку “Apply”. Пока плагин обрабатывает снимок, весь экран темнеет и в середине показывается статус процесса обработки (рис. 32). После завершения работы появляется сообщение с указанием успешности выполнения и площадью вычисленной опухоли (рис. 33). После этого в области просмотра снимка на снимок накладывается результат сегментации (рис. 34). Врач может задать различные типы масок. В программе предусмотрено четыре вида масок: смешивание, вырезание, смешение с затемнение исходного изображения и полное смешение.



Рисунок 29. Сегментация изображения

Рисунок 30. Сегментация изображения

Рисунок 31. Сегментация изображения

Рисунок 32. Сегментация изображения



Рисунок 33. Сегментация изображения

Рисунок 34. Сегментация изображения

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

Для того, чтобы получить доступ к плагинам, нужно в главном окне найти раздел “Plugins” и нажать на него (рис. 35). В открытом окне будет список всех плагинов, доступных на данный момент (рис. 36). Часть плагинов уже установлены в системе, а часть еще нет. Пользователь может удалить плагин или установить новый. Для каждого плагина приводится краткая информация о авторе, версии, типе плагина (сейчас предусматривается три вида плагина: фильтр, анализатор и определитель сегментов) и т.д.

Рисунок 35. Плагины



Рисунок 36. Плагины

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

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

- Загрузка и просмотр DICOM файлов;

- Настройка визуализации;

- Сегментация опухоли с помощью выбранного метода;

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

- Фильтрация и поиск снимков.



Заключение

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

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

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

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

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

2. Чтение информации из DICOM о снимке и выгрузка изображения;

3. Удаленный доступ к веб-сервисам DICOMWeb и DICOM;

4. Система плагинов на языке Python с возможностью написания дополнительных расширений на языках C/C++ и CUDA C++;

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

6. 2D-визуализация МРТ-снимков с возможностью применения фильтров, написанных на шейдерах;

7. Обработка МРТ-снимков и сегментация опухоли.

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

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



Список использованной литературы

1. Cancer Research UK. "Cancer statistics worldwid". [Электронный ресурс] //URL: http://www.cancerresearchuk.org/health-professional/cancer-statistics/worldwide-cancer#heading-Three (дата обращения: 04.02.2018).

2. Clang LLVM. "Clang" https://clang.llvm.org/get_started.html (дата обращения: 04.02.2018);

3. CMake. "CMake". https://cmake.org/ (дата обращения: 04.02.2018);

4. Cython. "Cython". http://cython.org/ (дата обращения: 04.02.2018);

5. D. Comaniciu, P. Meer. Mean shift: a robust approach toward feature space analysis // IEEE Transactions on Pattern Analysis and Machine Learning. - 2003. - с. 604-606;

6. Despotovic, B. Goosens, W. Philips. MRI Segmentation of the Human Brain: Challenges, Methods, and Applications. // Computational and Mathematical Methods in Medicine. - 2014. - с. 2 - 23.

7. Django. "Django". https://www.djangoproject.com/ (дата обращения: 04.02.2018);

8. Edoras. "Make - UNIX Utility". https://edoras.sdsu.edu/doc/make.html (дата обращения: 04.02.2018);

9. Jetbrains. "PyCharm". https://www.jetbrains.com/pycharm/ (дата обращения: 04.02.2018);

10. Jetbrains. "Webstorm". https://www.jetbrains.com/webstorm/ (дата обращения: 04.02.2018);

11. K. Karimi, N.G. Dickson, F. Hamze. A Performance Comparison of CUDA and OpenCL. // Cornell University Library. - 2010. - с. 1-10;

12. N. Otsu. A threshold selection from gray-level histograms. // IEEE Transactions on Systems, Man, and Cybernetics. - 1979. - с. 62 - 66.

13. Numpy. "Numpy". http://www.numpy.org/ (дата обращения: 04.02.2018);

14. NVIDIA Developer Blog. "An Even Easier Introduction to CUDA". [Электронный ресурс] //URL: https://devblogs.nvidia.com/even-easier-introduction-cuda/ (дата обращения: 04.02.2018).

15. NVIDIA Developer Blog. "How to Implement Performance Metrics in CUDA C/C++". [Электронный ресурс] //URL: https://devblogs.nvidia.com/how-implement-performance-metrics-cuda-cc/ (дата обращения: 04.02.2018).

16. NVIDIA Developer Blog. "Using Shared Memory in CUDA C/C++". [Электронный ресурс] //URL: https://devblogs.nvidia.com/using-shared-memory-cuda-cc/ (дата обращения: 04.02.2018).

17. Nvidia. "CUDA". https://docs.nvidia.com/cuda/ (дата обращения: 04.02.2018);

18. Nvidia. "Nsight". http://www.nvidia.com/object/nsight.html (дата обращения: 04.02.2018);

19. Nvidia. "NVCC". https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html (дата обращения: 04.02.2018);

20. OpenCV. "OpenCV". https://opencv.org/ (дата обращения: 04.02.2018);

21. Orthanc. Orthanc PACS server. [Электронный ресурс] //URL: https://www.orthanc-server.com/ (дата обращения: 04.02.2018);

22. OsiriX. OsiriX DICOM viewer. [Электронный ресурс] //URL: https://www.osirix-viewer.com/ (дата обращения: 04.02.2018);

23. PostgreSQL. "PostgreSQL". https://www.postgresql.org/ (дата обращения: 04.02.2018).

24. Pydicom. "Pydicom". https://pydicom.github.io/ (дата обращения: 04.02.2018);

25. Python. "Python". [Электронный ресурс] //URL: https://www.python.org/ (дата обращения: 04.02.2018);

26. R. Baert. The Gale encyclopedia of cancer: a guide to cancer and its treatment. - Фармингтон Хилс: Gale Group, 2002. - с. 93;

27. R. Farnoosh, B. Zarpak. Image segmentation using gaussian mixture model. // IUST International Journal of Engineering Science. - 2008. - с. 29-32;

28. S. Perreault, P. Hebert. Median filter in constant time // IEEE. - 2007. - с. 1-2;

29. Scikit-Learn. "Scikit-Learn". http://scikit-learn.org/stable/ (дата обращения: 04.02.2018);

30. T. Kanungo, N.S. Netanyahu, A.Y. Wu. An efficient k-Means clustering algorithm: analysis and implementation. // IEEE Transactions on Pattern Analysis and Machine Intelligence. - 2002. - с. 861 - 892.

31. TCIA. "Low grade glioma". [Электронный ресурс] //URL: https://wiki.cancerimagingarchive.net/display/Public/LGG-1p19qDeletion (дата обращения: 04.02.2018);

32. Tician. "PyCUDA". https://documen.tician.de/pycuda/ (дата обращения: 04.02.2018);

33. Tornado. "Tornado". http://www.tornadoweb.org/en/stable/ (дата обращения: 04.02.2018);

34. World health organization. "Cancer key facts". [Электронный ресурс] //URL: http://www.who.int/mediacentre/factsheets/fs297/en/ (дата обращения: 04.02.2018);

35. Z. Othman, H. Haron, M.R.A. Kadir. Comparison of Canny and Sobel Edge detection in MRI Images // Faculty of Computing. - 2013. - с. 1-2;

36. Д. Сандерс, Э. Кэндрот. Технология CUDA в примерах: введение в программирование графических процессоров. - М.: ДМК Пресс, 2011. - с. 227;

37. М.П. Димашова. Реализация алгоритма сегментации изображений Mean Shift на GPU // Нижегородский государственный университет им. Н.И. Лобачевского. - 2010. - с. 214-217.

Размещено на Allbest.ru