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