Currently job artifacts in CI/CD pipelines on LRZ GitLab never expire. Starting from Wed 26.1.2022 the default expiration time will be 30 days (GitLab default). Currently existing artifacts in already completed jobs will not be affected by the change. The latest artifacts for all jobs in the latest successful pipelines will be kept. More information: https://gitlab.lrz.de/help/user/admin_area/settings/continuous_integration.html#default-artifacts-expiration

Commit bac3287d authored by Michael Loipführer's avatar Michael Loipführer Committed by Tobias Lasser
Browse files

improve solver doc, add python bindings guide

parent 6dfb8cc5
......@@ -124,3 +124,17 @@ with `pre-commit install` such that they are run before each commit.
None of the commit hooks will change anything in your commit, they mearly check and error if
something is wrong.
## Building the Documentation
The [elsa documentation](https://ciip.in.tum.de/elsadocs/) is automatically built and deployed through the CI for each commit to master.
To build it locally the following packages are required: `sphinx doxygen` which should be available in
most major linux distributions or via pip. Additionally, the following sphinx extensions need to be installed via pip:
`sphinx-rtd-theme m2r2 breathe`.
Then simply build the documentation using ninja
```
mkdir -p build
cd build
cmake .. -GNinja
ninja docs
```
The docs should then be available at `build/docs/sphinx`.
......@@ -27,9 +27,12 @@ author = 'Tobias Lasser'
# Add any Sphinx extension module names here, as strings. They can be
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
# ones.
extensions = [ "breathe", "m2r2", "sphinxcontrib.katex"]
extensions = [ "breathe", "m2r2", "sphinxcontrib.katex", "sphinx.ext.autosectionlabel"]
source_suffix = ['.rst', '.md']
# Make sure the target is unique
autosectionlabel_prefix_document = True
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
......
Python Bindings
-----------------------------
Elsa also comes with python bindings for almost all aspects of the framework.
This short guide aims to give an introduction into some simple cases and explain how one can
easily translate C++ code into python code for faster prototyping and experimenting.
One major benefit that comes with the python bindings is that we are able to natively
use numpy arrays with our elsa data containers making it easy to work with other tools such as
matplotlib.
### Setup the python bindings
In case you installed elsa via a `make install` the bindings should work out of the box.
If you do not want to install elsa to your system or are developing locally make sure to
add the path to your build folder to the `PYTHONPATH` for python to be able to find and import
the binding code.
Once everything is set up simply
```python
import pyelsa as elsa
```
### 2D example
To give a short outline into the python usage of elsa we will recreate the 2D example of the [Quickstart](Quickstart) section
in python.
```python
import pyelsa as elsa
import numpy as np
import matplotlib.pyplot as plt
size = np.array([128, 128])
phantom = elsa.PhantomGenerator.createModifiedSheppLogan(size)
volume_descriptor = phantom.getDataDescriptor()
# generate circular trajectory
num_angles = 180
arc = 360
sino_descriptor = elsa.CircleTrajectoryGenerator.createTrajectory(
num_angles, phantom.getDataDescriptor(), arc, size[0] * 100, size[0])
# setup operator for 2d X-ray transform
projector = elsa.SiddonsMethod(volume_descriptor, sino_descriptor)
# simulate the sinogram
sinogram = projector.apply(phantom)
# setup reconstruction problem
problem = elsa.WLSProblem(projector, sinogram)
# solve the problem
solver = elsa.CG(problem)
n_iterations = 20
reconstruction = solver.solve(n_iterations)
# plot the reconstruction
plt.imshow(np.array(reconstruction), '2D Reconstruction')
plt.show()
```
As you can see the code is essentially equivalent to the C++ code shown in the Quickstart guide.
All the top level elsa modules normally available through
```cpp
#include "elsa.h"
```
are available under the top level `pyelsa` module.
All C++ functions and classes essentially have the same signatures.
One important aspect of the python bindings is, however, that in places where the C++ code would expect
`Eigen` vectors or matrices we can natively use `numpy` arrays as well as convert elsa `DataContainer` back to numpy
via
```python
data_container: elsa.DataContainer = ...
back_to_numpy = np.array(data_container)
```
This is important if we e.g. want to show a reconstruction image using matplotlib as it does not support the elsa
data containers.
### 3D reconstruction viewing
The tight integration with numpy and matplotlib also enables us to directly implement a 3D reconstruction viewer
within our experimentation code.
Let's take a simple 3D phantom reconstruction example using a CUDA projector to be more performant
```python
import pyelsa as elsa
import numpy as np
import matplotlib.pyplot as plt
size = np.array([64, 64, 64]) # 3D now
phantom = elsa.PhantomGenerator.createModifiedSheppLogan(size)
volume_descriptor = phantom.getDataDescriptor()
# generate circular trajectory
num_angles = 180
arc = 360
sino_descriptor = elsa.CircleTrajectoryGenerator.createTrajectory(
num_angles, phantom.getDataDescriptor(), arc, size[0] * 100, size[0])
# setup operator for 2d X-ray transform
projector = elsa.JosephsMethodCUDA(volume_descriptor, sino_descriptor)
# simulate the sinogram
sinogram = projector.apply(phantom)
# setup reconstruction problem
problem = elsa.WLSProblem(projector, sinogram)
# solve the problem
solver = elsa.CG(problem)
n_iterations = 20
reconstruction = np.array(solver.solve(n_iterations))
```
We can now implement a simple matlab viewer plugin to scroll through our 3D reconstruction using the mouse wheel as shown in
the [matplotlib documentation](https://matplotlib.org/stable/gallery/event_handling/image_slices_viewer.html)
```python
class IndexTracker:
def __init__(self, ax, X):
self.ax = ax
ax.set_title('use scroll wheel to navigate images')
self.X = X
rows, cols, self.slices = X.shape
self.ind = self.slices // 2
self.im = ax.imshow(self.X[:, :, self.ind])
self.update()
def on_scroll(self, event):
if event.button == 'up':
self.ind = (self.ind + 1) % self.slices
else:
self.ind = (self.ind - 1) % self.slices
self.update()
def update(self):
self.im.set_data(self.X[:, :, self.ind])
self.ax.set_ylabel('slice %s' % self.ind)
self.im.axes.figure.canvas.draw()
```
we then simply set up our viewer
```python
fig, ax = plt.subplots(1, 1)
tracker = IndexTracker(ax, reconstruction)
fig.canvas.mpl_connect('scroll_event', tracker.on_scroll)
plt.show()
```
and scroll through our 3D reconstruction.
......@@ -12,6 +12,7 @@
guides/quickstart-cxx.md
guides/admm-cxx.md
guides/python_bindings.md
.. toctree::
:caption: Modules
......@@ -23,7 +24,7 @@
modules/operators
modules/functionals
modules/problems
modules/solvers
modules/solvers/solvers
modules/projectors
modules/projectors_cuda
modules/proximity_operators
......
************
elsa solvers
************
.. _elsa-solvers-api:
***********
elsa solvers API
***********
.. contents:: Table of Contents
......
.. _elsa-solvers-choosing-a-solver:
******************
Choosing a Solver
******************
.. contents:: Table of Contents
A first example
===============
This guide assumes you have read the :ref:`Quickstart Guide <guides/quickstart-cxx:Quickstart>` or are already familiar with elsa.
A typical reconstruction pipeline in elsa looks like follows.
We would have a volume descriptor for some given 2d or 3d volume defining the
are to be reconstructed as well as a detector descriptor that describes the trajectory
of our measurement positions. Finally we also have a sinogram that contains the actual
measurement data.
.. code-block:: cpp
VolumeDescriptor volumeDescriptor;
DataContainer sinogram;
PlanarDetectorDescriptor sinoDescriptor;
SiddonsMethod projector(volumeDescriptor, sinoDescriptor);
To compute a reconstruction we would construct an optimisation problem, typically
a weighted least squares problem ``WLSProblem`` which we would then solve using an iterative
optimisation method (e.g. Gradient Descent) to get the actual reconstructed image.
.. code-block:: cpp
WLSProblem wlsProblem(projector, sinogram);
GradientDescent solver(problem);
auto reconstruction = solver.solve(20);
Since there are multiple options for optimisation methods available in elsa choosing the
right solver for the task is quite important.
Solver Options
===============
The following table should be regarded as an overview over what solvers exists in elsa and what their rough properties are.
The stated convergence speeds should be taken with a grain of salt as they are evaluated on general problems
(see the later sections) and can differ when specific solvers are applied to more specialized problems.
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Solver | Requirements on the problem | Convergence |
+==========================================================+=======================================================+=============+
| Gradient Descent (GD) | None | Slow |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Nesterov's Fast Gradient Method (FGM) | None | Medium |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Optimized Gradient Method (OGM) | None | Medium |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Separable Quadratic Surrogate Ordered Subsets (SQS OS) | None | Fast |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Conjugate Gradient Descent (CG) | Problem must be convertible to a quadric problem | Fast |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Iterative Shrinkage-Thresholding Algorithm (ISTA) | Requires a regularization term | Medium |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) | Requires a regularization term | Medium |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
| Alternating Direction Method of Multipliers (ADMM) | <WIP> | <WIP> |
+----------------------------------------------------------+-------------------------------------------------------+-------------+
For a more detailed description on the internals of each solvers also see the :ref:`API documentation <elsa-solvers-api>`.
To learn more about how to use ordered subsets and specifically SQS OS read the :ref:`Ordered Subsets Introduction <elsa-solvers-ordered-subsets>`.
Solver Convergence
==================
This section will take a closer look at the convergence speed of the different solvers for a simple phantom reconstruction
problem. The ADMM solver will be excluded here as it requires some special tweaking and cannot be compared to the other solvers
as easily.
To compare the convergence speed of the elsa solvers we set up a simple 2D 128 x 128 phantom reconstruction problem.
We run each solver for a number of iterations and plot the Mean Square Error of the reconstruction compared to the ground truth
over the number of iteration. The code containing these experiments can be found
`here <https://gitlab.lrz.de/IP/elsa/-/blob/master/examples/solver_experiments.py>`_.
The SQS solver is plotted both in its ordered subset mode and in normal mode where it is roughly equivalent to the Nesterov's Fast Gradient Method.
.. image:: convergence_all_solvers.png
:width: 900
:align: center
:alt: SolverConvergence
As we can see the main winners of this experiment are Conjugate Gradient Descent (CG) and SQS in ordered subset mode (SQS OS).
The only other solver coming somewhat close in terms of convergence speed by number of iterations is the Optimized Gradient Method (OGM).
Since the computational complexity of the solvers differs quite a bit we also need to have a look at how fast they converge
when measured in terms of run time.
Time complexity
==================
To measure the actual time-based convergence speed of the different solvers we again run the same experiment using
a sample 2D 128 x 128 phantom reconstruction. We run each solver for a number of iterations up to 50 and plot the
Mean Square Error of the reconstruction over the actual execution time of the optimization process.
Note: the irregularities stem from running the timing experiment on a normal multi-threaded desktop pc
.. image:: time_all_solvers.png
:width: 900
:align: center
:alt: SolverTime
This experiment now paints a slightly different picture of the actual convergence speeds than the previous one.
Again the SQS OS and CG have the fastest convergence speed but the main surprise comes when looking at the FGM and OGM
solvers which have a very high convergence speed when measured in the number of required iterations. When only measuring
the execution time both of those solver have a slow initial convergence speed but do catch up to the other medium speed
solvers after some time. This is due to the fact that both of these solvers have a more expensive implementation than
other first order methods (e.g. Gradient Descent) while having better convergence properties. OGM even is proven to reach
the theoretical maximum convergence speed for first order solvers.
As a rule of thumb we can conclude that using either CG or SQS OS for general reconstruction problems in elsa will
yield the best results. In case of more specialized problems using a solver that is specifically designed for such
problems might be more desirable though.
.. _elsa-solvers-ordered-subsets:
******************
Ordered Subsets
******************
.. contents:: Table of Contents
Introduction
===============
The elsa problem and solvers API supports solvers using ordered subsets to achieve faster initial convergence speeds.
Note: the only solver which currently implements ordered subsets is the Separable Quadratic Surrogate (SQS) method.
First order solvers using ordered subset function similarly to stochastic gradient descent. Instead of operating on the
full input data the solver computes the gradients needed for its update step only on a subset (or batch in SGD terminology) of the input data.
This process of reducing the amount of data used in the gradient computation can greatly increase initial convergence speeds.
One essential aspect of using OS methods for tomographic reconstruction problems is the process of sampling the initial
input data into subsets as one requirement of OS methods is the approximate equivalence of the gradient of the full problem with the
gradient of the ordered subsets. Since tomographic measurement data has some inherent structure stemming from the geometry of the
detector and the detector poses during measurement the subset sampling process should yield subsets whose measurement geometry
approximates that of the full measurement data.
To use ordered subsets in elsa we would start with our normal reconstruction input of a 2D or 3D volume descriptor,
a data container with the measurement data and a detector descriptor containing the measurement geometry and trajectory.
.. code-block:: cpp
VolumeDescriptor volumeDescriptor;
DataContainer sinogram;
PlanarDetectorDescriptor sinoDescriptor;
Instead of now directly instantiating a ``WLSProblem`` with the descriptors and a projector we first need to divide
our input data into subsets. This is done by instantiating a ``SubsetSampler`` and pass it the number of subsets we want.
Note: we also need to pass the type of detector descriptor as a template parameter due to some implementation internals.
.. code-block:: cpp
index_t nSubsets{4};
SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
volumeDescriptor, sinoDescriptor, nSubsets);
After instantiating our subset sampler we can then create a ``WLSSubsetProblem`` which we can pass to a solver which
supports ordered subsets. The solvers will check whether they received a ``SubsetProblem`` and automatically
run in ordered subsets mode. The API of a ``SubsetProblem`` also supports its usage as a normal ``Problem`` but will
not yield any benefits from its subset processing.
.. code-block:: cpp
WLSSubsetProblem<real_t> problem(
*subsetSampler.getProjector<SiddonsMethod<real_t>>(),
subsetSampler.getPartitionedData(sinogram),
subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>());
SQS solver{problem};
auto reconstruction = solver.solve(20);
Sampling Strategies
===================
The ``SubsetSampler`` API provides several strategies to generate subsets from input data to accommodate for
different measurement trajectories and detector geometries in order to keep the subset gradients approximately equal.
The desired sampling strategy can simply be passed as an argument to the ``SubsetSampler`` constructor.
Currently there are two strategies implemented available as enums in ``SubsetSampler::SamplingStrategy``:
+---------------------------+--------------------------------------------------------------------+
| Strategy | Description |
+===========================+====================================================================+
| ``ROUND_ROBIN`` | | a simple round robin based splitting (the default method), may |
| | | not work well for more complex 3D trajectories |
+---------------------------+--------------------------------------------------------------------+
| ``ROTATIONAL_CLUSTERING`` | | tries to assign data points to subsets based on their rotational |
| | | distance such that each subset approximates the geometry of the |
| | | full measurement data |
+---------------------------+--------------------------------------------------------------------+
Convergence
===========
A comparison between the convergence of SQS OS and most other solvers in elsa can be found in the :ref:`Solvers Overview <elsa-solvers-choosing-a-solver>`.
One very tunable parameter when using ordered subset solvers is the number of subsets. Generally the more subsets used
during the optimization the faster the initial convergence speed will be. The limit to how many subsets can and should
be used for a given problem is given by the number of measurement data points.
To get a rough understanding of how the SQS OS method performs when tweaking the number of subsets we evaluated a small
3D 32 x 32 x 32 phantom reconstruction problem. The code for these experiments can be found
`here <https://gitlab.lrz.de/IP/elsa/-/blob/master/examples/solver_experiments.py>`_.
.. image:: convergence_sqs.png
:width: 900
:align: center
:alt: SQSConvergence
As we can see increasing the number of subsets generally improves the convergence speed in terms of the number of iterations
needed. The same holds true for the actual execution time needed to achieve a certain accuracy as the next plots shows.
Note: the irregularities stem from running the timing experiment on a normal multi-threaded desktop pc
.. image:: time_sqs.png
:width: 900
:align: center
:alt: SQSTime
Increasing the number of subsets used can greatly improve the actual convergence speed of an ordered subsets based solver.
In practice some experimentation might be needed to find the maximum number of subsets usable for the given problem and
input data set, as a too small number of data points per subset might lead to unstable convergence.
\ No newline at end of file
*******
elsa solvers
*******
.. toctree::
:hidden:
choosing_a_solver.rst
ordered_subsets.rst
api.rst
This is the documentation of elsa's solvers module. It provides implementations of different iterative solvers for
general or more specific optimisation problems.
Get started by reading a quick overview on :ref:`which solver to choose <elsa-solvers-choosing-a-solver>` or jump
directly to the :ref:`API documentation <elsa-solvers-api>`.
You can read more about the following topics:
#. :ref:`Choose a solver <elsa-solvers-choosing-a-solver>` gives a basic overview over the different solvers and their strengths.
#. :ref:`Ordered Subsets <elsa-solvers-ordered-subsets>` gives an overview over ordered subsets and supported solvers.
#. :ref:`Solvers API <elsa-solvers-api>` provides a full documentation on the elsa solvers API.
......@@ -27,6 +27,7 @@ if(ELSA_BUILD_PYTHON_BINDINGS)
${ELSA_MODULE_TARGET_NAME} bind_${ELSA_MODULE_NAME}.cpp
${PROJECT_SOURCE_DIR}/tools/bindings_generation/hints/${ELSA_MODULE_NAME}_hints.cpp ${MODULE_SOURCES}
CircleTrajectoryGenerator.h
SphereTrajectoryGenerator.h
)
endif()
......
......@@ -27,7 +27,6 @@ namespace elsa
* @param[in] A the full system matrix of the whole WSL problem
* @param[in] b data vector
* @param[in] subsetAs the system matrices corresponding to each subset
* @param[in] subsetBs a data vector with a Block Descriptor containing the data vector vor
* each subset
*/
WLSSubsetProblem(const LinearOperator<data_t>& A, const DataContainer<data_t>& b,
......
......@@ -34,7 +34,6 @@ namespace elsa
*
* @param[in] volumeDescriptor of the problem
* @param[in] detectorDescriptor describes the geometry and trajectory of the measurements
* @param[in] sinogram is the actual measurement data
* @param[in] nSubsets is number of subsets that should be generated
* @param[in] samplingStrategy the strategy with which to sample the subsets
*/
......@@ -54,7 +53,9 @@ namespace elsa
DataContainer<data_t> getPartitionedData(const DataContainer<data_t>& sinogram);
/**
* @brief return the full projector
* @brief return
*
* @tparam projector_t the type of projector to instantiate
*/
template <typename Projector_t>
std::unique_ptr<LinearOperator<data_t>> getProjector()
......@@ -64,6 +65,8 @@ namespace elsa
/**
* @brief return a list of projectors that correspond to each subset
*
* @tparam projector_t the type of projector to instantiate
*/
template <typename Projector_t>
std::vector<std::unique_ptr<LinearOperator<data_t>>> getSubsetProjectors()
......
......@@ -36,7 +36,6 @@ namespace elsa
* function.
*
* @param[in] problem the problem that is supposed to be solved
* @param[in] stepSize the fixed step size to be used while solving
*/
GradientDescent(const Problem<data_t>& problem);
......
......@@ -33,7 +33,7 @@ namespace elsa
* in ordered subset mode.
*
* @param[in] problem the problem that is supposed to be solved
* @param[in] momentumAcceleration whether or not to enable momentum acceleration
* @param[in] momentumAcceleration whether to enable Nesterov's momentum acceleration
* @param[in] epsilon affects the stopping condition
*/
SQS(const Problem<data_t>& problem, bool momentumAcceleration = true,
......
from dataclasses import dataclass, field
from typing import List
import time
import numpy as np
import matplotlib.pyplot as plt
import pyelsa as elsa
class IndexTracker:
def __init__(self, ax, X):
self.ax = ax
ax.set_title('use scroll wheel to navigate images')
self.X = X
rows, cols, self.slices = X.shape
self.ind = self.slices // 2
self.im = ax.imshow(self.X[:, :, self.ind])
self.update()
def on_scroll(self, event):
if event.button == 'up':
self.ind = (self.ind + 1) % self.slices
else:
self.ind = (self.ind - 1) % self.slices
self.update()
def update(self):
self.im.set_data(self.X[:, :, self.ind])
self.ax.set_ylabel('slice %s' % self.ind)
self.im.axes.figure.canvas.draw()
@dataclass
class SolverTest:
solver_class: elsa.Solver
solver_name: str
uses_subsets: bool = False
n_subsets: int = 5
needs_lasso: bool = False
sampling_strategy: str = 'ROUND_ROBIN' # or 'ROTATIONAL_CLUSTERING'
extra_args: dict = field(default_factory=dict)
def mse(optimized: np.ndarray, original: np.ndarray) -> float:
size = original.size
diff = (original - optimized) ** 2
return np.sum(diff) / size
def instantiate_solvers(solvers: List[SolverTest], do_3d=False):
# generate the phantom
size = np.array([128, 128]) if not do_3d else np.array([32, 32, 32])
phantom = elsa.PhantomGenerator.createModifiedSheppLogan(size)
volume_descriptor = phantom.getDataDescriptor()
num_poses = 180
if do_3d:
# generate spherical trajectory
n_circles = 5
sino_descriptor = elsa.SphereTrajectoryGenerator.createTrajectory(
num_poses, phantom.getDataDescriptor(), n_circles, elsa.SourceToCenterOfRotation(size[0] * 100),
elsa.CenterOfRotationToDetector(size[0]))
else:
# generate circular trajectory
arc = 360
sino_descriptor = elsa.CircleTrajectoryGenerator.createTrajectory(
num_poses, phantom.getDataDescriptor(), arc, elsa.SourceToCenterOfRotation(size[0] * 100),
elsa.CenterOfRotationToDetector(size[0]))
# setup operator for 2d X-ray transform
projector = elsa.SiddonsMethod(volume_descriptor, sino_descriptor)