Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,940

ANALYSIS OF THE COMPUTATIONAL COMPLEXITY OF ALGORITHMS FOR SINGLING OUT AREAS OF INTEREST ON MRI DATA

Shustova M.V. 1
1 Ailamazyan Program Systems Institute
Today, cell therapy is one of the topical trends in regenerative medicine, and the use of magnetic resonance imaging is one of the main ways to carry out fundamental biomedical research. Analysis of the data of magnetic resonance imaging makes it possible to study areas of interest and the dynamics of their change. The combination of modern methods of magnetic resonance imaging and methods of scientific imaging makes it possible to build maps of the migration and distribution of stem cells after their transplantation, which gives an impetus to understanding the mechanisms of their action on the foci of the disease, for example, in ischemic stroke. However, the large volumes of information received from the tomograph slow down the speed of their interpretation by doctors in clinical and scientific research. To increase the efficiency of informational support of the researcher, methods and software have been developed for the automatic singling and visualization of areas of interest (clusters of mesenchymal stem cells and areas of ischemic lesion) on magnetic resonance imaging images. This paper provides descriptions of the developed algorithms in the pseudocode language. The presented methods and software are intended to be included in the information support system for the work of medical researchers. Estimates of their computational complexity are given. An experimental analysis of the performance of the developed algorithms is carried out.
image processing
computational complexity
magnetic resonance imaging
mesenchymal stem cells
ischemic stroke
algorithms
time complexity estimates

Одним из современных направлений исследований в области биомедицины является изучение влияния трансплантированных мезенхимальных стволовых клеток (МСК) на очаги ишемического инсульта, поскольку трансплантированные МСК являются эффективными терапевтическими агентами [1–3]. На данный момент остается до конца не исследованным механизм поведения и взаимодействия стволовых клеток с локальным микроокружением. Автоматизация анализа данных магнитно-резонансной томографии (МРТ) за счет методов научной визуализации предоставит возможность отслеживания МСК после трансплантации и построения карт их миграции, что является крайне важным для понимания механизмов воздействия МСК на очаги ишемического инсульта.

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

Материалы и методы исследования

Анализ алгоритма выделения мезенхимальных стволовых клеток

Разработанный алгоритм выделения скоплений МСК построен на одновременной обработке пары снимков МРТ (до введения МСК и после введения) и состоит из нескольких шагов [5]:

1) волновой алгоритм;

2) фильтрация по яркости;

3) удаление неинформативных объектов;

4) обработка T раз с разными наборами параметров.

Определим вычислительную сложность каждого шага алгоритма. Сложность волнового алгоритма известна и составляет O(MN) для матрицы размером M×N.

Фильтрация изображения по яркости размером M×N происходит по следующему правилу: рассматривается каждый пиксель изображения, и если яркость пикселя превышает заданный порог, то он удаляется. Для этого шага требуется O(MN) операций.

Удаление неинформативных объектов состоит из последовательных операций над двумя изображениями до и после введения стволовых клеток.

1. Выполняется маркировка оставшихся после предыдущих шагов объектов (рис. 1). Для обработки одного изображения размером M×N требуется O(MN) операций.

hustov1.tif

Рис. 1. Алгоритм маркировки объектов

2. После маркировки выполняется подсчет площади (количества пикселей) у каждого объекта. Количество требуемых операций – O(MN).

hustov2.tif

Рис. 2. Алгоритм фильтраций

3. Далее последовательно выполняются несколько фильтраций (рис. 2):

- фильтрация по количеству пикселей у объектов;

- удаление неравномерно пересекающихся объектов;

- фильтрация окном круглой формы.

Этот шаг выполняется не больше чем mCounterAfter раз, а поскольку mCounterAfter < M×N, то эта часть шага ограничена сверху O(MN) операциями.

Для улучшения показателей точности и полноты МРТ-снимок обрабатывается T раз с разными наборами параметров – соответственно, получается T масок. После этого выполняется дополнительная фильтрация: полученные изображения накладываются друг на друга и попиксельно считываются. Если рассматриваемый пиксель равен 255 как минимум на 80 % масок, то данный пиксель включается в результирующую маску. Требуемое количество операций: O(MN).

