STIR  6.2.0
line_distances.h
Go to the documentation of this file.
1 //
2 //
3 
12 /*
13  Copyright (C) 2005- 2005, Hammersmith Imanet Ltd
14  Copyright (C) 2016, University of Hull
15 
16  This file is part of STIR.
17 
18  SPDX-License-Identifier: Apache-2.0
19 
20  See STIR/LICENSE.txt for details.
21 */
22 
24 #include "stir/LORCoordinates.h"
25 #include <cmath>
26 #ifdef BOOST_NO_STDC_NAMESPACE
27 namespace std
28 {
29 using ::sqrt;
30 using ::fabs;
31 } // namespace std
32 #endif
33 
34 START_NAMESPACE_STIR
35 
44 template <class coordT>
45 inline coordT
47  const LORAs2Points<coordT>& line0,
48  const LORAs2Points<coordT>& line1)
49 {
50  /* Rationale:
51  parametrise points on the lines as
52  r0 + a0 d0, r1 + a1 d1
53  where r0 is a point on the line, and d0 is its direction
54 
55  compute distance squared between those points, i.e.
56  norm(a0 d0 - a1 d1 + r0-r1)^2
57  This can be written using inner products of all the vectors.
58 
59  Minimising it (by computing derivatives) results in a linear eq in a0, a1:
60  a0 d0d0 == a1 d0d1 + r10d0,
61  a0 d0d1 == a1 d1d1 + r10d1
62  which is easily solved.
63  The half-way point can then be found by using
64  (r0 + a0 d0 + r1 + a1 d1)/2
65  for the a0,a1 found.
66  */
67  const CartesianCoordinate3D<coordT>& r0 = line0.p1();
68  const CartesianCoordinate3D<coordT>& r1 = line1.p1();
69  const CartesianCoordinate3D<coordT> r10 = r1 - r0;
70  const CartesianCoordinate3D<coordT> d0 = line0.p2() - line0.p1();
71  const CartesianCoordinate3D<coordT> d1 = line1.p2() - line1.p1();
72  const coordT d0d0 = inner_product(d0, d0);
73  const coordT d0d1 = inner_product(d0, d1);
74  const coordT d1d1 = inner_product(d1, d1);
75  const coordT r10d0 = inner_product(r10, d0);
76  const coordT r10r10 = inner_product(r10, r10);
77 
78  const coordT eps = d0d0 * 10E-5; // small number for comparisons
79 
80  const coordT denom = square(d0d1) - d0d0 * d1d1;
81  coordT distance_squared;
82  if (std::fabs(denom) <= eps)
83  {
84  // parallel lines
85  const coordT a0 = r10d0 / d0d0;
86  result = r0 + d0 * a0;
87  distance_squared = r10r10 - square(r10d0) / d0d0;
88  }
89  else
90  {
91  const coordT r10d1 = inner_product(r10, d1);
92  const coordT a0 = (-d1d1 * r10d0 + d0d1 * r10d1) / denom;
93  const coordT a1 = (-d0d1 * r10d0 + d0d0 * r10d1) / denom;
94 
95  result = ((r0 + d0 * a0) + (r1 + d1 * a1)) / 2;
96  distance_squared = (d1d1 * square(r10d0) - 2 * d0d1 * r10d0 * r10d1 + d0d0 * square(r10d1)) / denom + r10r10;
97  }
98  if (distance_squared >= 0)
99  return std::sqrt(distance_squared);
100  else
101  {
102  if (-distance_squared < eps)
103  return 0;
104  else
105  {
106  assert(false);
107  return std::sqrt(distance_squared); // will return NaN
108  }
109  }
110 }
111 
116 template <class coordT>
117 inline coordT
119 {
120  /* Rationale:
121  parametrise points on the lines as
122  r0 + a0 d0
123  where r0 is a point on the line, and d0 is its direction
124 
125  compute distance squared between those points, i.e.
126  norm(a0 d0 + r0-r1)^2
127  This can be written using inner products of all the vectors:
128  a0^2 d0d0 - 2 a0 r10d0 + r10r10
129  Minimising it (by computing derivatives) results in a linear eq in a0, a1:
130  a0 d0d0 == r10d0
131  */
132  const CartesianCoordinate3D<coordT>& r0 = line.p1();
133  const CartesianCoordinate3D<coordT> r10 = r1 - r0;
134  const CartesianCoordinate3D<coordT> d0 = line.p2() - line.p1();
135  const coordT d0d0 = inner_product(d0, d0);
136  const coordT r10d0 = inner_product(r10, d0);
137  const coordT r10r10 = inner_product(r10, r10);
138 
139  // const coordT a0 = r10d0/d0d0;
140  // result = r0 + d0*a0;
141  const coordT distance_squared = r10r10 - square(r10d0) / d0d0;
142  if (distance_squared >= 0)
143  return std::sqrt(distance_squared);
144  else
145  {
146  if (-distance_squared < d0d0 * 10E-5)
147  return 0;
148  else
149  {
150  assert(false);
151  return std::sqrt(distance_squared); // will return NaN
152  }
153  }
154 }
155 
161 template <class coordT>
162 inline void
166 {
167 
168  const CartesianCoordinate3D<coordT> difference = p2 - p1;
169 
170  const CartesianCoordinate3D<coordT> r10 = r1 - p1;
171 
172  float inner_prod = inner_product(difference, difference);
173 
174  const float u = inner_product(r10, difference) / inner_prod;
175 
176  r1.x() = p1.x() + u * difference.x();
177  r1.y() = p1.y() + u * difference.y();
178  r1.z() = p1.z() + u * difference.z();
179 }
180 
181 END_NAMESPACE_STIR
a templated class for 3-dimensional coordinates.
Definition: CartesianCoordinate3D.h:52
defines various classes for specifying a line in 3 dimensions
STL namespace.
void project_point_on_a_line(const CartesianCoordinate3D< coordT > &p1, const CartesianCoordinate3D< coordT > &p2, CartesianCoordinate3D< coordT > &r1)
Project a point on a line.
Definition: line_distances.h:163
A class for LORs.
Definition: CListRecord.h:39
coordT inner_product(const BasicCoordinate< num_dimensions, coordT > &p1, const BasicCoordinate< num_dimensions, coordT > &p2)
compute sum_i p1[i] * p2[i]
Definition: BasicCoordinate.inl:408
NUMBER square(const NUMBER &x)
returns the square of a number, templated.
Definition: common.h:146
coordT distance_between_line_and_point(const LORAs2Points< coordT > &line, const CartesianCoordinate3D< coordT > &r1)
find the distance between a point and a line
Definition: line_distances.h:118
defines the stir::CartesianCoordinate3D<coordT> class
coordT coordinate_between_2_lines(CartesianCoordinate3D< coordT > &result, const LORAs2Points< coordT > &line0, const LORAs2Points< coordT > &line1)
find a point half-way between 2 lines and their distance
Definition: line_distances.h:46