Свойства регуляризованного алгоритма Гершберга-Папулиса в задаче веерной томографии Текст научной статьи по специальности «Математика»

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

Том 13, № 6, 2008

Свойства регуляризованного алгоритма Гершберга—Папулиса

в задаче веерной томографии*

В. В. Пикллов Институт теоретической и прикладной механики им. С.А. Христиановича СО РАН, Новосибирск, Россия e-mail: [email protected]

Д. И. КАЗАНЦЕВ Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: [email protected]

Problems of computational tomography with a small number of views require application of algorithms utilizing a priori information regarding the solution. Gerchberg— Papoulis algorithm (G—P), based on the theorem of a central slice, is known to be one of the best iterative methods of a few-view tomography for a parallel geometry of scanning. In a fan-beam tomography case such algorithm has not been yet investigated, because of there is no an analog of the central slice theorem for this case. In the present contribution, a recently developed generalization of the central slice theorem for fan-beam geometry is described. A new iterative G—P algorithm is developed based on the generalization. Two versions of the numerical scheme are investigated. An influence of the superimposed random noise on accuracy of reconstructions is studied; regularization criteria are developed which allow suppression of noise in measurements.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 07-07-00085, 07-01-00318, 05-08-50308, 05-02-16896), Интеграционного проекта СО РАН № 2006-35.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

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

1. Теория

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

Рассмотрим стандартную схему двумерной веерной томографии (рис. 1).

Здесь область определения неизвестной функции д(х, у) содержится внутри единичного круга х2 + у2 < 1. Источник веерных лучей Е движется по окружности:

где Б = ОЕ — расстояние от источника до начала координат в повернутой на угол в относительно системы (х,у) системе координат (в,р), а детекторы расположены на прямой линии Б1Б2 (рис. 1). Каждый луч ЕА в веере для фиксированного угла в можно

х(в) = -Б яп в, у (в) = Б сов в,

Р

Рис. 1. Формирование проекции на плоском детекторе в схеме веерной томографии [8]

охарактеризовать углом 7 к центральному лучу и расстоянием д до центра координат О. Угол в есть угол между осью х и нормалью к лучу, а угол между лучом и осью х определен как £. Ясно, что зарегистрированный сигнал на линейке детекторов О1О2 можно рассматривать и на прямой О1О2, которую назовем виртуальным детектором.

Введем следующее проективное ("деформирующее") преобразование в повернутой системе координат (в,р):

u

s/ (l - = s/Q;

(1)

v = p.

При таком преобразовании система веерных лучей в координатах (х,у) перейдет в систему прямых лучей, параллельных координатной оси V. Однако и исходное изображение д(х, у) будет деформировано в зависящее от угла в новое изображение д(Х(и, V), У (и, V)), где

x = X(u, v) = u ^ 1 — DD\ cos в — v sin в; y = Y(u, v) = u — — J sin в + v cos в■

(2)

Обозначим через и = и(х,у), V = V(х, у) обратную замену переменных в (2). Введем (г,ф) — полярные координаты точки (х,у), а также (д, в) — нормальные координаты веерного луча, проходящего через эту точку (рис. 1), в = в + 7, Р = г вт(ф — в). Отметим, что измеренный сигнал становится нулевым для |з| > О ^О"2 — 1. Тогда (см. также [5, 9]) веерная проекция может быть выражена следующим образом:

f(u)

8(q — U(x, y))g(x, y)dxdy

| J| 8(q — U)g(x,y)dudv.

(3)

—<x —<x

—<x —<x

Переход от координат (х,у) к (и^) осуществляется при помощи якобиана преобразования 3, который, используя (2), вычисляется как

J

dx/du dx/dv dy/du dy/dv

1 — — j cos в — sin в

v v \

1 — sin в cos в

1 — D = Q.

Окончательно уравнение (3) в новых координатах может быть записано как

1

fe(u) = / / |Q|8(q — U)g(x,y)dudv

g(X(u, v),Y(u, v))dv.

cos y (u)

При выводе уравнения (4) использованы следующие соотношения:

8(q — U) = 8(q — r cos(0 — ф)) = = 8(D sin y — cos y[r cos^ — ф) — tg y r sin(в — ф)]) =

(4)

8

cos

. r sin(в — ф)\ /n

u I 1 +-^—— — r cos^ — ф)

D

| cos

-8(uQ — r cos^ — ф))

1

|Q cos y(u)|

