CC BY
18
УДК 517.444, 519.22 О.В. Шестаков1
РЕГУЛЯРИЗАЦИЯ МЕТОДА РЕКОНСТРУКЦИИ ФУНКЦИИ ПО ЕЕ СФЕРИЧЕСКОМУ ПРЕОБРАЗОВАНИЮ РАДОНА*
В работе рассматривается задача реконструкции томографических изображений по проекционным данным в условиях математической модели, использующей сферическое преобразование Радона. Приводятся оценки точности реконструкции при использовании конечного числа сферических проекций.
Ключевые слова: сферическое преобразование Радона, томографические изображения, оценки точности.
1. Введение. Задача определения функции по ее сферическому преобразованию Радона возникает как в прикладной, так и в чисто теоретической области математики. За последнее время интерес к этой задаче возрос в связи с появлением новых методов медицинской диагностики, использующих термоакустичекую и фотоакустическую томографию (см. [1, 2]). В одной из первых работ, посвященных этой задаче, применяется разложение функции в ряд Фурье, коэффициенты которого находятся с помощью обращения преобразования Ганкеля (см. [3]). В работе [4] получена формула обращения, очень похожая на формулу обращения классического преобразования Радона.
На практике сферическое преобразование Радона известно лишь для конечного числа центров окружностей (сфер), и в этом случае задача обращения уже не имеет единственного решения даже в случае классического преобразования Радона, что приводит к так называемому парадоксу компьютерной томографии (см. [5]). Этот парадокс рассматривался в работах [5, 6], и были получены оценки близости между функциями, имеющими одинаковые или близкие проекции по конечному числу заданных направлений. В этой работе мы получим подобные оценки в случае использования сферического преобразования Радона.
2. Сферическое преобразование Радона. Пусть /(ж) — бесконечно дифференцируемая функция с компактным носителем. Далее, не ограничивая общности, будем считать, что носителем функции /(ж) является круг II единичного радиуса с центром в начале координат. Определим сферическое преобразование Радона на плоскости следующим образом:
= ^ / Др + Гв) йв. (1)
Здесь 51 — множество направлений, задаваемых единичными векторами в К2 с центром в начале координат, р € Б1 (т.е. р = (сов ср, втер), ср € [0,2ж)), а г е (0, сю). По аналогии с классическим преобразованием Радона будем называть интегральные преобразования вида (1) сферическими проекциями, т. е. сферическая проекция представляет собой функцию от г € (0, сю) при фиксированном р.
В работе [4] получены следующие формулы обращения.
Теорема 1. Пусть /(ж) € Со°(112) имеет носитель в круге 11. Тогда для х € I/
2
12
f(x) = 77-Да 2тг
rRf(p, г) log
■р\'
drdp,
(2)
или
'<*> = 2п
[r(Rfyr(p,r)yr log
drdp.
(3)
Здесь Ах — оператор Лапласа.
1 Факультет ВМК МГУ, доц., к.ф.-м.н.; ИПИ РАН, ст. науч. сотр., e-mail: oshestakovQcs.msu.su
* Работа выполнена при финансовой поддержке РФФИ (проекты 11-01-00515а и 11-01-12-26-офи-м), а также Министерства образования и науки РФ (ГК № 14.740.11.0996).
Функция Д/(р, г) естественным образом доопределяется для отрицательных г как четная функция. Однако если вместо этого доопределить ее нечетным образом, тогда из формул (2) и (3) можно получить следующую формулу обращения.
Следствие. Пусть выполнены условия теоремы 1, и по второй переменной Д/(р, г) доопределена как нечетная функция при отрицательных значениях г, тогда для х € 17
2
(4)
51 -2
Внутренний интеграл понимается в смысле главного значения.
Формула (4) очень похожа на стандартную формулу обращения классического преобразования Радона (см., например, [7]). Для дальнейшего перепишем ее в несколько ином виде. Поскольку функция /(ж) имеет носителем круг II, пределы внутреннего интегрирования можно продолжить от ^оо до оо. Таким образом,
/(ж) = ^ J \х-р\ р, \х-р\) ф,
где
ос
ОСр,-) = -
7г j 5 — г
Функция (¿(р, в) равна преобразованию Гильберта от функции г) по второй переменной. Пре-
образование Гильберта для функции д(г) определяется как
оо
Я„(.)=1
IX ,} 5 — Г
Если перейти к преобразованиям Фурье вида
оо
Ш = А= [ е-™д{г)йг,
—оо
то (см. [8])
я,(ш) = (5)
Воспользовавшись формулой (5) и свойством преобразования Фурье от производной, имеем
Я{р,ш) = И Е/(р,ш). Следовательно, формулу (4) можно переписать в виде
оо
№ = 1¡\х-р\-±= I е^-^^ЩХр^йшйр. (6)
—оо
3. Оценки точности реконструкции по конечному числу сферических проекций. Задача восстановления функции /(ж) по ее сферическим проекциям является некорректно поставленной, и, кроме того, на практике сферические проекции регистрируются только для конечного числа точек р, поэтому вместо точной формулы (6) нужно использовать ее регуляризованный вариант, в котором преобразование Гильберта Н заменяется на преобразование II,,. для которого
г
где множитель Wa(w) играет роль регуляризирующего окна. В результате вместо точной функции /(ж) получается приближенная функция /о-(ж), но метод реконструкции становится устойчивым к погрешностям. В дальнейшем мы будем использовать регуляризирующее окно вида
Wa(w) = e-aj2ij2/2.
Пусть Fu — класс всех неотрицательных функций /(ж) G Cq°(R2) с носителем в единичном круге U, таких, что
sup ||grad/^)|| ^ С. xeR2
Для регуляризованной формулы обращения в классе Fjj можно получить оценку точности реконструкции функции по конечному числу сферических проекций. Справедлива следующая теорема.
Теорема 2. Пусть N = 2п, где п — натуральное число, и Pi,-..,PN выбраны следуют,им образом:
Рз = {Vj, -1 )!(у) + I)1/2, j = 1,... ,П,
Рз = (1, vj-n)/(vj_n + 1)1/2, j = п + 1,..., 2п,
где
щ = соз(тт(2к — 1)/(2п)), к = 1,...,п.
Если f(x),g(x) G Fjj и сферические проекции функций /(ж) и д(ж) (доопределенные нечетным образом для отрицательных г согласно следствию из теоремы 1) совпадают для р, равного pj и —pj, j = 1,..., N, т. е.
RÎ(±Pj,r) = Rg{±P:h.r)-. J 1.....
тогда
supl/Дж) -g„(x)\ < (4^3+ ) —т. (7)
\J<TV J Я<7\,Л' / 1 ^ v " 1 о I о '
xeu \ 3 J na2
Доказательство. Используя регуляризованный вариант формулы (6), имеем
\U(x) -g<j{x)| =
1
2у/2тг
s1
—сю
сю
! \ш\е~ш'Л<Т'Л\к{р,ш)\ ёшёр, (8)
51 -оо
где ]г(р,ш) = К/(р,ш) — Кд(р,ш). Оценим \]г(р,ш)\ = \]г(р((р),ш)\. Для всех (р € [<£>з — ж/(2п), (р^ и (р € [<Рз, фз + 7г/(2п)], ;/' I...../V. используя формулу конечных приращений, имеем
\h(p((p),w)\ = \Цр((р),ш)\ - \Цр((рз),ш)\ < |Цр((р),ш) - Цр((рз),ш)| <
7г
< sup \к^(р(ср),ш)\\ср - cp:j\ < sup |/^,(p(¥>),w)| —.
^е[0,2тг) ^е[0,2тг) /п
Далее, учитывая тот факт, что функции /(ж) и д(ж) имеют носители в круге U, получаем
сю
-ос S1
2
1 Г Г Rf^
-2 s1
Следовательно, для всех р £ 51 имеем
2у/2тгС < —-—.
Подставляя эту оценку в (8), получаем
oc
С í Í 2 2 2С Í
\fa(x) - да(х)\ 5$ — / \х -р\ / |ш| е~ш а /2 duidp < ^ / |ж dp. (9)
S1 -оо 51
Оценим последний интеграл. Воспользуемся тем, что р = р(<р) = (cos (р, sin (р) и что подынтегральная функция имеет период 2я\ Имеем
2тг
J \x — p\dp = J \j\x\2 — 2 (x,p) + 1 dp = J \J I ж 12 — 2\x \ cos (p + 1 dip.
Подынтегральная функция зависит от х только через |ж|. При этом, поскольку х € U, |ж| ^ 1. При |ж| ^ 1 максимум выражения \х\2 — 2 |ж| cos р +1 достигается в точке |ж| = 0, если cos ^ 1/2,ивточке |ж| = 1, если cos < 1/2. Таким образом, для всех х € U
i 2-7г
2тг
2-7г з 3 27г
J \j\x\2 — 2 |ж| cos + 1 dp ^ J — 1 cos pdp + J dp + J dp = 4л/3 •
3
Подставляя эту оценку в (9), получаем (7). Теорема доказана.
На практике в силу несовершенства оборудования и наличия шума сферические проекции регистрируются с некоторой погрешностью. Если предположить, что погрешность не превосходит заданного уровня е, то можно получить оценку точности реконструкции с учетом этой погрешности. Справедлива следующая теорема.
Теорема 3. Пусть сферические проекции функций f(x),g(x) € Fy (доопределенные нечетным образом для отрицательных г согласно следствию из теоремы 1) отличаются не более чем на е (е > 0) для р, равного p:j и —pj, j = 1, • • •, N, где p:j те же, что и в предыдущей теореме, т. е.
sup |Rf(±pj,r) - Rg{±pj,r)j < e, j = 1,...,N, re R.
тогда
(2Ж \ / 2С \ 4^3+ —) I —^ + —j). (10)
3 J \7Ta2 no2 J
Доказательство. Поступая так же, как в теореме 2 при оценке \к{р,ш)\ = \1г(р((р),ш)\, для всех (р € [<fj — 7г/(2п), cpj] и (р € [<pj, <pj + ж/(2п)}, j = 1,..., N, имеем
\h(p(<p),u)\ = \h(p(<p),u)\ - \h(p(<pj),u)\ + \h(p(<pj),u)\ <
Далее
\h(p((pj),w)\ =
00
1
< Ih(p(v),w) - h(p((pj),u}) \ + \h(p((pj),u})\ < —-— + \h(p((pj),u})\ .
е~шг(Rf(p((pj),r) - Rg(p(<pj),r)) dr <
2
If 4e
< -= / \Rf(p(<pj),r) - Rg(p(<pj),r)I dr < — V2it J л/2п
-2
4e 2V2жС
-2
Следовательно, для всех р G S1
\1г(р,ш)\ <
67Г И
Подставляя эту оценку в (8) и действуя, как в теореме 2, получаем (10). Теорема доказана.
Правая часть в оценке (7) из теоремы 2 с ростом п убывает со скоростью 0(1/п), а правая часть в оценке (10) из теоремы 3 с ростом п и уменьшением е убывает со скоростью 0(е) + 0(1/п). Это означает, что использование регуляризованной формулы обращения (6) приводит к устойчивому методу реконструкции функции по конечному числу сферических проекций.
СПИСОК ЛИТЕРАТУРЫ
1. Wang X.D., Pang G., Ku Y. J., Xie X.Y., Stoica G., Wang L.-H. V. Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain // Nature Biotechnology. 2003. 21. P. 803-806.
2. Kruger R. A., Miller K.D., Reynolds H.E., Kiser W.L., Reinecke D.R., Kruger G. A. Breast cancer in vivo: contrast enhancement with thermoacoustic ct at 434 mhz-feasibility study // Radiology. 2000. 216. P. 279-283.
3. Norton S. J. Reconstruction of a two-dimensional reflecting medium over a circular domain: Exact solution // J. Acoust. Soc. Amer. 1980. 67. P. 1266-1273.
4. Finch D., Haltmeier M. Inversion of spherical means and the wave equation in even dimensions // SIAM. J. Appl. Math. 2007. 68. N 2. P. 392-412.
5. Khalfin L. A., Klebanov L.B. A solution of the computer tomography paradox and estimating the distances between the densities of measures with the same marginals // The Annals of Probability. 1994. 22. N 4. P. 22352241.
6. Шестаков О. В., Савенков Т. Ю. Оценка расстояния между плотностями вероятностных мер, имеющих близкие проекции // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2001. № 4. С. 44-46.
7. Хелгасон С. Преобразование Радона. М.: Мир, 1983.
8. Natterer F. Inversion of the attenuated Radon transform // Inverse Problems. 2001. 17. P. 113-119.
Поступила в редакцию 30.05.11
52
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2012. № 1
REGULARIZING RECONSTRUCTION OF FUNCTION FROM ITS SPHERICAL RADON TRANSFORM
Shestakov O. V.
This paper deals with the problem of reconstructing tomographic images in the frames of mathematical model based on spherical Radon transform. Some accuracy estimates for reconstruction when using finite number of spherical projections are obtained.
Keywords: spherical Radon transform, tomographic images, accuracy estimates.
