Как стать автором
Обновить
132.32
Smart Engines
Обработка изображений, распознавание в видеопотоке

Алгоритм реконструкции динамических томографических процессов

Уровень сложностиСредний
Время на прочтение17 мин
Количество просмотров513

Привет, Хабр! Мы уже рассказывали про наши успехи в рентгеновской томографии и разрабатываемой нами программе для томографической реконструкции Smart Tomo Engine, ознакомиться с которой можно у нас на сайте Smart Engines. В этом же посте мы хотим поделиться с вами деталями наших исследований в динамической или 4D томографии. Здесь, для исследования объекта, который менялся в процессе проведения измерений, нам пришлось разработать новый алгоритм томографической реконструкции и даже провести гидродинамическое моделирование. Но давайте обо всём по порядку.

Постановка задачи

В одном из предыдущих постов рассказывалось о методе томографической реконструкции объекта в случае неполных (недостоверных) данных на проекциях. Сейчас же мы представляем итеративный метод реконструкции для случая, когда изучаемый объект изменяется во времени и на каждый шаг по времени у нас есть только одно изображение объекта. Этот метод исследования называется времяразрешающей или 4D-томографией и позволяет восстанавливать и визуализировать внутреннюю структуру объектов, в которых происходят динамические процессы. Описанный алгоритм и эксперименты опубликованы в [1].

В качестве динамического процесса будет рассматриваться заполнение пористого объекта жидкостью. Однако все рассуждения справедливы и для других процессов при соблюдении перечисленных ниже условий.

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

Предлагаемый алгоритм реконструкции изображений динамических процессов основан на следующих предположениях:

  1. Мы предполагаем, что в процессе изменения исследуемого объекта количество материи в каждой пространственной точке не уменьшается со временем (процесс является монотонным). При рассмотрении случая заполнения образца жидкостью это означает, что в каждом элементарном объеме (вокселе) объекта количество жидкости не может уменьшаться. Это, например, справедливо при медленном заполнении пористой структуры жидкостью. Если жидкость зашла в пору, то пора уже никогда не будет пустой.

  2. Мы получаем только одно проекционное изображение объекта в каждый момент времени при регистрации динамического процесса (заполнение объекта жидкостью).

  3. Нам заранее известна структура той части объекта, которая заведомо не меняется при изучении объекта. Например, при исследовании заполнения объекта жидкостью мы предполагаем, что структура “пустого” объекта (матрикс) нам известна. Узнать эту структуру можно, например, сделав томографию объекта до начала его заполнения жидкостью. В дальнейшем мы предполагаем, что матрикс объекта будет неизменным, а поры - изменяемой областью, т. е. областью, в которой будет происходить динамический процесс. Это условие необязательно, но данные о неизменяемой структуре объекта улучшают сходимость алгоритма за счет уменьшения количества обрабатываемых пространственных точек.

  4. В описанном ниже эксперименте и моделировании используется приближение параллельной геометрии томографического измерения, т. е. источник находится достаточно далеко, чтобы считать падающие пучки излучения параллельными друг другу. В этом случае мы можем рассматривать 3D объект как набор горизонтальных 2D срезов, и в процессе работы алгоритма каждое такое пространственное сечение обрабатывается независимо от других. Алгоритм будет работать и в других геометриях измерения, однако это потребует больших вычислительных ресурсов. На рис. 1 показана схема получения томографических проекций в параллельном пучке.

Задача разрабатываемого алгоритма: на основе этих данных восстановить динамику процесса (заполнения объекта жидкостью) для каждой временной точки.

Рисунок 1: Схема получения томографических проекций в параллельном пучке [1].
Рисунок 1: Схема получения томографических проекций в параллельном пучке [1].

Описание алгоритма

Перед началом работы алгоритма исследуемый 3D объект разделяется на массив горизонтальных 2D срезов. Для каждого такого среза алгоритм реконструирует динамику изменения среза во времени. После завершения работы алгоритма все реконструированные срезы собираются в 3D-объекты для каждой временной точки. Набор таких 3D-объектов во времени дает желаемое 4D-изображение динамического процесса.

