Error_quadric.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): David Salinas
4  *
5  * Copyright (C) 2014 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10 
11 #ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_
12 #define GARLAND_HECKBERT_ERROR_QUADRIC_H_
13 
14 #include <boost/optional/optional.hpp>
15 
16 #include <vector>
17 #include <utility>
18 
19 template <typename Point> class Error_quadric {
20  private:
21  double coeff[10];
22 
23  public:
24  Error_quadric() {
25  clear();
26  }
27 
49  Error_quadric(const Point & p0, const Point & p1, const Point & p2) {
50  Point normal(unit_normal(p0, p1, p2));
51  double a = normal[0];
52  double b = normal[1];
53  double c = normal[2];
54  double d = -a * p0[0] - b * p0[1] - c * p0[2];
55  coeff[0] = a*a;
56  coeff[1] = a*b;
57  coeff[2] = a*c;
58  coeff[3] = a*d;
59  coeff[4] = b*b;
60  coeff[5] = b*c;
61  coeff[6] = b*d;
62  coeff[7] = c*c;
63  coeff[8] = c*d;
64  coeff[9] = d*d;
65 
66  double area_p0p1p2 = std::sqrt(squared_area(p0, p1, p2));
67  for (auto& x : coeff)
68  x *= area_p0p1p2;
69  }
70 
71  inline double squared_area(const Point& p0, const Point& p1, const Point& p2) {
72  // if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0
73  // then the squared area is = (u^2+v^2+w^2)/4
74  // with: u = x2 * y3 - x3 * y2;
75  // v = x3 * y1 - x1 * y3;
76  // w = x1 * y2 - x2 * y1;
77  Point p0p1(p1 - p0);
78  Point p0p2(p2 - p0);
79  double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1];
80  double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2];
81  double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0];
82  return 1. / 4. * (A * A + B * B + C * C);
83  }
84 
85  void clear() {
86  for (auto& x : coeff)
87  x = 0;
88  }
89 
90  Error_quadric& operator+=(const Error_quadric& other) {
91  if (this != &other) {
92  for (int i = 0; i < 10; ++i)
93  coeff[i] += other.coeff[i];
94  }
95  return *this;
96  }
97 
101  inline double cost(const Point& point) const {
102  double cost =
103  coeff[0] * point.x() * point.x() + coeff[4] * point.y() * point.y() + coeff[7] * point.z() * point.z()
104  + 2 * (coeff[1] * point.x() * point.y() + coeff[5] * point.y() * point.z() + coeff[2] * point.z() * point.x())
105  + 2 * (coeff[3] * point.x() + coeff[6] * point.y() + coeff[8] * point.z())
106  + coeff[9];
107  if (cost < 0) {
108  return 0;
109  } else {
110  return cost;
111  }
112  }
113 
114  inline double grad_determinant() const {
115  return
116  coeff[0] * coeff[4] * coeff[7]
117  - coeff[0] * coeff[5] * coeff[5]
118  - coeff[1] * coeff[1] * coeff[7]
119  + 2 * coeff[1] * coeff[5] * coeff[2]
120  - coeff[4] * coeff[2] * coeff[2];
121  }
122 
127  inline Point solve_linear_gradient(double det) const {
128  return Point({
129  (-coeff[1] * coeff[5] * coeff[8] + coeff[1] * coeff[7] * coeff[6] + coeff[2] * coeff[8] * coeff[4] -
130  coeff[2] * coeff[5] * coeff[6] - coeff[3] * coeff[4] * coeff[7] + coeff[3] * coeff[5] * coeff[5])
131  / det,
132  (coeff[0] * coeff[5] * coeff[8] - coeff[0] * coeff[7] * coeff[6] - coeff[5] * coeff[2] * coeff[3] -
133  coeff[1] * coeff[2] * coeff[8] + coeff[6] * coeff[2] * coeff[2] + coeff[1] * coeff[3] * coeff[7])
134  / det,
135  (-coeff[8] * coeff[0] * coeff[4] + coeff[8] * coeff[1] * coeff[1] + coeff[2] * coeff[3] * coeff[4] +
136  coeff[5] * coeff[0] * coeff[6] - coeff[5] * coeff[1] * coeff[3] - coeff[1] * coeff[2] * coeff[6])
137  / det
138  });
139  }
140 
146  boost::optional<Point> min_cost(double scale = 1) const {
147  // const double min_determinant = 1e-4 * scale*scale;
148  const double min_determinant = 1e-5;
149  boost::optional<Point> pt_res;
150  double det = grad_determinant();
151  if (std::abs(det) > min_determinant)
152  pt_res = solve_linear_gradient(det);
153  return pt_res;
154  }
155 
156  friend std::ostream& operator<<(std::ostream& stream, const Error_quadric& quadric) {
157  stream << "\n[ " << quadric.coeff[0] << "," << quadric.coeff[1] << "," << quadric.coeff[2] << "," <<
158  quadric.coeff[3] << ";\n";
159  stream << " " << quadric.coeff[1] << "," << quadric.coeff[4] << "," << quadric.coeff[5] << "," <<
160  quadric.coeff[6] << ";\n";
161  stream << " " << quadric.coeff[2] << "," << quadric.coeff[5] << "," << quadric.coeff[7] << "," <<
162  quadric.coeff[8] << ";\n";
163  stream << " " << quadric.coeff[3] << "," << quadric.coeff[6] << "," << quadric.coeff[8] << "," <<
164  quadric.coeff[9] << "]";
165  return stream;
166  }
167 };
168 
169 #endif // GARLAND_HECKBERT_ERROR_QUADRIC_H_
std::ostream & operator<<(std::ostream &os, const Permutahedral_representation< Vertex, OrderedSetPartition > &simplex)
Print a permutahedral representation to a stream.
Definition: Permutahedral_representation.h:173
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Fri Jan 14 2022 18:28:42 for GUDHIdev by Doxygen 1.9.1