# Denoise Undersampled Flow Phantom Data

The following compares divergence-free wavelet denoising with other methods on undersampled data reconstructed by ESPIRiT/SENSE with no regularization. Divergence-free wavelet with manually chosen thresholds is also compared

see demo_phantom.m for denoising performances on fully sampled data

## Contents

- Clear all and check path
- Load phantom data and set parameters
- Load Reconstructed Flow from undersampled data
- DivFree Wavelet with SureShrink and MAD sigma estimation
- DivFree Wavelet with SureShrink, MAD and random cycle spinning
- DivFree Wavelet with SureShrink and manually chosen thresholds
- Finite Difference Method
- DivFree Radial Basis Function
- Plot all

## 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 performance
    % 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 = 7000;                % Image threshold for segmentation
    imMask = (imMag>imThresh)*1;    % 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);

## Load Reconstructed Flow from undersampled data

    load 4dFlow_007_subsample

    % 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 good segmentation
    [vxN,vyN,vzN] = maskIM(imMask,vxN,vyN,vzN);

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

    % Calculate errors
    disp('Reconstructed Flow Error')
    [vNRMSE_Noise,vMagErr_Noise,angErr_Noise] = calcVelError(imMask,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)

    Reconstructed Flow Error
    PVNR: 			29.47dB
    NRMSE: 			0.033631
    vMag Error: 		0.019393
    Absolute Angle Error: 	0.024910

## 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 sureShrink w/ MAD')
    [vNRMSE_DFWsm,vMagErr_DFWsm,angErr_DFWsm] = calcVelError(imMask,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 sureShrink w/ MAD
    NRMSE: 			0.030420
    vMag Error: 		0.017794
    Absolute Angle Error: 	0.021827

## 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(imMask,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.029874
    vMag Error: 		0.017450
    Absolute Angle Error: 	0.021446

## DivFree Wavelet with SureShrink and manually chosen thresholds

    spins = 2;              % Number of cycle spinning per dimension
    isRandShift = 1;        % Use random shift
    minSize = 8*ones(1,3);  % Smallest wavelet level size
    dfThresh = 15;          % Threshold for divfree component
    nThresh = 20;           % Threshold for non-divfree component

    % Denoise
    [vxDFW,vyDFW,vzDFW] = dfwavelet_thresh_spin(vxN,vyN,vzN,minSize,res,dfThresh,nThresh,spins,isRandShift);

    % Plot
    figure,imshow_flow(imMask,vxDFW,vyDFW,vzDFW,vMax,[1,2,3]),
    title('Div Free Wavelet w/ manual threshold and random cycle spinning (Flow Field)','FontSize',15)
    figure,imshow_flowmag(imMask,vxDFW,vyDFW,vzDFW,vMax,[1,2,3]),
    title('Div Free Wavelet w/ manual threshold and random cycle spinning (Vel Mag)','FontSize',15)

    % Calculate errors
    vMagDFW = getVelMag(vxDFW,vyDFW,vzDFW);
    disp('DivFree Wavelet w/ Manual Thresh')
    [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 Thresh
    NRMSE: 			0.023951
    vMag Error: 		0.015891
    Absolute Angle Error: 	0.015785

## 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 measurements.

    % 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(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.038941
    vMag Error: 		0.023989
    Absolute Angle Error: 	0.029318

## 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(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.023357
    vMag Error: 		0.015485
    Absolute Angle Error: 	0.015471

## Plot all

Plotting results from all methods.

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

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