Project description

example workflow

MRIna: A library for MRI Noise Analysis

MRIna is a library for the analysis of reconstruction noise for Modeling noise random fields generated from [ArXiv]

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

The complete set of results from the above paper can be found at this link

Installation and documentation

You can install MRIna with pip (link to PyPI)

pip install mrina

For the documentation follow this link.

What you can do with MRIna.

The MRIna library provides the following functionalities.

It supports the generation of k-space undersampling masks of various types including Bernoulli, variable density triangular, variable density Gaussian, variable density exponential and Halton quasi-random sequences.
It supports arbitrary operators that implement a forward call (eval), and inverse call (adjoint), column restriction (colRestrict), shape and norm.
It supports various non-linear reconstruction methods including l1-norm minimization with iterative thresholding and stagewise orthogonal matching pursuit.
It provides a number of scripts to generate ensembles of synthetic, subsampled and noisy k-space images (4 complex images), to reconstruct image density and velocities, and to post-process to compute error patterns, correlation, MSE and relative errors.

Single-image examples

Original image		Wavelet transform
k-space mask		Noisy k-space measurements
Reconstruction: noiseless		Reconstruction: CS
Reconstruction: CSDEB		Reconstruction: stOMP

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 mrina import OperatorWaveletToFourier

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

Noiseless recovery using l1-norm minimization

from mrina 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 mrina 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

MRIna also provides scripts to automate:

the generation of noisy k-space signals.
linear and non-linear image reconstruction.
reconstructed images post-processing.