Line data Source code
1 : #include "PhantomGenerator.h"
2 : #include "EllipseGenerator.h"
3 : #include "Logger.h"
4 : #include "VolumeDescriptor.h"
5 : #include "CartesianIndices.h"
6 :
7 : #include <cmath>
8 : #include <stdexcept>
9 :
10 : namespace elsa
11 : {
12 : template <typename data_t>
13 : DataContainer<data_t> PhantomGenerator<data_t>::createCirclePhantom(IndexVector_t volumesize,
14 : data_t radius)
15 32 : {
16 32 : VolumeDescriptor dd(volumesize);
17 32 : DataContainer<data_t> dc(dd);
18 32 : dc = 0;
19 :
20 32 : const Vector_t<data_t> sizef = volumesize.template cast<data_t>();
21 32 : const auto center = (sizef.array() / 2).matrix();
22 :
23 69632 : for (auto pos : CartesianIndices(volumesize)) {
24 69632 : const Vector_t<data_t> p = pos.template cast<data_t>();
25 69632 : if ((p - center).norm() <= radius) {
26 12038 : dc(pos) = 1;
27 12038 : }
28 69632 : }
29 :
30 32 : return dc;
31 32 : }
32 :
33 : template <typename data_t>
34 : DataContainer<data_t> PhantomGenerator<data_t>::createRectanglePhantom(IndexVector_t volumesize,
35 : IndexVector_t lower,
36 : IndexVector_t upper)
37 8 : {
38 8 : VolumeDescriptor dd(volumesize);
39 8 : DataContainer<data_t> dc(dd);
40 8 : dc = 0;
41 :
42 9856 : for (auto pos : CartesianIndices(lower, upper)) {
43 9856 : dc(pos) = 1;
44 9856 : }
45 :
46 8 : return dc;
47 8 : }
48 :
49 : template <typename data_t>
50 : DataContainer<data_t> PhantomGenerator<data_t>::createModifiedSheppLogan(IndexVector_t sizes)
51 23 : {
52 : // sanity check
53 23 : if (sizes.size() < 2 || sizes.size() > 3)
54 0 : throw InvalidArgumentError(
55 0 : "PhantomGenerator::createModifiedSheppLogan: only 2d or 3d supported");
56 23 : if (sizes.size() == 2 && sizes[0] != sizes[1])
57 0 : throw InvalidArgumentError(
58 0 : "PhantomGenerator::createModifiedSheppLogan: 2d size has to be square");
59 23 : if (sizes.size() == 3 && (sizes[0] != sizes[1] || sizes[0] != sizes[2]))
60 0 : throw InvalidArgumentError(
61 0 : "PhantomGenerator::createModifiedSheppLogan: 3d size has to be cubed");
62 :
63 23 : Logger::get("PhantomGenerator")
64 23 : ->info("creating modified Shepp Logan phantom of size {}^{}", sizes[0], sizes.size());
65 :
66 23 : VolumeDescriptor dd(sizes);
67 23 : DataContainer<data_t> dc(dd);
68 23 : dc = 0;
69 :
70 23 : if (sizes.size() == 2) {
71 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(dc, 1.0,
72 16 : {scaleShift(dd, 0), scaleShift(dd, 0)},
73 16 : {scale(dd, 0.69f), scale(dd, 0.92f)}, 0);
74 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
75 16 : dc, -0.8f, {scaleShift(dd, 0), scaleShift(dd, -0.0184f)},
76 16 : {scale(dd, 0.6624f), scale(dd, 0.8740f)}, 0);
77 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
78 16 : dc, -0.2f, {scaleShift(dd, 0.22f), scaleShift(dd, 0)},
79 16 : {scale(dd, 0.11f), scale(dd, 0.31f)}, -18);
80 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
81 16 : dc, -0.2f, {scaleShift(dd, -0.22f), scaleShift(dd, 0)},
82 16 : {scale(dd, 0.16f), scale(dd, 0.41f)}, 18);
83 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
84 16 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.35f)},
85 16 : {scale(dd, 0.21f), scale(dd, 0.25)}, 0);
86 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
87 16 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.1f)},
88 16 : {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
89 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
90 16 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.1f)},
91 16 : {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
92 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
93 16 : dc, 0.1f, {scaleShift(dd, -0.08f), scaleShift(dd, -0.605f)},
94 16 : {scale(dd, 0.046f), scale(dd, 0.023f)}, 0);
95 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
96 16 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.606f)},
97 16 : {scale(dd, 0.023f), scale(dd, 0.023f)}, 0);
98 16 : EllipseGenerator<data_t>::drawFilledEllipse2d(
99 16 : dc, 0.1f, {scaleShift(dd, 0.06f), scaleShift(dd, -0.605f)},
100 16 : {scale(dd, 0.023f), scale(dd, 0.046f)}, 0);
101 16 : }
102 :
103 23 : if (sizes.size() == 3) {
104 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
105 7 : dc, 1.0, {scaleShift(dd, 0), scaleShift(dd, 0), scaleShift(dd, 0)},
106 7 : {scale(dd, 0.69f), scale(dd, 0.92f), scale(dd, 0.81f)}, 0, 0, 0);
107 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
108 7 : dc, -0.8f, {scaleShift(dd, 0), scaleShift(dd, -0.0184f), scaleShift(dd, 0)},
109 7 : {scale(dd, 0.6624f), scale(dd, 0.874f), scale(dd, 0.78f)}, 0, 0, 0);
110 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
111 7 : dc, -0.2f, {scaleShift(dd, 0.22f), scaleShift(dd, 0), scaleShift(dd, 0)},
112 7 : {scale(dd, 0.11f), scale(dd, 0.31f), scale(dd, 0.22f)}, -18, 0, 10);
113 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
114 7 : dc, -0.2f, {scaleShift(dd, -0.22f), scaleShift(dd, 0), scaleShift(dd, 0)},
115 7 : {scale(dd, 0.16f), scale(dd, 0.41f), scale(dd, 0.28f)}, 18, 0, 10);
116 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
117 7 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.35f), scaleShift(dd, -0.15f)},
118 7 : {scale(dd, 0.21f), scale(dd, 0.25f), scale(dd, 0.41f)}, 0, 0, 0);
119 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
120 7 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.1f), scaleShift(dd, 0.25f)},
121 7 : {scale(dd, 0.046f), scale(dd, 0.046f), scale(dd, 0.05f)}, 0, 0, 0);
122 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
123 7 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.1f), scaleShift(dd, 0.25f)},
124 7 : {scale(dd, 0.046f), scale(dd, 0.046f), scale(dd, 0.05f)}, 0, 0, 0);
125 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
126 7 : dc, 0.1f, {scaleShift(dd, -0.08f), scaleShift(dd, -0.605f), scaleShift(dd, 0)},
127 7 : {scale(dd, 0.046f), scale(dd, 0.023f), scale(dd, 0.05f)}, 0, 0, 0);
128 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
129 7 : dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.606f), scaleShift(dd, 0)},
130 7 : {scale(dd, 0.023f), scale(dd, 0.023f), scale(dd, 0.02f)}, 0, 0, 0);
131 7 : EllipseGenerator<data_t>::drawFilledEllipsoid3d(
132 7 : dc, 0.1f, {scaleShift(dd, 0.06f), scaleShift(dd, -0.605f), scaleShift(dd, 0)},
133 7 : {scale(dd, 0.023f), scale(dd, 0.046f), scale(dd, 0.02f)}, 0, 0, 0);
134 7 : }
135 :
136 23 : return dc;
137 23 : }
138 :
139 : template <typename data_t>
140 : index_t PhantomGenerator<data_t>::scale(const DataDescriptor& dd, data_t value)
141 530 : {
142 530 : return std::lround(
143 530 : value * static_cast<data_t>(dd.getNumberOfCoefficientsPerDimension()[0] - 1) / 2.0f);
144 530 : }
145 :
146 : template <typename data_t>
147 : index_t PhantomGenerator<data_t>::scaleShift(const DataDescriptor& dd, data_t value)
148 530 : {
149 530 : return std::lround(value
150 530 : * static_cast<data_t>(dd.getNumberOfCoefficientsPerDimension()[0] - 1)
151 530 : / 2.0f)
152 530 : + (dd.getNumberOfCoefficientsPerDimension()[0] / 2);
153 530 : }
154 :
155 : // ------------------------------------------
156 : // explicit template instantiation
157 : template class PhantomGenerator<float>;
158 : template class PhantomGenerator<double>;
159 :
160 : } // namespace elsa
|
