CC BY
37
Реконструкция томографических изображений с помощью пороговой обработки коэффициентов разложения проекционных данных по ортогональному базису
Ключевые слова: томография, вейвлеты, пороговая обработка, среднеквадратичный риск оценки, асимптотическая нормальность, преобразование Радона, устойчивый базис
Рассматривается задача реконструкции изображения по проекционным данным в модели с аддитивным шумом. Исследуются асимптотические свойства оценки риска при пороговой обработке коэффициентов разложения функции, описывающей проекционные данные, в ортогональный ряд- Приводятся условия, при которых имеет место асимптотическая нормальность оценки риска.
Кудрявцев A.A.,
МГУ им. М.В.Ломоносова, факультет ВМК, [email protected]
Шестаков О.В.,
МГУ им. М.В.Ломоносова, факультет ВМК;
Институт проблем информатики РАН, [email protected]
1. Введение
Томографические методы реконструкции изображений получили широкое распространение в самых разнообразных областях, включая медицину, биологию, физику плазмы, газовую динамику, геофизику, астрономию и радиолокацию (см., например, [1], [2] и [3]). Во многих видах
томографических экспериментов основную роль при реконструкции изображений играет преобразование (или оператор) Радона:
Rf(<p,s)= jf(x,y)dl =
V» (I)
00
= j/(.vcos^>-/sin q>,ssin <p + t cos(p)dt, s eR,<p e [0,2^),
где /(*,j/)eLr(R2) - функция с компактным носителем, описывающая искомое изображение, а ^ - прямая на
плоскости, по которой интегрируется функция f(x,y) (здесь s - расстояние от начала координат до прямой, а (р
- угол, образованный с осью X перпендикуляром, опущенным из начала координат на эту прямую). При фиксированном значении (р функцию Rf(tp,s) также называют проекцией функции f(x,y) под углом (р. Для преобразования Радона справедлива так называемая теорема о центральном сечении, играющая ключевую роль во многих методах реконсгрукции изображений по проекциям:
Rf(<p,(o) = yfl/г f Uocos(p,eos'm(p), (2)
где RJX<p,(o) - одномерное преобразование Фурье функции Rf(<p,s) по второй переменной, а /(«,,&),) - двумерное преобразование Фурье функции f(x,y) •
На практике проекционные данные регистрируются с
некоторой погрешностью (шумом), возникающей из-за несовершенства оборудования и различных случайных помех. И поскольку задача обращения преобразования Радона относится к классу так называемых некорректно поставленных задач, при реконструкции изображений приходится использовать методы регуляризации. В работах [4-7] рассматривается метод регуляризации, заключающийся в разложении функции изображения по двумерному ортонормированному вейвлет-базису с последующей пороговой обработкой коэффициентов разложения. При этом функция к/ раскладывается по так называемому вей глет-базису, вообще говоря не являющемуся ортонорми-рованным. Однако в силу «почти ортогональности» вейглетов среднеквадратичный риск метода оценивается при дополнительных предположениях, эквивалентных условию ортонормированности вейглет-базиса.
Методы анализа среднеквадратичного риска и его оценки в задаче реконструкции томографических изображений во многом опираются на тот факт, что преобразование Радона обладает многими свойствами линейных однородных операторов, которые рассматривались в работах [8-9]. В указанных работах рассматривается метод пороговой обработки коэффициентов вейглет-вейвлет разложения сигнала, и показывается, что при некоторых условиях регулярности оценка среднеквадратичного риска этого метода асимптотически нормальна. В данной работе мы представим метод обращения преобразования Радона, альтернативный методу вейвлет-вейглет разложения и исследуем свойства оценки среднеквадратичного риска этого метода.
Замечание 1 Везде в работе мы будем пользоваться следующими стандартными обозначениями:
математическое ожидание слу-
чайной величины ;
\А) индикатор события А;
Ф(х) стандартная нормальная функция
распределения;
ф(х) плотность стандартного нор-
мольного распределения;
</,£>
скалярное произведение функций /
слабая сходимость (сходимость по распределению).
2. Метод реконструкции и оценка риска.
Разложим функцию 1{/(<р,.ч) в ряд Фурье по переменной (р на [0,2л']- Поскольку функция /?/'(<-/>, .у) вещественна, мы будем использовать вещественный базис Фурье 1 £0%(п(р) 5Іп(«0>)|
[у/їтг УІЛ >/л
Для сокращения записи обозначим элементы этого базиса через \2а\'п_а, перенумеровав их соответствующим
образом. Имеем
«“О
2я
где Л/„(і)= |/?/'(<р,.9)2„(<ркЛр.
0
Рассмотрим вейвлет-базис {ш 4} 4с2, где
(//,,(*) = 2,п-1(/(2'х-к), а і//(д-) - некоторая вейвлет-функция. Индекс / называется масштабом, а индекс £ — сдвигом. Семейство {(/ А}у,б/ образует ортонормированный базис в Ь:(II) (см. [10]). Далее разложим функции Л/Дя) по базису
у.»} >.*6г •
/?/(«>,*) = X X ^Ч'іл)Ч',Л^гп{<р). (3)
»-0/,*б2
После применения к разложению (3) обратного преобразования К 1, получаем представление для функции /, которое является аналогом вейглет-вейвлет разложения, рассмотренного в работах [8], [9] и [11]:
Ах,у) = X X (4)
п°о)мг
где
, ч ЛТ^.К^Д') „ ц„.1г, ,ц
и^Лх’У)------------г---------. Д,„=||л [^Д
Ря.).к
Последовательность {//п . к} не образует
ортонормированную систему, однако если вейвлет-функция Ц/ удовлетворяет некоторым условиям гладкости, то последовательность {//„ ,} образует устойчивый базис, т.е. существуют такие константы 0 < А < В < оо, что
ЛХ Xе»./.* -
и-0 /.*• I 1"“®
п-0
(5)
для всех квадратично суммируемых последовательностей Иногда свойство (5) называют «почти
ортогональностью» (см. [7]). Функции иП1к по аналогии с
терминологией работ [5] и [11] назовем «вейглетами». Справедливо следующее утверждение, являющееся полным аналогом леммы 1 из работы [8].
Лемма 1 Пусть существуют такие константы А> 0. а> 1/2 и Ь>3/2. что
\ф((0^<,АЦ\\+\0^у{^а)Г1 (6)
для всех со є I*. тогда последовательность {и ,} образует устойчивый базис в Ь3^)-
Доказательство. Требуется показать, что существуют константы 0<Л<Й<оо такие, что для произвольной функции /■ є Ь3^) выполнено
М *Ш1 (7)
пш О >.*€2
(это соотношение эквивалентно условию (5)).
Определим оператор обратного проецирования, действующий на функцию %(<р,х):
2 я
К’я(х,у) = ^(<р,хсоь(р + у$\п(р)с!<р. ^
о
Оператор Я' является сопряженным к оператору Я (см.
[13]). Обозначим
“І.міХ’У) = Рп.)*К'[2»¥іЛх<У)-
Таким же образом, как в доказательстве леммы 1 из работы [8], можно убедиться, что (7) эквивалентно существованию таких положительных констант С, и С,, что
ІЕІ</.»иЛ’*с,и* <’>
п=О
И
<|0)
ч-оімг
Найдем константы Д (. Используя теорему о центральном сечении (2), имеем
РІІЛ =К’[^]|Г = ] '[2,у)к](х,у)\с1хс1у =
СО 00 _
= 11|^ 1[2пу/и\(ах,юг)|
| 2 ж» . і * 2
= — | ||Л х\2„уІХ ](л>соя<^),л>зіп ^)| \cc\doxltp = — Ц</;.д(<»)|'|<а|<Лг> =
4л
||^(®)|3|<4Л» = 2У||£||\
где через £ обозначена функция, для которой
2 4л
(П)
Таким образом, константы Д к не зависят от И и к и пропорциональны 2)П~. Докажем теперь (9). Имеем
£Д/> = X X \</'ип.м)\2 = X X 11Л*>У>«*.1Лх,У)1Ьф
л-0 ІМ7.
=хх
11/(л»,,йл)йлуд(<»,,«, )с/(о^(о2
=11
/г” О уДеЛ
■.я
| |/'(<У СОв^, <У БШ <0)м„ м («сое<р, ы$\п<рУос1(ос1<р
1
тХХ
(—^) я=0 у.Ас/
I
I /л/(«»,й>)Лмяуд(^
тХХт^г
(2л-) „ о / 4Д,^
ХЛЪ~А4г
(2л)" „ о,.*,/4Д"(
’^ЁХтз-
| |«л «9, «)г„ (<р)Ч>1к( (0)\(0\с/мЩ
о -®
2
90
р/„(<У)^уД®)Н</<У =
|л/,(ю)2 ‘аф(2 '(0)еМ2 \а^а
2||^|г гг,*/ л-2 =—УУ—
2)|^|2 ^,‘Г/л2/‘
|л/„(й)^(2^)Н'я-еыг/^
лгО/>с/
г7>+|
| * У^У„(&>+л2''1 м)|а>+л2''1 ^(2 'ео+2лп)с/а
Функция
Хл/,(<у + /г2у*'л;)|<у + л2у+|/и| ~ £(2~'ю + 2тн)
ме2
периодична с периодом л2/и. Следовательно,
. » >2^' _ _____________________________________________
и(Л=-туУУ I 1&(®+*2*«)|»+*2'*,|И| -«Г>®+2».)
2]|*| Я=0 ;сг О то. г
. _ Л" _______________
(1(0 =
“Пы^ХХ \ 1&(® + *2*,«)|© + *2*,Я|| "‘=(2 /<у+ 2л7и)х
7НЯ1 н=п ;.-7 ...7
^ю\ п=оуе2 о
хХ^>(, + -т2‘М/)Ь + /г2У*'/| ь(2 1 (о+ 2л1)с1 (о =
=—:ХХ \hfS(0\(^ &М0+л-*-' 'т\<0+^2’'£(2 'ы^(2 1 (0+2тп)с1(1).
ЦЩ ,го/л«г_»
Далее, поступая, как в доказательстве леммы 1 из работы [8], получаем, что при выполнении условия (6) для некоторой положительной константы с, выполнено
(У(Я<Хс,Ы2,
я*0
где gn((0)= к/п((0)|<и|1/2. Обозначим £(^,ю) = Д/(р,ю)|<у|,/2- По теореме о центральном сечении
2л оо 2я ев
И = I /|я(«».<»)|2Жос!<р= | ||я/(<р,«)| \a\doxlip =
О -<о 0 -се
2ж« , во оо ^
= 4л ||/(<ус°5<0,ю5т<0)| (А1м1<р = 4л | ||/(х,^)|2<£г<ф’=4л|/| =4л|/||\
О 0 -00-40
Кроме ТОГО,
\\g\f = | {НХя/ < ^ )Хя/„ (® )£,(? =
о I-О
= 1ХИН>^=Х1кГ-
^’(/)=£ х |</.«;.„>Г=ХХ|</.Д,,^[^,.]>Г =
я-0уд«2 /»=0 у.*
=хх
я=0 УД€/
-II
я-о )мг
00
-XX
Я“0 уДе?
=ХХИ/’^2»п*>Г =
я*Ч>у.Ае2
2ж оо ___________
№/: 1 |л/(«м)2»(р)^(5)*^
О -оо
2ж оо ______
||ь ||2г~ | | Я/(<?, <у )2„ (</> )у/ ,* ( <у )с/л*/<р
О -00
2» 00 _____________________
\ф.п 11к/((р,М)\(^ ~2п(<р)ц//Л((0)\с>\ 1"(1сос1(р
О -«о
|#|2Л }^Д»)Н1в<М®Н1в</®
=хх
я=* 0 уДс2
“II
я “О у >«2
И2>^ р/л(«)Н,я2
-«о
« _______________
|ф-^ ГЛ/„(«)^(2-'«)|«|1/2е'“*2 '</«
Следовательно, (9) выполнено. Докажем теперь (10).
=1X1
я**0 )Х^7.
Здесь через 7] обозначена функция, для которой {]{(!>) = ц/((о)\(о\ 12 • Далее, поступая, как при доказательстве (9), получаем, что найдется положительная константа с, такая, что
(/(/) £С2||/||2.
Лемма доказана.
Далее будем предполагать, что функция /, описывающая изображения, имеет компактный носитель и равномерно регулярна по Липшицу с некоторым показателем у >0 - Оказывается, что при этом функция Я/ равномерно регулярна по Липшицу с показателем / + \/2 (см. [13] и [14]). Также будем считать, что вейвлет-функция I// удовлетворяет условиям теоремы 6.3 из [10] и леммы 1. Таким образом, найдется такая константа А > 0, что
<л/,г„^,><^ (12)
при всех неотрицательных уег и всех к е2.
При практической реализации разложения (3) дискретные коэффициенты получаются умножением матрицы значений функции /{/' на ортогональную матрицу /■" и ортогональную матрицу IV, определяемые базисом Фурье ^а} и вейвлет-базисом \ц/ к \ • Причем, поскольку в
проекциях присутствует шум, рассматривается следующая модель измеряемых данных:
« = 1..................2',у = 1 2\ (13)
где J - некоторое положительное число («разрешение» изображения), х ~ наблюдаемые данные, Я -
преобразование Радона, / — истинная (незашумленная) функция изображения, а е ,, - независимые случайные погрешности, которые имеют одинаковое нормальное
распределение с нулевым средним и дисперсией а2. Таким образом, в силу ортогональности матриц преобразования, дискретные коэффициенты разложения наблюдаемых данных (13), которые мы обозначим через Уп)к*
описываются следующей моделью:
У.Л = К.и + «&. " = 0,...,2У-1,7 = 0,...,7 -1,* = 0,...,2; -1,
(14)
где а случайные погрешности
также независимы и имеют такое же распределение, как и
Поскольку в измерениях присутствует шум, необходимо использовать некоторые процедуры для его удаления. Мы рассмотрим процедуру пороговой обработки (см. [10]). Смысл пороговой обработки вейвлет-коэффициентов измеряемого сигнала заключается в удалении достаточно маленьких коэффициентов, которые считаются шумом. Будем использовать так называемую мягкую пороговую обработку с порогом т, зависящим от масштаба j. Для
подавления шума к каждому коэффициенту Уп . к применяется функция р1 (х), определяемая соотношением
рт (дг)=«я«(х)(|дг|-7';).
При такой пороговой обработке коэффициенты, которые по модулю меньше порога 7обнуляются, а абсолютные
величины остальных коэффициентов уменьшаются на величину порога. В результате функция сигнала / в разложении (4) (точнее, ее нормированная версия) оценивается следующим образом:
/=ХХХАм.*Рг/П.;.>,м-
/1-0 Jm0 кш0
Среднеквадратичная ошибка (риск) оценки / определяется как
г„ = е||/->
Так как система функций {и„1к} не является
ортонормированной, риск нельзя представить в виде суммы вкладов погрешностей отдельных коэффициентов, однако в силу устойчивости базиса {// (} найдутся такие константы
А и В, зависящие от I//, что А г, < га < Вг,, где
2 -и—12^-1
О =о(/,<т.Т)= ХХХДм.*Е<А-.у.»-Рт,(У..1лУ)\ (15) п-0 у-о»-о
(здесь Т = (7^,...,Г;|) - вектор порогов). Таким образом, с точностью до множителя величина гд асимптотически ведет себя так же, как , поэтому в дальнейшем мы будем рассматривать асимптотическое поведение Г,.
В выражении (15) присутствуют неизвестные величины ця , поэтому вычислить значение г3 нельзя. Однако его можно оценить. В каждом слагаемом если
\У„,к\>Т,< то вклад этого слагаемого в риск составляет р2п к(а~ + Т2), а если |уп т0 вклад составляет
fijj.Ml.jj, ■ Поскольку Y.Y2hk=a2+nljk, значение можно оценить разностью у2 , - а2 •
Таким образом, в качестве оценки риска можно использовать следующую величину:
о (/,ст. т)=х'ххк,^[С>-^7'у]-
п=0 j=0t=0
FIг, ajt] = (x-a- )l(|.t| <Т;) + (а2 + Tj )1(|дг| > Г;), (16) причем эта оценка является несмещенной (см. [ 10]).
3. Предельные теоремы для оценки риска.
По аналогии с работами [6], [8] и [9] мы рассмотрим ситуацию, когда для каждою масштаба у выбирается «универсальный» порог. Такой порог определяется соотношением 7)., = (Tyj2\n2l*J , поскольку на каждом
масштабе количество коэффициентов равно 2'*J, и получил название «универсальный», так как он не зависит от наблюдаемых данных и при выборе этого порога риск Kj
близок к минимальному (см. [12]).
Дтя оценки риска (16) справедлива следующая теорема, являющаяся аналогом теоремы 1 из [8].
Теорема 1 Пусть вейвлет-функция I// удовлетворяет условиям теоремы 6.3 из [10] и леммы I, а функция / е L‘(R) имеет компактный носитель и равномерно регулярна по Липшицу с показателем у > 0. тогда
■ ф(л-) при J-> оо, О 7)
где [)1 = ор2||^||:л/2/7 2?‘>, функция £ задана соотношением (11). а Ти=(Тил.........7^.,).
Зачастую дисперсия шума а2 не известна и ее также необходимо оценивать с помощью некоторой оценки <т2. При этом выражение (16) принимает вид
г,(/,«т. т) = ХХХД:.л*и^^.т-,].
п-0 у-0 к-0
(18)
и вместо порогов Tv i используются пороги
7) . = а^2\п2'*'' (обозначим вектор этих порогов через т, )•
Обычно дисперсия а2 оценивается по выборке снгнача (или изображения), однако ее можно оценить и по независимой выборке. Для этого следует произвести измерение пустого сигнала, тогда наблюдения будут представлять из себя чистый шум, по которому и оценивается а2 ■ Далее мы рассмотрим оба случая.
В случае, когда б2 строится по независимой выборке, для оценки риска (18) справедливо следующее утверждение.
Теорема 2 Пусть вейвлет-функция 4/ удовлетворяет условиям теоремы 6.3 из [10] и леммы /, а функция / е Ь2(1?) имеет компактный носи теп ь и равномерно регулярна по Липшицу с показателем у>0. Пусть <х2 - не зависящая от Уп (, оценка дисперсии <т2. для которой выполнено
р(2‘/(ст2 - <у2) < *)=> Фу(.х) при J-* оо, (19)
где ф (х) - функция распределения нормального закона с
нулевым средним и дисперсией . Тогда
(20)
(дг) при J
где ф .(х) — функция распредепения нормального закона с нулевым средним и дисперсией
18 а*
Если дисперсия оценивается по коэффициентам разложения у к и функция f удовлетворяет требуемым
условиям регулярности, то в силу (12) ее оценивают по половине всех коэффициентов на уровне j = J-1, так как эти коэффициенты фактически содержат только шум. Мы рассмотрим популярные оценки неизвестной дисперсии: <т2
— выборочная дисперсия, а„ - интерквартильный размах и ст., - абсолютное медианное отклонение от медианы. Эти
'V
оценки определяются следующими соотношениями:
&1 =-
-12 -1 I 2 —12*^— —I
I ZC-U- 1^-М
-0 к-0
пт 0 4*0
(21)
Y -У
- _ 1 2-/—1.(3/4) 2-/—1.( I /4)
СГ„-----------------------------------,
med I YnJ.l t -
med YeJ.u I
(22)
(23)
где К2 /_| (|11 и Ку м - выборочные квантили порядка 1/4 и 3/4, построенные по набору коэффициентов У„/.ц, О</7<2' -1, 0<к <2' 1 -1, £ - теоретическая квантиль
порядка 3/4 стандартного нормального распределения (4 = 0.6745), а обозначает выборочную медиану.
Теорема 3 Пусть вейвлет-функция у/ удовлетворяет условиям теоремы 6.3 из [10] и леммы 1, а функция / е Ь2(К) имеет компактный носи тень и равномерно регулярна по Липшицу с показатепем у > М2- Пусть оценка дисперсии <т~ определяется выражением (21). Тогда
р(/,ст. Т,.)-о(/,<т,Т,.) Dj
< х => Фт(-*) пР" J °°*
(24)
где Ф г (х) - функция распределения нормального закона с нулевым средним и дисперсией Т2 = 2/9.
Теорема 4 Пусть вейвлет-функция у/ удовлетворяет условиям теоремы 6.3 из [10] и леммы 1, а функция /є ) имеет компактный носитель и равномерно регулярна по Липшицу с показателем у>\. Пусть оценка
А 2
дисперсии О определяется выражением (22) или (23). Тогда
<х =>Фг(.г) при ./->00,
г, (/.<*, Т,.)-г,(/>.Т,.)
О.!
где Ф, (.г) — функция распределения нормального закона с нулевым средним и дисперсией
г;_ 7 4
Зб£,МСш))2 3
Доказательство теорем 1—4 полностью аналогично доказательству соответствующих теорем из работ [8, 9].
Замечание 2 Вектор порогов т, фактически используется для построения оценки функции /?/ . В работе [11] вместо вектора порогов Т, рассматривается единый для всех масштабов порог (в данном случае он равен Тк = а^41п22'1А который также назван «универсальным», и показывается, что при построении оценки функции f этот порог является асимптотически близким к оптимальному в минимаксном смысле. При рассмотрении порога Тк все результаты данной работы останутся неизменными.
Работа поддержана министерством образования и науки РФ (государственный контракт №П779).
Литература
1. Тихонов А.Н., Арсенин В.Я., Тимонов А.А. Математические задачи компьютерной томографии. - М.: Наука, 1987.
2. Троицкий И.Н. Статистическая теория томографии. - М.: Радио и связь, 1989.
3. Как А.С., Slaney М. Principles of Computerized Tomographic Imaging. - NY: IEEE Press, 1988.
4. Маркин A.B., Шесгаков O.B. Асимптотики оценки риска при пороговой обработке веивлет-вейглет коэффициентов в задаче томографии // Информатика и ее применения, 2010. Т.4. №2. С.36-45.
5. Donoho D. Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition // Applied and Computational Harmonic Analysis, 1995. Vol.2. P. 101-126.
6. Kolaczyk E.D. A wavelet shrinkage approach to tomographic image reconstruction // J. Am. Statist. Assoc., 1996. Vol.9l. P. 10791090.
7. Lee N. Wavelet-vaguelette decompositions and homogenous equations: PhD dissertation. Purdue University, 1997.
8. Кудрявцев A.A., Шестаков O.B. Асимптотика оценки риска при вейглет-вейвлет разложении наблюдаемого сигнала// T-Comm - Телекоммуникации и Транспорт, 2011. № 2. С.54—57.
