Add Reconstructors and API for easier MR reconstruction

examples/cartesian_reconstruction.py

@@ -17,11 +17,10 @@
 # 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.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
 import numpy as np
@@ -67,29 +66,21 @@
 # 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 reconstructor
+reconstructor = SingleChannelReconstructor(
+    kspace_data=kspace_data,
+    kspace_loc=kspace_loc,
+    uniform_data_shape=fourier_op.shape,
     wavelet_name="sym8",
-    samples=kspace_loc,
-    nb_scales=4,
     mu=2 * 1e-7,
+    nb_scale=4,
     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)
+    gradient_method='synthesis',
+    optimization_alg='fista',
+    verbose=1
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(num_iterations=200)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 # Calculate SSIM

examples/cartesian_reconstruction_3d.py

@@ -17,9 +17,8 @@
 # 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.utils import convert_mask_to_locations
+from mri.reconstructors import SingleChannelReconstructor
 from pysap.data import get_sample_data
 import pysap
 
@@ -70,29 +69,22 @@
 # 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 reconstructor
+reconstructor = SingleChannelReconstructor(
+    kspace_data=kspace_data,
+    kspace_loc=kspace_loc,
+    uniform_data_shape=fourier_op.shape,
     wavelet_name="sym8",
-    samples=kspace_loc,
-    nb_scales=2,
-    mu=2*1e-11,
+    mu=2 * 1e-11,
+    nb_scale=4,
     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)
+    gradient_method='synthesis',
+    optimization_alg='fista',
+    padding_mode='periodization',
+    verbose=1
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(num_iterations=200)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 # Calculate SSIM

examples/non_cartesian_reconstruction.py

@@ -16,11 +16,9 @@
 
 # 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.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
 
@@ -70,29 +68,22 @@
 # 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,
+# Setup the reconstructor
+reconstructor = SingleChannelReconstructor(
+    kspace_data=kspace_obs,
+    kspace_loc=kspace_loc,
+    uniform_data_shape=image.shape,
     wavelet_name="sym8",
-    samples=kspace_loc,
     mu=6 * 1e-7,
-    nb_scales=4,
+    nb_scale=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,
-    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)
+    gradient_method='synthesis',
+    optimization_alg='fista',
+    verbose=1
+)
+# Start Reconstruction
+x_final, costs, metrics = reconstructor.reconstruct(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

@@ -15,11 +15,9 @@
 
 # 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.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
 
@@ -82,29 +80,21 @@
 # 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,
+# Setup the reconstructor
+reconstructor = SingleChannelReconstructor(
+    kspace_data=kspace_obs,
+    kspace_loc=kspace_loc,
+    uniform_data_shape=image.shape,
     wavelet_name="sym8",
-    samples=kspace_loc,
     mu=6 * 1e-9,
-    nb_scales=4,
+    nb_scale=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,
-    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)
+    gradient_method='synthesis',
+    optimization_alg='fista',
+    verbose=1
+)
+x_final, costs, metrics = reconstructor.reconstruct(num_iterations=200)
 image_rec = pysap.Image(data=np.abs(x_final))
 # image_rec.show()
 recon_ssim = ssim(image_rec, image)