Анализ алгоритма выделения зоны ишемического поражения

Разработанный алгоритм выделения зоны ишемического поражения состоит из двух шагов [6]:

1) вычисление характерных признаков;

2) анализ полученных признаков с применением классификатора.

Для обработки изображений требуется определить два класса – в нашем случае это область здорового мозга («норма», класс Y1) и область ишемического поражения («поражение», класс Y2). Каждая область представляется набором небольших эталонных изображений одного из срезов мозга. Из этих изображений извлекаются признаки, необходимые для настройки классификатора.

Определим вычислительную сложность каждого шага алгоритма.

На первом шаге вычисляются характерные признаки. Всего задействовано три признака общего назначения и четыре текстурных признака Харалика. Признаки общего назначения вычисляются по следующим формулам [7]:

- математическое ожидание μ:

hustov01.wmf,

- стандартное отклонение:

hustov02.wmf,

- асимметричность уровня серого:

hustov03.wmf,

где M и N – размеры сканирующего окна, pij – значение яркости пикселя в i-й строке и j-м столбце сканирующего окна. Признаки общего назначения извлекаются отдельно для каждого канала RGB – таким образом, всего извлекается девять общих признаков.

Основные текстурные признаки Харалика вычисляются следующим образом [8]:

- энтропия: hustov04.wmf

- энергия: hustov05.wmf

- контраст: hustov06.wmf

- гомогенность: hustov07.wmf

где Co(i, j) – элемент матрицы смежности уровня серого (Gray-Level Co-occurrence Matrix) в i-й строке j-го столбца.

Алгоритм извлечения характерных признаков представлен на рис. 3, а. Этот этап алгоритма выделения ишемического поражения построен на обработке изображения размером M×N с помощью сканирующего окна, следовательно, он ограничен сверху O(MN) операциями.

hustov3a.tif

а) извлечение характерных признаков

hustov3b.tif

б) анализ извлеченных признаков

Рис. 3. Алгоритм извлечения признаков и их анализа с помощью классификатора на основе расстояния Евклида – Махаланобиса

На втором шаге для работы с перечисленными признаками применяется классификатор, использующий расстояние Евклида – Махаланобиса [9]. Расстояние между классом Y (например, классом «ишемическое повреждение») и точкой p (точка – вектор признаков, извлеченных из некоторой позиции сканирующего окна) вычисляется по формуле hustov08.wmf, где матрица hustov09.wmf. Ковариационная матрица CY для класса Y вычисляется по формуле

hustov10.wmf,

где hustov11.wmf – i-й признак k-й точки из класса Y, а hustov12.wmf – математическое ожидание значения этого признака в классе Y.

На рис. 3, б, приведен алгоритм анализа извлеченных признаков с помощью классификатора.

Этот этап алгоритма выделения ишемического поражения выполняется в том же цикле, что и предыдущий этап, следовательно, он также ограничен сверху O(MN) операциями.

Результаты исследования и их обсуждение

Проведем экспериментальный анализ полученной вычислительной сложности алгоритмов выделения МСК и зон ишемического поражения путем измерения времени выполнения их программных реализаций.

Анализ алгоритма выделения МСК

Алгоритм выделения МСК построен на одновременной обработке пары изображений, поэтому в качестве входных дан- ных использовались два снимка МРТ в разных разрешениях: от 64x44 пикселей до 512x352 пикселей. Для получения среднего времени обработки снимки обрабатывались 1000 раз. Программная реализация алгоритма запускалась на компьютере со следующей конфигурацией: Intel Core i3-8100 3.60 GHz, 16 GB RAM. В табл. 1 и рис. 4 приведены данные измерений времени работы алгоритма и аппроксимированные значения оценки.

hustov4.wmf

Рис. 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 видно, что подбором уточняющих коэффициентов удается достаточно точно приблизить теоретические оценки к экспериментальным результатам.

hustov5.wmf

Рис. 5. Графики роста (скорость работы алгоритма выделения зон ишемического поражения)

Заключение

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

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 17-29-07002-офи_м.