CC BY
28
Cloud of Science. 2018. T. 5. № 1 http:/ / cloudofscience.ru
Специальные процедуры для работы с объектами со сложной внутренней структурой по стеку КТ-сканов1
Э. П. Шурина*,**,***, Д. В. Добролюбова*,**,***, Е. И. Штанько*,***
*ФГБУН Институт нефтегазовой геологии и геофизики им. А. А. Трофимука Сибирского отделения Российской академии наук 630090, Новосибирск, пр-т Академика Коптюга, 3
**ФГБОУ ВО «Новосибирский государственный технический университет» 630073, Новосибирск, пр-т К Маркса, 20
***ФГАОУ ВО «Новосибирский национальный исследовательский государственный университет» 630090, Новосибирск, ул. Пирогова, 2
e-mail: [email protected]. ru
Аннотация. В работе предлагается алгоритм построения специальной иерархической структуры представления данных КТ- или МРТ-сканирования для искусственных или нативных сред со множеством сложных внутренних границ. Выполняется конвертация стека послойных изображений в градациях серого в сеточное разбиение, состоящее на верхнем уровне иерархии из полиэдральных (полигональных для 2D случаев) макроэлементов, а на микроуровне — из тетраэдров (треугольников для 2D). Полученные разбиения ориентированы на приложение в конечноэлементном многомасштабном моделировании и могут быть применены для решения широкого круга многофизичных задач. Ключевые слова: симплициальное разбиение, полигональное разбиение, искусственные и нативные среды.
1. Введение
В последние десятилетия различные методы визуализации внутренней структуры сложных объектов активно применяются практически во всех прикладных, научных и медицинских сферах. Среди таких методов наиболее распространены компьютерная томография (КТ), мультиспиральная компьютерная томография (МСКТ) и магниторезонансная томография (МРТ). Наиболее часто подобные технологии
1 Работа выполнена при поддержке комплексной программы фундаментальных исследований СО РАН «Междисциплинарные интеграционные исследования» на 2018-2020 гг. по теме «Экспериментальные исследования и математическое моделирование нативных и инженерных объектов с фазоизменяемы-ми материалами» и программы фундаментальных исследований Президиума РАН № 27 «Фундаментальные проблемы решения сложных практических задач с помощью суперкомпьютеров».
используются в медицине при диагностических исследованиях [1-3]. В материаловедении и геофизике КТ-установки широко применяются для исследования внутренней структуры искусственных [4, 5] и нативных [6] сред соответственно. В качестве инженерного приложения можно назвать, в частности, неразрушающий контроль качества и диагностику скрытых дефектов при производстве технических изделий [7], например полупроводников и изделий авиационной техники из полимерных материалов.
Визуализирующие методы исследования позволяют выполнить оцифровку внутренней структуры образца, которая может быть представлена в виде стека двумерных послойных изображений объекта.
На данный момент предлагается большое количество коммерческих и свободно распространяемых программных пакетов с широким спектром возможностей для обработки получаемого стека изображений. Они ориентированы в основном на построение трехмерной интерактивной модели с визуализацией внутренней структуры объекта и на конвертацию полученной трехмерной модели в формат, пригодный для вывода на 3Б-печать.
Среди предлагаемых пакетов можно выделить несколько основных, обладающих наиболее широким спектром возможностей. Avizo [8] является мощным коммерческим продуктом, ориентированным на такие прикладные области, как материаловедение, геофизика, медицина, контроль качества полупроводников. Avizo позволяет выполнять построение трехмерных моделей по двумерным сканам, производить развернутый статистический анализ внутренней структуры объекта, включающий в себя кластеризацию, исследование распределения включений, определение объемных концентраций и ориентации включений. В ряде случаев имеется возможность строить трехмерные сетки, однако для сравнительно простых внутренних структур, когда удается сгенерировать замкнутые поверхностные сетки без взаимных пересечений и касаний. Основным коммерческим конкурентом Avizo являются Simpleware [9] от Synopsys с аналогичным набором функций. Другой популярный коммерческий визуализатор Materialise Mimics [10] ориентирован главным образом на медицинские приложения и обработку данных МРТ и МСКТ.
Среди свободно распространяемых пакетов можно выделить 3 основных: 3D Slicer [11], ITKsnap [12] и FIJI [13] (на основе ImageJ). ITKsnap позволяет осуществлять только визуализацию и навигацию по трехмерной модели внутри своей пользовательской оболочки. Экспорт сеток не предусмотрен, пакет ориентирован на медицинские приложения. 3D Slicer предназначен, главным образом, для адаптации воссозданных трехмерных моделей для вывода на 3Б-печать. Наиболее широкие возможности предоставляет FIJI, который позволяет осуществлять сборку 3Б-модели из стека сканов и выполнять некоторую простейшую статистическую
обработку. Предусмотрен и вывод геометрии в формате, пригодном для сеточных генераторов, однако данная функция в полной мере не работает для объектов со сложной внутренней структурой, такой, например, как поровое пространство керна или трабекулярная структура кости.
Несмотря на большое количество коммерческих и свободно распространяемых программных продуктов, позволяющих обрабатывать данные послойного сканирования объектов, проблема остается актуальной и востребованной. Во многих пакетах заявляется возможность построения конечноэлементной сетки для полученных моделей, однако именно эта проблема в полной мере предлагаемыми пакетами не решается. Зачастую генерируются поверхностные сетки, требующие дальнейшей ручной доработки, так как они содержат большое количество вырожденных элементов и сеточных «дыр». Несмотря на то, что в визуализируемой пакетами трехмерной модели обычно присутствует разделение по подобластям с различными физическими свойствами, в явном виде на сеточном уровне выделить области с различными физическими характеристиками, особенно если речь идет о микровключениях, предлагаемые программные пакеты не позволяют.
Разработка научных программных продуктов налагает свои ограничения на формат выходных данных. В частности, если конечной целью является не визуализация внутренней структуры объекта, а решение задач математической физики в данной расчетной области, то возникают определенные требования, тесно связанные как с постановкой решаемой задачи, так и с применяемыми методами.
Для конечноэлементных методов в любом физическом приложении всегда является актуальной проблема генерации оптимальной сеточной дискретизации. Выбор алгоритмов построения конечноэлементных разбиений на основе данных КТ и МРТ сканирования обуславливается в первую очередь особенностями моделируемых сред и конкретной областью приложения. Следует также учитывать, что послойное снятие изображений предполагает высокое разрешение каждого из срезов, а также большое количество получаемых изображений с малым (порядка нескольких микрон) шагом дискретизации, как следствие, возникает необходимость обрабатывать и адаптировать для построения сетки большие объемы данных.
Наиболее распространена задача построения конечноэлементной модели на основе набора изображений в биомедицинских приложениях, где исследуемые среды имеют сложную и разномасштабную внутреннюю структуру. Для корректного учета структуры биологических тканей необходимо обеспечивать адаптивность генерируемых конечноэлементных разбиений.
В [1] Zhang Y. и соавторы на основе данных КТ- и МРТ-сканирования биологических тканей мозга, сердца и костной ткани осуществляют построение трехмерного тетраэдрального и гексаэдрального адаптивного конечноэлементного разбиения, которое
генерируется с помощью метода dual contouring и оптимизируется с помощью алгоритмов сглаживания и стягивания ребер.
В [2] предлагается алгоритм построения трехмерного конечноэлементного разбиения, основанный на методе шагающих тетраэдров (Marching tetrahedra). В качестве области приложения выступают биологические ткани, выполняется сегментация томографического изображения области грудины мыши. Авторы обращают внимание на то, что в многофазной среде, которую представляет собой живая ткань, необходимо правильно учитывать внутренние границы между органами, обладающими различными физическими свойствами. Это обеспечивается за счет дополнительных подразбиений первичной равномерной структурированной тетраэдральной сетки таким образом, чтобы каждый элемент результирующей сетки целиком принадлежал подобласти с постоянными физическими свойствами.
Построение конечноэлементных разбиений образцов губчатой костной ткани также требует тщательного учета всех особенностей внутренней структуры [3]. Для КТ-изображения тонких губчатых тканей особенную сложность представляет сохранение связности структур при построении конечноэлементной модели. Авторами используется гистерезисная фильтрация с динамически изменяющимися пороговыми значениями [14], позволяющая учесть сложную связность трабекулярных структур в расчетной области и ограниченная триангуляция Делоне для построения внутреннего конечноэлементного разбиения. Конечноэлементный анализ на основе построенной модели выполняется коммерческим программным продуктом ABAQUS [15].
Зачастую даже использование коммерческих программных пакетов требует достаточно трудоемкой предварительной подготовки входных данных. Предложенный в [16] алгоритм построения конечноэлементного разбиения основан на применении комбинации программных пакетов, таких как коммерческие программные продукты Mimics, MSC.Patran и Ansys, а также свободно распространяемый пакет Bonemat. В [16] отмечается, что представление входных данных в распространенном в медицинских измерениях формате DICOM требует их подготовки, при этом входные и выходные форматы используемого ПО могут также оказаться несовместимыми, что порождает необходимость дополнительной обработки данных. В материаловедении неразрушающие методы исследования получили широкое применение. В [4] предлагается многоступенчатый алгоритм построения ко-нечноэлементной модели углепластиков. Схожие плотности эпоксидных смол и карбоновых волокон, составляющих углепластики, обуславливают низкий контраст матрицы и включений на КТ-сканах, что может привести к неверному учету внутренней структуры. Стандартные методы обработки изображений не обеспечивают восстановление непрерывной структуры углеродных трубок, поэтому в [4] предла-
гается двухшаговый алгоритм, который позволяет в несколько этапов восстановить их внутреннюю структуру, для которой затем с помощью пакета ABAQUS строится конечноэлементное разбиение.
Предложенный алгоритм основан на Баейсовском выводе и учитывает наличие в образце кластеров углеродных трубок с различной ориентацией и другие особенности углепластиков. Отмечается, что разбиение модели на кластеры с различной ориентацией трубок необходимо, поскольку конечноэлементная сетка во всем образце будет содержать десятки миллионов элементов, что значительно усложнит численное моделирование.
Особенности структуры расчетной области определили подход к построению конечноэлементной модели теплозащитной системы многоцелевого командного модуля Орион [7]. Моделирование термозащитной системы требует серьезных вычислительных затрат, поэтому конечноэлементная модель строится только в части расчетной области. Особенности изображений, полученных в результате томографии, обусловили необходимость применения алгоритма шумоподавления, что позволило корректно учесть топологическую структуру области. В качестве фильтра использован фильтр анизотропной диффузии, который эффективно подавляет шум на изображении, но не сглаживает границы подобластей. Адаптивное разбиение выполнено пакетами ScanFE и Simpleware.
Наличие некоторой априорной информации о внутренней структуре исследуемой области позволяет упростить применяемые алгоритмы. В [5] при построении конечноэлементного разбиения на основе стека КТ-изображений образца бетона, содержащего небольшое число трещин, на этапе подготовки изображения использован алгоритм водораздела для сегментации изображения и оператор Кэнни для выделения контуров трещин. Построение трехмерного тетраэдрального разбиения расчетной области осуществляется алгоритмом Делоне по выделенным трехмерных изоповерхностям, являющимся трехмерными внутренними границами «скелет-трещина».
В работе выполняется сравнение двух подходов к построению сеточных разбиений для образца со сложной внутренней геометрией, близкой по своим свойствам к геологической среде.
Введенная специальная иерархическая сеточная структура ориентирована на применение в многомасштабном моделировании для решения многофизичных задач и определения эффективных характеристик нативных и искусственных сред. Во избежание резкого роста размерности сеточного разбиения, предлагается выполнять макро-элементную дискретизацию расчетной области. Предлагается использовать в качестве макроэлементов полигоны и полиэдры для 2D- и 3D-постановок соответственно, при-
чем для каждой из введенных макроподобластей осуществляется дополнительная внутренняя симплициальная триангуляция.
2. Постановка задачи
Для сохранения высокой информативности КТ- и МРТ-исследований сканирование выполняется в каждом слое с высоким разрешением и малым пространственным шагом по Z, особенно в инженерных приложениях. Разрешение снимка слоя может достигать порядка нескольких тысяч пикселей по каждой из осей X и У, причем количество слоев по Z может быть такого же порядка. При этом практически единственной информацией, которой может оперировать исследователь, является интенсивность цвета в пикселе изображения. Интенсивность цвета в данном случае представляет собой целочисленную величину, варьирующуюся в пределах от 0 до 255: рг(х,у) е [0,255].
В работе предлагается алгоритм конвертации данных КТ- либо МРТ-сканирования в сеточную форму для выполнения конечноэлементного многомасштабного моделирования. Разработанные программные комплексы позволяют подготовить макроэлементное разбиение образца, внутренние симплициальные разбиения с подразделением сеточных элементов на подобласти в соответствии с некоторым введенным критерием кластеризации цветовой интенсивности р.
Предполагается, что введенная кластеризация сопоставима с разбиением на физические подобласти с некоторой постоянной физической характеристикой, что позволяет однозначно идентифицировать скелет и включения различных типов. Основные этапы предлагаемого иерархического алгоритма приведены на рис. 1.
Алгоритм может применяться на уровне каждого послойного изображения и на уровне всего трехмерного образца. Покажем особенности его реализации для сечения образца.
3. Специальные процедуры исследования внутренней структуры отдельных сканов
Для исследования внутренней структуры на уровне микроразбиения в работе реализованы два различных алгоритма. Оба алгоритма в качестве подготовительного этапа реализуют построение неструктурированного симплициального разбиения
однородного макроэлемента М1 : Т = | /(}. где Ь — количество симплексов для макроэлемента М1. Алгоритм I включает в себя этапы, приведенные на рис. 2.
Рисунок 1. Иерархический алгоритм. Основные этапы
■ Для каждого /,67" определить значения
= 1, 5, где Р) -внутренние точки симплекса
Согласно выполненной ранее кластеризации определить С1ь к
= 1, Л! для каждого полученного значения рЛ(ЯД,/ = 1,5. Тогда П^ которому принадлежит большее число точек Pj
Сглаживание. Для r„ tt, Ii, lifMff-0, lf\t&<b, tf\lfi0, определить ii». k=l,N. Если 2 и более соседних симплексов принадлежат ih, отличной от определенной на этапе 2 для симплекса h области Ц,. то tiüilh. иначе ¡¡BiX/
Рисунок 2. Алгоритм I
Построение сеточного неструктурированного симплициального разбиения в алгоритме I выполняется свободно распространяемым пакетом Gmsh [17] в каждом макроэлементе.
Алгоритм II основан на определении границ раздела сред внутри расчетной области и последующем построении внутри каждого макроэлемента геометрической модели с помощью программного комплекса SALOME™ [18]. Затем выполняется адаптивная симплициальная дискретизация этой модели. Данный алгоритм состоит из этапов, показанных на рис. 3.
Построение внутренних границ St, i = 1, N подобластей в пределах макроэлемента Mi реализовано на основе алгоритма шагающих тетраэдров [2, 19]. При этом в качестве исходного разбиения выбирается неструктурированное симплициальное разбиение, введенное на макроэлементе на этапе 1. Граница раздела подобластей представляет собой трехмерную изоповерхность, для которой выполняется условие: Vxj е S PZ (Xj ) = Р , т. е. значение функции интенсивности цвета в любой точке изоповерхности постоянно и равно некоторому пороговому значению р .
Выделение внутренних
I
межфрагментарных границ Як JB
5,,/ = 1.М с помощью алгоритма шагающих т w
тетраэдров
Определение характеристик подобластей Лк,
* lpZk(n,)
ограниченных замкнутыми Pin ' Pout
контурами 5(1/ = , на В »
>2'(Пг)
основе выполненной ранее
кластеризации; щ
Построение адаптивного симплициального ЩЕШШ
разбиения для каждого
макроэлемента М^М на - cNL. * /
основе информации о границах раздела сред и характеристиках подобластей Пк
Рисунок 3. Алгоритм II
Основная идея алгоритма шагающих тетраэдров состоит в том, чтобы строить границу раздела подобластей в каждом тетраэдре исходного разбиения области. Считается, что граница раздела подобластей проходит через тетраэдр ^, если хотя
бы в одной из его вершин V] е ^, ] = 1,..., 4 функция р2к (V,-) принимает значение, отличное от значений в остальных его вершинах. Сечение тетраэдра изоповерхно-стью описывается одним из 16 возможных вариантов, при этом два из них тривиальны: значение функции р2к (V]) одинаково во всех узлах тетраэдра. Из 14 нетривиальных вариантов уникальны только два (рис. 4). Точки р е ^, в которых изо-поверхность пересекает ребра тетраэдра, вычисляются линейной интерполяцией: Р&' = VI + (Р - р(У ))(У2 - VI V (р( VI) - р(У2)), где VI и V 2 — вершины ребра тетраэдра.
а) б)
Рисунок 4. Варианты сечения тетраэдра изоповерхность: а) значение р2к (V]) отлично в одной вершине; б) значение р2к (V]) отлично в двух вершинах
Алгоритм шагающих тетраэдров позволяет построить топологически корректные изоповерхности, не содержащие «дыр». Однако в случае, когда границы подобластей пересекают границы макроэлементов, построенные в результате его применения поверхности (контуры в 2D) могут оказаться незамкнутыми. Для корректного учета пересечения межфрагментарных границ и границ полиэдрального макроэлемента алгоритм шагающих тетраэдров организован в рамках иерархического подхода, предполагающего построение замкнутых границ раздела сред на каждом уровне иерархии (рис. 5).
Рисунок 5. Иерархия выделения межфрагментарных границ
Для каждого элемента Kj межфрагментарной границы S вычисляются значения функции интенсивности цвета pZi и pz0lut по разные стороны от данной границы. На этапе 2 алгоритма II подобластям, разделенным построенной границей раздела, присваиваются значения функции интенсивности цвета, вычисленные осреднением по K ■. На основании проведенной ранее кластеризации делается вывод о
принадлежности фрагмента расчетной области той или иной подобласти Пк. Очевидно, что данный способ определения характеристик фрагментов подобластей в макроэлементе не дает корректного результата в расчетных областях, в которых отсутствует явное разделение на вмещающую среду и включения, а также в областях, содержащих сложные, вложенные одно в другое, включения. Для использования алгоритма II в таких расчетных областях требуется модификация данного этапа алгоритма.
Результатом выполнения этапов 1-3 алгоритма II является файл для Python-интерпретатора модуля GEOM программного пакета SALOME, в результате испол-
нения которого строится геометрическая модель расчетной области, а затем выполняется адаптивная симплициальная дискретизация.
Верификация предложенных алгоритмов выполняется на тестовой области с регулярной внутренней структурой (рис. 6), для которой можно точно определить такие характеристики, как длину границы раздела сред и концентрацию включений относительно общей площади образца. Так как разработанные алгоритмы ориентированы на формирование набора входных данных для решения задач определения эффективных характеристик сложных сред многомасштабными конечноэлемент-ными методами, то выбор предложенных критериев сравнения можно считать оптимальным, поскольку, согласно выполненным исследованиям [20], именно они оказывают основное влияние на значение эффективной характеристики среды.
а) 3Б б) 2Б
Рисунок 6. Эталонная область
Приведенная двумерная область (рис. 6б) соответствует ХГ-сечению трехмерного образца (рис. 6а) и может быть интерпретирована как сегментированное изображение КТ-сканирования (скан 84 из стека, состоящего из 170 изображений). Разрешение скана 733^733 пикселя. Для рассмотренной области общая длина границы раздела сред и концентрация включений относительно общей площади образца могут быть вычислены аналитически, исходя из того, что диаметры включений в данном сечении составляют 85 пикселей. В предположении, что область О целиком является единственным элементом макроразбиения, выполним исследование качества аппроксимации длины границы раздела сред й и площади включений на внутреннем микроразбиении предложенными алгоритмами I и II (табл. 1).
В табл. 1 N — размер сеточного разбиения, d — рассчитанная длина границы раздела сред, d * — точное значение длины границы, полученное для эталонного изображения, £ — рассчитанная концентрация включений в масштабах всего сечения, Б * — точное значение концентрации включений для эталонного скана.
На рис. 7 показаны восстановленные внутренние структуры 2D-образца: построенные сеточные разбиения и укрупненные изображения подобласти сеточного разбиения, содержащей одно включение, полученные алгоритмом I (со сглажива-
нием и без него) и алгоритмом II. Алгоритм II динамически изменяет исходное сеточное разбиение и строит адаптивную сетку. Алгоритм I позволяет получить неструктурированное, но не являющееся адаптивным, разбиение.
Таблица
N Алгоритм I Алгоритм 2
Без сглаживания Со сглаживанием
| й - ¿С | |5 - 5 *| | й - й | | 5 - 5 *| | й - й* | | 5 - 5 *|
й * 5 * й * 5 * й * 5 *
8 332 0.2977 0.0273 0.1057 0.0084 0.1313 0.0868
12 365 0.2928 0.0273 0.1076 0.0210 0.1533 0.0659
35 624 0.3752 0.0347 0.1114 0.0210 0.1152 0.0157
48 914 0.3176 0.0357 0.1184 0.0263 0.1305 0.0094
80 970 0.3457 0.0326 0.1211 0.0284 0.1740 0.0052
Для алгоритма I рост подробности сеточного разбиения приводит к увеличению погрешности аппроксимации как длины границы раздела сред, так и концентрации включений. Для алгоритма II также наблюдается увеличение погрешности аппроксимации длины границы раздела сред. Это связано с повышением уровня осцилляций функции интенсивности цвета р11 (х, у) вблизи границы раздела сред и активно проявляется в случае гладких границ.
а) алгоритм I без сглаживания
б) алгоритм I со сглаживанием
в) алгоритм II
Рисунок 7. Симплициальные разбиения восстановленной области
Далее рассмотрен образец среды, близкой по своей внутренней структуре к некоторой нативной (рис. 8а). Данная область содержит три типа включений, дифференцируемых на сканах по различной насыщенности цвета. Включения имеют произвольную форму и негладкие границы. Стек сканов содержит 400 изображений с разрешением 733^733. На рис. 8б приведен скан 245 для образца (рис. 8а).
Так как двумерная область моделирования (рис. 8б) имеет сложную внутреннюю структуру, для нее было простроено полигональное макроразбиение (рис. 8в). Обработка полигонов выполняется независимо, предложенные алгоритмы обладают свойством естественного параллелизма. Параллельная реализация алгоритмов позволяет значительно ускорить процесс распознавания внутренней структуры всего образца. Так, для алгоритма I время обработки скана уменьшается в 4.3 раза — со 101 секунды до 24 секунд2 для сеточных разбиений, содержащих 12 полигональных элементов и порядка 900 тыс. треугольных элементов микроразбиения во всем образце.
Р
а) область моделирования б) скан 245 в) макроразбиение
Рисунок 8. Расчетная область с нерегулярными включениями, близкая к нативной среде
а) алгоритм I
б) Алгоритм II
Рисунок 9. Восстановленная внутренняя структура области. Отдельные полигоны
! Вычисления выполнялись на ПК: процессор Intel Core i7-3930K 3.20 GHz (6 ядер); оперативная память 48 Гб.
Алгоритм I требует построения достаточно подробного внутреннего разбиения (рис. 9а около 80 тыс. симплексов). Алгоритм II позволяет сгенерировать более грубое разбиение, адаптивно сгущающееся к границе раздела сред (рис. 9б около 3.5 тыс. элементов). Однако даже в случае мелкого сеточного разбиения, при работе на отдельных макроэлементах, такая высокая подробность не представляет серьезной вычислительной сложности, несмотря на то, что в масштабах всего образца количество треугольников может достигать порядка миллиона.
4. Специальные процедуры исследования внутренней структуры всего образца
Для трехмерного образца с регулярно расположенными включениями (рис. 6а) были выполнены оценки объемной концентрации и площади поверхности включений (рис. 10).
а) алгоритм I б) алгоритм II
Рисунок 10. Восстановленная структура эталонной области
Аналогично двумерному случаю, погрешность аппроксимации концентрации и площади поверхности включений, восстановленных при помощи алгоритма I, возрастает при увеличении размерности исходного тетраэдрального разбиения. Так, при N ~ 2 ■ 105 погрешность аппроксимации объемной концентрации включений составляет 0.012 и 0.008 без сглаживания и со сглаживанием соответственно, при N «106 погрешности составляют соответственно 0.066 и 0.041. При использовании алгоритма II достаточно грубое исходное разбиение N = 43 032 позволяет восстановить объемную концентрацию включений c погрешностью 0.0590, а более подробное разбиение N = 215 613 — с погрешностью 0.0507.
Аппроксимация площади поверхности включений с увеличением подробности исходного разбиения ухудшается для обоих алгоритмов. Чтобы избежать такого эффекта, необходимо выполнять дополнительное сглаживание границы раздела сред, основанное на загрублении этой поверхности и удалении острых углов и «провалов».
Для алгоритма II повышение подробности исходного тетраэдрального разбиения приводит к увеличению времени обработки построенных границ раздела сред на этапе 3, так как они формируются из отдельных симплексов, для алгоритма I время увеличивается незначительно. Использование макроразбиения трехмерных образцов позволяет повысить точность аппроксимации внутренней структуры и скорость работы обоих алгоритмов, поскольку восстановление межфрагментарных границ и последующая обработка макроэлементов выполняются параллельно и независимо.
Для среды, приближенной к нативной (рис. 8а), восстановление внутренней структуры производилось на полиэдральном разбиении (рис. 11а). Для упрощения визуального восприятия приведем фрагменты восстановленных включений для одного полиэдра обоими алгоритмами (рис. 11 б—в).
а) полиэдральное разбиение расчетной области 8а
б) алгоритм I в) Алгоритм II
Рисунок 11. Восстановленная структура фрагмента расчетной области, принадлежащего одному из элементов макроразбиения
Подробность внутренних разбиений для алгоритмов I и II значительно различается, однако количество макроэлементов остается одинаковым. Требования к макроэлементной сетке налагаются исходя из геометрических особенностей рас-
четной области, а именно объемной концентрации включений, их расположения и кластеризации внутри образца.
5. Заключение
В работе реализован алгоритм конвертации стека послойных изображений натив-ной или искусственной среды со сложной внутренней структурой в специальную иерархическую сеточную структуру, ориентированную на решение широкого класса многофизичных задач многомасштабными методами. Алгоритм применим как для обработки отдельных сканов, так и для работы с трехмерным образцом. Идентификация внутренней структуры объекта на каждом из полиэдральных (полигональных для 2D случая) макроэлементов выполняется реализованными алгоритмами I и II, для которых выполнена верификация и исследование области применимости. Алгоритм I, основанный на сопоставлении значений интенсивности цвета изображения элементам нерегулярного симплициального разбиения, приводит к построению излишне мелких сеточных разбиений, с множественными осцилляция-ми на границе раздела сред. Алгоритм II показывает более высокую эффективность при оценке концентрации включений по отношению к объему/поверхности образца и длины границы раздела сред. Критическим параметром для алгоритма II является концентрация включений; для обработки областей со сложными вложенными включениями требуется модификация этапа 2.
Разработанные алгоритмы являются первым этапом численного определения эффективных характеристик сред, для которых выполняется неразрушающая визуализация внутренней структуры.
Литература
[1] Zhang Y., Bajaj C., Sohn B. S. 3D finite element meshing from imaging data // Computer methods in applied mechanics and engineering. 2005. Vol. 194, No. 48-49. P. 5083-5106.
[2] Cong A. et al. Geometrical modeling using multiregional marching tetrahedral for bioluminescence tomography // Medical Imaging 2005: Visualization, Image-Guided Procedures, and Display. International Society for Optics and Photonics, 2005. Vol. 5744. P. 756-764.
[3] Chen C. et al. Assessment of trabecular bone strength at in vivo ct imaging with space-variant hysteresis and finite element modelling // 2016 IEEE 13th International Symposium on. IEEE Biomedical Imaging (ISBI), 2016. P. 872-875.
[4] Sencu R. M., Yang Z., Wang Y., Withers P., Rau C., Parson A., Soutis C. Generation of Micro-scale Finite Element Models from Synchrotron X-ray CT Images for multidirectional Carbon Fibre Reinforced Composites // Composites Part A: Applied Science and Manufacturing, 2016. Vol. 91. P. 85-95.
[5] Wenjie N. et al. Three-dimensional finite element model generation based on CT image for concrete crack // Functional materials. 2016. No. 2. P. 326-330.
[6] Gelb J. et al. Multi-length Scale X-ray Microscopy: A Unique Solution for Digital Rock Physics // EAGE/FESM Joint Regional Conference Petrophysics Meets Geoscience. 2014.
[7] Abdul-Aziz A. et al. Material characterization and geometric segmentation of a composite structure using microfocus X-ray computed tomography image-based finite element modeling. — Cleveland, Ohio : National Aeronautics and Space Administration, Glenn Research Center, 2011.
[8] Avizo for Materials Science Thermo Fisher Scientific Inc. [Электронный ресурс]. URL: https://www.fei.com/software/avizo-for-materials-science
[9] Simpleware Software Solutions Synopsys, Inc. [Электронный ресурс]. URL: http://simpleware.com
[10] Materialise Mimics Overview [Электронный ресурс]. URL: http://www.materialise.com/ en/medical/software/mimics
[11] 3D Slicer Documentation/4.8 [Электронный ресурс]. URL: https://www.slicer.org/ wiki/Documentation/4.8
[12] Yushkevich P. A., Gerig G. ITK-SNAP: an intractive medical image segmentation tool to meet the need for expert-guided segmentation of complex medical images // IEEE pulse. 2017. Vol. 8, No. 4. P. 54-57.
[13] Schindelin J. et al. Fiji: an open-source platform for biological-image analysis // Nature methods. 2012. Vol. 9, No. 7. P. 676.
[14] Chen C. Finite element modeling of trabecular bone from multi-row detector CT imaging. — USA : The University of Iowa, 2014.
[15] Smith M. ABAQUS/Standard User's Manual, Version 6.9. — Providence, RI : Simulia, 2009.
[16] Jovanovic J. D., Jovanovic M. L. Finite element modeling of the vertebra with geometry and material properties retrieved from CT-Scan Data // Facta universitatis-series: Mechanical Engineering. 2010. Vol. 8, No. 1. P. 19-26.
[17] Geuzaine C., Remacle J.-F. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities // International Journal for Numerical Methods Engineering. 2009. Vol. 79, No. 11. P. 1309-1331.
[18] SALOME User's manual [Электронный ресурс]. URL: http://www.salome-platform.org/user-section/documentation/current-release/current-release-gui
[19] Doi A., Koide A. An efficient method of triangulating equi-valued surfaces by using tetrahe-dral cells // IEICE TRANSACTIONS on Information and Systems. 1991. Vol. 74, No. 1. P. 214-224.
[20] Штабель Н. В., Шурина Э. П., Эпов М. И., Михайлова Е. И., Кутищева А. Ю. Анализ электрофизических характеристик трехмерных объектов с микровключениями сложной
формы, возбуждаемых постоянным и переменным электромагнитным полем // Теоретические основы и конструирование численных алгоритмов для решения задач математической физики: XXI Всеросс. конференция и Молодежная школа-конференция, по-свящ. памяти К. И. Бабенко (Дюрсо, 5-11 сентября 2016): Тезисы докладов, 2016. С. 5758.
Авторы:
Элла Петровна Шурина — доктор технических наук, профессор, профессор кафедры Вычислительных технологий, Новосибирский государственный технический университет; главный научный сотрудник, Научно-исследовательская лаборатория моделирования и решения численных задач нефтегазовой отрасли, Новосибирский национальный исследовательский го-сударственный университет; главный научный сотрудник, Институт нефтегазовой геологии и геофизики им. А. А. Трофимука Сибирского отделения Российской академии наук
Дарья Владимировна Добролюбова — аспирант, кафедра Вычислительных технологий, Новосибирский государственный технический университет; младший научный сотрудник, Научно-исследовательская лаборатория моделирования и решения численных задач нефтегазовой отрасли, Новосибирский национальный исследовательский государственный университет; младший научный сотрудник, Институт нефтегазовой геологии и геофизики им. А. А. Трофимука Сибирского отделения Российской академии наук
Екатерина Игоревна Штанько — научный сотрудник, Институт нефтегазовой геологии и геофизики им. А. А. Трофимука Сибирского отделения Российской академии наук; научный сотрудник, Научно-исследовательская лаборатория моделирования и решения численных задач нефтегазовой отрасли, Новосибирский национальный исследовательский государственный университет
Special techniques for objects with complex inner structure based on a CT image sequence
£/~>l • * ** *** 7~~\ 7~~\ * 1 I * ** *** T-i /~>r i 1 * ***
. Shurina, , , D. Dobrolyubova >> , E. Shtanko,
*Trofimuk Institute of Petroleum Geology and Geophysics of Siberian Branch of Russian Academy of Sciences (IPGG SB RAS) 3, Koptug ave., Novosibirsk, Russia, 630090
Abstract. An algorithm for generating a special hierarchal data structure based on CT or MRI image of native and composite objects with complex inner structure is proposed. A CT grayscale image sequence is converted into the hierarchal finite element mesh, consisting of non-convex polyhedrons (polygons in 2D) as macro-scale elements and tetrahedrons (triangles in 2D) as micro-scale elements. Resulting meshes are designed for multiscale finite element modelling and can be used for a wide range of multi physical problems. Key words: simplicial mesh, polyhedral mesh, composite and native media.
References
[1] Zhang Y., Bajaj C., Sohn B. S. (2005) Computer methods in applied mechanics and engineering. 194(48-49):5083-5106.
[2] Cong A. et al. (2005) Geometrical modeling using multiregional marching tetrahedral for bioluminescence tomography. In book: Medical Imaging 2005: Visualization, Image-Guided Procedures, and Display. International Society for Optics and Photonics. Vol. 5744. p. 756764.
[3] Chen C. et al. (2016) Assessment of trabecular bone strength at in vivo ct imaging with space-variant hysteresis and finite element modelling // 2016 IEEE 13th International Symposium on. IEEE Biomedical Imaging (ISBI), , p. 872-875.
[4] Sencu R. M., Yang Z., Wang Y., Withers P., Rau C., Parson A., Soutis C. (2016) Composites Part A: Applied Science and Manufacturing, 91:85-95.
[5] Wenjie N. et al. (2016) Functional materials. 2:326-330.
**Novosibirsk State Technical University 20, Prospekt K. Marksa, Novosibirsk, Russia, 630073
**'Novosibirsk State University 2, Pirogova street, Novosibirsk, Russia, 630090
e-mail: [email protected]. ru
[6] Gelb J. et al. (2014) Multi-length Scale X-ray Microscopy: A Unique Solution for Digital Rock Physics. In book: EAGE/FESM Joint Regional Conference Petrophysics Meets Geosci-ence.
[7] Abdul-Aziz A. et al. (2011) Material characterization and geometric segmentation of a composite structure using microfocus X-ray computed tomography image-based finite element modeling. Cleveland, Ohio, National Aeronautics and Space Administration, Glenn Research Center.
[8] https://www.fei.com/software/avizo-for-materials-science
[9] http://simpleware.com
[ 10] http://www. materialise.com/en/medical/software/mimics
[11] https://www.slicer.org/wiki/Documentation/4.8
[12] Yushkevich P. A., Gerig G. (2017) IEEE pulse. 8(4):54-57.
[13] Schindelin J. et al. (2012) Nature methods. 9(7):676.
[14] Chen C. (2014) Finite element modeling of trabecular bone from multi-row detector CT imaging. USA, The University of Iowa.
[15] Smith M. (2009) ABAQUS/Standard User's Manual, Version 6.9. Providence, RI : Simulia.
[16] Jovanovic J. D., Jovanovic M. L. (2010) Facta universitatis-series: Mechanical Engineering. 8(1):19-26.
[17] Geuzaine C. and Remacle J.-F. (2009) International Journal for Numerical Methods Engineering. 79(11):1309—1331.
[18] http://www.salome-platform.org/user-section/documentation/current-release/current-release-gui
[19] Doi A., Koide A. (1991) IEICE TRANSACTIONS on Information and Systems. 74(1):214-224.
[20] Shtabel' N. V., Shurina E. P., Epov M. I., Mikhaylova Ye. I., Kutishcheva A. Yu. (2016) Analiz elektrofizicheskikh kharakteristik trekhmernykh ob"yektov s mikrovklyucheniyami slozhnoy formy, vozbuzhdayemykh postoyannym i peremennym elektromagnitnym polem. In book: Teoreticheskiye osnovy i konstruirovaniye chislennykh algoritmov dlya resheniya zadach ma-tematicheskoy fiziki: XXI Vseross. konferentsiya i Molodezhnaya shkola-konferentsiya, posvyashch. pamyati K. I. Babenko (Dyurso, 5-11 sentyabrya, 2016). p. 57-58. [In Rus]