8(u — u'),

(5)

1

' r cos(e - ф) s где и = -—-= — .Из геометр

os y (заметим, что в реальны

тан

итоге имеем, что модуль якобиана |q из-под знака дельта-функции.

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

f'eК) = J exp(-2niv

Таким образом новая теорема о ц

билинейной интерполяции ппояемонсг

обратной деформацией набора параллельных прямых с рис. 2, б. Несколько других примеров деформации томографических фантомов приведено в статье [8].

Итерационный алгоритм Г—П основан на применении теоремы о центральном сечении. Представим этот алгоритм согласно работе [10] в виде

д(п)(х,у) = -1~д(п-1\их,иу), п = 1, 2,...,

д(га)К,^) = ф(га-1Чд(п-1)(х,у), п =1,2,..., д(0)(х,у) = о, (7)

где F2 и F-1 — операторы прямого и обратного двумерного преобразования Фурье; ф("0, — соответственно операторы изменения п-го итерационного решения и его фурье-образа, которые отвечают за использование априорной информации. Ф(п) действует в области томограммы, Ф^га) действует в фурье-пространстве. В частности, в последний оператор входит процедура переноса фурье-спектра проекций, заданного на радиальных линиях в фурье-пространстве искомого изображения, на все двумерное пространство. Эта процедура является центральной частью алгоритма Г—П, учитывающей в нем измеренные данные — интегральные проекции.

На практике от объекта может быть получено только конечное число проекций и каждая проекция состоит из конечного числа отсчетов. Следовательно, фурье-образ экспериментально измеренных проекций / (ир) задан в наборе дискретных точек вдоль конечного числа радиальных линий. И для того чтобы осуществить обратное двумерное фурье-преобразование от функции / (^р), заданной в фурье-пространстве в полярной системе координат, необходимо определить значения соответствующего спектра в декартовой системе координат (ух,уу).

В алгоритме Г—П для параллельной геометрии (случай малоракурсной томографии) при переходе от полярной к декартовой системе координат хорошо зарекомендовал себя метод полосовой интерполяции [5]. Работа полосовой интерполяции осуществляется по следующему алгоритму: в частотной области вокруг каждой известной радиальной линии вводится полоса шириной 2Бщ, некоторый множитель 8s < 1 и период обновления ширины полосы равен т0. Внутри полосы осуществляются две одномерных интерполяции в текущий узел декартовой сетки: ближайшего соседа поперек полосы и линейной — вдоль полосы. В ходе работы алгоритма для той итерации с номером п, для которой выполнено условие шо^п,т0) = 0, обновляется значение ширины полосы по формуле

= 5га- . (8)

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

В дальнейшем в качестве оценки погрешности реконструкции вычислялась следующая норма отклонения точной томограммы от восстановленной Д,:

д,

\

N N

Е^ / „гее ,-.еха\ 2 — >

^^--100 %• (9)

Е Е (деха)2

г=1з=1

Здесь деха — точное решение; дгес — восстановленная томограмма.

Критерии останова итерационного алгоритма Г—П позволяют прервать итерационный процесс с ростом нормы невязки. Перечислим четыре исследованных здесь версии критерия останова по минимуму невязки:

1) по росту двух последних невязок;

2) по росту трех последних невязок;

3) по росту шести последних невязок;

4) по условию, когда ширина полосы становится меньше 1 пикселя, после чего включается 1-й критерий.

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

2. Численное моделирование для веерного алгоритма Г—П

Численная реализация веерного алгоритма Г—П была осуществлена двумя способами (модификации 1 и 2).

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

Деформирующее преобразование зададим "матрицей деформации" (1):

Т

Я-1 0 01

(10)

в которой Я = — . Посредством оператора Т осуществляется переход от системы

В.

координат (в,р) к новым координатам (и, у): ( ) = Т( ).

С помощью последовательного применения указанных операторов можно перейти к координатам (и, у), затем двумерное прямое преобразование Фурье дает возможность перейти к частотным координатам (ии, ^), в пространстве которых можно использовать теорему о центральном сечении, и далее, как и в алгоритме Г—П для параллельной геометрии, коэффициенты Фурье заменяются на известные, полученные применением одномерного фурье-преобразования к исходным проекционным данным (оператор Фf). Отметим, что подобная замена проводится в повернутой на угол в системе координат, в которой текущий угол веерного наблюдения всегда равен нулю. Как следствие, замена спектра в фурье-пространстве необходима лишь на оси абсцисс ии и требует лишь одномерной интерполяции. Обратное двумерное фурье-преобразование дает очередную оценку деформированной томограммы в координатах (и, у), и, чтобы вернуться к искомой функции д(х,у) (с учетом априорной информации через оператор ), необходимо осуществить обратную деформацию и поворот на угол — в.

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

