diff --git a/BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h b/BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h index 308da82f6d0..ecf70d29203 100644 --- a/BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h +++ b/BGL/include/CGAL/boost/graph/alpha_expansion_graphcut.h @@ -647,18 +647,143 @@ double alpha_expansion_graphcut (const InputGraph& input_graph, std::size_t vertex_i = get (vertex_index_map, vd); alpha_expansion.update(vertex_label_map, inserted_vertices, vd, vertex_i, alpha); } - } - } while(success); + + template + double min_cut(const InputGraph& input_graph, + EdgeCostMap edge_cost_map, + VertexLabelCostMap vertex_label_cost_map, + VertexLabelMap vertex_label_map, + const NamedParameters& np) + { + using parameters::choose_parameter; + using parameters::get_parameter; + + typedef boost::graph_traits GT; + typedef typename GT::edge_descriptor input_edge_descriptor; + typedef typename GT::vertex_descriptor input_vertex_descriptor; + + typedef typename GetInitializedVertexIndexMap::type VertexIndexMap; + VertexIndexMap vertex_index_map = CGAL::get_initialized_vertex_index_map(input_graph, np); + + typedef typename GetImplementationTag::type Impl_tag; + + // select implementation + typedef typename std::conditional + ::value, + Alpha_expansion_boost_adjacency_list_impl, + typename std::conditional + ::value, + Alpha_expansion_boost_compressed_sparse_row_impl, + Alpha_expansion_MaxFlow_impl>::type>::type + Alpha_expansion; + + typedef typename Alpha_expansion::Vertex_descriptor Vertex_descriptor; + + Alpha_expansion graph; + + // TODO: check this hardcoded parameter + const double tolerance = 1e-10; + + double min_cut = (std::numeric_limits::max)(); #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT - CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << - std::endl; - CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; - CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; + double vertex_creation_time, edge_creation_time, cut_time; + vertex_creation_time = edge_creation_time = cut_time = 0.0; #endif - return min_cut; -} + std::vector inserted_vertices; + inserted_vertices.resize(num_vertices(input_graph)); + + std::size_t number_of_labels = get(vertex_label_cost_map, *(vertices(input_graph).first)).size(); + + graph.clear_graph(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + Timer timer; + timer.start(); +#endif + + // For E-Data + // add every input vertex as a vertex to the graph, put edges to source & sink vertices + for (input_vertex_descriptor vd : CGAL::make_range(vertices(input_graph))) + { + std::size_t vertex_i = get(vertex_index_map, vd); + Vertex_descriptor new_vertex = graph.add_vertex(); + inserted_vertices[vertex_i] = new_vertex; + + double source_weight = get(vertex_label_cost_map, vd)[0]; + double sink_weight = get(vertex_label_cost_map, vd)[0]; + + graph.add_tweight(new_vertex, source_weight, sink_weight); + } +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + vertex_creation_time += timer.time(); + timer.reset(); +#endif + + // For E-Smooth + // add edge between every vertex, + for (input_edge_descriptor ed : CGAL::make_range(edges(input_graph))) + { + input_vertex_descriptor vd1 = source(ed, input_graph); + input_vertex_descriptor vd2 = target(ed, input_graph); + std::size_t idx1 = get(vertex_index_map, vd1); + std::size_t idx2 = get(vertex_index_map, vd2); + + double weight = get(edge_cost_map, ed); + + Vertex_descriptor v1 = inserted_vertices[idx1], + v2 = inserted_vertices[idx2]; + + std::size_t label_1 = get(vertex_label_map, vd1); + std::size_t label_2 = get(vertex_label_map, vd2); + graph.add_edge(v1, v2, weight, weight); + } + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + edge_creation_time += timer.time(); +#endif + + graph.init_vertices(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + timer.reset(); +#endif + + min_cut = graph.max_flow(); + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + cut_time += timer.time(); +#endif + + graph.get_labels(vertex_index_map, vertex_label_map, inserted_vertices, CGAL::make_range(vertices(input_graph))); + + //update labeling + for (auto vd : vertices(input_graph)) { + std::size_t idx = get(vertex_index_map, vd); + int label = graph.get_label(inserted_vertices[idx]); + put(vertex_label_map, vd, label); + } +/* + for (input_vertex_descriptor vd : CGAL::make_range(vertices(input_graph))) + { + std::size_t vertex_i = get(vertex_index_map, vd); + graph.update(vertex_label_map, inserted_vertices, vd, vertex_i, alpha); + }*/ + +#ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT + CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time << + std::endl; + CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl; + CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl; +#endif + + return min_cut; + } /// \cond SKIP_IN_MANUAL @@ -707,6 +832,23 @@ double alpha_expansion_graphcut (const std::vector +double min_cut(const std::vector >& edges, + const std::vector& edge_costs, + const std::vector >& cost_matrix, + std::vector& labels, + const AlphaExpansionImplementationTag&) +{ + internal::Alpha_expansion_old_API_wrapper_graph graph(edges, edge_costs, cost_matrix, labels); + + return min_cut(graph, + graph.edge_cost_map(), + graph.vertex_label_cost_map(), + graph.vertex_label_map(), + CGAL::parameters::vertex_index_map(graph.vertex_index_map()). + implementation_tag(AlphaExpansionImplementationTag())); +} /// \endcond }//namespace CGAL diff --git a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp index 6f8f9eea0cd..e287f1668a4 100644 --- a/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp +++ b/Kinetic_shape_reconstruction/examples/Kinetic_shape_reconstruction/kinetic_precomputed_shapes_example.cpp @@ -132,15 +132,11 @@ int main(const int argc, const char** argv) { std::back_inserter(output_faces), support_plane_idx, 6); assert(num_faces >= output_faces.size()); - int volume_level = -1; - const int num_volume_levels = ksr.number_of_volume_levels(); - assert(num_volume_levels > 0); - // Volumes. - const std::size_t num_volumes = ksr.number_of_volumes(volume_level); + const std::size_t num_volumes = ksr.number_of_volumes(); std::vector output_volumes; ksr.output_partition_volumes( - std::back_inserter(output_volumes), volume_level); + std::back_inserter(output_volumes)); assert(num_volumes == output_volumes.size()); // Support planes. @@ -163,7 +159,6 @@ int main(const int argc, const char** argv) { std::cout << "* number of faces: " << num_faces << std::endl; std::cout << "* number of volumes: " << num_volumes << std::endl; std::cout << "* number of support planes: " << num_support_planes << std::endl; - std::cout << "* number of volume levels: " << num_volume_levels << std::endl; std::cout << "* number of events: " << num_events << std::endl; // Export. diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h index 3dbf7b6603c..1ca0862d0d0 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR/conversions.h @@ -46,8 +46,8 @@ class Kinetic_traits_3 { using EK_to_IK = CGAL::Cartesian_converter; public: - Kinetic_traits_3(const bool use_hybrid_mode) : - m_use_hybrid_mode(use_hybrid_mode) { } + Kinetic_traits_3() : + m_use_hybrid_mode(false) { } inline const Point_2 intersection(const Line_2& t1, const Line_2& t2) const { diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/debug.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/debug.h index 3c9d98d2e59..60dd451e4f8 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/debug.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR/debug.h @@ -40,7 +40,7 @@ namespace KSR_3 { const std::tuple get_idx_color(std::size_t idx) { - CGAL::Random rand(idx); + CGAL::Random rand(static_cast(idx)); return std::make_tuple( static_cast(rand.get_int(32, 192)), static_cast(rand.get_int(32, 192)), @@ -118,24 +118,33 @@ void dump_2d_surface_mesh( std::vector map_vertices; map_vertices.clear(); + auto pvertices = data.pvertices(support_plane_idx); + //std::vector pts; + for (const auto pvertex : data.pvertices(support_plane_idx)) { if (map_vertices.size() <= pvertex.second) { map_vertices.resize(pvertex.second + 1); } + //pts.push_back(data.point_3(pvertex)); map_vertices[pvertex.second] = mesh.add_vertex(data.point_3(pvertex)); } + auto pfaces = data.pfaces(support_plane_idx); for (const auto pface : data.pfaces(support_plane_idx)) { vertices.clear(); + auto pvertices_of_pface = data.pvertices_of_pface(pface); for (const auto pvertex : data.pvertices_of_pface(pface)) { vertices.push_back(map_vertices[pvertex.second]); } CGAL_assertion(vertices.size() >= 3); const auto face = mesh.add_face(vertices); - CGAL_assertion(face != Mesh::null_face()); + if (face == Mesh::null_face()) { + std::cout << "could not dump mesh " << tag << std::endl; + return; + } std::tie(red[face], green[face], blue[face]) = - get_idx_color(support_plane_idx * (pface.second + 1)); + get_idx_color((support_plane_idx + 1) * (pface.second + 1)); } const std::string filename = (tag != std::string() ? tag + "-" : "") + "polygons.ply"; @@ -312,6 +321,7 @@ public: using FT = typename Traits::FT; using Point_2 = typename Traits::Point_2; using Point_3 = typename Traits::Point_3; + using Vector_3 = typename Traits::Vector_3; using Segment_2 = typename Traits::Segment_2; using Segment_3 = typename Traits::Segment_3; @@ -373,6 +383,54 @@ public: save(stream, file_name + ".xyz"); } + void export_regions_3( + const std::vector& points, + const std::vector& normals, + const std::vector& region, + const std::string file_name) const { + + if (points.size() != normals.size()) { + std::cout << "export_regions_3: number of points and normals are not equal" << std::endl; + return; + } + + if (points.size() != region.size()) { + std::cout << "export_regions_3: number of points and region indices are not equal" << std::endl; + return; + } + + std::stringstream stream; + initialize(stream); + + add_ply_header_regions(stream, points.size()); + + for (std::size_t i = 0; i < points.size(); ++i) { + stream << points[i] << " " << normals[i] << " " << region[i] << std::endl; + } + save(stream, file_name); + } + + void export_points_3( + const std::vector& points, + const std::vector& normals, + const std::string file_name) const { + + if (points.size() != normals.size()) { + std::cout << "export_regions_3: number of points and normals are not equal" << std::endl; + return; + } + + std::stringstream stream; + initialize(stream); + + add_ply_header_normals(stream, points.size()); + + for (std::size_t i = 0; i < points.size(); ++i) { + stream << points[i] << " " << normals[i] << " " << std::endl; + } + save(stream, file_name); + } + void export_segments_2( const std::vector& segments, const std::string file_name) const { @@ -516,7 +574,7 @@ public: } const Color get_idx_color(const std::size_t idx) const { - Random rand(idx); + Random rand(static_cast(idx)); const unsigned char r = rand.get_int(32, 192); const unsigned char g = rand.get_int(32, 192); const unsigned char b = rand.get_int(32, 192); @@ -566,16 +624,34 @@ private: const std::size_t size) const { stream << - "ply" + std::string(_NL_) + "" << - "format ascii 1.0" + std::string(_NL_) + "" << - "element vertex " << size << "" + std::string(_NL_) + "" << - "property double x" + std::string(_NL_) + "" << - "property double y" + std::string(_NL_) + "" << - "property double z" + std::string(_NL_) + "" << - "property double nx" + std::string(_NL_) + "" << - "property double ny" + std::string(_NL_) + "" << - "property double nz" + std::string(_NL_) + "" << - "end_header" + std::string(_NL_) + ""; + "ply" + std::string(_NL_) + "" << + "format ascii 1.0" + std::string(_NL_) + "" << + "element vertex " << size << "" + std::string(_NL_) + "" << + "property double x" + std::string(_NL_) + "" << + "property double y" + std::string(_NL_) + "" << + "property double z" + std::string(_NL_) + "" << + "property double nx" + std::string(_NL_) + "" << + "property double ny" + std::string(_NL_) + "" << + "property double nz" + std::string(_NL_) + "" << + "end_header" + std::string(_NL_) + ""; + } + + void add_ply_header_regions( + std::stringstream& stream, + const std::size_t size) const { + + stream << + "ply" << std::endl << + "format ascii 1.0" << std::endl << + "element vertex " << size << std::endl << + "property double x" << std::endl << + "property double y" << std::endl << + "property double z" << std::endl << + "property double nx" << std::endl << + "property double ny" << std::endl << + "property double nz" << std::endl << + "property int region" << std::endl << + "end_header" << std::endl; } void add_ply_header_mesh( @@ -605,7 +681,8 @@ void dump_volume( const DS& data, const std::vector& pfaces, const std::string file_name, - const bool use_colors = true) { + const bool use_colors = true, + const std::size_t volume_index = 0) { using Point_3 = typename DS::Kernel::Point_3; std::vector polygon; @@ -616,16 +693,18 @@ void dump_volume( polygons.reserve(pfaces.size()); Saver saver; + std::size_t i = 1; for (const auto& pface : pfaces) { const auto pvertices = data.pvertices_of_pface(pface); - const auto color = saver.get_idx_color(static_cast(pface.second)); + const auto color = saver.get_idx_color(volume_index); polygon.clear(); for (const auto pvertex : pvertices) { polygon.push_back(data.point_3(pvertex)); } if (use_colors) { colors.push_back(color); - } else { + } + else { colors.push_back(saver.get_idx_color(0)); } CGAL_assertion(polygon.size() >= 3); @@ -637,6 +716,43 @@ void dump_volume( saver.export_polygon_soup_3(polygons, colors, file_name); } +template +void dump_visi( + const DS& data, + const std::vector& pfaces, + const std::string &file_name, + const FT color) { + + using Point_3 = typename DS::Kernel::Point_3; + std::vector polygon; + std::vector< std::vector > polygons; + std::vector colors; + + colors.reserve(pfaces.size()); + polygons.reserve(pfaces.size()); + + const Color low(255, 255, 255); + const Color high(0, 0, 255); + + Saver saver; + for (const auto& pface : pfaces) { + const auto pvertices = data.pvertices_of_pface(pface); + polygon.clear(); + for (const auto pvertex : pvertices) { + polygon.push_back(data.point_3(pvertex)); + } + + colors.push_back(Color((1 - color) * low[0] + color * high[0], (1 - color) * low[1] + color * high[1], (1 - color) * low[2] + color * high[2], ((color > 0.5) ? 150 : 25))); + + CGAL_assertion(polygon.size() >= 3); + polygons.push_back(polygon); + } + CGAL_assertion(colors.size() == pfaces.size()); + CGAL_assertion(polygons.size() == pfaces.size()); + + saver.export_polygon_soup_3(polygons, colors, file_name); +} + template void dump_volumes(const DS& data, const std::string tag = std::string()) { @@ -687,6 +803,26 @@ void dump_polygon( saver.export_polygon_soup_3(polygons, name); } +template +void dump_polygons( + const DS& data, + const Polygon_2& input,// std::map< std::size_t, std::pair > polygons; + const std::string name) { + + using Kernel = typename DS::Kernel; + using Point_3 = typename Kernel::Point_3; + std::vector polygon; + std::vector< std::vector > polygons; + for (const auto& pair : input) { + for (const auto &point2d : pair.second.first) + polygon.push_back(data.to_3d(pair.first, point2d)); + polygons.push_back(polygon); + polygon.clear(); + } + Saver saver; + saver.export_polygon_soup_3(polygons, name); +} + template void dump_ifaces(const DS& data, const std::string tag = std::string()) { // write all polygons into a separate ply with support plane index and iface index @@ -727,6 +863,27 @@ void dump_pface( saver.export_polygon_soup_3(polygons, name); } +template +void dump_pfaces( + const DS& data, + const std::vector& pfaces, + const std::string name) { + + using Kernel = typename DS::Kernel; + using Point_3 = typename Kernel::Point_3; + std::vector< std::vector > polygons; + for (auto pface : pfaces) { + std::vector polygon; + for (const auto pvertex : data.pvertices_of_pface(pface)) { + polygon.push_back(data.point_3(pvertex)); + } + CGAL_assertion(polygon.size() >= 3); + polygons.push_back(polygon); + } + Saver saver; + saver.export_polygon_soup_3(polygons, name); +} + template void dump_pedge( const DS& data, @@ -745,13 +902,14 @@ void dump_info( const DS& data, const PFace& pface, const PEdge& pedge, - const std::vector& nfaces) { + const std::vector& nfaces, + const std::string &suffix) { std::cout << "DEBUG: number of found nfaces: " << nfaces.size() << std::endl; - dump_pface(data, pface, "face-curr"); - dump_pedge(data, pedge, "face-edge"); + dump_pface(data, pface, "face-curr-" + suffix); + dump_pedge(data, pedge, "face-edge-" + suffix); for (std::size_t i = 0; i < nfaces.size(); ++i) { - dump_pface(data, nfaces[i], "nface-" + std::to_string(i)); + dump_pface(data, nfaces[i], "nface-" + std::to_string(i) + "-" + suffix); } } @@ -816,6 +974,29 @@ void dump_cdt( out.close(); } +template +void dump(const InputRange input_range, PointMap point_map, NormalMap normal_map, Regions regions, const std::string &file_name) { + std::vector region_index(input_range.size(), -1); + + for (std::size_t r = 0; r < regions.size(); r++) { + for (std::size_t i = 0; i < regions[r].size(); i++) { + CGAL_assertion(regions[r][i] < input_range.size()); + region_index[regions[r][i]] = r; + } + } + + std::vector::value_type> pts(input_range.size()); + std::vector::value_type> normals(input_range.size()); + + for (std::size_t i = 0; i < input_range.size(); i++) { + pts[i] = get(point_map, i); + normals[i] = get(normal_map, i); + } + + Saver::value_type::R> saver; + saver.export_regions_3(pts, normals, region_index, file_name); +} + } // namespace KSR_3 } // namespace CGAL diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h index f69cdf5362a..6c9cad59b4f 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR/parameters.h @@ -25,14 +25,10 @@ struct Parameters_3 { unsigned int n = 0; // n subdivisions, not implemented yet FT enlarge_bbox_ratio = FT(11) / FT(10); // ratio to enlarge bbox - FT distance_tolerance = FT(5) / FT(10); // distance tolerance between planes + FT distance_tolerance = FT(0.005) / FT(10); // distance tolerance between planes bool reorient = false; // true - optimal bounding box, false - axis aligned - // true - apply exact interesections while all other computations are inexact, - // works only with EPICK as input kernel, switched off currently, see conversions.h - bool use_hybrid_mode = false; - // All files are saved in the current build directory. bool verbose = true; // print basic verbose information bool debug = false; // print all steps and substeps + export initial and final configurations diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR/utils.h b/Kinetic_shape_reconstruction/include/CGAL/KSR/utils.h index eb2584aad76..b76c9d36413 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR/utils.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR/utils.h @@ -92,17 +92,17 @@ point_3_from_point_2(const Point_2& point_2) { // Global tolerance. template static FT tolerance() { - return FT(1) / FT(100000); + return 0;// FT(1) / FT(100000); } template static FT point_tolerance() { - return tolerance(); + return 0;// tolerance(); } template static FT vector_tolerance() { - return FT(99999) / FT(100000); + return 0;// FT(99999) / FT(100000); } // Normalize vector. @@ -111,7 +111,7 @@ inline const Vector_d normalize(const Vector_d& v) { using Traits = typename Kernel_traits::Kernel; using FT = typename Traits::FT; const FT dot_product = CGAL::abs(v * v); - CGAL_assertion(dot_product != FT(0)); + //CGAL_assertion(dot_product != FT(0)); return v / static_cast(CGAL::sqrt(CGAL::to_double(dot_product))); } diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h index 04c30159d08..d9a7c42754a 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Data_structure.h @@ -181,6 +181,7 @@ public: struct Volume_cell { std::vector pfaces; + std::vector pface_oriented_outwards; std::vector neighbors; std::set pvertices; std::size_t index = std::size_t(-1); @@ -191,6 +192,9 @@ public: FT outside = FT(0); FT weight = FT(0); + FT inside_count = FT(0); + FT outside_count = FT(0); + void add_pface(const PFace& pface, const int neighbor) { pfaces.push_back(pface); neighbors.push_back(neighbor); @@ -253,6 +257,7 @@ public: m_reconstructed_model.clear(); } + void precompute_iedge_data() { for (std::size_t i = 0; i < number_of_support_planes(); ++i) { diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h index 9e2c5ee66d3..6d5bcd6f5e8 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Graphcut.h @@ -17,8 +17,9 @@ // CGAL includes. #include -#define CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE +//#define CGAL_DO_NOT_USE_BOYKOV_KOLMOGOROV_MAXFLOW_SOFTWARE #include +#include #include #include #include @@ -43,9 +44,11 @@ namespace KSR_3 { using FT = typename Kernel::FT; using Point_3 = typename Kernel::Point_3; using Triangle_2 = typename Kernel::Triangle_2; + using Triangle_3 = typename Kernel::Triangle_3; using Indices = std::vector; using Data_structure = KSR_3::Data_structure; + using Mesh = typename Data_structure::Mesh; using Volume_cell = typename Data_structure::Volume_cell; using PFace = typename Data_structure::PFace; @@ -60,7 +63,7 @@ namespace KSR_3 { PFace pface; FT weight = FT(0); std::pair neighbors; - bool is_boundary = false; + std::size_t support_plane; }; Graphcut( @@ -70,14 +73,14 @@ namespace KSR_3 { m_beta(graphcut_beta) { } - void compute(std::vector& volumes) { + void compute(std::vector& volumes, std::size_t inliers) { if (volumes.size() == 0) return; std::vector wrappers; create_pface_wrappers(wrappers); - compute_weights(wrappers); + compute_weights(wrappers, inliers); compute_weights(volumes); std::vector< std::pair > edges; @@ -97,6 +100,7 @@ namespace KSR_3 { private: const Data_structure& m_data; const FT m_beta; + FT m_total_area; void create_pface_wrappers( std::vector& wrappers) const { @@ -109,7 +113,7 @@ namespace KSR_3 { const auto pfaces = m_data.pfaces(i); for (const auto pface : pfaces) { wrapper.pface = pface; - wrapper.is_boundary = (i < 6) ? true : false; + wrapper.support_plane = i; CGAL_assertion(pface_neighbors.find(pface) != pface_neighbors.end()); const auto& pair = pface_neighbors.at(pface); wrapper.neighbors = pair; @@ -120,13 +124,23 @@ namespace KSR_3 { } void compute_weights( - std::vector& wrappers) const { + std::vector& wrappers, std::size_t inliers) { - FT sum = FT(0); + FT min = 10000000000; + FT max = 0; + + std::vector inside, boundary; + + m_total_area = 0; for (auto& wrapper : wrappers) { auto& weight = wrapper.weight; const auto& pface = wrapper.pface; + if (wrapper.neighbors.first < 0 || wrapper.neighbors.second < 0) { + boundary.push_back(pface); + } + else inside.push_back(pface); + Delaunay_2 tri; const auto pvertices = m_data.pvertices_of_pface(pface); for (const auto pvertex : pvertices) { @@ -144,12 +158,22 @@ namespace KSR_3 { fit->vertex(2)->point()); weight += triangle.area(); } - sum += weight; + + m_total_area += weight; + min = (std::min)(min, weight); + max = (std::max)(max, weight); } - CGAL_assertion(sum > FT(0)); + std::cout << "total area: " << m_total_area << " min: " << min << " max: " << max << " mean: " << (m_total_area / wrappers.size()) << std::endl; + + dump_pfaces(m_data, inside, "inside_pfaces.ply"); + dump_pfaces(m_data, boundary, "boundary_pfaces.ply"); + std::cout << inside.size() << " inside faces" << std::endl; + std::cout << boundary.size() << " boundary faces" << std::endl; + + CGAL_assertion(m_total_area > FT(0)); for (auto& wrapper : wrappers) { - wrapper.weight /= sum; + wrapper.weight = 2.0 * inliers * wrapper.weight / m_total_area; } } @@ -159,6 +183,8 @@ namespace KSR_3 { FT sum = FT(0); const Converter converter; + std::size_t index = 0; + for (auto& volume : volumes) { auto& weight = volume.weight; const auto& pvertices = volume.pvertices; @@ -176,12 +202,15 @@ namespace KSR_3 { weight += tet.volume(); } sum += weight; + index++; } CGAL_assertion(sum > FT(0)); for (auto& volume : volumes) { volume.weight /= sum; } + + std::cout << volumes.size() << " volumes" << std::endl; } void set_graph_edges( @@ -189,8 +218,11 @@ namespace KSR_3 { std::vector< std::pair >& edges, std::vector& edge_costs) const { + std::map, std::size_t> edges2index; + edges.clear(); edge_costs.clear(); + std::size_t internal = 0, external = 0; for (const auto& wrapper : wrappers) { const FT edge_weight = wrapper.weight; @@ -199,26 +231,61 @@ namespace KSR_3 { const int idx1 = neighbors.first; const int idx2 = neighbors.second; + std::size_t id1, id2; + + bool intern = false; + // Boundary edges. CGAL_assertion(idx1 >= 0 || idx2 >= 0); - if (idx1 < 0 && idx2 >= 0) { continue; } - if (idx2 < 0 && idx1 >= 0) { continue; } + if( (idx1 < 0 && idx2 >= 0) || (idx2 < 0 && idx1 >= 0)) { + if (idx1 < 0) { + id1 = wrapper.support_plane; + id2 = static_cast(idx2) + 6; + } + else { + id1 = static_cast(idx1) + 6; + id2 = wrapper.support_plane; + } + } + else { + intern = true; + // Internal edges. + CGAL_assertion(idx1 >= 0); + id1 = static_cast(idx1) + 6; + CGAL_assertion(idx2 >= 0); + id2 = static_cast(idx2) + 6; + } - // Internal edges. - CGAL_assertion(idx1 >= 0); - const std::size_t id1 = static_cast(idx1); - CGAL_assertion(idx2 >= 0); - const std::size_t id2 = static_cast(idx2); + if (id2 < id1) { + std::size_t tmp = id1; + id1 = id2; + id2 = tmp; + } + + std::pair idx(id1, id2); + + auto it = edges2index.find(idx); + if (it == edges2index.end()) { + edges2index[idx] = edges.size(); + edges.push_back(std::make_pair(id1, id2)); + edge_costs.push_back(compute_edge_cost(edge_weight)); + if (intern) + internal++; + else + external++; + } + else { + edge_costs[edges2index[idx]] += compute_edge_cost(edge_weight); + } - edges.push_back(std::make_pair(id1, id2)); - edge_costs.push_back(compute_edge_cost(edge_weight)); } + std::cout << internal << " internal " << external << " external edges" << std::endl; } double compute_edge_cost(const FT edge_weight) const { - CGAL_assertion(m_beta >= FT(0) && m_beta <= FT(1)); - CGAL_assertion(edge_weight >= FT(0) && edge_weight <= FT(1)); + //CGAL_assertion(m_beta >= FT(0) && m_beta <= FT(1)); + //CGAL_assertion(edge_weight >= FT(0) && edge_weight <= FT(1)); return CGAL::to_double(m_beta * edge_weight); } @@ -228,26 +295,97 @@ namespace KSR_3 { cost_matrix.clear(); cost_matrix.resize(2); - cost_matrix[0].resize(volumes.size()); - cost_matrix[1].resize(volumes.size()); + cost_matrix[0].resize(volumes.size() + 6); + cost_matrix[1].resize(volumes.size() + 6); - for (std::size_t i = 0; i < volumes.size(); ++i) { + FT min_in = 100000000, max_in = 0, mean_in = 0; + FT min_out = 100000000, max_out = 0, mean_out = 0; + int min_in_count = 100000000, max_in_count = -100000000; + int min_out_count = 100000000, max_out_count = -100000000; + + // Searching minimum values + for (std::size_t i = 0; i < volumes.size(); i++) { const auto& volume = volumes[i]; - - const FT in = volume.inside; - const FT out = volume.outside; - - CGAL_assertion(in >= FT(0) && in <= FT(1)); - CGAL_assertion(out >= FT(0) && out <= FT(1)); - CGAL_assertion((in + out) == FT(1)); - - const FT face_weight = volume.weight; - const double cost_in = get_face_cost(in , face_weight); - const double cost_out = get_face_cost(out, face_weight); - - cost_matrix[0][i] = cost_in; - cost_matrix[1][i] = cost_out; + min_in_count = (std::min)(min_in_count, static_cast(volume.inside_count)); + max_in_count = (std::max)(max_in_count, static_cast(volume.inside_count)); + min_out_count = (std::min)(min_out_count, static_cast(volume.outside_count)); + max_out_count = (std::max)(max_out_count, static_cast(volume.outside_count)); } + + std::cout << "min in count: " << min_in_count << " max in count: " << max_in_count << std::endl; + std::cout << "min out count: " << min_out_count << " max out count: " << max_out_count << std::endl; + + int minimum = (min_in_count < min_out_count) ? min_in_count : min_out_count; + + // Setting preferred outside label for bbox plane nodes + // Order: + // 0 zmin + // 1 ymin + // 2 xmax + // 3 ymax + // 4 xmin + // 5 zmax + + for (std::size_t i = 0; i < 6; i++) { + cost_matrix[0][i] = 100000000; + } + cost_matrix[0][0] = -100000000; + + for (std::size_t i = 0; i < 6; i++) { + cost_matrix[1][i] = -cost_matrix[0][i]; + } + + // Using new data term + if (true) { + for (std::size_t i = 0; i < volumes.size(); i++) { + cost_matrix[0][i + 6] = (volumes[i].outside_count - volumes[i].inside_count) * (1.0 - m_beta); + cost_matrix[1][i + 6] = (volumes[i].inside_count - volumes[i].outside_count) * (1.0 - m_beta); + //std::cout << i << ": " << cost_matrix[0][i + 6] << " " << cost_matrix[1][i + 6] << std::endl; + mean_in += cost_matrix[0][i + 6]; + mean_out += cost_matrix[1][i + 6]; + } + } + else { + for (std::size_t i = 0; i < volumes.size(); ++i) { + const auto& volume = volumes[i]; + + const FT in = volume.inside; + const FT out = volume.outside; + + CGAL_assertion(in >= FT(0) && in <= FT(1)); + CGAL_assertion(out >= FT(0) && out <= FT(1)); + CGAL_assertion((in + out) == FT(1)); + + const FT face_weight = volume.weight; + const double cost_in = get_face_cost(in, face_weight); + const double cost_out = get_face_cost(out, face_weight); + + cost_matrix[0][i + 6] = cost_in; + cost_matrix[1][i + 6] = cost_out; + + min_in = (std::min)(min_in, cost_in); + max_in = (std::max)(max_in, cost_in); + mean_in += cost_in; + + min_out = (std::min)(min_out, cost_out); + max_out = (std::max)(max_out, cost_out); + mean_out += cost_out; + + min_in_count = (std::min)(min_in_count, volume.inside_count); + max_in_count = (std::max)(max_in_count, volume.inside_count); + min_out_count = (std::min)(min_out_count, volume.outside_count); + max_out_count = (std::max)(max_out_count, volume.outside_count); + + //std::cout << volume.index << " in: " << cost_in << " before: " << in << " " << volume.inside_count << std::endl; + //std::cout << " out: " << cost_out << " before " << out << " " << volume.outside_count << std::endl; + } + } + + mean_in /= FT(volumes.size()); + mean_out /= FT(volumes.size()); + + //std::cout << "min in: " << min_in << " max in: " << max_in << " mean: " << mean_in << std::endl; + //std::cout << "min out: " << min_out << " max out: " << max_out << " mean: " << mean_out << std::endl; } double get_face_cost( @@ -266,14 +404,18 @@ namespace KSR_3 { std::vector& labels) const { labels.clear(); - labels.resize(volumes.size()); + labels.resize(volumes.size() + 6); + + for (std::size_t i = 0; i < 6; i++) { + labels[i] = 1; + } for (std::size_t i = 0; i < volumes.size(); ++i) { const auto& volume = volumes[i]; if (volume.visibility == Visibility_label::INSIDE) { - labels[i] = 0; + labels[i + 6] = 0; } else { - labels[i] = 1; + labels[i + 6] = 1; } } } @@ -283,18 +425,68 @@ namespace KSR_3 { const std::vector& edge_costs, const std::vector< std::vector >& cost_matrix, std::vector& labels) const { + std::vector tmp = labels; + + std::cout << "beta" << m_beta << std::endl; + + double min = 10000000000; + double max = -min; + for (std::size_t i = 0; i < edge_costs.size(); i++) { + min = (std::min)(edge_costs[i], min); + max = (std::max)(edge_costs[i], max); + } + + std::cout << "edge costs" << std::endl; + std::cout << "min: " << min << std::endl; + std::cout << "max: " << max << std::endl; + + min = 1000000000; + max = -min; + + for (std::size_t i = 6; i < cost_matrix[0].size(); i++) { + min = (std::min)(cost_matrix[0][i], min); + max = (std::max)(cost_matrix[0][i], max); + } + + std::cout << "label costs" << std::endl; + std::cout << "min: " << min << std::endl; + std::cout << "max: " << max << std::endl; + +/* + CGAL::min_cut( + edges, edge_costs, cost_matrix, labels, CGAL::parameters::implementation_tag(CGAL::Alpha_expansion_MaxFlow_tag())); + + bool difference = false; + for (std::size_t i = 0; i < labels.size(); i++) { + if (tmp[i] != labels[i]) { + difference = true; + break; + } + } + std::cout << "Labels changed: " << difference << std::endl; + */ CGAL::alpha_expansion_graphcut( - edges, edge_costs, cost_matrix, labels); + edges, edge_costs, cost_matrix, labels, CGAL::parameters::implementation_tag(CGAL::Alpha_expansion_MaxFlow_tag())); + + bool difference = false; + for (std::size_t i = 0; i < labels.size(); i++) { + if (tmp[i] != labels[i]) { + difference = true; + break; + } + } + + std::cout << "Labels changed: " << difference << std::endl; } void apply_new_labels( const std::vector& labels, std::vector& volumes) const { - CGAL_assertion(volumes.size() == labels.size()); - for (std::size_t i = 0; i < labels.size(); ++i) { - const std::size_t label = labels[i]; + CGAL_assertion((volumes.size() + 6) == labels.size()); + for (std::size_t i = 0; i < volumes.size(); ++i) { + const std::size_t label = labels[i + 6]; auto& volume = volumes[i]; if (label == 0) { diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h index 1eaf98ea668..6220a8d45a2 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Reconstruction.h @@ -202,6 +202,7 @@ public: template bool detect_planar_shapes( + const std::string& file_name, const NamedParameters& np) { if (m_verbose) { @@ -216,12 +217,13 @@ public: create_approximate_walls(np); create_approximate_roofs(np); } else { - create_all_planar_shapes(np); + create_all_planar_shapes(np, file_name); } CGAL_assertion(m_planes.size() == m_polygons.size()); CGAL_assertion(m_polygons.size() == m_region_map.size()); - if (m_debug) dump_polygons("detected-planar-shapes"); + //if (m_debug) + dump_polygons("detected-planar-shapes"); if (m_polygons.size() == 0) { if (m_verbose) std::cout << "* no planar shapes found" << std::endl; @@ -230,6 +232,37 @@ public: return true; } + template + bool planar_shapes_from_map( + RegionMap ®ion_map, + const NamedParameters& np + ) { + + if (m_verbose) { + std::cout << std::endl << "--- PLANAR SHAPES from map: " << std::endl; + } + m_planes.clear(); + m_polygons.clear(); + m_region_map.clear(); + + std::vector > regions; + + for (std::size_t i = 0; i < m_input_range.size(); i++) { + if (region_map[i] < 0) + continue; + if (regions.size() <= region_map[i]) + regions.resize(region_map[i] + 1); + regions[region_map[i]].push_back(i); + } + + for (const auto& region : regions) { + const auto plane = fit_plane(region); + const std::size_t shape_idx = add_planar_shape(region, plane); + CGAL_assertion(shape_idx != std::size_t(-1)); + m_region_map[shape_idx] = region; + } + } + template bool regularize_planar_shapes( const NamedParameters& np) { @@ -310,7 +343,7 @@ public: std::cout << std::endl << "--- COMPUTING THE MODEL: " << std::endl; } - if (m_data.number_of_volumes(-1) == 0) { + if (m_data.number_of_volumes() == 0) { if (m_verbose) std::cout << "* no volumes found, skipping" << std::endl; return false; } @@ -323,7 +356,7 @@ public: CGAL_assertion(m_data.volumes().size() > 0); visibility.compute(m_data.volumes()); - // if (m_debug) dump_volumes("visibility/visibility"); + dump_visibility("visibility/visibility", pface_points); if (m_verbose) { std::cout << "done" << std::endl; @@ -334,8 +367,8 @@ public: parameters::get_parameter(np, internal_np::graphcut_beta), FT(1) / FT(2)); Graphcut graphcut(m_data, beta); - graphcut.compute(m_data.volumes()); - // if (m_debug) dump_volumes("graphcut/graphcut"); + graphcut.compute(m_data.volumes(), visibility.inliers()); + dump_volumes("graphcut/graphcut"); if (m_verbose) { std::cout << "done" << std::endl; @@ -596,14 +629,14 @@ private: } template - void create_all_planar_shapes(const NamedParameters& np) { + void create_all_planar_shapes(const NamedParameters& np, const std::string& file_name) { if (m_free_form_points.size() < 3) { if (m_verbose) std::cout << "* no points found, skipping" << std::endl; return; } if (m_verbose) std::cout << "* getting planar shapes using region growing" << std::endl; - const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_free_form_points); + const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_free_form_points, file_name); if (m_verbose) std::cout << "* found " << num_shapes << " approximate walls" << std::endl; } @@ -615,7 +648,7 @@ private: return; } if (m_verbose) std::cout << "* getting walls using region growing" << std::endl; - const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_boundary_points); + const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_boundary_points, ""); if (m_verbose) std::cout << "* found " << num_shapes << " approximate walls" << std::endl; } @@ -627,17 +660,22 @@ private: return; } if (m_verbose) std::cout << "* getting roofs using region growing" << std::endl; - const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_interior_points); + const std::size_t num_shapes = compute_planar_shapes_with_rg(np, m_interior_points, ""); if (m_verbose) std::cout << "* found " << num_shapes << " approximate roofs" << std::endl; } template std::size_t compute_planar_shapes_with_rg( const NamedParameters& np, - const std::vector& input_range) { + const std::vector& input_range, + const std::string& file_name) { std::vector< std::vector > regions; apply_region_growing_3(np, input_range, regions); + + if (file_name.size() > 0) + dump(m_input_range, m_point_map_3, m_normal_map_3, regions, file_name); + for (const auto& region : regions) { const auto plane = fit_plane(region); const std::size_t shape_idx = add_planar_shape(region, plane); @@ -717,7 +755,7 @@ private: std::size_t num_shapes = 0; if (wall_points.size() >= 3) { - num_shapes += compute_planar_shapes_with_rg(np, wall_points); + num_shapes += compute_planar_shapes_with_rg(np, wall_points, ""); // dump_polygons("walls-1"); } @@ -1238,6 +1276,19 @@ private: } } + void dump_visibility(const std::string file_name, const std::map &pface_points) { + for (const auto& volume : m_data.volumes()) { + std::size_t sample_count = 0; + for (auto pface : volume.pfaces) { + const auto indices = pface_points.at(pface); + sample_count += indices.size(); + } + dump_visi(m_data, volume.pfaces, + file_name + "-" + std::to_string(volume.index) + "-" + std::to_string(volume.inside_count) + "-" + + std::to_string(volume.outside_count) + "-" + std::to_string(sample_count), (volume.inside_count)/(volume.inside_count + volume.outside_count)); + } + } + void dump_model(const std::string file_name) { std::vector polygon; diff --git a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h index 4d8f2765cd1..9a65833ddc6 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h +++ b/Kinetic_shape_reconstruction/include/CGAL/KSR_3/Visibility.h @@ -70,18 +70,44 @@ namespace KSR_3 { m_pface_points(pface_points), m_point_map_3(point_map_3), m_normal_map_3(normal_map_3), - m_num_samples(100), + m_num_samples(0), m_random(0) { CGAL_assertion(m_pface_points.size() > 0); + m_inliers = 0; + for (const auto pair : m_pface_points) { + m_inliers += pair.second.size(); + } + std::cout << "inliers: " << m_inliers << std::endl; } void compute(std::vector& volumes) const { CGAL_assertion(volumes.size() > 0); if (volumes.size() == 0) return; + std::size_t i = 0; + FT xmin, ymin, zmin; + xmin = ymin = zmin = 10000000; + FT xmax, ymax, zmax; + xmax = ymax = zmax = -xmin; for (auto& volume : volumes) { estimate_volume_label(volume); + const Point_3& c = volume.centroid; + xmin = (std::min)(xmin, c.x()); + xmax = (std::max)(xmax, c.x()); + ymin = (std::min)(ymin, c.y()); + ymax = (std::max)(ymax, c.y()); + zmin = (std::min)(zmin, c.z()); + zmax = (std::max)(zmax, c.z()); + i++; } + std::cout << "Sampled " << m_num_samples << " for data term" << std::endl; + std::cout << "x: " << xmin << " " << xmax << std::endl; + std::cout << "y: " << ymin << " " << ymax << std::endl; + std::cout << "z: " << zmin << " " << zmax << std::endl; + } + + std::size_t inliers() const { + return m_inliers; } private: @@ -89,7 +115,8 @@ namespace KSR_3 { const std::map& m_pface_points; const Point_map_3& m_point_map_3; const Vector_map_3& m_normal_map_3; - const std::size_t m_num_samples; + mutable std::size_t m_num_samples; + mutable std::size_t m_inliers; Random m_random; void estimate_volume_label(Volume_cell& volume) const { @@ -97,10 +124,8 @@ namespace KSR_3 { const auto stats = estimate_in_out_values(volume); CGAL_assertion(stats.first >= FT(0) && stats.first <= FT(1)); CGAL_assertion(stats.second >= FT(0) && stats.second <= FT(1)); - CGAL_assertion( - CGAL::abs(stats.first + stats.second - FT(1)) < KSR::tolerance()); - if (stats.first >= FT(1) / FT(2)) { + if (volume.inside_count > volume.outside_count) { volume.visibility = Visibility_label::INSIDE; } else { volume.visibility = Visibility_label::OUTSIDE; @@ -113,12 +138,14 @@ namespace KSR_3 { } const std::pair estimate_in_out_values( - const Volume_cell& volume) const { + Volume_cell& volume) const { std::size_t in = 0, out = 0; std::vector samples; create_samples(volume, samples); compute_stats( volume, samples, in, out); + volume.inside_count = in; + volume.outside_count = out; if (in == 0 && out == 0) { in = 1; out = 1; } @@ -194,6 +221,9 @@ namespace KSR_3 { const Point_3& query, std::size_t& in, std::size_t& out) const { + std::vector inside, outside; + std::vector insideN, outsideN; + bool found = false; const auto& pfaces = volume.pfaces; for (const auto& pface : pfaces) { @@ -202,6 +232,8 @@ namespace KSR_3 { if (indices.size() == 0) continue; found = true; + m_num_samples += indices.size(); + for (const std::size_t index : indices) { const auto& point = get(m_point_map_3 , index); const auto& normal = get(m_normal_map_3, index); @@ -209,15 +241,27 @@ namespace KSR_3 { const Vector_3 vec(point, query); const FT dot_product = vec * normal; if (dot_product < FT(0)) { + inside.push_back(point); + insideN.push_back(normal); in += 1; - } else { + } + else { + outside.push_back(point); + outsideN.push_back(normal); out += 1; } } } - if (!found) { - out += 1; return false; + if (volume.index != -1) { + std::ofstream vout("visibility/" + std::to_string(volume.index) + "-query.xyz"); + vout.precision(20); + vout << query << std::endl; + vout.close(); + + Saver saver; + saver.export_points_3(inside, insideN, "visibility/" + std::to_string(volume.index) + "-inside.ply"); + saver.export_points_3(outside, outsideN, "visibility/" + std::to_string(volume.index) + "-outside.ply"); } return true; } diff --git a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h index 2dc01b89d6c..76cf6c6fd17 100644 --- a/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h +++ b/Kinetic_shape_reconstruction/include/CGAL/Kinetic_shape_reconstruction_3.h @@ -242,168 +242,16 @@ public: template< typename InputRange, - typename PolygonMap, + typename PointMap, + typename VectorMap, + typename SemanticMap, typename NamedParameters> - bool partition_by_faces( - const InputRange& input_range, - const PolygonMap polygon_map, - const NamedParameters& np) { - - Timer timer; - m_parameters.k = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::k_intersections), 1); - m_parameters.n = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::n_subdivisions), 0); - m_parameters.enlarge_bbox_ratio = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::enlarge_bbox_ratio), FT(11) / FT(10)); - m_parameters.distance_tolerance = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::distance_tolerance), FT(5) / FT(10)); - m_parameters.reorient = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::reorient), false); - m_parameters.use_hybrid_mode = parameters::choose_parameter( - parameters::get_parameter(np, internal_np::use_hybrid_mode), false); - - std::cout.precision(20); - if (input_range.size() == 0) { - CGAL_warning_msg(input_range.size() > 0, - "WARNING: YOUR INPUT IS EMPTY! RETURN WITH NO CHANGE!"); - return false; - } - - if (m_parameters.n != 0) { - CGAL_assertion_msg(false, "TODO: IMPLEMENT KINETIC SUBDIVISION!"); - if (m_parameters.n > 3) { - CGAL_warning_msg(m_parameters.n <= 3, - "WARNING: DOES IT MAKE SENSE TO HAVE MORE THAN 64 INPUT BLOCKS? SETTING N TO 3!"); - m_parameters.n = 3; - } - } - - if (m_parameters.enlarge_bbox_ratio < FT(1)) { - CGAL_warning_msg(m_parameters.enlarge_bbox_ratio >= FT(1), - "WARNING: YOU SET ENLARGE_BBOX_RATIO < 1.0! THE VALID RANGE IS [1.0, +INF). SETTING TO 1.0!"); - m_parameters.enlarge_bbox_ratio = FT(1); - } - - if (m_parameters.verbose) { - const unsigned int num_blocks = std::pow(m_parameters.n + 1, 3); - const std::string is_reorient = (m_parameters.reorient ? "true" : "false"); - - std::cout << std::endl << "--- PARTITION OPTIONS: " << std::endl; - std::cout << "* number of intersections k: " << m_parameters.k << std::endl; - std::cout << "* number of subdivisions per bbox side: " << m_parameters.n << std::endl; - std::cout << "* number of subdivision blocks: " << num_blocks << std::endl; - std::cout << "* enlarge bbox ratio: " << m_parameters.enlarge_bbox_ratio << std::endl; - std::cout << "* reorient: " << is_reorient << std::endl; - std::cout << "* hybrid mode: " << m_parameters.use_hybrid_mode << std::endl; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- INITIALIZING PARTITION:" << std::endl; - } - - // Initialization. - timer.reset(); - timer.start(); - m_data.clear(); - - Initializer initializer(m_data, m_parameters); - initializer.initialize(input_range, polygon_map); - - timer.stop(); - const double time_to_initialize = timer.time(); - - if (m_parameters.verbose) { - std::cout << std::endl << "* initialization (sec.): " << time_to_initialize << std::endl; - std::cout << "INITIALIZATION SUCCESS!" << std::endl << std::endl; - } - - if (m_parameters.k == 0) { // for k = 0, we skip propagation - CGAL_warning_msg(m_parameters.k > 0, - "WARNING: YOU SET K TO 0! THAT MEANS NO PROPAGATION! THE VALID VALUES ARE {1,2,...}. INTERSECT AND RETURN!"); - return false; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- RUNNING THE QUEUE:" << std::endl; - std::cout << "* propagation started" << std::endl; - } - - // FacePropagation. - timer.reset(); - timer.start(); - std::size_t num_queue_calls = 0; - FacePropagation propagation(m_data, m_parameters); - std::tie(num_queue_calls, m_num_events) = propagation.propagate(); - timer.stop(); - const double time_to_propagate = timer.time(); - - if (m_parameters.verbose) { - std::cout << "* propagation finished" << std::endl; - std::cout << "* number of queue calls: " << num_queue_calls << std::endl; - std::cout << "* number of events handled: " << m_num_events << std::endl; - } - - if (m_parameters.verbose) { - std::cout << std::endl << "--- FINALIZING PARTITION:" << std::endl; - } - - // Finalization. - timer.reset(); - timer.start(); - if (m_parameters.debug) dump(m_data, "jiter-final-a-result"); - - Finalizer finalizer(m_data, m_parameters); - finalizer.clean(); - - if (m_parameters.verbose) std::cout << "* checking final mesh integrity ..."; - CGAL_assertion(m_data.check_integrity(true, true, true)); - if (m_parameters.verbose) std::cout << " done" << std::endl; - - if (m_parameters.debug) dump(m_data, "jiter-final-b-result"); - // std::cout << std::endl << "CLEANING SUCCESS!" << std::endl << std::endl; - // exit(EXIT_SUCCESS); - - if (m_parameters.verbose) std::cout << "* getting volumes ..." << std::endl; - finalizer.create_polyhedra(); - timer.stop(); - const double time_to_finalize = timer.time(); - if (m_parameters.verbose) { - std::cout << "* found all together " << m_data.number_of_volumes(-1) << " volumes" << std::endl; - } - // std::cout << std::endl << "CREATING VOLUMES SUCCESS!" << std::endl << std::endl; - // exit(EXIT_SUCCESS); - - for (std::size_t i = 0; i < m_data.number_of_support_planes(); i++) { - dump_2d_surface_mesh(m_data, i, "final-surface-mesh-" + std::to_string(i)); - } - - // Timing. - if (m_parameters.verbose) { - std::cout << std::endl << "--- TIMING (sec.):" << std::endl; - } - const double total_time = - time_to_initialize + time_to_propagate + time_to_finalize; - if (m_parameters.verbose) { - std::cout << "* initialization: " << time_to_initialize << std::endl; - std::cout << "* propagation: " << time_to_propagate << std::endl; - std::cout << "* finalization: " << time_to_finalize << std::endl; - std::cout << "* total time: " << total_time << std::endl; - } - return true; - } - - template< - typename InputRange, - typename PointMap, - typename VectorMap, - typename SemanticMap, - typename NamedParameters> bool reconstruct( const InputRange& input_range, const PointMap point_map, const VectorMap normal_map, const SemanticMap semantic_map, + const std::string file_name, const NamedParameters& np) { using Reconstruction = KSR_3::Reconstruction< @@ -411,7 +259,8 @@ public: Reconstruction reconstruction( input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); - bool success = reconstruction.detect_planar_shapes(np); + + bool success = reconstruction.detect_planar_shapes(file_name, np); if (!success) { CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); return false; @@ -441,6 +290,57 @@ public: return success; } + template< + typename InputRange, + typename PointMap, + typename VectorMap, + typename SemanticMap, + typename RegionMap, + typename NamedParameters> + bool reconstruct( + const InputRange& input_range, + const PointMap point_map, + const VectorMap normal_map, + const SemanticMap semantic_map, + const RegionMap region_map, + const NamedParameters& np) { + + using Reconstruction = KSR_3::Reconstruction< + InputRange, PointMap, VectorMap, SemanticMap, Kernel>; + + Reconstruction reconstruction( + input_range, point_map, normal_map, semantic_map, m_data, m_parameters.verbose, m_parameters.debug); + + bool success = reconstruction.planar_shapes_from_map(region_map, np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, DETECTING PLANAR SHAPES FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + //success = reconstruction.regularize_planar_shapes(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, REGULARIZATION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = partition( + reconstruction.planar_shapes(), reconstruction.polygon_map(), np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, PARTITION FAILED!"); + return false; + } + // exit(EXIT_SUCCESS); + + success = reconstruction.compute_model(np); + if (!success) { + CGAL_assertion_msg(false, "ERROR: RECONSTRUCTION, COMPUTING MODEL FAILED!"); + return false; + } + return success; + } + /******************************* ** STATISTICS ** ********************************/