LSST Applications  21.0.0+3c14b91618,21.0.0+9f51b1e3f7,21.0.0-1-ga51b5d4+6691386486,21.0.0-10-g2408eff+1a328412bf,21.0.0-10-g560fb7b+52fd22d7b4,21.0.0-10-g8d1d15d+2f9043cae0,21.0.0-10-gcf60f90+d15de71c48,21.0.0-11-g25eff31+d43066e4ef,21.0.0-15-g490e301a+a676f0d5cf,21.0.0-2-g103fe59+4758c8ef83,21.0.0-2-g1367e85+a9f57e981a,21.0.0-2-g45278ab+9f51b1e3f7,21.0.0-2-g5242d73+a9f57e981a,21.0.0-2-g7f82c8f+1bcc828e4f,21.0.0-2-g8f08a60+e6fd6d9ff9,21.0.0-2-ga326454+1bcc828e4f,21.0.0-2-gde069b7+66c51b65da,21.0.0-2-gecfae73+251b9830c3,21.0.0-2-gfc62afb+a9f57e981a,21.0.0-20-g09baf175d+e1e7d1c708,21.0.0-3-g357aad2+854c3902c3,21.0.0-3-g4be5c26+a9f57e981a,21.0.0-3-g65f322c+feaa1990e9,21.0.0-3-g7d9da8d+3c14b91618,21.0.0-3-ge02ed75+bda8df9b93,21.0.0-4-g591bb35+bda8df9b93,21.0.0-4-g65b4814+52fd22d7b4,21.0.0-4-g88306b8+656365ce3f,21.0.0-4-gccdca77+86bf7a300d,21.0.0-4-ge8a399c+b99e86088e,21.0.0-5-gd00fb1e+6a0dc09319,21.0.0-51-gd3b42663+bf9f0d1c8b,21.0.0-6-g2d4f3f3+9f51b1e3f7,21.0.0-7-g04766d7+ee2ae02087,21.0.0-7-g98eecf7+3609eddee2,21.0.0-7-gde99fe0+bda8df9b93,21.0.0-8-gd70ce6f+79e45e76e4,master-gac4afde19b+bda8df9b93,w.2021.10
LSST Data Management Base Package
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
This file declares a class for arbitrary precision integers.
int m
Definition: SpanSet.cc:49
BigInteger is an arbitrary precision signed integer class.
Definition: BigInteger.h:45
BigInteger & multiplyPow2(unsigned n)
multiplyPow2 multiplies this integer by 2ⁿ.
Definition: BigInteger.cc:214
int getSign() const
getSign returns -1, 0 or 1 if this integer is negative, zero or positive.
Definition: BigInteger.h:68
void negate()
negate multiplies this integer by -1.
Definition: BigInteger.h:103
BigInteger & add(BigInteger const &b)
add adds b to this integer.
Definition: BigInteger.cc:156
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
static UnitVector3d Z()
Definition: UnitVector3d.h:101
static UnitVector3d Y()
Definition: UnitVector3d.h:97
static UnitVector3d X()
Definition: UnitVector3d.h:93
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
double x() const
Definition: Vector3d.h:66
double y() const
Definition: Vector3d.h:68
double z() const
Definition: Vector3d.h:70
T fabs(T... args)
T frexp(T... args)
int orientationZ(UnitVector3d const &b, UnitVector3d const &c)
orientationZ(b, c) is equivalent to orientation(UnitVector3d::Z(), b, c).
Definition: orientation.cc:237
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
int orientationX(UnitVector3d const &b, UnitVector3d const &c)
orientationX(b, c) is equivalent to orientation(UnitVector3d::X(), b, c).
Definition: orientation.cc:227
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
int orientationY(UnitVector3d const &b, UnitVector3d const &c)
orientationY(b, c) is equivalent to orientation(UnitVector3d::Y(), b, c).
Definition: orientation.cc:232
A base class for image defects.
T sort(T... args)
This file declares functions for orienting points on the sphere.
BigInteger * mantissa
Definition: orientation.cc:40
int exponent
Definition: orientation.cc:41
table::Key< int > b
table::Key< int > a