Для описания этапов алгоритма введем следующие обозначения (см. рис. 2): N - количество зарегистрированных проекционных изображений объекта. Под проекцией мы понимаем линеаризованное, значение показаний детектора. В предлагаемом подходе N также равно количеству углов поворота объекта и количеству временных точек в эксперименте. I_n - это проекция объекта в динамическом эксперименте, полученная при угле поворота a_n , в момент времени t_n, где n=1\ldots N.

Рисунок 2: Схема получения экспериментальных данных и их обозначения:   - временные точки;  и  - соответствующие углы поворота и проекционные данные, соответственно;  - экспериментальная синограмма. В верхнем ряду показаны состояния одного среза объекта в ходе эксперимента в указанные моменты времени. Справа показана экспериментальная синограмма . Цифрами 1-4 на первом (левом) изображении среза объекта и на экспериментальной синограмме обозначены поры и их проекционные следы, соответственно. Ниже, под состояниями среза объекта, показаны соответствующие проекционные изображения состояния среза. Наблюдается постепенное, неравномерное заполнение пор жидкостью от момента времени .до момента времени .Чем светлее пора, тем больше в ней жидкости (коэффициент поглощения выше)  [1].
Рисунок 2: Схема получения экспериментальных данных и их обозначения: t_1<t_j<t_i<t_N - временные точки; a_1,a_j,a_i,a_N и I_1,I_j,I_i,I_N - соответствующие углы поворота и проекционные данные, соответственно; S_e - экспериментальная синограмма. В верхнем ряду показаны состояния одного среза объекта в ходе эксперимента в указанные моменты времени. Справа показана экспериментальная синограмма S_e. Цифрами 1-4 на первом (левом) изображении среза объекта и на экспериментальной синограмме обозначены поры и их проекционные следы, соответственно. Ниже, под состояниями среза объекта, показаны соответствующие проекционные изображения состояния среза. Наблюдается постепенное, неравномерное заполнение пор жидкостью от момента времени t_1.до момента времени t_N..Чем светлее пора, тем больше в ней жидкости (коэффициент поглощения выше) [1].

Набор проекций для каждого пространственного среза объекта образует экспериментальную синограмму S_e. Срез имеет толщину в один воксель; воксели в срезе параметризуются набором двумерных векторов \overline{r}. Этот набор делится на две части: неизменяемые, соответствующие твердой структуре (матриксу), и изменяемые, соответствующие заполняемым порам. Набор реконструкций \hat{R} состоит из N элементов R_n, которые являются реконструкциями состояния объекта во все моменты времени t_n. R_0 является реконструкцией начального состояния объекта. Обозначим реконструированные значения коэффициента поглощения (количество жидкости в вокселе) в данном вокселе \overline{r} в момент времени t_n через R_n(\overline{r}) и R_0(\overline{r}) для начального состояния.

Исходя из нашего предположения о процессе (количество материала в любой точке не уменьшается), условие R_n(\overline{r})\geq R_{n-1}(\overline{r}) должно выполняться для любого n и любого \overline{r}. Определим L_2-норму для единичной реконструкции \|R_n\|_2:

\|R_n\|_2=\left( \sum_{ch.a}R_n^2(\overline{r}) \right)^{1/2}

и всего множества реконструкций \| \hat{R} \|_2:

\|\hat{R} \|_2=\left(\frac{1}{N} \sum_{ch.a}\sum_{n=1}^{N}R_n^2(\overline{r}) \right)^{1/2},

где ch.a обозначает изменяемую область объекта.

Последовательность шагов алгоритма (см рис. 3):

Шаг 0. Инициализация данных.

Создаём N одинаковых синограмм S_n\leftarrow S_e(n=1\ldots N), соответствующие начальному состоянию объекта (без жидкости). В ходе дальнейших шагов алгоритма синограммы S_n (n=1\ldots N) будут меняться.

Шаг 1. Использование экспериментальных данных.

