В этом примере показано, как анализировать 3D данные с помощью 3D аналитического инструмента вейвлета, и как отобразить lowpass и компоненты высокой передачи вдоль данного среза. Пример фокусируется на изображениях магнитного резонанса.
Ключевая возможность этого анализа должна отследить оптимальное, или по крайней мере хорошая, основанная на вейвлете разреженность изображения, которое является самым низким процентом коэффициентов преобразования, достаточных для реконструкции диагностического качества.
Чтобы проиллюстрировать это, мы сохраняем приближение 3D MRI, чтобы показать сокращение сложности. Результат может быть улучшен, если изображения были преобразованы и восстановлены от самого большого, преобразовывают коэффициенты, где определение качества оценено медицинскими специалистами.
Мы будем видеть, что Вейвлет преобразовывает для изображений мозга, позволяет эффективным и точным реконструкциям включающие только 5-10% коэффициентов.
Во-первых, загрузите 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 данных на уровне 3.
n = 3; % Decomposition Level w = 'sym4'; % Near symmetric wavelet WT = wavedec3(X,n,w); % Multilevel 3D wavelet decomposition.
Восстановите от коэффициентов приближения и детали для уровней 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
Восстановленные приближения и детали вдоль 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 исходного массива 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 визуализации с помощью анимации через различные срезы.