Merge pull request #233 from cabouman/prerelease

Release 0.2.8, adding support for fan beam

Author: cabouman

README.md

@@ -1,7 +1,7 @@
 # svmbir
 
-*Python code for fast parallel-beam MBIR (Model Based Iterative Reconstruction)  
-This is a python wrapper around HPImaging's supervoxel C code, [sv-mbirct](https://github.com/HPImaging/sv-mbirct).*
+*Python code for fast MBIR (Model Based Iterative Reconstruction)  
+This is a python wrapper for High Performance Imaging's supervoxel C code, [HPImaging/sv-mbirct](https://github.com/HPImaging/sv-mbirct).*
 
 Full documentation is available at: https://svmbir.readthedocs.io
 

demo/demo_fanbeam.py

@@ -0,0 +1,50 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import svmbir
+
+"""
+Fan beam demo
+"""
+
+# Simulated sinogram parameters
+geometry = 'fan'
+dist_source_detector = 512.0
+magnification = 2.0
+num_views = 512
+num_channels = 512
+angles = np.linspace(-np.pi, np.pi, num_views, endpoint=False)
+
+# Reconstruction parameters
+img_size = 256
+sharpness = 0.0
+T = 1.0
+p = 1.2
+snr_db = 40.0
+
+# Generate phantom with a single slice
+phantom = svmbir.phantom.gen_shepp_logan(img_size,img_size)
+phantom = np.expand_dims(phantom, axis=0)
+sino = svmbir.project(phantom, angles, num_channels, geometry=geometry, dist_source_detector=dist_source_detector, magnification=magnification)
+
+# Perform MBIR reconstruction
+recon = svmbir.recon(sino, angles, num_rows=img_size, num_cols=img_size, T=T, p=p, sharpness=sharpness, snr_db=snr_db, geometry=geometry, dist_source_detector=dist_source_detector, magnification=magnification)
+
+# Compute Normalized Root Mean Squared Error
+nrmse = svmbir.phantom.nrmse(recon[0], phantom[0])
+
+#plt.ion()
+vmin = 1.0
+vmax = 1.2
+plt.figure(); plt.imshow(phantom[0],vmin=vmin,vmax=vmax,cmap='gray'); plt.colorbar()
+plt.title('Shepp Logan Phantom')
+
+plt.figure(); plt.imshow(np.squeeze(sino),cmap='gray'); plt.colorbar()
+plt.title('Sinogram')
+
+plt.figure(); plt.imshow(recon[0],vmin=vmin,vmax=vmax,cmap='gray'); plt.colorbar()
+plt.title(f'Reconstruction, nmrse={nrmse:.3f}')
+
+print("Close figures to continue")
+plt.show()
+

dev_scripts/test_pypi.sh

@@ -9,7 +9,7 @@
 # Set these accordingly:
 
 python_versions=("3.7" "3.8" "3.9")
-svmbir_version=0.2.7
+svmbir_version=0.2.8
 
 echo "*********************************************************"
 echo "**** Test install "
@@ -22,7 +22,7 @@ for pyv in ${python_versions[@]}; do
     conda activate sv${pyv}
 
     echo "*********** Installing svmbir ${svmbir_version} *************"
-    #pip install -i https://test.pypi.org/simple/ svmbir==$svmbir_version
+    #pip install -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple svmbir==$svmbir_version
     pip install svmbir==$svmbir_version
     pip install -r ../demo/requirements_demo.txt
 

docs/README.rst

@@ -1,14 +1,13 @@
-**svmbir** stands for Super-Voxel Model-Based Iterative Reconstruction.
-svmbir is an easy-to-use python package for fast parallel-beam reconstruction of tomography data using model-based priors.
+**svmbir** (super-voxel model-based iterative reconstruction) is an easy-to-use python package for fast iterative reconstruction of tomography data using model-based priors.
+The `svmbir <https://github.com/cabouman/svmbir>`_ package is a python interface to HPImaging's optimized 
+C implementation `[sv-mbirct] <https://github.com/HPImaging/sv-mbirct>`_
+of the super-voxel algorithm :cite:`wang2016high` :cite:`wang2017massively`.
 
 
 Features
 --------
-* Easy-to-use python code for fast parallel-beam MBIR (Model Based Iterative Reconstruction)
-
-* Interface to HPImaging's optimized C implementation `[sv-mbirct] <https://github.com/HPImaging/sv-mbirct>`_ of the super-voxel algorithm :cite:`wang2016high` :cite:`wang2017massively`
-
-* Supports MBIR reconstruction with Bayesian and Plug-and-Play prior models.
+* Supports MBIR reconstruction with Bayesian and Plug-and-Play prior models
+* Supports parallel-beam and fan-beam geometries
 
 
 System Requirements
@@ -23,11 +22,11 @@ Optional System Requirements
 Fastest reconstruction can be obtained with,
 
 * Intel-based CPU(s) supporting AVX2,AVX512
