Commit d7555c08 authored by David Frank's avatar David Frank
Browse files

Rework quickstart guide to refer to code snippets from example

parent 83523fef
Pipeline #834975 passed with stages
in 9 minutes and 51 seconds
Quickstart
----------
##########
To demonstrate our framework, we're working through a simple exampled for Computed Tomography.
We'll be using the C++ interface of our framework. Please note, that this is a brief peak at the
code usage, not an introduction to tomographic reconstruction.
### 2D example
2D example
**********
To make the examples short and less noisy, we first include all the headers of elsa and
pull the elsa namespace into scope.
```cpp
#include "elsa.h"
using namespace elsa;
```
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon header begin
:end-before: simplerecon header end
Next we'll create the template to put all the code inside:
.. code-block:: cpp
// The code from above
#include "elsa.h"
using namespace elsa
void example2d() {
// From now all the code goes here
}
int main() {
example2d();
}
The first thing, we have to do is set up some phantom data. Working with phantom data is quite
useful for many applications, as we have a ground truth and compute error norms.
```cpp
const auto size = IndexVector_t::Constant(2, 1, 256);
const auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
```
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon phantom create
:end-before: simplerecon phantom create
Here we settle on a 2D Shepp-Logan phantom. To have a look at the phantom, you can write it to
disk using the `PGM <https://de.wikipedia.org/wiki/Portable_Anymap>`_ format, like this:
Here we settle on a 2D Shepp-Logan phantom of size 256x256. To have a look at the phantom,
you can write it to disk using the [PGM](https://de.wikipedia.org/wiki/Portable_Anymap) format,
like this:
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon phantom write
:end-before: simplerecon phantom write
```cpp
PGM::write(phantom, "phantom.pgm");
```
But if you change the file extension to `edf`, you can write it - well - to the `edf` file format.
`EDF <https://en.Wikipedia.org/wiki/European_Data_Format>`_ files can be read back in to elsa, PGM
files can not. PCM files are just meant to as a quick and easy way to debug reconstructions. If you
plan to further use the output, you should use EDF files.
From here, you can write all :code:`DataContainer`'s, like the sinogram or the reconstruction, we will
not explicitly add it here in the code, but you are encouraged to write the results to disc hand
have a look.
Next, we have to define a trajectory. The trajectory defines the position of the source and the
detector. In our case, this is a fan-beam geometry setup.
```cpp
const index_t numAngles = 180;
const index_t arc = 360;
const real_t distance = size(0);
const auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
numAngles, phantom.getDataDescriptor(), arc, distance * 100.0f, distance);
```
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon trajectory
:end-before: simplerecon trajectory
The trajectory has 180 positions (`numAngles`) over the whole circle (`arc`). The two distances are
the distance from the source to the center of the volume and the distance from the center of the volume
to the detector.
to the detector. In this case we make the distance from the center to detector and source dependent
on the size of the volume. The values provided here are okay defaults.
It's usually a good idea to move the source further away from the volume, and the detector closer to it.
This way also the edges of the phantom are covered nicely, and no artifacts appear during
reconstruction (but don't trust me, try it!).
Now, we need to create the sinogram. The sinogram is the result of projecting the phantom.
```cpp
const auto volumeDescriptor = phantom.getDataDescriptor();
SiddonsMethod projector(volumeDescriptor, *sinoDescriptor);
auto sinogram = projector.apply(phantom);
// Also write it to disk and have a look at it!
PGM::write(phantom, "sinogram.pgm");
```
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon sinogram
:end-before: simplerecon sinogram
The `projector` is the approximation of the Radon transform. Our current projectors are implemented
based on ray traversals.
With this setup, we finally have everything to setup the tomographic problem.
```cpp
const WLSProblem wlsProblem(projector, sinogram);
const CG solver(wlsProblem);
```
Without getting to into the math behind it (checkout the [paper](#citation), if you want to dive
into it), this sets up an optimization problem. We want to find an solution to this problem
In this case, we'll find it it using the
[Conjugate gradient method](https://en.wikipedia.org/wiki/Conjugate_gradient_method), or CG for short.
It's an iterative algorithm to solve. As this is still a quite simple example, we don't need
to let it run for too many iterations.
```cpp
const index_t noIterations{20};
const auto reconstruction = solver.solve(noIterations);
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon solver
:end-before: simplerecon solver
PGM::write(reconstruction, "reconstruction.pgm");
```
Without getting to into the math behind it (checkout the
`paper <https://doi.org/10.1117/12.2534833>`_ , if you want to dive into it), this sets up an
optimization problem. We want to find an solution to this problem In this case, we'll find it it
using the `Conjugate gradient method <https://en.wikipedia.org/wiki/Conjugate_gradient_method>`_,
or CG for short. It's an iterative algorithm to solve. As this is still a quite simple example, we
don't need to let it run for too many iterations.
Now, you have the reconstruction! In the best case, this should already look quite similar to the
original phantom. We can have a look at the difference and the L2-Norm:
```cpp
PGM::write(phantom - reconstruction, "reconstruction.pgm");
infoln("L2-Norm of the phantom: {:f.5}", phantom.l2Norm());
infoln("L2-Norm of the reconstruction: {:f.5}", reconstruction.l2Norm());
```
### 3D example
If you want to run a 3D reconstruction, all you need to change, is the initial size at the beginning.
```cpp
const auto size = IndexVector_t::Constant(3, 1, 64);
const auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
```
In this case, we create a 64x64x64 sized phantom, as it's already 4 times larger than the 2D 256x256
phantom, and therefore, will also take longer to compute.
With the current capabilities of the framework, you'll also have to remove or comment out the `PGM::write`
calls. We don't support 3D phantoms to be printed with it. If you still want to look at it, you can
replace the calls with `EDF::write`, which is better suited.
### CUDA projectors
To speed up the reconstruction up, we'll leverage CUDA in the next step. Be sure that you've set it up
and have a CUDA capable card.
elsa by default checks if you've got a working CUDA compiler, therefore, you should be good to go.
If something went wrong, pass CMake the flag: `ELSA_BUILD_CUDA_PROJECTORS=ON`. This should enable
the CUDA projectors.
In the above example, the `SiddondsMethod` projector was used. We'll replace it with the
`JosephsMethodCUDA`:
```cpp
JosephsMethodCUDA projector(volumeDescriptor, *sinoDescriptor);
```
.. literalinclude:: ../../examples/simple_recon2d.cpp
:language: cpp
:start-after: simplerecon analysis
:end-before: simplerecon analysis
The reconstruction should now be quite a bit faster.
......@@ -10,7 +10,7 @@
:caption: Guides
:maxdepth: 2
guides/quickstart-cxx.md
guides/quickstart-cxx.rst
guides/python_bindings.md
modules/solvers/choosing_a_solver.rst
guides/admm-cxx.md
......
# Create binaries for examples in bin/examples
# Create binaries for examples in bin/examples
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin/examples)
function(add_example example_name example_source)
add_executable(${example_name} ${example_source})
target_link_libraries(${example_name} PRIVATE elsa::all)
target_compile_features(${example_name} PUBLIC cxx_std_17)
set_quickvec_options(${example_name} ${example_source})
add_dependencies(build-examples ${example_name})
endfunction()
# define a target for all examples
add_custom_target(build-examples)
# build a simple 2D reconstruction
add_example(simple_recon2d simple_recon2d.cpp)
# build the 2d example program
add_executable(example2d example2d.cpp)
target_link_libraries(example2d PRIVATE elsa::all)
target_compile_features(example2d PUBLIC cxx_std_17)
set_quickvec_options(example2d example2d.cpp)
add_dependencies(build-examples example2d)
add_example(example2d example2d.cpp)
# build the 2d limited angle trajectory example program
add_executable(limited_angle_example2d limited_angle_example2d.cpp limited_angle_example2d.cpp)
target_link_libraries(limited_angle_example2d PRIVATE elsa::all)
target_compile_features(limited_angle_example2d PUBLIC cxx_std_17)
set_quickvec_options(limited_angle_example2d limited_angle_example2d.cpp limited_angle_example2d.cpp)
add_dependencies(build-examples limited_angle_example2d)
add_example(limited_angle_example2d limited_angle_example2d.cpp)
# build the 2d example admm program
add_executable(example2d_admm example2d_admm.cpp)
target_link_libraries(example2d_admm PRIVATE elsa::all)
target_compile_features(example2d_admm PUBLIC cxx_std_17)
set_quickvec_options(example2d_admm example2d_admm.cpp)
add_dependencies(build-examples example2d_admm)
add_example(example2d_admm example2d_admm.cpp)
# build the shearlet example program
add_executable(shearlet_example example2d_shearlet.cpp)
target_link_libraries(shearlet_example PRIVATE elsa::all)
target_compile_features(shearlet_example PUBLIC cxx_std_17)
set_quickvec_options(shearlet_example example2d_shearlet.cpp)
add_dependencies(build-examples shearlet_example)
# add_executable(exampleMl exampleMl.cpp)
# target_link_libraries(exampleMl elsa::all)
# target_compile_features(exampleMl PUBLIC cxx_std_17)
# add_dependencies(build-examples exampleMl)
# add_executable(example_MNIST example_MNIST.cpp)
# target_link_libraries(example_MNIST elsa::all)
# target_compile_features(example_MNIST PUBLIC cxx_std_17)
# add_dependencies(build-examples example_MNIST)
add_example(shearlet_example example2d_shearlet.cpp)
if(ELSA_BUILD_CUDA_PROJECTORS)
include(CheckLanguage)
......@@ -49,17 +34,9 @@ if(ELSA_BUILD_CUDA_PROJECTORS)
enable_language(CUDA)
# build the 3d example program
add_executable(example3d example3d.cpp)
target_link_libraries(example3d PRIVATE elsa::all)
target_compile_features(example3d PUBLIC cxx_std_17)
set_quickvec_options(example3d example3d.cpp)
add_dependencies(build-examples example3d)
add_example(example3d example3d.cpp)
# build the GPU projector speed test program
add_executable(speed_test speed_test.cpp)
target_link_libraries(speed_test PRIVATE elsa::all)
target_compile_features(speed_test PUBLIC cxx_std_17)
set_quickvec_options(speed_test speed_test.cpp)
add_dependencies(build-examples speed_test)
add_example(speed_test speed_test.cpp)
endif()
endif()
/// [simplerecon header begin]
#include "elsa.h"
using namespace elsa;
/// [simplerecon header end]
void example2d()
{
/// [simplerecon phantom create]
IndexVector_t size({{128, 128}});
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
/// [simplerecon phantom create]
/// [simplerecon phantom write]
io::write(phantom, "2dphantom.pgm");
/// [simplerecon phantom write]
/// [simplerecon trajectory]
const index_t numAngles = 180;
const index_t arc = 360;
const auto distance = static_cast<real_t>(size[0]);
const auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
numAngles, phantom.getDataDescriptor(), arc, distance * 100.0f, distance);
/// [simplerecon trajectory]
/// [simplerecon sinogram]
const auto& phanDescriptor = dynamic_cast<const VolumeDescriptor&>(phantom.getDataDescriptor());
SiddonsMethod projector(phanDescriptor, *sinoDescriptor);
// simulate the sinogram
auto sinogram = projector.apply(phantom);
/// [simplerecon sinogram]
// write the sinogram out
io::write(sinogram, "2dsinogram.pgm");
// setup reconstruction problem
/// [simplerecon solver]
WLSProblem wlsProblem(projector, sinogram);
CG solver(wlsProblem);
index_t iterations{20};
auto reconstruction = solver.solve(iterations);
/// [simplerecon solver]
// write the reconstruction out
io::write(reconstruction, "2dreconstruction.pgm");
/// [simplerecon analysis]
DataContainer diff = phantom - reconstruction;
io::write(diff, "2ddiff.pgm");
Logger::get("example")->info("L2-Norm of the phantom: {:.4f}", phantom.l2Norm());
Logger::get("example")->info("L2-Norm of the phantom: {:.4f}", reconstruction.l2Norm());
/// [simplerecon analysis]
}
int main()
{
try {
example2d();
} catch (std::exception& e) {
std::cerr << "An exception occurred: " << e.what() << "\n";
}
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment