From 14f55cf58e9be168e6be635ddafebc6c86cc7eea Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 31 May 2017 16:48:08 +0000 Subject: Convert more examples into utilities git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2494 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 54fa703f0489cba059a8df973738118f6a092f70 --- src/Contraction/example/CMakeLists.txt | 5 +- src/Contraction/example/Garland_heckbert.cpp | 195 --------------------- .../example/Garland_heckbert/Error_quadric.h | 182 ------------------- src/Contraction/utilities/CMakeLists.txt | 7 + src/Contraction/utilities/Garland_heckbert.cpp | 195 +++++++++++++++++++++ .../utilities/Garland_heckbert/Error_quadric.h | 182 +++++++++++++++++++ 6 files changed, 385 insertions(+), 381 deletions(-) delete mode 100644 src/Contraction/example/Garland_heckbert.cpp delete mode 100644 src/Contraction/example/Garland_heckbert/Error_quadric.h create mode 100644 src/Contraction/utilities/CMakeLists.txt create mode 100644 src/Contraction/utilities/Garland_heckbert.cpp create mode 100644 src/Contraction/utilities/Garland_heckbert/Error_quadric.h (limited to 'src/Contraction') diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt index b2b38dea..5cd32daf 100644 --- a/src/Contraction/example/CMakeLists.txt +++ b/src/Contraction/example/CMakeLists.txt @@ -3,10 +3,8 @@ project(Contraction_examples) add_executable(RipsContraction Rips_contraction.cpp) -add_executable(GarlandHeckbert Garland_heckbert.cpp) -target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) -target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) +target_link_libraries(RipsContraction ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) add_test(NAME Contraction_example_tore3D_0.2 COMMAND $ @@ -18,4 +16,3 @@ add_test(NAME Contraction_example_tore3D_0.2 COMMAND $. - * - */ - - -#ifndef GARLAND_HECKBERT_H_ -#define GARLAND_HECKBERT_H_ - -#include -#include -#include -#include - -#include -#include - -#include "Garland_heckbert/Error_quadric.h" - -struct Geometry_trait { - typedef Point_d Point; -}; - -/** - * The vertex stored in the complex contains a quadric. - */ -struct Garland_heckbert_traits - : public Gudhi::skeleton_blocker::Skeleton_blocker_simple_geometric_traits { - public: - struct Garland_heckbert_vertex : public Simple_geometric_vertex { - Error_quadric quadric; - }; - typedef Garland_heckbert_vertex Graph_vertex; -}; - -using Complex = Gudhi::skeleton_blocker::Skeleton_blocker_geometric_complex< Garland_heckbert_traits >; -using EdgeProfile = Gudhi::contraction::Edge_profile; -using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor; - -/** - * How the new vertex is placed after an edge collapse : here it is placed at - * the point minimizing the cost of the quadric. - */ -class GH_placement : public Gudhi::contraction::Placement_policy { - Complex& complex_; - - public: - typedef Gudhi::contraction::Placement_policy::Placement_type Placement_type; - - GH_placement(Complex& complex) : complex_(complex) { (void)complex_; } - - Placement_type operator()(const EdgeProfile& profile) const override { - auto sum_quad(profile.v0().quadric); - sum_quad += profile.v1().quadric; - - boost::optional min_quadric_pt(sum_quad.min_cost()); - if (min_quadric_pt) - return Placement_type(*min_quadric_pt); - else - return profile.p0(); - } -}; - -/** - * How much cost an edge collapse : here the costs is given by a quadric - * which expresses a squared distances with triangles planes. - */ -class GH_cost : public Gudhi::contraction::Cost_policy { - Complex& complex_; - - public: - typedef Gudhi::contraction::Cost_policy::Cost_type Cost_type; - - GH_cost(Complex& complex) : complex_(complex) { (void)complex_; } - - Cost_type operator()(EdgeProfile const& profile, boost::optional const& new_point) const override { - Cost_type res; - if (new_point) { - auto sum_quad(profile.v0().quadric); - sum_quad += profile.v1().quadric; - res = sum_quad.cost(*new_point); - } - return res; - } -}; - -/** - * Visitor that is called at several moment. - * Here we initializes the quadrics of every vertex at the on_started call back - * and we update them when contracting an edge (the quadric become the sum of both quadrics). - */ -class GH_visitor : public Gudhi::contraction::Contraction_visitor { - Complex& complex_; - - public: - GH_visitor(Complex& complex) : complex_(complex) { (void)complex_; } - - // Compute quadrics for every vertex v - // The quadric of v consists in the sum of quadric - // of every triangles passing through v weighted by its area - - void on_started(Complex & complex) override { - for (auto v : complex.vertex_range()) { - auto & quadric_v(complex[v].quadric); - for (auto t : complex.triangle_range(v)) { - auto t_it = t.begin(); - const auto& p0(complex.point(*t_it++)); - const auto& p1(complex.point(*t_it++)); - const auto& p2(complex.point(*t_it++)); - quadric_v += Error_quadric(p0, p1, p2); - } - } - } - - /** - * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. - */ - void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement) - override { - profile.v0().quadric += profile.v1().quadric; - } -}; - -int main(int argc, char *argv[]) { - if (argc != 4) { - std::cerr << "Usage " << argv[0] << - " input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n"; - return EXIT_FAILURE; - } - - Complex complex; - typedef Complex::Vertex_handle Vertex_handle; - - // load the points - Gudhi::skeleton_blocker::Skeleton_blocker_off_reader off_reader(argv[1], complex); - if (!off_reader.is_valid()) { - std::cerr << "Unable to read file:" << argv[1] << std::endl; - return EXIT_FAILURE; - } - - if (!complex.empty() && !(complex.point(Vertex_handle(0)).dimension() == 3)) { - std::cerr << "Only points of dimension 3 are supported." << std::endl; - return EXIT_FAILURE; - } - - std::cout << "Load complex with " << complex.num_vertices() << " vertices" << std::endl; - - int num_contractions = atoi(argv[3]); - - boost::timer::auto_cpu_timer t; - - // constructs the contractor object with Garland Heckbert policies. - Complex_contractor contractor(complex, - new GH_cost(complex), - new GH_placement(complex), - Gudhi::contraction::make_link_valid_contraction(), - new GH_visitor(complex)); - - std::cout << "Contract " << num_contractions << " edges" << std::endl; - contractor.contract_edges(num_contractions); - - std::cout << "Final complex has " << - complex.num_vertices() << " vertices, " << - complex.num_edges() << " edges and " << - complex.num_triangles() << " triangles." << std::endl; - - // write simplified complex - Gudhi::skeleton_blocker::Skeleton_blocker_off_writer off_writer(argv[2], complex); - - return EXIT_SUCCESS; -} - -#endif // GARLAND_HECKBERT_H_ - - - - diff --git a/src/Contraction/example/Garland_heckbert/Error_quadric.h b/src/Contraction/example/Garland_heckbert/Error_quadric.h deleted file mode 100644 index e7dafaa0..00000000 --- a/src/Contraction/example/Garland_heckbert/Error_quadric.h +++ /dev/null @@ -1,182 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): David Salinas - * - * Copyright (C) 2014 INRIA Sophia Antipolis-M�diterran�e (France) - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - * - */ - -#ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_ -#define GARLAND_HECKBERT_ERROR_QUADRIC_H_ - -#include - -#include -#include - -template class Error_quadric { - private: - double coeff[10]; - - public: - Error_quadric() { - clear(); - } - - /** - * Quadric corresponding to the L2 distance to the plane. - * - * According to the notation of Garland Heckbert, they - * denote a quadric symetric matrix as : - * Q = [ q11 q12 q13 q14] - * [ q12 q22 q23 q24] - * [ q13 q23 q33 q34] - * [ q14 q24 q34 q44] - * - * which is represented by a vector with 10 elts that - * are denoted ci for clarity with : - * Q = [ c0 c1 c2 c3 ] - * [ c1 c4 c5 c6 ] - * [ c2 c5 c7 c8 ] - * [ c3 c6 c8 c9 ] - * - * The constructor return the quadrics that represents - * the squared distance to the plane defined by triangle p0,p1,p2 - * times the area of triangle p0,p1,p2. - */ - Error_quadric(const Point & p0, const Point & p1, const Point & p2) { - Point normal(unit_normal(p0, p1, p2)); - double a = normal[0]; - double b = normal[1]; - double c = normal[2]; - double d = -a * p0[0] - b * p0[1] - c * p0[2]; - coeff[0] = a*a; - coeff[1] = a*b; - coeff[2] = a*c; - coeff[3] = a*d; - coeff[4] = b*b; - coeff[5] = b*c; - coeff[6] = b*d; - coeff[7] = c*c; - coeff[8] = c*d; - coeff[9] = d*d; - - double area_p0p1p2 = std::sqrt(squared_area(p0, p1, p2)); - for (auto& x : coeff) - x *= area_p0p1p2; - } - - inline double squared_area(const Point& p0, const Point& p1, const Point& p2) { - // if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0 - // then the squared area is = (u^2+v^2+w^2)/4 - // with: u = x2 * y3 - x3 * y2; - // v = x3 * y1 - x1 * y3; - // w = x1 * y2 - x2 * y1; - Point p0p1(p1 - p0); - Point p0p2(p2 - p0); - double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1]; - double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2]; - double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0]; - return 1. / 4. * (A * A + B * B + C * C); - } - - void clear() { - for (auto& x : coeff) - x = 0; - } - - Error_quadric& operator+=(const Error_quadric& other) { - if (this != &other) { - for (int i = 0; i < 10; ++i) - coeff[i] += other.coeff[i]; - } - return *this; - } - - /** - * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert - */ - inline double cost(const Point& point) const { - double cost = - coeff[0] * point.x() * point.x() + coeff[4] * point.y() * point.y() + coeff[7] * point.z() * point.z() - + 2 * (coeff[1] * point.x() * point.y() + coeff[5] * point.y() * point.z() + coeff[2] * point.z() * point.x()) - + 2 * (coeff[3] * point.x() + coeff[6] * point.y() + coeff[8] * point.z()) - + coeff[9]; - if (cost < 0) { - return 0; - } else { - return cost; - } - } - - inline double grad_determinant() const { - return - coeff[0] * coeff[4] * coeff[7] - - coeff[0] * coeff[5] * coeff[5] - - coeff[1] * coeff[1] * coeff[7] - + 2 * coeff[1] * coeff[5] * coeff[2] - - coeff[4] * coeff[2] * coeff[2]; - } - - /** - * Return the point such that it minimizes the gradient of the quadric. - * Det must be passed with the determinant value of the gradient (should be non zero). - */ - inline Point solve_linear_gradient(double det) const { - return Point({ - (-coeff[1] * coeff[5] * coeff[8] + coeff[1] * coeff[7] * coeff[6] + coeff[2] * coeff[8] * coeff[4] - - coeff[2] * coeff[5] * coeff[6] - coeff[3] * coeff[4] * coeff[7] + coeff[3] * coeff[5] * coeff[5]) - / det, - (coeff[0] * coeff[5] * coeff[8] - coeff[0] * coeff[7] * coeff[6] - coeff[5] * coeff[2] * coeff[3] - - coeff[1] * coeff[2] * coeff[8] + coeff[6] * coeff[2] * coeff[2] + coeff[1] * coeff[3] * coeff[7]) - / det, - (-coeff[8] * coeff[0] * coeff[4] + coeff[8] * coeff[1] * coeff[1] + coeff[2] * coeff[3] * coeff[4] + - coeff[5] * coeff[0] * coeff[6] - coeff[5] * coeff[1] * coeff[3] - coeff[1] * coeff[2] * coeff[6]) - / det - }); - } - - /** - * returns the point that minimizes the quadric. - * It inverses the quadric if its determinant is higher that a given threshold . - * If the determinant is lower than this value the returned value is uninitialized. - */ - boost::optional min_cost(double scale = 1) const { - // const double min_determinant = 1e-4 * scale*scale; - const double min_determinant = 1e-5; - boost::optional pt_res; - double det = grad_determinant(); - if (std::abs(det) > min_determinant) - pt_res = solve_linear_gradient(det); - return pt_res; - } - - friend std::ostream& operator<<(std::ostream& stream, const Error_quadric& quadric) { - stream << "\n[ " << quadric.coeff[0] << "," << quadric.coeff[1] << "," << quadric.coeff[2] << "," << - quadric.coeff[3] << ";\n"; - stream << " " << quadric.coeff[1] << "," << quadric.coeff[4] << "," << quadric.coeff[5] << "," << - quadric.coeff[6] << ";\n"; - stream << " " << quadric.coeff[2] << "," << quadric.coeff[5] << "," << quadric.coeff[7] << "," << - quadric.coeff[8] << ";\n"; - stream << " " << quadric.coeff[3] << "," << quadric.coeff[6] << "," << quadric.coeff[8] << "," << - quadric.coeff[9] << "]"; - return stream; - } -}; - -#endif // GARLAND_HECKBERT_ERROR_QUADRIC_H_ diff --git a/src/Contraction/utilities/CMakeLists.txt b/src/Contraction/utilities/CMakeLists.txt new file mode 100644 index 00000000..77aa99d5 --- /dev/null +++ b/src/Contraction/utilities/CMakeLists.txt @@ -0,0 +1,7 @@ +cmake_minimum_required(VERSION 2.6) +project(Contraction_utilities) + +add_executable(GarlandHeckbert Garland_heckbert.cpp) +target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY} ${Boost_SYSTEM_LIBRARY}) + +install(TARGETS GarlandHeckbert DESTINATION bin) diff --git a/src/Contraction/utilities/Garland_heckbert.cpp b/src/Contraction/utilities/Garland_heckbert.cpp new file mode 100644 index 00000000..8b5a6a6c --- /dev/null +++ b/src/Contraction/utilities/Garland_heckbert.cpp @@ -0,0 +1,195 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-M�diterran�e (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + */ + + +#ifndef GARLAND_HECKBERT_H_ +#define GARLAND_HECKBERT_H_ + +#include +#include +#include +#include + +#include +#include + +#include "Garland_heckbert/Error_quadric.h" + +struct Geometry_trait { + typedef Point_d Point; +}; + +/** + * The vertex stored in the complex contains a quadric. + */ +struct Garland_heckbert_traits + : public Gudhi::skeleton_blocker::Skeleton_blocker_simple_geometric_traits { + public: + struct Garland_heckbert_vertex : public Simple_geometric_vertex { + Error_quadric quadric; + }; + typedef Garland_heckbert_vertex Graph_vertex; +}; + +using Complex = Gudhi::skeleton_blocker::Skeleton_blocker_geometric_complex< Garland_heckbert_traits >; +using EdgeProfile = Gudhi::contraction::Edge_profile; +using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor; + +/** + * How the new vertex is placed after an edge collapse : here it is placed at + * the point minimizing the cost of the quadric. + */ +class GH_placement : public Gudhi::contraction::Placement_policy { + Complex& complex_; + + public: + typedef Gudhi::contraction::Placement_policy::Placement_type Placement_type; + + GH_placement(Complex& complex) : complex_(complex) { (void)complex_; } + + Placement_type operator()(const EdgeProfile& profile) const override { + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + + boost::optional min_quadric_pt(sum_quad.min_cost()); + if (min_quadric_pt) + return Placement_type(*min_quadric_pt); + else + return profile.p0(); + } +}; + +/** + * How much cost an edge collapse : here the costs is given by a quadric + * which expresses a squared distances with triangles planes. + */ +class GH_cost : public Gudhi::contraction::Cost_policy { + Complex& complex_; + + public: + typedef Gudhi::contraction::Cost_policy::Cost_type Cost_type; + + GH_cost(Complex& complex) : complex_(complex) { (void)complex_; } + + Cost_type operator()(EdgeProfile const& profile, boost::optional const& new_point) const override { + Cost_type res; + if (new_point) { + auto sum_quad(profile.v0().quadric); + sum_quad += profile.v1().quadric; + res = sum_quad.cost(*new_point); + } + return res; + } +}; + +/** + * Visitor that is called at several moment. + * Here we initializes the quadrics of every vertex at the on_started call back + * and we update them when contracting an edge (the quadric become the sum of both quadrics). + */ +class GH_visitor : public Gudhi::contraction::Contraction_visitor { + Complex& complex_; + + public: + GH_visitor(Complex& complex) : complex_(complex) { (void)complex_; } + + // Compute quadrics for every vertex v + // The quadric of v consists in the sum of quadric + // of every triangles passing through v weighted by its area + + void on_started(Complex & complex) override { + for (auto v : complex.vertex_range()) { + auto & quadric_v(complex[v].quadric); + for (auto t : complex.triangle_range(v)) { + auto t_it = t.begin(); + const auto& p0(complex.point(*t_it++)); + const auto& p1(complex.point(*t_it++)); + const auto& p2(complex.point(*t_it++)); + quadric_v += Error_quadric(p0, p1, p2); + } + } + } + + /** + * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. + */ + void on_contracting(EdgeProfile const &profile, boost::optional< Point > placement) + override { + profile.v0().quadric += profile.v1().quadric; + } +}; + +int main(int argc, char *argv[]) { + if (argc != 4) { + std::cerr << "Usage " << argv[0] << + " input.off output.off N to load the file input.off, contract N edges and save the result to output.off.\n"; + return EXIT_FAILURE; + } + + Complex complex; + typedef Complex::Vertex_handle Vertex_handle; + + // load the points + Gudhi::skeleton_blocker::Skeleton_blocker_off_reader off_reader(argv[1], complex); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file:" << argv[1] << std::endl; + return EXIT_FAILURE; + } + + if (!complex.empty() && !(complex.point(Vertex_handle(0)).dimension() == 3)) { + std::cerr << "Only points of dimension 3 are supported." << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Load complex with " << complex.num_vertices() << " vertices" << std::endl; + + int num_contractions = atoi(argv[3]); + + boost::timer::auto_cpu_timer t; + + // constructs the contractor object with Garland Heckbert policies. + Complex_contractor contractor(complex, + new GH_cost(complex), + new GH_placement(complex), + Gudhi::contraction::make_link_valid_contraction(), + new GH_visitor(complex)); + + std::cout << "Contract " << num_contractions << " edges" << std::endl; + contractor.contract_edges(num_contractions); + + std::cout << "Final complex has " << + complex.num_vertices() << " vertices, " << + complex.num_edges() << " edges and " << + complex.num_triangles() << " triangles." << std::endl; + + // write simplified complex + Gudhi::skeleton_blocker::Skeleton_blocker_off_writer off_writer(argv[2], complex); + + return EXIT_SUCCESS; +} + +#endif // GARLAND_HECKBERT_H_ + + + + diff --git a/src/Contraction/utilities/Garland_heckbert/Error_quadric.h b/src/Contraction/utilities/Garland_heckbert/Error_quadric.h new file mode 100644 index 00000000..e7dafaa0 --- /dev/null +++ b/src/Contraction/utilities/Garland_heckbert/Error_quadric.h @@ -0,0 +1,182 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): David Salinas + * + * Copyright (C) 2014 INRIA Sophia Antipolis-M�diterran�e (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + */ + +#ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_ +#define GARLAND_HECKBERT_ERROR_QUADRIC_H_ + +#include + +#include +#include + +template class Error_quadric { + private: + double coeff[10]; + + public: + Error_quadric() { + clear(); + } + + /** + * Quadric corresponding to the L2 distance to the plane. + * + * According to the notation of Garland Heckbert, they + * denote a quadric symetric matrix as : + * Q = [ q11 q12 q13 q14] + * [ q12 q22 q23 q24] + * [ q13 q23 q33 q34] + * [ q14 q24 q34 q44] + * + * which is represented by a vector with 10 elts that + * are denoted ci for clarity with : + * Q = [ c0 c1 c2 c3 ] + * [ c1 c4 c5 c6 ] + * [ c2 c5 c7 c8 ] + * [ c3 c6 c8 c9 ] + * + * The constructor return the quadrics that represents + * the squared distance to the plane defined by triangle p0,p1,p2 + * times the area of triangle p0,p1,p2. + */ + Error_quadric(const Point & p0, const Point & p1, const Point & p2) { + Point normal(unit_normal(p0, p1, p2)); + double a = normal[0]; + double b = normal[1]; + double c = normal[2]; + double d = -a * p0[0] - b * p0[1] - c * p0[2]; + coeff[0] = a*a; + coeff[1] = a*b; + coeff[2] = a*c; + coeff[3] = a*d; + coeff[4] = b*b; + coeff[5] = b*c; + coeff[6] = b*d; + coeff[7] = c*c; + coeff[8] = c*d; + coeff[9] = d*d; + + double area_p0p1p2 = std::sqrt(squared_area(p0, p1, p2)); + for (auto& x : coeff) + x *= area_p0p1p2; + } + + inline double squared_area(const Point& p0, const Point& p1, const Point& p2) { + // if (x1,x2,x3) = p1-p0 and (y1,y2,y3) = p2-p0 + // then the squared area is = (u^2+v^2+w^2)/4 + // with: u = x2 * y3 - x3 * y2; + // v = x3 * y1 - x1 * y3; + // w = x1 * y2 - x2 * y1; + Point p0p1(p1 - p0); + Point p0p2(p2 - p0); + double A = p0p1[1] * p0p2[2] - p0p1[2] * p0p2[1]; + double B = p0p1[2] * p0p2[0] - p0p1[0] * p0p2[2]; + double C = p0p1[0] * p0p2[1] - p0p1[1] * p0p2[0]; + return 1. / 4. * (A * A + B * B + C * C); + } + + void clear() { + for (auto& x : coeff) + x = 0; + } + + Error_quadric& operator+=(const Error_quadric& other) { + if (this != &other) { + for (int i = 0; i < 10; ++i) + coeff[i] += other.coeff[i]; + } + return *this; + } + + /** + * @return The quadric quost defined by the scalar product v^T Q v where Q is the quadratic form of Garland/Heckbert + */ + inline double cost(const Point& point) const { + double cost = + coeff[0] * point.x() * point.x() + coeff[4] * point.y() * point.y() + coeff[7] * point.z() * point.z() + + 2 * (coeff[1] * point.x() * point.y() + coeff[5] * point.y() * point.z() + coeff[2] * point.z() * point.x()) + + 2 * (coeff[3] * point.x() + coeff[6] * point.y() + coeff[8] * point.z()) + + coeff[9]; + if (cost < 0) { + return 0; + } else { + return cost; + } + } + + inline double grad_determinant() const { + return + coeff[0] * coeff[4] * coeff[7] + - coeff[0] * coeff[5] * coeff[5] + - coeff[1] * coeff[1] * coeff[7] + + 2 * coeff[1] * coeff[5] * coeff[2] + - coeff[4] * coeff[2] * coeff[2]; + } + + /** + * Return the point such that it minimizes the gradient of the quadric. + * Det must be passed with the determinant value of the gradient (should be non zero). + */ + inline Point solve_linear_gradient(double det) const { + return Point({ + (-coeff[1] * coeff[5] * coeff[8] + coeff[1] * coeff[7] * coeff[6] + coeff[2] * coeff[8] * coeff[4] - + coeff[2] * coeff[5] * coeff[6] - coeff[3] * coeff[4] * coeff[7] + coeff[3] * coeff[5] * coeff[5]) + / det, + (coeff[0] * coeff[5] * coeff[8] - coeff[0] * coeff[7] * coeff[6] - coeff[5] * coeff[2] * coeff[3] - + coeff[1] * coeff[2] * coeff[8] + coeff[6] * coeff[2] * coeff[2] + coeff[1] * coeff[3] * coeff[7]) + / det, + (-coeff[8] * coeff[0] * coeff[4] + coeff[8] * coeff[1] * coeff[1] + coeff[2] * coeff[3] * coeff[4] + + coeff[5] * coeff[0] * coeff[6] - coeff[5] * coeff[1] * coeff[3] - coeff[1] * coeff[2] * coeff[6]) + / det + }); + } + + /** + * returns the point that minimizes the quadric. + * It inverses the quadric if its determinant is higher that a given threshold . + * If the determinant is lower than this value the returned value is uninitialized. + */ + boost::optional min_cost(double scale = 1) const { + // const double min_determinant = 1e-4 * scale*scale; + const double min_determinant = 1e-5; + boost::optional pt_res; + double det = grad_determinant(); + if (std::abs(det) > min_determinant) + pt_res = solve_linear_gradient(det); + return pt_res; + } + + friend std::ostream& operator<<(std::ostream& stream, const Error_quadric& quadric) { + stream << "\n[ " << quadric.coeff[0] << "," << quadric.coeff[1] << "," << quadric.coeff[2] << "," << + quadric.coeff[3] << ";\n"; + stream << " " << quadric.coeff[1] << "," << quadric.coeff[4] << "," << quadric.coeff[5] << "," << + quadric.coeff[6] << ";\n"; + stream << " " << quadric.coeff[2] << "," << quadric.coeff[5] << "," << quadric.coeff[7] << "," << + quadric.coeff[8] << ";\n"; + stream << " " << quadric.coeff[3] << "," << quadric.coeff[6] << "," << quadric.coeff[8] << "," << + quadric.coeff[9] << "]"; + return stream; + } +}; + +#endif // GARLAND_HECKBERT_ERROR_QUADRIC_H_ -- cgit v1.2.3