Анализ вейвлета для 3D данных

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

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

Чтобы проиллюстрировать это, мы сохраняем приближение 3D MRI, чтобы показать сокращение сложности. Результат может быть улучшен, если изображения были преобразованы и восстановлены от самого большого, преобразовывают коэффициенты, где определение качества оценено медицинскими специалистами.

Мы будем видеть, что Вейвлет преобразовывает для изображений мозга, позволяет эффективным и точным реконструкциям включающие только 5-10% коэффициентов.

Загрузите и отобразите 3D данные MRI

Во-первых, загрузите wmri.mat файл, который создается из набора данных MRI, который идет с MATLAB®.

load wmri

Мы теперь отображаем некоторые срезы вдоль Z-ориентации исходных мозговых данных.

map = pink(90);
idxImages = 1:3:size(X,3);
figure('DefaultAxesXTick',[],'DefaultAxesYTick',[],...
    'DefaultAxesFontSize',8,'Color','w')
colormap(map)
for k = 1:9
    j = idxImages(k);
    subplot(3,3,k)
    image(X(:,:,j))
    xlabel(['Z = ' int2str(j)])
    if k==2
        title('Some slices along the Z-orientation of the original brain data')
    end
end

Figure contains 9 axes. Axes 1 contains an object of type image. Axes 2 with title Some slices along the Z-orientation of the original brain data contains an object of type image. Axes 3 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image. Axes 7 contains an object of type image. Axes 8 contains an object of type image. Axes 9 contains an object of type image.

Мы теперь переключаемся на срез Y-ориентации.

perm = [1 3 2];
XP = permute(X,perm);
figure('DefaultAxesXTick',[],'DefaultAxesYTick',[],...
    'DefaultAxesFontSize',8,'Color','w')
colormap(map)
for k = 1:9
    j = idxImages(k);
    subplot(3,3,k)
    image(XP(:,:,j))
    xlabel(['Y = ' int2str(j)])
    if k==2
    	title('Some slices along the Y-orientation')
    end
end

Figure contains 9 axes. Axes 1 contains an object of type image. Axes 2 with title Some slices along the Y-orientation contains an object of type image. Axes 3 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image. Axes 7 contains an object of type image. Axes 8 contains an object of type image. Axes 9 contains an object of type image.

clear XP

Многоуровневое 3D разложение вейвлета

Вычислите разложение вейвлета 3D данных на уровне 3.

n = 3;                   % Decomposition Level 
w = 'sym4';              % Near symmetric wavelet
WT = wavedec3(X,n,w);    % Multilevel 3D wavelet decomposition.

Многоуровневая 3D реконструкция вейвлета

Восстановите от коэффициентов приближения и детали для уровней 1 - 3.

A = cell(1,n);
D = cell(1,n);
for k = 1:n
    A{k} = waverec3(WT,'a',k);   % Approximations (low-pass components)
    D{k} = waverec3(WT,'d',k);   % Details (high-pass components)
end

Проверяйте на совершенную реконструкцию.

err = zeros(1,n);
for k = 1:n
    E = double(X)-A{k}-D{k};
    err(k) = max(abs(E(:)));
end
disp(err)
   1.0e-09 *

    0.3043    0.3043    0.3043

Отобразите компоненты Lowpass и Высокой Передачи

Восстановленные приближения и детали вдоль Z-ориентации отображены ниже.

nbIMG = 6;
idxImages_New = [1 7 10 16 19 25];
for ik = 1:nbIMG
    j = idxImages_New(ik);
    figure('DefaultAxesXTick',[],'DefaultAxesYTick',[],...
        'DefaultAxesFontSize',8,'Color','w')
    colormap(map)
    for k = 1:n
        labstr = [int2str(k) ' - Z = ' int2str(j)];
        subplot(2,n,k)
        image(A{k}(:,:,j))
        xlabel(['A' labstr])
        if k==2
        	title(['Approximations and details at level 1 to 3 - Slice = ' num2str(j)])
        end
        subplot(2,n,k+n)
        imagesc(abs(D{k}(:,:,j)))
        xlabel(['D' labstr])
    end
end

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 1 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 7 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 10 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 16 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 19 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

Figure contains 6 axes. Axes 1 contains an object of type image. Axes 2 contains an object of type image. Axes 3 with title Approximations and details at level 1 to 3 - Slice = 25 contains an object of type image. Axes 4 contains an object of type image. Axes 5 contains an object of type image. Axes 6 contains an object of type image.

3D отображение исходных данных и приближения на уровне 2

Размер 3D исходного массива X (128 x 128 x 27) = 442368. Мы можем использовать 3D отображение, чтобы показать его.

figure('DefaultAxesXTick',[],'DefaultAxesYTick',[],...
        'DefaultAxesFontSize',8,'Color','w')
XR = X;
Ds = smooth3(XR);
hiso = patch(isosurface(Ds,5),'FaceColor',[1,.75,.65],'EdgeColor','none');
hcap = patch(isocaps(XR,5),'FaceColor','interp','EdgeColor','none');
colormap(map)
daspect(gca,[1,1,.4])
lightangle(305,30); 
fig = gcf;
fig.Renderer = 'zbuffer';
lighting phong
isonormals(Ds,hiso)
hcap.AmbientStrength = .6;
hiso.SpecularColorReflectance = 0;
hiso.SpecularExponent = 50;
ax = gca;
ax.View = [215,30];
ax.Box = 'On';
axis tight
title('Original Data')

Figure contains an axes. The axes with title Original Data contains 2 objects of type patch.

3D массив коэффициентов приближения на уровне 2, размер которого (22 x 22 x 9) = 4356, меньше 1% размер исходных данных. С этими коэффициентами мы можем восстановить A2, приближение на уровне 2, который является своего рода сжатием исходного 3D массива. A2 можно также показать с помощью 3D отображения.

figure('DefaultAxesXTick',[],'DefaultAxesYTick',[],...
        'DefaultAxesFontSize',8,'Color','w')
XR = A{2};
Ds = smooth3(XR);
hiso = patch(isosurface(Ds,5),'FaceColor',[1,.75,.65],'EdgeColor','none');
hcap = patch(isocaps(XR,5),'FaceColor','interp','EdgeColor','none');
colormap(map)
daspect(gca,[1,1,.4])
lightangle(305,30); 
fig = gcf;
fig.Renderer = 'zbuffer';
lighting phong
isonormals(Ds,hiso)
hcap.AmbientStrength = .6;
hiso.SpecularColorReflectance = 0;
hiso.SpecularExponent = 50;
ax = gca;
ax.View = [215,30];
ax.Box = 'On';
axis tight
title('Approximation at level 2')

Figure contains an axes. The axes with title Approximation at level 2 contains 2 objects of type patch.

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

В этом примере показано, как использовать 3D функции вейвлета, чтобы анализировать 3D данные. Приложение Wavelet Analyzer, waveletAnalyzer, позволяет вам выполнить те же шаги более легко, такие как симуляция 3D визуализации с помощью анимации через различные срезы.

Для просмотра документации необходимо авторизоваться на сайте