-* Intel ICC compiler (in "Parallel Studio XE", and now free "oneAPI")
+* Intel ICC compiler (included with Intel's free "oneAPI" toolkit)
 
 We also recommend:
 
-* Installation using conda environment
+* Installation to a conda virtual environment
 
 License
 -------

docs/requirements.txt

@@ -1,10 +1,10 @@
+numpy==1.21.*
+psutil>=5.8
+Pillow
+Cython
+ruamel.yaml
+pytest
 sphinx>=2.2
 sphinxcontrib-bibtex==1.0
-ipython>=6.3.1
-numpy~=1.20.1
-matplotlib~=3.3.4
-ruamel.yaml~=0.16.12
 sphinx-rtd-theme~=0.5.1
-psutil~=5.8.0
-Pillow~=8.4.0
-Cython~=0.29.21
+ipython>=6.3.1

docs/source/geom-fan.jpg

@@ -2,17 +2,36 @@
 Overview
 ========
 
-**svmbir** is a Python implementation of the Super-Voxel Model Based Iterative Reconstruction (MBIR) algorithm :cite:`wang2016high` :cite:`wang2017massively` for fast reconstruction of parallel beam 3D data.
+**svmbir** is a Python implementation of the Super-Voxel Model Based Iterative Reconstruction (MBIR) algorithm :cite:`wang2016high` :cite:`wang2017massively` for fast reconstruction of parallel-beam and fan-beam tomography data (3D).
 The code performs Bayesian reconstruction of tomographic data, so it is particularly well-suited for sparse view reconstruction from noisy data.
 It also has hooks to support Plug-and-Play prior models that can dramatically improve image quality :cite:`venkatakrishnan2013plug` :cite:`sreehari2016plug`.
+The reconstruction engine for *svmbir* is written and optimized in C 
+`[HPImaging/sv-mbirct] <https://github.com/HPImaging/sv-mbirct>`_ .
 
 **How does it work?**
 
-The super-voxel code can be 100x to 1000x faster than conventional code because it reorganizes operations in a way that is much better matched to a computer's cache structure. To do this, it precomputes the ``system matrix`` that describes the geometry of the tomography system and stores it in a file in encoded form. Whenever you do a reconstruction with a new geometry, the ``svmbir`` package automatically detects the new geometry, precomputes a new system matrix, and stores it in a library for future use. By default, the system matrices are stored in a subdirectory of your personal ``.cache`` directory in your home directory. So if they are taking too much space on your disk, you can remove these files at any time, and they will just be recomputed when needed.
+The super-voxel code can be 100x to 1000x faster than conventional MBIR code because it reorganizes operations in a way that is much better matched to a computer's cache structure. 
+Part of this involves precomputing a *system matrix* that models the system geometry, and encoding it in a layout that facilitates parallelization and reduces the required fetches from memory during reconstruction. 
+When system matrices are computed, they are stored to disk and will be automatically loaded whenever the same geometry is subsequently encountered. 
 
-**Troubleshooting**
+**Geometry**
+
+**svmbir** currently supports *parallel-beam* and *fan-beam* (equiangular, see below) imaging geometries.
+
+.. figure:: geom-fan.jpg
+   :width: 50%
+   :alt: fan beam geometry
+   :align: center
+
+   Equiangular fan-beam geometry
+
+.. figure:: geom-parallel.jpg
+   :width: 50%
+   :alt: parallel beam geometry
+   :align: center
+
+   Parallel-beam geometry
 
-*Cached System Matrix Problems:* Rare updates to the software package could include changes to the encoding of the system matrix and result in existing pre-computed matrix files being incompatible after the update. The last such update was on (2020-12-02). To remove the outdated files, delete the cache directory located at ``~/.cache/svmbir``. The package will regenerate the system matrices as needed at the time of reconstruction.
 
 *View Angle Ordering:* It is common practice to collect view data using techniques such as the "golden ratio" method in which the view angles are not collected in monotonically increasing order on the interval :math:`[0,2\pi)`. While ``svmbir`` will produce the correct reconstruction regardless of view ordering, its reconstruction speed will be substantially degraded when the views are not in monotone order. In this case, we highly recommend that users reorder the sinogram views using the provided ``sino_sort`` function. The ``sino_sort``  function first wraps the view angles modulo :math:`2\pi`, and then sorts the views to be in monotonically increasing order by view angle.
 
@@ -37,3 +56,11 @@ Using this convention, the 3D array, ``image``, will be in units of photons/AU.
 .. math::
 
     \mbox{image in photons/mm} = \frac{ \mbox{image in photons/ALU} }{ 5 \mbox{mm} / \mbox{ALU}}
+
+**Matrix caching**
+
+When system matrices are computed, they are stored to disk and will be automatically loaded whenever the same geometry is subsequently encountered. 
+By default, the system matrices are stored in the subfolder ``~/.cache/svmbir/sysmatrix`` of your home directory.
+The matrix files can be removed at any time, and should be periodically cleaned out to reduce disk use.
+Occasionally, updates to the software package include changes to the encoding of the system matrix, in which case the the cached matrix files should also be cleaned out to avoid incompatibility.
+

docs/source/geom-parallel.jpg

@@ -2,7 +2,7 @@
 Theory
 ======
 
-Super-Voxel Model-Based Iterative Reconstruction (SVMBIR) is a fast algorithm for computing MBIR reconstructions from parallel beam tomographic data.
+Super-Voxel Model-Based Iterative Reconstruction (SVMBIR) is a fast algorithm for computing MBIR reconstructions from tomographic data.
 MBIR reconstruction works by solving the following optimization problem
 
 .. math::
@@ -22,7 +22,7 @@ The forward model term has the form,
 
 where :math:`y` is the sinogram data, 
 where :math:`x` is the unknown image to be reconstructed, 
-:math:`A` is the linear projection operator for the parallel-beam geometry, 
+:math:`A` is the linear projection operator for the specified imaging geometry, 
 :math:`\Lambda` is the diagonal matrix of sinogram weights, :math:`\Vert y \Vert_\Lambda^2 = y^T \Lambda y`, and :math:`\sigma_y` is a parameter controling the assumed standard deviation of the measurement noise.
 
 These quantities correspond to the following python variables:

docs/source/overview.rst

@@ -1,4 +1,4 @@
-__version__ = '0.2.7'
+__version__ = '0.2.8'
 from .svmbir import *
 from .phantom import *
 __all__ = ['_svmbir_lib_path','_clear_cache','sino_sort','auto_sigma_x','auto_sigma_p','auto_sigma_y', 'calc_weights','project','backproject','recon']