Документация

Около инвариантности сдвига двойного древовидного DWT

DWT страдает от отклонения сдвига, означая, что маленькие сдвиги во входном сигнале или изображении могут вызвать существенные изменения в распределении энергии сигнала/изображения через шкалы в коэффициентах DWT. DT-CWT является приблизительно инвариантом сдвига.

Чтобы продемонстрировать это на тестовом сигнале, создайте два переключенных импульса дискретного времени 128 выборок в длине. Один сигнал имеет модульный импульс в демонстрационных 60, в то время как другой сигнал имеет модульный импульс в демонстрационных 64. Оба сигнала ясно имеют модульную энергию (норма).

kronDelta1 = zeros(128,1);
kronDelta1(60) = 1;
kronDelta2 = zeros(128,1);
kronDelta2(64) = 1;

Получите DWT и двойной древовидный DWT двух сигналов вниз к уровню 3 с вейвлетом и масштабирующимися фильтрами длины 14. Извлеките коэффициенты детали уровня 3 для сравнения.

J = 3;
dwt1 = dddtree('dwt',kronDelta1,J,'sym7');
dwt2 = dddtree('dwt',kronDelta2,J,'sym7');
dwt1Cfs = dwt1.cfs{J};
dwt2Cfs = dwt2.cfs{J};

dt1 = dddtree('cplxdt',kronDelta1,J,'dtf3');
dt2 = dddtree('cplxdt',kronDelta2,J,'dtf3');
dt1Cfs = dt1.cfs{J}(:,:,1)+1i*dt1.cfs{J}(:,:,2);
dt2Cfs = dt2.cfs{J}(:,:,1)+1i*dt2.cfs{J}(:,:,2);

Постройте абсолютные значения DWT и коэффициентов DT-CWT для двух сигналов на уровне 3 и вычислите энергию (нормы в квадрате) коэффициентов.

figure;
subplot(1,2,1)
stem(abs(dwt1Cfs),'markerfacecolor',[0 0 1]);
title(['DWT Squared 2-norm = ' num2str(norm(dwt1Cfs,2)^2,3)],...
    'fontsize',10);
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title(['DWT Squared 2-norm = ' num2str(norm(dwt2Cfs,2)^2,3)],...
    'fontsize',10);

figure;
subplot(1,2,1)
stem(abs(dt1Cfs),'markerfacecolor',[0 0 1]);
title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt1Cfs,2)^2,3)],...
    'fontsize',10);
subplot(1,2,2)
stem(abs(dwt2Cfs),'markerfacecolor',[0 0 1])
title(['Dual-tree DWT Squared 2-norm = ' num2str(norm(dt2Cfs,2)^2,3)],...
    'fontsize',10);

Обратите внимание, что четыре выборки переключаются на нижний регистр, сигнал вызвал изменение на почти 6,5% в энергии уровня 3 коэффициенты вейвлета DWT. Однако двойные древовидные коэффициенты вейвлета показывают только изменение на 0,3% в энергии.

Чтобы продемонстрировать утилиту аппроксимированной инвариантности сдвига в действительных данных, мы анализируем электрокардиограмму (ECG) сигнал. Интервал выборки для сигнала ECG является 1/180 секундами. Данные взяты из Percival & Walden (2000), p.125 (данные, первоначально обеспеченные Уильямом Константином и На Reinhall, Вашингтонский университет). Для удобства мы берем данные, чтобы запуститься в t=0.

load wecg
dt = 1/180;
t = 0:dt:(length(wecg)*dt)-dt;
figure;
plot(t,wecg)
xlabel('Seconds'); ylabel('Millivolts');

Большой положительный peaks на расстоянии приблизительно в 0,7 секунды является волнами R сердечного ритма. Во-первых, анализируйте сигнал с помощью критически выбранного DWT и постройте исходный сигнал наряду с коэффициентами вейвлета уровня 2 и уровня 3. Коэффициенты уровня 2 и уровня 3 были выбраны, потому что волны R изолируются наиболее заметно в тех шкалах для данной выборки уровня.

J = 6;
dtDWT1 = dddtree('dwt',wecg,J,'farras');
details = zeros(2048,3);
details(2:4:end,2) = dtDWT1.cfs{2};
details(4:8:end,3) = dtDWT1.cfs{3};
subplot(311)
stem(t,details(:,2),'Marker','none','ShowBaseline','off')
title('Level 2'); ylabel('mV');
subplot(312)
stem(t,details(:,3),'Marker','none','ShowBaseline','off')
title('Level 3'); ylabel('mV');
subplot(313)
plot(t,wecg); title('Original Signal');
xlabel('Seconds'); ylabel('mV');

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

