Решается задача реконструкции внешних и внутренних границ органов, имеющих сложную пространственную структуру (мозг, сложные кости) по данным томографии. Для разрешения неоднозначностей реконструкция выполняется в ортогональных плоскостях. Рассматриваются 2 подхода. В первом, производится трехмерная сегментация методом разрастания регионов, и решается задача реконструкции двумерных контуров в двух ортогональных сечениях. Во втором - алгоритм реконструирует контуры в плоскостях, ортогональных осям координат, и объединяет их в граф поверхности. В обоих алгоритмах, полученный набор контуров триангулируется. За счет независимой обработки плоскостей, в параллельной реализации удалось добиться достаточно высокого ускорения. Ключевые слова: реконструкция, томография, сегментация, контур, неоднозначность.
На вход алгоритма поступают данные компьютерной томограммы, представляющие собой множество точек со значением поля. Точки размещаются в узлах трехмерной прямоугольной сетки. С помощью этих точек, используя интерполяцию, можно вычислить значение поля в произвольной точке. Совокупность точек с одинаковым значением поля образуют изоповерхность. Необходимо выполнить ее полигональную реконструкцию с информацией о смежных вершинах.
Идея первого подхода заключается в следующем. Сначала производится сегментация какоголибо органа. В данной работе применяется трехмерная сегментация с использованием простейшего метода разрастания регионов и метода активных контуров. В результате этого получаем набор вокселей, представляющий сегментированный объект. Затем, используя полученный набор вокселей, получаем контуры на каждом слое путем рассечения вокселей плоскостью, параллельной плоскости XY. Полученные контуры отдельно для каждого слоя аппроксимируются с помощью периодических B-сплайнов (Рис. 1).
Таким образом, получаем достаточно точное приближение для контура на каждом слое
томограммы. B-сплайны параметризованы по длине, поэтому в процессе триангуляции генерируются
треугольники из соответствующих точек B-сплайнов. В данном подходе приходится иметь дело с
неоднозначностями. Обычно, неоднозначности в контурном анализе разделяют на 3 группы.
Вопервых – это проблема нахождения когерентных контуров. Имея несколько контуров на
предыдущем и последующем слое, необходимо определить какие, контуры соединять друг с другом.
Вовторых, необходимо найти соответственные точки на когерентных контурах. В-третьих,
необходимо разрешить проблему ветвления в случае, если, число контуров меняется от слоя к слою.
Решение для всех этих проблем предлагается, например, в статьях [
Рис. 1. Результат аппроксимации контуров для всех слоев Аналогичные операции производятся в ортогональном направлении. Из-за использования аппроксимации для приближения контура, в общем случае, нельзя получить точек пересечения ортогональных контуров. Поэтому, для того, чтобы получить сетку, учитываются точки уже сгенерированных контуров для сплайн-аппроксимации в ортогональном направлении (Рис. 2). Рис. 2. Результат аппроксимации контуров для всех слоев в двух ортогональных направлениях 3. Реконструкция методом объединения контуров в ортогональных плоскостях
Первым этапом реконструкции является получение набора точек. В данном подходе, будем рассматривать поверхность, интерполирующую данные точки, то есть каждая вычисленная точка будет являться вершиной поверхности. Рассматривая произвольный набор точек в пространстве, возникает ряд неоднозначностей, наиболее важной из которых является определение соседних вершин.
Для разрешения этих неоднозначностей алгоритм Marching Cubes предлагает принять следующие условия: • В пространстве рассматривается прямоугольная сетка; • Ребро ячейки сетки пересекает только одну изоповерхность; • Вершины изоповерхности находятся на узлах и ребрах сетки; • Вершины, находящиеся на грани ячейки, являются соседними (их можно соединить линией).
Данные свойства дают возможность применить несколько иной подход к реконструкции –
вначале реконструировать кривые (или контуры) исходной поверхности. Данный подход является
широко распространенным, и, например, хорошо описан в [
Алгоритм реконструкции состоит из трех этапов. На каждом из них вычисляются вершины и выявляются связи между ними, основываясь на установленном критерии. В процессе работы алгоритма последовательно рассматриваются 2D сетки, ортогональные осям. Для оси X получаем последовательность контуров (Рис. 3, слева) и множество вершин на ребрах сеток. Затем рассматриваются 2D сетки, ортогональные оси Z. Множество вершин на сетках, полученное на предыдущем этапе, переносится на текущие сетки, и вычисляются новые вершины. После выявления связей получаем последовательность контуров (Рис. 3, в середине). В итоге, рассматриваются 2D сетки, ортогональные оси Y. Вычисление новых вершин не требуется, и остается перенести полученные на предыдущих этапах вершины на текущие сетки. Получаем последовательность контуров (Рис.3, справа).
Рис. 3. Результат работы алгоритма Последовательность направлений может быть любой. Но только выбор по трем направлениям гарантирует единственность получаемой модели. Сетки в пределах одного направления обрабатываются независимо, что дает возможность распараллелить вычисления.
Таким образом, основными этапами алгоритма являются реконструкция сетки поверхности и триангуляция полученной сетки. 4. Использованные инструменты и реализация
Для реализации сегментации и выделения контуров использовались библиотеки ITK и VTK.
Данные библиотеки с открытым исходным кодом и распространяются под лицензией Apache 2.0.
Лицензия Apache дает пользователю право использовать программное обеспечение для любых
целей, свободно распространять, изменять, и распространять измененные копии, за исключением
названия. Библиотека ITK содержит большое количество реализаций для различных методов
сегментации, причем сегментация может быть произведена в пространстве любой размерности.
Библиотека работает по принципу конвейера, то есть данные для сегментации и любых других
операций загружаются в память порциями, размер которых определяется пользователем. Таким образом,
можно работать с данными очень больших объемов и производить параллельные вычисления.
Например, для сегментации методом разрастания регионов достаточно выбрать начальную точку
алгоритма и указать, что сегментация будет произведена в трехмерном пространстве. Для
сегментации методом активных контуров необходимо задать начальный контур. Для анализа была взята
томограмма мозга, которую можно бесплатно получить в Интернете в рамках проекта brainweb
[
Рис. 4. Результат работы алгоритма.
После выделения точечных контуров, они аппроксимируются с использованием B-сплайнов.
Для этой задачи используется библиотека Geometric Tools [
B-сплайны параметризованы по длине. Для простых структур, у которых форма и размеры контура не сильно отличаются от слоя к слою, этого оказывается достаточно для нахождения соответственных точек и триангуляции.
Библиотека VTK служит для визуализации результатов обработки изображений. К сожалению, на данный момент, она не содержит механизмов триангуляции параллельных контуров, и, тем более, разрешения неоднозначностей.
Реализацию второго подхода рассмотрим на примере трехмерной прямоугольной сетки 6x5x5 с известными значениями поля в узлах. Для построения модели сначала рассматриваются 2D сетки, ортогональные оси Z.
Для текущей 2D сетки, введем нумерацию узлов (совпадает с нумерацией вертикальных ребер) и горизонтальных ребер (Рис. 6).
Рис. 6. Нумерация узлов Для каждого ребра значения поля в инцидентных узлах сравниваются с изозначением. Если в одном узле значение больше, а в другом – меньше изозначения, то через это ребро проходит вершина поверхности. В этом случае интерполируется положение вершины (Рис. 7). Возможно, изозначение совпадет со значением поля в узле сетки. В этом случае положение вершины будет точно на узле.
Рис. 7. Интерполированные положения вершин. Таким образом, алгоритм работает с тремя типами вершин: 1) Вершины, лежащие на горизонтальных линиях сетки; 2) Вершины, лежащие на вертикальных линиях сетки; 3) Вершины, лежащие в узлах сетки. Вместе c типом вершины, соответствующий номер ребра определяет ее положение на сетке. В дальнейшем, номер будет использоваться при вычислении соседних вершин.
Для определения соседних вершин предлагается использовать критерий: связывать вершины, находящиеся на ближайших ребрах (Рис. 7). Для горизонтальных и вертикальных вершин можно выделить три возможных комбинации размещения соседей (Рис. 8).
Рис. 8. Вычисление связей для вершин. Для каждой вершины обозначена окрестность, на границе которой возможные вершины являются соседними.
Обращение к соответствующим ребрам для проверки наличия на них соседних вершин осуществляется через введенную нумерацию вертикальных (узловых) и горизонтальных вершин. Например, на рисунках 9 и 10 изображен способ вычисления номеров ребер для вершины на ребре J.
Рис. 9. Вычисление номеров нижнего и верхнего ребра, W – ширина сетки (в ребрах) Рис. 10. Вычисление номеров для левых и правых ребер, R = J / W Рис. 11. Построенный контур Таким образом, строим последовательность контуров, ортогональных Z (Рис. 12).
Рис. 12. Последовательность контуров, ортогональных оси Z Можно заметить, что вертикальные и узловые вершины, принадлежащие сеткам, ортогональным Z, будут принадлежать и сеткам, ортогональным X. Эти вершины (вместе с информацией об их связях) переносятся на сетки, ортогональные X (предполагается использование индекса вершины или указателя), и там, аналогично, будут являться вертикальными или узловыми. Для переноса, выполняется обращение к соответствующим ребрам с помощью нумерации. Также, вычисляются новые вершины (На рисунке 13 – горизонтальные). После этого применяется описанный ранее алгоритм вычисления соседних вершин к каждой сетке. Рис. 13. Последовательность контуров, ортогональных оси X На последнем этапе вычислять новые вершины не требуется – они вычислены на предыдущих этапах. Их необходимо правильно перенести – вершины, вычисленные на первом этапе – останутся горизонтальными, а вершины со второго этапа – были горизонтальными, но станут вертикальными. Аналогично, к каждой сетке применяется алгоритм вычисления связей (Рис. 14) Рис. 14. Последовательность контуров, ортогональных оси Y В полученном графе, степень вершины и количество ребер в каждой внутренней грани, как правило, ограничены константой. Это дает возможность триангулировать данный граф, например, посредством перебора инцидентных граней для каждой вершины. 5. Результаты
Для первого подхода, на четырехъядерном процессоре AMD A-10 5745M 2.10GHz с 12 Гб оперативной памяти были получены следующие результаты.
Таблица 1. Результаты сегментации Число итераций лельных вычислений высока, было получено небольшое ускорение. В дальнейшем, планируется распараллелить этот процесс с помощью OpenMP, что должно повысить производительность.
Таблица 2. Результаты: выделение контуров, аппроксимация и триангуляция (для всех слоев) Ускорение, полученное для данной фазы, согласуется с теоретическими оценками для доли параллельного кода около 75% (закон Амдала). В дальнейшем, планируется оптимизировать алгоритм для повышения доли параллельных вычислений.
Для второго подхода, на двухъядерном процессоре Intel®Core™ i5 750 2.66 GHz с 4 Гб оперативной памяти были получены следующие результаты.
Таблица 3. Таблица производительности для алгоритма вычисления сетки контуров Долю последовательного кода в данном алгоритме можно оценить отношением объема томограммы к числу полученных вершин. Таким образом, из-за небольшого количества вершин в изоповерхности по сравнению с количеством обрабатываемых ребер, и как следствие, малой долей последовательного кода, теоретическое ускорение, по закону Амдала, почти равно числу процессоров. Однако, на практике, из-за неравномерной загрузки и наличия кэш-промахов в функции вычисления связей, полученное ускорение оказалось меньше теоретического.
Далее показаны таблица и график для сравнения с реализацией Marching Cubes [
В данной статье представлены два подхода, применяющиеся для реконструкции поверхностей
произвольного топологического типа. Вначале реконструировались контуры поверхности.
Показано, как наличие ортогональных плоскостей позволяет уточнять текущие и строить новые контуры.
В рамках первого подхода, трехмерная сегментация органов и автоматическое выделение контуров
производились с помощью библиотеки ITK. С целью получения гладкого параметризованного
представления выполнена аппроксимация контуров на каждом слое периодическими B-сплайнами
третьей степени. Во втором подходе, реконструкция поверхности выполнялась за счет построения
ее графа. Вычисление координат вершин поверхности производилось аналогично алгоритму
Marching Cubes [
В процессе реконструкции, плоскости обрабатываются независимо. Это позволило получить достаточно высокие ускорения в параллельных версиях алгоритмов. В дальнейшем, после получения замкнутой поверхности в рамках первого подхода, будет произведена оценка нормалей во всех точках, участвовавших в триангуляции. Затем будет применена технология криволинейных PN-треугольников для уточнения поверхности. Используя средства библиотеки ITK, можно рассчитывать площади построенных поверхностей и их объемы. Также планируется использовать аппроксимацию всего набора контуров с помощью В-сплайновых поверхностей.
Для структуры, полученной во втором подходе, будет разработан более эффективный
алгоритм триангуляции, позволяющий получить приближение поверхности криволинейными
PNтреугольниками [
In this paper we solve the reconstruction problem of external and internal borders of organs with a complex spatial structure (brain, complex fracture). Reconstruction performed in orthogonal planes to resolve ambiguities. First approach is the following: at the beginning we perform the 3-dimensional segmentation using region growing algorithm. Then we solve the 2dimensional contour reconstruction problem for the two orthogonal cuts. In the second case, the algorithm reconstructs contours in planes which are orthogonal to the coordinate axis, and combines them into graph of surface. Further, in the both approaches, resulting set of contours are triangulating. Due to the independent processing of the planes in the parallel implementation achieved a sufficiently high acceleration.