X \ \ D / \ I m I S

g{s,p) = g(i-1) yjJ , g (u,v) = g\T

gD(vu,vv) = Ф/F2[gD(u,v)], gf (u, v) = Ф^-1^v^)], I rp-U u l l l - l s

gф(s,p) = gф ^T ^^^^ , gi(x,y) = gф . (П)

В данной реализации алгоритма все операции поворота реализованы в виде двумерной интерполяции, а прямая и обратная деформации изображения осуществляются при помощи построчной одномерной интерполяции, одинаковой для каждого ракурса (операторы T и T-1), так как при соответствующем повороте фокальная точка всегда оказывается расположенной на оси ординат v = p. Подчеркнем еще раз, что формулы (11) (в отличие от формул (7)) относятся лишь к одной итерации и описывают набор операций для обработки i-й проекции под углом Д.

Результаты, полученные с помощью алгоритма в модификации 1, представлены на рис. 3, где в качестве элементарного фантома рассмотрена модель № 224 (сдвинутая и повернутая гауссиана) пакета программ по вычислительной томографии TopasMicro [5].

Модификация 1 продемонстрировала хорошую сходимость итерационного процесса, но дала низкую точность реконструкции за счет большого числа интерполяций на каждом ракурсе (рис. 3, а). Увеличение точности реконструкции путем применения интерполяции более высокого порядка, а именно одномерной интерполяции эрмитовыми сплайнами [11], позволило повысить точность реконструкции. Результаты эксперимента с зашумленными данными и одномерным сглаживанием проекций (рис. 3, б) показали, что применение процедуры сглаживания в модификации 1 не приносит улучшения, а напротив, вносит искажения в решение. Чтобы уменьшить высокие погрешности модификации 1, необходимо сократить общее число интерполяций. Суть следующей модификации заключается в замене двух интерполяций, которые отвечают за поворот и деформацию томограммы, на одну двумерную, в которую одновременно входят и поворот, и деформация.

Модификация 2. Отметим, что процедуру деформации томограммы можно выполнить сразу для истинного положения фокальной точки, без ее предварительного переноса на ось ординат, как это делается в модификации 1 алгоритма. В этом случае операции двумерного фурье-преобразования происходят в пространстве неповернутой деформированной томограммы (x ,y ), которое связано с исходным пространством матричными соотношениями

/ //;\

7' \ / /" \ / i" \ / i"

-1

у,/= ^КуУ и=^МУ! (12)

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

8 10 12 14 16 18 П

8 10 12 14 16 18 п

Рис. 3. Зависимость ошибок реконструкции от числа итераций для модели № 224. Модификация алгоритма № 1 с параметрами К = 13, В = 1.5, N = N = 128; а — одномерные интерполяции: линейная (кривые 1-3), эрмитовыми сплайнами (4-6) с параметрами полосовой интерполяции = 1.8, 5,3 = 0.8, то = 3 (1, 4), то = 1 (2, 5), = 0.9 (3, 6); б — параметры для кривых (1-8): = 1.8, = 0.8, то = 1; интерполяции: линейная (1-3), эрмитовыми сплайнами (4-8); уровень случайного шума в проекциях: к = 0% (кривые 1, 4), к = 3% (7), к = 3%, одномерная фильтрация проекций (8), к = 10% (2, 5), к = 10%, одномерная фильтрация (3, 6); в — точная томограмма, модель № 224; г — результат итоговой реконструкции для кривой 6 на а

а

в

г

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

На рис. 5 продемонстрированы результаты вычислительного эксперимента с применением двумерной интерполяции бикубическими сплайнами [12] (для модификации 2 алгоритма Г—П) к фантомам, различающимся по степени гладкости: модель № 224 (гауссиана, гладкая модель), 230 (парабола, частично гладкая), 235 (прямоугольник и ступенька в эллипсе, разрывная).

