Denoise Simulated Couette flow
The following compares divergence-free wavelet denoising with existing methods on a simulated Couette flow field.
Contents
- Clear all and check path
- Generate flow data and set parameters
- Add Noise
- DivFree Wavelet with Manual Threshold
- DivFree Wavelet with manual threshold and partial cycle spinning
- DivFree Wavelet with SureShrink and MAD sigma estimation
- DivFree Wavelet with SureShrink, MAD and random cycle spinning
- Finite Difference Method
- DivFree Radial Basis Function
- Plot all
- Visualize DivFree Wavelet Coefficients
- Plot Velocity profiles
Clear all and check path
close all clc clear if ~exist('dfwavelet_thresh','file') error('Cannot find dfwavelet functions. run setPath!'); end
Generate flow data and set parameters
N = 64; % Image width Rout = N/2; % Radius of outer cylinder Rin = N/4; % Radius of inner cylinder Win = 10; % vMax = 1; % Maximum velocity % Generate Couette flow field [vx,vy,vz,imMag] = genCouette(vMax,N,Rout,Rin,Win); % Set parameters vMag = sqrt(vx.^2+vy.^2+vz.^2); % Velocity magnitude vMax = max(vMag(:)); % Maximum speed FOV = size(vMag); % Field of view N = FOV(3); % Number of slices (for plotting) res = [1,1,1]; % Relative resolution ph0 = zeros(FOV); % Reference phase in phase contrast imMask = imMag; % Segmentation mask % Plot figure,imshow_flow(imMag,vx,vy,vz,vMax,[1,2,3]),title('Original Flow Field','FontSize',14); figure,imshow_flowmag(imMag,vx,vy,vz,vMax,[1,2,3]),title('Original Velocity Magnitude','FontSize',14); pause(1);


Add Noise
sigma = 0.2*max(imMag(:)); % Noise standard deviation gamma = 26753; % Gyromagnetic ratio of proton * 2pi M1 = pi/(2*vMax*gamma); % First moment % Transform velocity to 5-point method phase encodes [s1,s2,s3,s4,s5] = fivePointPC(vx,vy,vz,imMag,ph0,M1); % Add noise s1 = s1+sigma/sqrt(2)*(randn(FOV)+1i*randn(FOV)); s2 = s2+sigma/sqrt(2)*(randn(FOV)+1i*randn(FOV)); s3 = s3+sigma/sqrt(2)*(randn(FOV)+1i*randn(FOV)); s4 = s4+sigma/sqrt(2)*(randn(FOV)+1i*randn(FOV)); s5 = s5+sigma/sqrt(2)*(randn(FOV)+1i*randn(FOV)); % Transform phase encodes back to velocity [vxN,vyN,vzN,imMagN] = invfivePointPC(s1,s2,s3,s4,s5,M1); % Crop velocities using image mask [vxN,vyN,vzN] = maskIM(imMask,vxN,vyN,vzN); % Plot figure,imshow_flow(imMag,vxN,vyN,vzN,vMax,[1,2,3]),title('Noisy Flow Field','FontSize',14); figure,imshow_flowmag(imMag,vxN,vyN,vzN,vMax,[1,2,3]),title('Noisy Velocity Magnitude','FontSize',14); % Calculate errors in region of interest disp('Noisy Flow Error') [vNRMSE_Noise,vMagErr_Noise,angErr_Noise] = calcVelError(imMask,vx,vy,vz,vxN,vyN,vzN); PVNR = 8*log10(1/vNRMSE_Noise); fprintf('PVNR: \t\t\t%.2fdB\nNRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',PVNR,vNRMSE_Noise,vMagErr_Noise,angErr_Noise); pause(1);
Noisy Flow Error PVNR: 8.83dB NRMSE: 0.078838 vMag Error: 0.047772 Absolute Angle Error: 0.094471


