www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/wavelet3ddemo.m
%% Wavelet Analysis for 3D Data % This example shows how to analyze 3D data using the three-dimensional % wavelet analysis tool, and how to display low-pass and high-pass % components along a given slice. The example focuses on magnetic resonance % images. % % A key feature of this analysis is to track the optimal, or at least a % good, wavelet-based sparsity of the image which is the lowest percentage % of transform coefficients sufficient for diagnostic-quality % reconstruction. % % To illustrate this, we keep the approximation of a 3D MRI to show the % complexity reduction. The result can be improved if the images were % transformed and reconstructed from the largest transform coefficients % where the definition of the quality is assessed by medical specialists. % % We will see that Wavelet transform for brain images allows efficient and % accurate reconstructions involving only 5-10% of the coefficients. % Copyright 2009-2014 The MathWorks, Inc. %% Load and Display 3D MRI Data % First, load the |wmri.mat| file which is built from the MRI data set that % comes with MATLAB(R). load wmri %% % We now display some slices along the Z-orientation of the original brain % data. 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 %% % We now switch to the Y-orientation slice. 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 %% Multilevel 3D Wavelet Decomposition % Compute the wavelet decomposition of the 3D data at level 3. n = 3; % Decomposition Level w = 'sym4'; % Near symmetric wavelet WT = wavedec3(X,n,w); % Multilevel 3D wavelet decomposition. %% Multilevel 3D Wavelet Reconstruction % Reconstruct from coefficients the approximations and details for levels 1 to 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 %% % Check for perfect reconstruction. err = zeros(1,n); for k = 1:n E = double(X)-A{k}-D{k}; err(k) = max(abs(E(:))); end disp(err) %% Display Low-Pass and High-Pass Components % The reconstructed approximations and details along the Z-orientation are % displayed below. 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 Display of Original Data and Approximation at Level 2 % The size of the 3D original array X is *(128 x 128 x 27) = 442368*. We % can use a 3D display to show it. 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'); %% % The 3D array of the coefficients of approximation at level 2, whose size % is *(22 x 22 x 9) = 4356*, is less than 1% the size of the original data. % With these coefficients, we can reconstruct A2, the approximation at % level 2, which is a kind of compression of the original 3D array. % A2 can also be shown using a 3D display. 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'); %% Summary % This example shows how to use 3D wavelet functions to analyze 3D data. The % Wavelet Analyzer app, |waveletAnalyzer|, lets you perform the same steps % more easily, such as simulating 3D visualization using an animation % across different slices.