11#ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_
12#define GARLAND_HECKBERT_ERROR_QUADRIC_H_
14#include <boost/optional/optional.hpp>
19template <
typename Po
int>
class Error_quadric {
49 Error_quadric(
const Point & p0,
const Point & p1,
const Point & p2) {
50 Point normal(unit_normal(p0, p1, p2));
54 double d = -a * p0[0] - b * p0[1] - c * p0[2];
66 double area_p0p1p2 = std::sqrt(squared_area(p0, p1, p2));
71 inline double squared_area(
const Point& p0,
const Point& p1,
const Point& p2) {
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);
90 Error_quadric& operator+=(
const Error_quadric& other) {
92 for (
int i = 0; i < 10; ++i)
93 coeff[i] += other.coeff[i];
101 inline double cost(
const Point& point)
const {
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())
114 inline double grad_determinant()
const {
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];
127 inline Point solve_linear_gradient(
double det)
const {
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])
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])
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])
146 boost::optional<Point> min_cost(
double scale = 1)
const {
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);
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] <<
"]";
std::ostream & operator<<(std::ostream &os, const Permutahedral_representation< Vertex, OrderedSetPartition > &simplex)
Print a permutahedral representation to a stream.
Definition: Permutahedral_representation.h:173