DivFree Wavelet with Manual Threshold
To denoise, we apply the soft-threshold operation on the divergence-free component and the non-divergence-free component separately
minSize = 8*ones(1,3); % Smallest wavelet level size dfThresh = 0.05; % Threshold for divergence-free components nThresh = 0.08; % Threshold for non-divergence-free components % Forward DivFree Wavelet Transform [wcdf1,wcdf2,wcn,numLevels,wcSizes] = dfwavelet_forward(vxN,vyN,vzN,minSize,res); % Soft Threshold [wcdf1,wcdf2,wcn] = dfsoft_thresh(wcdf1,wcdf2,wcn,minSize,res,FOV,dfThresh,nThresh); % Inverse DivFree Wavelet Transform [vxDFW,vyDFW,vzDFW] = dfwavelet_inverse(wcdf1,wcdf2,wcn,minSize,res,FOV); % Alternatively, you can use the following function to do the same things % as above: % [vxDFW,vyDFW,vzDFW] = dfwavelet_thresh(vxN,vyN,vzN,minSize,res,dfthresh,nthresh); % Plot figure,imshow_flow(imMag,vxDFW,vyDFW,vzDFW,vMax,[1,2,3]) title('Div Free Wavelet w/ Manual Thresholds (Flow Field)','FontSize',14) figure,imshow_flowmag(imMag,vxDFW,vyDFW,vzDFW,vMax,[1,2,3]) title('Div Free Wavelet w/ Manual Thresholds (Vel Mag)','FontSize',14) % Calculate errors disp('DivFree Wavelet w/ Manual Thresholds') [vNRMSE_DFW,vMagErr_DFW,angErr_DFW] = calcVelError(imMask,vx,vy,vz,vxDFW,vyDFW,vzDFW); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFW,vMagErr_DFW,angErr_DFW); pause(1);
DivFree Wavelet w/ Manual Thresholds NRMSE: 0.048441 vMag Error: 0.037032 Absolute Angle Error: 0.023616


DivFree Wavelet with manual threshold and partial cycle spinning
To remove the blocking artifacts, we do partial cycle spinning Here we do 2^3=8 random shifts
spins = 2; % Number of cycle spinning per dimension isRandShift = 1; % Use random shift minSize = 8*ones(1,3); % Smallest wavelet level size dfThresh = 0.05; % Threshold for divergence-free components nThresh = 0.08; % Threshold for non-divergence-free components % Initialize velocities vxDFW_s = zeros(FOV); vyDFW_s = zeros(FOV); vzDFW_s = zeros(FOV); % Cycle Spinning for i = 1:spins for j = 1:spins for k = 1:spins % Generate spins if isRandShift shifts = randi(8,1,3) - 1; else shifts = [i,j,k]; end % Circular shift vxN_spin = circshift(vxN, shifts ); vyN_spin = circshift(vyN, shifts ); vzN_spin = circshift(vzN, shifts ); % Forward DivFree Wavelet Transform [wcdf1_s,wcdf2_s,wcn_s,numLevels,wcSizes] = dfwavelet_forward(vxN_spin,vyN_spin,vzN_spin,minSize,res); % Soft Threshold [wcdf1_s,wcdf2_s,wcn_s] = dfsoft_thresh(wcdf1_s,wcdf2_s,wcn_s,minSize,res,FOV,dfThresh,nThresh); % Inverse DivFree Wavelet Transform [vxN_spin,vyN_spin,vzN_spin] = dfwavelet_inverse(wcdf1_s,wcdf2_s,wcn_s,minSize,res,FOV); % Circular unshift vxN_spin = circshift(vxN_spin, -shifts ); vyN_spin = circshift(vyN_spin, -shifts ); vzN_spin = circshift(vzN_spin, -shifts ); % Add back vxDFW_s = vxDFW_s + vxN_spin; vyDFW_s = vyDFW_s + vyN_spin; vzDFW_s = vzDFW_s + vzN_spin; end end end vxDFW_s = vxDFW_s/spins^3; vyDFW_s = vyDFW_s/spins^3; vzDFW_s = vzDFW_s/spins^3; % Alternatively, you can use the following function to do the same things % as above: % [vxDFW_s,vyDFW_s,vzDFW_s] = dfwavelet_thresh_spin(vxN,vyN,vzN,minSize,res,dfthresh,nthresh,spins,isRandShift); % Plot figure,imshow_flow(imMag,vxDFW_s,vyDFW_s,vzDFW_s,vMax,[1,2,3]) title('Div Free Wavelet w/ Manual Thresholds and Partial Cycle Spinning (Flow Field)','FontSize',14) figure,imshow_flowmag(imMag,vxDFW,vyDFW,vzDFW,vMax,[1,2,3]) title('Div Free Wavelet w/ Manual Thresholds and Partial Cycle Spinning(Vel Mag)','FontSize',14) % Calculate errors disp('DivFree Wavelet w/ Manual Thresholds and Partial Cycle Spinning') [vNRMSE_DFW_s,vMagErr_DFW_s,angErr_DFW_s] = calcVelError(imMask,vx,vy,vz,vxDFW_s,vyDFW_s,vzDFW_s); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFW_s,vMagErr_DFW_s,angErr_DFW_s); pause(1);
DivFree Wavelet w/ Manual Thresholds and Partial Cycle Spinning NRMSE: 0.032196 vMag Error: 0.024539 Absolute Angle Error: 0.014880


