diff options
Diffstat (limited to 'src/Contraction')
24 files changed, 1988 insertions, 0 deletions
diff --git a/src/Contraction/doc/COPYRIGHT b/src/Contraction/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Contraction/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + +Author(s): Vincent Rouvreau + +Copyright (C) 2015 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT.
\ No newline at end of file diff --git a/src/Contraction/doc/SO3_rips.png b/src/Contraction/doc/SO3_rips.png Binary files differnew file mode 100644 index 00000000..60452f86 --- /dev/null +++ b/src/Contraction/doc/SO3_rips.png diff --git a/src/Contraction/doc/SO3_simplified.png b/src/Contraction/doc/SO3_simplified.png Binary files differnew file mode 100644 index 00000000..f70a1903 --- /dev/null +++ b/src/Contraction/doc/SO3_simplified.png diff --git a/src/Contraction/doc/SO3points.png b/src/Contraction/doc/SO3points.png Binary files differnew file mode 100644 index 00000000..0362d98f --- /dev/null +++ b/src/Contraction/doc/SO3points.png diff --git a/src/Contraction/doc/so3.png b/src/Contraction/doc/so3.png Binary files differnew file mode 100644 index 00000000..e66acae1 --- /dev/null +++ b/src/Contraction/doc/so3.png diff --git a/src/Contraction/doc/so3.svg b/src/Contraction/doc/so3.svg new file mode 100644 index 00000000..adea3f38 --- /dev/null +++ b/src/Contraction/doc/so3.svg @@ -0,0 +1,209 @@ +<?xml version="1.0" encoding="UTF-8" standalone="no"?> +<!-- Created with Inkscape (http://www.inkscape.org/) --> + +<svg + xmlns:dc="http://purl.org/dc/elements/1.1/" + xmlns:cc="http://creativecommons.org/ns#" + xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" + xmlns:svg="http://www.w3.org/2000/svg" + xmlns="http://www.w3.org/2000/svg" + xmlns:xlink="http://www.w3.org/1999/xlink" + xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd" + xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape" + width="744.09448819" + height="1052.3622047" + id="svg2" + version="1.1" + inkscape:version="0.48.4 r9939" + sodipodi:docname="so3.svg" + inkscape:export-filename="/home/dsalinas/Documents/CodeSVN/gudhi_depot/trunk/src/Contraction/doc/so3.png" + inkscape:export-xdpi="200.20428" + inkscape:export-ydpi="200.20428"> + <defs + id="defs4"> + <marker + inkscape:stockid="Arrow1Lend" + orient="auto" + refY="0.0" + refX="0.0" + id="Arrow1Lend" + style="overflow:visible;"> + <path + id="path3888" + d="M 0.0,0.0 L 5.0,-5.0 L -12.5,0.0 L 5.0,5.0 L 0.0,0.0 z " + style="fill-rule:evenodd;stroke:#000000;stroke-width:1.0pt;" + transform="scale(0.8) rotate(180) translate(12.5,0)" /> + </marker> + </defs> + <sodipodi:namedview + id="base" + pagecolor="#ffffff" + bordercolor="#666666" + borderopacity="1.0" + inkscape:pageopacity="0.0" + inkscape:pageshadow="2" + inkscape:zoom="2.8" + inkscape:cx="302.8754" + inkscape:cy="816.37285" + inkscape:document-units="px" + inkscape:current-layer="layer1" + showgrid="false" + inkscape:window-width="2560" + inkscape:window-height="1523" + inkscape:window-x="0" + inkscape:window-y="0" + inkscape:window-maximized="1" /> + <metadata + id="metadata7"> + <rdf:RDF> + <cc:Work + rdf:about=""> + <dc:format>image/svg+xml</dc:format> + <dc:type + rdf:resource="http://purl.org/dc/dcmitype/StillImage" /> + <dc:title></dc:title> + </cc:Work> + </rdf:RDF> + </metadata> + <g + inkscape:label="Layer 1" + inkscape:groupmode="layer" + id="layer1"> + <image + y="175.32289" + x="87.6194" + id="image3026" + xlink:href="file:///user/dsalinas/home/Documents/CodeSVN/gudhi_depot/trunk/src/Contraction/doc/SO3points.png" + height="107.55493" + width="121.70161" /> + <image + y="174.31145" + x="250.86069" + id="image3037" + xlink:href="file:///user/dsalinas/home/Documents/CodeSVN/gudhi_depot/trunk/src/Contraction/doc/SO3_rips.png" + height="107.95626" + width="121.70161" /> + <image + y="174.31216" + x="415.46198" + id="image3048" + xlink:href="file:///user/dsalinas/home/Documents/CodeSVN/gudhi_depot/trunk/src/Contraction/doc/SO3_simplified.png" + height="107.85593" + width="122.0026" /> + <rect + style="color:#000000;fill:none;stroke:#000000;stroke-width:0.10033109;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none;stroke-dashoffset:0;marker:none;visibility:visible;display:inline;overflow:visible;enable-background:accumulate" + id="rect3080-2" + width="27.49983" + height="26.862415" + x="160.74173" + y="176.24422" /> + <image + y="124.78581" + x="168.92697" + id="image3077" + xlink:href="file:///user/dsalinas/home/Documents/CodeSVN/gudhi_depot/trunk/src/Contraction/doc/zoom.png" + height="59.596668" + width="60.499645" /> + <rect + style="color:#000000;fill:none;stroke:#000000;stroke-width:0.20066218;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none;stroke-dashoffset:0;marker:none;visibility:visible;display:inline;overflow:visible;enable-background:accumulate" + id="rect3080" + width="61.215229" + height="59.796326" + x="168.33478" + y="124.91287" /> + <path + style="fill:none;stroke:#000000;stroke-width:0.20066218;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none;marker-end:url(#Arrow1Lend)" + d="m 163.59862,176.17779 c -2.13689,-4.88892 1.1683,-8.41755 4.60541,-9.6451" + id="path3879" + inkscape:connector-curvature="0" + sodipodi:nodetypes="cc" /> + <text + xml:space="preserve" + style="font-size:6.70418215px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans" + x="151.0036" + y="300.44409" + id="text4507" + sodipodi:linespacing="125%"><tspan + sodipodi:role="line" + id="tspan4509" + x="152.07077" + y="300.44409" + style="text-align:center;text-anchor:middle">Point cloud sampling SO3 </tspan><tspan + sodipodi:role="line" + x="151.0036" + y="308.82431" + id="tspan4513" + style="text-align:center;text-anchor:middle">(points are in R but projected into R</tspan><tspan + sodipodi:role="line" + x="151.0036" + y="317.20456" + id="tspan4515" + style="text-align:center;text-anchor:middle">for vizualization)</tspan><tspan + sodipodi:role="line" + x="151.0036" + y="325.58478" + id="tspan4511" + style="text-align:center;text-anchor:middle" /></text> + <text + xml:space="preserve" + style="font-size:4.58914995px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans" + x="141.32632" + y="304.69067" + id="text4517" + sodipodi:linespacing="125%"><tspan + sodipodi:role="line" + id="tspan4519" + x="141.32632" + y="304.69067">9 </tspan></text> + <text + xml:space="preserve" + style="font-size:4.58914995px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans" + x="210.88516" + y="304.76022" + id="text4521" + sodipodi:linespacing="125%"><tspan + sodipodi:role="line" + id="tspan4523" + x="210.88516" + y="304.76022">3</tspan></text> + <text + xml:space="preserve" + style="font-size:6.70418215px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans" + x="309.4176" + y="300.58682" + id="text4507-8" + sodipodi:linespacing="125%"><tspan + sodipodi:role="line" + x="309.4176" + y="300.58682" + id="tspan4515-4" + style="text-align:center;text-anchor:middle">Rips complex built uppon these points</tspan><tspan + sodipodi:role="line" + x="309.4176" + y="308.96704" + style="text-align:center;text-anchor:middle" + id="tspan4599">20 millions simplices</tspan><tspan + sodipodi:role="line" + x="309.4176" + y="317.34729" + id="tspan4511-3" + style="text-align:center;text-anchor:middle" /></text> + <text + xml:space="preserve" + style="font-size:6.70418215px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans" + x="476.61395" + y="300.4592" + id="text4507-8-0" + sodipodi:linespacing="125%"><tspan + sodipodi:role="line" + x="476.61395" + y="300.4592" + id="tspan4511-3-9" + style="text-align:center;text-anchor:middle">Simplicial complex obtained after simplification</tspan><tspan + sodipodi:role="line" + x="476.61395" + y="308.83942" + style="text-align:center;text-anchor:middle" + id="tspan4601">714 simplices</tspan></text> + </g> +</svg> diff --git a/src/Contraction/doc/sphere_contraction_representation.png b/src/Contraction/doc/sphere_contraction_representation.png Binary files differnew file mode 100644 index 00000000..edf37bf3 --- /dev/null +++ b/src/Contraction/doc/sphere_contraction_representation.png diff --git a/src/Contraction/doc/zoom.png b/src/Contraction/doc/zoom.png Binary files differnew file mode 100644 index 00000000..38d2b520 --- /dev/null +++ b/src/Contraction/doc/zoom.png diff --git a/src/Contraction/example/CMakeLists.txt b/src/Contraction/example/CMakeLists.txt new file mode 100644 index 00000000..f0dc885d --- /dev/null +++ b/src/Contraction/example/CMakeLists.txt @@ -0,0 +1,17 @@ +project(Contraction_examples) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable(RipsContraction Rips_contraction.cpp) + + add_executable(GarlandHeckbert Garland_heckbert.cpp) + target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY}) + + add_test(NAME Contraction_example_tore3D_0.2 COMMAND $<TARGET_FILE:RipsContraction> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "0.2") + # TODO(DS) : These tests are too long under Windows + #add_test(NAME Contraction_example_sphere_0.2 COMMAND $<TARGET_FILE:RipsContraction> + # "${CMAKE_SOURCE_DIR}/data/points/sphere3D_2646.off" "0.2") + #add_test(NAME Contraction_example_SO3_0.3 COMMAND $<TARGET_FILE:RipsContraction> + # "${CMAKE_SOURCE_DIR}/data/points/SO3_10000.off" "0.3") + +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Contraction/example/Garland_heckbert.cpp b/src/Contraction/example/Garland_heckbert.cpp new file mode 100644 index 00000000..9c0b5205 --- /dev/null +++ b/src/Contraction/example/Garland_heckbert.cpp @@ -0,0 +1,179 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + + +#ifndef GARLAND_HECKBERT_H_ +#define GARLAND_HECKBERT_H_ + +#include <gudhi/Point.h> +#include <gudhi/Edge_contraction.h> +#include <gudhi/Skeleton_blocker.h> +#include <gudhi/Off_reader.h> + +#include <iostream> + +#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<Geometry_trait> { + public: + struct Garland_heckbert_vertex : public Simple_geometric_vertex { + Error_quadric<Geometry_trait::Point> 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<Complex>; +using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor<Complex>; + +/** + * 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<EdgeProfile> { + Complex& complex_; + + public: + typedef Gudhi::contraction::Placement_policy<EdgeProfile>::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<Point> 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<EdgeProfile> { + Complex& complex_; + + public: + typedef Gudhi::contraction::Cost_policy<EdgeProfile>::Cost_type Cost_type; + + GH_cost(Complex& complex) : complex_(complex) { (void)complex_; } + + Cost_type operator()(EdgeProfile const& profile, boost::optional<Point> 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<EdgeProfile> { + 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<Point>(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<Complex> 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]); + + // 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<EdgeProfile>(), + 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<Complex> 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 new file mode 100644 index 00000000..49250d7a --- /dev/null +++ b/src/Contraction/example/Garland_heckbert/Error_quadric.h @@ -0,0 +1,169 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GARLAND_HECKBERT_ERROR_QUADRIC_H_ +#define GARLAND_HECKBERT_ERROR_QUADRIC_H_ + +#include <boost/optional/optional.hpp> + +#include <vector> +#include <utility> + +template <typename Point> 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<Point> min_cost(double scale = 1) const { + // const double min_determinant = 1e-4 * scale*scale; + const double min_determinant = 1e-5; + boost::optional<Point> 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/example/Rips_contraction.cpp b/src/Contraction/example/Rips_contraction.cpp new file mode 100644 index 00000000..b5ce06c1 --- /dev/null +++ b/src/Contraction/example/Rips_contraction.cpp @@ -0,0 +1,86 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ +#include <gudhi/Edge_contraction.h> +#include <gudhi/Skeleton_blocker.h> +#include <gudhi/Off_reader.h> +#include <gudhi/Point.h> +#include <gudhi/Clock.h> + +#include <iostream> + +struct Geometry_trait { + typedef Point_d Point; +}; + +using Complex_geometric_traits = Gudhi::skeleton_blocker::Skeleton_blocker_simple_geometric_traits<Geometry_trait>; +using Complex = Gudhi::skeleton_blocker::Skeleton_blocker_geometric_complex< Complex_geometric_traits >; +using Profile = Gudhi::contraction::Edge_profile<Complex>; +using Complex_contractor = Gudhi::contraction::Skeleton_blocker_contractor<Complex>; + + +template<typename ComplexType> +void build_rips(ComplexType& complex, double offset) { + if (offset <= 0) return; + auto vertices = complex.vertex_range(); + for (auto p = vertices.begin(); p != vertices.end(); ++p) + for (auto q = p; ++q != vertices.end(); /**/) { + if (squared_dist(complex.point(*p), complex.point(*q)) < 4 * offset * offset) + complex.add_edge_without_blockers(*p, *q); + } +} + +int main(int argc, char *argv[]) { + if (argc != 3) { + std::cerr << "Usage " << argv[0] << " ../../../data/meshes/SO3_10000.off 0.3 to load the file " << + "../../data/SO3_10000.off and contract the Rips complex built with paremeter 0.3.\n"; + return -1; + } + + Complex complex; + + // load only the points + Gudhi::skeleton_blocker::Skeleton_blocker_off_reader<Complex> off_reader(argv[1], complex, true); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file:" << argv[1] << std::endl; + return EXIT_FAILURE; + } + + std::cout << "Build the Rips complex with " << complex.num_vertices() << " vertices" << std::endl; + + build_rips(complex, atof(argv[2])); + + Gudhi::Clock contraction_chrono("Time to simplify and enumerate simplices"); + + std::cout << "Initial complex has " << + complex.num_vertices() << " vertices and " << + complex.num_edges() << " edges" << std::endl; + + Complex_contractor contractor(complex, + new Gudhi::contraction::Edge_length_cost<Profile>, + Gudhi::contraction::make_first_vertex_placement<Profile>(), + Gudhi::contraction::make_link_valid_contraction<Profile>(), + Gudhi::contraction::make_remove_popable_blockers_visitor<Profile>()); + contractor.contract_edges(); + + std::cout << "Counting final number of simplices \n"; + unsigned num_simplices = std::distance(complex.complex_simplex_range().begin(), complex.complex_simplex_range().end()); + + std::cout << "Final complex has " << + complex.num_vertices() << " vertices, " << + complex.num_edges() << " edges, " << + complex.num_blockers() << " blockers and " << + num_simplices << " simplices" << std::endl; + + std::cout << contraction_chrono; + + return EXIT_SUCCESS; +} + + diff --git a/src/Contraction/include/gudhi/Contraction/Edge_profile.h b/src/Contraction/include/gudhi/Contraction/Edge_profile.h new file mode 100644 index 00000000..0e914de9 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/Edge_profile.h @@ -0,0 +1,119 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_EDGE_PROFILE_H_ +#define CONTRACTION_EDGE_PROFILE_H_ + +#include <ostream> +#include <cassert> + +namespace Gudhi { + +namespace contraction { + +template<typename GeometricSimplifiableComplex> class Edge_profile { + public: + typedef GeometricSimplifiableComplex Complex; + typedef typename Complex::GT GT; + + typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle; + typedef typename GeometricSimplifiableComplex::Root_vertex_handle Root_vertex_handle; + + + typedef typename GeometricSimplifiableComplex::Edge_handle Edge_handle; + typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex; + typedef typename GeometricSimplifiableComplex::Graph_edge Graph_edge; + typedef typename GeometricSimplifiableComplex::Point Point; + + Edge_profile(GeometricSimplifiableComplex& complex, Edge_handle edge) : complex_(complex), edge_handle_(edge), + v0_(complex_.first_vertex(edge_handle_)), v1_(complex_.second_vertex(edge_handle_)) { + assert(complex_.get_address(complex_[edge_handle_].first())); + assert(complex_.get_address(complex_[edge_handle_].second())); + assert(complex_.contains_edge(v0_handle(), v1_handle())); + assert(v0_handle() != v1_handle()); + } + + virtual ~Edge_profile() { } + + GeometricSimplifiableComplex& complex() const { + return complex_; + } + + Edge_handle edge_handle() const { + return edge_handle_; + } + + Graph_edge& edge() const { + return complex_[edge_handle_]; + } + + Graph_vertex& v0() const { + return complex_[v0_handle()]; + } + + Graph_vertex& v1() const { + return complex_[v1_handle()]; + } + + Vertex_handle v0_handle() const { + return v0_; + // Root_vertex_handle root = complex_[edge_handle_].first(); + // assert(complex_.get_address(root)); + // return *complex_.get_address(root); + } + + Vertex_handle v1_handle() const { + return v1_; + // Root_vertex_handle root = complex_[edge_handle_].second(); + // assert(complex_.get_address(root)); + // return *complex_.get_address(root); + } + + const Point& p0() const { + return complex_.point(v0_handle()); + } + + const Point& p1() const { + return complex_.point(v1_handle()); + } + + friend std::ostream& operator<<(std::ostream& o, const Edge_profile& v) { + return o << "v0:" << v.v0_handle() << " v1:" << v.v1_handle(); + } + + private: + GeometricSimplifiableComplex& complex_; + + Edge_handle edge_handle_; + + Vertex_handle v0_; + + Vertex_handle v1_; +}; + +template<typename EdgeProfile> class Edge_profile_factory { + public: + typedef typename EdgeProfile::Edge_handle Edge_handle_; + typedef typename EdgeProfile::Complex Complex_; + + virtual EdgeProfile make_profile( + Complex_& complex, + Edge_handle_ edge) const { + return EdgeProfile(complex, edge); + } + + virtual ~Edge_profile_factory() { } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_EDGE_PROFILE_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Contraction_visitor.h b/src/Contraction/include/gudhi/Contraction/policies/Contraction_visitor.h new file mode 100644 index 00000000..243bc51c --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Contraction_visitor.h @@ -0,0 +1,79 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_CONTRACTION_VISITOR_H_ +#define CONTRACTION_POLICIES_CONTRACTION_VISITOR_H_ + +#include <gudhi/Contraction/Edge_profile.h> +#include <boost/optional.hpp> + +namespace Gudhi { + +namespace contraction { + +/** + *@class Contraction_visitor + *@brief Interface for a visitor of the edge contraction process. + *@ingroup contr + */ +template <typename EdgeProfile> +class Contraction_visitor { // : public Dummy_complex_visitor<typename EdgeProfile::Vertex_handle> { + public: + // typedef typename ComplexType::GeometryTrait GT; + typedef EdgeProfile Profile; + typedef double FT; + typedef typename Profile::Complex Complex; + typedef Complex ComplexType; + typedef typename ComplexType::Point Point; + + virtual ~Contraction_visitor() { } + + /** + * @brief Called before the edge contraction process starts. + */ + virtual void on_started(ComplexType & complex) { } + + /** + * @brief Called when the algorithm stops. + */ + virtual void on_stop_condition_reached() { } + + /** + * @brief Called during the collecting phase (when a cost is assigned to the edges), for each edge collected. + */ + virtual void on_collected(const Profile &profile, boost::optional< FT > cost) { } + + /** + * @brief Called during the processing phase (when edges are contracted), for each edge that is selected. + */ + virtual void on_selected(const Profile &profile, boost::optional< FT > cost, int initial_count, int current_count) { } + + /** + * @brief Called when an edge is about to be contracted and replaced by a vertex whose position is *placement. + */ + virtual void on_contracting(const Profile &profile, boost::optional< Point > placement) { } + + /** + * @brief Called when after an edge has been contracted onto a new point placement. + * A possibility would to remove popable blockers at this point for instance. + */ + virtual void on_contracted(const Profile &profile, boost::optional< Point > placement) { } + + /** + * @brief Called for each selected edge which cannot be contracted because the ValidContractionPredicate is false + */ + virtual void on_non_valid(const Profile &profile) { } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_CONTRACTION_VISITOR_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Cost_policy.h b/src/Contraction/include/gudhi/Contraction/policies/Cost_policy.h new file mode 100644 index 00000000..97114794 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Cost_policy.h @@ -0,0 +1,41 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_COST_POLICY_H_ +#define CONTRACTION_POLICIES_COST_POLICY_H_ + +#include <boost/optional.hpp> + +namespace Gudhi { + +namespace contraction { + +/** + *@brief Policy to specify the cost of contracting an edge. + *@ingroup contr + */ +template< typename EdgeProfile> +class Cost_policy { + public: + typedef typename EdgeProfile::Point Point; + typedef typename EdgeProfile::Graph_vertex Graph_vertex; + + typedef boost::optional<double> Cost_type; + + virtual Cost_type operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const = 0; + + virtual ~Cost_policy() { } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_COST_POLICY_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Dummy_valid_contraction.h b/src/Contraction/include/gudhi/Contraction/policies/Dummy_valid_contraction.h new file mode 100644 index 00000000..27a4dc7a --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Dummy_valid_contraction.h @@ -0,0 +1,37 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_DUMMY_VALID_CONTRACTION_H_ +#define CONTRACTION_POLICIES_DUMMY_VALID_CONTRACTION_H_ + +#include <gudhi/Contraction/policies/Valid_contraction_policy.h> + +namespace Gudhi { + +namespace contraction { + +/** + *@brief Policy that accept all edge contraction. + */ +template< typename EdgeProfile> +class Dummy_valid_contraction : public Valid_contraction_policy<EdgeProfile> { + public: + typedef typename EdgeProfile::Point Point; + + bool operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) { + return true; + } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_DUMMY_VALID_CONTRACTION_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Edge_length_cost.h b/src/Contraction/include/gudhi/Contraction/policies/Edge_length_cost.h new file mode 100644 index 00000000..97589385 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Edge_length_cost.h @@ -0,0 +1,44 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_EDGE_LENGTH_COST_H_ +#define CONTRACTION_POLICIES_EDGE_LENGTH_COST_H_ + +#include <gudhi/Contraction/policies/Cost_policy.h> + +namespace Gudhi { + +namespace contraction { + +/** + * @brief return a cost corresponding to the squared length of the edge + */ +template< typename EdgeProfile> +class Edge_length_cost : public Cost_policy<EdgeProfile> { + public: + typedef typename Cost_policy<EdgeProfile>::Cost_type Cost_type; + typedef typename EdgeProfile::Point Point; + + Cost_type operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const override { + double res = 0; + auto p0_coord = profile.p0().begin(); + auto p1_coord = profile.p1().begin(); + for (; p0_coord != profile.p0().end(); p0_coord++, p1_coord++) { + res += (*p0_coord - *p1_coord) * (*p0_coord - *p1_coord); + } + return res; + } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_EDGE_LENGTH_COST_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/First_vertex_placement.h b/src/Contraction/include/gudhi/Contraction/policies/First_vertex_placement.h new file mode 100644 index 00000000..005b80e0 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/First_vertex_placement.h @@ -0,0 +1,40 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_FIRST_VERTEX_PLACEMENT_H_ +#define CONTRACTION_POLICIES_FIRST_VERTEX_PLACEMENT_H_ + +#include <gudhi/Contraction/policies/Placement_policy.h> + +namespace Gudhi { + +namespace contraction { + +/** + * @brief Places the contracted point onto the first point of the edge + */ +template< typename EdgeProfile> +class First_vertex_placement : public Placement_policy<EdgeProfile> { + public: + typedef typename EdgeProfile::Point Point; + typedef typename EdgeProfile::Edge_handle Edge_handle; + + typedef typename Placement_policy<EdgeProfile>::Placement_type Placement_type; + + Placement_type operator()(const EdgeProfile& profile) const override { + return Placement_type(profile.p0()); + } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_FIRST_VERTEX_PLACEMENT_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h b/src/Contraction/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h new file mode 100644 index 00000000..2e7ea481 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Link_condition_valid_contraction.h @@ -0,0 +1,44 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_LINK_CONDITION_VALID_CONTRACTION_H_ +#define CONTRACTION_POLICIES_LINK_CONDITION_VALID_CONTRACTION_H_ + +#include <gudhi/Contraction/policies/Valid_contraction_policy.h> +#include <gudhi/Debug_utils.h> + + +namespace Gudhi { + +namespace contraction { + +/** + *@brief Policy that only accept edges verifying the link condition (and therefore whose contraction preserving homotopy type). + *@ingroup contr + */ +template< typename EdgeProfile> +class Link_condition_valid_contraction : public Valid_contraction_policy<EdgeProfile> { + public: + typedef typename EdgeProfile::Edge_handle Edge_handle; + typedef typename EdgeProfile::Point Point; + // typedef typename EdgeProfile::Edge_handle Edge_handle; + + bool operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const override { + Edge_handle edge(profile.edge_handle()); + DBGMSG("Link_condition_valid_contraction:", profile.complex().link_condition(edge)); + return profile.complex().link_condition(edge); + } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_LINK_CONDITION_VALID_CONTRACTION_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Middle_placement.h b/src/Contraction/include/gudhi/Contraction/policies/Middle_placement.h new file mode 100644 index 00000000..7dcf708b --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Middle_placement.h @@ -0,0 +1,39 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_MIDDLE_PLACEMENT_H_ +#define CONTRACTION_POLICIES_MIDDLE_PLACEMENT_H_ + +#include <gudhi/Contraction/policies/Placement_policy.h> + +namespace Gudhi { + +namespace contraction { + +template< typename EdgeProfile> +class Middle_placement : public Placement_policy<EdgeProfile> { + public: + typedef typename EdgeProfile::Point Point; + typedef typename EdgeProfile::Edge_handle Edge_handle; + typedef typename EdgeProfile::Graph_vertex Graph_vertex; + + typedef typename Placement_policy<EdgeProfile>::Placement_type Placement_type; + + Placement_type operator()(const EdgeProfile& profile) const override { + // todo compute the middle + return Placement_type(profile.p0()); + } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_MIDDLE_PLACEMENT_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Placement_policy.h b/src/Contraction/include/gudhi/Contraction/policies/Placement_policy.h new file mode 100644 index 00000000..5f97d6a7 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Placement_policy.h @@ -0,0 +1,39 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_PLACEMENT_POLICY_H_ +#define CONTRACTION_POLICIES_PLACEMENT_POLICY_H_ + +#include <boost/optional.hpp> + +namespace Gudhi { + +namespace contraction { + +/** + *@brief Policy to specify where the merged point had to be placed after an edge contraction. + *@ingroup contr + */ +template< typename EdgeProfile> +class Placement_policy { + public: + typedef typename EdgeProfile::Point Point; + typedef boost::optional<Point> Placement_type; + + virtual Placement_type operator()(const EdgeProfile& profile) const = 0; + + virtual ~Placement_policy() { } +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // CONTRACTION_POLICIES_PLACEMENT_POLICY_H_ diff --git a/src/Contraction/include/gudhi/Contraction/policies/Valid_contraction_policy.h b/src/Contraction/include/gudhi/Contraction/policies/Valid_contraction_policy.h new file mode 100644 index 00000000..413c5bd6 --- /dev/null +++ b/src/Contraction/include/gudhi/Contraction/policies/Valid_contraction_policy.h @@ -0,0 +1,39 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONTRACTION_POLICIES_VALID_CONTRACTION_POLICY_H_ +#define CONTRACTION_POLICIES_VALID_CONTRACTION_POLICY_H_ + +namespace Gudhi { + +namespace contraction { + +/** + *@brief Policy to specify if an edge contraction is valid or not. + *@ingroup contr + */ +template< typename EdgeProfile> +class Valid_contraction_policy { + public: + typedef typename EdgeProfile::Point Point; + typedef typename EdgeProfile::Edge_handle Edge_handle; + typedef typename EdgeProfile::Graph_vertex Graph_vertex; + + virtual bool operator()(const EdgeProfile& profile, const boost::optional<Point>& placement) const = 0; + + virtual ~Valid_contraction_policy() { } +}; + +} // namespace contraction + +} // namespace Gudhi + + +#endif // CONTRACTION_POLICIES_VALID_CONTRACTION_POLICY_H_ diff --git a/src/Contraction/include/gudhi/Edge_contraction.h b/src/Contraction/include/gudhi/Edge_contraction.h new file mode 100644 index 00000000..6058d64b --- /dev/null +++ b/src/Contraction/include/gudhi/Edge_contraction.h @@ -0,0 +1,219 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef EDGE_CONTRACTION_H_ +#define EDGE_CONTRACTION_H_ + + +#include <gudhi/Skeleton_blocker_contractor.h> +#include <gudhi/Contraction/policies/Edge_length_cost.h> +#include <gudhi/Contraction/policies/First_vertex_placement.h> +#include <gudhi/Contraction/policies/Valid_contraction_policy.h> +#include <gudhi/Contraction/policies/Dummy_valid_contraction.h> +#include <gudhi/Contraction/policies/Link_condition_valid_contraction.h> +#include <gudhi/Debug_utils.h> + +namespace Gudhi { + +namespace contraction { + + +/** \defgroup contr Edge contraction + +\author David Salinas + +\section edgecontractionintroduction Introduction + +The purpose of this package is to offer a user-friendly interface for edge contraction simplification of huge simplicial complexes. +It uses the \ref skbl data-structure whose size remains small during simplification +of most used geometrical complexes of topological data analysis such as the Rips or the Delaunay complexes. In practice, the +size of this data-structure is even much lower than the total number of simplices. + +The edge contraction operation consists in identifying two vertices of a simplicial complex. +A lot of algorithms have been developed in computer graphics that allows to reduce efficiently the size of 2-dimensional simplicial complexes + while preserving its geometry (for instance see \cite Garland \cite Lindstrom). +These approaches can be extended to higher-dimensional simplicial complexes. +The main advantage of using the Skeleton-Blocker data structure for edge contraction is that when the number of blockers is small, +the operations needed for edge contraction algorithms have polynomial complexity regarding the size the graph. +Therefore, the simplification can be done without enumerating the set of simplices that is often non tracktable in high-dimension and is then very efficient +(sub-linear with regards to the number of simplices in practice). + +A typical application of this package is homology group computation. It is illustrated in the next figure where a Rips complex is built uppon a set of high-dimensional points and +simplified with edge contractions. +It has initially a big number of simplices (around 20 millions) but simplifying it to a much reduced form with only 15 vertices (and 714 simplices) takes only few seconds on a desktop machine (see the example bellow). +One can then compute homology group with a simplicial complex having very few simplices instead of running the homology algorithm on the much bigger initial set of +simplices which would take much more time and memory. + + +\image html "so3.png" "Edge contraction illustration" + +\section Design + +This class design is policy based and heavily inspired from the similar edge collapse package of CGAL http://doc.cgal.org/latest/Surface_mesh_simplification/index.html (which is however restricted to 2D triangulations). + + +\subsection Policies + +Four policies can be customized in this package: + +\li Cost_policy: specify how much cost an edge contraction of a given edge. The edge with lowest cost is iteratively picked and contracted if valid. +\li Valid_contraction_policy: specify if a given edge contraction is valid. For instance, this policy can check the link condition which ensures that the homotopy type is preserved afer the edge contraction. +\li Placement_policy: every time an edge is contracted, its points are merge to one point specified by this policy. This may be the middle of the edge of some more sophisticated point such as the minimum of a cost as in +\cite Garland. + + +\subsection Visitor + +A visitor which implements the class Contraction_visitor gets called at several key moments +during the simplification: + +\li when the algorithms starts or ends, +\li when an edge is seen for the first time, +\li when an edge is considered for an edge contraction, +\li before and after an edge is contracted. + +This allows to implements most of edge contraction based algorithm with this +package without having to change the main simplification source code. + + +\section Performance + +The next figure shows the computation time to reduce a random 2-sphere to a single tetrahedron with +this package and with the CGAL equivalent package. +Despite this package is able to deal with \a arbitrary simplicial complexes (any dimensions, manifold or non manifold), +it is still \a 65% times faster than the CGAL package which is focused on 2-manifold. +The main reason is that few blockers appears during the simplification and hence, +the algorithm only have to deal with the graph and not higher-dimensional simplices +(in this case triangles). However, we recall that higher-dimensional simplices are \a implicitely +stored in the \ref skbl data-structure. Hence, one has to store +simplices in an external map if some information needs to be associated with them (information that could be a filtration value or +an orientation for instance). + + +\image html "sphere_contraction.png" "Time in seconds to simplify random 2-spheres to a tetrahedron" width=10cm + +\section Example + + +This example loads points from an OFF file and builds the Rips complex with an user provided parameter. It then simplifies the built Rips complex +while ensuring its homotopy type is preserved during the contraction (edge are contracted only when the link condition is valid). + + \code{.cpp} +#include <boost/timer/timer.hpp> +#include <iostream> +#include <gudhi/Edge_contraction.h> +#include <gudhi/Skeleton_blocker.h> +#include <gudhi/Off_reader.h> + + +using namespace std; +using namespace Gudhi; +using namespace skeleton_blocker; +using namespace contraction; + +struct Geometry_trait{ + typedef std::vector<double> Point; +}; + +typedef Geometry_trait::Point Point; +typedef Skeleton_blocker_simple_geometric_traits<Geometry_trait> Complex_geometric_traits; +typedef Skeleton_blocker_geometric_complex< Complex_geometric_traits > Complex; +typedef Edge_profile<Complex> Profile; +typedef Skeleton_blocker_contractor<Complex> Complex_contractor; + +template<typename Point> +double eucl_distance(const Point& a,const Point& b){ + double res = 0; + auto a_coord = a.begin(); + auto b_coord = b.begin(); + for(; a_coord != a.end(); a_coord++, b_coord++){ + res += (*a_coord - *b_coord) * (*a_coord - *b_coord); + } + return sqrt(res); +} + +template<typename ComplexType> +void build_rips(ComplexType& complex, double offset){ + if (offset<=0) return; + auto vertices = complex.vertex_range(); + for (auto p = vertices.begin(); p != vertices.end(); ++p) + for (auto q = p; ++q != vertices.end(); ) + if (eucl_distance(complex.point(*p), complex.point(*q)) < 2 * offset) + complex.add_edge_without_blockers(*p,*q); +} + +int main (int argc, char *argv[]) +{ + if (argc!=3){ + std::cerr << "Usage "<<argv[0]<<" ../../data/SO3_10000.off 0.3 to load the file ../../data/SO3_10000.off and contract the Rips complex built with paremeter 0.3.\n"; + return -1; + } + + Complex complex; + + // load the points + Skeleton_blocker_off_reader<Complex> off_reader(argv[1],complex,true); + if(!off_reader.is_valid()){ + std::cerr << "Unable to read file:"<<argv[1]<<std::endl; + return EXIT_FAILURE; + } + std::cout << "Build the Rips complex"<<std::endl; + + build_rips(complex,atof(argv[2])); + + boost::timer::auto_cpu_timer t; + + std::cout << "Initial complex has "<< + complex.num_vertices()<<" vertices and "<< + complex.num_edges()<<" edges"<<std::endl; + + Complex_contractor contractor(complex, + new Edge_length_cost<Profile>, + contraction::make_first_vertex_placement<Profile>(), + contraction::make_link_valid_contraction<Profile>(), + contraction::make_remove_popable_blockers_visitor<Profile>()); + contractor.contract_edges(); + + std::cout << "Counting final number of simplices \n"; + unsigned num_simplices = std::distance(complex.star_simplex_range().begin(),complex.star_simplex_range().end()); + + std::cout << "Final complex has "<< + complex.num_vertices()<<" vertices, "<< + complex.num_edges()<<" edges, "<< + complex.num_blockers()<<" blockers and "<< + num_simplices<<" simplices"<<std::endl; + + + std::cout << "Time to simplify and enumerate simplices:\n"; + + return EXIT_SUCCESS; +} +} + \endcode + +\verbatim +./example/Contraction/RipsContraction ../../data/SO3_10000.off 0.3 +[ 50%] [100%] Built target SkeletonBlockerIteration +Built target RipsContraction +Build the Rips complex +Initial complex has 10000 vertices and 195664 edges +Counting final number of simplices +Final complex has 15 vertices, 101 edges, 165 blockers and 714 simplices +Time to simplify and enumerate simplices: + 3.166621s wall, 3.150000s user + 0.010000s system = 3.160000s CPU (99.8%) +\endverbatim + +*/ +/** @} */ // end defgroup +} // namespace contraction + +} // namespace Gudhi + +#endif // EDGE_CONTRACTION_H_ diff --git a/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h new file mode 100644 index 00000000..c2b3157c --- /dev/null +++ b/src/Contraction/include/gudhi/Skeleton_blocker_contractor.h @@ -0,0 +1,576 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): David Salinas + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SKELETON_BLOCKER_CONTRACTOR_H_ +#define SKELETON_BLOCKER_CONTRACTOR_H_ + +#include <gudhi/Contraction/Edge_profile.h> +#include <gudhi/Contraction/policies/Cost_policy.h> +#include <gudhi/Contraction/policies/Edge_length_cost.h> +#include <gudhi/Contraction/policies/Placement_policy.h> +#include <gudhi/Contraction/policies/First_vertex_placement.h> +#include <gudhi/Contraction/policies/Valid_contraction_policy.h> +#include <gudhi/Contraction/policies/Dummy_valid_contraction.h> // xxx remove +#include <gudhi/Contraction/policies/Link_condition_valid_contraction.h> // xxx remove +#include <gudhi/Contraction/policies/Contraction_visitor.h> + +#include <gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h> +#include <gudhi/Debug_utils.h> + +// todo remove the queue to be independent from cgal +#include <CGAL/Modifiable_priority_queue.h> +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <boost/scoped_array.hpp> +#include <boost/scoped_ptr.hpp> + +#include <memory> +#include <cassert> +#include <list> +#include <utility> // for pair +#include <vector> + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +namespace Gudhi { + +namespace contraction { + +template <class Profile> +Placement_policy<Profile>* make_first_vertex_placement() { + return new First_vertex_placement<Profile>(); +} + +template <class Profile> +Valid_contraction_policy<Profile>* make_link_valid_contraction() { + return new Link_condition_valid_contraction<Profile>(); +} + +/** + *@brief Visitor to remove popable blockers after an edge contraction. + */ +template <class Profile> +class Contraction_visitor_remove_popable : public Contraction_visitor<Profile> { + public: + typedef typename Profile::Point Point; + typedef typename Profile::Complex Complex; + typedef typename Complex::Vertex_handle Vertex_handle; + + void on_contracted(const Profile &profile, boost::optional< Point > placement) override { + profile.complex().remove_all_popable_blockers(profile.v0_handle()); + } +}; + +template <class Profile> +Contraction_visitor<Profile>* make_remove_popable_blockers_visitor() { + return new Contraction_visitor_remove_popable<Profile>(); +} + +/** + *@class Skeleton_blocker_contractor + *@brief Class that allows to contract iteratively edges of a simplicial complex. + *@ingroup contr + * + * @details The simplification algorithm consists in iteratively picking the + * edge with lowest cost and performing an edge contraction if the contraction is valid. + * This class is policy based (and much inspired from the edge collapse package of CGAL http://doc.cgal.org/latest/Surface_mesh_simplification/index.html). + * + * Policies that can be changed are : + * - the cost policy : how much cost an edge contraction + * - the placement policy : where will be placed the contraction point + * - the valid contraction policy : is the contraction valid. For instance, it can be + * a topological condition (link condition) or a geometrical condition (normals messed up). + * + */ +template<class GeometricSimplifiableComplex, class EdgeProfile = Edge_profile<GeometricSimplifiableComplex>> +class Skeleton_blocker_contractor : private skeleton_blocker::Dummy_complex_visitor< +typename GeometricSimplifiableComplex::Vertex_handle> { + GeometricSimplifiableComplex& complex_; + + public: + typedef typename GeometricSimplifiableComplex::Graph_vertex Graph_vertex; + typedef typename GeometricSimplifiableComplex::Vertex_handle Vertex_handle; + typedef typename GeometricSimplifiableComplex::Simplex Simplex; + typedef typename GeometricSimplifiableComplex::Root_vertex_handle Root_vertex_handle; + typedef typename GeometricSimplifiableComplex::Graph_edge Graph_edge; + typedef typename GeometricSimplifiableComplex::Edge_handle Edge_handle; + typedef typename GeometricSimplifiableComplex::Point Point; + + typedef EdgeProfile Profile; + + + typedef Cost_policy<Profile> Cost_policy_; + typedef Placement_policy<Profile> Placement_policy_; + typedef Valid_contraction_policy<Profile> Valid_contraction_policy_; + typedef Contraction_visitor<EdgeProfile> Contraction_visitor_; + typedef Edge_profile_factory<EdgeProfile> Edge_profile_factory_; + + + + typedef boost::optional<double> Cost_type; + typedef boost::optional<Point> Placement_type; + + typedef size_t size_type; + + typedef Skeleton_blocker_contractor Self; + + private: + struct Compare_id { + Compare_id() : algorithm_(0) { } + + Compare_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } + + bool operator()(Edge_handle a, Edge_handle b) const { + return algorithm_->get_undirected_edge_id(a) < algorithm_->get_undirected_edge_id(b); + } + + Self const* algorithm_; + }; + + struct Compare_cost { + Compare_cost() : algorithm_(0) { } + + Compare_cost(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } + + bool operator()(Edge_handle a, Edge_handle b) const { + // NOTE: A cost is an optional<> value. + // Absent optionals are ordered first; that is, "none < T" and "T > none" for any defined T != none. + // In consequence, edges with undefined costs will be promoted to the top of the priority queue and popped out + // first. + return algorithm_->get_data(a).cost() < algorithm_->get_data(b).cost(); + } + + Self const* algorithm_; + }; + + struct Undirected_edge_id : boost::put_get_helper<size_type, Undirected_edge_id> { + typedef boost::readable_property_map_tag category; + typedef size_type value_type; + typedef size_type reference; + typedef Edge_handle key_type; + + Undirected_edge_id() : algorithm_(0) { } + + Undirected_edge_id(Self const* aAlgorithm) : algorithm_(aAlgorithm) { } + + size_type operator[](Edge_handle e) const { + return algorithm_->get_undirected_edge_id(e); + } + + Self const* algorithm_; + }; + + typedef CGAL::Modifiable_priority_queue<Edge_handle, Compare_cost, Undirected_edge_id> PQ; + typedef typename PQ::handle pq_handle; + + + // An Edge_data is associated with EVERY edge in the complex (collapsible or not). + // It relates the edge with the PQ-handle needed to update the priority queue + // It also relates the edge with a policy-based cache + + class Edge_data { + public: + Edge_data() : PQHandle_(), cost_() { } + + Cost_type const& cost() const { + return cost_; + } + + Cost_type & cost() { + return cost_; + } + + pq_handle PQ_handle() const { + return PQHandle_; + } + + bool is_in_PQ() const { + return PQHandle_ != PQ::null_handle(); + } + + void set_PQ_handle(pq_handle h) { + PQHandle_ = h; + } + + void reset_PQ_handle() { + PQHandle_ = PQ::null_handle(); + } + + private: + pq_handle PQHandle_; + Cost_type cost_; + }; + typedef Edge_data* Edge_data_ptr; + typedef boost::scoped_array<Edge_data> Edge_data_array; + + int get_undirected_edge_id(Edge_handle edge) const { + return complex_[edge].index(); + } + + const Edge_data& get_data(Edge_handle edge) const { + return edge_data_array_[get_undirected_edge_id(edge)]; + } + + Edge_data& get_data(Edge_handle edge) { + return edge_data_array_[get_undirected_edge_id(edge)]; + } + + Cost_type get_cost(const Profile & profile) const { + return (*cost_policy_)(profile, get_placement(profile)); + } + + Profile create_profile(Edge_handle edge) const { + if (edge_profile_factory_) + return edge_profile_factory_->make_profile(complex_, edge); + else + return Profile(complex_, edge); + } + + void insert_in_PQ(Edge_handle edge, Edge_data& data) { + data.set_PQ_handle(heap_PQ_->push(edge)); + ++current_num_edges_heap_; + } + + void update_in_PQ(Edge_handle edge, Edge_data& data) { + data.set_PQ_handle(heap_PQ_->update(edge, data.PQ_handle())); + } + + void remove_from_PQ(Edge_handle edge, Edge_data& data) { + data.set_PQ_handle(heap_PQ_->erase(edge, data.PQ_handle())); + --current_num_edges_heap_; + } + + boost::optional<Edge_handle> pop_from_PQ() { + boost::optional<Edge_handle> edge = heap_PQ_->extract_top(); + if (edge) + get_data(*edge).reset_PQ_handle(); + return edge; + } + + private: + /** + * @brief Collect edges. + * + * Iterates over all edges of the simplicial complex and + * 1) inserts them in the priority queue sorted according to the Cost policy. + * 2) set the id() field of each edge + */ + void collect_edges() { + // + // Loop over all the edges in the complex in the heap + // + size_type size = complex_.num_edges(); + DBG("Collecting edges ..."); + DBGMSG("num edges :", size); + + edge_data_array_.reset(new Edge_data[size]); + + heap_PQ_.reset(new PQ(size, Compare_cost(this), Undirected_edge_id(this))); + + std::size_t id = 0; + + // xxx do a parralel for + for (auto edge : complex_.edge_range()) { + complex_[edge].index() = id++; + Profile const& profile = create_profile(edge); + Edge_data& data = get_data(edge); + data.cost() = get_cost(profile); + ++initial_num_edges_heap_; + insert_in_PQ(edge, data); + if (contraction_visitor_) contraction_visitor_->on_collected(profile, data.cost()); + } + + DBG("Edges collected."); + } + + bool should_stop(double lCost, const Profile &profile) const { + return false; + } + + boost::optional<Point> get_placement(const Profile& profile) const { + return (*placement_policy_)(profile); + } + + bool is_contraction_valid(Profile const& profile, Placement_type placement) const { + if (!valid_contraction_policy_) return true; + return (*valid_contraction_policy_)(profile, placement); + } + + + public: + /** + * \brief Contract edges. + * + * While the heap is not empty, it extracts the edge with the minimum + * cost in the heap then try to contract it. + * It stops when the Stop policy says so or when the number of contractions + * given by 'num_max_contractions' is reached (if this number is positive). + */ + void contract_edges(int num_max_contractions = -1) { + DBG("\n\nContract edges"); + int num_contraction = 0; + + bool unspecified_num_contractions = (num_max_contractions == -1); + // + // Pops and processes each edge from the PQ + // + boost::optional<Edge_handle> edge; + while ((edge = pop_from_PQ()) && ((num_contraction < num_max_contractions) || (unspecified_num_contractions))) { + Profile const& profile = create_profile(*edge); + Cost_type cost(get_data(*edge).cost()); + if (contraction_visitor_) contraction_visitor_->on_selected(profile, cost, 0, 0); + + DBGMSG("\n\n---- Pop edge - num vertices :", complex_.num_vertices()); + + if (cost) { + DBGMSG("sqrt(cost):", std::sqrt(*cost)); + if (should_stop(*cost, profile)) { + if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached(); + DBG("should_stop"); + break; + } + Placement_type placement = get_placement(profile); + if (is_contraction_valid(profile, placement) && placement) { + DBG("contraction_valid"); + contract_edge(profile, placement); + ++num_contraction; + } else { + DBG("contraction not valid"); + if (contraction_visitor_) contraction_visitor_->on_non_valid(profile); + } + } else { + DBG("uncomputable cost"); + } + } + if (contraction_visitor_) contraction_visitor_->on_stop_condition_reached(); + } + + bool is_in_heap(Edge_handle edge) const { + if (heap_PQ_->empty()) { + return false; + } else { + return edge_data_array_[get_undirected_edge_id(edge)].is_in_PQ(); + } + } + + bool is_heap_empty() const { + return heap_PQ_->empty(); + } + + /** + * @brief Returns an Edge_handle and a Placement_type. This pair consists in + * the edge with the lowest cost in the heap together with its placement. + * The returned value is initialized iff the heap is non-empty. + */ + boost::optional<std::pair<Edge_handle, Placement_type > > top_edge() { + boost::optional<std::pair<Edge_handle, Placement_type > > res; + + if (!heap_PQ_->empty()) { + auto edge = heap_PQ_->top(); + Profile const& profile = create_profile(edge); + Placement_type placement = get_placement(profile); + res = std::make_pair(edge, placement); + DBGMSG("top edge:", complex_[edge]); + } + return res; + } + + /** + * @brief Constructor with default policies. + * + * @details The default cost, placement, valid and visitor policies + * are respectively : the edge length, the first point, the link condition + */ + Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex) + : complex_(complex), + cost_policy_(new Edge_length_cost<Profile>), + placement_policy_(new First_vertex_placement<Profile>), + valid_contraction_policy_(new Link_condition_valid_contraction<Profile>), + contraction_visitor_(new Contraction_visitor_()), + edge_profile_factory_(0), + initial_num_edges_heap_(0), + current_num_edges_heap_(0) { + complex_.set_visitor(this); + if (contraction_visitor_) contraction_visitor_->on_started(complex); + collect_edges(); + } + + /** + * @brief Constructor with customed policies. + * @remark Policies destruction is handle by the class with smart pointers. + */ + Skeleton_blocker_contractor(GeometricSimplifiableComplex& complex, + Cost_policy_ *cost_policy, + Placement_policy_ * placement_policy = new First_vertex_placement<Profile>, + Valid_contraction_policy_ * valid_contraction_policy = + new Link_condition_valid_contraction<Profile>, + Contraction_visitor_* contraction_visitor = new Contraction_visitor_(), + Edge_profile_factory_* edge_profile_factory = NULL) : + complex_(complex), + cost_policy_(cost_policy), + placement_policy_(placement_policy), + valid_contraction_policy_(valid_contraction_policy), + contraction_visitor_(contraction_visitor), + edge_profile_factory_(edge_profile_factory), + initial_num_edges_heap_(0), + current_num_edges_heap_(0) { + complex_.set_visitor(this); + if (contraction_visitor) contraction_visitor->on_started(complex); + collect_edges(); + } + + ~Skeleton_blocker_contractor() { + complex_.set_visitor(0); + } + + private: + void contract_edge(const Profile& profile, Placement_type placement) { + if (contraction_visitor_) contraction_visitor_->on_contracting(profile, placement); + + assert(complex_.contains_vertex(profile.v0_handle())); + assert(complex_.contains_vertex(profile.v1_handle())); + assert(placement); + + profile.complex().point(profile.v0_handle()) = *placement; + + // remark : this is not necessary since v1 will be deactivated + // profile.complex().point(profile.v1_handle()) = *placement; + + complex_.contract_edge(profile.v0_handle(), profile.v1_handle()); + + assert(complex_.contains_vertex(profile.v0_handle())); + assert(!complex_.contains_vertex(profile.v1_handle())); + + update_changed_edges(); + + // the visitor could do something as complex_.remove_popable_blockers(); + if (contraction_visitor_) contraction_visitor_->on_contracted(profile, placement); + } + + private: + // every time the visitor's method on_changed_edge is called, it adds an + // edge to changed_edges_ + std::vector< Edge_handle > changed_edges_; + + /** + * @brief we update the cost and the position in the heap of an edge that has + * been changed + */ + inline void on_changed_edge(Vertex_handle a, Vertex_handle b) override { + boost::optional<Edge_handle> ab(complex_[std::make_pair(a, b)]); + assert(ab); + changed_edges_.push_back(*ab); + } + + void update_changed_edges() { + // xxx do a parralel for + DBG("update edges"); + + // sequential loop + for (auto ab : changed_edges_) { + // 1-get the Edge_handle corresponding to ab + // 2-change the data in mEdgeArray[ab.id()] + // 3-update the heap + Edge_data& data = get_data(ab); + Profile const& profile = create_profile(ab); + data.cost() = get_cost(profile); + if (data.is_in_PQ()) { + update_in_PQ(ab, data); + } else { + insert_in_PQ(ab, data); + } + } + changed_edges_.clear(); + } + + + private: + void on_remove_edge(Vertex_handle a, Vertex_handle b) override { + boost::optional<Edge_handle> ab((complex_[std::make_pair(a, b)])); + assert(ab); + Edge_data& lData = get_data(*ab); + if (lData.is_in_PQ()) { + remove_from_PQ(*ab, lData); + } + } + + private: + /** + * @brief Called when the edge 'ax' has been added while the edge 'bx' + * is still there but will be removed on next instruction. + * We assign the index of 'bx' to the edge index of 'ax' + */ + void on_swaped_edge(Vertex_handle a, Vertex_handle b, Vertex_handle x) override { + boost::optional<Edge_handle> ax(complex_[std::make_pair(a, x)]); + boost::optional<Edge_handle> bx(complex_[std::make_pair(b, x)]); + assert(ax && bx); + complex_[*ax].index() = complex_[*bx].index(); + } + + private: + /** + * @brief Called when a blocker is removed. + * All the edges that passes through the blocker may be edge-contractible + * again and are thus reinserted in the heap. + */ + void on_delete_blocker(const Simplex * blocker) override { + // we go for all pairs xy that belongs to the blocker + // note that such pairs xy are necessarily edges of the complex + // by definition of a blocker + + // todo uniqument utile pour la link condition + // laisser a l'utilisateur ? booleen update_heap_on_removed_blocker? + Simplex blocker_copy(*blocker); + for (auto x = blocker_copy.begin(); x != blocker_copy.end(); ++x) { + for (auto y = x; ++y != blocker_copy.end();) { + auto edge_descr(complex_[std::make_pair(*x, *y)]); + assert(edge_descr); + Edge_data& data = get_data(*edge_descr); + Profile const& profile = create_profile(*edge_descr); + data.cost() = get_cost(profile); + + // If the edge is already in the heap + // its priority has not changed. + // If the edge is not present, we reinsert it + // remark : we could also reinsert the edge + // only if it is valid + if (!data.is_in_PQ()) { + insert_in_PQ(*edge_descr, data); + } + } + } + } + + + private: + std::shared_ptr<Cost_policy_> cost_policy_; + std::shared_ptr<Placement_policy_> placement_policy_; + std::shared_ptr<Valid_contraction_policy_> valid_contraction_policy_; + std::shared_ptr<Contraction_visitor_> contraction_visitor_; + + // in case the user wants to do something special when the edge profile + // are created (for instance add some info) + std::shared_ptr<Edge_profile_factory_> edge_profile_factory_; + Edge_data_array edge_data_array_; + + boost::scoped_ptr<PQ> heap_PQ_; + int initial_num_edges_heap_; + int current_num_edges_heap_; +}; + +} // namespace contraction + +} // namespace Gudhi + +#endif // SKELETON_BLOCKER_CONTRACTOR_H_ |
