CC BY
29
УДК 61:001.89+61:347.77:57.08+57.08
К.С. Бразовский, канд. мед. наук, В. П. Демкин, д-р физ.-мат. наук,
Сибирский государственный медицинский университет, г. Томск Я. С. Пеккер, канд. техн. наук,
Национальный исследовательский томский государственный университет
Трехмерная реконструкция поверхности биологических объектов с использованием высокопроизводительного вычислительного кластера
Ключевые слова: трехмерная реконструкция, высокопроизводительные вычисления, медицнская визуализация. Key words: 3D raconstruction, high performance computiong, medical imaging.
Трехмерная реконструкция сложных поверхно- использовании профессиональных видеопроцес-
стей биологических объектов играет важную роль соров визуализация такой модели особой слож-
в медицинской визуализирующей диагностике, ности не представляет, однако при решении задач
так как позволяет не только корректно отображать математического моделирования и создания фи-
органы и ткани, но и создавать аппроксимацию зического аналога с помощью стереолитографии
границ анатомических структур с помощью про- избыточная детализация объекта создает боль-
В статье рассмотрены вопросы параллельной реализации алгоритма реконструкции поверхности биологических объектов. Предложен способ разделения изображения для метода деформируемой поверхности. Алгоритм с блочным распределением трехмерного изображения по вычислительным узлам был реализован на вычислительном кластере СКИФ-Сибирь. Вычислительный эксперимент был проведен на 16 вычислительных узлах и показал ускорение более 8. Работа выполнена в рамках федеральной целевой программы «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007—2013 годы», научно-исследовательских работ по лоту шифр 2011-1.4-514-111 «Проведение проблемно-ориентированных поисковых исследований в области информационно-телекоммуникационных систем для решения задач технологической платформы „Биоиндустрия и биоресурсы — БиоТех2030"», государственный контракт № 07.514.11.4054.
стых геометрических объектов. Геометрическая модель может быть использована в дальнейшем для математического моделирования механических, электрических и других свойств биологического объекта, а также для создания реальных физических моделей с помощью метода стереолитографии.
Введение
В основе реконструкции поверхности трехмерного объекта лежит принцип определения границ раздела отдельных органов и тканей. Выделенные границы служат в качестве ориентиров для формального описания поверхности. Существует несколько широко распространенных методов выделения поверхностей, среди которых особенно выделяется метод marching cubes, алгоритм «шагающих кубов» [1]. Достоинством метода является высокая скорость работы и естественный параллелизм вычислительной процедуры, что делает его пригодным для эффективной реализации на специализированной аппаратуре, в частности на графических процессорах (GPU, graphic processor unit). Несмотря на то что метод изначально разрабатывался для обработки медицинских изображений, он имеет ряд недостатков именно с точки зрения трехмерной реконструкции поверхности биообъектов. При обработке трехмерных медицинских изображений высокого разрешения поверхность получается слишком сложной, содержащей несколько десятков миллионов треугольников. При
Рис. 1 Реконструированная поверхность головы с помощью метода «шагающих кубов». Количество поверхностных треугольных конечных элементов более 2 млн
шие проблемы, связанные с размерностью задачи. В частности, модель поверхности головы человека (рис. 1), построенная на основе магнитно-резонансного трехмерного изображения с разрешением 1 мм по каждой оси, содержит более 2 млн треугольников. Реконструкция выполнена методом «шагающих кубов», реализованных в программе Slicer 3D [2].
Визуальное качество реконструкции достаточно для качественного анализа изображений, однако при создании трехмерной математической модели на основе данной поверхности приблизительная (предвычисленная) размерность модели составляет более 200 млн трехмерных конечных элементов в форме тетраэдра. Кроме того, большинство алгоритмов объемной реконструкции предполагают, в явном или неявном виде, ограниченную кривизну поверхности и ее неразрывность. При создании модели методом «шагающих кубов» в общем случае эти ограничения выполнить невозможно, и приходится применять дополнительные алгоритмы восстановления и улучшения поверхности (алгоритмы «латания дыр», filling holes) [3].
Метод деформируемых поверхностей
Альтернативным способом реконструкции является метод деформируемых поверхностей [4], который показал хорошие результаты при реконструкции медицинских изображений, в том числе с очень сложной поверхностью [5].
Метод деформируемых поверхностей появился как обобщение двумерного метода деформируемого контура для трехмерного пространства, однако увеличение размерности привело к появлению ряда специфических проблем метода:
1) исходная и реконструируемая поверхности должны быть топологически эквивалентны (обе должны быть замкнутыми или незамкнутыми, одно-связными или нет); тип деформируемой поверхности выбирают с учетом анатомических особенностей реконструируемой поверхности;
2) требуется определение ряда формальных характеристик деформируемой поверхности («упругость», «эластичность»), которые характеризуют гладкость реконструированной поверхности и скорость сходимости решения;
3) важнейшую роль для качественной реконструкции методом деформируемых поверхностей играет выбор оператора воздействия, аналога механической силы, сжимающего или растягивающего поверхность.
Для метода «шагающих кубов» не требуется таких исходных данных для реконструкции, что и обуславливает его большую универсальность и более широкое применение. Однако в случае обработки медицинских изображений мы можем пожертвовать универсальностью для достижения гибкости в настройке параметров реконструированной поверхности (степень детализации, гладкость, предельная кривизна).
Поверхность в трехмерном пространстве может быть представлена в виде функции трех переменных:
S(u, v) = x(u, v), y(u, v), z(u, v),
(1)
где 5 — поверхность; х, у, х — декартовы координаты поверхности; и, V — в общем случае криволинейные координаты на поверхности.
Необходимость взаимного преобразования координат обусловлена тем, что трехмерные медицинские изображения представлены в декартовых координатах, а вычисления удобнее проводить в координатах реконструируемой поверхности.
Идея метода деформируемых поверхностей основана на аппроксимации существующей поверхности биологического объекта 5о^, для которой не существует приемлемого формального описания, поверхностью с легко вычисляемыми параметрами:
м
5(и, V) = 5о^(х, у, х); 5(и, V) = ^ 5'(иь, vi), (2)
I=1
где 5(и, V) — аппроксимирующая поверхность; 8'(и^ Vi) — элемент поверхности; М — количество элементов в аппроксимирующей поверхности.
Нахождение поверхности биообъекта на исходном изображении является нетривиальной задачей и может быть выполнено различными способами. В частности, она может быть получена методом «шагающих кубов», но практически такой подход применяется редко вследствие больших вычислительных затрат. Чаще всего граница на трехмерном изображении определяется как набор точек с определенным градиентом яркости. Существует большое количество хорошо исследованных операторов, выделяющих края, в том числе и на трехмерном изображении, в частности, лапласиан гауссова фильтра. С точки зрения программной
реализации оператор выделения границы — это цифровой фильтр следующего вида:
п1 п 2 п3 F(х, у, х) = ^ ^ ^ а^/(х + i, у + ], х + к),
i=-n1]=-п2к=—п3 (3)
где F(x, у, х) — преобразованное изображение с выделенной границей; а^ — коэффициенты фильтра; 1(х, у, х) — исходное изображение; п1, п2, пз — порядок фильтра вдоль соответствующей оси. Как правило, для удобства выбирают п1 = п2 = п3, но в общем случае порядок фильтра может быть разным по разным осям. Предполагается, что координаты изображения х, у, х дискретны.
Процесс аппроксимации может быть описан уравнением движением точек поверхности под воздействием двух сил: внешней деформирующей силы и внутренней силой упругости границы биологического объекта:
М + х аЩ + ЬХ(г) = Fext (X, г), (4)
йг2
йг
где М — масса поверхности; X — коэффициент сглаживания (демпфирования); Ь — коэффициент упругости поверхности, Fext — сила сопротивления границы биологического объекта.
При выполнении условия (2) инерция поверхности равна нулю, что соответствует стационарному состоянию, поэтому в расчетах мы можем принять М = 0. Кроме того, при реконструкции статических изображений сила сопротивления Fext будет зависеть только от координат. В результате процесс реконструкции поверхности сводится к численному решению неоднородного нелинейного относительно координат дифференциального уравнения (4).
Очень важную роль играет выбор вычислительной схемы решения уравнения (4). В нашем случае не требуется высокой точности, поэтому подход на основе линеаризации вполне себя оправдывает. Заменив в выражении (4) первую производную простейшим разностным аналогом, получим систему линейных алгебраических уравнений для вычисления малого приращения координат:
Х (г0) = Х0;
X(^) = X(^—1) + aЬX—1) + pFext (X), i = 1, ..., ы,
(5)
где i — номер итерации, ^ - 1 — временной шаг; N — максимально разрешенное количество итераций; а < 1, Р < 1 — коэффициенты влияния сил: коэффициент Р определяет скорость деформации (чем больше, тем быстрее деформируется исходная поверхность), коэффициент а определяет гладкость поверхности; Fext — сила сопротивления поверхности, вычисляемая с помощью (4). Итерации (5) прекращают при достижении максимально разрешенного количества, либо когда приращение
координат становится меньше заданной малой величины: е << 1.
Для уменьшения размерности системы (5) исходную деформируемую поверхность разбивают на простые конечные элементы с аналитически вычисляемой функцией формы, чаще всего на плоские треугольники. В этом случае матрица системы (5) является разреженной и максимальное количество ненулевых элементов в каждой строке X] равно числу соседних с X] узлов.
Для решения системы (5) на каждом шаге можно использовать как прямые, так и итерационные методы. Прямые методы эффективны только при сравнительно небольшом количестве узлов деформируемой поверхности. Поскольку в реальных изображениях поверхность состоит из нескольких десятков тысяч конечных элементов, то более походящими являются итерационные методы решения (5). В своей работе мы использовали метод сопряженных градиентов. На рис. 2 приведена реконструированная поверхность головы человека. Вычисления проводились на основе того же изображения, которое использовалось для получения поверхности на рис. 1, однако реконструированная поверхность содержит всего 10 240 поверхностных конечных элементов, что более чем на два порядка меньше по сравнению с реконструкцией на основе алгоритма «шагающих кубов» при сравнимой погрешности аппроксимации.
Регулируя параметр а можно изменять степень сглаживания поверхности, а уменьшая параметр Р можно увеличить точность аппроксимации ценой увеличения времени реконструкции.
Существенным недостатком метода деформируемых поверхностей является кубическая зависимость между количеством узлов в деформируемой поверхности и количеством арифметических операций. При большом количестве узлов, когда требуется точная аппроксимации сложной поверхности, время счета может стать критическим фактором.
Кроме того, по сравнению с методом «шагающих кубов» использование деформируемых поверхно-
Рис. 2
Реконструированная поверхность головы с помощью метода деформируемой поверхности. Количество поверхностных треугольных конечных элементов 10240
стеи для реконструкции трехмерных медицинских изображений требует значительных вычислительных ресурсов и существенно сложнее в части программной реализации и отладки кода.
Параллельная реализация этого алгоритма также затруднена в связи с тем, что на каждом шаге итерации необходимо вычислять новую матрицу системы линейных уравнений (5) на основе исходного изображения и деформируемой поверхности, что возможно при хранении исходного изображения на каждом узле параллельной вычислительной системы. Однако это приведет к неоправданному расходу оперативной памяти либо существенно замедлит вычисления при хранении изображения на внешнем носителе. В этой связи весьма актуальна задача разработки эффективного параллельного алгоритма реконструкции методом деформируемых поверхностей. Наиболее очевидный вариант — выполнять параллельное решение системы линейных уравнений (5) с помощью библиотечных процедур. Узким местом такого алгоритма будет необходимость пересылки на каждом шаге итерации большого объема данных вычислительным узлам. Проведенные нами исследования на вычислительном кластере СКИФ-Сибирь Национального исследовательского томского государственного университета показали, что существенного прироста производительности по сравнению с последовательной реализацией не наблюдается. Основные затраты времени связаны с формированием матрицы системы уравнений (5) (последовательная процедура, выполняемая управляющим процессом), ее рассылкой по вычислительным узлам, ожиданием решения и вычислением нового значения вектора координат (также последовательная процедура, выполняемая управляющим процессом).
Для уменьшения расчетного времени и эффективного использования ресурсов высокопроизводительного вычислительного кластера мы предлагаем подход к распараллеливанию задачи, который успешно оправдал себя при решении задачи медицинской визуализации на основе электроимпедансной томографии [6]: использование естественной блочной структуры трехмерного медицинского изображения. Был разработан следующий параллельный алгоритм деформируемой поверхности:
1) исходное изображение разделяется на трехмерные блоки (рис. 3); количество блоков определяется количеством доступных вычислительных узлов: по одному блоку на узел; для каждого блока создается список соседних блоков (от 2 до 6) с указанием номера вычислительного узла, на котором находится соседний блок;
2) формируется исходная недеформированная поверхность, выполняются ее триангуляция и разделение на фрагменты соответственно выделенным блокам исходного изображения; в каждый фрагмент входят только те конечные элементы, координаты вершин которых лежат в пределах блока изображения;
Рис. 3
Распределение блоков трехмерного изображения между вычислительными узлами
3) на каждом вычислительном узле проводятся решение системы (5) и вычисление нового значения вектора координат вершин;
4) в том случае, когда в результате коррекции координат одна вершина выходит за границы блока, процесс запрашивает данные для текущей итерации у соседнего блока; в том случае, когда две или три вершины покидают границы блока, весь треугольник целиком передается соседнему блоку.
Данный алгоритм позволяет уменьшить интенсивность обмена между узлами вычислительного кластера и эффективно использовать присущий алгоритму деформируемой поверхности естественный параллелизм.
Реализация алгоритма была выполнена на языке программирования C++, использовались оптимизирующий компилятор Intel, библиотеки MKL и MPI. При тестировании на упомянутом выше магнитно-резонансном изображении было показано, что использование 16 вычислительных узлов дает максимальное ускорение по сравнению с последовательной реализацией в пределах от 8 до 9, в среднем 8,7.
Заключение
Разработан параллельный алгоритм для трехмерной реконструкции поверхности биологических объектов, который обеспечивает эффективное использование вычислительных ресурсов высокопроизводительного кластера. Реализация алгоритма выполнена с максимальной независимостью от аппаратуры, при этом удалось достичь существенного прироста производительности за счет оптимизации количества обменов между вычислительными узлами.
| Литература |
1. William E. Lorensen, Harvey E. Cline: Marching Cubes: A high resolution 3D surface construction algorithm // In: Computer Graphics. 1984. Vol. 21, N 4.
2. The NA-MIC Kit: ITK, VTK, Pipelines, Grids and 3D Slicer as an Open Platform for the Medical Image Computing Community / S. Pieper, B. Lorensen, W. Schroeder, R. Kiki-nis // Proceedings of the 3rd IEEE International Symposium on Biomedical Imaging: From Nano to Macro. 2006. Vol. 1. Р. 698-701.
3. Filling Holes in Triangular Meshes by Curve Unfolding / A. Brunton, S. Wuhrer, C. Shu [et al.] // Shape Modeling and
6.
Applications.IEEE International Conference. 2009. Р. 66-72. [Электронный ресурс, дата обращения 15.03.2012] http:// hdl.handle.net/1721.1/53749.
Staib H. L. Model-based Deformable Surface Finding for Medical Images //IEEE transaction on medical imaging. 1996. Vol. 5. N 5. Р. 720-731.
Law M., Chung A. A Deformable Surface Model for Vascular Segmentation // Medical Image Computing and ComputerAssisted Intervention (MICCAI). 12th Interntional Conference. London, 2009. Р. 59-67.
Пеккер Я. С., Бразовский К. С., Фокин А. В. Применение высокопроизводительных параллельных вычислений для решения задач электроимпедансной томографии // Биотехносфера. 2010. № 5-6. С. 25-29.
УДК 612.172.4+303.717 Е. С. Якушенко, аспирант
Санкт-Петербургский государственный электротехнический университет «ЛЭТИ»
Программа моделирования ЭКГ в среде LabVIEW
Ключевые слова: имитатор ЭКГ, вариабельность сердечного ритма, LabVIEW. Key words: ECG imitator, heart rate variability, LabVIEW.
Представлен программно-аппаратный комплекс, формирующий искусственную электрокардиограмму. В основе программы лежат алгоритмы формирования морфологии QPR5T-комплекса и вариабельности сердечного ритма. Комплекс способен отображать сигнал ЭКГ на дисплее компьютера, выводить его через ЦАП и записывать в файл. Есть возможность работы с записями ЭКГ в формате М1Т-В1Н.
Введение
В настоящее время существует множество различных методик, аппаратных и программных средств для регистрации и обработки сигналов электрокардиограммы (ЭКГ). Извлечение полезной клинической информации из реальной ЭКГ требует надежных технологий обработки сигнала. Они включают обнаружение К-зубца, измерение частоты сердечных сокращений и частоты дыхания по ЭКГ, определение параметров вариабельности сердечного ритма (ВСР) —УЪГ, ЪГ, ИГ и ЬЕ/ИЕ и т. д. [1]
Существует необходимость сравнения различных методов обработки кардиосигналов между
собой для выбора оптимального и наиболее подходящего алгоритма для решения поставленных задач. Также бывает трудно понять, как будут меняться эффективность и надежность этих методов в различных клинических условиях при широком диапазоне уровня шума и частот дискретизации [1].
Такая оценка возможна при наличии реальных сигналов ЭКГ [1], но не всегда имеется возможность снимать сигнал ЭКГ с пациента. Также реальные сигналы могут не обладать необходимыми характеристиками. Данная задача может быть решена путем создания программного комплекса, который способен был бы генерировать реалистичный искусственный сигнал ЭКГ. Разработка такой системы и является целью данной работы.
Методы решения
В основе программы лежит модель формирования морфологии PQRST-комплекса, описанная в работах [1, 2]. Модель основана на системе из трех простых дифференциальных уравнений (1). Она генерирует траекторию в трехмерном пространстве с координатами (х, у, х). В течение каждого сокращения сердца ЭКГ имеет квазипериодичную волновую форму, представленную зубцами Р, Q,
