PhantomGenerator.cpp 7.52 KB
Newer Older
1
2
#include "PhantomGenerator.h"
#include "EllipseGenerator.h"
3
#include "Logger.h"
4
#include "VolumeDescriptor.h"
5
#include "CartesianIndices.h"
6
7
8
9
10
11

#include <cmath>
#include <stdexcept>

namespace elsa
{
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
    template <typename data_t>
    DataContainer<data_t> PhantomGenerator<data_t>::createCirclePhantom(IndexVector_t volumesize,
                                                                        data_t radius)
    {
        VolumeDescriptor dd(volumesize);
        DataContainer<data_t> dc(dd);
        dc = 0;

        const Vector_t<data_t> sizef = volumesize.template cast<data_t>();
        const auto center = (sizef.array() / 2).matrix();

        for (auto pos : CartesianIndices(volumesize)) {
            const Vector_t<data_t> p = pos.template cast<data_t>();
            if ((p - center).norm() <= radius) {
                dc(pos) = 1;
            }
        }

        return dc;
    }

    template <typename data_t>
    DataContainer<data_t> PhantomGenerator<data_t>::createRectanglePhantom(IndexVector_t volumesize,
                                                                           IndexVector_t lower,
                                                                           IndexVector_t upper)
    {
        VolumeDescriptor dd(volumesize);
        DataContainer<data_t> dc(dd);
        dc = 0;

        for (auto pos : CartesianIndices(lower, upper)) {
            dc(pos) = 1;
        }

        return dc;
    }

49
50
51
52
53
    template <typename data_t>
    DataContainer<data_t> PhantomGenerator<data_t>::createModifiedSheppLogan(IndexVector_t sizes)
    {
        // sanity check
        if (sizes.size() < 2 || sizes.size() > 3)
54
            throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
55
                "PhantomGenerator::createModifiedSheppLogan: only 2d or 3d supported");
56
        if (sizes.size() == 2 && sizes[0] != sizes[1])
57
            throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
58
                "PhantomGenerator::createModifiedSheppLogan: 2d size has to be square");
59
        if (sizes.size() == 3 && (sizes[0] != sizes[1] || sizes[0] != sizes[2]))
60
            throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
61
                "PhantomGenerator::createModifiedSheppLogan: 3d size has to be cubed");
62

Jens Petit's avatar
Jens Petit committed
63
64
        Logger::get("PhantomGenerator")
            ->info("creating modified Shepp Logan phantom of size {}^{}", sizes[0], sizes.size());
65

66
        VolumeDescriptor dd(sizes);
67
        DataContainer<data_t> dc(dd);
68
        dc = 0;
69
70

        if (sizes.size() == 2) {
Jens Petit's avatar
Jens Petit committed
71
72
            EllipseGenerator<data_t>::drawFilledEllipse2d(dc, 1.0,
                                                          {scaleShift(dd, 0), scaleShift(dd, 0)},
73
                                                          {scale(dd, 0.69f), scale(dd, 0.92f)}, 0);
Jens Petit's avatar
Jens Petit committed
74
            EllipseGenerator<data_t>::drawFilledEllipse2d(
75
76
                dc, -0.8f, {scaleShift(dd, 0), scaleShift(dd, -0.0184f)},
                {scale(dd, 0.6624f), scale(dd, 0.8740f)}, 0);
Jens Petit's avatar
Jens Petit committed
77
            EllipseGenerator<data_t>::drawFilledEllipse2d(
78
79
                dc, -0.2f, {scaleShift(dd, 0.22f), scaleShift(dd, 0)},
                {scale(dd, 0.11f), scale(dd, 0.31f)}, -18);
Jens Petit's avatar
Jens Petit committed
80
            EllipseGenerator<data_t>::drawFilledEllipse2d(
81
82
                dc, -0.2f, {scaleShift(dd, -0.22f), scaleShift(dd, 0)},
                {scale(dd, 0.16f), scale(dd, 0.41f)}, 18);
Jens Petit's avatar
Jens Petit committed
83
            EllipseGenerator<data_t>::drawFilledEllipse2d(
84
85
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.35f)},
                {scale(dd, 0.21f), scale(dd, 0.25)}, 0);
