LSSTApplications  19.0.0+10,19.0.0+80,19.0.0-1-g20d9b18+31,19.0.0-1-g49a97f9+4,19.0.0-1-g8c57eb9+31,19.0.0-1-g9a028c0+13,19.0.0-1-ga72da6b+3,19.0.0-1-gb77924a+15,19.0.0-1-gbfe0924+66,19.0.0-1-gc3516c3,19.0.0-1-gefe1d0d+49,19.0.0-1-gf8cb8b4+3,19.0.0-14-g7511ce4+6,19.0.0-16-g3dc1a33c+6,19.0.0-17-g59f0e8a+4,19.0.0-17-g9c22e3c+9,19.0.0-18-g35bb99870+2,19.0.0-19-g2772d4a+9,19.0.0-2-g260436e+53,19.0.0-2-g31cdcee,19.0.0-2-g9675b69+3,19.0.0-2-gaa2795f,19.0.0-2-gbcc4de1,19.0.0-2-gd6f004e+6,19.0.0-2-gde8e5e3+5,19.0.0-2-gff6972b+15,19.0.0-21-ge306cd8,19.0.0-21-gf122e69+2,19.0.0-3-g2d43a51+2,19.0.0-3-gf3b1435+6,19.0.0-4-g56feb96,19.0.0-4-gac56cce+17,19.0.0-49-gce872c1+1,19.0.0-5-g66946eb+6,19.0.0-5-gd8897ba+6,19.0.0-51-gfc4a647e,19.0.0-7-g686a884+5,19.0.0-7-g886f805+5,19.0.0-8-g305ff64,w.2020.17
LSSTDataManagementBasePackage
utils.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 
26 #include "lsst/sphgeom/utils.h"
27 
28 #include <cmath>
29 
31 
32 
33 namespace lsst {
34 namespace sphgeom {
35 
37  Vector3d const & a,
38  Vector3d const & b,
39  Vector3d const & n)
40 {
41  Vector3d vxn = v.cross(n);
42  if (vxn.dot(a) > 0.0 && vxn.dot(b) < 0.0) {
43  // v is in the lune defined by the half great circle passing through
44  // n and a and the half great circle passing through n and b, so p
45  // is in the interior of the great circle segment from a to b. The
46  // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
47  // and ‖v‖ ‖n‖ cos θ = ‖v × n‖. The desired squared chord length is
48  // 4 sin²(θ/2).
49  double s = std::fabs(v.dot(n));
50  double c = vxn.getNorm();
51  double theta = (c == 0.0) ? 0.5 * PI : std::atan(s / c);
52  double d = std::sin(0.5 * theta);
53  return 4.0 * d * d;
54  }
55  return 4.0;
56 }
57 
59  Vector3d const & a,
60  Vector3d const & b,
61  Vector3d const & n)
62 {
63  Vector3d vxn = v.cross(n);
64  if (vxn.dot(a) < 0.0 && vxn.dot(b) > 0.0) {
65  // v is in the lune defined by the half great circle passing through
66  // n and -a and the half great circle passing through n and -b, so p
67  // is in the interior of the great circle segment from a to b. The
68  // angle θ between p and v satisfies ‖v‖ ‖n‖ sin θ = |v·n|,
69  // and ‖v‖ ‖n‖ cos θ = -‖v × n‖. The desired squared chord length is
70  // 4 sin²(θ/2).
71  double s = std::fabs(v.dot(n));
72  double c = - vxn.getNorm();
73  double d = std::sin(0.5 * std::atan2(s, c));
74  return 4.0 * d * d;
75  }
76  return 0.0;
77 }
78 
80  UnitVector3d const & v1,
81  UnitVector3d const & v2)
82 {
83  // For the details, see:
84  //
85  // The centroid and inertia tensor for a spherical triangle
86  // John E. Brock
87  // 1974, Naval Postgraduate School, Monterey Calif.
88  //
89  // https://openlibrary.org/books/OL25493734M/The_centroid_and_inertia_tensor_for_a_spherical_triangle
90 
91  Vector3d x01 = v0.robustCross(v1); // twice the cross product of v0 and v1
92  Vector3d x12 = v1.robustCross(v2);
93  Vector3d x20 = v2.robustCross(v0);
94  double s01 = 0.5 * x01.normalize(); // sine of the angle between v0 and v1
95  double s12 = 0.5 * x12.normalize();
96  double s20 = 0.5 * x20.normalize();
97  double c01 = v0.dot(v1); // cosine of the angle between v0 and v1
98  double c12 = v1.dot(v2);
99  double c20 = v2.dot(v0);
100  double a0 = (s12 == 0.0 && c12 == 0.0) ? 0.0 : std::atan2(s12, c12);
101  double a1 = (s20 == 0.0 && c20 == 0.0) ? 0.0 : std::atan2(s20, c20);
102  double a2 = (s01 == 0.0 && c01 == 0.0) ? 0.0 : std::atan2(s01, c01);
103  return 0.5 * (x01 * a2 + x12 * a0 + x20 * a1);
104 }
105 
106 }} // namespace lsst::sphgeom
double normalize()
normalize scales this vector to have unit norm and returns its norm prior to scaling.
Definition: Vector3d.cc:41
T atan(T... args)
This file declares miscellaneous utility functions.
table::Key< int > b
table::Key< int > a
Vector3d robustCross(UnitVector3d const &v) const
a.robustCross(b) is (b + a).cross(b - a) - twice the cross product of a and b.
Definition: UnitVector3d.h:161
Vector3d is a vector in ℝ³ with components stored in double precision.
Definition: Vector3d.h:44
T atan2(T... args)
double getMinSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector closest to v that lies on the plane with normal n in the direction of the cr...
Definition: utils.cc:36
T sin(T... args)
This file declares a class for representing unit vectors in ℝ³.
Vector3d getWeightedCentroid(UnitVector3d const &v0, UnitVector3d const &v1, UnitVector3d const &v2)
getWeightedCentroid returns the center of mass of the given spherical triangle (assuming a uniform ma...
Definition: utils.cc:79
A base class for image defects.
constexpr double PI
Definition: constants.h:36
double getMaxSquaredChordLength(Vector3d const &v, Vector3d const &a, Vector3d const &b, Vector3d const &n)
Let p be the unit vector furthest from v that lies on the plane with normal n in the direction of the...
Definition: utils.cc:58
T fabs(T... args)
Vector3d cross(Vector3d const &v) const
cross returns the cross product of this vector and v.
Definition: Vector3d.h:101
UnitVector3d is a unit vector in ℝ³ with components stored in double precision.
Definition: UnitVector3d.h:55
double dot(Vector3d const &v) const
dot returns the inner product of this unit vector and v.
Definition: UnitVector3d.h:152
double dot(Vector3d const &v) const
dot returns the inner product of this vector and v.
Definition: Vector3d.h:73
double getNorm() const
getNorm returns the L2 norm of this vector.
Definition: Vector3d.h:81