DivFree Wavelet with SureShrink and MAD sigma estimation
Here, we use Median Absolute Deviation to estimate noise std and then use SureShrink to find the optimal threshold that minimizes MSE
minSize = 8*ones(1,3); % Smallest wavelet level size percentZero = sum(imMask(:)==0)/numel(imMask); % Percentage of masked elements % Forward DivFree Wavelet Transform [wcdf1_sm,wcdf2_sm,wcn_sm,numLevels,wcSizes] = dfwavelet_forward(vxN,vyN,vzN,minSize,res); % Estimate noise standard deviation in velocity est_sigma = getMADsigma(wcn_sm,minSize,res,FOV,percentZero); % SureShrink [wcdf1_sm,wcdf2_sm,wcn_sm] = dfSUREshrink(wcdf1_sm,wcdf2_sm,wcn_sm,minSize,res,FOV,percentZero,est_sigma); % Inverse DivFree Wavelet Transform [vxDFW_sm,vyDFW_sm,vzDFW_sm] = dfwavelet_inverse(wcdf1_sm,wcdf2_sm,wcn_sm,minSize,res,FOV); % Alternatively, you can use the following function to do the same things % as above: % [vxDFW_sm,vyDFW_sm,vzDFW_sm] = dfwavelet_thresh_SURE_MAD(vxN,vyN,vzN,minSize,res); % Plot figure,imshow_flow(imMag,vxDFW_sm,vyDFW_sm,vzDFW_sm,vMax,[1,2,3]) title('Div Free Wavelet w/ SureShrink (Flow Field)','FontSize',14) figure,imshow_flowmag(imMag,vxDFW_sm,vyDFW_sm,vzDFW_sm,vMax,[1,2,3]) title('Div Free Wavelet w/ SureShrink (Vel Mag)','FontSize',14) % Calculate errors disp('DivFree Wavelet w/ SureShrink and MAD') [vNRMSE_DFW_sm,vMagErr_DFW_sm,angErr_DFW_sm] = calcVelError(imMask,vx,vy,vz,vxDFW_sm,vyDFW_sm,vzDFW_sm); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFW_sm,vMagErr_DFW_sm,angErr_DFW_sm); pause(1);
DivFree Wavelet w/ SureShrink and MAD NRMSE: 0.034796 vMag Error: 0.027149 Absolute Angle Error: 0.014671


