GteApprGreatCircle3.h
Go to the documentation of this file.
1 // David Eberly, Geometric Tools, Redmond WA 98052
2 // Copyright (c) 1998-2017
3 // Distributed under the Boost Software License, Version 1.0.
4 // http://www.boost.org/LICENSE_1_0.txt
5 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
6 // File Version: 3.0.0 (2016/06/19)
7 
8 #pragma once
9 
13 
14 namespace gte
15 {
16 
17 // Least-squares fit of a great circle to unit-length vectors (x,y,z) by
18 // using distance measurements orthogonal (and measured along great circles)
19 // to the proposed great circle. The inputs akPoint[] are unit length. The
20 // returned value is unit length, call it N. The fitted great circle is
21 // defined by Dot(N,X) = 0, where X is a unit-length vector on the great
22 // circle.
23 template <typename Real>
25 {
26 public:
27  void operator()(int numPoints, Vector3<Real> const* points,
28  Vector3<Real>& normal);
29 };
30 
31 
32 // In addition to the least-squares fit of a great circle, the input vectors
33 // are projected onto that circle. The sector of smallest angle (possibly
34 // obtuse) that contains the points is computed. The endpoints of the arc
35 // of the sector are returned. The returned endpoints A0 and A1 are
36 // perpendicular to the returned normal N. Moreover, when you view the
37 // arc by looking at the plane of the great circle with a view direction
38 // of -N, the arc is traversed counterclockwise starting at A0 and ending
39 // at A1.
40 template <typename Real>
42 {
43 public:
44  void operator()(int numPoints, Vector3<Real> const* points,
45  Vector3<Real>& normal, Vector3<Real>& arcEnd0,
46  Vector3<Real>& arcEnd1);
47 };
48 
49 
50 template <typename Real>
52  Vector3<Real> const* points, Vector3<Real>& normal)
53 {
54  // Compute the covariance matrix of the vectors.
55  Real covar00 = (Real)0, covar01 = (Real)0, covar02 = (Real)0;
56  Real covar11 = (Real)0, covar12 = (Real)0, covar22 = (Real)0;
57  for (int i = 0; i < numPoints; i++)
58  {
59  Vector3<Real> diff = points[i];
60  covar00 += diff[0] * diff[0];
61  covar01 += diff[0] * diff[1];
62  covar02 += diff[0] * diff[2];
63  covar11 += diff[1] * diff[1];
64  covar12 += diff[1] * diff[2];
65  covar22 += diff[2] * diff[2];
66  }
67 
68  Real invNumPoints = ((Real)1) / static_cast<Real>(numPoints);
69  covar00 *= invNumPoints;
70  covar01 *= invNumPoints;
71  covar02 *= invNumPoints;
72  covar11 *= invNumPoints;
73  covar12 *= invNumPoints;
74  covar22 *= invNumPoints;
75 
76  // Solve the eigensystem.
78  std::array<Real, 3> eval;
79  std::array<std::array<Real, 3>, 3> evec;
80  es(covar00, covar01, covar02, covar11, covar12, covar22, false, +1,
81  eval, evec);
82  normal = evec[0];
83 }
84 
85 template <typename Real>
87  Vector3<Real> const* points, Vector3<Real>& normal,
88  Vector3<Real>& arcEnd0, Vector3<Real>& arcEnd1)
89 {
90  // Get the least-squares great circle for the vectors. The circle is on
91  // the plane Dot(N,X) = 0. Generate a basis from N.
92  Vector3<Real> basis[3]; // { N, U, V }
93  ApprGreatCircle3<Real>()(numPoints, points, basis[0]);
95 
96  // The vectors are X[i] = u[i]*U + v[i]*V + w[i]*N. The projections
97  // are P[i] = (u[i]*U + v[i]*V)/sqrt(u[i]*u[i] + v[i]*v[i]). The great
98  // circle is parameterized by C(t) = cos(t)*U + sin(t)*V. Compute the
99  // angles t in [-pi,pi] for the projections onto the great circle. It
100  // is not necesarily to normalize (u[i],v[i]), instead computing
101  // t = atan2(v[i],u[i]).
102  std::vector<std::array<Real, 3>> items(numPoints); // item (u, v, angle)
103  for (int i = 0; i < numPoints; ++i)
104  {
105  items[i][0] = Dot(basis[1], points[i]);
106  items[i][1] = Dot(basis[2], points[i]);
107  items[i][2] = atan2(items[i][1], items[i][0]);
108  }
109  std::sort(items.begin(), items.end(),
110  [](std::array<Real, 3> const& item0, std::array<Real, 3> const& item1)
111  {
112  return item0[2] < item1[2];
113  }
114  );
115 
116  // Locate the pair of consecutive angles whose difference is a maximum.
117  // Effectively, we are constructing a cone of minimum angle that contains
118  // the unit-length vectors.
119  int numPointsM1 = numPoints - 1;
120  Real maxDiff = (Real)GTE_C_TWO_PI + items[0][2] - items[numPointsM1][2];
121  int end0 = 0, end1 = numPointsM1;
122  for (int i0 = 0, i1 = 1; i0 < numPointsM1; i0 = i1++)
123  {
124  Real diff = items[i1][2] - items[i0][2];
125  if (diff > maxDiff)
126  {
127  maxDiff = diff;
128  end0 = i1;
129  end1 = i0;
130  }
131  }
132 
133  normal = basis[0];
134  arcEnd0 = items[end0][0] * basis[1] + items[end0][1] * basis[2];
135  arcEnd1 = items[end1][0] * basis[1] + items[end1][1] * basis[2];
136  Normalize(arcEnd0);
137  Normalize(arcEnd1);
138 }
139 
140 
141 }
void operator()(int numPoints, Vector3< Real > const *points, Vector3< Real > &normal)
GLfixed GLfixed GLint GLint GLfixed points
Definition: glext.h:4927
DualQuaternion< Real > Dot(DualQuaternion< Real > const &d0, DualQuaternion< Real > const &d1)
Real Normalize(GVector< Real > &v, bool robust=false)
Definition: GteGVector.h:454
Real ComputeOrthogonalComplement(int numInputs, Vector2< Real > *v, bool robust=false)
Definition: GteVector2.h:123
#define GTE_C_TWO_PI
Definition: GteConstants.h:20
void operator()(int numPoints, Vector3< Real > const *points, Vector3< Real > &normal, Vector3< Real > &arcEnd0, Vector3< Real > &arcEnd1)


geometric_tools_engine
Author(s): Yijiang Huang
autogenerated on Thu Jul 18 2019 03:59:59