Моделирование артефактов и методы их фильтрации в рентгеновской компьютерной томографии Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 519.6

МОДЕЛИРОВАНИЕ АРТЕФАКТОВ И МЕТОДЫ ИХ ФИЛЬТРАЦИИ В РЕНТГЕНОВСКОЙ КОМПЬЮТЕРНОЙ ТОМОГРАФИИ

В.В. Ласьков, Е.Н. Симонов

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

Ключевые слова: компьютерная томография, артефакты, метод фильтрации.

Введение

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

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

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

Из-за особенностей процесса реконструкции изображения по проекциям и его сложности каждое измерение в проекции усредняется («размазывается») в виде прямой линии на изображении, и ошибка в измерении проекции не локализуется, как в случае традиционной рентгенологии.

Артефакты могут сильно искажать томографическое изображение и приводить в итоге к постановке ошибочного диагноза.

1. Основные типы артефактов

Согласно монографии [1], в общем случае артефакты изображений КТ можно классифицировать по четырем основным категориям: строки, затенения, кольца и полосы. На рис. 1 представлены томографические изображения, содержащие артефакты различных типов.

Строки часто появляются как яркие прямые линии (не обязательно параллельные) на изображении. Они могут быть как светлыми, так и темными. Во многих случаях светлые и темные строки появляются парами из-за структуры «ядра» фильтрации алгоритма реконструкции.

Артефакты затенения появляются возле областей объекта с высоким контрастом. Например, они появляются либо возле воздушных карманов, либо в области мягких тканей, около которых имеются плотные структуры. Артефакты затенения создают непредсказуемые скачки чисел КТ (чисел Хаунсфилда) на изображении и могут привести к неправильной диагностике, если они не были ранее идентифицированы.

Рис. 1. Различные типы артефактов томографических изображений

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

Было проведено исследование, целью которого являлось изучение возможности применения «классических» методов фильтрации [2] для подавления артефактов томографического изображения: линейная пространственная фильтрация с различным размером маски фильтрации (метод скользящего среднего, метод взвешенного среднего), нелинейная пространственная фильтрация с различным размером маски фильтрации (медианная, метод первой производной, метод второй производной, эквализация, комбинированные фильтры), адаптивная пространственная фильтрация (адаптивные локальные фильтры, адаптивные медианные фильтры, рекурсивные фильтры, адаптивные оптимальные фильтры), частотная фильтрация (фурье-фильтрация, вейвлет-фильтрация).

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

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

Для проведения моделирования был разработан универсальный реконструктор томографического изображения методом обратного проецирования с фильтрацией сверткой на произвольное количество ракурсов облучения (проекций) и единичных детекторов (отсчетов, измерений в проекции), для любой конфигурации сечения объекта [3, 4]. Реконструктор позволяет решать как прямую задачу томографии - определение проекций, так и обратную - восстановление изображения по проекциям.

Прямая задача томографии - вычисление интенсивности каждого луча, прошедшего через объект. Лучи распространяются в объекте вдоль прямой линии l, определяемой уравнением:

х cos 0 + y sin 0- s = 0, где s - расстояние от соответствующего луча до начала координат.

Интенсивность луча на выходе из объекта равна интегралу от искомой функции u (х, y) вдоль траектории луча l :

да да

p ( s, 0)= I I u ( х, y )5( х cos 0 + y sin 0- s ) dxdy,

—да —да

где 5(*) - дельта-функция.

p (s, 0) называется радоновским образом, а преобразование - преобразованием Радона.

Проекционные данные для «карандашного» рентгеновского луча на основе закона Бугера -Ламберта - Бера рентгеновского излучения:

(I (^, ер

p ( 5,6) = - ln

(1)

где I (s, 0) и /0 - интенсивность рентгеновского излучения, соответственно, после объекта исследования и до объекта.

Проекционные данные p(s,0) представляются обычно в виде цифровой матрицы М х N

(рис. 2), где М - количество ракурсов (углов) облучения объекта исследования, обычно М = 600^1200, N - количество единичных детекторов в системе детектирования (отсчетов), обычно N = 512^1024; Д0 и As - соответственно, размер единичного ракурса и детектора, j и i -соответственно, их индексы.

Обратная задача томографии - восстановление u ( x, y) при известных проекционных данных

Р(s0) .

Проекция p(s,0) функции двух переменных u (x,y) для каждого угла 0 представляет одномерную функцию. Ее можно преобразовать в двумерную, зафиксировав угол 0 и выполнив обратное проецирование по всей плоскости (x, y) в соответствии с выражением

p0( x, y) = p (x cos 0 + y sin 0,0).

Далее складываются все обратные проекции p0 (x, y) для 0 < 0 < n . В результате получается суммарное изображение, которое используется в качестве оценки функции плотности u (x, y) . Суммарное изображение определяется соотношением

u ( x, y) = j p ( x cos 0 + y sin 0,0) d 0,

0

n да

u (x, y) = j J g ( s'- s ) p ( s, 0) dsd 0,

или

0 —да

да

g(s )=4П2 J \w\K (w)

(2)

где K(w) - функция окна, s' = xcos 0 + y sin 0.

Томографическое изображение u (x, y) представляется обычно в виде цифровой матрицы

Мх х My (рис. 3) значений серо-белого (256 градаций), где Мх, My - количество пикселей изображения, соответственно, по оси Х и Y, обычно матрица Мх х My равна 512 х 512 или 1024 х 1024 пикселей.

Рис. 2. Матрица проекционных данных p( s, 0)

Рис. 3. Матрица томографического изображения

и(я, j)

0

Учитывая, что обратная задача томографии относится к классу некорректных задач, для ее регуляризации используют функции окна К(ц>), для которых получают различные сворачивающие функции g (5), например, функции Ромачандрана - Лакшминараянана, Шеппа - Логана,

Ханна, Хемминга, Римана и т. д.

Для прямоугольного окна:

[1, при 1^ ^ wпр,

К>М = [ р

10, при \w\ > wпр.

Для этого случая сворачивающая функция примет вид:

1 7 1^1 ■

g(s) = ^ [ \w\K0(w)elWSdw = — [ \w\elW5dw. 4 л2 4л2

Такой интеграл уже был рассмотрен выше, он равен следующему выражению:

^пр 5

1 “р W 1 2 (^„5^

g(5) =----2 I W ^^5^ = —5)-----------------------sin

2л о 2л 5

2 2 л5

2

(3)

Эту функцию называют сворачивающей функцией Ромачандрана - Лакшминараянана.

л

ДТ

при k = 0 :

Ее можно представить в дискретном виде. Если принять 5 = кД£ , wпр = —, тогда

g (0) = lim

4 ’ к^0

1

(

л

-Нт

2Д12 k^о

2л2 Д£2 к

sin (лк )

лк

sin (лк ) -

1

_2,_2*й2^п I 2

л2к 2ДГ

лк

Л

1

-Нт

4Д£2 к^о

С . (лк smI —

1 2

лк

2

2

2Д£2 4Д£2 4Д£2

при к четном:

sin (лк) = sin | — | = 0,

лк

g (к М ) = -

1

2Д£ 2лк при к нечетном: sin (лк) = 0,

-2 )=

(лк ) -

л2к 2Д£2

^т2 I — I = 0;

g (к Д1 ) =

1

2Д1 лк

-sin (лк ) -

1

л2к 2Д£2

Sin

лк

2

л2к2Д£2 '

(4)

Реконструктор - программа, разработанная с применением формул (1)-(4), для моделирования проекционных данных и реконструкции по ним изображения (обратная задача томографии), а также для моделирования изображения и получения по нему проекционных данных (прямая задача).

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

На рис. 4 представлены «идеальные» проекционные данные и «идеальное» томографическое изображение, восстановленное по этим проекционным данным [5].

1

1

1

1

2

1

о □ О О

а) б) в) г)

Рис. 4: а - «идеальные» проекционные данные; б - «идеальное» томографическое изображение, востанов-ленное по этим проекционным данным; в - трехмерная визуализация «идеальных» проекционных данных; г - трехмерная визуализация «идеального» томографического изображения, восстановленного по этим проекционным данным

1.1. Артефакты в виде полос

Причина появления артефактов в виде полос - однократные сбои единичного канала детекторной системы во время измерения проекций.

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

а) б) в) г)

Рис. 5: а - матрица проекционных данных с единичными сбоями канала; б - томографическое изображение, восстановленное по этим проекционным данным; в - трехмерная визуализация матрицы проекционных данных; г - трехмерная визуализация восстановленного изображения

Семь сбоев отобразились в семь прямых линий на реконструируемом изображении.

1.2. Артефакты в виде колец

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

На рис. 6 показаны проекционные данные с неисправностью единичного канала и реконструкция томографического изображения по этим проекционным данным.

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