Из эксперимента известно, что n-й ряд синограмм S_e соответствует углу проекции a_n, a n-й ряд в каждой синограмме S_n должен быть равен экспериментальным данным, соответствующим углу проекции a_n. Для каждого n=1\dots N мы заменяем n-ю строку в синограмме S_n на n-ю строку в экспериментальной синограмме S_e. Остальные N-1 оставшихся строк остаются неизменными.

Шаг 2а. Реконструкция.

Вычисляем набор реконструкций \hat{R} с помощью алгебраического метода SIRT для каждой обновленной синограммы S_n.

Шаг 2b. Использование начальных условий.

Для каждого \overline{r} неизменяемой области заменяем значение R_n(\overline{r}) на R_0(\overline{r}), если известна реконструкция начального состояния объекта.

Шаг 3. Использование предварительной (априорной) информации о процессе.

Для каждого n и каждой точки \overline{r} изменяемой области объекта, проверяем условие неубывания количества материала. Если оно не выполняется, скорректируем соответствующие значения для каждого R_n, и сохраним их в переменной \hat{R}_n^{iter}, которая содержит R_n^{iter} и будет использоваться на следующем шаге алгоритма:

R_n^{iter}(\overline{r})=\min{\left(R_n(\overline{r}), R_{n+1}(\overline{r}) \right)}.

Шаг 4. Повтор итерации или выход. Рассчитываем синограммы S_n по каждому из \hat{R}^{iter} и вернуться к Шагу 1. Если результат шагов 1-3 приводит к недостаточному изменению набора реконструкций (изменение нормы L_2 на величину меньше, чем 10^{-5}), мы выходим из алгоритма, т.к. от итерации к итерации уже ничего не меняется. Незначительное изменение набора реконструкций означает, что замены рядов в синограммах (Шаг 1) и корректировки за счет информации о процессе (Шаг 2) не оказывают заметного влияния на результат. Таким образом, рассчитанные синограммы S_n близки к экспериментальной синограмме S_e и являются описанием наблюдаемого процесса.

Рисунок 3: Иллюстрация этапов алгоритма. На шаге 2 для точек 1 и 3 не выполняется условие неубывающего коэффициента поглощения, что было исправлено на шаге 3 [1].
Рисунок 3: Иллюстрация этапов алгоритма. На шаге 2 для точек 1 и 3 не выполняется условие неубывающего коэффициента поглощения, что было исправлено на шаге 3 [1].

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

Вход: экспериментальная двумерная синограмма S_e, каждая n\left( n=1\ldots N \right) строка S_e^n соответствует экспериментальным данным под углом a_n (в момент времени t_n), R(\overline{\mathbf{r}}) - двумерная реконструкция пустой пористой структуры, \overline{\mathbb{r}}_{\text{unchangeable}} - область на реконструкции, где значения реконструкции не должны меняться из-за предварительной информации об образце. \varepsilon - точность.
Предварительная информация об объекте: Концентрация жидкости не может уменьшаться \mathbf{R}_n(\overline{\mathbf{r}})\leq \mathbf{R}_{n+1}(\overline{\mathbf{r}}), для n=1\ldots N-1
Выходные данные: \mathbf{R}_n(\overline{\mathbf{r}}) - 2D-реконструкции объекта в каждый момент времени \mathbf{t}_nt_n для n=1…N
Инициализация:
Создаём N синограмм для итерационной процедуры S_n \leftarrow S_e для n=1\ldots N
Создаём N реконструкций для расчета первой итерации набора реконструкций R^n=\text{SIRT}(S_n)
Шаг 1: Использование экспериментальных данных
do:
Обновление n-й строки в каждой синограмме S_n для приведения ее в соответствие с экспериментальными данными

for n in 1...N:

S_n^n \leftarrow S_e^n

Шаг 2а: Реконструкция.

for n in 1...N:

R_n^{prev}\leftarrow R_n
R_n\leftarrow \text{SIRT}(S_n)

Шаг 2b: Использование начальных условий для каждой \overline{r} неизменной области \overline{\textbf{r}}_{\text{unchangeable}}.

for n in 1...N:

R_n(\overline{r}) \leftarrow R_0(\overline{r}), \forall \overline{r} \in \overline{\mathbf{r}}_{\text{unchangeable}}
Шаг 3: Использование предварительной (априорной) информации о процессе.
for n in 1...N-1:

R_n^{iter}(\overline{r}) \leftarrow \min(R_n(\overline{r}),R_{n+1}(\overline{r}))

for n in 1...N:

R_n(\overline{r}) \leftarrow R_n^{iter}(\overline{r})

S_n = CalculateSinogram(R_n).

Шаг 4. Повтор итерации или выход.

while \|R_n-R_n^{prev} \|_2>\varepsilon \: \forall n=1 \ldots N

return R_n \forall n=1\ldots N

Примечание 1. На этапе 2а могут использоваться и другие методы томографической реконструкции, например, FBP. Выбор конкретного метода зависит от параметров эксперимента (отношение сигнал/шум, общее количество проекций и т.д.). Алгебраический метод SIRT лучше работает при низком соотношении сигнал/шум, в то время как FBP является быстрым при достаточно хороших данных.
Замечание 2. Если начальное состояние объекта неизвестно, то шаг 2b пропускается.
Замечание 3. Наши расчеты показали, что алгоритм либо сходится за разумное число итераций, либо не сходится вообще. Таким образом, критерий выхода можно заменить достижением достаточно большого числа итераций (например, 10 000) и контролем изменения нормы L_2 на последней итерации. Ниже будет показана эквивалентность этих условий.

Проверка работы алгоритма

Как говорилось выше, одной из областей применения предложенного алгоритма является изучение динамики движения жидкости в пористых средах. Для проверки алгоритма было проведено два эксперимента по заполнению пористой структуры, в случае, когда известна информация о структуре матрикса.

Во-первых, был использован модельный объект с искусственно сгенерированными данными. В этом случае структура матрикса была известна из модели.

Второй эксперимент связан с физическим объектом. В нём мы наблюдаем подъем масла в вертикальном капилляре. Простая геометрия позволяет сравнить результаты томографии с теоретической моделью из гидродинамики, а также с визуальными наблюдениями.

В следующих подразделах приводятся подробности этих экспериментов.

Виртуальный эксперимент с модельным объектом

Для виртуального эксперимента мы создаем объект: двумерную модель пористой структуры (один горизонтальный срез трехмерного объекта; см. рис. 4). Мы будем называть этот объект моделью. Он состоит из 1024 (32x32x1) вокселей. Часть из них образует неизменную область, представляющую матрицу, остальные - поры. Поры расположены в центре объекта, чтобы избежать артефактов томографической реконструкции, и образуют изменяемую область. Они пусты при t=0. На рис. 4 матрикс и поры показаны белым и черным цветом соответственно.

Рисунок 4: Модель одного среза пористого объекта. Поры (пустоты), которые могут быть заполнены жидкостью, показаны черным цветом [1].
Рисунок 4: Модель одного среза пористого объекта. Поры (пустоты), которые могут быть заполнены жидкостью, показаны черным цветом [1].

Далее мы моделируем заполнение пространства жидкостью. Для этого генерируем состояние модельного объекта в 100 временных точках с равным интервалом между ними: в любой момент времени t_n для каждого вокселя \overline{r} в изменяемой области, мы устанавливаем число c_n(\overline{r}) от 0 до 1, представляющее собой коэффициент поглощения (который соответствует количеству жидкости). Для фиксированного \overline{r} последовательность c_n(\overline{r}) монотонно растет с увеличением n. На рис. 5 показан модельный объект для нескольких временных точек; цвет пор меняется от черного к белому, пока жидкость заполняет поры, а коэффициент поглощения растет. Мы вычисляем проекции модельного объекта и используем их в качестве входной информации для нашего алгоритма.

Рисунок 5: Модель распределения жидкости в вокселях пористого объекта в разные моменты времени: (а) точка 0, (б) точка 50 и (в) точка 100. Пустоты показаны тёмным цветом, по мере наполнения жидкостью они светлеют [1].
Рисунок 5: Модель распределения жидкости в вокселях пористого объекта в разные моменты времени: (а) точка 0, (б) точка 50 и (в) точка 100. Пустоты показаны тёмным цветом, по мере наполнения жидкостью они светлеют [1].

