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