Point2.cpp
Go to the documentation of this file.
1 /* ----------------------------------------------------------------------------
2 
3  * GTSAM Copyright 2010, Georgia Tech Research Corporation,
4  * Atlanta, Georgia 30332-0415
5  * All Rights Reserved
6  * Authors: Frank Dellaert, et al. (see THANKS for the full author list)
7 
8  * See LICENSE for the license information
9 
10  * -------------------------------------------------------------------------- */
11 
18 #include <gtsam/geometry/Point2.h>
19 #include <cmath>
20 #include <iostream>
21 
22 using namespace std;
23 
24 namespace gtsam {
25 
26 /* ************************************************************************* */
28  double r = std::sqrt(p.x() * p.x() + p.y() * p.y());
29  if (H) {
30  if (std::abs(r) > 1e-10)
31  *H << p.x() / r, p.y() / r;
32  else
33  *H << 1, 1; // really infinity, why 1 ?
34  }
35  return r;
36 }
37 
38 /* ************************************************************************* */
39 double distance2(const Point2& p, const Point2& q, OptionalJacobian<1, 2> H1,
41  Point2 d = q - p;
42  if (H1 || H2) {
43  Matrix12 H;
44  double r = norm2(d, H);
45  if (H1) *H1 = -H;
46  if (H2) *H2 = H;
47  return r;
48  } else {
49  return d.norm();
50  }
51 }
52 
53 /* ************************************************************************* */
54 // Math inspired by http://paulbourke.net/geometry/circlesphere/
55 std::optional<Point2> circleCircleIntersection(double R_d, double r_d,
56  double tol) {
57 
58  double R2_d2 = R_d*R_d; // Yes, RD-D2 !
59  double f = 0.5 + 0.5*(R2_d2 - r_d*r_d);
60  double h2 = R2_d2 - f*f; // just right triangle rule
61 
62  // h^2<0 is equivalent to (d > (R + r) || d < (R - r))
63  // Hence, there are only solutions if >=0
64  if (h2<-tol) return {}; // allow *slightly* negative
65  else if (h2<tol) return Point2(f,0.0); // one solution
66  else return Point2(f,std::sqrt(h2)); // two solutions
67 }
68 
69 /* ************************************************************************* */
71  std::optional<Point2> fh) {
72 
73  list<Point2> solutions;
74  // If fh==std::nullopt, there are no solutions, i.e., d > (R + r) || d < (R - r)
75  if (fh) {
76  // vector between circle centers
77  Point2 c12 = c2-c1;
78 
79  // Determine p2, the point where the line through the circle
80  // intersection points crosses the line between the circle centers.
81  Point2 p2 = c1 + fh->x() * c12;
82 
83  // If h == 0, the circles are touching, so just return one point
84  if (fh->y()==0.0)
85  solutions.push_back(p2);
86  else {
87  // determine the offsets of the intersection points from p
88  Point2 offset = fh->y() * Point2(-c12.y(), c12.x());
89 
90  // Determine the absolute intersection points.
91  solutions.push_back(p2 + offset);
92  solutions.push_back(p2 - offset);
93  }
94  }
95  return solutions;
96 }
97 
98 /* ************************************************************************* */
99 list<Point2> circleCircleIntersection(Point2 c1, double r1, Point2 c2,
100  double r2, double tol) {
101 
102  // distance between circle centers.
103  double d = distance2(c1, c2);
104 
105  // centers coincide, either no solution or infinite number of solutions.
106  if (d<1e-9) return list<Point2>();
107 
108  // Calculate f and h given normalized radii
109  double _d = 1.0/d, R_d = r1*_d, r_d=r2*_d;
110  std::optional<Point2> fh = circleCircleIntersection(R_d,r_d);
111 
112  // Call version that takes fh
113  return circleCircleIntersection(c1, c2, fh);
114 }
115 
116 Point2Pair means(const std::vector<Point2Pair> &abPointPairs) {
117  const size_t n = abPointPairs.size();
118  if (n == 0) throw std::invalid_argument("Point2::mean input Point2Pair vector is empty");
119  Point2 aSum(0, 0), bSum(0, 0);
120  for (const Point2Pair &abPair : abPointPairs) {
121  aSum += abPair.first;
122  bSum += abPair.second;
123  }
124  const double f = 1.0 / n;
125  return {aSum * f, bSum * f};
126 }
127 
128 /* ************************************************************************* */
129 ostream &operator<<(ostream &os, const gtsam::Point2Pair &p) {
130  os << p.first << " <-> " << p.second;
131  return os;
132 }
133 
134 } // namespace gtsam
Vector2 Point2
Definition: Point2.h:32
int n
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate offset
Definition: BFloat16.h:88
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set set set surface set nocontour set clabel set mapping cartesian set nohidden3d set cntrparam order set cntrparam linear set cntrparam levels auto set cntrparam points set size set set xzeroaxis lt lw set x2zeroaxis lt lw set yzeroaxis lt lw set y2zeroaxis lt lw set tics in set ticslevel set tics set mxtics default set mytics default set mx2tics default set my2tics default set xtics border mirror norotate autofreq set ytics border mirror norotate autofreq set ztics border nomirror norotate autofreq set nox2tics set noy2tics set timestamp bottom norotate set rrange [*:*] noreverse nowriteback set trange [*:*] noreverse nowriteback set urange [*:*] noreverse nowriteback set vrange [*:*] noreverse nowriteback set xlabel matrix size set x2label set timefmt d m y n H
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
static const double r2
EIGEN_DEVICE_FUNC const Scalar & q
traits
Definition: chartTesting.h:28
Point2Pair means(const std::vector< Point2Pair > &abPointPairs)
Calculate the two means of a set of Point2 pairs.
Definition: Point2.cpp:116
ostream & operator<<(ostream &os, const gtsam::Point2Pair &p)
Definition: Point2.cpp:129
static const double r1
ofstream os("timeSchurFactors.csv")
double distance2(const Point2 &p, const Point2 &q, OptionalJacobian< 1, 2 > H1, OptionalJacobian< 1, 2 > H2)
distance between two points
Definition: Point2.cpp:39
std::pair< Point2, Point2 > Point2Pair
Definition: Point2.h:35
list< Point2 > circleCircleIntersection(Point2 c1, double r1, Point2 c2, double r2, double tol)
Intersect 2 circles.
Definition: Point2.cpp:99
float * p
static Point3 p2
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const G double tol
Definition: Group.h:86
#define abs(x)
Definition: datatypes.h:17
2D Point
double norm2(const Point2 &p, OptionalJacobian< 1, 2 > H)
Distance of the point from the origin, with Jacobian.
Definition: Point2.cpp:27


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:14