Для оценки влияния количества проекций на качество реконструкции были создали наборы из N проекций I_N для N=100,50,25 с учетом каждой первой, второй и четвертой временной точки. Углы проекций равномерно распределены по окружности для каждого случая. Сравнение результатов для различных значений N дает оценки влияния количества проекций.

Анализ данных

Сначала был проанализирован наиболее полный набор данных для N=100, а затем та же процедура была применена для худшего временного разрешения. В этом подразделе результаты соответствуют большому значению N, если иное не указано явно.
В качестве метода реконструкции (шаг 2 описанного алгоритма) использовался алгебраический метод SIRT из пакета программ ASTRA toolbox. Выполнялось 10 000 итераций нашего алгоритма. На рис. 6 показаны промежуточные реконструированные изображения для середины всего временного интервала (50-й момент времени). К 1000-й итерации алгоритма полученное изображение находится в хорошем согласии с истиной (среднеквадратичные отклонения восстановленных значений от заданных в модели для итераций 1, 200 и 1000 составили 0,0322, 0,0095 и 0,0049 соответственно). На рис. 7 показана зависимость нормы реконструкции от номера итерации (рис. 7а) и изменение разницы между нормами последовательных итераций (рис. 7б). В данном случае общее число итераций 1000 представляется разумным компромиссом между точностью реконструкции и временем работы алгоритма.

Рисунок 6: Результаты итеративной реконструкции заполнения вокселя жидкостью в момент времени 50 : (a) 1 итерация, (b) 200 итераций, (c) 1000 итераций, и (d) значение модели (заложенное в модель значение) [1].
Рисунок 6: Результаты итеративной реконструкции заполнения вокселя жидкостью в момент времени 50 : (a) 1 итерация, (b) 200 итераций, (c) 1000 итераций, и (d) значение модели (заложенное в модель значение) [1].
Рисунок 7: (a) Зависимость  -нормы набора реконструкций от номера итерации (красная пунктирная линия обозначает  -норму объекта моделирования); (b) изменение разницы между нормами последовательных итераций;  [1].
Рисунок 7: (a) Зависимость L_2-нормы набора реконструкций от номера итерации (красная пунктирная линия обозначает L_2-норму объекта моделирования); (b) изменение разницы между нормами последовательных итераций; N=100 [1].

На рис. 8 показаны результаты реконструкции для различных значений N (шагов по времени), представляющие зависимости коэффициента поглощения от времени для трех выбранных точек (обведены кружками на вставке в верхнем левом углу рис. 8). Отклонения реконструированных данных от истины уменьшаются при увеличении количества шагов по времени N. Наилучшая реконструкция наблюдается для точки, выделенной фиолетовым цветом (на вставке в верхнем левом углу рис. 8). Немонотонный результат для зависимости коэффициента поглощения от времени на промежуточных итерациях является препятствием для распространения коррекции между итерациями. Можно заметить, что реконструированные зависимости не являются монотонными и даже показывают значения поглощения меньше нуля и больше единицы, что не соответствует заложенной модели объекта. Это может быть связано с тем, что реализованная в алгоритме коррекция (шаг 3) рассматривает все пары реконструкций для последовательных временных шагов независимо, но при этом не производит глобального упорядочивания. В результате малые значения, ошибочно полученные для больших n, не влияют на результат для меньших n. Несмотря на это, разработанный алгоритм на тестовых данных даёт ошибку меньше 0.2, при поглощении в интервале от 0 до 1, тогда как реконструкция традиционными методами на этом объекте даёт ошибку до 0.9.

Рисунок 8: Сравнение томографической реконструкции модельного объекта (см. изображение в левом верхнем углу) для различного количества проекций (). Сплошные линии - профили грунтовых истин в пикселях, отмеченных на модели соответствующим цветом. Пунктирные линии - реконструкция при разном количестве проекций () [1].
Рисунок 8: Сравнение томографической реконструкции модельного объекта (см. изображение в левом верхнем углу) для различного количества проекций (N). Сплошные линии - профили грунтовых истин в пикселях, отмеченных на модели соответствующим цветом. Пунктирные линии - реконструкция при разном количестве проекций (N=100,50,25) [1].

