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

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

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

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

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 и вычислите энергию (в квадрате $${\ell }^2$$ нормы) коэффициентов.

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 [3], 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(3,1,1)
stem(t,details(:,2),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
subplot(3,1,2)
stem(t,details(:,3),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
subplot(3,1,3)
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(3,1,1)
stem(t,real(details(:,2)),'Marker','none','ShowBaseline','off')
title('Level 2')
ylabel('mV')
subplot(3,1,2)
stem(t,real(details(:,3)),'Marker','none','ShowBaseline','off')
title('Level 3')
ylabel('mV')
subplot(3,1,3)
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=7×3 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=7×3 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';'2-D separable wavelets (sym4) -- LH, HL, HH'})

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

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

Фильтры применились к столбцам, и строки могут быть от той же пары фильтра, $$\{h_0,h_1\}$$ или $$\{g_0,g_1\}$$, или от различных пар фильтра, $$\{h_0,g_1\},\{g_0,h_1\}$$. Следующий код показывает, что ориентация этих 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 2-D 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(1,2,1)
imagesc(dtImage)
axis equal
axis off
colormap gray
title('Complex Oriented Dual-Tree')
subplot(1,2,2)
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(1,2,1)
imagesc(dtImage)
axis equal
axis off
colormap gray
title('DT-CWT')
subplot(1,2,2)
imagesc(dwtImage)
axis equal
axis off
colormap gray
title('DWT')

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

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

Из-за способности изолировать отличные ориентации в отдельных поддиапазонах, двойной древовидный DWT часто может превзойти стандартный отделимый DWT по характеристикам в приложениях как шумоподавление изображений. Чтобы продемонстрировать это, используйте функцию помощника helperCompare2DDenoising. Функция помощника загружает изображение и добавляет нулевой средний белый Гауссов шум с $$\sigma =\mathrm{}25$$. Для предоставленной пользователями области значений порогов функция сравнивает шумоподавление с помощью мягкой пороговой обработки в критически произведенном 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.

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

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

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

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

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

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

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

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

helperVisualize3D('DWT',7);

Смешивание ориентаций, наблюдаемых в 2D DWT, является еще более явным в 3-D. Так же, как в 2D случае, смешивание ориентаций в 3-D приводит к значительному вызову или блокированию артефактов. Чтобы продемонстрировать это, исследуйте 3-D 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 случаю направленная селективность 3-D 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. Мы продемонстрировали, как эти свойства могут привести к улучшенной производительности в анализе сигнала, представлении ребер в изображениях и объемах, и шумоподавлении объема и изображении.