LSSTApplications  16.0-1-gce273f5+15,16.0-1-gf77f410+8,16.0-10-g230e10e+8,16.0-10-g5a5abec+4,16.0-10-gc1446dd+8,16.0-11-g5cf5345,16.0-12-g1dc09ba+2,16.0-12-g726f8f3+6,16.0-13-g4c33ca5+8,16.0-13-g5f6c0b1+8,16.0-13-gd9b1b71+8,16.0-14-g71e547a+4,16.0-17-g0bdc215,16.0-17-g6a7bfb3b+8,16.0-18-g95848a16+2,16.0-2-g0febb12+12,16.0-2-g839ba83+46,16.0-2-g9d5294e+35,16.0-3-g404ea43+7,16.0-3-gbc759ec+6,16.0-3-gcfd6c53+33,16.0-4-g03cf288+24,16.0-4-g13a27c5+10,16.0-4-g2cc461c+3,16.0-4-g5f3a788+11,16.0-4-g8a0f11a+30,16.0-4-ga3eb747+1,16.0-43-gaa65a4573+4,16.0-5-g1991253+8,16.0-5-g865efd9+8,16.0-5-gb3f8a4b+40,16.0-5-gd0f1235+4,16.0-6-g316b399+8,16.0-6-gf0acd13+27,16.0-6-gf9cb114+9,16.0-7-ga8e1655+4,16.0-8-g23bbf3f+1,16.0-8-g4dec96c+21,16.0-8-gfd407c0,w.2018.39
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)
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.
Definition: cameraGeom.dox:3
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)
solver_t * s
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
This file declares miscellaneous utility functions.
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