MBIRJAX v0.3.2 update (#32)

* Convert from cos_sin_angles to angles. * Working version of new install. * Update install to work with readthedocs (I hope). * Update docs and doc installation. * Improve developer docs. * Add gpu memory and timing scripts. * Improve and revise docs * Minor correction * buzzard_dev (#10) * Refactor num_iterations as a method input rather than a class parameter. * Implement voxel_batch_size for forward projection. * Implement voxel_batch_size for back projection. * Add initialization feature to recon * Prepare for merge of add_prox. * Improve error message. --------- Co-authored-by: Charles Bouman <Charles.Bouman@gmail.com> * Add proximal map (#12) * Possibly buggy but functioning commit of refactored VCD functions * Partially tested initial implementation of prox map * Prox bug causes nan when initial condition is perfect * Update gpu evaluatation for new interface. * Add epsilon to avoid alpha=nan. * Pull auto_set_regularization_params() out of prox() * Update demo script --------- Co-authored-by: Greg Buzzard <buzzard@purdue.edu> Co-authored-by: gbuzzard <54102356+gbuzzard@users.noreply.github.com> * Remove unused attributes: proximal_map_flag, max_resolutions, initialization * Add batch size management to recon_test.py * Add advanced features to documentation * Fix typo. * Add prox_map to docs * Convert to use public instance variables in Projectors. * Refactor to use recon_shape. * Refactor to use recon_shape and projector_params. * Remove prox_recon parameter * Improve documentation of Projectors. * Improve documentation of Projectors. * Basic Conebeam (#16) * Refactor parallel_beam.py for clarity * Remove references to svmbir. * More renaming in parallel_beam.py * Add partial implementation of cone_beam.py * Update incomplete version of cone_beam.py * Updates to incomplete version cone beam projector * Updates to incomplete version cone beam projector * Progress towards working version. * More progress towards working version. * More progress towards working version. * More progress towards working version. * More progress towards working version. * More progress towards working version. * Working version - still needs to be tested in multiple cases. * Add cone beam demo and experiment scripts * Refactor per-pixel magnification to allow for infinite source-detector distance. * Include more options to exercise conebeam model. * Remove reshape from test and demo scripts * Bug fix to get the correct system matrix row entries. * Add link to powerpoint slides. * Rename variables and functions for clarity * Rename various variables and methods. * Rename voxel to pixel as needed. * Fix naming problems and restore vcd_figs_for_abst.py to working condition. * Minor improvement to vcd_figs_for_abst.py script * Add variable intensity window to slice_viewer and add 'modified' to the function that generates the modified Shepp-Logan phantom. * Partial progress towards improving efficiency in conebeam. * Update cvpr scripts * Partial progress towards improving efficiency in conebeam. * Rename delta_voxel_xy, delta_voxel_z, and delta_pixel_recon to delta_voxel. * First implementation of forward vertical fan beam projection. * Fix bugs in sphinx docs * Add TomographyModel parameter definitions * Refactor sphinx docs names * Add documentation for cone beam class * Fix bug in conebeam backprojection and simplify parameter handling. --------- Co-authored-by: gbuzzard <54102356+gbuzzard@users.noreply.github.com> Co-authored-by: Greg Buzzard <buzzard@purdue.edu> * Efficient Cone Beam and Parameter Handling (#18) * Partial implementation of horizontal fan beam projection. * Add cone beam developer sphinx doc * Working refactored forward projection for cone beam. * Remove unused code. * Change from magnification to source_iso_dist. * Refactor magnification from a parameter to a function. * Additional comment. * Working version of vertical fan beam back projection. * Working version of full refactored back projection. * Working version of full refactored back projection. * Fix Hessian calculation for cone beam. * Add cone beam figure * Remove logo * Add logo * Add pdf, jpg, pptx to gitignore * Fix bug in docs. * Change interface for basic geometry class to operate on a batch of pixels and one view for forward and back projections. Improve docstrings. * Refactor to vmap over pixels in vertical back projection. * Move p to the list of geometry parameters. * Include default batch sizes. Are you happy? * Refactor top level data flow in forward projector. * Add projectors. Yes :-) * Add cone beam tests and pytest support * Refactor top level data flow in back projector. * Add vmap over slices in horizontal fan beam. * Add warnings when tests fail. * Update pixel and voxel batch sizes and handling. Make source_detector_distance a multiple of num_det_channels. * Change auto_set_recon_size to be geometry specific. * Refactored auto_set_sigma functions * Refactor TomographModel to put parameter handling in a superclass. Perform some code cleanup. * Fix docs to reflect ParameterHandler. * Change p to psf_radius and include sinogram viewer. * Add psf_radius support * Add psf_radius for both parallel and cone geometry * Add ability to save and load to file. * Check that loaded parameters match the type of the loading class. * Update documentation. * Edit documentation --------- Co-authored-by: Greg Buzzard <buzzard@purdue.edu> * Add bug-fixes and PyPI support (#20) * Refactor to separate out the calculation of quantities needed for projection. Improve documentation. * Provide a display default when vmin=vmax. * Provide a display default when vmin=vmax. * Allow the slice label to be set. * Fix bug to prevent zooming in sliders. * Add cpu memory stats. * Refactor for joint pan and zoom with two images. * Minor modification to pyproject.toml * Add PyPI dependencies --------- Co-authored-by: Greg Buzzard <buzzard@purdue.edu> * Update step size in vcd and other various improvements. (#26) * Update PyPI instructions. * Correct docstring, tweak auto_regularization, improve demos, and update version number * Update gitignore for jupyter notebooks and add LLNL experiments directory * Minor bug fixes. * Add files to reconstruct and view nersc data. * Fix bug in backproject and improve parameter return for recon. * Fix bug in default weights. * Fix step size in vcd updates. Clean up docs. Change default partitions. Fix bug in row and channel offets in cone and parallel. Combine 2 demos into one. Include more tests. Fix bug in default weights. Use reflected boundary conditions for qggmrf top and bottom slices. * Fix minor bug. * Update code for cvpr abstract to use the original partition sequences. * Remove lbnl files for PR. * Update version number --------- Co-authored-by: Charles Bouman <Charles.Bouman@gmail.com> * Improve package maintenance instructions. * Refactor get_delta as a separate function and include a test for it. * Refactor to use apply_map_in_batches. * Refactor to use sum_function_in_batches and rename apply_map_in batches to concatenate_function_in_batches. * Label regularization parameters as static for more efficient compiling and execution. * Reorganize qggmrf names and functions in preparation for refactoring to use vmap. * Fix error in doc related to recon_shape * Add parameter checking to raise an exception if the user tries to set an invalid parameter. * Add parameter checking to raise an exception if the user tries to set an invalid parameter. * Relax the tolerances a little to account for randomness in vcd partition selection. * Add cube phantom. * Refactor qggmrf gradient and hessian to use vmap. * Update qggmrf test. * Remove randomness to ensure reproducibility. * Partially working version - not fully debugged: Refactor sum and concatenate in batches to work on tuples of inputs and outputs and include tests of these functions. Implement infrastructure to return prior cost and to choose optimal alpha using prior gradient and hessian. * Add QR code generation script * Update tests. * Include cost function for full recon and test for qggmrf gradient correctness. * Start to remove cost function. * Finish refactoring prior cost. * Change default number of iterations. * Update cvpr figures. * Include test for surrogate hessian. * Update displayed total loss. * Rename cost to loss, sigma_p to sigma_prox. * Move qggmrf functions to separate file. * Move qggmrf functions to separate file. * Working version of alpha using approximate prior quadratic. * Remove deprecated keywords. * Change version number. * Buzzard dev (#31) * Fix transparancy of corners of QR code * Put names on variables. * Add spectral plots for cvpr. * Update docs to remove **kwargs from constructor. * Include blue noise sampling and updated partition sequence for demo. * Change default granularity sequence and handle case of more subsets than pixels. * Allow display of 2D images. * Minor changes to experiments. * Update and correct documentation * Update docs and output messages. * Update docs and output messages. * Make random samples the default and blue noise the alternate. --------- Co-authored-by: Charles Bouman <Charles.Bouman@gmail.com> * Include nrms update values. * Include nrms update values. * Update version number to v0.3.2 --------- Co-authored-by: Greg Buzzard <buzzard@purdue.edu> Co-authored-by: gbuzzard <54102356+gbuzzard@users.noreply.github.com>
cabouman · Jun 24, 2024 · 0e75606 · 0e75606
1 parent dc6ea60
commit 0e75606
Show file tree

Hide file tree

Showing 14 changed files with 766 additions and 86 deletions.
diff --git a/dev_scripts/QR-code/create-QR-code.py b/dev_scripts/QR-code/create-QR-code.py
@@ -20,10 +20,10 @@
 qr.make(fit=True)
 
 # Create an image from the QR code instance
-img_qr = qr.make_image(fill='black', back_color='white').convert('RGB')
+img_qr = qr.make_image(fill='black', back_color='white').convert('RGBA')
 
 # Open the logo image
-logo = Image.open(logo_path)
+logo = Image.open(logo_path).convert('RGBA')
 
 # Resize the logo
 logo_width = 400
@@ -37,11 +37,11 @@
 total_height = qr_height + logo_height + padding  # Add some space between QR and logo
 
 # Create a new image with an off-white background
-bg_color = (235, 235, 230)  # Off-white
-final_image = Image.new('RGB', (qr_width, total_height), bg_color)
+bg_color = (235, 235, 230, 255)  # Off-white with full opacity
+final_image = Image.new('RGBA', (qr_width, total_height), bg_color)
 
 # Paste the QR code onto the final image
-final_image.paste(img_qr, (0, 0))
+final_image.paste(img_qr, (0, 0), mask=img_qr)
 
 # Paste the logo onto the final image
 logo_pos = ((qr_width - logo_width) // 2, qr_height + padding)  # Reduced padding
@@ -54,7 +54,7 @@
 draw.rounded_rectangle([(0, 0), final_image.size], radius, fill=255)
 
 # Apply rounded corners mask to the final image
-rounded_final_image = Image.new('RGB', final_image.size)
+rounded_final_image = Image.new('RGBA', final_image.size, (0, 0, 0, 0))
 rounded_final_image.paste(final_image, (0, 0), mask=mask)
 
 # Save the final image

diff --git a/docs/source/_static/new_model_template.py b/docs/source/_static/new_model_template.py
@@ -11,7 +11,12 @@ class TemplateModel(TomographyModel):
     This class inherits all methods and properties from the :ref:`TomographyModelDocs` and may override some
     to suit the specific geometry.
 
+    Parameters not included in the constructor can be set using the set_params method of :ref:`TomographyModelDocs`.
+    Refer to :ref:`TomographyModelDocs` documentation for a detailed list of possible parameters.
+
     Items to change for a particular geometry are highlighted with TODO.
+    It is not recommended to use **kwargs in the constructor since doing so complicates checking for valid parameter
+    names in set_params.
 
     Args:
         sinogram_shape (tuple):
@@ -22,10 +27,6 @@ class TemplateModel(TomographyModel):
             These are view-independent scalar parameters that are required for the geometry and are not already included in the parent class.
         view_dependent_vec1, view_dependent_vec2 (jnp.ndarray):
             These are view-dependent parameter vectors each with length = number of views.
-        **kwargs (dict):
-            Additional keyword arguments that are passed to the :ref:`TomographyModelDocs` constructor. These can
-            include settings and configurations specific to the tomography model such as noise models or image dimensions.
-            Refer to :ref:`TomographyModelDocs` documentation for a detailed list of possible parameters.
     """
 
     # TODO: Adjust the signature as needed for a particular geometry and update the docstring to match.

diff --git a/docs/source/usr_parameters.rst b/docs/source/usr_parameters.rst
@@ -1,35 +1,31 @@
 .. _ParametersDocs:
 
 
-==================
-Primary Parameters
-==================
+===============
+Base Parameters
+===============
 
-The ``TomographyModel`` provides the basic interface for all specific geometries for tomographic projection
-and reconstruction.
+The following documents the base parameters used by the :ref:`TomographyModelDocs` class.
+Any of these parameters can be modified with :func:`TomographyModel.set_params`.
 
-Parameter Documentation
------------------------
+Parameters that are specific to particular geometries are documented in the geometry's documentation.
 
-The following documents basic TomographyModel class parameters that are commonly used in reconstruction.
-Other parameters may be used for specific geometries.
-
-
-Basic Reconstruction Parameters
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Reconstruction Parameters
+^^^^^^^^^^^^^^^^^^^^^^^^^
 
 sharpness : float
     Specifies the sharpness of the reconstruction. Defaults to 0.0. Larger values produce sharper images. Smaller values produce softer images.
 
-snr_db : float
-    Specifies the assumed SNR of sinogram measurements in dB. Defaults to 30.0. Larger values produce sharper images.
-
 verbose : int
     Larger values produce more status information. Defaults to 0 for silent operation.
 
 
-Basic Geometry Parameters
-^^^^^^^^^^^^^^^^^^^^^^^^^
+Geometry Parameters
+^^^^^^^^^^^^^^^^^^^
+
+recon_shape : tuple (num_rows, num_cols, num_slices)
+    Array size of reconstruction. This is set automatically and is available from :meth:`get_params('recon_shape')`.
+    It is recommended to use :func:`set_params` to increase this by a factor of 10-15% when the object extends beyond the field of view.
 
 delta_det_channel : float
     Spacing between detector channels in ALU. Defaults to 1.0.
@@ -40,8 +36,8 @@ delta_det_row : float
 det_channel_offset : float
     Assumed offset between center of rotation and center of detector between detector channels in ALU. Defaults to 0.0.
 
-magnification : float
-    Ratio of (source to detector distance)/(source to iso distance). Defaults to 1.0.
+det_row_offset : float
+    Assumed offset in rows of the source-to-detector line with center of detector in ALU. Defaults to 0.0.
 
 delta_voxel : float
     Spacing between voxels in ALU. Defaults to 1.0.
@@ -57,15 +53,5 @@ sigma_p : float
     Proximal map parameter. Defaults to 1.0.
 
 
-Memory Allocation Parameters
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-pixel_batch_size : int
-    Maximum number of pixels (i.e., voxel cylinders) processed simultaneously.
-
-view_batch_size : int
-    Maximum number of views processed simultaneously.
-
-
 
 
diff --git a/experiments/cvpr-2024/spectral_convergence.py b/experiments/cvpr-2024/spectral_convergence.py
@@ -0,0 +1,94 @@
+import numpy as np
+import jax.numpy as jnp
+import matplotlib.pyplot as plt
+import mbirjax
+
+if __name__ == "__main__":
+    """
+    This is a script to develop, debug, and tune the vcd reconstruction with a parallel beam projector
+    """
+    # Set parameters
+    num_views = 256
+    num_det_rows = 1
+    num_det_channels = 256
+    start_angle = 0
+    end_angle = np.pi
+    sharpness = 0.0
+    num_phantom_slices = 10
+    phantom_shape = (num_det_channels, num_det_channels, num_phantom_slices)
+
+    # Initialize sinogram
+    sinogram_shape = (num_views, num_det_rows, num_det_channels)
+    angles = jnp.linspace(start_angle, np.pi, num_views, endpoint=False)
+
+    # Set up parallel beam model
+    parallel_model = mbirjax.ParallelBeamModel(sinogram_shape, angles)
+
+    # Generate 3D Shepp Logan phantom
+    phantom = mbirjax.generate_3d_shepp_logan_low_dynamic_range(phantom_shape)
+    center_slice = num_phantom_slices // 2
+    phantom = phantom[:, :, center_slice:center_slice+1]
+
+    # Generate synthetic sinogram data
+    print('Creating sinogram')
+    sinogram = parallel_model.forward_project(phantom)
+
+    # Generate weights array
+    weights = parallel_model.gen_weights(sinogram / sinogram.max(), weight_type='transmission_root')
+
+    # Set reconstruction parameter values
+    parallel_model.set_params(sharpness=sharpness, verbose=1)
+    granularity = np.array([1, 2, 64, 512])
+    parallel_model.set_params(granularity=granularity)
+
+    # vcd_partition_sequence = [0, 1, 2, 3, 1, 2, 3, 2, 3, 3, 0, 1, 2, 3, 1, 2, 3, 2, 3, 3]
+    vcd_partition_sequence = [0, 1, 2, 3, 3, 2, 2, 2, 3, 3]  # 2, 3, 2, 3, 3, 0, 1, 2, 3, 1, 2, 3, 2, 3, 3]
+    gd_partition_sequence =  [0, ]
+    icd_partition_sequence = [3, ]
+
+    num_iterations = 5
+
+    sequences = [vcd_partition_sequence, gd_partition_sequence, icd_partition_sequence]
+    residual_fft = []
+    residual_recon = []
+    all_recons = []
+
+    # ##########################
+    # Perform reconstructions
+    for sequence in sequences:
+        parallel_model.set_params(partition_sequence=sequence)
+        cur_residual = np.zeros((num_iterations,) + phantom.shape[0:2])
+        cur_recon = np.zeros((num_iterations,) + phantom.shape[0:2])
+        recon = None
+        for iteration in range(num_iterations):
+            recon, cur_recon_params = parallel_model.recon(sinogram, weights=weights, num_iterations=iteration + 1,
+                                                               first_iteration=iteration, init_recon=recon,
+                                                               compute_prior_loss=True)
+            cur_residual[iteration] = recon[:, :, 0] - phantom[:, :, 0]
+            cur_recon[iteration] = recon[:, :, 0]
+
+        all_recons.append(cur_recon)
+        residual_recon.append(cur_residual)
+        cur_residual_fft = np.fft.fftn(cur_residual, axes=(1, 2))
+        cur_residual_fft = np.fft.fftshift(cur_residual_fft, axes=(1, 2))
+        residual_fft.append(20 * np.log10(1e-6 + np.abs(cur_residual_fft)))
+
+    vcd_recons, gd_recons, icd_recons = all_recons
+    vcd_residual, gd_residual, icd_residual = residual_recon
+    vcd_residual_fft, gd_residual_fft, icd_residual_fft = residual_fft
+
+    mbirjax.slice_viewer(vcd_recons, gd_recons, slice_axis=0, vmin=0, vmax=1, title='Recon, vcd left, gd right', slice_label='Iteration')
+    mbirjax.slice_viewer(vcd_residual, gd_residual, slice_axis=0, vmin=-1, vmax=1, title='Residual recon, vcd left, gd right', slice_label='Iteration')
+    mbirjax.slice_viewer(vcd_residual_fft, gd_residual_fft, title='Residual |FFT| in dB\nVCD left, GD right', slice_axis=0, slice_label='Iteration', vmin=0, vmax=60)
+
+    mbirjax.slice_viewer(vcd_recons, icd_recons, slice_axis=0, vmin=0, vmax=1, title='Recon, vcd left, icd right', slice_label='Iteration')
+    mbirjax.slice_viewer(vcd_residual, icd_residual, slice_axis=0, vmin=-1, vmax=1, title='Residual recon, vcd left, icd right', slice_label='Iteration')
+    mbirjax.slice_viewer(vcd_residual_fft, icd_residual_fft, title='Residual |FFT| in dB\nVCD left, ICD right', slice_axis=0, slice_label='Iteration', vmin=0, vmax=60)
+
+    # fm_rmse_vcd = cur_recon_params.fm_rmse
+    # prior_loss_vcd = cur_recon_params.prior_loss
+    # default_partition_sequence = parallel_model.get_params('partition_sequence')
+    # partition_sequence = mbirjax.gen_partition_sequence(default_partition_sequence, num_iterations=num_iterations)
+    # granularity_sequence_vcd = granularity[partition_sequence]
+
+    a = 0
diff --git a/experiments/cvpr-2024/spectral_response.py b/experiments/cvpr-2024/spectral_response.py
@@ -0,0 +1,126 @@
+import numpy as np
+import jax.numpy as jnp
+import matplotlib.pyplot as plt
+import mbirjax
+import mbirjax.parallel_beam
+from scipy.sparse.linalg import svds, aslinearoperator, LinearOperator
+import jax
+
+if __name__ == "__main__":
+    """
+    This is a script to develop, debug, and tune the vcd reconstruction with a parallel beam projector
+    """
+    view_batch_size = None
+    pixel_batch_size = None
+    with jax.experimental.enable_x64(True):  # Finite difference requires 64 bit arithmetic
+
+        # Initialize sinogram
+        num_views = 128
+        num_det_rows = 1
+        num_det_channels = 128
+        start_angle = 0
+        end_angle = jnp.pi
+        sinogram = jnp.zeros((num_views, num_det_rows, num_det_channels))
+        angles = jnp.linspace(start_angle, jnp.pi, num_views, endpoint=False)
+
+        # Initialize a random key
+        seed_value = np.random.randint(1000000)
+        key = jax.random.PRNGKey(seed_value)
+
+        # Set up parallel beam model
+        parallel_model = mbirjax.ParallelBeamModel(sinogram.shape, angles)
+
+        # Generate phantom
+        recon_shape = parallel_model.get_params('recon_shape')
+        num_recon_rows, num_recon_cols, num_recon_slices = recon_shape[:3]
+        phantom = mbirjax.gen_cube_phantom(recon_shape)
+
+        # Generate indices of pixels
+        num_subsets = 1
+        full_indices = mbirjax.gen_pixel_partition(recon_shape, num_subsets)[0]
+
+        # Generate sinogram data
+        voxel_values = phantom.reshape((-1,) + recon_shape[2:])[full_indices]
+
+        parallel_model.set_params(view_batch_size=view_batch_size, pixel_batch_size=pixel_batch_size)
+
+        print('Starting forward projection')
+        sinogram = parallel_model.sparse_forward_project(voxel_values, full_indices)
+
+        # Determine resulting number of views, slices, and channels and image size
+        print('Sinogram shape: {}'.format(sinogram.shape))
+
+        # Get the vector of indices
+        indices = jnp.arange(num_recon_rows * num_recon_cols)
+
+        sinogram = jnp.array(sinogram)
+        indices = jnp.array(indices)
+
+        # Run once to finish compiling
+        print('Starting back projection')
+        bp = parallel_model.sparse_back_project(sinogram, indices)
+        print('Recon shape: ({}, {}, {})'.format(num_recon_rows, num_recon_cols, num_recon_slices))
+
+        # ##########################
+        # Test the adjoint property
+        # Get a random 3D phantom to test the adjoint property
+        key, subkey = jax.random.split(key)
+        x = jax.random.uniform(subkey, shape=bp.shape)
+        key, subkey = jax.random.split(key)
+        y = jax.random.uniform(subkey, shape=sinogram.shape)
+
+        # Do a forward projection, then a backprojection
+        x = x.reshape((-1, num_recon_slices))[indices]
+        Ax = parallel_model.sparse_forward_project(x, indices)
+        Aty = parallel_model.sparse_back_project(y, indices)
+
+        # Calculate <Aty, x> and <y, Ax>
+        Aty_x = jnp.sum(Aty * x)
+        y_Ax = jnp.sum(y * Ax)
+
+        assert(np.allclose(Aty_x, y_Ax))
+        print("Adjoint property holds for random x, y <y, Ax> = <Aty, x>: {}".format(np.allclose(Aty_x, y_Ax)))
+
+        def Ax_flat(local_x):
+            local_x = np.reshape(local_x, x.shape)
+            ax_flat = parallel_model.sparse_forward_project(local_x, indices)
+            ax_flat = np.array(ax_flat.flatten())
+            return ax_flat
+
+        assert(np.allclose(Ax.flatten(), Ax_flat(x.flatten())))
+        print('Ax_flat matches forward projection')
+
+        def Aty_flat(y):
+            local_y = np.reshape(y, Ax.shape)
+            aty_flat = parallel_model.sparse_back_project(local_y, indices)
+            aty_flat = np.array(aty_flat.flatten())
+            return aty_flat
+
+        assert(np.allclose(Aty.flatten(), Aty_flat(y.flatten())))
+        print('Aty_flat matches back projection')
+
+        Ax_linear_operator = LinearOperator(matvec=Ax_flat, rmatvec=Aty_flat, shape=(Ax.size, x.size))
+
+        Ax_lo = Ax_linear_operator(x.flatten())
+        Aty_lo = Ax_linear_operator.rmatvec(np.array(y).flatten())
+
+        assert(np.allclose(Ax.flatten(), Ax_lo))
+        assert(np.allclose(Aty.flatten(), Aty_lo))
+        print('Linear operator matches known projectors')
+
+        num_sing_values = 15  # num_views * num_det_channels
+        sing_vects = True
+        if sing_vects:
+            u, s, vh = svds(Ax_linear_operator, k=num_sing_values, tol=1e-6, return_singular_vectors=True, solver='propack')
+            vh = vh.reshape(num_sing_values, num_det_channels, num_det_channels)
+            vh = vh[::-1, :, :]
+            mbirjax.slice_viewer(vh, slice_axis=0)
+            u = u.reshape((num_det_channels, num_det_channels, num_sing_values))
+            u = u[:, :, ::-1]
+            mbirjax.slice_viewer(u, slice_axis=2)
+        else:
+            s = svds(Ax_linear_operator, k=num_sing_values, tol=1e-6, return_singular_vectors=False,
+                     solver='propack')
+        plt.plot(np.sort(s)[::-1], '.')
+        plt.show()
+        a = 0