Denoise Flow Phantom Data with incorrect segmentation

The following compares divergence-free wavelet denoising with existing methods on a flow phantom with incorrect segmentation.

see demo_phantom.m for performances with correct segmentation

Contents

Clear all and check path

close all
clc
clear

if ~exist('dfwavelet_thresh','file')
    error('Cannot find dfwavelet functions. run setPath!');
end

Load phantom data and set parameters

load 4dFlow_007


% Eddy Current Correction
calib = [127-100,128+100;127-40,128+40;32-20,32+20];
[vxEddy,vyEddy,vzEddy] = eddyCorr(vx,vy,vz,calib);
vx = vx-vxEddy;
vy = vy-vyEddy;
vz = vz-vzEddy;


% Crop to pipe for faster runtime
% Comment the following to denoise the entire flow data
imMag = imMag(20:250,1:64,17:49);
vx = vx(20:250,1:64,17:49);
vy = vy(20:250,1:64,17:49);
vz = vz(20:250,1:64,17:49);
vxEddy = vxEddy(20:250,1:64,17:49);
vyEddy = vyEddy(20:250,1:64,17:49);
vzEddy = vzEddy(20:250,1:64,17:49);


% Set parameters
imThresh = 2000;                % Image threshold for segmentation
imMask = (imMag>imThresh)*1;    % Segmentation mask
imCorrThresh = 7000;            % Correct image threshold
imCorrMask = (imMag>imCorrThresh)*1;    % Correct segmentation mask
vMag = getVelMag(vx,vy,vz);     % Velocity magnitude
vMax = 1.0127e+03;              % Maximum speed
FOV = size(vMag);               % Field of view
N = FOV(3);                     % Number of slices (for plotting)
ph0 = zeros(FOV);               % Reference phase in phase contrast


% Plot
figure,imshow_flow(imMask,vx,vy,vz,vMax,[1,2,3]),title('Original Flow Field','FontSize',14);
figure,imshow_flowmag(imMask,vx,vy,vz,vMax,[1,2,3]),title('Original Velocity Magnitude','FontSize',14);
pause(1);

Get Noisy Data

load 4dFlow_007_noisy


% Crop to pipe for faster performance
% Comment the following to denoise the entire flow data
vxN = vxN(20:250,1:64,17:49);
vyN = vyN(20:250,1:64,17:49);
vzN = vzN(20:250,1:64,17:49);


% Eddy Current Correction
vxN = vxN-vxEddy;
vyN = vyN-vyEddy;
vzN = vzN-vzEddy;


% Mask flow data with incorrect segmentation
[vxN,vyN,vzN] = maskIM(imMask,vxN,vyN,vzN);


% Plot
figure,imshow_flow(imMask,vxN,vyN,vzN,vMax,[1,2,3]),title('Noisy Flow Field','FontSize',14);
figure,imshow_flowmag(imMask,vxN,vyN,vzN,vMax,[1,2,3]),title('Noisy Velocity Magnitude','FontSize',14);


% Calculate errors
disp('Noisy Flow Error')
[vNRMSE_Noise,vMagErr_Noise,angErr_Noise] = calcVelError(imCorrMask,vx,vy,vz,vxN,vyN,vzN);
PVNR = 20*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: 			33.44dB
NRMSE: 			0.021276
vMag Error: 		0.012291
Absolute Angle Error: 	0.008179

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

% Denoise
[vxDFWsm,vyDFWsm,vzDFWsm] = dfwavelet_thresh_SURE_MAD(vxN,vyN,vzN,minSize,res);


% Plot
figure,imshow_flow(imMask,vxDFWsm,vyDFWsm,vzDFWsm,vMax,[1,2,3])
title('Div Free Wavelet w/ SureShrink (Flow Field)','FontSize',14)
figure,imshow_flowmag(imMask,vxDFWsm,vyDFWsm,vzDFWsm,vMax,[1,2,3])
title('Div Free Wavelet w/ SureShrink (Vel Mag)','FontSize',14)


% Calculate errors
disp('DivFree Wavelet w/ SureShrink and MAD')
[vNRMSE_DFWsm,vMagErr_DFWsm,angErr_DFWsm] = calcVelError(imCorrMask,vx,vy,vz,vxDFWsm,vyDFWsm,vzDFWsm);
fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFWsm,vMagErr_DFWsm,angErr_DFWsm);
pause(1);
DivFree Wavelet w/ SureShrink and MAD
NRMSE: 			0.017372
vMag Error: 		0.010133
Absolute Angle Error: 	0.007159

