Add Reconstructors and API for easier MR reconstruction

.circleci/config.yml

@@ -28,20 +28,25 @@ jobs:
           command: |
             python3 -m venv venv
             . venv/bin/activate
+            pushd ../
+            git clone https://github.com/CEA-COSMIC/pysap
+            cd pysap
+            python setup.py install
+            echo "export PATH=$PATH:$PWD/build/temp.linux-x86_64-3.6/extern/bin" >> $BASH_ENV
+            popd
             pip install -r doc/requirements.txt
+            pip list > doc/state.txt
       - save_cache:
           paths:
             - ./venv
-          key: v1-dependencies-{{ checksum "doc/requirements.txt" }}
+          key: v1-dependencies-{{ checksum "doc/state.txt" }}
       - run:
           name: Running CircleCI tests
           command: |
             . venv/bin/activate
             pip install Cython
-            pip install git+https://github.com/CEA-COSMIC/pysap.git
             pip install -e .
             python .circleci/test_script.py
       - store_artifacts:
           path: test-reports
           destination: test-reports
-

examples/cartesian_reconstruction.py

@@ -15,15 +15,16 @@
 """
 
 # Package import
-from modopt.math.metrics import ssim
-from mri.operators import FFT
-from mri.optimizers import fista
-from mri.reconstruct.utils import generate_operators
+from mri.operators import FFT, WaveletN
 from mri.operators.utils import convert_mask_to_locations
-from pysap.data import get_sample_data
+from mri.reconstructors import SingleChannelReconstructor
 import pysap
+from pysap.data import get_sample_data
 
 # Third party import
+from modopt.opt.proximity import SparseThreshold
+from modopt.opt.linear import Identity
+from modopt.math.metrics import ssim
 import numpy as np
 
 # Loading input data
@@ -67,29 +68,23 @@
 # We now want to refine the zero order solution using a FISTA optimization.
 # The cost function is set to Proximity Cost + Gradient Cost
 
-# Generate operators
-gradient_op, linear_op, prox_op, cost_op = generate_operators(
-    data=kspace_data,
-    wavelet_name="sym8",
-    samples=kspace_loc,
-    nb_scales=4,
-    mu=2 * 1e-7,
-    fourier_type='cartesian',
-    uniform_data_shape=None,
-    gradient_space="synthesis",
-    padding_mode="periodization")
-
-# Start the FISTA reconstruction
-max_iter = 200
-x_final, costs, metrics = fista(
-    gradient_op,
-    linear_op,
-    prox_op,
-    cost_op,
-    lambda_init=1,
-    max_nb_of_iter=max_iter,
-    atol=1e-4,
-    verbose=1)
+# Setup the operators
+linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
+regularizer_op = SparseThreshold(Identity(), 2 * 1e-7, thresh_type="soft")
+# Setup Reconstructor
+reconstructor = SingleChannelReconstructor(
+    fourier_op=fourier_op,
+    linear_op=linear_op,
+    regularizer_op=regularizer_op,
+    gradient_formulation='synthesis',
+    verbose=1,
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(
+    kspace_data=kspace_data,
+    optimization_alg='fista',
+    num_iterations=200,
+)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 # Calculate SSIM

examples/cartesian_reconstruction_3d.py

@@ -16,14 +16,15 @@
 
 # Package import
 from modopt.math.metrics import ssim
-from mri.operators import FFT
-from mri.optimizers import fista
-from mri.reconstruct.utils import generate_operators
+from mri.operators import FFT, WaveletN
 from mri.operators.utils import convert_mask_to_locations
+from mri.reconstructors import SingleChannelReconstructor
 from pysap.data import get_sample_data
 import pysap
 
 # Third party import
+from modopt.opt.linear import Identity
+from modopt.opt.proximity import SparseThreshold
 import numpy as np
 
 # Loading input data and convert it into a single channel using Sum-Of-Squares
@@ -70,29 +71,28 @@
 # We now want to refine the zero order solution using a FISTA optimization.
 # The cost function is set to Proximity Cost + Gradient Cost
 
-# Generate operators
-gradient_op, linear_op, prox_op, cost_op = generate_operators(
-    data=kspace_data,
+# Setup the operators
+linear_op = WaveletN(
     wavelet_name="sym8",
-    samples=kspace_loc,
-    nb_scales=2,
-    mu=2*1e-11,
-    fourier_type='cartesian',
-    uniform_data_shape=None,
-    gradient_space="synthesis",
-    padding_mode="periodization")
-
-# Start the FISTA reconstruction
-max_iter = 200
-x_final, costs, metrics = fista(
-    gradient_op,
-    linear_op,
-    prox_op,
-    cost_op,
-    lambda_init=1,
-    max_nb_of_iter=max_iter,
-    atol=1e-4,
-    verbose=1)
+    nb_scales=4,
+    dim=3,
+    padding_mode="periodization",
+)
+regularizer_op = SparseThreshold(Identity(), 2 * 1e-11, thresh_type="soft")
+# Setup Reconstructor
+reconstructor = SingleChannelReconstructor(
+    fourier_op=fourier_op,
+    linear_op=linear_op,
+    regularizer_op=regularizer_op,
+    gradient_formulation='synthesis',
+    verbose=1,
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(
+    kspace_data=kspace_data,
+    optimization_alg='fista',
+    num_iterations=200,
+)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 # Calculate SSIM

examples/non_cartesian_reconstruction.py

@@ -15,17 +15,17 @@
 """
 
 # Package import