Рис. 4. Зависимость ошибок реконструкции от числа итераций для модели № 224. Модификация алгоритма №2, K = 13, D = 1.5, Nx = Ny = 128, ôs = 0.8. Билинейная интерполяция, параметры кривых: 1 — Sw = 5.5, m0 = 3; 2 — Sw = 5.5, m0 = 1; 3 — Sw = 1.8, m0 = 3; 4 — Sw = 1.8, m0 = 1; 5 — Sw = 0.9, m0 = 3; б — Sw = 0.9, m0 = 1

Рис. 5. Зависимость ошибок реконструкции от числа итераций для различных моделей. Модификация алгоритма № 2 с использованием бикубической интерполяции для моделей 224 (кривая 1), 230 (2) и 235 (3) соответственно; кривая (4) — модель 235 с процедурой очистки

По результатам, отраженным на рис. 5, видно, что разрывные, сложные модели восстанавливаются значительно хуже гладких. Гладкие модели (кривые 1, 2) дают при реконструкции быструю сходимость алгоритма и существенно меньшую ошибку реконструкции (3.5%). При использовании априорной информации о положительности и пространственной ограниченности томограмм, удается значительно улучшить реконструкцию. К разрывной модели № 235 применена "очистка" [5], вытекающая из свойства положительности томограммы, в основе такой "очистки" лежит процедура обратного проецирования. На выбранном ракурсе зануляются те участки томограммы, где обратное проецирование дало нулевое значение (из-за нулевых участков проекций). После применения процедуры очистки ошибка реконструкции снизилась (кривая 4 на рис. 5).

Аналогично эксперименту с модификацией 1 алгоритма (рис. 3, б) по влиянию шумов в проекционных данных на точность реконструкции, подобное влияние исследовалось и для модификации 2 алгоритма. Результаты этого эксперимента представлены на рис. 6, и они значительно отличаются от аналогичных результатов для модификации 1. На случайных шумах, меньших 5%, веерный алгоритм Г—П устойчив к шуму и обладает регуляризирующим свойством, однако на более высоких шумах необходимо сглаживание проекций. Это продемонстрировали реконструкции на рис. 7, б, в. Здесь показано, что применение сглаживания к проекциям позволяет не только избавиться от артефактов возле объекта (рис. 7, б), но и улучшить структуру самой томограммы (рис. 7, в). Демонстрируется также положительное влияние процедуры "очистки" на примере сложной разрывной модели № 235. На рис. 7, г можно видеть значительное количество низкоамплитудных артефактов, сопровождающих реконструкцию этой сложной модели. С помощью "очистки" удается в значительной мере избавиться от их присутствия на реконструкции (рис. 7, <9).

Рис. 6. Зависимость ошибок реконструкции от числа итераций и уровня шума для модели № 224; модификация алгоритма № 2, бикубическая интерполяция с параметрами К = 13, В = 1.5, N = Му = 128, Бы = 1.8, 53 = 0.8, то = 1; уровень шума к = 0% (кривая 1), к = 10% (2), к = 10%, одномерная фильтрация проекций (3)

0

I

0.9

0.5

■ 1.2

1

0.6 0

-1

! 0

%

-1 0 1

г

