Одним из современных направлений исследований в области биомедицины является изучение влияния трансплантированных мезенхимальных стволовых клеток (МСК) на очаги ишемического инсульта, поскольку трансплантированные МСК являются эффективными терапевтическими агентами [1–3]. На данный момент остается до конца не исследованным механизм поведения и взаимодействия стволовых клеток с локальным микроокружением. Автоматизация анализа данных магнитно-резонансной томографии (МРТ) за счет методов научной визуализации предоставит возможность отслеживания МСК после трансплантации и построения карт их миграции, что является крайне важным для понимания механизмов воздействия МСК на очаги ишемического инсульта.
За последние несколько лет в мире было разработано большое количество профессиональных систем и различных программных средств обработки и визуализации данных МРТ. Практически все медицинские программы имеют открытые исходные коды для анализа и обработки изображений и распространяются без лицензионных ограничений. Однако при ближайшем рассмотрении оказывается, что имеющихся программных средств недостаточно для решения конкретных специализированных задач, в частности прижизненной визуализации трансплантированных МСК в головном мозге реципиентов в условиях ишемического инсульта [4]. Для реализации автоматического выделения скоплений МСК и зон ишемического поражения были разработаны (с участием автора) соответствующие методы и алгоритмы. Настоящая работа посвящена подробному описанию этих алгоритмов и оценке их вычислительной сложности.
Материалы и методы исследования
Анализ алгоритма выделения мезенхимальных стволовых клеток
Разработанный алгоритм выделения скоплений МСК построен на одновременной обработке пары снимков МРТ (до введения МСК и после введения) и состоит из нескольких шагов [5]:
1) волновой алгоритм;
2) фильтрация по яркости;
3) удаление неинформативных объектов;
4) обработка T раз с разными наборами параметров.
Определим вычислительную сложность каждого шага алгоритма. Сложность волнового алгоритма известна и составляет O(MN) для матрицы размером M×N.
Фильтрация изображения по яркости размером M×N происходит по следующему правилу: рассматривается каждый пиксель изображения, и если яркость пикселя превышает заданный порог, то он удаляется. Для этого шага требуется O(MN) операций.
Удаление неинформативных объектов состоит из последовательных операций над двумя изображениями до и после введения стволовых клеток.
1. Выполняется маркировка оставшихся после предыдущих шагов объектов (рис. 1). Для обработки одного изображения размером M×N требуется O(MN) операций.
Рис. 1. Алгоритм маркировки объектов
2. После маркировки выполняется подсчет площади (количества пикселей) у каждого объекта. Количество требуемых операций – O(MN).
Рис. 2. Алгоритм фильтраций
3. Далее последовательно выполняются несколько фильтраций (рис. 2):
- фильтрация по количеству пикселей у объектов;
- удаление неравномерно пересекающихся объектов;
- фильтрация окном круглой формы.
Этот шаг выполняется не больше чем mCounterAfter раз, а поскольку mCounterAfter < M×N, то эта часть шага ограничена сверху O(MN) операциями.
Для улучшения показателей точности и полноты МРТ-снимок обрабатывается T раз с разными наборами параметров – соответственно, получается T масок. После этого выполняется дополнительная фильтрация: полученные изображения накладываются друг на друга и попиксельно считываются. Если рассматриваемый пиксель равен 255 как минимум на 80 % масок, то данный пиксель включается в результирующую маску. Требуемое количество операций: O(MN).
Анализ алгоритма выделения зоны ишемического поражения
Разработанный алгоритм выделения зоны ишемического поражения состоит из двух шагов [6]:
1) вычисление характерных признаков;
2) анализ полученных признаков с применением классификатора.
Для обработки изображений требуется определить два класса – в нашем случае это область здорового мозга («норма», класс Y1) и область ишемического поражения («поражение», класс Y2). Каждая область представляется набором небольших эталонных изображений одного из срезов мозга. Из этих изображений извлекаются признаки, необходимые для настройки классификатора.
Определим вычислительную сложность каждого шага алгоритма.
На первом шаге вычисляются характерные признаки. Всего задействовано три признака общего назначения и четыре текстурных признака Харалика. Признаки общего назначения вычисляются по следующим формулам [7]:
- математическое ожидание μ:
,
- стандартное отклонение:
,
- асимметричность уровня серого:
,
где M и N – размеры сканирующего окна, pij – значение яркости пикселя в i-й строке и j-м столбце сканирующего окна. Признаки общего назначения извлекаются отдельно для каждого канала RGB – таким образом, всего извлекается девять общих признаков.
Основные текстурные признаки Харалика вычисляются следующим образом [8]:
- энтропия:
- энергия:
- контраст:
- гомогенность:
где Co(i, j) – элемент матрицы смежности уровня серого (Gray-Level Co-occurrence Matrix) в i-й строке j-го столбца.
Алгоритм извлечения характерных признаков представлен на рис. 3, а. Этот этап алгоритма выделения ишемического поражения построен на обработке изображения размером M×N с помощью сканирующего окна, следовательно, он ограничен сверху O(MN) операциями.
а) извлечение характерных признаков
б) анализ извлеченных признаков
Рис. 3. Алгоритм извлечения признаков и их анализа с помощью классификатора на основе расстояния Евклида – Махаланобиса
На втором шаге для работы с перечисленными признаками применяется классификатор, использующий расстояние Евклида – Махаланобиса [9]. Расстояние между классом Y (например, классом «ишемическое повреждение») и точкой p (точка – вектор признаков, извлеченных из некоторой позиции сканирующего окна) вычисляется по формуле , где матрица . Ковариационная матрица CY для класса Y вычисляется по формуле
,
где – i-й признак k-й точки из класса Y, а – математическое ожидание значения этого признака в классе Y.
На рис. 3, б, приведен алгоритм анализа извлеченных признаков с помощью классификатора.
Этот этап алгоритма выделения ишемического поражения выполняется в том же цикле, что и предыдущий этап, следовательно, он также ограничен сверху O(MN) операциями.
Результаты исследования и их обсуждение
Проведем экспериментальный анализ полученной вычислительной сложности алгоритмов выделения МСК и зон ишемического поражения путем измерения времени выполнения их программных реализаций.
Анализ алгоритма выделения МСК
Алгоритм выделения МСК построен на одновременной обработке пары изображений, поэтому в качестве входных дан- ных использовались два снимка МРТ в разных разрешениях: от 64x44 пикселей до 512x352 пикселей. Для получения среднего времени обработки снимки обрабатывались 1000 раз. Программная реализация алгоритма запускалась на компьютере со следующей конфигурацией: Intel Core i3-8100 3.60 GHz, 16 GB RAM. В табл. 1 и рис. 4 приведены данные измерений времени работы алгоритма и аппроксимированные значения оценки.
Рис. 4. Графики роста (скорость работы алгоритма выделения МСК)
Таблица 1
Экспериментальные данные измерения времени работы алгоритма выделения МСК
Разрешение, px |
Площадь, px |
Масштаб, % |
Среднее время, мс |
Аппроксимированное значение оценки, мс |
64x44 |
2816 |
100 |
187,093 |
175,2634 |
96x66 |
6336 |
225 |
332,48 |
280,8634 |
128x88 |
11264 |
400 |
371,381 |
428,7034 |
160x110 |
17600 |
625 |
474,037 |
618,7834 |
192x132 |
25344 |
900 |
656,502 |
851,1034 |
224x154 |
34496 |
1225 |
861,543 |
1125,6634 |
256x176 |
45056 |
1600 |
1114,689 |
1442,4634 |
288x198 |
57024 |
2025 |
1493,989 |
1801,5034 |
320x220 |
70400 |
2500 |
1962,928 |
2202,7834 |
352x242 |
85184 |
3025 |
2441,073 |
2646,3034 |
384x264 |
101376 |
3600 |
2635,321 |
3132,0634 |
416x286 |
118976 |
4225 |
3178,921 |
3660,0634 |
448x308 |
137984 |
4900 |
3790,305 |
4230,3034 |
480x330 |
158400 |
5625 |
4536,728 |
4842,7834 |
512x352 |
180224 |
6400 |
5386,007 |
5497,5034 |
Анализ алгоритма выделения зон ишемического поражения
В качестве входных данных использовался снимок МРТ в разных разрешениях: от 64x45 пикселей до 512x360 пикселей. Для получения среднего времени обработки снимки обрабатывались 1000 раз. Программная реализация алгоритма запускалась на компьютере со следующей конфигурацией: Intel Core i3-8100 3.60 GHz, 16 GB RAM. В табл. 2 и рис. 5 приведены данные измерений времени работы алгоритма и аппроксимированные значения оценки.
Таблица 2
Экспериментальные данные измерения времени работы алгоритма выделения зон ишемического поражения
Разрешение, px |
Площадь, px |
Масштаб, % |
Среднее время, мс |
Аппроксимированное значение оценки, мс |
64x45 |
2880 |
100 |
37,706 |
43,2 |
96x68 |
6528 |
226,7 |
89,931 |
97,92 |
128x90 |
11520 |
400 |
161,107 |
172,8 |
160x113 |
18080 |
627,8 |
256,721 |
271,2 |
192x135 |
25920 |
900 |
373,501 |
388,8 |
224x158 |
35392 |
1228,9 |
513,973 |
530,88 |
256x180 |
46080 |
1600 |
674,432 |
691,2 |
288x203 |
58464 |
2030 |
856,569 |
876,96 |
320x225 |
72000 |
2500 |
1058,62 |
1080 |
352x248 |
87296 |
3031,1 |
1280,49 |
1309,44 |
384x270 |
103680 |
3600 |
1532,801 |
1555,2 |
416x293 |
121888 |
4232,2 |
1806,694 |
1828,32 |
448x315 |
141120 |
4900 |
2095,11 |
2116,8 |
480x338 |
162240 |
5633,3 |
2413,322 |
2433,6 |
512x360 |
184320 |
6400 |
2746,128 |
2764,8 |
Из рис. 4 и 5 видно, что подбором уточняющих коэффициентов удается достаточно точно приблизить теоретические оценки к экспериментальным результатам.
Рис. 5. Графики роста (скорость работы алгоритма выделения зон ишемического поражения)
Заключение
Проведенный анализ разработанных алгоритмов выделения зон интереса врача-исследователя показал, что полученные результаты в виде оценок их вычислительной сложности согласуются с результатами экспериментального анализа. Разработанные алгоритмы эффективно реализуются на типовом компьютере, что позволяет использовать их в составе исследовательского медицинского комплекса, предназначенного для врачей, занимающихся анализом МРТ-данных в лабораторных условиях.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 17-29-07002-офи_м.