dtcplx1 = dddtree('cplxdt',wecg,J,'dtf3');
details = zeros(2048,3);
details(2:4:end,2) = dtcplx1.cfs{2}(:,1,1)+1i*dtcplx1.cfs{2}(:,1,2);
details(4:8:end,3) = dtcplx1.cfs{3}(:,1,1)+1i*dtcplx1.cfs{3}(:,1,2);
subplot(311)
stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off')
title('Level 2'); ylabel('mV');
subplot(312)
stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off')
title('Level 3'); ylabel('mV');
subplot(313)
plot(t,wecg); title('Original Signal');
xlabel('Seconds'); ylabel('mV');

И критически выбранный и двойной древовидный вейвлет преобразовывает, локализуют важную функцию формы волны ECG к подобным шкалам.

Важное приложение вейвлетов в 1D сигналах должно получить дисперсионный анализ шкалой. Это выдерживает обосновать, что этот дисперсионный анализ не должен быть чувствителен к циклическим сдвигам во входном сигнале. К сожалению, дело обстоит не так с критически выбранным DWT. Чтобы продемонстрировать это, мы циркулярный сдвиг сигнал ECG 4 выборками, анализирует непереключенные и переключенные сигналы с критически выбранным DWT и вычисляет распределение энергии через шкалы.

wecgShift = circshift(wecg,4);
dtDWT2 = dddtree('dwt',wecgShift,J,'farras');
sigenrgy = norm(wecg,2)^2;
enr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT1.cfs,'uni',0));
enr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,dtDWT2.cfs,'uni',0));
levels = {'D1';'D2';'D3';'D4';'D5';'D6';'A6'};
enr1 = enr1(:);
enr2 = enr2(:);
table(levels,enr1,enr2,'VariableNames',{'Level','enr1','enr2'})

ans =

  7x3 table

    Level     enr1      enr2 
    _____    ______    ______

    'D1'     4.1994    4.1994
    'D2'      8.425     8.425
    'D3'     13.381    10.077
    'D4'     7.0612    10.031
    'D5'     5.4606    5.4436
    'D6'     3.1273    3.4584
    'A6'     58.345    58.366

Обратите внимание на то, что коэффициенты вейвлета на уровнях 3 и 4 показывают приблизительно 3%-е изменения в энергии между исходным и переключенным сигналом. Затем, мы повторяем, что этот анализ с помощью комплексного двойного древовидного вейвлета преобразовывает.

dtcplx2 = dddtree('cplxdt',wecgShift,J,'dtf3');
cfs1 = cellfun(@squeeze,dtcplx1.cfs,'uni',0);
cfs2 = cellfun(@squeeze,dtcplx2.cfs,'uni',0);
cfs1 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs1,'uni',0);
cfs2 = cellfun(@(x) complex(x(:,1),x(:,2)),cfs2,'uni',0);
dtenr1 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs1,'uni',0));
dtenr2 = cell2mat(cellfun(@(x)(norm(x,2)^2/sigenrgy)*100,cfs2,'uni',0));
dtenr1 = dtenr1(:);
dtenr2 = dtenr2(:);
table(levels,dtenr1,dtenr2, 'VariableNames',{'Level','dtenr1','dtenr2'})

ans =

  7x3 table

    Level    dtenr1    dtenr2
    _____    ______    ______

    'D1'      4.936     4.936
    'D2'     6.6691    6.6691
    'D3'     12.682    12.611
    'D4'     8.3891    8.4808
    'D5'     5.8868    5.8791
    'D6'      3.053    3.0415
    'A6'     58.384    58.382

Двойное древовидное преобразование производит сопоставимый дисперсионный анализ шкалой для исходного сигнала и его циркулярной переключенной версии.

Направленная селективность в обработке изображений

Стандартная реализация DWT в 2D использовании отделимая фильтрация столбцов и строки изображения. LH, HL и вейвлеты HH для Добечиса меньше всего - асимметричный вейвлет фазы с 4 исчезающими моментами (sym4) показывают в следующей фигуре.

