Повышение контрастности малоракурсных томограмм, полученных алгебраическими алгоритмами реконструкции Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 14, № 3, 2009

Повышение контрастности малоракурсных томограмм, полученных алгебраическими алгоритмами реконструкции*

А. В. Лихачев Учреждение Российской академии наук Институт автоматики и электрометрии СО РАН, Новосибирск, Россия

e-mail: [email protected]

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

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

Введение

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

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

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 06-01-81000 и № 07-07-00251-а).

© ИВТ СО РАН, 2009.

двух переменных набор интегралов от нее вдоль всевозможных прямых, лежащих в

плоскости [1, 2]. Оно записывается следующим образом:

+^

/(<р,1) = (Яд)(р,/) =| д(р + агс1ап(р//),^/2 + р2)<!р , (1)

—те

где 1,р £ Я1, р £ [0;2п]. Функция одной переменной /^(1) = /(р, 1)|^=со^ называется одномерной параллельной проекцией функции д(ф,т), далее для краткости будем иногда называть ее просто проекцией. Согласно (1), значение проекции /(/) в некоторой точке есть интеграл от функции д(ф,т) вдоль прямой, заданной углом наклона р и расстоянием до начала координат I. Следует отметить, что для регистрации параллельной проекции в описанной выше схеме измерений источник излучения должен быть унесен на бесконечность. В реальной ситуации такое невозможно. Тем не менее это не умаляет практического значения проведенных в работе исследований. Во-первых, при достаточно большом удалении источника от объекта погрешность, обусловленная непараллельностью лучей, вдоль которых производятся измерения, незначительна. Во-вторых, может быть проведена пересортировка проекционных данных, в результате которой будет получен набор проекций /^(1) (см., например, [1]).

Для двумерного преобразования Радона (1) имеет место формула обращения, которая служит основой для построения алгоритмов томографической реконструкции [1-5]. Эта формула может быть записана в различных эквивалентных формах. Ниже приведена та из них, по которой был построен один из алгоритмов, используемых в работе:

1 717 )

д(Ф,т) = 4Л1 (У 1и )ехр сОй(Р - ^) ^ ' (2)

0 <—те /

где через ) обозначен фурье-образ проекции. Внешний интеграл по углу в вы-

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

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

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

1. Алгоритм реконструкции

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

gW = ф (A-1g(n-1) ,gFB) , n = 1, 2 ... (3)

Здесь А-1 — оператор, соответствующий итерации основного алгоритма, т.е. основной алгоритм имеет вид g(n) = A-1g(n-1); Ф — оператор, переводящий две функции в одну (см. формулу (6)); gFB — томограмма, полученная алгоритмом фильтрации и обратного проецирования.

Для приведенных в настоящей работе результатов основным алгоритмом является ART (Algebraic Reconstruction Technique) [1], широко используемый при решении задач томографии. Сходные результаты были получены и тогда, когда в качестве основного алгоритма использовались SIRT [1] и MART [6]. Согласно [1] (k + 1)-й шаг для ART дается формулой

- /ai(k)

д(к+1) = д(к) + А/г(к) ^ (4)

В уравнении (4) аг(к) — г(к)-я строка матрицы, полученной в результате дискретизации задачи томографии, эта матрица называется проецирующей; /г(к) — правая часть г(к)-го уравнения системы, т. е. показание отдельного детектора; (о, о) и || о || — скалярное произведение и эвклидова норма в конечномерном векторном пространстве соответственно; А — параметр релаксации. Проецирующая матрица имеет размерность I х 3, где I — полное число измерений (т. е. число детекторов в линейке, умноженное на число положений источника); 3 — число элементов (пикселов), на которые разбита область реконструкции. Элемент проецирующей матрицы а^ часто определяют как длину пересечения луча, соединяющего отдельный детектор и источник (положение детектора на линейке и поворот источника определяются индексом г) с ]-м пикселом. Зависимость индекса г от шага алгоритма к имеет вид г(к) = [к(шо^) + 1], таким образом, строки проецирующей матрицы перебираются циклически. Итерация, которая имеется в виду в формуле (3), состоит из I шагов алгоритма (4), т. е. в процессе итерации перебираются все проекционные данные. В [1] было показано, что итерационный процесс (4) сходится для любого начального приближения д(0) € &, если 0 < А < 2.