(

I

-1 0 1

1.2

0.6

в

е

Рис. 7. Результаты реконструкции модели № 224 (а—в) и модели № 235 (г — е); алгоритм Г— П в модификации 2 с параметрами К = 13, В = 1.5, М = Му = 128, ё8 = 0.8; а—в — реконструкции с использованием бикубической интерполяции, для уровня случайных шумов в проекциях 0% (а), 10% (б) и 10% (со сглаживанием проекций) (в). Ошибка реконструкции на 20-й итерации А, = 3.7% (а), 13% (б), 12.3% (в). Параметры реконструкции модели № 235: бикубическая интерполяция, п = 8, А, = 22.3% (г); бикубическая интерполяция, очистка томограммы, п = 12, А, =21% (д); билинейная интерполяция, очистка томограммы, М = Му = 512, Мр = 128, п = 20, А, = 19.4% (е)

Время счета одной итерации для модели № 224 при различных размерностях томограммы

Размерность (Nx = Ny, Np = 128) Время счета (одна итерация), с Ai, %

128 6 16.7

512 86 7.1

1024 295 4.5

Примечание. Эксперимент произведен на модификации 2 веерного алгоритма Г—П с использованием билинейной интерполяции.

Результаты веерного алгоритма Г—П модификации 2 с бикубической интерполяцией являются наилучшими на данный момент, но основная проблема заключается в большом времени счета: одна итерация данного алгоритма в программной среде Matlab выполняется около 80 мин. Были предприняты попытки обойти эту проблему увеличением размерности томограммы (без увеличения числа лучей в веере) и использованием билинейной интерполяции вместо бикубической. Полученные результаты оказались весьма успешными, так как на размерностях (Nx = Ny = 512 и 1024) удалось приблизиться к точности бикубической интерполяции, сократив время счета в десятки раз.

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

Таким образом, имея неизменное количество детекторов, в данном случае Np = 128, можно достаточно точно восстановить томограмму и на больших размерностях. В таблице приведена зависимость времени счета при увеличении размерности томограммы с использованием билинейной интерполяции в модификации 2 на примере модели № 224. Численное моделирование проведено на компьютере с процессором AMD Athlon 3000+ (2.01 ГГц). Результаты в таблице необходимо сравнить с наилучшими результатами, полученными при помощи бикубической интерполяции. Как уже упоминалось, бикубическая интерполяция на размерности Nx = Ny = 128 требует порядка 80 мин на одну итерацию, что в 800 раз больше, чем билинейная интерполяция на той же размерности. Билинейная интерполяция имеет весьма большую ошибку реконструкции на размерности Nx = Ny = 128 (Ai = 16.7%), однако при увеличении размерности до Nx = Ny = 512 удается уменьшить погрешность в 2 раза (A1 = 7.1 %). При этом время счета по сравнению с размерностью Nx = Ny = 128 увеличилось в 15 раз, но относительно бикубической интерполяции оно все равно осталось в 55 раз меньше. При увеличении размерности томограммы до Nx = Ny = 1024 время счета уменьшилось в 16 раз по сравнению с бикубической интерполяцией, а ошибка реконструкции сравнима с ошибкой при использовании в алгоритме бикубической интерполяции на размерности Nx = Ny = 128 (A1 = 4.5%). В подтверждение данного вывода проведена реконструкция модели № 235 модификацией 2 алгоритма Г—П с билинейной интерполяцией, но на размерностях Nx = Ny = 512, Np = 128 (рис. 7, е). После применения процедуры очистки видно, что результат реконструкции и погрешность почти идентичны восстановлению этой же модели с использованием бикубической интерполяции (рис. 7, д), однако реконструкция на рис. 7, е получена в десятки раз быстрее.

Заключение

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

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

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

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

[1] Defrise M., De Mol C. A regularized iterative algorithm for limited-angle inverse Radon transform // Optica Acta. 1983. Vol. 30, N 4. P. 403-408.

[2] Sato T., Norton S.J., Linzer M.J., Ikeda U., Hirama M. Tomographic image reconstruction from limited projections using iterative revisions in image and transform spaces // Applied Optics. 1981. Vol. 20, N 3. P. 395-399.

[3] Вишняков Г.Н., Гильмлн Г.А., Левин Г.Г. Восстановление томограмм при ограниченном числе проекций. Итерационные методы // Оптика и спектроскопия. 1985. Т. 58, № 2. C. 406-413.

[4] Kak A.C., Slaney M. Principles of computerized tomographic imaging. N.Y.: IEEE Press, 1988.

[5] Пикалов В.В., Мельникова Т.С. Томография плазмы. Новосибирск: Наука, 1995.

[6] Гуварени Н.М. Вычислительные методы и алгоритмы малоракурсной компьютерной томографии. Киев: Наукова думка, 1997.

[7] Pickalov V.V., Kazantzev D.I., Ayupova N.B., Golubyatnikoy V.P. Considerations on iterative algorithms for fan-beam tomography scheme // Proc. 4-th World Congress on Industrial Process Tomography. Aizu, Japan, 5-8 September 2005. Vol. 2. P. 687-690.

[8] Пикалов В.В., Казанцев Д.И., Голубятников В.П. Обобщение теоремы о центральном сечении на задачу веерной томографии // Вычисл. методы и программирование. 2006. Т. 7, № 2. С. 180-184.

[9] Minerbo G. Maximum entropy reconstruction from cone-beam projection data // Comput. Biol. Med. 1979. Vol. 9, N 1. P. 29-37.

[10] Пикалов В.В., Лихачев А.В. Применение метода Гершберга—Папулиса в трехмерной доплеровской томографии // Вычисл. методы и программирование. 2004. Т. 5, № 2. С. 27-34.

[11] Завьялов Ю.С., Квлсов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980.

[12] Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980.

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