В этом примере показано, как анализировать 3D данные с помощью инструмента трехмерного вейвлет-анализа и как отображать низкочастотные и высокочастотные компоненты вдоль данного среза. Пример посвящен магнитно-резонансным изображениям.
Ключевой особенностью этого анализа является отслеживание оптимальной или, по меньшей мере, хорошей, основанной на вейвлет разреженности изображения, которая является самым низким процентом коэффициентов преобразования, достаточным для восстановления диагностического качества.
Чтобы проиллюстрировать это, мы сохраняем приближение 3D МРТ, чтобы показать снижение сложности. Результат может быть улучшен, если изображения были преобразованы и реконструированы из самых больших коэффициентов преобразования, где определение качества оценивается медицинскими специалистами.
Мы увидим, что вейвлет-преобразование для изображений мозга позволяет эффективно и точно реконструировать только 5-10% коэффициентов.
Сначала загрузите wmri.mat файл, созданный из набора данных МРТ, поставляемого с 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.3043
Реконструированные приближения и детали вдоль 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 визуализацию с помощью анимации для различных фрагментов.