mrina 0.1.0

Library for MRI noise analysis

MRINA: A library for MRI Noise Analysis

Modeling noise random fields generated from [ArXiv]

Lauren Hensley Partin, Daniele E. Schiavazzi and Carlos A. Sing-Long Collao

Result images

---

The mrina library provides functionalities to

---

Examples for single-image processing

Original Image	k-space

Read grayscale image

import cv2
im = cv2.imread('city.png', cv2.IMREAD_GRAYSCALE)/255.0

Read grayscale image

from mrina import generateSamplingMask

# Select an undesampling ratio
delta = 0.25
# Generate an undersampling mask
omega = generateSamplingMask(im.shape, delta, 'bernoulli')
# Verify the undersampling ratio
nsamp = np.sum((omega == 1).ravel())/np.prod(omega.shape)
print('Included frequencies: %.1f%%' % (nsamp*100))

Read grayscale image

import pywt

waveName = 'haar'
waveMode = 'zero'
wim = pywt.coeffs_to_array(pywt.wavedec2(im, wavelet=waveName, mode=waveMode))[0]
plt.figure(figsize=(8,8))
plt.imshow(np.log(np.abs(wim)+1.0e-5), cmap='gray')
plt.axis('off')
plt.show()

Initialize a WaveletToFourier operator and generate undersampled k-space measurements

from maps import OperatorWaveletToFourier

A = OperatorWaveletToFourier(im.shape, samplingSet=omega[0], waveletName=waveName)
yim = A.eval(wim, 1)

Noiseless recovery using l1-norm minimization

from solver_l1_norm import RecoveryL1NormNoisy

# Recovery - for low values of eta it is better to use SoS-L1Ball
wimrec_cpx, _ = RecoveryL1NormNoisy(0.01, yim, A, disp=True, method='SoS-L1Ball')
# The recovered coefficients could be complex!
imrec_cpx = A.getImageFromWavelet(wimrec_cpx)
imrec = np.abs(imrec_cpx)

Generate k-space noise

# Target SNR
SNR = 50
# Signal power. The factor 2 accounts for real/imaginary parts
yim_pow = la.norm(yim.ravel()) ** 2 / (2 * yim.size)
# Noise st. dev.
sigma = np.sqrt(yim_pow / SNR)
# Noisy measurements
y = yim + sigma * (np.random.normal(size=yim.shape) + 1j * np.random.normal(size=yim.shape))

Image recovery with l1-norm minimization

# Parameter eta
eta = np.sqrt(2 * y.size) * sigma
# Recovery
wimrec_noisy_cpx, _ = RecoveryL1NormNoisy(eta, y, A, disp=True, disp_method=False, method='BPDN')
# The recovered coefficients could be complex!
imrec_noisy = np.abs(A.getImageFromWavelet(wimrec_noisy_cpx))

Estimator debiasing

# Support of noisy solution
wim_supp = np.where(np.abs(wimrec_noisy_cpx) > 1E-4 * la.norm(wimrec_noisy_cpx.ravel(), np.inf), True, False)
# Restriction of the operator
Adeb = A.colRestrict(wim_supp)
# Solve least-squares problem
lsqr = lsQR(Adeb)  
lsqr.solve(y[Adeb.samplingSet])
wimrec_noisy_cpx_deb = np.zeros(Adeb.wavShape,dtype=np.complex)
wimrec_noisy_cpx_deb[Adeb.basisSet] = lsqr.x[:]
# The recovered coefficients could be complex!
imrec_noisy_deb = np.abs(Adeb.getImageFromWavelet(wimrec_noisy_cpx_deb))

Image recovery with stOMP

from solver_omp import lsQR,OMPRecovery
# Recovery
wimrec_noisy_cpx, _ = OMPRecovery(A, y)
# The recovered coefficients could be complex!
imrec_noisy_cpx = A.getImageFromWavelet(wimrec_noisy_cpx)
imrec_noisy = np.abs(imrec_noisy_cpx)

Script functionalities

Sample generation

  python -m mrina.genSamples --fromdir KSPACEDIR \
                             --repetitions REALIZATIONS \
                             --origin IMGNAME \
                             --dest RECDIR \
                             --utype SAMPTYPE \
                             --urate PVAL \
                             --noisepercent NOISEVAL

• KSPACEDIR -
• REALIZATIONS -
• IMGNAME -
• RECDIR -
• SAMPTYPE -
• PVAL -
• NOISEVAL -

Image recovery

  python -m mrina.recover --noisepercent $NOISEVAL \
                          --urate $PVAL \
                          --utype $SAMPTYPE \
                          --repetitions $REALIZATIONS \
                          --numprocesses $PROCESSES \
                          --fromdir $KSPACEDIR \
                          --recdir $RECDIR \
                          --maskdir $PATTERNDIR \
                          --method $SOLVERMODE \
                          --wavelet $WAVETYPE \
                          --savevels

Post-processing - Saving reconstructed images

  python -m mrina.saveimgs --numsamples 100 \
                           --maindir ./ \
                           --recdir ./CS/ \
                           --maskdir ./ \
                           --outputdir ./OUT/01_imgs/ \
                           --savetrue \
                           --savemask \
                           --saverec \
                           --savenoise \
                           --savelin \
                           --usetrueasref \
                           --printlevel 1 \
                           --savel in \
                           --limits 0.0 1.0 0.0 0.9733396442983887 0.0 1.0 0.0 0.9733215407741752 \
                           --fluidmaskfile ia_mask.npy

For additional information on the script input parameters, type

python -m mrina.gen_samples --help

Core Dependencies

• Python 3.6.5
• PyWavelets 1.1.1
• Numpy 1.18.1
• Scipy 1.1.0
• Matplotlib 3.1.0
• [Cython]
• [opencv]

Citation

Find this useful or like this work? Cite us with:

Add Paper once published...