CC BY
36
УДК 621.391:519.2
ВЕЙВЛЕТ - ТОМОГРАФИЯ В УСЛОВИЯХ ШУМА
И.А. Патрикеев, П.Г, Фрик (Пермь) Abstracts
A view on a Radon transform as on a singular wavelet transform allows to reconstruct only needed scales from the projections data with noise. The main properties of the algorithm are shown on axisymmetrical example for good obviousness only, however algorithm works out for arbitrary 2D function.
Двумерные вейвлеты. Семейство вейвлетов на плоскости И2 генерируется из материнской функции g путем применения двумерного сдвига Т и масштабирования Б, а в случае анизотропного вейвлета еще и поворота П [1]:
(Tbxbyg)
х х -bx
= g 1 СҐ 'С 1
(Dag)
X -2 а *х
.У. = а g -1
-а У.
(1)
(Пф8)
= g
coscp - sm<p x
|_апф coscp J| у
Вейвлет-преобразование определяется как свертка произвольной функции Дх,у) и вейвлета с параметрами Ьх,Ьу,а,<р.
Wf(bx, by ,а, ф) = JJ f(x,y)(Tbx’by DaQcpg)(x, y)dxdy.
(2)
R
Таким образом, двумерное вейвлет-преобразование отображает функцию Дх,у) в четырехмерное пространство. Для центрально-симметричной функции размерность вейвлет-пространства сокращается до трех:
Wf(b, а, ф) = JJ f(p)(TbDaQ<pg)(x, y)dxdy, (3)
где введены обозначения
R2
Р = ц/х2 +у2, b = /ьх2 +Ьу^. (4)
Преобразование Радона. Моделирование процесса получения проекционных данных традиционно строится на основе аппарата преобразования Радона, которое для пространства с размерностью п=2 представляет собой множество интегралов вдоль всевозможных секущих прямых s с параметрами р и ср:
ОС
Rf(p,(p) - Jf(pcoscp - ssin<p,psin(p + scoscp)ds. (5)
—00
Учитывая, что интеграл от f вдоль прямой равен интегралу по всей плоскости от произведения f на 5-функцию, сосредоточенную на этой прямой, предыдущее уравнение
преобразуется к виду:
Rf(p,<p)= jjf(x,y)5(x cosip + у sin ф - p)dxdy, (6)
R2
где аргументом 5-функции является левая часть уравнения секущей прямой
X COS ф + у БШф - р - 0. (7)
Поскольку любая прямая может быть получена из прямой х=0 путем сдвига и поворота, можно определить преобразование Радона в терминах вейвлет-преобразования [2]:
Яґ(р,ф)= |Jf(x,y)(TpQ(p5)(x)dxdy. (8)
н2
При этом ядро преобразования Радона рассматривается как сингулярный вейвлет. В отличие от регулярных вейвлетов среднее значение 5-функции не равно нулю. Кроме того, регулярные вейвлеты локализованы одновременно в пространстве сигналов и в частотной области, в то время как сингулярные вейвлеты сосредоточены в бесконечно малой области и, следовательно, имеют бесконечно широкий спектр.
Учитывая особенности 5-функции как сингулярного вейвлета, устанавливается связь между сингулярным вейвлет-преобразованием функции ґ и ее преобразованием Радона:
\\Т(Ьх,Ьу ,а,ф) = -Я^р'.ф), (9)
д.
где р ’ определяется по формуле, описывающей известную в томографии операцию «обратной проекции» [3],
р'= Ьх СОБф + Ьу БІПф (Ю)
Отметим, что с точностью до масштабного коэффициента вид сингулярного вейвлет-преобразования одинаков при любом а.
Для центрально-симметричной функции Радон-образ не зависит от (р и выражение (8) упрощается подстановкой х=р:
*^(Р)= | Г( р, У)с1у,
—00
а сингулярное вейвлет-преобразование вычисляется по формуле
(П)
\ЩЪ,а, ф) = — Ш’(Ъ собф). а
(12)
Для иллюстрации алгоритма в качестве теста рассматривается двумерная цен-трально-симмегричная функция
Г(х, у) - ехр
(д/х2 +у2 -Мо]
¿от
(13)
с параметрами Мо=0,5, сг=0,05.
Преобразование Радона данной тестовой функции вычислялось численно и дополнялось аддитивным белым шумом с параметром ш, задаваемым как отношение максимального значения амплшуды шума к среднему значению Радон-образа.
На рис. 1 показано сингулярное вейвлет-преобразование тестовой функции при а=1, вычисленное по ее Радон-образу без шума. По горизонтали - радиальная координата Ь, по вертикали - угловая координата ф (координата ф изменяется от 0 до п/2 вследствие угловой симметрии) Светлые участки соответствуют большим значениям, темные - меньшим.
-1 -0.5 О 0.5 1
ь
Рис. 1. Сингулярное вейвлет-преобразование тестовой функции
Вейвлет-синтез. Реконструкция функции по ее сингулярному вейвлет-образу осуществляется сверткой с синтезирующим вейвлетом, интегрированием по углу ф и по масштабу а. Выбор синтезирующего вейвлета Ь производится, исходя из условия [1,2].
ггё(г)Ку)^ ,2
И—ГГ?—dy = (27t) ,
R2
\у\
(14)
где ¿(у), Ь(у) - двумерные Фурье-спектры анализирующего и синтезирующего вейвлетов соответственно. Данное условие позволяет скомпенсировать сингулярность анализирующего вейвлета выбором подходящего регулярного вейвлета на этапе синтеза.
Условие (14) упрощается при подстановке в качестве анализирующего вейвлета 5-функции, сосредоточенной на прямой
00 л
г ау - _
] 777Ь(у е) = 2я, (15)
где ё = (соэф.япф), Ь - одномерный Фурье-спектр синтезирующего вейвлета, у - радиальная частота:
ОС
Ь(у)= |ь(р)е”1у,Мр. (16)
В результате свертки с h(pj по анализирующему вейвлет-образу Wf находится синтезирующий вейвлет-образ W'f:
27t со cog Л
W'f(b,a) = ]^ф j" Wf(p,a^)a-1hi---------- ----—Jdp. (17)
0 -00
На рис.2 показаны синтезирующие вейвлет-образы W’f тестовой функции при различных уровнях шума w, полученные сверткой с мексиканской шляпой h(p)=(l-p2)exp(-p2/2). Сравнивая синтезирующие вейвлет-образы тестовой функции (см. рис.2) можно отметить, что влияние шума проявляется в виде характерной структуры в области малых масштабов. Соответственно, при интегрировании по а необходимо найти оптимальную границу а,пш, отсекающую область высокочастотного шума
В общем случае алгоритм, обеспечивающий точное восстановление функции по проекционным данным без шума [2],
f(p)= lim fep (18)
г.->0,р->со
должен быть для обеспечения устойчивости в условиях шума преобразован к виду
f(p)= lim fe р, (19)
K_>,amin (W),P max
(W)
ь
Рис.2. Вейвлет-образ \¥' тестовой функции при различных уровнях шума \у где V/ - уровень шума, а е и р - пределы интегрирования в формуле
\ Рг ёа
(2°)
В последней формуле С- нормирующий коэффициент. Для мексиканской шляпы, использовавшейся в качестве синтезирующего вейвлета, С=2к.
Устойчивость реконструкции. Вследствие ограничения пределов формула (19) обеспечивает только приближенное восстановление Др) по тем масштабам, где влияние шума невелико. С уменьшением уровня шума пределы интегрирования расширяются и функция восстанавливается все более точно, стремясь к истинному значению.
Наглядно а„т, а™* можно определить по графику энергетического вейвлет-спектра, показывающему распределение энергии сигнала по масштабам:
ОО
ЕГ(а)= |^Т(Ь,а)|2сЗЬ. (21)
—00
На рис. 3 представлен энергетический спектр тестовой функции (13), вычисленный по 256 пространственным отсчетам для 16 масштабных коэффициентов. Е(п - график-спектра без шума, Б!?« - спектр функции, искаженный влиянием шума \у=20%, ЕДо -спектр с шумом 40% (по горизонтальной оси - а в логарифмическом масштабе).
Рис.З. Энергетический вейвлет-спектр тестовой функции
Видно, что шум искажаег спектр сигнала преимущественно в области малых а, поэтому для устойчивого восстановления тестовой функции достаточно выбрать нижнюю границу агат в соответствии с уровнем шума. Таким образом, обеспечение устойчивости в данном случае обеспечивается «фильтром больших масштабов»
Если требуется реконструкция ограниченного масштабного диапазона, то пределы интегрирования агат и атах выбираются таким образом, чтобы одновременно обеспечить устойчивость и подавить ненужные масштабы. В этом случае работает «полосовой фильтр». В общем случае для восстановления отдельных масштабов соответствующим образом выбираются пределы в формуле (20).
На рис.4 показаны результаты восстановления тестовой функции по 256 дискретным отсчетам Радон-образа с шумом при 64 масштабных коэффициентах, изменяющихся в пределах от атах=10 до атт=0,02 (верхний график) и атш: 0,01 (нижний график) (Г - исходная функция, ^ - результат восстановления без шума, Гзо - \у=20%, 1'40 -\у=40%).
Особенности поведения восстановленной функции вблизи точки р=0 объясняются когерентным усилением влияния шума (только для центральносимметричных функций) на масштабах порядка ширины минимального вейвлета. Дальнейшее уменьшении ап™ приводит к потере устойчивости даже при восстановлении тестовой функции по Радон-образу без шума (« 0%) вследствие дискретности проекционных данных
р
-1 -0.5 0 0-5 1
V
Рис.4. Реконструкция тестовой функции при различных значениях а™,,.
Выводы. Рассмотренный алгоритм вейвлет-томографии позволяет провести одновременно с восстановлением по проекционным данным фильтрацию шума и тем самым обеспечить устойчивость решения. При этом реконструкция только интересующих масштабов осуществляется выбором границ а,™ и a™*, причем в наиболее распространенном на практике случае высокочастотного шума достаточно определить только а^га. Общий алгоритм предполагает поиск весовой функции для каждого масштаба а в зависимости от распределения шума.
Библиографический список
1. Holschneider М. Wavelets: an analysis tool. Oxford Univ. Press, 1995, pp.423.
2. Holschneider M. Inverse Radon transform through inverse wavelet transform. Preprint P2364 CPT90, C.N.R.S. , Fev.1990.
3 Наттерер Ф. Математические аспекты компьютерной томографии. - М.: Мир, 1990. -288с.
