Line data Source code
1 : #include "Ellipsoid.h"
2 :
3 : namespace elsa::phantoms
4 : {
5 :
6 : template <typename data_t>
7 : Ellipsoid<data_t>::Ellipsoid(data_t amplit, elsa::phantoms::Vec3i center,
8 : Vec3X<data_t> halfAxis, Vec3X<data_t> eulers)
9 : : _amplit{amplit}, _center{center}, _halfAxis{halfAxis}, _eulers{eulers}
10 99 : {
11 99 : bSqrcSqr =
12 99 : _halfAxis[INDEX_B] * _halfAxis[INDEX_B] * _halfAxis[INDEX_C] * _halfAxis[INDEX_C];
13 99 : aSqrcSqr =
14 99 : _halfAxis[INDEX_A] * _halfAxis[INDEX_A] * _halfAxis[INDEX_C] * _halfAxis[INDEX_C];
15 99 : aSqrbSqr =
16 99 : _halfAxis[INDEX_A] * _halfAxis[INDEX_A] * _halfAxis[INDEX_B] * _halfAxis[INDEX_B];
17 99 : aSqrbSqrcSqr = aSqrbSqr * _halfAxis[INDEX_C] * _halfAxis[INDEX_C];
18 :
19 99 : rotated = (std::abs(eulers[INDEX_PHI]) + std::abs(eulers[INDEX_THETA])
20 99 : + std::abs(eulers[INDEX_PSI]))
21 99 : > 0;
22 :
23 99 : if (rotated) {
24 18 : fillRotationMatrix(eulers, rot);
25 18 : rot.transposeInPlace();
26 18 : }
27 99 : };
28 :
29 : template <typename data_t>
30 : bool Ellipsoid<data_t>::isInEllipsoid(const Vec3i& idx) const
31 51376 : {
32 51376 : return (double(idx[INDEX_X] * idx[INDEX_X]) * bSqrcSqr
33 51376 : + double(idx[INDEX_Y] * idx[INDEX_Y]) * aSqrcSqr
34 51376 : + double(idx[INDEX_Z] * idx[INDEX_Z]) * aSqrbSqr)
35 51376 : <= aSqrbSqrcSqr;
36 51376 : };
37 :
38 : template <typename data_t>
39 : bool Ellipsoid<data_t>::isInEllipsoid(const Vec3X<data_t>& idx) const
40 99422 : {
41 :
42 99422 : return (idx[INDEX_X] * idx[INDEX_X] * bSqrcSqr + idx[INDEX_Y] * idx[INDEX_Y] * aSqrcSqr
43 99422 : + idx[INDEX_Z] * idx[INDEX_Z] * aSqrbSqr)
44 99422 : <= aSqrbSqrcSqr;
45 99422 : };
46 :
47 : template <Blending b, typename data_t>
48 : void rasterizeNoRotation(Ellipsoid<data_t>& el, VolumeDescriptor& dd, DataContainer<data_t>& dc)
49 79 : {
50 :
51 79 : Vec3i idx(3);
52 79 : idx << 0, 0, 0;
53 :
54 79 : auto strides = dd.getProductOfCoefficientsPerDimension();
55 :
56 : // global offests for a fast memory reuse in the for loops
57 79 : index_t zOffset = 0;
58 79 : index_t yOffset = 0;
59 :
60 79 : index_t zOffsetNeg = 0;
61 79 : index_t yOffsetNeg = 0;
62 :
63 79 : index_t xOffset = 0;
64 79 : index_t xOffsetNeg = 0;
65 :
66 79 : Vec3i _center = el.getCenter();
67 79 : data_t _amplit = el.getAmplitude();
68 :
69 : // Slicing on z axis to lower cache misses
70 423 : for (; el.isInEllipsoid(idx); idx[INDEX_Z]++) {
71 344 : zOffset = (idx[INDEX_Z] + _center[INDEX_Z]) * strides[INDEX_Z];
72 344 : zOffsetNeg = (-idx[INDEX_Z] + _center[INDEX_Z]) * strides[INDEX_Z];
73 :
74 3768 : for (; el.isInEllipsoid(idx); idx[INDEX_Y]++) {
75 3424 : yOffset = (idx[INDEX_Y] + _center[INDEX_Y]) * strides[INDEX_Y];
76 3424 : yOffsetNeg = (-idx[INDEX_Y] + _center[INDEX_Y]) * strides[INDEX_Y];
77 :
78 47185 : for (; el.isInEllipsoid(idx); idx[INDEX_X]++) {
79 43761 : xOffset = (idx[INDEX_X] + _center[INDEX_X]) * strides[INDEX_X];
80 43761 : xOffsetNeg = (-idx[INDEX_X] + _center[INDEX_X]) * strides[INDEX_X];
81 :
82 43761 : blend<b>(dc, xOffset + yOffset + zOffset, _amplit);
83 :
84 : // Voxel in ellipsoids can be mirrored 8 times if they are not on a mirror
85 : // plane. Depending on the voxels location they are mirrored or not.
86 : // This exclude prevents double increment of the same voxel on the mirror plane.
87 :
88 43761 : if (idx[INDEX_X] != 0) {
89 40337 : blend<b>(dc, xOffsetNeg + yOffset + zOffset, _amplit);
90 40337 : }
91 43761 : if (idx[INDEX_Y] != 0) {
92 40931 : blend<b>(dc, xOffset + yOffsetNeg + zOffset, _amplit);
93 40931 : }
94 43761 : if (idx[INDEX_Z] != 0) {
95 41100 : blend<b>(dc, xOffset + yOffset + zOffsetNeg, _amplit);
96 41100 : }
97 :
98 43761 : if (idx[INDEX_X] != 0 && idx[INDEX_Y] != 0) {
99 37851 : blend<b>(dc, xOffsetNeg + yOffsetNeg + zOffset, _amplit);
100 37851 : }
101 :
102 43761 : if (idx[INDEX_X] != 0 && idx[INDEX_Z] != 0) {
103 37990 : blend<b>(dc, xOffsetNeg + yOffset + zOffsetNeg, _amplit);
104 37990 : }
105 :
106 43761 : if (idx[INDEX_Y] != 0 && idx[INDEX_Z] != 0) {
107 38542 : blend<b>(dc, xOffset + yOffsetNeg + zOffsetNeg, _amplit);
108 38542 : }
109 :
110 43761 : if (idx[INDEX_X] != 0 && idx[INDEX_Y] != 0 && idx[INDEX_Z] != 0) {
111 35697 : blend<b>(dc, xOffsetNeg + yOffsetNeg + zOffsetNeg, _amplit);
112 35697 : }
113 43761 : };
114 3424 : idx[INDEX_X] = 0;
115 3424 : }
116 344 : idx[INDEX_Y] = 0;
117 344 : }
118 79 : };
119 :
120 : template <typename data_t>
121 : data_t Ellipsoid<data_t>::getRoundMaxHalfWidth() const
122 20 : {
123 20 : data_t max = _halfAxis.colwise().maxCoeff()[0];
124 20 : return std::ceil(max);
125 20 : };
126 :
127 : /**
128 : * Bounding Box in object space
129 : */
130 : template <typename data_t>
131 : std::array<data_t, 6>
132 : boundingBox(const data_t maxHalfAxis, const Vec3i& _center,
133 : const IndexVector_t& ncpd /*numberOfCoefficientsPerDimension zero based*/)
134 20 : {
135 :
136 20 : data_t minX, minY, minZ, maxX, maxY, maxZ;
137 :
138 20 : minX = std::max(-maxHalfAxis, data_t(-_center[INDEX_X]));
139 20 : minY = std::max(-maxHalfAxis, data_t(-_center[INDEX_Y]));
140 20 : minZ = std::max(-maxHalfAxis, data_t(-_center[INDEX_Z]));
141 :
142 20 : maxX = std::min(data_t(ncpd[INDEX_X] - 1 - _center[INDEX_X]), maxHalfAxis);
143 20 : maxY = std::min(data_t(ncpd[INDEX_Y] - 1 - _center[INDEX_Y]), maxHalfAxis);
144 20 : maxZ = std::min(data_t(ncpd[INDEX_Z] - 1 - _center[INDEX_Z]), maxHalfAxis);
145 :
146 20 : return {minX, minY, minZ, maxX, maxY, maxZ};
147 20 : };
148 :
149 : template <Blending b, typename data_t>
150 : void rasterizeWithClipping(Ellipsoid<data_t>& el, VolumeDescriptor& dd,
151 : DataContainer<data_t>& dc, MinMaxFunction<data_t> clipping)
152 20 : {
153 20 : const Vec3i _center = el.getCenter();
154 20 : data_t _amplit = el.getAmplitude();
155 :
156 20 : auto strides = dd.getProductOfCoefficientsPerDimension();
157 :
158 20 : auto [minX, minY, minZ, maxX, maxY, maxZ] = clipping(boundingBox<data_t>(
159 20 : el.getRoundMaxHalfWidth(), _center, dd.getNumberOfCoefficientsPerDimension()));
160 :
161 20 : index_t minXSchifted = index_t(minX) + _center[INDEX_X];
162 20 : index_t minYSchifted = index_t(minY) + _center[INDEX_Y];
163 20 : index_t minZSchifted = index_t(minZ) + _center[INDEX_Z];
164 :
165 20 : Vec3X<data_t> idx(3);
166 20 : Vec3i idx_shifted(3);
167 :
168 20 : idx << minX, minY, minZ;
169 20 : idx_shifted << minXSchifted, minYSchifted, minZSchifted;
170 :
171 20 : Vec3X<data_t> rotated(3);
172 :
173 20 : index_t xOffset = 0;
174 20 : index_t yOffset = 0;
175 20 : index_t zOffset = 0;
176 :
177 267 : for (; idx[INDEX_Z] <= maxZ; idx[INDEX_Z]++, idx_shifted[INDEX_Z]++) {
178 247 : zOffset = idx_shifted[INDEX_Z] * strides[INDEX_Z];
179 4902 : for (; idx[INDEX_Y] <= maxY; idx[INDEX_Y]++, idx_shifted[INDEX_Y]++) {
180 4655 : yOffset = idx_shifted[INDEX_Y] * strides[INDEX_Y];
181 104077 : for (; idx[INDEX_X] <= maxX; idx[INDEX_X]++, idx_shifted[INDEX_X]++) {
182 99422 : xOffset = idx_shifted[INDEX_X] * strides[INDEX_X];
183 99422 : rotated = el.getInvRotationMatrix() * idx;
184 99422 : if (el.isInEllipsoid(rotated)) {
185 13991 : blend<b>(dc, xOffset + yOffset + zOffset, _amplit);
186 13991 : }
187 99422 : }
188 : // reset X
189 4655 : idx[INDEX_X] = minX;
190 4655 : idx_shifted[INDEX_X] = minXSchifted;
191 4655 : }
192 : // reset Y
193 247 : idx[INDEX_Y] = minY;
194 247 : idx_shifted[INDEX_Y] = minYSchifted;
195 247 : }
196 20 : };
197 :
198 : template <typename data_t>
199 18 : auto noClipping = [](std::array<data_t, 6> minMax) { return minMax; };
200 :
201 : template <Blending b, typename data_t>
202 : void rasterize(Ellipsoid<data_t>& el, VolumeDescriptor& dd, DataContainer<data_t>& dc)
203 97 : {
204 97 : if (el.isRotated()) {
205 18 : rasterizeWithClipping<b, data_t>(el, dd, dc, noClipping<data_t>);
206 79 : } else {
207 79 : rasterizeNoRotation<b, data_t>(el, dd, dc);
208 79 : }
209 97 : };
210 :
211 : // ------------------------------------------
212 : // explicit template instantiation
213 : template class Ellipsoid<float>;
214 : template class Ellipsoid<double>;
215 :
216 : template void rasterize<Blending::ADDITION, float>(Ellipsoid<float>& el, VolumeDescriptor& dd,
217 : DataContainer<float>& dc);
218 : template void rasterize<Blending::ADDITION, double>(Ellipsoid<double>& el, VolumeDescriptor& dd,
219 : DataContainer<double>& dc);
220 :
221 : template void rasterizeWithClipping<Blending::ADDITION, float>(Ellipsoid<float>& el,
222 : VolumeDescriptor& dd,
223 : DataContainer<float>& dc,
224 : MinMaxFunction<float> clipping);
225 : template void rasterizeWithClipping<Blending::ADDITION, double>(
226 : Ellipsoid<double>& el, VolumeDescriptor& dd, DataContainer<double>& dc,
227 : MinMaxFunction<double> clipping);
228 :
229 : template void rasterize<Blending::OVERWRITE, float>(Ellipsoid<float>& el, VolumeDescriptor& dd,
230 : DataContainer<float>& dc);
231 : template void rasterize<Blending::OVERWRITE, double>(Ellipsoid<double>& el,
232 : VolumeDescriptor& dd,
233 : DataContainer<double>& dc);
234 :
235 : template void rasterizeWithClipping<Blending::OVERWRITE, float>(Ellipsoid<float>& el,
236 : VolumeDescriptor& dd,
237 : DataContainer<float>& dc,
238 : MinMaxFunction<float> clipping);
239 : template void rasterizeWithClipping<Blending::OVERWRITE, double>(
240 : Ellipsoid<double>& el, VolumeDescriptor& dd, DataContainer<double>& dc,
241 : MinMaxFunction<double> clipping);
242 :
243 : } // namespace elsa::phantoms
|