Jens Petit's avatar
Jens Petit committed
86
            EllipseGenerator<data_t>::drawFilledEllipse2d(
87
88
89
90
91
92
93
94
95
96
97
98
99
100
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.1f)},
                {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
            EllipseGenerator<data_t>::drawFilledEllipse2d(
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.1f)},
                {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
            EllipseGenerator<data_t>::drawFilledEllipse2d(
                dc, 0.1f, {scaleShift(dd, -0.08f), scaleShift(dd, -0.605f)},
                {scale(dd, 0.046f), scale(dd, 0.023f)}, 0);
            EllipseGenerator<data_t>::drawFilledEllipse2d(
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.606f)},
                {scale(dd, 0.023f), scale(dd, 0.023f)}, 0);
            EllipseGenerator<data_t>::drawFilledEllipse2d(
                dc, 0.1f, {scaleShift(dd, 0.06f), scaleShift(dd, -0.605f)},
                {scale(dd, 0.023f), scale(dd, 0.046f)}, 0);
101
102
103
        }

        if (sizes.size() == 3) {
Jens Petit's avatar
Jens Petit committed
104
105
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
                dc, 1.0, {scaleShift(dd, 0), scaleShift(dd, 0), scaleShift(dd, 0)},
106
                {scale(dd, 0.69f), scale(dd, 0.92f), scale(dd, 0.81f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
107
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
108
109
                dc, -0.8f, {scaleShift(dd, 0), scaleShift(dd, -0.0184f), scaleShift(dd, 0)},
                {scale(dd, 0.6624f), scale(dd, 0.874f), scale(dd, 0.78f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
110
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
111
112
                dc, -0.2f, {scaleShift(dd, 0.22f), scaleShift(dd, 0), scaleShift(dd, 0)},
                {scale(dd, 0.11f), scale(dd, 0.31f), scale(dd, 0.22f)}, -18, 0, 10);
Jens Petit's avatar
Jens Petit committed
113
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
114
115
                dc, -0.2f, {scaleShift(dd, -0.22f), scaleShift(dd, 0), scaleShift(dd, 0)},
                {scale(dd, 0.16f), scale(dd, 0.41f), scale(dd, 0.28f)}, 18, 0, 10);
Jens Petit's avatar
Jens Petit committed
116
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
117
118
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.35f), scaleShift(dd, -0.15f)},
                {scale(dd, 0.21f), scale(dd, 0.25f), scale(dd, 0.41f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
119
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
120
121
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.1f), scaleShift(dd, 0.25f)},
                {scale(dd, 0.046f), scale(dd, 0.046f), scale(dd, 0.05f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
122
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
123
124
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.1f), scaleShift(dd, 0.25f)},
                {scale(dd, 0.046f), scale(dd, 0.046f), scale(dd, 0.05f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
125
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
126
127
                dc, 0.1f, {scaleShift(dd, -0.08f), scaleShift(dd, -0.605f), scaleShift(dd, 0)},
                {scale(dd, 0.046f), scale(dd, 0.023f), scale(dd, 0.05f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
128
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
129
130
                dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.606f), scaleShift(dd, 0)},
                {scale(dd, 0.023f), scale(dd, 0.023f), scale(dd, 0.02f)}, 0, 0, 0);
Jens Petit's avatar
Jens Petit committed
131
            EllipseGenerator<data_t>::drawFilledEllipsoid3d(
132
133
                dc, 0.1f, {scaleShift(dd, 0.06f), scaleShift(dd, -0.605f), scaleShift(dd, 0)},
                {scale(dd, 0.023f), scale(dd, 0.046f), scale(dd, 0.02f)}, 0, 0, 0);
134
135
136
137
138
139
140
141
        }

        return dc;
    }

    template <typename data_t>
    index_t PhantomGenerator<data_t>::scale(const DataDescriptor& dd, data_t value)
    {
142
143
        return std::lround(
            value * static_cast<data_t>(dd.getNumberOfCoefficientsPerDimension()[0] - 1) / 2.0f);
144
145
146
147
148
    }

    template <typename data_t>
    index_t PhantomGenerator<data_t>::scaleShift(const DataDescriptor& dd, data_t value)
    {
149
150
151
        return std::lround(value
                           * static_cast<data_t>(dd.getNumberOfCoefficientsPerDimension()[0] - 1)
                           / 2.0f)
Jens Petit's avatar
Jens Petit committed
152
               + (dd.getNumberOfCoefficientsPerDimension()[0] / 2);
153
154
155
156
157
158
159
    }

    // ------------------------------------------
    // explicit template instantiation
    template class PhantomGenerator<float>;
    template class PhantomGenerator<double>;

160
} // namespace elsa