Проверка работы алгоритма на реальных экспериментальных данных

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

Экспериментальная установка

Эксперимент проводился на описанном выше лабораторном томографе “ТОМАС” в ФНИЦ “Кристаллография и фотоника” РАН. В качестве источника использовалась рентгеновская трубка с молибденовым анодом. Для выделения характеристической линии использовался кристалл- монохроматор из пиролитического графита. Энергия зондирующего излучения составляла 17,5 кэВ. Расстояние от источника до образца составляло 1250 мм, а от образца до детектора - 15 мм. Соотношение этих расстояний позволяло работать в приближении параллельной томографической схемы (см. рис. 1). Регистрация излучения осуществлялась с помощью детектора XIMEA xiRAY11 с размером пикселя 9 × 9 мкм и полем зрения 36 × 24 мм. Характерное время для полного томографического эксперимента составляет 0,5-2 часа, в зависимости от уровня поглощения исследуемого объекта.

Рисунок 9: Исследуемый капилляр, установленный в рентгеновской микротомографической установке “ТОМАС” (слева) и держатель капилляра с емкостью для рабочей жидкости (справа) [1].
Рисунок 9: Исследуемый капилляр, установленный в рентгеновской микротомографической установке “ТОМАС” (слева) и держатель капилляра с емкостью для рабочей жидкости (справа) [1].

Характерное время эксперимента и размер поля зрения детектора определяют диапазоны размеров капилляра и вязкость рабочей жидкости. Для того, чтобы время проведения динамического эксперимента было сравнимо со временем традиционной томографии мы подобрали параметры капилляра и масла таким образом, чтобы время заполнения капилляра составило 30 мин: использовали стеклянный капилляр высотой 33,4 мм с внутренним каналом диаметром 0,2 мм и, заполнявшийся силиконовом маслом (полиметилсилоксан) PMS-1000 вязкостью 1000 сСт, поверхностным натяжением 0,019 Н/м, и контактным углом со стеклом 0 градусов.

Капилляр был установлен на специальный держатель с емкостью для рабочей жидкости (см. рис. 9). Протокол проведения эксперимента состоял из двух частей. Сначала проводилось томографическое исследование пустого капилляра (без жидкости). Затем в держатель капилляра залили масло, и во второй части эксперимента мы повторили томографическое исследование во время подъема жидкости в капилляре. В каждом томографическом сканировании мы записали 444 проекции с шагом поворота образца 0,45 градуса. Время экспозиции каждой проекции составляло 1,5 секунды, а пауза между проекциями, связанная с поворотом образца и инициализацией детектора, составляла 3,1 секунды. Общая продолжительность одного сканирования составила 34 минуты. Таким образом, для каждого угла поворота образца мы получили проекции как пустого (рис. 10а), так и частично заполненного капилляра (рис. 10б). Для каждой пары изображений мы вычисляли логарифм отношения нормализованного сигнала в каждой точке (разность изображений), как показано на рисунке рис. 10c. Полученные разности использовались в качестве исходных данных для работы алгоритма. Вычисления проводились с помощью реализации алгоритма на языке Python на процессоре Intel Core i7- 5930K, 64Gb RAM, Nvidia GTX 980Ti. В качестве метода реконструкции использовался алгебраический метод SIRT. Для каждого среза томограммы было выполнено 100 итераций предложенного алгоритма, так как дальнейшие итерации не улучшают качество изображения из-за шума в экспериментальных данных. Расчет каждого среза занял 20 минут.

Рисунок 10: Теневые проекции пустых (a) и частично заполненных (b) капилляров и разница между этими изображениями (c). Заполненная область соответствует темным точкам на панели (c) [1].
Рисунок 10: Теневые проекции пустых (a) и частично заполненных (b) капилляров и разница между этими изображениями (c). Заполненная область соответствует темным точкам на панели (c) [1].