Томограмма дРБ реконструировалась регуляризующим алгоритмом фильтрации и обратного проецирования, разработанным в [5]. В формуле обращения (2) фильтрация производится в частотной области: фурье-образ проекции умножается на модуль частоты. При ее численной реализации возникают трудности, связанные с наличием бесконечных пределов интегрирования, поскольку спектр проекции, вообще говоря, неограничен. Поэтому при построении алгоритмов на основании уравнения (2) вместо | используют некоторую функцию ), обладающую при V ^ то нужными свойствами. В [5] одной из предложенных в качестве функций П^) была следующая:

О (. ) [ ^1 ехр(-а^|2) ^ < ^ (5)

°аМ=10, 21V |й> 1, (5)

где h — шаг сетки, задаваемой на проекции. Параметр регуляризации а вычислялся в [5] исходя из критерия невязки [7].

На рис. 1, а представлен математический фантом, использованный в настоящей работе в процессе вычислительного эксперимента. На рис. 1, б и в показаны его томограммы, полученные по 25 проекциям. Томограмма на рис. 1, б получена алгоритмом ART, а томограмма на рис. 1, в — алгоритмом фильтрации и обратного проецирования с регуляризующим фильтром (5). Видно, что томограммы на рис. 1, б и в существенно отличаются друг от друга. Решение, полученное посредствам ART, более размыто, нежели решение, изображенное на рис. 1, в, однако последнее имеет более интенсивные артефакты. Цель настоящей работы — получить такую томограмму, которая в местах разрывов имела бы такое же качество, как томограмма, реконструированная алгоритмами фильтрации и обратного проецирования, а в остальной части представляла бы собой решение, полученное алгебраическим алгоритмом.

Исходя из сформулированной цели, оператор Ф в формуле (3) строили следующим

Рис. 1. Математический фантом (а); томограммы, число проекций 25, уровень шума £ = 0 %: б — ART, в — фильтрация и обратное проецирование, г — предлагаемый комбинированный алгоритм (3)

образом:

Ф(01,52)(*,у) = | g1 И' |glИ - g |/ |g |- (6)

В выражении (6) д — среднее значение, взятое по большой области томограммы ниже поясним, как выбирать эту область. Через д1(ж,у) обозначена средняя величина томограммы в окрестности точки (ж, у), через е — положительная пороговая константа. Таким образом, если величина томограммы в окрестности некоторой точки незначительно отличается от окружающего фона, то оператор Ф не изменяет в ней томограммы. Если же имеет место обратная ситуация, то оператор Ф заменяет величину д^ж,^) на значение томограммы д2 в этой же точке.

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

2. Вычислительный эксперимент

Эффективность предлагаемого комбинированного алгоритма была показана в процессе вычислительного эксперимента. При этом использовался математический фантом, изображенный на рис. 1, а. Он состоит из слабо вытянутого эллипса (отношение осей 10:11) c амплитудой, равной 1.0, на фоне которого расположено восемь прямоугольников. Четыре из них имитируют трещины, значение фантома в их пределах равно нулю; остальные — вкрапления более плотных материалов, их амплитуды равны 2.0, 1.9, 1.8 и 1.7.

Математический фантом и томограммы вычислялись на квадратной сетке 1025 х 1025 узлов. Проекционные данные моделировались как набор одномерных проекций, распределенных равномерно в интервале от 0 до 180°. На каждой проекции задавалась сетка в 1025 узлов. Реконструкция проводилась как в отсутствие случайного шума в проекционных данных, так и при его наличии. Шум предполагался гауссовым со средним, равным нулю, и переменной дисперсией, составлявшей в каждой точке проекции £ процентов от ее значения. В расчетах варьировались число проекций, уровень шума, а также величина порога e в выражении (6). Значение g вычислялось по центральной части (размером 500 х 500 пикселов) томограммы, полученной алгебраическим алгоритмом. Средняя величина рТО^у-) в окрестности узла (ж^у) вычислялась путем усреднения в окне 3 х 3.

