LSSTApplications  19.0.0-14-gb0260a2+f7b33f63ed,20.0.0+028ca7aeb5,20.0.0+578ca8d785,20.0.0+95fe209e70,20.0.0+d1f8484762,20.0.0+ec0eb203f4,20.0.0-1-g253301a+95fe209e70,20.0.0-1-g2b7511a+ce6f728470,20.0.0-1-g5b95a8c+36cabaa438,20.0.0-1-gedffbd8+08e3982e2f,20.0.0-11-gd9dafd18+ee11e5b92d,20.0.0-16-gfab17e72e+f119078cbc,20.0.0-19-gb3c8180+76b8b4580d,20.0.0-2-g4dae9ad+b9c65d2db1,20.0.0-2-g61b8584+62d102034e,20.0.0-2-gb780d76+f45b7d88f4,20.0.0-2-ged6426c+c16e0e6c56,20.0.0-2-gf072044+95fe209e70,20.0.0-2-gf1f7952+b9c65d2db1,20.0.0-22-gdf434b7+b9c65d2db1,20.0.0-23-g10eeb28+266feae6ef,20.0.0-25-g3dcad98+7edba9c73c,20.0.0-3-g1653f94+62d102034e,20.0.0-3-g4cc78c6+1c0c6df15e,20.0.0-3-g8f21e14+cae0a2a29f,20.0.0-3-gbd60e8c+eb59cfae10,20.0.0-3-gbecbe05+a2681a7a86,20.0.0-36-g96d5af2b+330d9ff9b7,20.0.0-4-g97dc21a+7edba9c73c,20.0.0-4-gb4befbc+0084b8bcb8,20.0.0-4-gfea843c+f45b7d88f4,20.0.0-5-gdfe0fee+88049eac2d,20.0.0-6-g64f541c+f45b7d88f4,20.0.0-6-g9a5b7a1+db91213b54,20.0.0-67-g32d6278+0d86c2e7cd,20.0.0-9-g4aef684+e9f4309ef0,w.2020.44
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)
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
mantissa
BigInteger * mantissa
Definition: orientation.cc:40
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
exponent
int exponent
Definition: orientation.cc:41
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
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