После обработки данных с помощью предложенного алгоритма мы получили реконструированные изображения горизонтальных срезов объекта во всех временных точках (рис. 11). Последующее соединение срезов в единый объект позволило получить 4D изображение капиллярного подъема.

Рисунок 11: Изображения столба жидкости в капилляре, полученные из проекционных данных и реконструированные алгоритмом срезов объекта на высоте 2,4 мм (красный), 9 мм (зеленый) и 16 мм (синий) для временных точек (a) 380 с, (b) 1120 с и (c) 1500 с. Заполненная область соответствует темным точкам [1].
Рисунок 11: Изображения столба жидкости в капилляре, полученные из проекционных данных и реконструированные алгоритмом срезов объекта на высоте 2,4 мм (красный), 9 мм (зеленый) и 16 мм (синий) для временных точек (a) 380 с, (b) 1120 с и (c) 1500 с. Заполненная область соответствует темным точкам [1].

Сравнение с теорией и моделированием

Чтобы оценить точность реконструкции динамического процесса по представленному алгоритму, необходимо сравнить реконструированное 4D-изображение с экспериментальными данными. В случае простого объекта (например, капилляра) независимые результаты могут быть получены из теневой проекции или расчетов.

Различия между теневыми проекциями пустого и частично заполненного капилляра обеспечивают четкие границы жидкости и дают точность определения высоты столба жидкости порядка нескольких пикселей (~10 мкм). Теоретическая зависимость высоты столба жидкости ℎ над уровнем жидкости вдали от капилляра от времени ? следует из баланса между подъемной силой поверхностного натяжения, тянущей вниз силой тяжести и замедляющим вязким трением. Это выражается формулой:

\frac{t-t_0}{T}=-\left( 1+ \frac{h_0}{H} \right) \ln \left( 1-\frac{h}{H} \right)-\frac{h+h_0}{H},

где H=\frac{2 \sigma \cos\theta}{\rho g R}, T=\frac{16 \mu \sigma \cos \theta}{\rho^2 g^2 R^3}- опорная высота и время. Они определяются коэффициентом поверхностного натяжения и контактным углом \sigma, \theta, радиусом капилляра R, плотностью и вязкостью жидкости \rho и \mu, а также глубиной погружения нижнего края капилляра h_0; t_0 свободный параметр, позволяющий учесть начальный интервал.

При томографической реконструкции в параллельной схеме мы считаем, что горизонтальные слои могут быть восстановлены независимо друг от друга. Будем считать, что слой (воксель) на высоте h заполненным, когда средний коэффициент поглощения вокселя внутри капилляра достигает заданного значения.

Для сравнения томографической реконструкции с теорией и моделированием необходимо определиться, как определять по томографической реконструкции уровень поднятия жидкости. Уровень заполнения пикселя в теории может составлять от 0% до 100%. Однако при томографических исследованиях всегда присутствуют шумы, из-за которых крайние значения диапазона не достигаются. Нами были выбраны 3 критерия, при достижении которых можно считать, что воксель наполнился жидкостью - 85%, 90% и 95%. Сравнения скорости поднятия жидкости по теории, численному моделированию и томографии сравнивались именно по этим трём критериям заполнения.

Рисунок 12: Сравнение результатов визуального определения высоты подъема капилляров с расчетными значениями и с данными томографической реконструкции [1].
Рисунок 12: Сравнение результатов визуального определения высоты подъема капилляров с расчетными значениями и с данными томографической реконструкции [1].

На рис. 12 показано сравнение экспериментальных данных визуальных наблюдений, расчета с использованием безынерционного режима подъема капилляров и реконструкции динамики процесса по предложенному алгоритму. Для последнего используются три описанных выше критерия заполнения.

Время выполнения критерия заполнения менее чувствительно к процентному содержанию на нижнем слое, поскольку жидкость заполняет нижние слои достаточно быстро. Для верхних же слоев мы имеем небольшое количество проекций с заполненным слоем, и отношение сигнал/шум хуже, чем для нижних слоев.