Значение параметра релаксации Л для алгоритма АИ.Т (4) равнялось 1.0. Поскольку цель применения предлагаемого комбинированного алгоритма — получение томограмм с максимально выделенной мелкомасштабной структурой, величина параметра регуляризации а для фильтра (5) бралась ниже той, которая получалась исходя из критерия невязки. Для реконструкции по незашумленным данным значение а равнялось 0.00005, а для зашумленных проекций — 0.0001.

При оценке результатов реконструкции использовались как качественные, так и количественные критерии. К первым относятся визуальное сравнение полученных изображений и сравнение профилей вдоль 470-й и 605-й строк. Одна из них проходит через "трещины", вторая — через прямоугольники с более высокой по сравнению с фоном амплитудой.

Среднеквадратичные отклонения полученных томограмм от математического фантома вычислялись в прямоугольных областях вблизи неоднородностей, поскольку именно здесь точность реконструкции представляет наибольший интерес. Эти области показаны на рис. 1, а. Ниже при описании результатов вычислительного эксперимента область слева (ту, которая охватывает "трещины") будем обозначать через , а область справа — через Среднеквадратичные отклонения для них будем обозначать через А 1, А2 и вычислять согласно формуле

А,-

а № - ^

£ [9з )

мА 2

1, 2.

(7)

В формуле (7) суммирование ведется по пикселам, принадлежащим соответствующим областям.

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

а

Рис. 2. Зависимости ошибок Д1 (а) и Д2 (б) от числа проекций: кривая 1 — ЛИТ; 2 — фильтрация и обратное проецирование; 3 — комбинированный алгоритм

рис. 2, б — ошибка Д2. На рисунке кривые 1, 2 и 3 относятся, соответственно, к алгоритмам ART, фильтрации и обратного проецирования с регуляризующим фильтром (5) и алгоритму (3), (6), полученному путем их комбинации. Значения ошибок для ART и алгоритма (3), (6) приводятся после 10-й итерации. При таком количестве итераций приближенно достигается минимум зависимостей Д^п) и Д2(п) (см. ниже — рис. 4). По рис. 2 можно заключить, что при небольшом числе проекций (M < 50) предлагаемый алгоритм в области неоднородностей дает лучший результат, нежели ART или обратное проецирование с фильтрацией. Однако при увеличении числа проекций точность реконструкции всеми алгоритмами (относительно ошибок Д1 и Д2) становится примерно одинаковой. Для результатов, представленных на рис. 2, восстановление производилось по незашумленным данным.

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

а б

Рис. 3. Зависимости ошибок Ai (а) и Д2 (б) от уровня шума: кривая 1 — ART, 2 — фильтрация и обратное проецирование, 3 — комбинированный алгоритм

Рис. 4. Зависимости ошибок Ai (кривые 1 и 2) и Д2 (кривые 3 и 4) от номера итерации (25 проекций): кривые 1 и 3 — ART, кривые 2 и 4 — комбинированный алгоритм

О сходимости итерационных процессов можно судить по рис. 4, где приведены зависимости ошибок Д1 (кривые 1 и 2) и Д2 (кривые 3 и 4) от номера итерации. К алгоритму ART относятся кривые 1 и 3, а к предлагаемому комбинированному алгоритму — кривые 2 и 4. Реконструкция производилась по 25 проекциям. Как следует из рисунка, сходимость обоих алгоритмов имеет одинаковый характер. Существенное уменьшение

200 400 600

д

1.5

0.5

0

200 400 600 800 1000 0

в

1 5

0.5

1000 о

200 400 600

г

и

200 400 600

е

1000

1000

а

Рис. 5. Профили изображенных на рис. 1 томограмм: вдоль 470-й (а, в, д) и 605-й (б, г, е) строк; а и б — ART; в и г — фильтрация и обратное проецирование, д и е — комбинированный алгоритм (3)

