Library for MRI noise analysis
Project description
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
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.
Sample generation
python -m mrina.genSamples --fromdir $KSPACEDIR \
--repetitions $REALIZATIONS \
--origin $IMGNAME \
--dest $RECDIR \
--utype $SAMPTYPE \
--urate $PVAL \
--noisepercent $NOISEVAL
For additional information on the script input parameters, type
python -m mrina.gen_samples --help
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
For additional information on the script input parameters, type
python -m mrina.recover --help
Post-processing - Saving reconstructed images
python -m mrina.saveimgs --numsamples $REALIZATIONS \
--maindir $MAINDIR \
--recdir $RECDIR \
--maskdir $MASKDIR \
--outputdir $OUTDIR \
--savetrue \
--savemask \
--saverec \
--savenoise \
--savelin \
--usetrueasref \
--printlevel $PRINTLEV \
--savelin \
--limits $LIMITS \
--fluidmaskfile $FMFILE
For additional information on the script input parameters, type
python -m mrina.saveimgs --help
Post-processing - Computing correlations
python -m mrina.correlation --numsamples $REALIZATIONS \
--numpts 50 \
--recdir ./CS/ \
--ptsdir ./ \
--vencdir ./ \
--maindir ./ \
--usefluidmask \
--printlevel 1
For additional information on the script input parameters, type
python -m mrina.correlation --help
Post-processing - Plot correlations
python -m mrina.plot_corr --noise 0.1 0.01 0.05 0.3 \
--uval 0.75 0.25 0.5 \
--utype vardengauss bernoulli \
--method cs csdebias omp \
--wavelet haar db8 \
--numsamples 100 \
--numpts 50 \
--dir ./ \
--outputdir ./OUT/02_corr/ \
--usefluidmask \
--printlevel 1
For additional information on the script input parameters, type
python -m mrina.plot_corr --help
Post-processing - Compute MSE and relative errors
python -m mrina.plot_mse --noise 0.1 0.01 0.05 0.3 \
--uval 0.75 0.25 0.5 \
--utype vardengauss bernoulli \
--method cs csdebias omp \
--wavelet haar db8 \
--numsamples 100 \
--numpts 50 \
--dir ./ \
--outputdir ./OUT/03_mse/ \
--maskdir ./ \
--usecompleximgs \
--addlinearrec \
--usetrueimg \
--usefluidmask \
--fluidmaskfile ia_mask.npy \
--printlevel 1 \
--percstring 1
For additional information on the script input parameters, type
python -m mrina.plot_mse --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...
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Source Distribution
Built Distributions
Hashes for mrina-0.1.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | 53c4090870b38fb1f8b05f72dc5fd61ea9da961d6f5a2d700fd2411b1f8c71bb |
|
MD5 | e96b508c8feccec23d7c421d9035f2cf |
|
BLAKE2b-256 | 15ccf3923a9dcc86d6efffe57b85561229c71b11feb7195bcf33fc9279b376ca |
Hashes for mrina-0.1.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
Algorithm | Hash digest | |
---|---|---|
SHA256 | f721d4c398920289de281bb958f1a76786474247811b63a98e1b40949fe45591 |
|
MD5 | b8e0c038b502bdcabb354e39072deed3 |
|
BLAKE2b-256 | 4258027a52f1a4cd4b0de3f671a73c9cd6af747436f9da0c8dd049016c2bed16 |