figure;
J = 5;
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(L,3*L);
dt = dddtree2('dwt',Y,J,'sym4');
dt.cfs{J}(N/3,N/2,1) = 1;
dt.cfs{J}(N/2,N/2+N,2) = 1;
dt.cfs{J}(N/2,N/2+2*N,3) = 1;
dwtImage = idddtree2(dt);
imagesc(dwtImage); axis xy; axis off;
title({'Critically Sampled DWT';'2D separable wavelets (sym4) -- LH, HL, HH'});

Обратите внимание на то, что LH и вейвлеты HL имеют ясные горизонтальные и вертикальные ориентации соответственно. Однако вейвлет HH на ультраправых смесях и +45 и-45 направлений степени, производя артефакт шахматной доски. Это смешивание ориентаций происходит из-за использования отделимых фильтров с действительным знаком. HH отделимый фильтр с действительным знаком имеет полосы пропускания во всех четырех высокочастотных углах 2D плоскости частоты.

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

Фильтры применились к столбцам, и строки могут быть от той же пары фильтра, или, или от различных пар фильтра. Следующий код показывает, что ориентация этих 12 вейвлетов, соответствующих действительным и мнимым частям комплекса, ориентировала двойной древовидный DWT.

J = 4;
L = 3*2^(J+1);
N = L/2^J;
Y = zeros(2*L,6*L);
wt = dddtree2('cplxdt',Y,J,'dtf3');
wt.cfs{J}(N/2,N/2+0*N,2,2,1) = 1;
wt.cfs{J}(N/2,N/2+1*N,3,1,1) = 1;
wt.cfs{J}(N/2,N/2+2*N,1,2,1) = 1;
wt.cfs{J}(N/2,N/2+3*N,1,1,1) = 1;
wt.cfs{J}(N/2,N/2+4*N,3,2,1) = 1;
wt.cfs{J}(N/2,N/2+5*N,2,1,1) = 1;
wt.cfs{J}(N/2+N,N/2+0*N,2,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+1*N,3,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+2*N,1,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+3*N,1,1,2) = 1;
wt.cfs{J}(N/2+N,N/2+4*N,3,2,2) = 1;
wt.cfs{J}(N/2+N,N/2+5*N,2,1,2) = 1;
waveIm = idddtree2(wt);
imagesc(waveIm); axis off;
title('Complex Dual-Tree 2D Wavelets');

Верхняя строка предыдущей фигуры показывает, что шесть направленных вейвлетов действительного ориентированного двойного древовидного вейвлета преобразовывают. Вторая строка показывает мнимые части. Вместе действительные и мнимые части формируются, комплекс ориентировался, двойной древовидный вейвлет преобразовывают. Действительные и мнимые части ориентированы в том же направлении. Можно использовать dddtree2 с опцией 'realdt', чтобы получить действительный ориентированный двойной древовидный DWT, который использует только действительные части. Используя действительный ориентированный двойной древовидный вейвлет преобразовывают, можно достигнуть направленной селективности, но вы не получаете полную выгоду от использования аналитических вейвлетов, таких как аппроксимированная инвариантность сдвига.

Представление ребра в двух измерениях

Аппроксимированная аналитичность и выборочная ориентация комплексных двойных древовидных вейвлетов обеспечивают наилучшее решение по стандартному 2D DWT в представлении ребер в изображениях. Чтобы проиллюстрировать это, мы анализируем тестовые изображения с ребрами, состоящими из строки и особенностей кривой в нескольких направлениях с помощью критически выбранного 2D DWT, и 2D комплекс ориентировался, двойное дерево преобразовывают. Во-первых, анализируйте изображение восьмиугольника, который состоит из особенностей строки.

load woctagon;
figure;
imagesc(woctagon); colormap gray;
title('Original Image'); axis equal; axis off;

Анализируйте изображение вниз к уровню 4 и восстановите приближение изображений на основе коэффициентов детали уровня 4.

dtcplx = dddtree2('cplxdt',woctagon,4,'dtf3');
dtDWT = dddtree2('dwt',woctagon,4,'farras');

dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));

dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));

dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(121)
imagesc(dtImage); axis equal; axis off; colormap gray;
title('Complex Oriented Dual-Tree');
subplot(122)
imagesc(dwtImage); axis equal; axis off; colormap gray;
title('DWT')

Затем, анализируйте восьмиугольник с гиперболическими ребрами. Ребра в гиперболическом восьмиугольнике являются особенностями кривой.

load woctagonHyperbolic;
figure;
imagesc(woctagonHyperbolic); colormap gray;
title('Octagon with Hyperbolic Edges'); axis equal; axis off;

