Вейвлет-томография в условиях шума

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.

Двумерные вейвлеты. Семейство вейвлетов на плоскости R² генерируется из материнской функции g путем применения двумерного сдвига Т и масштабирования D, а в случае анизотропного вейвлета еще и поворота Q [1]:

№)

СО8ф sin ср

- 51ПфТх cos о J у

Вейвлет-преобразование определяется как свертка произвольной функции Дх,у) и вейвлета с параметрами b_x,b_y,a,

Wf(bx, Ьу ,а,ф) = ff f(x,y)(Tb’!’byDaQtpg)(x, y)dxdy.

R2

Таким образом, двумерное вейвлет-преобразование отображает функцию Дх,у) в четырехмерное пространство. Для центрально-симметричной функции f(p) размерность вейвлет-пространства сокращается до трех:

Wf(b, а, ф) = JJ fXpXThDWgXx, y)dxdy,                       (3)

R² где введены обозначения

Р = Vx2 +У2-

b = ^bx2 + by2.

Преобразование Радона. Моделирование процесса получения проекционных данных традиционно строится на основе аппарата преобразования Радона, которое для пространства с размерностью п=2 представляет собой множество интегралов вдоль всевозможных секущих прямых 5 с параметрами р и ф:

Rf(p,

—00

Учитывая, что интеграл от f вдоль прямой равен интегралу по всей плоскости от произведения f на 5-функцию, сосредоточенную на этой прямой, предыдущее уравнение преобразуется к виду:

КЦр,ф) = JJ f(x, у)5(хсобф + у зшф - p)dxdy,(6)

R2

где аргументом 5-функции является левая часть уравнения секущей прямой х cos ф + у 5Щф - р - 0.(7)

Поскольку любая прямая может быть получена из прямой х=0 путем сдвига и поворота, можно определить преобразование Радона в терминах вейвлет-преобразования [2]:

КГ(р,ф)= JJ f(x,y)(TpQ(p5)(x)dxdy.(8)

R2

При этом ядро преобразования Радона рассматривается как сингулярный вейвлет-. В отличие от регулярных вейвлетов среднее значение 5-функции не равно нулю. Кроме того, регулярные вейвлеты локализованы одновременно в пространстве сигналов и в частотной области, в то время как сингулярные вейвлеты сосредоточены в бесконечно малой области и, следовательно, имеют бесконечно широкий спектр.

Учитывая особенности 5-функции как сингулярного вейвлета, устанавливается связь между сингулярным вейвлет-преобразованием функции f и ее преобразованием Радона:

Wf(bx,by ,а,ф) = ^КГ(р',ф),                                     (9)

где р ’ определяется по формуле, описывающей известную в томографии операцию «обратной проекции» [3], р'= Ьх СО8ф + Ьу 8Щф                                          (10)

Отметим, что с точностью до масштабного коэффициента вид сингулярного вейвлет-преобразования одинаков при любом а.

Для центрально-симметричной функции Радон-образ не зависит от ф и выражение (8) упрощается подстановкой х=р:

Rf(p)= ff(PO')dy,                                            (И)

а сингулярное вейвлет-преобразование вычисляется по формуле

Wf(b,a,

Для иллюстрации алгоритма в качестве теста рассматривается двумерная центрально-симметричная функция

f(x,y) -exp<- -

х2+у2-

(В)

с параметрами Мо=0,5, о=0,05.

Преобразование Радона данной тестовой функции вычислялось численно и дополнялось аддитивным белым шумом с параметром w, задаваемым как отношение максимального значения амплитуды шума к среднему значению Радон-образа.

На рис.1 показано сингулярное вейвлет-преобразование тестовой функции при а=1, вычисленное по ее Радон-образу без шума. По горизонтали - радиальная координата Ь, по вертикали - угловая координата (р (координата ср изменяется от 0 до л/2 вследствие угловой симметрии) Светлые участки соответствуют большим значениям, темные - меньшим.

-I                       -0.5                        0                        0.5                         1

Рис. 1. Сингулярное вейвлет-преобразование тестовой функции

Вейвлет-синтез. Реконструкция функции по ее сингулярному вейвлет-образу осуществляется сверткой с синтезирующим вейвлетом, интегрированием по углу <р и по масштабу а. Выбор синтезирующего вейвлета h производится, исходя из условия [1,2].

ГГё<У)Ь(у) .. .2

JJ—~2—бу = (2я) ,                                      (14)

R2

где g(y),h(y)- двумерные Фурье-спектры анализирующего и синтезирующего вейвлетов соответственно. Данное условие позволяет скомпенсировать сингулярность анализирующего вейвлета выбором подходящего регулярного вейвлета на этапе синтеза.

Условие (14) упрощается при подстановке в качестве анализирующего вейвлета 5-функции, сосредоточенной на прямой

7 dy - _

]-^Ь(уе) = 2я,                                         (15)

-ооМ где e = (cos(p,sin

h(y)= Jh(p)e”1YPdp                                            (16)

—00

В результате свертки с h(p) по анализирующему вейвлет-образу Wf находится синтезирующий вейвлет-образ W'f:

2л со г г              л | о - pcosm)

W'f(b,a) = j dcp J Wf(p,a,

1 hl----------jdp.                            (17)

0   -00

На рис.2 показаны синтезирующие вейвлет-образы W'f тестовой функции при различных уровнях шума w, полученные сверткой с мексиканской шляпой h(p)=(l-p²)exp(-p²/2). Сравнивая синтезирующие вейвлет-образы тестовой функции (см. рис.2) можно отметить, что влияние шума проявляется в виде характерной структуры в области малых масштабов. Соответственно, при интегрировании по а необходимо найти оптимальную границу а,™, отсекающую область высокочастотного шума

В общем случае алгоритм, обеспечивающий точное восстановление функции по проекционным данным без шума [2],

f(p)= lim fep                                               (18)

г.->0,р->со должен быть для обеспечения устойчивости в условиях шума преобразован к виду

f(p)=         lim fe р,                                    (19)

^а1п^п (ч ),р —*amax(w)

-1                    -оз                   о                    оз                    1

Рис.2. Вейвлет-образ W' тестовой функции при различных уровнях шума w

где w - уровень шума, а е и р - пределы интегрирования в формуле lr      da

В последней формуле С- нормирующий коэффициент. Для мексиканской шляпы, использовавшейся в качестве синтезирующего вейвлета, С=2л.

Устойчивость реконструкции. Вследствие ограничения пределов формула (19) обеспечивает только приближенное восстановление f(p) по тем масштабам, где влияние шума невелико. С уменьшением уровня шума пределы интегрирования расширяются и функция восстанавливается все более точно, стремясь к истинному значению. .

Наглядно a„in, а™ можно определить по графику энергетического вейвлет-спектра, показывающему распределение энергии сигнала по масштабам:

Ef(a)= JjW f(b,a)j2 db (21)

На рис. 3 представлен энергетический спектр тестовой функции (13), вычисленный по 256 пространственным отсчетам для 16 масштабных коэффициентов. Ef₀ - график спектра без шума, ЕЕо - спектр функции, искаженный влиянием шума w=20%, Ef^ -спектр с шумом 40% (по горизонтальной оси - а в логарифмическом масштабе).

$20

$40

Рис.З. Энергетический вейвлет-спектр тестовой функции

Видно, что шум искажает спектр сигнала преимущественно в области малых а, поэтому для устойчивого восстановления тестовой функции достаточно выбрать нижнюю границу а_пмп в соответствии с уровнем шума. Таким образом, обеспечение устойчивости в данном случае обеспечивается «фильтром больших масштабов»

Если требуется реконструкция ограниченного масштабного диапазона, то пределы интегрирования агаи> и апих выбираются таким образом, чтобы одновременно обеспечить устойчивость и подавить ненужные масштабы. В этом случае работает «полосовой фильтр». В общем случае для восстановления отдельных масштабов соответствующим образом выбираются пределы в формуле (20).

На рис.4 показаны результаты восстановления тестовой функции по 256 дискретным отсчетам Радон-образа с шумом при 64 масштабных коэффициентах, изменяющихся в пределах от атах=10 до атт=0,02 (верхний график) и атт=0,01 (нижний график) (f - исходная функция, f0 - результат восстановления без шума, f20 - w=20%, f40 -w=40%).

Особенности поведения восстановленной функции вблизи точки р=0 объясняются когерентным усилением влияния шума (только для центральносимметричных функций) на масштабах порядка ширины минимального вейвлета. Дальнейшее уменьшении а„ип приводит к потере устойчивости даже при восстановлении тестовой функции по Ра-дон-образу без шума (w-0%) вследствие дискретности проекционных данных

Рис 4. Реконструкция тестовой функции при различных значениях а,,,,,,.

Выводы. Рассмотренный алгоритм вейвлет-томографии позволяет провести одновременно с восстановлением по проекционным данным фильтрацию шума и тем самым обеспечить устойчивость решения. При этом реконструкция только интересующих масштабов осуществляется выбором границ a_miri и a^x, причем в наиболее распространенном на практике случае высокочастотного шума достаточно определить только amm. Общий алгоритм предполагает поиск весовой функции для каждого масштаба а в зависимости от распределения шума.