Проекционные данные на рис. 7, а, б содержат неисправности каналов детекторной системы, распространяющиеся на все ракурсы облучения. Проекционные данные на рис. 7, а содержат положительные постоянные ошибки, проекционные данные на рис. 7, б содержат отрицательные постоянные ошибки. Реконструкции представлены на рис. 7, д, е.

д) е) ж) з)

Рис. 6: а, б - матрицы проекционных данных с неисправным каналом; в, г - трехмерные визуализации матриц проекционных данных; д, е - томографические изображения, восстановленные по этим проекционным данным; ж, з - трехмерные визуализации восстановленных изображений

д) е) ж) з)

Рис. 7: а, б - матрицы проекционных данных с неисправностями каналов; в, г - трехмерные визуализации матриц проекционных данных; д, е - томографические изображения, восстановленные по этим проекционным данным; ж, з - трехмерные визуализации восстановленных изображений

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

Также следует отметить, что амплитуда колец (значение яркости) уменьшается по мере отдаления от центральных каналов.

2. Предлагаемый алгоритм

Был разработан метод фильтрации проекционных данных, содержащий следующие этапы:

1. Сжатие

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

Значение К выбирается таким образом, чтобы расстояние между двумя сжатыми проекциями соответствовало углу разворота примерно 3° или 120 проекциям на один полный разворот.

2. Высокочастотная фильтрация

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

l=-L/ 2

3. Низкочастотная фильтрация 1

Применение низкочастотной фильтрации в проекционном направлении позволяет уменьшить амплитуду коротких структур и сохранить амплитуду длинных колец.

L 2

Lowpass1 (j,i)= ^ Highpass(j-1,/')• lp(l). (7)

4. Блок дифференцирования

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

Diff (j,i) = Lowpass1( j,i) - Lowpass1 (j -1,i ). (8)

5. Блок принятия решений

Решение, принадлежит ли точка кольцу, принимается по двум пороговым значениям.

Если амплитуда значения дифференцированного сигнала Diff (j, i) не превышает порога градиента S 4 < 0, предполагается, что точка данных входит в кольцо, амплитуда которого была рассчитана по проекционной матрице Lowpass1(y,k) . В ином случае это значение устанавливается равным 0.

Если амплитуда значения проекционной матрицы Lowpass1 ( y, k) превышает амплитудное

пороговое значение S3 > 0 , то значение ограничивается ±S3 . Этот шаг предпринимается для уменьшения эффекта неточных коррекций.

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

Dec( j,i) = Lowpass1 ( j,i), если Diff ( j,i)< S4;

Пороговое значение 53 не зависит от номера канала, а значение £4 зависит: для области центральных каналов применяется пороговое значение 54/, для периферийных каналов - пороговое значение 54а. Устанавливается равенство 54/ = 2*54а, так как условия обнаружения для внешних каналов жестче, чем в области центральных каналов. Это связано с тем, что внутренняя область проекций обычно содержит данные структур более высоких частот.

(5)

L 2

Highpass (j,i) = £ Comp (j, i -1 )• hp (l).

(6)

Dec (j, i) = 0, Dec (j, i) = S 3, Dec (j, i) = -S 3,

если Diff ( j, i )> S 4; если Lowpass1( j,i)> S3; если Lowpass1(j,i)<-S3.

(9)

6. Низкочастотная фильтрация 2

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

L/ 2

Lowpass2 (j, i )= ^ Dec (j -1, i )• lp2 (l). (10)

l=-L/ 2

7. Декомпрессия

Если в начале метода проводилось сжатие, то следует провести декомпрессию как обратную процедуру по отношению к сжатию.

8. Коррекция

Данные, содержащиеся в Lowpass2 (j,i) , представляют корректирующую проекционную матрицу, которая на данном шаге вычитается из исходной проекционной матрицы. В этом случае каждая корректирующая проекция Lowpass2 (jcomp,*) применяется к K входных проекций:

Summ [( j -1) * K +1,i] = P[( j -1) * K +1,i] - Lowpass2( j,i). (11)

На рис. 8 представлены проекционные данные, отфильтрованные разработанным методом (формулы (5)—(11)), и восстановленные изображения.

и) к) л) м)

Рис. 8: а, б - изображения с единичными ошибками каналов; в, г - изображения с многократной ошибкой единичного канала; д-з - изображения с неисправным каналом; и-м - изображения с множественными неисправными каналами