Снова, анализируйте изображение вниз к уровню 4 и восстановите приближение изображений на основе коэффициентов детали уровня 4 и для стандартного 2D DWT и для комплекса, ориентированного на двойной древовидный DWT.

dtcplx = dddtree2('cplxdt',woctagonHyperbolic,4,'dtf3');
dtDWT = dddtree2('dwt',woctagonHyperbolic,4,'farras');

dtcplx.cfs{1} = zeros(size(dtcplx.cfs{1}));
dtcplx.cfs{2} = zeros(size(dtcplx.cfs{2}));
dtcplx.cfs{3} = zeros(size(dtcplx.cfs{3}));
dtcplx.cfs{5} = zeros(size(dtcplx.cfs{5}));

dtDWT.cfs{1} = zeros(size(dtDWT.cfs{1}));
dtDWT.cfs{2} = zeros(size(dtDWT.cfs{2}));
dtDWT.cfs{3} = zeros(size(dtDWT.cfs{3}));
dtDWT.cfs{5} = zeros(size(dtDWT.cfs{5}));

dtImage = idddtree2(dtcplx);
dwtImage = idddtree2(dtDWT);
subplot(121)
imagesc(dtImage); axis equal; axis off; colormap gray;
title('DT-CWT');
subplot(122)
imagesc(dwtImage); axis equal; axis off; colormap gray;
title('DWT')

Обратите внимание на то, что звонящие артефакты, очевидные в 2D критически, выбрали DWT, отсутствуют в 2D DT-CWT обоих изображений. DT-CWT более искренне воспроизводит особенности кривой и строка.

Отобразите шумоподавление

Из-за способности изолировать отличные ориентации в отдельных поддиапазонах, двойной древовидный DWT часто может превзойти стандартный отделимый DWT по характеристикам в приложениях как шумоподавление изображений. Чтобы продемонстрировать это, используйте функцию помощника helperCompare2DDenoising. Функция помощника загружает изображение и добавляет нулевой средний белый Гауссов шум с. Для предоставленной пользователями области значений порогов функция сравнивает шумоподавление с помощью мягкой пороговой обработки для критически выбранного DWT, действительного ориентированного двойного древовидного DWT, и комплекс ориентировал двойной древовидный DWT. Для каждого порогового значения отображены среднеквадратичная (RMS) ошибка и пиковое отношение сигнал-шум (PSNR).

numex = 3;
helperCompare2DDenoising(numex,0:2:100,'PlotMetrics');

И ориентированное действительное и комплекс ориентировались, двойные древовидные DWTs превосходят стандартный DWT по характеристикам по ошибке RMS и PSNR.

Затем, получите изображения denoised для порогового значения 25, который равен стандартному отклонению аддитивного шума.

numex = 3;
helperCompare2DDenoising(numex,25,'PlotImage');

С пороговым значением, равным стандартному отклонению аддитивного шума, ориентировался комплекс, двойное дерево преобразовывают, обеспечивает PSNR на почти 4 дБ выше, чем стандартный 2D DWT.

Направленная селективность в 3D

Звонящие артефакты, наблюдаемые с отделимым DWT в двух измерениях, усилены при расширении анализа вейвлета к более высоким размерностям. DT-CWT позволяет вам поддержать направленную селективность в 3D с минимальным сокращением. В 3D в двойном древовидном преобразовании существует 28 поддиапазонов вейвлета.

Чтобы продемонстрировать направленную селективность 3D двойного древовидного вейвлета преобразовывают, визуализируют пример 3D изоповерхности и 3D двойного дерева и отделимых вейвлетов DWT. Во-первых, визуализируйте действительные и мнимые части отдельно двух двойных древовидных поддиапазонов.

helperVisualize3D('Dual-Tree',28,'separate');
helperVisualize3D('Dual-Tree',25,'separate');

Красный фрагмент графика изоповерхности указывает на положительное отклонение вейвлета от нуля, в то время как синий обозначает отрицательное отклонение. Можно ясно видеть направленную селективность на пробеле действительных и мнимых частей двойных древовидных вейвлетов. Теперь визуализируйте один из двойных древовидных поддиапазонов с действительными и мнимыми графиками, построенными вместе как одна изоповерхность.

helperVisualize3D('Dual-Tree',25,'real-imaginary');