-from mri.operators import NonCartesianFFT
-from mri.optimizers import fista
-from mri.reconstruct.utils import generate_operators
-from mri.operators.utils import convert_locations_to_mask
-from mri.parallel_mri.extract_sensitivity_maps import \
+from mri.operators import NonCartesianFFT, WaveletUD2
+from mri.operators.utils import convert_locations_to_mask, \
     gridded_inverse_fourier_transform_nd
+from mri.reconstructors import SingleChannelReconstructor
 import pysap
 from pysap.data import get_sample_data
 
 # Third party import
 from modopt.math.metrics import ssim
+from modopt.opt.linear import Identity
+from modopt.opt.proximity import SparseThreshold
 import numpy as np
 
 # Loading input data
@@ -70,29 +70,25 @@
 # We now want to refine the zero order solution using a FISTA optimization.
 # The cost function is set to Proximity Cost + Gradient Cost
 
-# Generate operators
-gradient_op, linear_op, prox_op, cost_op = generate_operators(
-    data=kspace_obs,
-    wavelet_name="sym8",
-    samples=kspace_loc,
-    mu=6 * 1e-7,
-    nb_scales=4,
-    fourier_type='non-cartesian',
-    nfft_implementation='cpu',
-    uniform_data_shape=image.shape,
-    gradient_space="synthesis")
-
-# Start the FISTA reconstruction
-max_iter = 200
-x_final, costs, metrics = fista(
-    gradient_op=gradient_op,
+# Setup the operators
+linear_op = WaveletUD2(
+    wavelet_id=24,
+    nb_scale=4,
+)
+regularizer_op = SparseThreshold(Identity(), 6 * 1e-7, thresh_type="soft")
+# Setup Reconstructor
+reconstructor = SingleChannelReconstructor(
+    fourier_op=fourier_op,
     linear_op=linear_op,
-    prox_op=prox_op,
-    cost_op=cost_op,
-    lambda_init=1.0,
-    max_nb_of_iter=max_iter,
-    atol=1e-4,
-    verbose=1)
+    regularizer_op=regularizer_op,
+    gradient_formulation='synthesis',
+    verbose=1,
+)
+x_final, costs, metrics = reconstructor.reconstruct(
+    kspace_data=kspace_obs,
+    optimization_alg='fista',
+    num_iterations=200,
+)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 recon_ssim = ssim(image_rec, image)

examples/non_cartesian_reconstruction_stack_3d.py

@@ -14,17 +14,17 @@
 """
 
 # Package import
-from mri.operators import Stacked3DNFFT
-from mri.optimizers import fista
-from mri.reconstruct.utils import generate_operators
-from mri.operators.utils import convert_locations_to_mask, get_stacks_fourier
-from mri.parallel_mri.extract_sensitivity_maps import \
-    gridded_inverse_fourier_transform_stack
+from mri.operators import Stacked3DNFFT, WaveletN
+from mri.operators.utils import convert_locations_to_mask, \
+    gridded_inverse_fourier_transform_stack, get_stacks_fourier
+from mri.reconstructors import SingleChannelReconstructor
 import pysap
 from pysap.data import get_sample_data
 
 # Third party import
 from modopt.math.metrics import ssim
+from modopt.opt.linear import Identity
+from modopt.opt.proximity import SparseThreshold
 import numpy as np
 
 # Loading input data
@@ -82,29 +82,23 @@
 # The cost function is set to Proximity Cost + Gradient Cost
 
 # TODO get the right mu operator
-# Generate operators
-gradient_op, linear_op, prox_op, cost_op = generate_operators(
-    data=kspace_obs,
-    wavelet_name="sym8",
-    samples=kspace_loc,
-    mu=6 * 1e-9,
-    nb_scales=4,
-    fourier_type='stack',
-    nfft_implementation='cpu',
-    uniform_data_shape=image.shape,
-    gradient_space="synthesis")
-
-# Start the FISTA reconstruction
-max_iter = 10
-x_final, costs, metrics = fista(
-    gradient_op=gradient_op,
+# Setup the operators
+linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
+regularizer_op = SparseThreshold(Identity(), 6 * 1e-9, thresh_type="soft")
+# Setup Reconstructor
+reconstructor = SingleChannelReconstructor(
+    fourier_op=fourier_op,
     linear_op=linear_op,
-    prox_op=prox_op,
-    cost_op=cost_op,
-    lambda_init=1.0,
-    max_nb_of_iter=max_iter,
-    atol=1e-4,
-    verbose=1)
+    regularizer_op=regularizer_op,
+    gradient_formulation='synthesis',
+    verbose=1,
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(
+    kspace_data=kspace_obs,
+    optimization_alg='fista',
+    num_iterations=10,
+)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 recon_ssim = ssim(image_rec, image)