STIR 6.4.0
LORCoordinates.inl
Go to the documentation of this file.
1//
2//
3
13/*
14 Copyright (C) 2004- 2013, Hammersmith Imanet Ltd
15 Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics
16 Copyright (C) 2018, University College London
17 This file is part of STIR.
18
19 SPDX-License-Identifier: Apache-2.0
20
21 See STIR/LICENSE.txt for details
22*/
23
24#include "stir/modulo.h"
25#include "stir/Succeeded.h"
26#include <algorithm>
27START_NAMESPACE_STIR
28
29/********************************************************************
30 The simple and boring constructors
31********************************************************************/
32template <class coordT>
33PointOnCylinder<coordT>::PointOnCylinder()
34 : _psi(0) // set to value to avoid assert
35{}
36
37template <class coordT>
38PointOnCylinder<coordT>::PointOnCylinder(const coordT z, const coordT psi)
39 : _z(z),
40 _psi(psi)
41{}
42
43template <class coordT>
44LORInCylinderCoordinates<coordT>::LORInCylinderCoordinates(const coordT radius)
45 : _radius(radius)
46{}
47
48template <class coordT>
49LORInCylinderCoordinates<coordT>::LORInCylinderCoordinates(const PointOnCylinder<coordT>& p1,
50 const PointOnCylinder<coordT>& p2,
51 const coordT radius)
52 : _p1(p1),
53 _p2(p2),
54 _radius(radius)
55{}
56
57template <class coordT>
58LORAs2Points<coordT>::LORAs2Points()
59{}
60
61template <class coordT>
62LORAs2Points<coordT>::LORAs2Points(const CartesianCoordinate3D<coordT>& p1, const CartesianCoordinate3D<coordT>& p2)
63 : _p1(p1),
64 _p2(p2)
65{}
66
67template <class coordT>
68LORInAxialAndSinogramCoordinates<coordT>::LORInAxialAndSinogramCoordinates(const coordT radius)
69 : LORCylindricalCoordinates_z_and_radius<coordT>(radius),
70 _phi(0),
71 _s(0),
72 _swapped(false) // set _phi,_s to value to avoid assert
73{
74 check_state();
75}
76
77template <class coordT>
78LORInAxialAndSinogramCoordinates<coordT>::LORInAxialAndSinogramCoordinates(
79 const coordT z1, const coordT z2, const coordT phi, const coordT s, const coordT radius, const bool swapped)
80 : LORCylindricalCoordinates_z_and_radius<coordT>(z1, z2, radius),
81 _phi(to_0_2pi(phi)),
82 _s(s),
83 _swapped(swapped)
84{
85 if (_phi >= _PI)
86 {
87 _phi -= coordT(_PI);
88 _s = -_s;
89 std::swap(private_base_type::_z1, private_base_type::_z2);
90 _swapped = !swapped;
91 }
92 check_state();
93}
94
95template <class coordT>
96LORInAxialAndNoArcCorrSinogramCoordinates<coordT>::LORInAxialAndNoArcCorrSinogramCoordinates(const coordT radius)
97 : LORCylindricalCoordinates_z_and_radius<coordT>(radius),
98 _phi(0),
99 _beta(0),
100 _swapped(false) // set _phi,_beta to value to avoid assert
101
102{
103 check_state();
104}
105
106template <class coordT>
107LORInAxialAndNoArcCorrSinogramCoordinates<coordT>::LORInAxialAndNoArcCorrSinogramCoordinates(
108 const coordT z1, const coordT z2, const coordT phi, const coordT beta, const coordT radius, const bool swapped)
109 : LORCylindricalCoordinates_z_and_radius<coordT>(z1, z2, radius),
110 _phi(to_0_2pi(phi)),
111 _beta(beta),
112 _swapped(swapped)
113{
114 if (_phi >= _PI)
115 {
116 _phi -= coordT(_PI);
117 _beta = -_beta;
118 std::swap(private_base_type::_z1, private_base_type::_z2);
119 _swapped = !swapped;
120 }
121 check_state();
122}
123
124/********************************************************************
125 more constructors
126 These convert from one class to another
127********************************************************************/
128template <class coordT>
129LORInCylinderCoordinates<coordT>::LORInCylinderCoordinates(const LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& na_coords)
130 : _radius(na_coords.radius())
131{
132 _p1.z() = na_coords.z1();
133 _p2.z() = na_coords.z2();
134 _p1.psi() = to_0_2pi(na_coords.phi() + na_coords.beta());
135 _p2.psi() = to_0_2pi(na_coords.phi() - na_coords.beta() + static_cast<coordT>(_PI));
136 if (na_coords.is_swapped())
137 std::swap(_p1, _p2);
138 check_state();
139}
140
141template <class coordT>
142LORInCylinderCoordinates<coordT>::LORInCylinderCoordinates(const LORInAxialAndSinogramCoordinates<coordT>& coords)
143 : _radius(coords.radius())
144{
145 _p1.z() = coords.z1();
146 _p2.z() = coords.z2();
147 _p1.psi() = to_0_2pi(coords.phi() + coords.beta());
148 _p2.psi() = to_0_2pi(coords.phi() - coords.beta() + static_cast<coordT>(_PI));
149 if (coords.is_swapped())
150 std::swap(_p1, _p2);
151 check_state();
152}
153
154template <class coordT>
155static inline void
156get_sino_coords(
157 coordT& _z1, coordT& _z2, coordT& _phi, coordT& _beta, bool& _swapped, const LORInCylinderCoordinates<coordT>& cyl_coords)
158{
159 _beta = to_0_2pi((cyl_coords.p1().psi() - cyl_coords.p2().psi() + static_cast<coordT>(_PI)) / 2);
160 if (_beta > _PI)
161 _beta -= static_cast<coordT>(2 * _PI);
162 _phi = to_0_2pi((cyl_coords.p1().psi() + cyl_coords.p2().psi() - static_cast<coordT>(_PI)) / 2);
163
164 /* beta is now between -Pi and Pi, phi between 0 and 2Pi
165 Now bring into standard range and set z accordingly */
166 if (_phi < static_cast<coordT>(_PI))
167 {
168 if (_beta >= static_cast<coordT>(_PI) / 2)
169 {
170 _beta = static_cast<coordT>(_PI) - _beta;
171 _z2 = cyl_coords.p1().z();
172 _z1 = cyl_coords.p2().z();
173 _swapped = true;
174 }
175 else if (_beta < -static_cast<coordT>(_PI) / 2)
176 {
177 _beta = -static_cast<coordT>(_PI) - _beta;
178 _z2 = cyl_coords.p1().z();
179 _z1 = cyl_coords.p2().z();
180 _swapped = false;
181 }
182 else
183
184 {
185 _z1 = cyl_coords.p1().z();
186 _z2 = cyl_coords.p2().z();
187 _swapped = false;
188 }
189 }
190 else
191 {
192 _phi -= static_cast<coordT>(_PI);
193 assert(_phi >= 0);
194 if (_beta >= static_cast<coordT>(_PI) / 2)
195 {
196 _beta -= static_cast<coordT>(_PI);
197 _z1 = cyl_coords.p1().z();
198 _z2 = cyl_coords.p2().z();
199 _swapped = false;
200 }
201 else if (_beta < -static_cast<coordT>(_PI) / 2)
202 {
203 _beta += static_cast<coordT>(_PI);
204 _z1 = cyl_coords.p1().z();
205 _z2 = cyl_coords.p2().z();
206 _swapped = true;
207 }
208 else
209 {
210 _beta *= -1;
211 _z2 = cyl_coords.p1().z();
212 _z1 = cyl_coords.p2().z();
213 _swapped = true;
214 }
215 }
216 assert(_phi >= 0);
217 assert(_phi <= _PI);
218 assert(_beta >= -_PI / 2);
219 assert(_beta <= _PI / 2);
220}
221
222template <class coordT>
223LORInAxialAndNoArcCorrSinogramCoordinates<coordT>::LORInAxialAndNoArcCorrSinogramCoordinates(
224 const LORInCylinderCoordinates<coordT>& cyl_coords)
225 : LORCylindricalCoordinates_z_and_radius<coordT>(cyl_coords.radius())
226{
227#ifndef NDEBUG
228 // set these to prevent assert breaking in check_state()
229 _phi = 0;
230 _beta = 0;
231#endif
232 get_sino_coords(z1(), z2(), _phi, _beta, _swapped, cyl_coords);
233 check_state();
234}
235
236template <class coordT>
237LORInAxialAndSinogramCoordinates<coordT>::LORInAxialAndSinogramCoordinates(const LORInCylinderCoordinates<coordT>& cyl_coords)
238 : LORCylindricalCoordinates_z_and_radius<coordT>(cyl_coords.radius())
239{
240 coordT beta;
241#ifndef NDEBUG
242 // set these to prevent assert breaking in check_state()
243 _phi = 0;
244 _s = 0;
245#endif
246 get_sino_coords(z1(), z2(), _phi, beta, _swapped, cyl_coords);
247 _s = this->_radius * sin(beta);
248 check_state();
249}
250
251template <class coordT>
252LORInAxialAndSinogramCoordinates<coordT>::LORInAxialAndSinogramCoordinates(
253 const LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& coords)
254 : LORCylindricalCoordinates_z_and_radius<coordT>(coords.z1(), coords.z2(), coords.radius()),
255 _phi(coords.phi()),
256 _s(coords.s()),
257 _swapped(coords.is_swapped())
258{
259 check_state();
260}
261
262template <class coordT>
263LORInAxialAndNoArcCorrSinogramCoordinates<coordT>::LORInAxialAndNoArcCorrSinogramCoordinates(
264 const LORInAxialAndSinogramCoordinates<coordT>& coords)
265 : LORCylindricalCoordinates_z_and_radius<coordT>(coords.z1(), coords.z2(), coords.radius()),
266 _phi(coords.phi()),
267 _beta(coords.beta()),
268 _swapped(coords.is_swapped())
269{
270 check_state();
271}
272
273#if 0
274template <class coordT>
275LORInAxialAndSinogramCoordinates<coordT>::
276LORInAxialAndSinogramCoordinates(const LORAs2Points<coordT>& coords)
277 : LORInAxialAndSinogramCoordinates(LORInCylinderCoordinates<coordT>(coords))
278{}
279
280template <class coordT>
281LORInAxialAndNoArcCorrSinogramCoordinates<coordT>::
282LORInAxialAndNoArcCorrSinogramCoordinates(const LORAs2Points<coordT>& coords)
283 : LORInAxialAndNoArcCorrSinogramCoordinates(LORInCylinderCoordinates<coordT>(coords))
284{}
285#endif
286
287template <class coordT>
288LORAs2Points<coordT>::LORAs2Points(const LORInCylinderCoordinates<coordT>& cyl_coords)
289{
290 // make sure the return values are in STIR coordinates
291 p1().z() = cyl_coords.p1().z();
292 p1().y() = -cyl_coords.radius() * cos(cyl_coords.p1().psi());
293 p1().x() = cyl_coords.radius() * sin(cyl_coords.p1().psi());
294
295 p2().z() = cyl_coords.p2().z();
296 p2().y() = -cyl_coords.radius() * cos(cyl_coords.p2().psi());
297 p2().x() = cyl_coords.radius() * sin(cyl_coords.p2().psi());
298}
299
300template <class coordT>
301LORAs2Points<coordT>::LORAs2Points(const LORInAxialAndSinogramCoordinates<coordT>& coords)
302{
303 *this = LORInCylinderCoordinates<coordT>(coords);
304}
305
306template <class coordT>
307LORAs2Points<coordT>::LORAs2Points(const LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& coords)
308{
309 *this = LORInCylinderCoordinates<coordT>(coords);
310}
311
312template <class coordT1, class coordT2>
313Succeeded
314find_LOR_intersections_with_cylinder(LORAs2Points<coordT1>& intersection_coords,
315 const LORAs2Points<coordT2>& coords,
316 const double radius)
317{
318 const CartesianCoordinate3D<coordT2>& c1 = coords.p1();
319 const CartesianCoordinate3D<coordT2>& c2 = coords.p2();
320
321 const CartesianCoordinate3D<coordT2> d = c2 - c1;
322 /* parametrisation of LOR is
323 c = l*d+c1
324 l has to be such that c.x^2 + c.y^2 = R^2
325 i.e.
326 (l*d.x+c1.x)^2+(l*d.y+c1.y)^2==R^2
327 l^2*(d.x^2+d.y^2) + 2*l*(d.x*c1.x + d.y*c1.y) + c1.x^2+c2.y^2-R^2==0
328 write as a*l^2+2*b*l+e==0
329 l = (-b +- sqrt(b^2-a*e))/a
330 argument of sqrt simplifies to
331 R^2*(d.x^2+d.y^2)-(d.x*c1.y-d.y*c1.x)^2
332 */
333 const double a = square(d.x()) + square(d.y());
334 const double b = d.x() * c1.x() + d.y() * c1.y();
335 const double e = square(c1.x()) + square(c1.y()) - square(radius);
336 double argsqrt = square(b) - a * e;
337 if (argsqrt <= 0)
338 {
339 return Succeeded::no; // LOR is outside detector radius
340 }
341 const coordT2 root = static_cast<coordT2>(sqrt(argsqrt));
342
343 auto l1 = static_cast<coordT2>((-b - root) / a);
344 auto l2 = static_cast<coordT2>((-b + root) / a);
345
346 intersection_coords.p1() = d * l1 + c1;
347 intersection_coords.p2() = d * l2 + c1;
348 assert(fabs(square(intersection_coords.p1().x()) + square(intersection_coords.p1().y()) - square(radius))
349 < square(radius) * 10.E-5);
350 assert(fabs(square(intersection_coords.p2().x()) + square(intersection_coords.p2().y()) - square(radius))
351 < square(radius) * 10.E-5);
352 return Succeeded::yes;
353}
354
355template <class coordT1, class coordT2>
356Succeeded
357find_LOR_intersections_with_cylinder(LORInCylinderCoordinates<coordT1>& cyl_coords,
358 const LORAs2Points<coordT2>& cart_coords,
359 const double radius)
360{
361 LORAs2Points<coordT1> intersection_coords;
362 if (find_LOR_intersections_with_cylinder(intersection_coords, cart_coords, radius) == Succeeded::no)
363 return Succeeded::no;
364
365 const CartesianCoordinate3D<coordT1>& c1 = intersection_coords.p1();
366 const CartesianCoordinate3D<coordT1>& c2 = intersection_coords.p2();
367 cyl_coords.reset(static_cast<float>(radius));
368
369 cyl_coords.p1().psi() = from_min_pi_plus_pi_to_0_2pi(static_cast<coordT1>(atan2(c1.x(), -c1.y())));
370 cyl_coords.p2().psi() = from_min_pi_plus_pi_to_0_2pi(static_cast<coordT1>(atan2(c2.x(), -c2.y())));
371 cyl_coords.p1().z() = static_cast<coordT1>(c1.z());
372 cyl_coords.p2().z() = static_cast<coordT1>(c2.z());
373
374 return Succeeded::yes;
375}
376
377template <class coordT1, class coordT2>
378Succeeded
379find_LOR_intersections_with_cylinder(LORInAxialAndNoArcCorrSinogramCoordinates<coordT1>& lor,
380 const LORAs2Points<coordT2>& cart_coords,
381 const double radius)
382{
383 LORInCylinderCoordinates<coordT1> cyl_coords;
384 if (find_LOR_intersections_with_cylinder(cyl_coords, cart_coords, radius) == Succeeded::no)
385 return Succeeded::no;
386 lor = cyl_coords;
387 return Succeeded::yes;
388}
389
390template <class coordT1, class coordT2>
391Succeeded
392find_LOR_intersections_with_cylinder(LORInAxialAndSinogramCoordinates<coordT1>& lor,
393 const LORAs2Points<coordT2>& cart_coords,
394 const double radius)
395{
396 LORInCylinderCoordinates<coordT1> cyl_coords;
397 if (find_LOR_intersections_with_cylinder(cyl_coords, cart_coords, radius) == Succeeded::no)
398 return Succeeded::no;
399 lor = cyl_coords;
400 return Succeeded::yes;
401}
402
403template <class coordT>
404Succeeded
405LORAs2Points<coordT>::change_representation(LORInCylinderCoordinates<coordT>& lor, const double radius) const
406{
407 return find_LOR_intersections_with_cylinder(lor, *this, radius);
408}
409
410template <class coordT>
411Succeeded
412LORAs2Points<coordT>::change_representation(LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& lor, const double radius) const
413{
414 return find_LOR_intersections_with_cylinder(lor, *this, radius);
415}
416
417template <class coordT>
418Succeeded
419LORAs2Points<coordT>::change_representation(LORInAxialAndSinogramCoordinates<coordT>& lor, const double radius) const
420{
421 return find_LOR_intersections_with_cylinder(lor, *this, radius);
422}
423
424template <class coordT>
425Succeeded
426LORAs2Points<coordT>::get_intersections_with_cylinder(LORAs2Points<coordT>& lor, const double radius) const
427{
428 return find_LOR_intersections_with_cylinder(lor, *this, radius);
429}
430
431template <class coordT>
432Succeeded
433LORAs2Points<coordT>::change_representation_for_block(LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& lor,
434 const double radius) const
435{
436 const CartesianCoordinate3D<coordT>& c1 = this->p1();
437 const CartesianCoordinate3D<coordT>& c2 = this->p2();
438
439 // To check if LOR is inside the detector
440 const CartesianCoordinate3D<coordT> d = c2 - c1;
441 const double dxy2 = (square(d.x()) + square(d.y()));
442 const double argsqrt = (square(radius) * dxy2 - square(d.x() * c1.y() - d.y() * c1.x()));
443
444 LORInCylinderCoordinates<coordT> cyl_coords;
445 cyl_coords.reset(static_cast<float>(radius));
446
447 cyl_coords.p1().psi() = from_min_pi_plus_pi_to_0_2pi(static_cast<coordT>(atan2(c1.x(), -c1.y())));
448 cyl_coords.p2().psi() = from_min_pi_plus_pi_to_0_2pi(static_cast<coordT>(atan2(c2.x(), -c2.y())));
449 cyl_coords.p1().z() = static_cast<coordT>(c1.z());
450 cyl_coords.p2().z() = static_cast<coordT>(c2.z());
451 lor = cyl_coords;
452
453 if (argsqrt <= 0)
454 return Succeeded::no; // LOR is outside detector radius
455 else
456 return Succeeded::yes;
457}
458
459#define DEFINE_LOR_GET_FUNCTIONS(TYPE) \
460 template <class coordT> \
461 Succeeded TYPE<coordT>::change_representation(LORInCylinderCoordinates<coordT>& lor, const double radius) const \
462 { \
463 lor = *this; \
464 return lor.set_radius(static_cast<coordT>(radius)); \
465 } \
466 \
467 template <class coordT> \
468 Succeeded TYPE<coordT>::change_representation(LORInAxialAndNoArcCorrSinogramCoordinates<coordT>& lor, const double radius) \
469 const \
470 { \
471 lor = *this; \
472 return lor.set_radius(static_cast<coordT>(radius)); \
473 } \
474 \
475 template <class coordT> \
476 Succeeded TYPE<coordT>::change_representation(LORInAxialAndSinogramCoordinates<coordT>& lor, const double radius) const \
477 { \
478 lor = *this; \
479 return lor.set_radius(static_cast<coordT>(radius)); \
480 } \
481 \
482 template <class coordT> \
483 Succeeded TYPE<coordT>::get_intersections_with_cylinder(LORAs2Points<coordT>& lor, const double radius) const \
484 { \
485 LORAs2Points<coordT> actual_lor = *this; \
486 if (find_LOR_intersections_with_cylinder(lor, actual_lor, radius) == Succeeded::no) \
487 return Succeeded::no; \
488 return Succeeded::yes; \
489 }
490
491DEFINE_LOR_GET_FUNCTIONS(LORInCylinderCoordinates)
492DEFINE_LOR_GET_FUNCTIONS(LORInAxialAndNoArcCorrSinogramCoordinates)
493DEFINE_LOR_GET_FUNCTIONS(LORInAxialAndSinogramCoordinates)
494
495#undef DEFINE_LOR_GET_FUNCTIONS
496
497END_NAMESPACE_STIR
Succeeded.h
Declaration of class stir::Succeeded.
stir::CartesianCoordinate3D
a templated class for 3-dimensional coordinates.
Definition CartesianCoordinate3D.h:53
stir::LORAs2Points
A class for LORs.
Definition LORCoordinates.h:296
stir::LORAs2Points::change_representation_for_block
Succeeded change_representation_for_block(LORInAxialAndNoArcCorrSinogramCoordinates< coordT > &, const double radius) const
Calculate intersections with block. Used in: ProjDataInfoBlocksOnCylindrical::get_LOR.
Definition LORCoordinates.inl:433
stir::LORCylindricalCoordinates_z_and_radius
An internal class for LORs. Do not use.
Definition LORCoordinates.h:118
stir::LORInAxialAndNoArcCorrSinogramCoordinates
A class for LORs.
Definition LORCoordinates.h:511
stir::LORInAxialAndNoArcCorrSinogramCoordinates::is_swapped
bool is_swapped() const override
Return if the LOR direction is opposite from normal.
Definition LORCoordinates.h:571
stir::LORInAxialAndSinogramCoordinates
A class for LORs.
Definition LORCoordinates.h:367
stir::LORInCylinderCoordinates
A class for LORs.
Definition LORCoordinates.h:185
stir::Succeeded
a class containing an enumeration type that can be used by functions to signal successful operation o...
Definition Succeeded.h:44
stir::find_LOR_intersections_with_cylinder
Succeeded find_LOR_intersections_with_cylinder(LORInCylinderCoordinates< coordT1 > &, const LORAs2Points< coordT2 > &, const double radius)
Given an LOR, find its intersection with a (infintely long) cylinder.
Definition LORCoordinates.inl:357
_PI
#define _PI
The constant pi to high precision.
Definition common.h:141
stir::square
NUMBER square(const NUMBER &x)
returns the square of a number, templated.
Definition common.h:154
modulo.h
defines stir::modulo() and related functions
stir::from_min_pi_plus_pi_to_0_2pi
FloatOrDouble from_min_pi_plus_pi_to_0_2pi(const FloatOrDouble phi)
A function to convert an angle from one range to another.
Definition modulo.h:118
stir::to_0_2pi
FloatOrDouble to_0_2pi(const FloatOrDouble phi)
Convert angle to standard range.
Definition modulo.h:138
Generated on Wed May 27 2026 17:04:05 for STIR by doxygen 1.13.2