OGS
AnalyticalGeometry-impl.h
Go to the documentation of this file.
1 
11 namespace GeoLib
12 {
13 template <typename InputIterator>
14 std::pair<Eigen::Vector3d, double> getNewellPlane(InputIterator pnts_begin,
15  InputIterator pnts_end)
16 {
17  Eigen::Vector3d plane_normal({0, 0, 0});
18  Eigen::Vector3d centroid({0, 0, 0});
19  for (auto i=std::prev(pnts_end), j=pnts_begin; j!=pnts_end; i = j, ++j) {
20  auto &pt_i = *(*i);
21  auto &pt_j = *(*j);
22  plane_normal[0] += (pt_i[1] - pt_j[1])
23  * (pt_i[2] + pt_j[2]); // projection on yz
24  plane_normal[1] += (pt_i[2] - pt_j[2])
25  * (pt_i[0] + pt_j[0]); // projection on xz
26  plane_normal[2] += (pt_i[0] - pt_j[0])
27  * (pt_i[1] + pt_j[1]); // projection on xy
28 
29  centroid += Eigen::Map<Eigen::Vector3d const>(pt_j.getCoords());
30  }
31 
32  plane_normal.normalize();
33  auto n_pnts(std::distance(pnts_begin, pnts_end));
34  assert(n_pnts > 2);
35  double d = 0.0;
36  if (n_pnts > 0)
37  {
38  d = centroid.dot(plane_normal) / n_pnts;
39  }
40  return std::make_pair(plane_normal, d);
41 }
42 
43 template <class T_POINT>
44 std::pair<Eigen::Vector3d, double> getNewellPlane(
45  const std::vector<T_POINT*>& pnts)
46 {
47  return getNewellPlane(pnts.begin(), pnts.end());
48 }
49 
50 template <class T_POINT>
51 std::pair<Eigen::Vector3d, double> getNewellPlane(
52  const std::vector<T_POINT>& pnts)
53 {
54  Eigen::Vector3d plane_normal({0, 0, 0});
55  Eigen::Vector3d centroid({0, 0, 0});
56  std::size_t n_pnts(pnts.size());
57  for (std::size_t i = n_pnts - 1, j = 0; j < n_pnts; i = j, j++) {
58  plane_normal[0] += (pnts[i][1] - pnts[j][1])
59  * (pnts[i][2] + pnts[j][2]); // projection on yz
60  plane_normal[1] += (pnts[i][2] - pnts[j][2])
61  * (pnts[i][0] + pnts[j][0]); // projection on xz
62  plane_normal[2] += (pnts[i][0] - pnts[j][0])
63  * (pnts[i][1] + pnts[j][1]); // projection on xy
64 
65  centroid += Eigen::Map<Eigen::Vector3d const>(pnts[j].getCoords());
66  }
67 
68  plane_normal.normalize();
69  double const d = centroid.dot(plane_normal) / n_pnts;
70  return std::make_pair(plane_normal, d);
71 }
72 
73 
74 template <typename InputIterator>
75 void rotatePoints(Eigen::Matrix3d const& rot_mat, InputIterator pnts_begin,
76  InputIterator pnts_end)
77 {
78  for (auto it = pnts_begin; it != pnts_end; ++it)
79  {
80  Eigen::Map<Eigen::Vector3d>((*it)->getCoords()) =
81  rot_mat * Eigen::Map<Eigen::Vector3d const>((*it)->getCoords());
82  }
83 }
84 
85 template <typename InputIterator1, typename InputIterator2>
86 Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin,
87  InputIterator1 p_pnts_end,
88  InputIterator2 r_pnts_begin,
89  InputIterator2 r_pnts_end)
90 {
91  assert(std::distance(p_pnts_begin, p_pnts_end) > 2);
92 
93  // compute the plane normal
94  auto const [plane_normal, d] =
95  GeoLib::getNewellPlane(p_pnts_begin, p_pnts_end);
96 
97  // rotate points into x-y-plane
98  Eigen::Matrix3d const rot_mat = computeRotationMatrixToXY(plane_normal);
99  rotatePoints(rot_mat, r_pnts_begin, r_pnts_end);
100 
101  for (auto it = r_pnts_begin; it != r_pnts_end; ++it)
102  {
103  (*(*it))[2] = 0.0; // should be -= d but there are numerical errors
104  }
105 
106  return rot_mat;
107 }
108 
109 template <typename P>
110 void rotatePoints(Eigen::Matrix3d const& rot_mat, std::vector<P*> const& pnts)
111 {
112  rotatePoints(rot_mat, pnts.begin(), pnts.end());
113 }
114 
115 } // end namespace GeoLib
116 
void rotatePoints(Eigen::Matrix3d const &rot_mat, InputIterator pnts_begin, InputIterator pnts_end)
Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const &n)
Eigen::Matrix3d rotatePointsToXY(InputIterator1 p_pnts_begin, InputIterator1 p_pnts_end, InputIterator2 r_pnts_begin, InputIterator2 r_pnts_end)
std::pair< Eigen::Vector3d, double > getNewellPlane(InputIterator pnts_begin, InputIterator pnts_end)