Line data Source code
1 : #include "DataContainer.h"
2 : #include "IdenticalBlocksDescriptor.h"
3 : #include "Math.hpp"
4 : #include "SphericalCoefficientsDescriptor.h"
5 : #include "TypeCasts.hpp"
6 : #include "elsaDefines.h"
7 : #include "SphericalHarmonicsTransform.h"
8 : #include "XGIDetectorDescriptor.h"
9 :
10 : namespace elsa::axdt
11 : {
12 :
13 : namespace detail
14 : {
15 :
16 : template <typename data_t>
17 : inline std::array<data_t, 15> exactCoeffsMalecki(const DirVec<data_t> ray,
18 : const DirVec<data_t> sensitivity)
19 719 : {
20 719 : using math::sq;
21 :
22 719 : const data_t lx = ray[0], ly = ray[1], lz = ray[2];
23 719 : const data_t tx = sensitivity[0], ty = sensitivity[1], tz = sensitivity[2];
24 :
25 : // Thanks, Mathematica!
26 719 : data_t h_00 =
27 719 : (4 * sqrt(pi<data_t>)
28 719 : * (1 + sq(ty) - 2 * lx * lz * tx * tz + sq(tz) - 2 * ly * ty * (lx * tx + lz * tz)
29 719 : - sq(ly) * (-1 + 2 * sq(ty) + sq(tz)) - sq(lz) * (-1 + sq(ty) + 2 * sq(tz))))
30 719 : / 15.;
31 :
32 719 : data_t h_2_2 = (4 * sqrt(pi<data_t> / 15.)
33 719 : * (2 * tx * ((2 + sq(lz)) * ty - ly * lz * tz)
34 719 : + lx * (-2 * lz * ty * tz + ly * (-3 + 2 * sq(tz)))))
35 719 : / 7.;
36 :
37 719 : data_t h_2_1 = (4 * sqrt(pi<data_t> / 15.)
38 719 : * (2 * sq(ly) * ty * tz + 2 * ty * (lx * lz * tx + (-3 + sq(lz)) * tz)
39 719 : + ly * (2 * lx * tx * tz + lz * (1 + 2 * sq(ty) + 2 * sq(tz)))))
40 719 : / 7.;
41 :
42 719 : data_t h_20 =
43 719 : (2 * sqrt(pi<data_t> / 5.)
44 719 : * (-1 - 4 * sq(ly) - 7 * sq(lz) + 8 * lx * ly * tx * ty
45 719 : + 4 * (-1 + 2 * sq(ly) + sq(lz)) * sq(ty) - 4 * lz * (lx * tx + ly * ty) * tz
46 719 : + 2 * (7 + 2 * sq(ly) - 2 * sq(lz)) * sq(tz)))
47 719 : / 21.;
48 :
49 719 : data_t h_21 = (-4 * sqrt(pi<data_t> / 15.)
50 719 : * (2 * tx * (-(ly * lz * ty) + (2 + sq(ly)) * tz)
51 719 : + lx * (lz * (-3 + 2 * sq(ty)) - 2 * ly * ty * tz)))
52 719 : / 7.;
53 :
54 719 : data_t h_22 =
55 719 : (-2 * sqrt(pi<data_t> / 15.)
56 719 : * (-1 + 8 * sq(ty) + 4 * lx * lz * tx * tz - 4 * ly * lz * ty * tz + 2 * sq(tz)
57 719 : + sq(ly) * (-6 + 4 * sq(tz)) + sq(lz) * (-5 + 4 * sq(ty) + 4 * sq(tz))))
58 719 : / 7.;
59 :
60 719 : data_t h_4_4 =
61 719 : (4 * sqrt(pi<data_t> / 35.)
62 719 : * ((-1 + 2 * sq(ly) + sq(lz)) * tx * ty + lx * ly * (-1 + 2 * sq(ty) + sq(tz))))
63 719 : / 3.;
64 :
65 719 : data_t h_4_3 = (-2 * sqrt((2 * pi<data_t>) / 35.)
66 719 : * (2 * sq(ly) * ty * tz + ty * (-2 * lx * lz * tx + (-1 + sq(lz)) * tz)
67 719 : + ly * (-2 * lx * tx * tz + lz * (-1 + 2 * sq(ty) + sq(tz)))))
68 719 : / 3.;
69 :
70 719 : data_t h_4_2 = (-4 * sqrt(pi<data_t> / 5.)
71 719 : * (tx * ((-1 + 3 * sq(lz)) * ty + 4 * ly * lz * tz)
72 719 : + lx * (4 * lz * ty * tz + ly * (-1 + 3 * sq(tz)))))
73 719 : / 21.;
74 :
75 719 : data_t h_4_1 =
76 719 : (-2 * sqrt((2 * pi<data_t>) / 5.)
77 719 : * (2 * sq(ly) * ty * tz + ty * (2 * lx * lz * tx + (1 - 5 * sq(lz)) * tz)
78 719 : + ly * (2 * lx * tx * tz + lz * (1 + 2 * sq(ty) - 5 * sq(tz)))))
79 719 : / 21.;
80 :
81 719 : data_t h_40 =
82 719 : (2 * sqrt(pi<data_t>)
83 719 : * (-3 + 2 * sq(ly) + 7 * sq(lz) - 4 * lx * ly * tx * ty
84 719 : - 2 * (-1 + 2 * sq(ly) + sq(lz)) * sq(ty) + 16 * lz * (lx * tx + ly * ty) * tz
85 719 : + (7 - 2 * sq(ly) - 19 * sq(lz)) * sq(tz)))
86 719 : / 105.;
87 :
88 719 : data_t h_41 = (2 * sqrt((2 * pi<data_t>) / 5.)
89 719 : * (tx * (-2 * ly * lz * ty + (-3 + 2 * sq(ly) + 7 * sq(lz)) * tz)
90 719 : + lx * (-2 * ly * ty * tz + lz * (-3 + 2 * sq(ty) + 7 * sq(tz)))))
91 719 : / 21.;
92 :
93 719 : data_t h_42 =
94 719 : (4 * sqrt(pi<data_t> / 5.)
95 719 : * (1 - sq(ty) - 4 * lx * lz * tx * tz + 4 * ly * lz * ty * tz - 2 * sq(tz)
96 719 : + sq(ly) * (-1 + 3 * sq(tz)) + sq(lz) * (-2 + 3 * sq(ty) + 3 * sq(tz))))
97 719 : / 21.;
98 :
99 719 : data_t h_43 = (-2 * sqrt((2 * pi<data_t>) / 35.)
100 719 : * (tx * (2 * ly * lz * ty + (-1 + 2 * sq(ly) + sq(lz)) * tz)
101 719 : + lx * (2 * ly * ty * tz + lz * (-1 + 2 * sq(ty) + sq(tz)))))
102 719 : / 3.;
103 :
104 719 : data_t h_44 = (-2 * sqrt(pi<data_t> / 35.)
105 719 : * (-4 * lx * ly * tx * ty + 2 * sq(ly) * (-1 + 2 * sq(ty) + sq(tz))
106 719 : + (-1 + sq(lz)) * (-1 + 2 * sq(ty) + sq(tz))))
107 719 : / 3.;
108 :
109 719 : return {h_00, h_2_2, h_2_1, h_20, h_21, h_22, h_4_4, h_4_3,
110 719 : h_4_2, h_4_1, h_40, h_41, h_42, h_43, h_44};
111 719 : }
112 : } // namespace detail
113 :
114 : template <typename data_t>
115 : DataContainer<data_t> exactWeightingFunction(const XGIDetectorDescriptor& detectorDescriptor,
116 : const Symmetry symmetry, const index_t degree)
117 36 : {
118 :
119 36 : const index_t poses = detectorDescriptor.getNumberOfCoefficientsPerDimension()[2];
120 :
121 : // setup complete block descriptor
122 36 : SphericalCoefficientsDescriptor resultDescriptor{detectorDescriptor, symmetry, degree};
123 :
124 : // setup factors block
125 : // Result has shape (detX×detY)×#poses×#coeffs
126 : // Set to zero as odd weights are zero and appear in Symmetry::regular
127 36 : auto weights = zeros<data_t>(resultDescriptor);
128 :
129 36 : if (detectorDescriptor.isParallelBeam()) {
130 :
131 36 : #pragma omp parallel for
132 36 : for (index_t pose = 0; pose < poses; ++pose) {
133 :
134 : // obtain geometry object for current image
135 0 : const Geometry camera = detectorDescriptor.getGeometryAt(pose);
136 :
137 : // beam direction: last column of rotation matrix (scaling is forced to be
138 : // isotropic, z-direction/viewing axis)
139 0 : const DirVec<data_t> l =
140 0 : camera.getRotationMatrix().transpose().col(2).cast<data_t>();
141 :
142 : // sensitivity direction: first column of rotation matrix (scaling is forced to
143 : // be isotropic, x-direction/horizontal axis)
144 0 : const DirVec<data_t> t =
145 0 : (camera.getRotationMatrix().transpose() * detectorDescriptor.getSensDir())
146 0 : .cast<data_t>();
147 :
148 0 : const auto coeffs = detail::exactCoeffsMalecki(l, t);
149 :
150 700 : if (symmetry == Symmetry::regular) {
151 : // Truncation at l = 4 is exact
152 700 : const auto lMax = std::min(degree, index_t{4});
153 700 : auto j = 0;
154 2114 : for (int l = 0; l <= lMax; l += 2) {
155 6658 : for (int m = -l; m <= l; ++m) {
156 5244 : weights.getBlock(l * l + l + m).slice(pose).fill(coeffs[j]);
157 5244 : j++;
158 5244 : }
159 1414 : }
160 :
161 >1844*10^16 : } else if (symmetry == Symmetry::even) {
162 80 : const auto coeffCount = std::min(
163 80 : SphericalCoefficientsDescriptor::coefficientCount(symmetry, degree),
164 80 : index_t{15});
165 1273 : for (auto j = 0; j < coeffCount; ++j) {
166 1193 : weights.getBlock(j).slice(pose).fill(coeffs[j]);
167 1193 : }
168 80 : }
169 0 : }
170 36 : } else {
171 0 : const index_t detX = detectorDescriptor.getNumberOfCoefficientsPerDimension()[0];
172 0 : const index_t detY = detectorDescriptor.getNumberOfCoefficientsPerDimension()[1];
173 :
174 0 : #pragma omp parallel for
175 0 : for (index_t pose = 0; pose < poses; ++pose) {
176 : // set index to base for currently processed image and direction
177 0 : index_t idx = pose * detX * detY;
178 :
179 : // obtain geometry object for current image
180 0 : const auto camera = detectorDescriptor.getGeometryAt(pose);
181 :
182 : // sensitivity direction: first column of rotation matrix (scaling is
183 : // forced to be isotropic, x-direction/horizontal axis)
184 0 : const DirVec<data_t> t =
185 0 : (camera.getRotationMatrix().transpose() * detectorDescriptor.getSensDir())
186 0 : .cast<data_t>();
187 :
188 : // traverse image
189 0 : for (index_t y = 0; y < detY; ++y) {
190 0 : for (index_t x = 0; x < detX; ++x) {
191 : // sampling direction: set above, normalized by design
192 : // beam direction: get ray as provided by projection matrix
193 :
194 0 : const IndexVector_t coord{{x, y, pose}};
195 :
196 0 : const auto ray = detectorDescriptor.computeRayFromDetectorCoord(coord);
197 :
198 0 : const DirVec<data_t> l = ray.direction().cast<data_t>();
199 :
200 0 : const auto coeffs = detail::exactCoeffsMalecki(l, t);
201 :
202 0 : if (symmetry == Symmetry::regular) {
203 : // Truncation at l = 4 is exact
204 0 : auto lMax = std::min(degree, index_t{4});
205 0 : auto j = 0;
206 0 : for (int l = 0; l <= lMax; l += 2) {
207 0 : for (int m = -l; m <= l; ++m) {
208 : // As result has shape (detX×detY)×#poses×#coeffs, the j-th
209 : // block corresponds to h_j, which is the same for all
210 : // pixels in the detector at a given pose, hence we fill
211 : // slice(pose)
212 0 : weights.getBlock(l * l + l + m)[idx] = coeffs[j];
213 0 : j++;
214 0 : }
215 0 : }
216 :
217 0 : } else if (symmetry == Symmetry::even) {
218 0 : const auto coeffCount = std::min(
219 0 : SphericalCoefficientsDescriptor::coefficientCount(symmetry, degree),
220 0 : index_t{15});
221 0 : for (auto j = 0; j < coeffCount; ++j) {
222 0 : weights.getBlock(j)[idx] = coeffs[j];
223 0 : }
224 0 : }
225 0 : idx++;
226 0 : }
227 0 : }
228 0 : }
229 0 : }
230 :
231 36 : return weights;
232 36 : }
233 :
234 : template <typename data_t>
235 : DataContainer<data_t> evaluateMalecki(const XGIDetectorDescriptor& detectorDescriptor,
236 : const DirVecList<data_t>& directions)
237 28 : {
238 28 : using namespace axdt;
239 :
240 : // obtain number of sampling directions, number of measured images
241 28 : const index_t numDirs = directions.size();
242 28 : const index_t poses = detectorDescriptor.getNumberOfCoefficientsPerDimension()[2];
243 :
244 : // set up the block descriptor
245 28 : auto directionSpace = IdenticalBlocksDescriptor{numDirs, detectorDescriptor};
246 :
247 : // Shape (detX×detY)×#poses×#dirs
248 28 : auto sampledWeightFunction = empty<data_t>(directionSpace);
249 :
250 28 : if (detectorDescriptor.isParallelBeam()) {
251 :
252 1036 : for (index_t i = 0; i < numDirs; ++i) {
253 : // obtain sampling direction
254 1008 : DirVec<data_t> e = directions[asUnsigned(i)];
255 1008 : if (abs(e.norm() - 1) > 1e-5)
256 0 : throw std::invalid_argument("direction vector at index " + std::to_string(i)
257 0 : + " not normalized");
258 :
259 1008 : #pragma omp parallel for
260 1008 : for (index_t pose = 0; pose < poses; ++pose) {
261 :
262 : // obtain geometry object for current image
263 0 : const Geometry camera = detectorDescriptor.getGeometryAt(pose);
264 :
265 : // allocate the grating's sensitivity vector
266 0 : DirVec<data_t> t;
267 :
268 : // allocate helper objects for ray computation
269 0 : DirVec<data_t> l;
270 :
271 : // sampling direction: set above, normalized by design
272 :
273 : // beam direction: last column of rotation matrix (scaling is forced to be
274 : // isotropic, z-direction/viewing axis)
275 : // TODO: Why do we take the last row?
276 0 : l = camera.getRotationMatrix().row(2).transpose().cast<data_t>();
277 :
278 : // sensitivity direction: first column of rotation matrix (scaling is forced
279 : // to be isotropic, x-direction/horizontal axis)
280 0 : t = (camera.getRotationMatrix().transpose() * detectorDescriptor.getSensDir())
281 0 : .cast<data_t>();
282 :
283 : // compute the sensitivty function (|s × e|<e,t>)²
284 0 : auto h = math::sq(l.cross(e).norm() * e.dot(t));
285 :
286 : // h is the same for all pixels in the detector
287 : // we can therefore fill the whole slice (corresponding to a slab of shape
288 : // (detX×detY) with the same value
289 0 : sampledWeightFunction.getBlock(i).slice(pose).fill(h);
290 0 : }
291 1008 : }
292 28 : } else {
293 :
294 0 : const index_t px = detectorDescriptor.getNumberOfCoefficientsPerDimension()[0];
295 0 : const index_t py = detectorDescriptor.getNumberOfCoefficientsPerDimension()[1];
296 :
297 0 : for (index_t i = 0; i < numDirs; ++i) {
298 : // obtain sampling direction
299 0 : DirVec<data_t> e = directions[asUnsigned(i)];
300 :
301 0 : if (abs(e.norm() - 1) > 1e-5)
302 0 : throw std::invalid_argument("direction vector at index " + std::to_string(i)
303 0 : + " not normalized");
304 :
305 0 : #pragma omp parallel for
306 0 : for (index_t pose = 0; pose < poses; ++pose) {
307 : // set index to base for currently processed image and direction
308 0 : index_t idx = pose * px * py;
309 :
310 : // obtain geometry object for current image
311 0 : auto camera = detectorDescriptor.getGeometryAt(pose);
312 :
313 : // allocate the grating's sensitivity vector
314 0 : DirVec<data_t> t;
315 :
316 : // allocate helper objects for ray computation
317 0 : DirVec<data_t> s;
318 0 : IndexVector_t pt(3);
319 :
320 : // traverse image
321 0 : for (index_t y = 0; y < py; ++y) {
322 0 : for (index_t x = 0; x < px; ++x) {
323 : // sampling direction: set above, normalized by design
324 : // beam direction: get ray as provided by projection matrix
325 0 : pt << x, y, pose;
326 0 : auto ray = detectorDescriptor.computeRayFromDetectorCoord(pt);
327 0 : s = ray.direction().cast<data_t>();
328 :
329 : // sensitivity direction: first column of rotation matrix (scaling
330 : // is forced to be isotropic, x-direction/horizontal axis)
331 0 : t = (camera.getRotationMatrix().transpose()
332 0 : * detectorDescriptor.getSensDir())
333 0 : .cast<data_t>();
334 :
335 : // compute the factor: (|sxe|<e,t>)^2
336 : // apply the factor (the location is x/y/n, but there is no need to
337 : // compute that)
338 0 : const auto h = math::sq(s.cross(e).norm() * e.dot(t));
339 :
340 0 : sampledWeightFunction.getBlock(i)[idx++] = h;
341 0 : }
342 0 : }
343 0 : }
344 0 : }
345 0 : }
346 :
347 28 : return sampledWeightFunction;
348 28 : }
349 :
350 : template <typename data_t>
351 : DataContainer<data_t>
352 : numericalWeightingFunction(const XGIDetectorDescriptor& detectorDescriptor,
353 : const Symmetry symmetry, const index_t degree,
354 : const DirVecList<data_t>& directions)
355 28 : {
356 : // Shape of sampled weighting function: (detX×detY)×#poses×#dirs
357 28 : auto sampledMalecki = evaluateMalecki(detectorDescriptor, directions);
358 :
359 : // Rescale in order as to obtain the correct normalization for the legacy tests
360 28 : sampledMalecki *= sqrt(4 * pi<data_t> / directions.size());
361 :
362 : // Result has shape (detX×detY)×#poses×#coeffs
363 28 : SphericalCoefficientsDescriptor coefficientDescriptor{detectorDescriptor, symmetry, degree};
364 :
365 28 : SphericalHarmonicsTransform<data_t> evaluate{coefficientDescriptor, directions};
366 :
367 28 : return evaluate.applyAdjoint(sampledMalecki);
368 28 : }
369 :
370 : template DataContainer<float> axdt::exactWeightingFunction(const XGIDetectorDescriptor&,
371 : const Symmetry, const index_t);
372 : template DataContainer<double> axdt::exactWeightingFunction(const XGIDetectorDescriptor&,
373 : const Symmetry, const index_t);
374 :
375 : template DataContainer<float> axdt::numericalWeightingFunction(const XGIDetectorDescriptor&,
376 : const Symmetry, const index_t,
377 : const DirVecList<float>&);
378 :
379 : template DataContainer<double> axdt::numericalWeightingFunction(const XGIDetectorDescriptor&,
380 : const Symmetry, const index_t,
381 : const DirVecList<double>&);
382 :
383 : } // namespace elsa::axdt
|
