CC BY
23
Вычислительные технологии
Том 19, № 2, 2014
ROI-томография по данным, содержащим шум с переменной дисперсией
А. В. Лихачёв
Институт автоматики и электрометрии СО РАН, Новосибирск, Россия
e-mail: [email protected]
Лихачёв А.В. ROI-томография по данным, содержащим шум с переменной дисперсией // Вычисл. технологии. 2014. Т. 19, № 2. С. 62-75.
Разработаны два метода решения задачи ROI-томографии, в которой область, рассматриваемая в качестве предмета исследования, расположена в центральной части изображения, а проекционные данные содержат шум, имеющий наибольшую интенсивность на периферии. В первом методе проекции умножаются на функцию окна, во втором — делятся на несколько частей, которые отдельно фильтруются, после чего суммируются с весами. Получены оценки усиления шумовой составляющей. Вычислительный эксперимент показал, что ошибка реконструкции с ростом уровня шума для классического подхода, состоящего в непосредственной реализации формулы обращения преобразования Радона, растёт значительно быстрее, чем для предлагаемых методов.
Ключевые слова: ROI-томография, нестационарный шум.
Likhachov A.V. ROI-tomography from the data distorted with variable dispersion noise // Comput. Technologies. 2014. Vol. 19, No. 2. P. 62-75.
Two methods are developed to solve the following problem of ROI-tomography: the region of interest is placed in the central part of the image, while the projection data contain noise that has its maximal intensity on the periphery. In one of these methods, a multiplication of the projections by a function of window is applied. Using another method, one divides each projection into several parts, which are separately filtered and multiplied by the weights, and then the obtained functions are summarized. The estimations of noise component amplification are derived. Numerical simulations show that a classical approach, which consists in the immediate implementation of the inverse Radon transform formula, provides much faster growing of the reconstruction error when the noise level increases compared to the proposed methods.
Key words: ROI-tomography, non-stationary noise.
Введение
В приложениях предметом исследования часто является некоторая, возможно небольшая, область внутри объекта, остальные же его части интереса не представляют. Такие задачи относятся к ROI-томографии (от английского region of interest). Их изучению посвящено большое число работ [1-5]. Сходными вопросами занимается также локальная томография. Этот термин используется в случае, когда для определения свойств функции в некоторой точке (например, её значения или наличия разрыва) требуют-
ся только данные, зарегистрированные на лучах, пересекающих сколь угодно малую окрестность данной точки [2, 5].
Будем рассматривать двумерный случай. Примем лучевое приближение, согласно которому реконструкция внутренней структуры объекта сводится к нахождению функции двух переменных, ниже она обозначается как g(x,y), по её интегралам вдоль прямых линий [6]. Основой многих алгоритмов лучевой томографии является преобразование Радона [6, 7], которое записывается следующим образом [6]:
сю
{Rg)ip, ф) = f (p, ф) = J g (у^ + Т2, ф + arctgp^j di. (i)
— с
Здесь f (p, ф) — значение интеграла от g(x,y) вдоль прямой, заданной расстоянием до начала координат p и углом ф между вектором нормали и положительным направлением оси X. Функцию f (p, ф) называют проекционными данными.
Для понимания трудностей, возникающих в задачах ROI-томографии, обратимся к формуле обращения преобразования (1), которая может быть представлена в виде [7]
п сю
g(x,y) = 2 J J luf (u, ф) exp(iu(x cos ф + y sin ip))dwd<p, (2)
0 с
где f(u, ф) — преобразование Фурье от f (p, ф) по первой переменной. Зафиксируем угол ф. Назовём функцию fv(p) = f (p,ф) |^=const одномерной проекцией. Формула (2) определяет следующий порядок операций. Каждая проекция фильтруется путём умножения её спектра на модуль частоты, в результате чего получается функция f'F(p). В англоязычной литературе эта процедура известна как гатр-фильтрация. Затем производится обратное проецирование — интегрирование отфильтрованных данных по углу.
Предположим, что посредством формулы (2) вычисляется функция g(x,y) в точке (x0,y0). Обратное проецирование суммирует значения (p0(ф)), где p0(ф) = y0 cos ф — x0 sin ф. Таким образом, оно является локальным в том смысле, что для вычисления g(x0,y0) используются только данные, зарегистрированные вдоль прямых, проходящих через (x0,y0). Подробнее об этой процедуре см. [6]. Однако функции (p) получаются путём гатр-фильтрации одномерных проекций, которая, очевидно, требует всех значений ftp(p). Следовательно, чтобы восстановить искомую характеристику объекта хотя бы в одной точке, нужно иметь проекционные данные от всего объекта, т. е. в целом формула обращения (2) является нелокальной.
Пусть функция g(x,y) имеет носитель supp(g). Требуется найти её в области D С supp(g), при этом по какой-либо причине не возможно или не желательно использовать все проекционные данные. Тогда из приведённого выше рассуждения следует, что по формуле (2) точное решение задачи получить нельзя. В этом случае для оценки g(x,y) в области, представляющей предмет исследования, проекции обрабатываются каким-либо фильтром высоких частот, имеющим локальные свойства. В работе [8] для этой цели предложен оператор Лапласа. Позднее было показано, что таким способом можно точно восстановить места разрывов функции g(x,y), хотя сами её значения оказываются неверными [9, 10]. Для ряда приложений, например для дефектоскопии, этой информации оказывается вполне достаточно, чтобы провести диагностику объекта.
Другой подход к проблеме восстановления положения разрывов разрабатывается в [11, 12]. Здесь используются интегро-дифференциальные операторы первого порядка,
что имеет определённые преимущества. Однако в этих работах рассматривается более общая задача томографии, в частности, учитываются распределение фотонов по энергии и рассеяние зондирующего излучения. Поэтому для её решения необходимо много проекционных данных, зарегистрированных посредством сложной системы измерений.
Трудность применения оператора Лапласа состоит в том, что он по сравнению с гашр-фильтром в большей степени усиливает высокочастотную компоненту шумовой составляющей, поскольку в частотной области ему соответствует умножение на ш2. До некоторой степени эта проблема может быть решена при помощи сглаживания одномерных проекций регуляризующими сплайнами [13, 14]. Ситуация осложняется, если шум является нестационарным. В этом случае крайне затруднительно найти единое для всех данных значение параметра регуляризации.
Таким образом, локальная томография нуждается в дальнейших исследованиях, в том числе по разработке новых алгоритмов реконструкции объекта в заданной области. В настоящей работе рассматривается одна из задач, представляющая существенный практический интерес.
1. Постановка задачи
Не уменьшая общности, можно полагать, что функция д(х,у), описывающая распределение какого-либо параметра внутри объекта, имеет компактный носитель, заключённый в единичном круге. При этом каждая одномерная проекция определена на отрезке [—1; 1]. Если объект полностью сканируется во всех выбранных направлениях, то для его реконструкции могут быть использованы алгоритмы, реализующие формулу обращения двумерного преобразования Радона (2).
Предположим, что проекции содержат нестационарный аддитивный шум, причём его дисперсия монотонно возрастает с увеличением модуля аргумента. Требуется оценить функцию д(х,у) в круге радиуса г0 < 1 с центром в начале координат. Несмотря на то что эта задача может быть решена средствами обычной, нелокальной, томографии, применение её методов может привести к значительным ошибкам. Действительно, информация об объекте в заданной области содержится только в центральных частях проекций, где уровень шума наиболее низкий. Однако благодаря нелокальной фильтрации на результат реконструкции будет влиять интенсивный шум с периферии. При этом очевидно, что точность метода понизится.
И,01-томография в рассматриваемом случае представляется более перспективной. Однако хорошо изученная локальная фильтрация посредством оператора Лапласа, с одной стороны, не обеспечивает получения оценок значений функции д(х,у), с другой — приводит к сильному росту шумовой составляющей. В связи с этим в работе предложены и исследованы два новых метода решения сформулированной задачи, которые были названы метод окна и метод деления.
2. Предлагаемые методы 2.1. Метод окна
Поскольку область, представляющая предмет исследования, в каждом направлении проецируется на отрезок [—г0; г0], полезная информация заключена только в его пределах. Остальные данные используются для того, чтобы отделить её от информации
о других частях объекта. Идея метода состоит в том, чтобы произвести эту операцию посредством окна. Таким образом, из рассмотрения будут исключены измерения, наиболее искажённые шумом. В работе реализован следующий алгоритм. Одномерные проекции умножаются на некоторую чётную функцию, равную единице в начале координат и не возрастающую на интервале [0; 1]. Далее они фильтруются гашр-фильтром. Затем проводится процедура обратного проецирования.
Применение окна искажает спектры проекций. При этом существенную роль играет его форма. В частности, гладкие окна вносят меньшую ошибку, чем окна, содержащие разрывы. Систематическое изложение данного вопроса можно найти, например, в [15]. В проведённом вычислительном эксперименте использовались следующие известные функции: окно Хэмминга
Шн(р) = 0.54 + 0.46 ес8 (пр/ртах),
(3)
окно Блэкмена
окно Парзена
Шв (р) = 0.42 + 0.5 осе (пр/ртах) + 0.08 осе (2пр/ртах),
Шр (р) =
1 - 6(р/ртах)2 + 6|р/ртах|3, 0 < |р| < 0.5ртах,
0.5р тах < |р| < рп
|р| > ртах.
2(1 - |р/ртах|)3
0,
В формулах (3)-(5) ртах > 0 и -ртах < р < рП
2.2. Метод деления
Данный метод основан на подходе, аналогичном предложенному Бартлеттом для вычисления периодограмм случайных сигналов [16]. В этом случае проекция представляется в виде суммы К > 1 слагаемых:
К/2-1 К/2-1
/^(р) = X] = X /(р)п2/к
к=-К/2 к=-К/2
р
2^ + 1 К
/(р)
(К-1)/2
Е
к=-(К-1)/2
(К-1)/2
/^(^,р) = X Шп1/к (р - К)
к=-(К-1)/2 ^ '
(7)
Формулы (6) и (7) записаны для чётных и нечётных К соответственно. Через Щ/к(р) обозначен прямоугольный импульс единичной амплитуды, имеющий ширину 2/К, центр которого расположен в начале координат. Таким образом, &-е слагаемое в (6) равно одномерной проекции при р € [2^/К; 2(& + 1)/К] (для формулы (7) при р € [(2& — 1)/К; (2^ + 1)/К]) и нулю в остальной области. Отметим, что равенства (6) и (7) не выполняются на границах интервалов.
Каждая из функций /(&,р) фильтруется посредством гашр-фильтра. Оценка фильтрованной проекции вычисляется путём суммирования полученных функций /^(&,р)
с весами. Вес должен быть тем меньше, чем выше уровень зашумления на к-м отрезке. В работе он определяется как
Здесь а — параметр, принимающий значения 0 < а < 1, О(р) — переменная дисперсия шума. Интегрирование в (8) ведётся по к-му отрезку. В знаменателе подкоренного выражения стоит наибольшее для рассматриваемых отрезков значение интеграла, по условию задачи достигаемое при минимальном или максимальном к.
Преимуществом метода деления является то, что в нём используются все имеющиеся данные. При этом измерения, сильнее искажённые случайным шумом, входят в результат с меньшим весом. Однако применение численно реализованной фильтрации приводит к нарушению равенств (6) и (7). Кроме того, сложение с весами, хотя и подавляет шум, но, с другой стороны, вносит дополнительные ошибки.
3. Оценки усиления шумовой составляющей в проекциях от области, представляющей предмет исследования
3.1. Общие замечания
В современных алгоритмах томографической реконструкции, реализующих формулу обращения (2), фурье-образ одномерной проекции умножается не на модуль частоты, а на произведение |ш|Ф(ш), где Ф(ш) — передаточная функция фильтра низких частот. Низкочастотная фильтрация требуется как для повышения устойчивости алгоритма, так и для подавления боковых лепестков спектра проекции, возникающих в результате дискретного преобразования Фурье. В настоящее время используется большое число различных функций Ф(ш). Наиболее распространённые из них приведены в монографиях [6, 7]. Там же дан подробный анализ их свойств. Здесь мы только отметим, что Ф(ш) является чётной функцией частоты, т.е. Ф(ш) = Ф(-ш). Выбор того или иного фильтра зависит от условий задачи. Введём обозначение 0,(ш) = |ш|Ф(ш). Очевидно, что 0,(ш), так же как и Ф(ш), — чётная функция.
Назовём данные, зарегистрированные в отсутствие шума, точными и обозначим через ¡0(р, ф), а результат их фильтрации посредством умножения спектра на 0,(ш) — через ¡^ (р, ф). Пусть вместо точных проекций известны функции ¡р (р) = ¡0>1р (р) + ^ (р). Здесь (р) — нестационарный случайный сигнал с математическим ожиданием тр (р), дисперсией (р) и корреляционной функцией Яр(р\,р2). Будем считать, что статистические свойства шума во всех направлениях сканирования одинаковы. Следовательно, функции тр(р), О(р) и Яр(р\,р2) не зависят от угла ф, поэтому нижний индекс в их обозначениях может быть опущен. Пусть т(р) = 0, тогда О(р) = Я(р\,р2)|Р1=Р2=Р. Предположим далее
т. е. шум является некоррелированным.
После фильтрации получим (р) = (р) + ^ (р). Определим среднюю ошибку, содержащуюся в отфильтрованных проекционных данных от интересующей нас области,
(8)
Я(рг,р2) = В(рг)6(р1 - р2),
обусловленную наличием шума:
N
N
П=1
\
N
М{/(р, ^п) - /Т(р, Ы)2}Ф = - £
-го
п=1
\
М {(£Т (р,^,))2}ф. (10)
-го
Здесь N — количество ракурсов наблюдения, М{о} — математическое ожидание. Нетрудно видеть, что гашр-фильтрация и другие операции над проекциями, входящие в предлагаемые методы, оставляют шум центрированным, а его дисперсию — изотропной по направлениям, обозначим её через (р). Тогда (10) можно записать в виде
1
N
N
п=1
\
го
(р)ф
-го
\
го
(р)ф.
-го
Для различных методов решения рассматриваемой задачи функции (р) будут отличаться друг от друга. Ниже в разделах 3.2 — 3.4 получены их выражения.
3.2. Реализация формулы обращения
В данном случае фильтр П(р) применяется непосредственно к имеющимся зашумлён-ным одномерным проекциям. Ниже для краткости такой метод оценки фантома внутри области, представляющей интерес, будет называться классическим. Дисперсию (р) определим по корреляционной функции, которая при воздействии на случайный процесс линейного оператора Т преобразуется по правилу Т1Т2Д(р1,р2), где Т1 и Т2 означают оператор Т, действующий по переменным р1 и р2 соответственно. Воспользуемся представлением в пространстве Фурье. Учитывая (9), имеем
-те -те те
Др1Жр1,р2) ехр (-г(^1р1 + ^2р2))ф1ф2 Д(р1) ехр (-¿(ш^ + ^2р1))^р1 = ¿(ш + ^2).
:12)
Преобразованием является умножение на функцию П(ш). Таким образом, получим
КТ (ш1, ш2) = Я(ш1,ш2 )П(ш1)П(Ш2) = + ш2)П(ш1)П(ш2). (13)
Возьмём двумерное обратное преобразование Фурье от (13):
те те
Лт(р1,р2) = 4П2 J ! ^Ош + ^2)^(^1)1^(^2) ехр (¿(ш^ + ш2р2))^ш^ш2. (14)
-те -те
Произведём в (14) замену переменных ш = + ш2, ш' = ш2. Поскольку при р1 = р2 левая часть (14) есть дисперсия шума после фильтрации, запишем
БТ (р)
4п2
¿(ш) I / П(ш)П(ш - I ехр (гшр)^ш.
:15)
1
Внутренний интеграл в (15) представляет свёртку фурье-образа фильтра с самим собой. Согласно известной теореме он равен 2п(П2)(ш). Тогда
те те
Вр= В(ш)(№)(ш)ехр(1шр)с1ш = J О(р')П2(р - р')ф. (16)
— те —те
В итоге для ошибки (11) имеем следующее выражение:
0 =
\
г о
I (Б * П2)(р)ф. (17)
го
Здесь и далее запись (о * о) означает свёртку.
3.3. Метод окна
Обозначим через Ш(р) функцию окна. При умножении на неё одномерной проекции корреляционная функция шума будет
Яш(Р1,Р2) = Ш(рг)Ш(Р2)Я(Р1,Р2) = Ш2(р1)Б(р1)8(р1 - р2). (18)
Таким образом, имеем некоррелированный шум с дисперсией (р) = Ш2(р)Б(р). Используя результаты раздела 3.2, получим
0
\
у ((Ш2Б) * П2)(р)ф. (19)
—го
3.4. Метод деления
При определении ошибки 0 для метода деления будем считать, что носитель проекции разбит на нечётное число интервалов. Тогда для дисперсии шума имеет место представление, аналогичное формуле (7):
(К—1)/2 , .
Б(р)= £ Б(р)П2/Л р - - . (20)
к=—(К—1)/2 ^ '
Проведём гашр-фильтрацию слагаемых в правой части (20) и просуммируем их с весами (8). Учитывая (17), получим выражение
(К—1)/2
БП(Р)= ^ ^к(Бп,к * П2)(р). (21)
к=—(К—1)/2
Здесь введено обозначение Бп,к (р) = Б(р)Щ/К (р — 2к/К).
Таким образом, после фильтрации ошибка в проекциях от интересующей нас области, обусловленная случайным шумом, будет
0
где БП определяется равенством (21)
\
/ ОП Шр, (22)
—го
4. Результаты вычислительного эксперимента
Предложенные методы были исследованы в процессе вычислительного эксперимента. Математический фантом представлен на рис. 1, а. Радиус рассматриваемой области г0 равнялся 0.35. На рис. 1, а её граница помечена штриховой линией. Изображение фантома внутри данной области приведено на рис. 1, б. Расчёты производились на квадратной сетке 1025 х 1025 узлов. Ракурсы наблюдения были равномерно распределены по углу в интервале от 0 до п, их число N составляло 180. На каждом направлении было 1025 отсчётов. Значения случайной функции £<(р), моделирующей шум, в каждой точке проекции представляли собой нормально распределённую случайную величину с нулевым математическим ожиданием и дисперсией Д(р) = а/2|р|, где / — среднее значение проекционных данных в отсутствии шума.
Точность реконструкции оценивалась посредством нормированной среднеквадратичной ошибки А, которая рассчитывалась по формуле
A=,/I>j - 9i )2/Х (go,j )2, (23)
у j j
где g0j и gj — значения фантома и результата реконструкции в j-м узле соответственно. Суммирование ведётся по узлам, заключённым внутри заданной области.
Использовалась реализация гашр-фильтра, предложенная Шеппом и Логаном [17]. В настоящее время она широко применяется в различных приложениях томографии. Фильтр Шеппа — Логана имеет вид
^Lx п/2 - P^max Sin (P^ max) /а a
(п/2)2 - (p^max)2 ' ( )
где wmax — максимальная по модулю частота, ограничивающая ширину полосы пропускания. Для приведённых здесь результатов wmax = 0.7wN, где = п/h — циклическая частота Найквиста, h — шаг сетки, задаваемой на проекции. Такая величина wmax была выбрана по результатам серии реконструкций по данным, содержащим шум различного
Рис. 1. Математический фантом с отмеченной областью, представляющей предмет исследования (о); изображение фантома внутри этой области (б)
ОД
уровня. При меньших штах среднеквадратичная ошибка (23) была ниже, однако мелкая структура фантома на полученных томограммах оказалась сильно размытой. Следует отметить, что для разных значений штах рассматриваемые в работе зависимости в целом имели сходный характер.
Вычислительный эксперимент показал, что по мере повышения уровня шума ошибка реконструкции рассматриваемой области А, полученная обоими предлагаемыми методами, становится существенно ниже той, которую даёт непосредственное применение алгоритма Шеппа — Логана, реализующего формулу обращения двумерного преобразования Радона (2). При этом ошибка в проекционных данных в, определённая согласно (11), для алгоритма Шеппа — Логана растёт значительно быстрее. На рис. 2 представлены зависимости в и А от параметра а, определяющего дисперсию шума. Кривая 1 на рисунке соответствует классическому методу, кривая 2 — методу окна, используется функция Хэмминга (3) с ртах = 1.3го, кривая 3 — методу деления, число отрезков К равно 7, величина а в формуле для вычисления весовых коэффициентов (8) принималась такой, при которой ошибка А была минимальной.
150
0
120
90
60
30
'1
>
А. г—-*-**
г
а
0.2
0.4 0.6
5
0.8
1.6
1.2
0.8
0.4
2
3
а
0.2
0.4
0.6
0.8
Рис. 2. Зависимости ошибки В (а) и ошибки А (б) от параметра а, определяющего дисперсию шума О(р) = а/2|р|; 1 — классический метод, 2 — метод окна, 3 — метод деления
а
1.1
1.03
0.96
/1
У2
у* !
0 75
1.5
2.5
Рис. 3. Зависимость ошибки А от относительной ширины окна в = ртах/го, дисперсия шума соответствует а = 0.5; 1 — окно Хэмминга, 2 — окно Блэкмена, 3 — окно Парзена
0.82 0.808 0.796 0.784 0.772 0.76
1.05 0.99 0.93 0.87 0.81 0.75
1
д
Л
2
р г»-"-к*
1
II ■
ю
19 28
5
37
К
46
1У/
С 2
а
0.2
0.4
0.6
0.8
а
Рис. 4. Зависимость ошибки А от числа интервалов К: 1 — а = 0.4, 2 — а = 0.6 (а) и от параметра а: 1 — К = 7, 2 — К = 21 (б), дисперсия шума соответствует а = 0.5
Влияние формы окна иллюстрирует рис. 3. Кривые на рисунке — зависимости Д(в), где в = Ртах/го — отношение ширины окна к размеру области интереса. Проекционные данные содержат шум с дисперсией Б(р) = 0.5/2|р|. Видно, что все зависимости имеют по одному минимуму, значения Д в которых примерно одинаковы. Для окна Хэмминга минимум достигается при в = 1.3, для окна Блэкмена — при в = 1.6, для окна Пар-зена — при в = 1.8. С увеличением в ошибка для окна Блэкмена растёт быстрее. Это объясняется тем, что на границах его функция окна не падает до нуля.
Результаты исследования метода деления представлены на рис. 4. Шум в проекционных данных такой же, как и для рис. 3. На рис. 4, а приведена зависимость ошибки Д от числа интервалов К, на которые разбит носитель проекции. Кривые 1 и 2 соответствуют значениям а = 0.4 и а = 0.6 в формуле (8). Как следует из рисунка, при рассматриваемом уровне шума и выбранных а число К слабо влияет на точность реконструкции. Разброс ошибки Д для разных К не превышает 7% от её величины. При К < 20 графики имеют осциллирующий характер, что, по-видимому, обусловлено формой фантома. При К > 30 значения Д практически не изменяются.
На рис. 4, б показана зависимость ошибки Д от параметра а. Из приведённых графиков можно заключить, что с увеличением а от 0 до 0.5 ошибка незначительно уменьшается, а кривые проходят близко друг к другу. Дальнейший рост а приводит к росту Д,
Рис. 5. Результаты восстановления области, представляющей интерес для исследования: классический метод, А = 0.482 (а), метод окна, А = 0.673 (б), метод деления, А = 0.688 (в), метод Вайнберга (г). Дисперсия шума а = 0.05
а б
Рис. 6. Результаты восстановления области, представляющей интерес для исследования: классический метод, А = 1.910 (а); метод окна, А = 0.786 (б); метод деления, А = 0.784 (в); метод Вайнберга (г). Дисперсия шума а = 0.5
который составляет 25-27% на интервале а Е [0.6; 1.0]. При этом значения А несколько ниже для К = 21. В вычислительном эксперименте также получено, что для более высокого уровня шума минимум зависимости А(а) достигается при больших а.
На рис. 5, 6 даны изображения рассматриваемой области, полученные различными методами. Шум в проекционных данных определяется значением параметра а, равным соответственно 0.05 и 0.5. Таким образом, дисперсия шума в каждой точке проекции для рис. 6 в 10 раз больше, чем для рис. 5. Рисунки 5, а, 6, а соответствуют классическому методу, 5, б, 6, б — методу окна, при этом использовалось окно Хэмминга с в = 1.3, 5, в, 6, в — методу деления, К = 7, а = 0.1 (рис. 5, в) и а = 0.4 (рис. 6, в). Наконец, на рис. 5, г, 6, г приведены изображения, полученные методом Вайнберга. Для них одномерные проекции сглаживались фильтром скользящего среднего с шириной, равной пяти шагам сетки.
Заключение
В работе предложены два новых метода И,01-томографии по данным, искажённым нестационарным случайным шумом. Один из них заключается в выделении при помощи окна участков проекций, содержащих информацию об области, представляющей
предмет для исследования. В рамках другого метода проекции делятся на несколько частей, которые фильтруются по отдельности, после чего полученные функции суммируются с весами, зависящими от степени их зашумлённости.
Для разработанных методов, а также для алгоритма, реализующего формулу обращения (классический метод), получены оценки усиления шумовой составляющей. Проведённые вычисления показали, что в последнем случае средняя дисперсия шума в 1.6-1.8 раза больше, чем для метода окна, и в 2.2-2.7 раза больше, чем для метода деления. Установлено, что с усилением шума среднеквадратичная ошибка реконструкции для классического метода растёт быстрее. В частности, при максимальном уровне шума, который использовался в вычислительном эксперименте, она была больше примерно в 3 раза.
Список литературы
[1] Sahiner B., Yagle A.E. Region-of-interest tomography using exponential radial sampling // IEEE Trans. Image Process. 1995. Vol. 4, No. 8. P. 1120-1127.
[2] Sourbelle K., Lauritsch G., Tam K.C. et al. Performance evaluation of local ROI algorithms for exact ROI reconstruction in spiral cone-beam computed tomography // IEEE Trans. Nucl. Sci. 2001. Vol. 48, No. 3, pt II. P. 697-702.
[3] Lauritsch G., Sourbelle K., Tam K.C. Optimization of derivative kernels for exact cone-beam ROI reconstruction in spiral computed tomography // Ibid. 2002. Vol. 49, No. 3, Pt I, P. 728-732.
[4] Das P.C., Sastry Ch.S. Region-of-interest tomography using a composite Fourier-wavelet algorithm // Numer. Functional Analysis and Optimizat. 2002. Vol. 23, No. 7-8. P. 757-777.
[5] Zeng G.L., Gagnon D., Natterer F. et al. Local tomography property of residual minimization reconstruction with planar integral data // IEEE Trans. Nucl. Sci. 2003. Vol. 50, No. 5. P. 1590-1594.
[6] Хермен Г.Т. Восстановление изображений по проекциям. Основы реконструктивной томографии: Пер. с англ. М.: Мир, 1983. 349 c.
[7] Наттерер Ф. Математические аспекты компьютерной томографии: Пер. с англ. М.: Мир, 1990. 288 с.
[8] Вайнберг Э.И., Казак И.А., Курозаев В.П. Реконструкция внутренней пространственной структуры объектов по интегральным проекциям в реальном масштабе времени // Докл. АН СССР. 1981. Т. 257, № 1. С. 89-94.
[9] Kuchment P., Lancaster K., Mogilevskaya L. On local tomography // Inverse Problems. 1995. Vol. 11, No. 3. P. 571-589.
[10] Katsevich A.I., Ramm A.G. New methods for finding values of the jumps of a function from its local tomographic data // Ibid. 1995. Vol. 11, No. 5. P. 1005-1023.
[11] Аниконов Д.С., Ковтанюк А.Е., Прохоров Н.В. Использование уравнения переноса в томографии. М.: Логос, 2000. 224 c.
[12] Аниконов Д.С., Балакина Е.Ю. Полихроматический индикатор неоднородности неизвестной среды для задачи рентгеновской томографии // Сибирский математический журнал. 2012. Т. 53, № 4. С. 721-740.
[13] Василенко В.А. Сплайн-функции: Теория, алгоритмы, программы. Новосибирск: Наука, 1983. 214 c.
[14] Пикалов В.В., Преображенский Н.Г. Реконструктивная томография в газодинамике и физике плазмы. Новосибирск: Наука, 1987. 230 c.
[15] Хемминг Р.В. Цифровые фильтры: Пер. с англ. М.: Недра, 1987. 221 с.
[16] Bartlett M.S. Smoothing periodograms from time series with continuous spectra // Nature. 1948. Vol. 161. P. 686-687.
[17] Shepp L.A., Logan B.F. The Fourier reconstruction of a head section // IEEE Trans. Nucl. Sci. 1974. Vol. 21, No. 3. P. 21-43.
Поступила в 'редакцию 20 декабря 2013 г., с доработки — 4 марта 2014 г.