Если заполняемый объект имеет более сложную форму (например, пористый керамический фильтр или керн горной породы), визуально определить количественные характеристики заполнения пористой среды жидкостью по данным проекции может быть сложно. Приближенные аналитические методы расчета потока также будут неприменимы. Для сравнения можно использовать, например, зависимость массы жидкости внутри образца от времени и оценить положение границы раздела, рассчитанное по томографической реконструкции.

Обсуждение полученного результата

Как видно из представленных результатов, предложенный алгоритм томографической реконструкции позволяет восстанавливать структуру объектов, изменяющихся во времени по монотонному закону. При этом временное разрешение потенциально ограничено только временем экспозиции одной проекции. Предложенный метод реконструкции приводит к уменьшению общего времени эксперимента по сравнению с проведением полного томографического измерения в каждый момент времени и, соответственно, к снижению накопленной дозы облучения исследуемого объекта. Это также позволяет применять данный метод для изучения динамических процессов на лабораторных источниках рентгеновского излучения, где эксперименты по 4D-томографии с реконструкцией квазистатических состояний по большому набору проекций сильно ограничены в своих возможностях из-за относительно большого времени экспозиции каждой проекции. В то же время применение предложенного алгоритма к синхротронным источникам позволит изучать чрезвычайно быстрые и непериодические процессы.

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

Таким образом продемонстрировано, что представленный алгоритм математической реконструкции 4D томографии позволяет реконструировать динамические процессы при наличии только одного проекционного изображения на момент времени, при этом процесс изменения объекта является неубывающим по интенсивности вокселей реконструированных изображений. Мы доказали реализуемость предложенного алгоритма, используя данные двух экспериментов: виртуального со сложной структурой пористого объекта, когда входные данные генерировались искусственно и априорно была известна истина, и реального с данными, полученными при рентгеновской томографии простого объекта. Дизайн последнего эксперимента позволил проводить независимый мониторинг процесса заполнения по теневым снимкам и простому моделированию. Экспериментальная проверка путем изучения процесса заполнения капилляра жидкостью продемонстрировала возможность практического применения алгоритма. Примером применения предложенного метода является исследование потоков жидкости в пористых объектах. Поскольку визуального анализа восстановленных изображений пористых объектов уже недостаточно, были также разработаны методы автоматической бинаризации томографических изображений пористых структур, используемых в качестве маски в нашем алгоритме [2]. Мы также понимаем, что очень важен компромисс между скоростью вращения и соотношением сигнал/шум при мониторинге динамических процессов. Наши эксперименты показали, что мониторинг быстрых процессов, когда скорость процесса не позволяет измерить более одной томографической проекции, потенциально возможен. Другим преимуществом такого подхода является возможность его применения в медицинских исследованиях из-за чрезвычайно низких доз облучения при КТ. Сверхмалое время экспозиции приводит к крайне плохому соотношению сигнал/шум, однако в настоящее время ведется разработка новых подходов для работы с томографическими данными, собранными в условиях сверхмалых доз. В настоящее время мы работаем над усовершенствованием предложенного алгоритма 4D томографии. Мы планируем адаптировать его к таким динамическим процессам, где изменения в исследуемых объектах могут быть немонотонными.

Список литературы
  1. Buzmakov A., Krivonosov Y., Grigoriev M., Mogilevskiy E., Chukalina M., Nikolaev D., Asadchikov V. Iterative algorithm for 4d tomography reconstruction using a single projection per time step // IEEE Access. — 2022. — Vol. 10. — P. 46963–46974.

  2. Chukalina M. V., Khafizov A. V., Kokhan V. V., Buzmakov A. V., Senin R. A., Uvarov V. I., Grigoriev M. V. Algorithm for post-processing of tomography images to calculate the dimension-geometric features of porous structures // Computer Optics. — 2021. — Vol. 45, no. 1. — P. 110–121.

Теги:
Хабы:
Всего голосов 2: ↑2 и ↓0+2
Комментарии0

Публикации

Информация

Сайт
smartengines.ru
Дата регистрации
Дата основания
Численность
51–100 человек
Местоположение
Россия