simple_recon2d.cpp 1.95 KB
Newer Older
1
2
3
/// [simplerecon header begin]
#include "elsa.h"

4
5
#include <iostream>

6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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";
    }
}