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