Предыдущий график демонстрирует, что действительные и мнимые части являются переключенными версиями друг друга на пробеле. Это отражает то, что мнимая часть комплексного вейвлета является аппроксимированным Гильбертовым преобразованием действительной части. Затем, визуализируйте изоповерхность действительного ортогонального вейвлета в 3D для сравнения.

helperVisualize3D('DWT',7);

Смешивание ориентаций, наблюдаемых в 2D DWT, является еще более явным в 3D. Так же, как в 2D случае, смешивание ориентаций в 3D приводит к значительному вызову или блокированию артефактов. Чтобы продемонстрировать это, исследуйте 3D DWT и детали вейвлета DT-CWT сферического объема. Сфера 64 64 64.

load sphr;
[A,D] = dualtree3(sphr,2,'excludeL1');
A = zeros(size(A));
sphrDTCWT = idualtree3(A,D);
montage(reshape(sphrDTCWT,[64 64 1 64]),'DisplayRange',[])
title('DT-CWT Level 2 Details');

Сравните предыдущий график против деталей второго уровня на основе отделимого DWT.

sphrDEC = wavedec3(sphr,2,'sym4','mode','per');
sphrDEC.dec{1} = zeros(size(sphrDEC.dec{1}));
for kk = 2:8
    sphrDEC.dec{kk} = zeros(size(sphrDEC.dec{kk}));
end
sphrrecDWT = waverec3(sphrDEC);
figure;
montage(reshape(sphrrecDWT,[64 64 1 64]),'DisplayRange',[])
title('DWT Level 2 Details');

Увеличьте масштаб изображений и в DT-CWT и в монтаже DWT, и вы будете видеть, как видный блокирующиеся артефакты в деталях DWT сравниваются с DT-CWT.

Шумоподавление объема

Подобно 2D случаю направленная селективность 3D DT-CWT часто приводит к улучшениям шумоподавления объема.

Чтобы продемонстрировать это, рассмотрите набор данных MRI, состоящий из 16 срезов. Гауссов шум со стандартным отклонением 10 был добавлен к исходному набору данных. Отобразите шумный набор данных.

load MRI3D;
montage(reshape(noisyMRI,[128 128 1 16]),'DisplayRange',[]);

Обратите внимание на то, что исходный ОСШ до шумоподавления составляет приблизительно 11 дБ.

20*log10(norm(origMRI(:),2)/norm(origMRI(:)-noisyMRI(:),2))

ans =

   11.2997

Denoise набор данных MRI вниз к уровню 4 с помощью и DT-CWT и DWT. Подобные длины фильтра вейвлета используются в обоих случаях. Постройте получившийся ОСШ как функцию порога. Отобразите результаты denoised и для DT-CWT и для DWT, полученного в лучшем ОСШ.

[imrecDTCWT,imrecDWT] = helperCompare3DDenoising(origMRI,noisyMRI);
figure;
montage(reshape(imrecDTCWT,[128 128 1 16]),'DisplayRange',[]);
title('DT-CWT Denoised Volume');
figure;
montage(reshape(imrecDWT,[128 128 1 16]),'DisplayRange',[]);
title('DWT Denoised Volume');

Сводные данные

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

Дополнительные материалы для чтения

Chen, H. & Kingsbury, N.G. "Efficient registration of nonrigid 3-D
bodies", IEEE Transactions on Image Processing, Vol. 21, No. 1, Jan
2012, pp. 262-272.

Kingsbury, N.G. "Complex Wavelets for Shift Invariant Analysis and
Filtering of Signals". Journal of Applied and Computational Harmonic
Analysis. Vol 10, Number 3, May 2001, pp. 234-253.

Percival, D.B. and A.T. Walden. "Wavelet Methods for Time Series
Analysis", Cambridge University Press, 2000.

Selesnick, I., Baraniuk, R.G., and N.G. Kingsbury. "The Dual-Tree
Complex Wavelet Transform." IEEE Signal Processing Magazine. Vol. 22,
Number 6, November, 2005, pp. 123-151.

Selesnick, I. "The Double Density DWT". Wavelets in Signal and Image
Analysis: From Theory to Practice (A.A Petrosian, F.G. Meyer, eds.),
Norwell, MA: Kluwer Academic Publishers, 2001, pp. 39-66.

Selesnick, I. "The Double-Density Dual-Tree Wavelet Transform". IEEE
Transactions on Signal Processing. Vol. 52, Number 5, May 2004, pp.
1304-1314.