DivFree Wavelet with SureShrink, MAD and random cycle spinning
We do partial cycle spinning with SureShrink to further improve the performance
spins = 2; % Number of cycle spinning per dimension isRandShift = 1; % Use random shift minSize = 8*ones(1,3); % Smallest wavelet level size % Initialize velocities vxDFW_sms = zeros(FOV); vyDFW_sms = zeros(FOV); vzDFW_sms = zeros(FOV); % Cycle Spinning for i = 1:spins for j = 1:spins for k = 1:spins % Generate spins if isRandShift shifts = randi(8,1,3) - 1; else shifts = [i,j,k]; end % Circular shift vxN_spin = circshift(vxN, shifts ); vyN_spin = circshift(vyN, shifts ); vzN_spin = circshift(vzN, shifts ); % Forward DivFree Wavelet Transform [wcdf1_sms,wcdf2_sms,wcn_sms,numLevels,wcSizes] = dfwavelet_forward(vxN_spin,vyN_spin,vzN_spin,minSize,res); % Estimate noise standard deviation in velocity est_sigma = getMADsigma(wcn_sms,minSize,res,FOV,percentZero); % SureShrink [wcdf1_sms,wcdf2_sms,wcn_sms] = dfSUREshrink(wcdf1_sms,wcdf2_sms,wcn_sms,minSize,res,FOV,percentZero,est_sigma); % Inverse DivFree Wavelet Transform [vxN_spin,vyN_spin,vzN_spin] = dfwavelet_inverse(wcdf1_sms,wcdf2_sms,wcn_sms,minSize,res,FOV); % Circular unshift vxN_spin = circshift(vxN_spin, -shifts ); vyN_spin = circshift(vyN_spin, -shifts ); vzN_spin = circshift(vzN_spin, -shifts ); % Add back vxDFW_sms = vxDFW_sms + vxN_spin; vyDFW_sms = vyDFW_sms + vyN_spin; vzDFW_sms = vzDFW_sms + vzN_spin; end end end vxDFW_sms = vxDFW_sms/spins^3; vyDFW_sms = vyDFW_sms/spins^3; vzDFW_sms = vzDFW_sms/spins^3; % Alternatively, you can use the following function to do the same things % as above: % [vxDFW_sms,vyDFW_sms,vzDFW_sms] = dfwavelet_thresh_SURE_MAD_spin(vxN,vyN,vzN,minSize,res,spins,isRandShift); % Plot figure,imshow_flow(imMag,vxDFW_sms,vyDFW_sms,vzDFW_sms,vMax,[1,2,3]) title('Div Free Wavelet w/ SureShrink and Partial Cycle Spinning (Flow Field)','FontSize',14) figure,imshow_flowmag(imMag,vxDFW_sms,vyDFW_sms,vzDFW_sms,vMax,[1,2,3]) title('Div Free Wavelet w/ SureShrink and Partial Cycle Spinning (Vel Mag)','FontSize',14) % Calculate errors disp('DivFree Wavelet w/ SureShrink, MAD and Partial Cycle Spinning') [vNRMSE_DFW_sms,vMagErr_DFW_sms,angErr_DFW_sms] = calcVelError(imMask,vx,vy,vz,vxDFW_sms,vyDFW_sms,vzDFW_sms); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFW_sms,vMagErr_DFW_sms,angErr_DFW_sms); pause(1);
DivFree Wavelet w/ SureShrink, MAD and Partial Cycle Spinning NRMSE: 0.022368 vMag Error: 0.018292 Absolute Angle Error: 0.006091


Finite Difference Method
The following implements finite difference method denoising as described in: Song SM, Pelc NJ., et al. JMRI 1993 Noise reduction in three-dimensional phase-contrast MR velocity measurementsl.
% FDM Denoise [vxFDM,vyFDM,vzFDM] = fdmDenoise(vxN,vyN,vzN,res); % Plot figure,imshow_flow(imMag,vxFDM,vyFDM,vzFDM,vMax,[1,2,3]) title('Div Free FDM','FontSize',14) figure,imshow_flowmag(imMag,vxFDM,vyFDM,vzFDM,vMax,[1,2,3]) title('Div Free FDM','FontSize',14) % Calculate errors disp('Finite Difference Method Error') [vNRMSE_FDM,vMagErr_FDM,angErr_FDM] = calcVelError(imMask,vx,vy,vz,vxFDM,vyFDM,vzFDM); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_FDM,vMagErr_FDM,angErr_FDM); pause(1);
Finite Difference Method Error NRMSE: 0.075755 vMag Error: 0.047040 Absolute Angle Error: 0.080547