ошибок приходится на первые 4-5 итераций. Далее процессы практически стабилизируются. Таким образом, для достижения лучшего качества реконструкции, как для алгоритма ART, так и для комбинированного алгоритма, требуется одинаковое число итераций (порядка 10, согласно проведенным расчетам). Поэтому разница во времени реконструкции этими алгоритмами будет зависеть только от разницы во времени, затрачиваемом на одну итерацию. Последнюю можно оценить как tps/n+т, где n — число итераций; tpB — время, необходимое для реконструкции томограммы gFB алгоритмом обратного проецирования с фильтрацией (очевидно, что томограмму gFB достаточно вычислить один раз в начале итерационного процесса); т — время фильтрации оператором Ф. В проведенном вычислительном эксперименте были получены следующие результаты для полного времени реконструкции по 25 проекциям: обратное проецирование с фильтрацией — 1.1 с, ART — 17.3, комбинированный алгоритм — 21.4 с (для двух последних алгоритмов приведено время для 10 итераций). Реконструкция производилась на персональном компьютере с тактовой частотой процессора 2.0 ГГц.

Рис. 6. Томограммы: число проекций 50; уровень шума £ = 3%; а — ART; б — фильтрация и обратное проецирование; в — комбинированный алгоритм (3)

Относительно величины порога е для оператора (6) были сделаны следующие выводы. В проведенном вычислительном эксперименте томограммы лучшего качества были получены при е = 0.05... 0.15. С увеличением уровня шума, а также с уменьшением числа проекций значение е целесообразно увеличить. Однако не рекомендуется делать его больше одной трети от абсолютной разности между амплитудами включения и фона, деленной на амплитуду фона.

Примеры полученных томограмм и их сечений представлены на рис. 1, 5 и 6. На рис. 1 показаны математический фантом (рис. 1, а) и три его томограммы (рис. 1, б—г), полученные по 25 свободным от шума проекциям. Томограммы на рис. 1, б—г относятся соответственно к алгоритмам ART, фильтрации и обратного проецирования с регуля-ризующим фильтром (5) и алгоритму (3), (6), полученному путем их комбинации. На рис. 5 изображены профили этих томограмм вдоль 470-й (а, в, д) и 605-й (б, г, е) строк (а, б — ART, в, г — фильтрация и обратное проецирование, д, е — комбинированный алгоритм).

На рис. 6 даны результаты реконструкции по 50 зашумленным проекциям: а — ART, б — фильтрация и обратное проецирование, в — комбинированный алгоритм. Уровень шума £ составлял 3%. Следует обратить внимание на то, что "размытые" томограммы были получены алгоритмом ART без применения каких-либо дополнительных процедур сглаживания, причем такой вид томограмм имел место при реконструкции как по незашумленным, так и по зашумленным проекциям.

Заключение

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

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

Список литературы

[1] Хермен Г.Т. Восстановление изображений по проекциям. Основы реконструктивной томографии: Пер. с англ. М.: Мир, 1983.

[2] Наттерер Ф. Математические аспекты компьютерной томографии: Пер. с англ. М.: Мир, 1990.

[3] Тихонов А.Н., Арсенин В.Я., Тимонов А.А. Математические задачи компьютерной томографии: М.: Наука, 1987.

[4] Лаврентьев М.М., Зеркаль С.М., Трофимов О.Е. Численное моделирование в томографии и условно-корректные задачи. Новосибирск: Изд-во ИДМИ НГУ, 1999.

[5] Лихачев А.В. Регуляризующая фильтрация проекций в алгоритмах двумерной томографии // Сиб. ЖВМ. 2008. Т. 11, № 2. С. 187-200.

[6] Censor Y. Finite series-expansion reconstruction methods // Proc. IEEE. 1983. Vol. 71. P. 409-419.

[7] Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач: М.: Наука, 1986.

Поступила в редакцию 4 мая 2007 г., в переработанном виде — 14 января 2009 г.