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

Скрипт Open Live Script

В этом примере показано, как анализировать 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

Мы теперь переключаемся на срез 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

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.3044

Отобразите компоненты 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

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');

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');

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

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

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