DivFree Radial Basis Function
The following implements divergence-free radial basis function denoising as described in: Busch J, Kozerke S., et al. MRM 812 Construction of divergence-free velocity fields from cine 3D phase-contrast flow measurements.
radius = 3; % Radius of kernel nIter = 15; % Number of iterations for lsqr % Plot during iterations, if on, does gradient descent instead of lsqr doplot = 0; % RBF Denoise [vxRBF,vyRBF,vzRBF] = rbfDenoise(vxN,vyN,vzN,imMask,radius,res,nIter,doplot); % Plot figure,imshow_flow(imMag,vxRBF,vyRBF,vzRBF,vMax,[1,2,3]), title('Div Free RBF','FontSize',14) figure,imshow_flowmag(imMag,vxRBF,vyRBF,vzRBF,vMax,[1,2,3]), title('Div Free RBF','FontSize',14) % Calculate Error disp('Radial Basis Function Error') [vNRMSE_RBF,vMagErr_RBF,angErr_RBF] = calcVelError(imMask,vx,vy,vz,vxRBF,vyRBF,vzRBF); fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_RBF,vMagErr_RBF,angErr_RBF); pause(1);
RBF Denoising... done Radial Basis Function Error NRMSE: 0.034087 vMag Error: 0.023883 Absolute Angle Error: 0.037344


Plot all
Plotting results from all methods. Notice how DFW has less errors near edges than other methods
figure,imshow_flowmag(... cat(1,cat(2,imMask,imMask,imMask,imMask),cat(2,imMask,imMask,imMask,imMask)),... cat(1,cat(2,vx,vxN,vxDFW,vxDFW_s),cat(2,vxDFW_sm,vxDFW_sms,vxFDM,vxRBF)),... cat(1,cat(2,vy,vyN,vyDFW,vyDFW_s),cat(2,vyDFW_sm,vyDFW_sms,vyFDM,vyRBF)),... cat(1,cat(2,vz,vzN,vzDFW,vzDFW_s),cat(2,vzDFW_sm,vzDFW_sms,vzFDM,vzRBF)),... vMax,[1,2,3]); title('(Top row) Original,Noisy, DFW/manual, DFW/manual/spin (Bottom Row) DFW/SURE, DFW/SURE/spin, FDM,RBF','FontSize',14) pause(1);

Visualize DivFree Wavelet Coefficients
% Wavelet coefficients for original minSize = 8*ones(1,3); [wcdf1,wcdf2,wcn,numLevels,wcSizes] = dfwavelet_forward(vx,vy,vz,minSize,res); wcdf1_tile = wcTile(wcdf1,numLevels,wcSizes); wcdf2_tile = wcTile(wcdf2,numLevels,wcSizes); wcn_tile = wcTile(wcn,numLevels,wcSizes); % Wavelet coefficients for noisy [wcdf1N,wcdf2N,wcnN,numLevels,wcSizes] = dfwavelet_forward(vxN,vyN,vzN,minSize,res); wcdf1N_tile = wcTile(wcdf1N,numLevels,wcSizes); wcdf2N_tile = wcTile(wcdf2N,numLevels,wcSizes); wcnN_tile = wcTile(wcnN,numLevels,wcSizes); % Wavelet coefficeints for div free wavelet w/ cycle spinning [wcdf1_sms,wcdf2_sms,wcn_sms,numLevels,wcSizes] = dfwavelet_forward(vxDFW_sms,vyDFW_sms,vzDFW_sms,minSize,res); wcdf1_sms_tile = wcTile(wcdf1_sms,numLevels,wcSizes); wcdf2_sms_tile = wcTile(wcdf2_sms,numLevels,wcSizes); wcn_sms_tile = wcTile(wcn_sms,numLevels,wcSizes); % Plot figure,imshow3f(abs(cat(2,wcn_tile,wcnN_tile,wcn_sms_tile)),[1,2,3],round(N/8)),title('Non-divfree wavelet subbands for Original, Noisy and DFW','FontSize',14); pause(1);

Plot Velocity profiles
figure,plot3f(cat(2,vx,vxN,vxDFW_sms),[1,2,3]), title('Velocity profiles for Original,Noisy, and DFW/spin','FontSize',14);
