LSSTApplications  19.0.0-14-gb0260a2+9346bf5579,20.0.0+34a42eae2c,20.0.0+4d97b31663,20.0.0+5a87225079,20.0.0+8558dd3f48,20.0.0+9180b0bcc6,20.0.0+b290a576ab,20.0.0+b2ea66fa67,20.0.0+bba7c37fb9,20.0.0+cd847060a9,20.0.0+d138450326,20.0.0+d8493377e7,20.0.0+dcf29472a8,20.0.0+ef162136e0,20.0.0+f45b7d88f4,20.0.0-1-g10df615+6305e2b088,20.0.0-1-g253301a+dcf29472a8,20.0.0-1-g498fb60+ff88705a28,20.0.0-1-g4d801e7+d83096fe1b,20.0.0-1-g8a53f90+2817c06967,20.0.0-1-gc96f8cb+bba7c37fb9,20.0.0-1-gd1c87d7+2817c06967,20.0.0-1-gdb27ee5+52b05b0b7e,20.0.0-12-ga81c59a+61094d0bf4,20.0.0-18-g08fba245+aea2d85f7a,20.0.0-2-gec03fae+3bc057fb2a,20.0.0-28-gb33ccd1+1ae6d82017,20.0.0-3-gd2e950e+f45b7d88f4,20.0.0-3-gdd5c15c+990b4320db,20.0.0-4-g4a2362f+f45b7d88f4,20.0.0-5-gac0d578b1+6c871ee35c,20.0.0-5-gfcebe35+988ee452db,20.0.0-6-g01203fff+883dccf1c0,20.0.0-7-g3c4151b+a8ac49de8d,20.0.0-8-gc2abeef+bba7c37fb9,20.0.0-9-gabd0d4c+52b05b0b7e,w.2020.33
LSSTDataManagementBasePackage
orientation.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2014-2015 AURA/LSST.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <https://www.lsstcorp.org/LegalNotices/>.
21  */
22 
25 
27 
28 #include <algorithm>
29 
31 
32 
33 namespace lsst {
34 namespace sphgeom {
35 
36 namespace {
37 
38 // `BigFloat` is an exact floating point type.
39 struct BigFloat {
40  BigInteger * mantissa;
41  int exponent;
42 
43  BigFloat() : mantissa(0), exponent(0) {}
44  BigFloat(BigInteger * m) : mantissa(m), exponent(0) {}
45 };
46 
47 // `computeProduct` computes the product of 3 doubles exactly and stores the
48 // result in p.
49 void computeProduct(BigFloat & p, double d0, double d1, double d2) {
50  // This constant (2^53) is used to scale the fractions returned by
51  // std::frexp and turn them into integer mantissas.
52  static double const SCALE = 9007199254740992.0;
53  uint32_t buf[2];
54  BigInteger i(buf, sizeof(buf) / sizeof(uint32_t));
55  // Unpack the 3 input doubles into integer mantissas and exponents.
56  int e0 = 0;
57  int e1 = 0;
58  int e2 = 0;
59  double m0 = std::frexp(d0, &e0) * SCALE;
60  double m1 = std::frexp(d1, &e1) * SCALE;
61  double m2 = std::frexp(d2, &e2) * SCALE;
62  // Compute the product of the 3 input doubles using exact arithmetic.
63  p.mantissa->setTo(static_cast<int64_t>(m0));
64  i.setTo(static_cast<int64_t>(m1));
65  p.mantissa->multiply(i);
66  i.setTo(static_cast<int64_t>(m2));
67  p.mantissa->multiply(i);
68  // Finally, adjust the exponent of the result to compensate for the 3
69  // multiplications by 2^53 performed above.
70  p.exponent = e0 + e1 + e2 - 3 * 53;
71 }
72 
73 } // unnamed namespace
74 
75 
77  Vector3d const & b,
78  Vector3d const & c)
79 {
80  // Product mantissa storage buffers.
81  uint32_t mantissaBuffers[6][6];
82  // Product mantissas.
83  BigInteger mantissas[6] = {
84  BigInteger(mantissaBuffers[0], 6),
85  BigInteger(mantissaBuffers[1], 6),
86  BigInteger(mantissaBuffers[2], 6),
87  BigInteger(mantissaBuffers[3], 6),
88  BigInteger(mantissaBuffers[4], 6),
89  BigInteger(mantissaBuffers[5], 6)
90  };
91  BigFloat products[6] = {
92  BigFloat(&mantissas[0]),
93  BigFloat(&mantissas[1]),
94  BigFloat(&mantissas[2]),
95  BigFloat(&mantissas[3]),
96  BigFloat(&mantissas[4]),
97  BigFloat(&mantissas[5])
98  };
99  // An accumulator and its storage.
100  uint32_t accumulatorBuffer[512];
101  BigInteger accumulator(accumulatorBuffer,
102  sizeof(accumulatorBuffer) / sizeof(uint32_t));
103  // Compute the products in the determinant. Performing all multiplication
104  // up front means that each product mantissa occupies at most 3*53 bits.
105  computeProduct(products[0], a.x(), b.y(), c.z());
106  computeProduct(products[1], a.x(), b.z(), c.y());
107  computeProduct(products[2], a.y(), b.z(), c.x());
108  computeProduct(products[3], a.y(), b.x(), c.z());
109  computeProduct(products[4], a.z(), b.x(), c.y());
110  computeProduct(products[5], a.z(), b.y(), c.x());
111  mantissas[1].negate();
112  mantissas[3].negate();
113  mantissas[5].negate();
114  // Sort the array of products in descending exponent order.
115  std::sort(products, products + 6, [](BigFloat const & a, BigFloat const & b) {
116  return a.exponent > b.exponent;
117  });
118  // First, initialize the accumulator to the product with the highest
119  // exponent, then add the remaining products. Prior to each addition, we
120  // must shift the accumulated value so that its radix point lines up with
121  // the the radix point of the product to add.
122  //
123  // More precisely, at each step we have an accumulated value A·2ʲ and a
124  // product P·2ᵏ, and we update the accumulator to equal (A·2ʲ⁻ᵏ + P)·2ᵏ.
125  // Because the products were sorted beforehand, j ≥ k and 2ʲ⁻ᵏ is an
126  // integer.
127  accumulator = *products[0].mantissa;
128  for (int i = 1; i < 6; ++i) {
129  accumulator.multiplyPow2(products[i - 1].exponent - products[i].exponent);
130  accumulator.add(*products[i].mantissa);
131  }
132  return accumulator.getSign();
133 }
134 
136  UnitVector3d const & b,
137  UnitVector3d const & c)
138 {
139  // This constant is a little more than 5ε, where ε = 2^-53. When multiplied
140  // by the permanent of |M|, it gives an error bound on the determinant of
141  // M. Here, M is a 3x3 matrix and |M| denotes the matrix obtained by
142  // taking the absolute value of each of its components. The derivation of
143  // this proceeds in the same manner as the derivation of the error bounds
144  // in section 4.3 of:
145  //
146  // Adaptive Precision Floating-Point Arithmetic
147  // and Fast Robust Geometric Predicates,
148  // Jonathan Richard Shewchuk,
149  // Discrete & Computational Geometry 18(3):305–363, October 1997.
150  //
151  // available online at http://www.cs.berkeley.edu/~jrs/papers/robustr.pdf
152  static double const relativeError = 5.6e-16;
153  // Because all 3 unit vectors are normalized, the maximum absolute value of
154  // any vector component, cross product component or dot product term in
155  // the calculation is very close to 1. The permanent of |M| must therefore
156  // be below 3 + c, where c is some small multiple of ε. This constant, a
157  // little larger than 3 * 5ε, is an upper bound on the absolute error in
158  // the determinant calculation.
159  static double const maxAbsoluteError = 1.7e-15;
160  // This constant accounts for floating point underflow (assuming hardware
161  // without gradual underflow, just to be conservative) in the computation
162  // of det(M). It is a little more than 14 * 2^-1022.
163  static double const minAbsoluteError = 4.0e-307;
164 
165  double bycz = b.y() * c.z();
166  double bzcy = b.z() * c.y();
167  double bzcx = b.z() * c.x();
168  double bxcz = b.x() * c.z();
169  double bxcy = b.x() * c.y();
170  double bycx = b.y() * c.x();
171  double determinant = a.x() * (bycz - bzcy) +
172  a.y() * (bzcx - bxcz) +
173  a.z() * (bxcy - bycx);
174  if (determinant > maxAbsoluteError) {
175  return 1;
176  } else if (determinant < -maxAbsoluteError) {
177  return -1;
178  }
179  // Expend some more effort on what is hopefully a tighter error bound
180  // before falling back on arbitrary precision arithmetic.
181  double permanent = std::fabs(a.x()) * (std::fabs(bycz) + std::fabs(bzcy)) +
182  std::fabs(a.y()) * (std::fabs(bzcx) + std::fabs(bxcz)) +
183  std::fabs(a.z()) * (std::fabs(bxcy) + std::fabs(bycx));
184  double maxError = relativeError * permanent + minAbsoluteError;
185  if (determinant > maxError) {
186  return 1;
187  } else if (determinant < -maxError) {
188  return -1;
189  }
190  // Avoid the slow path when any two inputs are identical or antipodal.
191  if (a == b || b == c || a == c || a == -b || b == -c || a == -c) {
192  return 0;
193  }
194  return orientationExact(a, b, c);
195 }
196 
197 
198 namespace {
199 
200  inline int _orientationXYZ(double ab, double ba) {
201  // Calling orientation() with a first argument of (1,0,0), (0,1,0) or
202  // (0,0,1) corresponds to computing the sign of a 2x2 determinant
203  // rather than a 3x3 determinant. The corresponding error bounds
204  // are also tighter.
205  static double const relativeError = 1.12e-16; // > 2^-53
206  static double const maxAbsoluteError = 1.12e-16; // > 2^-53
207  static double const minAbsoluteError = 1.0e-307; // > 3 * 2^-1022
208 
209  double determinant = ab - ba;
210  if (determinant > maxAbsoluteError) {
211  return 1;
212  } else if (determinant < -maxAbsoluteError) {
213  return -1;
214  }
215  double permanent = std::fabs(ab) + std::fabs(ba);
216  double maxError = relativeError * permanent + minAbsoluteError;
217  if (determinant > maxError) {
218  return 1;
219  } else if (determinant < -maxError) {
220  return -1;
221  }
222  return 0;
223  }
224 
225 }
226 
227 int orientationX(UnitVector3d const & b, UnitVector3d const & c) {
228  int o = _orientationXYZ(b.y() * c.z(), b.z() * c.y());
229  return (o != 0) ? o : orientationExact(UnitVector3d::X(), b, c);
230 }
231 
232 int orientationY(UnitVector3d const & b, UnitVector3d const & c) {
233  int o = _orientationXYZ(b.z() * c.x(), b.x() * c.z());
234  return (o != 0) ? o : orientationExact(UnitVector3d::Y(), b, c);
235 }
236 
237 int orientationZ(UnitVector3d const & b, UnitVector3d const & c) {
238  int o = _orientationXYZ(b.x() * c.y(), b.y() * c.x());
239  return (o != 0) ? o : orientationExact(UnitVector3d::Z(), b, c);
240 }
241 
242 }} // namespace lsst::sphgeom
lsst::sphgeom::BigInteger::negate
void negate()
negate multiplies this integer by -1.
Definition: BigInteger.h:103
std::fabs
T fabs(T... args)
orientation.h
This file declares functions for orienting points on the sphere.
lsst::sphgeom::UnitVector3d::x
double x() const
Definition: UnitVector3d.h:145
lsst::sphgeom::orientation
int orientation(UnitVector3d const &a, UnitVector3d const &b, UnitVector3d const &c)
orientation computes and returns the orientations of 3 unit vectors a, b and c.
Definition: orientation.cc:135
BigInteger.h
This file declares a class for arbitrary precision integers.
std::sort
T sort(T... args)
lsst::sphgeom::Vector3d
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
lsst::sphgeom::UnitVector3d::Y
static UnitVector3d Y()
Definition: UnitVector3d.h:97
std::frexp
T frexp(T... args)
mantissa
BigInteger * mantissa
Definition: orientation.cc:40
lsst::sphgeom::UnitVector3d
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
lsst::sphgeom::Vector3d::z
double z() const
Definition: Vector3d.h:70
lsst::sphgeom::BigInteger
BigInteger is an arbitrary precision signed integer class.
Definition: BigInteger.h:45
lsst::sphgeom::orientationX
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
Definition: orientation.cc:227
lsst::sphgeom::UnitVector3d::Z
static UnitVector3d Z()
Definition: UnitVector3d.h:101
b
table::Key< int > b
Definition: TransmissionCurve.cc:467
lsst::sphgeom::UnitVector3d::z
double z() const
Definition: UnitVector3d.h:149
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::sphgeom::orientationY
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
Definition: orientation.cc:232
lsst::sphgeom::UnitVector3d::y
double y() const
Definition: UnitVector3d.h:147
lsst::sphgeom::Vector3d::y
double y() const
Definition: Vector3d.h:68
a
table::Key< int > a
Definition: TransmissionCurve.cc:466
lsst::sphgeom::orientationExact
int orientationExact(Vector3d const &a, Vector3d const &b, Vector3d const &c)
orientationExact computes and returns the orientations of 3 vectors a, b and c, which need not be nor...
Definition: orientation.cc:76
lsst::sphgeom::UnitVector3d::X
static UnitVector3d X()
Definition: UnitVector3d.h:93
exponent
int exponent
Definition: orientation.cc:41
lsst::sphgeom::orientationZ
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
Definition: orientation.cc:237
lsst::sphgeom::Vector3d::x
double x() const
Definition: Vector3d.h:66
m
int m
Definition: SpanSet.cc:49