Можно сделать вывод о том, что метод позволяет обнаруживать и подавлять ошибки каналов, распространяющиеся на одну или несколько проекций (вплоть до угла разворота в 30°) и проявляющиеся как артефакты колец на восстановленных изображениях.

Заключение

Проведенные исследования по моделированию артефактов томографического изображения показали, что все без исключения рассмотренные артефакты существенным образом ухудшают качество получаемого изображения.

Получены признаки проявления артефактов для различных состояний детектирующей системы КТ.

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

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

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

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

Литература

1. Computed Tomography: Principles, Design, Artifacts, and Recent Advances. SPIE Press Monograph Vol. PM114 by Jiang Hsieh, 2003.

2. Гонсалес, Р. Цифровая обработка изображений: пер. с англ. / Р. Гонсалес, Р. Вудс. - М.: Техносфера, 2005. - 1070 с.

3. Симонов, Е.Н. Физико-математические основы проектирования томографических рентгеновских компьютерных комплексов / Е.Н. Симонов. - М. : Российская Академия Естествознания, 2011. - 410 с.

4. Симонов, Е.Н. Томографические измерительные информационные системы. Рентгеновская компьютерная томография / Е.Н. Симонов. - М. : НИЯУ МИФИ, 2011. - 440 с.

5. Симонов, Е.Н. Физика визуализации изображений в рентгеновской компьютерной томографии /Е.Н. Симонов. - Челябинск: Издат. центр ЮУрГУ, 2013. -505 с. - В печати.

Ласьков Вячеслав Валерьевич, соискатель, инженер-программист ФГУП «РФЯЦ-ВНИИТФ им. акад. Е.И. Забабахина»; snzst86@gmail.com.

Симонов Евгений Николаевич, д-р техн. наук, профессор кафедры радиотехники, ЮжноУральский государственный университет, филиал в г. Кыштыме; e.n.simonov@yandex.ru.

Bulletin of the South Ural State University Series “Computer Technologies, Automatic Control, Radio Electronics”

2013, vol. 13, no. 3, pp. 13-22

COMPUTED TOMOGRAPHY ARTIFACTS ANALYSIS, SIMULATION AND REDUCTION

V.V. Las’kov, FSUE “RFNC-VNIITF”, Russian Federation, snzst86@gmail.com,

E.N. Simonov, Branch of South Ural State University in Kyshtym, Russian Federation, e.n. simonov@yandex. ru

The analysis of the artifacts of projection data and images in the X-ray computed tomography has been done. The mechanism for modeling artifacts, which developed a universal software reconstructor tomographic images using filtered back projection convolution into several angles of exposure (projections) and single detector (counts and measurements in projection) with different functions for any configuration section of the object, has been proposed. Reconstructor allows to solve direct problem - defining projections and inverse problem - image reconstruction from projections. Modeling artifacts allows developing filters to suppress artifacts and assess the effectiveness of filters used. An effective method to reduce artifacts in tomographic images based on filtering the projection data has been proposed. The developed method allows to take into account peculiarities of obtaining projection data from the detectors to conduct optimal filtering distorted projections from the matrix of the projection data and, eventually, to suppress many types of artifacts in the tomographic image.

Keywords: computed tomography, artifacts, filtering procedure.

References

1. Computed Tomography: Principles, Design, Artifacts, and Recent Advances. SPIE Press Monograph Vol. PM114 by Jiang Hsieh, 2003.

2. Gonzalez R., Woods R. Tsifrovaya obrabotka izobrazheniy. Translated from English [Digital image processing]. Moscow, Technosfera, 2005. 1070 p.

3. Simonov E.N. Fiziko-matematicheskie osnovy proektirovania tomograficheskikh rentegenovskikh komp'yuternykh kompleksov [Physical and Mathematical Basis of the Design of X-ray Tomographic Computer Systems]. Moscow, Russian Academiya of Natural History, 2011. 410 p.

4. Simonov E.N. Tomograficheskie izmeritelnye informatsionnye sistemy. Rentgenoskaya komp'yu-ternaya tomografiya [Tomographic Measuring Information System. X-ray Computed Tomography]. Moscow NRNU MiFi, 2011. 440 p.

5. Simonov E.N. Fizika vizyalizatsii izobrazheniy v rentgenovskoy komp'yuternoy tomografii [Physics of Imaging Visualization in X-ray Computer Tomography]. Chelyabinsk, Publishing Center of SUSU, 2013. 505 p. In print.

Поступила в редакцию 7 апреля 2013 г.