DivFree Wavelet with SureShrink, MAD and random 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


% Denoise
[vxDFWsms,vyDFWsms,vzDFWsms] = dfwavelet_thresh_SURE_MAD_spin(vxN,vyN,vzN,minSize,res,spins,isRandShift);


% Plot
figure,imshow_flow(imMask,vxDFWsms,vyDFWsms,vzDFWsms,vMax,[1,2,3])
title('Div Free Wavelet w/ SureShrink and Partial Cycle Spinning (Flow Field)','FontSize',14)
figure,imshow_flowmag(imMask,vxDFWsms,vyDFWsms,vzDFWsms,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_DFWsms,vMagErr_DFWsms,angErr_DFWsms] = calcVelError(imCorrMask,vx,vy,vz,vxDFWsms,vyDFWsms,vzDFWsms);
fprintf('NRMSE: \t\t\t%f\nvMag Error: \t\t%f\nAbsolute Angle Error: \t%f\n\n',vNRMSE_DFWsms,vMagErr_DFWsms,angErr_DFWsms);
pause(1);
DivFree Wavelet w/ SureShrink, MAD and Partial Cycle Spinning
NRMSE: 			0.014120
vMag Error: 		0.008493
Absolute Angle Error: 	0.004810

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.

% Denoise
[vxFDM,vyFDM,vzFDM] = fdmDenoise(vxN,vyN,vzN,res);


% Plot
figure,imshow_flow(imMask,vxFDM,vyFDM,vzFDM,vMax,[1,2,3])
title('Div Free FDM','FontSize',14)
figure,imshow_flowmag(imMask,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(imCorrMask,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.037992
vMag Error: 		0.022113
Absolute Angle Error: 	0.026409

DivFree Radial Basis Function

The following implements divergence-free radial basis function denoising as described in: Busch J, Kozerke S., et al. MRM 2012 Construction of divergence-free velocity fields from cine 3D phase-contrast flow measurements.

radius = 4;     % Radius of kernel
nIter = 20;     % Number of iterations for lsqr

% Plot during iterations, if on, does gradient descent instead of lsqr
doplot = 0;


% Denoise
% Takes a while to run
[vxRBF,vyRBF,vzRBF] = rbfDenoise(vxN,vyN,vzN,imMask,radius,res,nIter,doplot);


% Plot
figure,imshow_flow(imMask,vxRBF,vyRBF,vzRBF,vMax,[1,2,3]),
title('Div Free RBF','FontSize',14)
figure,imshow_flowmag(imMask,vxRBF,vyRBF,vzRBF,vMax,[1,2,3]),
title('Div Free RBF','FontSize',14)


% Calculate errors
disp('Radial Basis Function Error')
[vNRMSE_RBF,vMagErr_RBF,angErr_RBF] = calcVelError(imCorrMask,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.031683
vMag Error: 		0.019790
Absolute Angle Error: 	0.021085

Plot all

Plotting results from all methods.

figure,imshow_flowmag(...
    cat(2,imMask,imMask,imMask,imMask,imMask,imMask),...
    cat(2,vx,vxN,vxDFWsm,vxDFWsms,vxFDM,vxRBF),...
    cat(2,vy,vyN,vyDFWsm,vyDFWsms,vyFDM,vyRBF),...
    cat(2,vz,vzN,vzDFWsm,vzDFWsms,vzFDM,vzRBF),...
    vMax,[1,2,3]);

title('Original,Noisy,DFW,DFW/spin,FDM and RBF','FontSize',14)
pause(1)

Visualize DivFree Wavelet Coefficients

% Wavelet coefficients for original
minSize = 8*ones(1,3);
[vx,vy,vz] = maskIM(imMask,vx,vy,vz);
[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
[wcdf1sms,wcdf2sms,wcnsms,numLevels,wcSizes] = dfwavelet_forward(vxDFWsms,vyDFWsms,vzDFWsms,minSize,res);
wcdf1sms_tile = wcTile(wcdf1sms,numLevels,wcSizes);
wcdf2sms_tile = wcTile(wcdf2sms,numLevels,wcSizes);
wcnsms_tile = wcTile(wcnsms,numLevels,wcSizes);


% Plot
figure,imshow3f(abs(cat(2,wcn_tile,wcnN_tile,wcnsms_tile)),[1,2,3],round(N/4)),title('Non-divfree wavelet subbands for Original, Noisy and DFW','FontSize',14);