diff --git a/.gitattributes b/.gitattributes index 9ae0185ba62..40509b80fa1 100644 --- a/.gitattributes +++ b/.gitattributes @@ -4016,11 +4016,13 @@ Surface_mesh_segmentation/example/Surface_mesh_segmentation/segmentation_example Surface_mesh_segmentation/example/Surface_mesh_segmentation/segmentation_example_sdf.cpp -text Surface_mesh_segmentation/example/Surface_mesh_segmentation/simple_example.cpp -text Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h -text +Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h -text Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_traversal_traits.h -text Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h -text Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h -text Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h -text Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h -text +Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h -text Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/copyright -text Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/license.txt -text Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/maintainer -text diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h index 60ca6de6486..4734a14aff5 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h @@ -24,6 +24,7 @@ #include #include +#include #include #include @@ -32,12 +33,15 @@ #include #include #include - +#include //AF: macros must be prefixed with "CGAL_" //IOY: done #define CGAL_NORMALIZATION_ALPHA 5.0 #define CGAL_CONVEX_FACTOR 0.08 - +#define CGAL_DEFAULT_NUMBER_OF_CLUSTERS 5 +#define CGAL_DEFAULT_SMOOTHING_LAMBDA 23.0 +#define CGAL_DEFAULT_CONE_ANGLE (2.0 / 3.0) * CGAL_PI +#define CGAL_DEFAULT_NUMBER_OF_RAYS 25 //IOY: these are going to be removed at the end (no CGAL_ pref) #define ACTIVATE_SEGMENTATION_DEBUG #ifdef ACTIVATE_SEGMENTATION_DEBUG @@ -55,7 +59,11 @@ namespace CGAL * It is a connector class which uses soft clustering and graph cut in order to segment meshes. * All preprocessing and postprocessing issues are handled here. */ -template > +template < +class Polyhedron, + class FacetIndexMap = + boost::associative_property_map > + > class Surface_mesh_segmentation { //type definitions @@ -65,81 +73,110 @@ public: typedef typename Polyhedron::Facet Vertex; typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::Point_3 Point; - typedef typename Polyhedron::Facet_iterator Facet_iterator; - typedef typename Polyhedron::Facet_handle Facet_handle; - typedef typename Polyhedron::Halfedge_handle Halfedge_handle; - typedef typename Polyhedron::Edge_iterator Edge_iterator; - typedef typename Polyhedron::Vertex_handle Vertex_handle; - typedef typename SDFCalculation::Parameters SDF_Parameters; + typedef typename Polyhedron::Edge_const_iterator Edge_const_iterator; + typedef typename Polyhedron::Halfedge_const_iterator Halfedge_const_iterator; + typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator; + typedef typename Polyhedron::Vertex_const_iterator Vertex_const_iterator; + + typedef internal::SDF_calculation SDF_calculation; protected: typedef typename Kernel::Plane_3 Plane; - typedef std::map Face_value_map; - typedef std::map Face_center_map; - typedef std::map Face_segment_map; +// member variables +public: // going to be protected ! + const Polyhedron& mesh; - template - struct Compare_second_element { - bool operator()(const ValueTypeName& v1, const ValueTypeName& v2) const { - return v1.second < v2.second; - } - }; + FacetIndexMap facet_index_map; -//member variables -public: - Polyhedron* mesh; - - Face_value_map sdf_values; - Face_center_map centers; - Face_segment_map segments; + std::vector sdf_values; + std::vector centers; + std::vector segments; int number_of_centers; double smoothing_lambda; + std::map facet_index_map_internal; - internal::Expectation_maximization fitter;/**< going to be removed */ +// member functions -//member functions public: - Surface_mesh_segmentation(Polyhedron* mesh): mesh(mesh) { -#ifdef SEGMENTATION_PROFILE - profile("profile.txt"); -#endif + Surface_mesh_segmentation(const Polyhedron& mesh) + : mesh(mesh), facet_index_map(facet_index_map_internal) { + int facet_index = 0; + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); + ++facet_it, ++facet_index) { + boost::put(facet_index_map, facet_it, facet_index); + } } - void calculate_sdf_values(SDF_Parameters parameters = SDF_Parameters()) { - SEG_DEBUG(CGAL::Timer t) + Surface_mesh_segmentation(Polyhedron* mesh, FacetIndexMap facet_index_map) + : mesh(mesh), facet_index_map(facet_index_map) { + } + + void calculate_sdf_values(double cone_angle = CGAL_DEFAULT_CONE_ANGLE, + int number_of_ray = CGAL_DEFAULT_NUMBER_OF_RAYS) { + typedef std::map internal_map; + internal_map facet_value_map_internal; + boost::associative_property_map sdf_pmap( + facet_value_map_internal); + + SDF_calculation(cone_angle, number_of_ray).calculate_sdf_values(mesh, sdf_pmap); + + sdf_values = std::vector(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + get(sdf_values, facet_it) = boost::get(sdf_pmap, facet_it); + } + + check_zero_sdf_values(); + smooth_sdf_values_with_bilateral(); + linear_normalize_sdf_values(); + } + + template + void calculate_sdf_values(SDFPropertyMap sdf_pmap, double cone_angle, + int number_of_ray) { + SEG_DEBUG(Timer t) SEG_DEBUG(t.start()) - sdf_values.clear(); - SDFCalculation(parameters).calculate_sdf_values(*mesh, sdf_values); + SDF_calculation(cone_angle, number_of_rays).calculate_sdf_values(mesh, + sdf_pmap); + + sdf_values = std::vector(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + get(sdf_values, facet_it) = boost::get(sdf_pmap, facet_it); + } SEG_DEBUG(std::cout << t.time() << std::endl) check_zero_sdf_values(); smooth_sdf_values_with_bilateral(); - normalize_sdf_values(); + linear_normalize_sdf_values(); + + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + boost::put(sdf_pmap, facet_it, get(sdf_values, facet_it)); + } SEG_DEBUG(std::cout << t.time() << std::endl) } - void partition(int number_of_centers = 5, double smoothing_lambda = 23.0) { + template + void partition(FacetSegmentMap segment_pmap, + int number_of_centers = CGAL_DEFAULT_NUMBER_OF_CLUSTERS, + double smoothing_lambda = CGAL_DEFAULT_SMOOTHING_LAMBDA) { this->number_of_centers = number_of_centers; this->smoothing_lambda = smoothing_lambda; centers.clear(); - - std::vector sdf_vector; - sdf_vector.reserve(sdf_values.size()); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it) { - sdf_vector.push_back(sdf_values[facet_it]); - } + // log normalize sdf values + normalize_sdf_values(); // soft clustering using GMM-fitting initialized with k-means - fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, - internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); + internal::Expectation_maximization fitter(number_of_centers, sdf_values, + internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); std::vector labels; fitter.fill_with_center_ids(labels); @@ -156,21 +193,62 @@ public: // apply graph cut internal::Alpha_expansion_graph_cut gc(edges, edge_weights, probability_matrix, labels); - - std::vector::iterator center_it = labels.begin(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } + centers = labels; // assign a segment id for each facet assign_segments(); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + boost::put(segment_pmap, facet_it, get(segments, facet_it)); + } } -public: - double calculate_dihedral_angle_of_edge(const Halfedge_handle& edge) const { - Facet_handle f1 = edge->facet(); - Facet_handle f2 = edge->opposite()->facet(); + template + void partition(SDFPropertyMap sdf_pmap, FacetSegmentMap segment_pmap, + int number_of_centers = CGAL_DEFAULT_NUMBER_OF_CLUSTERS, + double smoothing_lambda = CGAL_DEFAULT_SMOOTHING_LAMBDA) { + sdf_values = std::vector(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + get(sdf_values, facet_it) = boost::get(sdf_pmap, facet_it); + } + + this->number_of_centers = number_of_centers; + this->smoothing_lambda = smoothing_lambda; + centers.clear(); + // log normalize sdf values + normalize_sdf_values(); + // soft clustering using GMM-fitting initialized with k-means + internal::Expectation_maximization fitter(number_of_centers, sdf_values, + internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); + + std::vector labels; + fitter.fill_with_center_ids(labels); + + std::vector > probability_matrix; + fitter.fill_with_probabilities(probability_matrix); + log_normalize_probability_matrix(probability_matrix); + + // calculating edge weights + std::vector > edges; + std::vector edge_weights; + calculate_and_log_normalize_dihedral_angles(edges, edge_weights); + + // apply graph cut + internal::Alpha_expansion_graph_cut gc(edges, edge_weights, probability_matrix, + labels); + centers = labels; + // assign a segment id for each facet + assign_segments(); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + boost::put(segment_pmap, facet_it, get(segments, facet_it)); + } + } + +protected: + double calculate_dihedral_angle_of_edge(Edge_const_iterator& edge) const { + Facet_const_iterator f1 = edge->facet(); + Facet_const_iterator f2 = edge->opposite()->facet(); const Point& f2_v1 = f2->halfedge()->vertex()->point(); const Point& f2_v2 = f2->halfedge()->next()->vertex()->point(); @@ -204,14 +282,14 @@ public: return angle; } - double calculate_dihedral_angle_of_edge_2(const Halfedge_handle& edge) const { + double calculate_dihedral_angle_of_edge_2(Edge_const_iterator& edge) const { const Point& a = edge->vertex()->point(); const Point& b = edge->prev()->vertex()->point(); const Point& c = edge->next()->vertex()->point(); const Point& d = edge->opposite()->next()->vertex()->point(); // As far as I check: if, say, dihedral angle is 5, this returns 175, // if dihedral angle is -5, this returns -175. - double n_angle = CGAL::Mesh_3::dihedral_angle(a, b, c, d) / 180.0; + double n_angle = Mesh_3::dihedral_angle(a, b, c, d) / 180.0; bool concave = n_angle > 0; double angle = 1 + ((concave ? -1 : +1) * n_angle); @@ -220,8 +298,8 @@ public: } return angle; - //Facet_handle f1 = edge->facet(); - //Facet_handle f2 = edge->opposite()->facet(); + //Facet_const_iterator f1 = edge->facet(); + //Facet_const_iterator f2 = edge->opposite()->facet(); // //const Point& f2_v1 = f2->halfedge()->vertex()->point(); //const Point& f2_v2 = f2->halfedge()->next()->vertex()->point(); @@ -240,8 +318,8 @@ public: //const Point& f1_v1 = f1->halfedge()->vertex()->point(); //const Point& f1_v2 = f1->halfedge()->next()->vertex()->point(); //const Point& f1_v3 = f1->halfedge()->prev()->vertex()->point(); - //Vector f1_normal = CGAL::unit_normal(f1_v1, f1_v2, f1_v3); - //Vector f2_normal = CGAL::unit_normal(f2_v1, f2_v2, f2_v3); + //Vector f1_normal = unit_normal(f1_v1, f1_v2, f1_v3); + //Vector f2_normal = unit_normal(f2_v1, f2_v2, f2_v3); // //double dot = f1_normal * f2_normal; //if(dot > 1.0) { dot = 1.0; } @@ -257,37 +335,46 @@ public: } void normalize_sdf_values() { - typedef typename Face_value_map::iterator fv_iterator; - Compare_second_element comparator; + typedef std::vector::iterator fv_iterator; std::pair min_max_pair = - CGAL::min_max_element(sdf_values.begin(), sdf_values.end(), comparator, - comparator); + min_max_element(sdf_values.begin(), sdf_values.end()); - double max_value = min_max_pair.second->second, - min_value = min_max_pair.first->second; + double max_value = *min_max_pair.second, min_value = *min_max_pair.first; double max_min_dif = max_value - min_value; - for(fv_iterator pair_it = sdf_values.begin(); pair_it != sdf_values.end(); - ++pair_it) { - double linear_normalized = (pair_it->second - min_value) / max_min_dif; + for(fv_iterator it = sdf_values.begin(); it != sdf_values.end(); ++it) { + double linear_normalized = (*it - min_value) / max_min_dif; double log_normalized = log(linear_normalized * CGAL_NORMALIZATION_ALPHA + 1) / log(CGAL_NORMALIZATION_ALPHA + 1); - pair_it->second = log_normalized; + *it = log_normalized; + } + } + + void linear_normalize_sdf_values() { + typedef std::vector::iterator fv_iterator; + std::pair min_max_pair = + min_max_element(sdf_values.begin(), sdf_values.end()); + + double max_value = *min_max_pair.second, min_value = *min_max_pair.first; + double max_min_dif = max_value - min_value; + for(fv_iterator it = sdf_values.begin(); it != sdf_values.end(); ++it) { + *it = (*it - min_value) / max_min_dif; } } void smooth_sdf_values() { - Face_value_map smoothed_sdf_values; - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - Facet_handle f = pair_it->first; - typename Facet::Halfedge_around_facet_circulator facet_circulator = - f->facet_begin(); + std::vector smoothed_sdf_values(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + typename Facet::Halfedge_around_facet_const_circulator facet_circulator = + facet_it->facet_begin(); double total_neighbor_sdf = 0.0; do { - total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; - } while( ++facet_circulator != f->facet_begin()); + total_neighbor_sdf += get(sdf_values, facet_circulator->opposite()->facet()); + } while( ++facet_circulator != facet_it->facet_begin()); + total_neighbor_sdf /= 3.0; - smoothed_sdf_values[f] = (pair_it->second + total_neighbor_sdf) / 2.0; + get(smoothed_sdf_values, facet_it) = (get(sdf_values, + facet_it) + total_neighbor_sdf) / 2.0; } sdf_values = smoothed_sdf_values; } @@ -298,23 +385,22 @@ public: const int iteration = 1; for(int i = 0; i < iteration; ++i) { - Face_value_map smoothed_sdf_values; - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - Facet_handle facet = pair_it->first; - std::map neighbors; - get_neighbors_by_vertex(facet, neighbors, window_size); + std::vector smoothed_sdf_values(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + std::map neighbors; + get_neighbors_by_vertex(facet_it, neighbors, window_size); double total_sdf_value = 0.0; double total_weight = 0.0; - for(typename std::map::iterator it = neighbors.begin(); - it != neighbors.end(); ++it) { + for(typename std::map::iterator it = + neighbors.begin(); it != neighbors.end(); ++it) { double weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0), 2))); // window_size => 2*sigma - total_sdf_value += sdf_values[it->first] * weight; + total_sdf_value += get(sdf_values, it->first) * weight; total_weight += weight; } - smoothed_sdf_values[facet] = total_sdf_value / total_weight; + get(smoothed_sdf_values, facet_it) = total_sdf_value / total_weight; } sdf_values = smoothed_sdf_values; } @@ -326,17 +412,16 @@ public: const int iteration = 1; for(int i = 0; i < iteration; ++i) { - Face_value_map smoothed_sdf_values; - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { + std::vector smoothed_sdf_values(mesh.size_of_facets());; + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { //Find neighbors and put their sdf values into a list - Facet_handle facet = pair_it->first; - std::map neighbors; - get_neighbors_by_vertex(facet, neighbors, window_size); + std::map neighbors; + get_neighbors_by_vertex(facet_it, neighbors, window_size); std::vector sdf_of_neighbors; - for(typename std::map::iterator it = neighbors.begin(); - it != neighbors.end(); ++it) { - sdf_of_neighbors.push_back(sdf_values[it->first]); + for(typename std::map::iterator it = + neighbors.begin(); it != neighbors.end(); ++it) { + sdf_of_neighbors.push_back(get(sdf_values, it->first)); } // Find median. double median_sdf = 0.0; @@ -351,7 +436,7 @@ public: } else { median_sdf = sdf_of_neighbors[half_neighbor_count]; } - smoothed_sdf_values[facet] = median_sdf; + get(smoothed_sdf_values, facet_it) = median_sdf; } sdf_values = smoothed_sdf_values; } @@ -366,51 +451,51 @@ public: const int iteration = 1; for(int i = 0; i < iteration; ++i) { - Face_value_map smoothed_sdf_values; - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - Facet_handle facet = pair_it->first; - std::map neighbors; - get_neighbors_by_vertex(facet, neighbors, window_size); + std::vector smoothed_sdf_values(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + //Facet_handle facet = facet_it; + std::map neighbors; + get_neighbors_by_vertex(facet_it, neighbors, window_size); double total_sdf_value = 0.0, total_weight = 0.0; - double current_sdf_value = sdf_values[facet]; + double current_sdf_value = get(sdf_values, facet_it); // calculate deviation for range weighting. double deviation = 0.0; - for(typename std::map::iterator it = neighbors.begin(); - it != neighbors.end(); ++it) { - deviation += std::pow(sdf_values[it->first] - current_sdf_value, 2); + for(typename std::map::iterator it = + neighbors.begin(); it != neighbors.end(); ++it) { + deviation += std::pow(get(sdf_values, it->first) - current_sdf_value, 2); } deviation = std::sqrt(deviation / neighbors.size()); if(deviation == 0.0) { deviation = std::numeric_limits::epsilon(); //this might happen } - for(typename std::map::iterator it = neighbors.begin(); + for(std::map::iterator it = neighbors.begin(); it != neighbors.end(); ++it) { double spatial_weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0), 2))); // window_size => 2*sigma - double domain_weight = exp(-0.5 * (std::pow((sdf_values[it->first] - - current_sdf_value) / (std::sqrt(2.0)*deviation), 2))); + double domain_weight = exp(-0.5 * (std::pow( (get(sdf_values, + it->first) - current_sdf_value) / (std::sqrt(2.0)*deviation), 2))); double weight = spatial_weight * domain_weight; - total_sdf_value += sdf_values[it->first] * weight; + total_sdf_value += get(sdf_values, it->first) * weight; total_weight += weight; } - smoothed_sdf_values[facet] = total_sdf_value / total_weight; + get(smoothed_sdf_values, facet_it) = total_sdf_value / total_weight; } sdf_values = smoothed_sdf_values; } } - void get_neighbors_by_edge(const Facet_handle& facet, - std::map& neighbors, int max_level) { - typedef std::pair facet_level_pair; + void get_neighbors_by_edge(Facet_const_iterator& facet, + std::map& neighbors, int max_level) { + typedef std::pair facet_level_pair; std::queue facet_queue; facet_queue.push(facet_level_pair(facet, 0)); while(!facet_queue.empty()) { const facet_level_pair& pair = facet_queue.front(); bool inserted = neighbors.insert(pair).second; if(inserted && pair.second < max_level) { - typename Facet::Halfedge_around_facet_circulator facet_circulator = + typename Facet::Halfedge_around_facet_const_circulator facet_circulator = pair.first->facet_begin(); do { facet_queue.push(facet_level_pair(facet_circulator->opposite()->facet(), @@ -421,20 +506,20 @@ public: } } - void get_neighbors_by_vertex(const Facet_handle& facet, - std::map& neighbors, int max_level) { - typedef std::pair facet_level_pair; + void get_neighbors_by_vertex(Facet_const_iterator& facet, + std::map& neighbors, int max_level) { + typedef std::pair facet_level_pair; std::queue facet_queue; facet_queue.push(facet_level_pair(facet, 0)); while(!facet_queue.empty()) { const facet_level_pair& pair = facet_queue.front(); bool inserted = neighbors.insert(pair).second; if(inserted && pair.second < max_level) { - Facet_handle facet = pair.first; - Halfedge_handle edge = facet->halfedge(); + Facet_const_iterator facet = pair.first; + Halfedge_const_iterator edge = facet->halfedge(); do { // loop on three vertices of the facet - Vertex_handle vertex = edge->vertex(); - typename Facet::Halfedge_around_vertex_circulator vertex_circulator = + Vertex_const_iterator vertex = edge->vertex(); + typename Facet::Halfedge_around_vertex_const_circulator vertex_circulator = vertex->vertex_begin(); do { // for each vertex loop on neighbor vertices facet_queue.push(facet_level_pair(vertex_circulator->facet(), pair.second + 1)); @@ -444,34 +529,35 @@ public: facet_queue.pop(); } } + void check_zero_sdf_values() { // If there is any facet which has no sdf value, assign average sdf value of its neighbors - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - if(pair_it->second == 0.0) { - typename Facet::Halfedge_around_facet_circulator facet_circulator = - pair_it->first->facet_begin(); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + if(get(sdf_values, facet_it) == 0.0) { + typename Facet::Halfedge_around_facet_const_circulator facet_circulator = + facet_it->facet_begin(); double total_neighbor_sdf = 0.0; do { - total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; - } while( ++facet_circulator != pair_it->first->facet_begin()); - pair_it->second = total_neighbor_sdf / 3.0; + total_neighbor_sdf += get(sdf_values, facet_circulator->opposite()->facet()); + } while( ++facet_circulator != facet_it->facet_begin()); + get(sdf_values, facet_it) = total_neighbor_sdf / 3.0; } } // Find minimum sdf value other than 0 double min_sdf = (std::numeric_limits::max)(); - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - if(pair_it->second < min_sdf && pair_it->second != 0.0) { - min_sdf = pair_it->second; + for(std::vector::iterator it = sdf_values.begin(); + it != sdf_values.end(); ++it) { + if(*it < min_sdf && *it != 0.0) { + min_sdf = *it; } } // If still there is any facet which has no sdf value, assign minimum sdf value. // This is meaningful since (being an outlier) 0 sdf values might effect normalization & log extremely. - for(typename Face_value_map::iterator pair_it = sdf_values.begin(); - pair_it != sdf_values.end(); ++pair_it) { - if(pair_it->second == 0.0) { - pair_it->second = min_sdf; + for(std::vector::iterator it = sdf_values.begin(); + it != sdf_values.end(); ++it) { + if(*it == 0.0) { + *it = min_sdf; } } } @@ -495,19 +581,11 @@ public: std::vector >& edges, std::vector& edge_weights) const { const double epsilon = 1e-5; - //assign an id for every facet (facet-id) - std::map facet_indices; - int index = 0; - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++index) { - facet_indices.insert(std::pair(facet_it, index)); - } - //edges and their weights. pair stores facet-id pairs (see above) (may be using CGAL::Triple can be more suitable) - for(Edge_iterator edge_it = mesh->edges_begin(); edge_it != mesh->edges_end(); - ++edge_it) { - int index_f1 = facet_indices[edge_it->facet()]; - int index_f2 = facet_indices[edge_it->opposite()->facet()]; + //edges and their weights. pair stores facet-id pairs (see above) (may be using boost::tuple can be more suitable) + for(Edge_const_iterator edge_it = mesh.edges_begin(); + edge_it != mesh.edges_end(); ++edge_it) { + int index_f1 = boost::get(facet_index_map, edge_it->facet()); + int index_f2 = boost::get(facet_index_map, edge_it->opposite()->facet()); edges.push_back(std::pair(index_f1, index_f2)); double angle = calculate_dihedral_angle_of_edge(edge_it); @@ -522,240 +600,47 @@ public: } void assign_segments() { - segments.clear(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { - segments.insert(std::pair(facet_it, -1)); - } + segments = std::vector(mesh.size_of_facets(), -1); int segment_id = 0; - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { - if(segments[facet_it] == -1) { + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + if(get(segments, facet_it) == -1) { depth_first_traversal(facet_it, segment_id); segment_id++; } } } - void depth_first_traversal(const Facet_handle& facet, int segment_id) { - segments[facet] = segment_id; - typename Facet::Halfedge_around_facet_circulator facet_circulator = + void depth_first_traversal(Facet_const_iterator& facet, int segment_id) { + get(segments, facet) = segment_id; + typename Facet::Halfedge_around_facet_const_circulator facet_circulator = facet->facet_begin(); double total_neighbor_sdf = 0.0; do { - Facet_handle neighbor = facet_circulator->opposite()->facet(); - if(segments[neighbor] == -1 && centers[facet] == centers[neighbor]) { + Facet_const_iterator neighbor = facet_circulator->opposite()->facet(); + if(get(segments, neighbor) == -1 + && get(centers, facet) == get(centers, neighbor)) { depth_first_traversal(neighbor, segment_id); } } while( ++facet_circulator != facet->facet_begin()); } - /** - * Going to be removed - */ - void apply_GMM_fitting() { - centers.clear(); - std::vector sdf_vector; - sdf_vector.reserve(sdf_values.size()); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it) { - sdf_vector.push_back(sdf_values[facet_it]); - } - SEG_DEBUG(CGAL::Timer t) - SEG_DEBUG(t.start()) - //internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10); - - fitter = internal::Expectation_maximization(number_of_centers, sdf_vector); - //fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 40); - SEG_DEBUG(std::cout << "GMM fitting time: " << t.time() << std::endl) - std::vector center_memberships; - fitter.fill_with_center_ids(center_memberships); - std::vector::iterator center_it = center_memberships.begin(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } - } - /** - * Going to be removed - */ - void apply_K_means_clustering() { - centers.clear(); - std::vector sdf_vector; - sdf_vector.reserve(sdf_values.size()); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it) { - sdf_vector.push_back(sdf_values[facet_it]); - } - internal::K_means_clustering clusterer(number_of_centers, sdf_vector); - std::vector center_memberships; - clusterer.fill_with_center_ids(center_memberships); - std::vector::iterator center_it = center_memberships.begin(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } - //center_memberships_temp = center_memberships; //remove - } - /** - * Going to be removed - */ - void apply_GMM_fitting_with_K_means_init() { - centers.clear(); - std::vector sdf_vector; - sdf_vector.reserve(sdf_values.size()); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it) { - sdf_vector.push_back(sdf_values[facet_it]); - } - std::vector center_memberships; - fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, - internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); - center_memberships.clear(); - fitter.fill_with_center_ids(center_memberships); - std::vector::iterator center_it = center_memberships.begin(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } - } - /** - * Going to be removed - */ - void apply_GMM_fitting_and_K_means() { - centers.clear(); - std::vector sdf_vector; - sdf_vector.reserve(sdf_values.size()); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it) { - sdf_vector.push_back(sdf_values[facet_it]); - } - internal::Expectation_maximization gmm_random_init(number_of_centers, - sdf_vector); - - internal::K_means_clustering k_means(number_of_centers, sdf_vector); - std::vector center_memberships; - k_means.fill_with_center_ids(center_memberships); - internal::Expectation_maximization gmm_k_means_init(number_of_centers, - sdf_vector, - internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); - - if(gmm_k_means_init.final_likelihood > gmm_random_init.final_likelihood) { - fitter = gmm_k_means_init; - } else { - fitter = gmm_random_init; - } - center_memberships.clear(); - fitter.fill_with_center_ids(center_memberships); - std::vector::iterator center_it = center_memberships.begin(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } - } - /** - * Going to be removed - */ - void apply_graph_cut() { - - std::vector > edges; - std::vector edge_weights; - calculate_and_log_normalize_dihedral_angles(edges, edge_weights); - std::vector > probability_matrix; - fitter.fill_with_probabilities(probability_matrix); - - std::vector labels; - fitter.fill_with_center_ids(labels); - - log_normalize_probability_matrix(probability_matrix); - ////////////////////////////////////////////////////////////// - // FOR READING FROM MATLAB, GOING TO BE REMOVED -// std::ifstream cc("D:/GSoC/Matlab/ccount.txt"); -// cc >> number_of_centers; -// std::vector > probability_matrix(number_of_centers, std::vector(sdf_values.size(), 0.0)); -// read_center_ids("D:/GSoC/Matlab/result.txt"); -// read_probabilities("D:/GSoC/Matlab/probs.txt", probability_matrix); -// int f = 0; -// for (Facet_iterator facet_it = mesh->facets_begin(); facet_it != mesh->facets_end(); facet_it++, f++) - //{ - // for(int i = 0; i < number_of_centers; ++i) - // { - // double probability = probability_matrix[i][f]; - // probability += 1e-8; -// probability = (std::min)(probability, 1.0); -// probability = -log(probability); -// probability_matrix[i][f] = (std::max)(probability, std::numeric_limits::epsilon()); - // } - //} -// -// std::vector labels; -// for(Facet_iterator facet_it = mesh->facets_begin(); facet_it != mesh->facets_end(); -// ++facet_it) -// { -// labels.push_back(centers[facet_it]); -// } - ////////////////////////////////////////////////////////////// - - //apply graph cut - internal::Alpha_expansion_graph_cut gc(edges, edge_weights, probability_matrix, - labels); - - std::vector::iterator center_it = labels.begin(); - centers.clear(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); - ++facet_it, ++center_it) { - centers.insert(std::pair(facet_it, (*center_it))); - } - } - - /** - * Going to be removed - */ - void select_cluster_number() { - int min_cluster_count = 3; - int max_cluster_count = 5; - int range = max_cluster_count - min_cluster_count + 1; - std::vector distortions(range+1); - for(int i = min_cluster_count -1; i <= max_cluster_count; ++i) { - number_of_centers = i; - apply_GMM_fitting_with_K_means_init(); - double distortion = fitter.calculate_distortion(); - distortions[i-(min_cluster_count -1)] = std::pow(distortion, -0.5); - } - double max_jump = 0.0; - for(int i = 1; i < range+1; ++i) { - double jump = distortions[i] - distortions[i-1]; - if(jump > max_jump) { - max_jump = jump; - number_of_centers = i + min_cluster_count - 1; - } - } - for(int i = 0; i < range + 1; ++i) { - std::cout << "d: " << distortions[i] << std::endl; - } - std::cout << "noc: " << number_of_centers << std::endl; - //apply_GMM_fitting_and_K_means_init(); + template + T& get(std::vector& data, const Facet_const_iterator& facet) { + return data[ boost::get(facet_index_map, facet) ]; } +public: /** * Going to be removed */ void write_sdf_values(const char* file_name) { std::ofstream output(file_name); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { - output << sdf_values[facet_it] << std::endl; + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + output << get(sdf_values, facet_it) << std::endl; } output.close(); } @@ -764,31 +649,29 @@ public: */ void read_sdf_values(const char* file_name) { std::ifstream input(file_name); - sdf_values.clear(); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { + sdf_values = std::vector(mesh.size_of_facets()); + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { double sdf_value; input >> sdf_value; - sdf_values.insert(std::pair(facet_it, sdf_value)); + get(sdf_values, facet_it) = sdf_value; } } /** * Going to be removed */ void read_center_ids(const char* file_name) { - std::ifstream input(file_name); - centers.clear(); - int max_center = 0; - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { - int center_id; - input >> center_id; - centers.insert(std::pair(facet_it, center_id)); - if(center_id > max_center) { - max_center = center_id; - } - } - number_of_centers = max_center + 1; + //std::ifstream input(file_name); + //centers.clear(); + //int max_center = 0; + //for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); ++facet_it) + //{ + // int center_id; + // input >> center_id; + // centers.insert(std::pair(facet_it, center_id)); + // if(center_id > max_center) { max_center = center_id; } + //} + //number_of_centers = max_center + 1; } /** @@ -812,9 +695,9 @@ public: void write_segment_ids(const char* file_name) { assign_segments(); std::ofstream output(file_name); - for(Facet_iterator facet_it = mesh->facets_begin(); - facet_it != mesh->facets_end(); ++facet_it) { - output << segments[facet_it] << std::endl; + for(Facet_const_iterator facet_it = mesh.facets_begin(); + facet_it != mesh.facets_end(); ++facet_it) { + output << get(segments, facet_it) << std::endl; } output.close(); } @@ -837,7 +720,7 @@ public: use_minimum_segment = true; multiplier_for_segment = 1.0 + i; - CGAL::Timer t; + Timer t; t.start(); calculate_sdf_values(); @@ -877,7 +760,7 @@ public: use_minimum_segment = false; multiplier_for_segment = 1.0 + i * 0.2; - CGAL::Timer t; + Timer t; t.start(); calculate_sdf_values(); @@ -920,6 +803,10 @@ public: #undef CGAL_NORMALIZATION_ALPHA #undef CGAL_CONVEX_FACTOR +#undef CGAL_DEFAULT_NUMBER_OF_CLUSTERS +#undef CGAL_DEFAULT_SMOOTHING_LAMBDA +#undef CGAL_DEFAULT_CONE_ANGLE +#undef CGAL_DEFAULT_NUMBER_OF_RAYS #ifdef SEG_DEBUG #undef SEG_DEBUG diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h new file mode 100644 index 00000000000..e688e5047af --- /dev/null +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h @@ -0,0 +1,142 @@ +// Copyright (c) 2009 INRIA Sophia-Antipolis (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/features/gsoc2012-segmentation-iyaz/AABB_tree/include/CGAL/AABB_polyhedron_triangle_primitive.h $ +// $Id: AABB_polyhedron_triangle_primitive.h 67117 2012-01-13 18:14:48Z lrineau $ +// +// +// Author(s) : Stéphane Tayeb, Pierre Alliez +// +//****************************************************************************** +// File Description : +// +//****************************************************************************** + +#ifndef CGAL_AABB_POLYHEDRON_TRIANGLE_PRIMITIVE_H_ +#define CGAL_AABB_POLYHEDRON_TRIANGLE_PRIMITIVE_H_ + +namespace CGAL +{ + +/** +* @class AABB_polyhedron_triangle_primitive +* +* +*/ +template +class AABB_polyhedron_triangle_primitive +{ +public: + /// AABBPrimitive types + typedef typename GeomTraits::Point_3 Point; + typedef typename GeomTraits::Triangle_3 Datum; + typedef typename Polyhedron::Facet_const_iterator Id; + + /// Self + typedef AABB_polyhedron_triangle_primitive Self; + + /// Constructors + AABB_polyhedron_triangle_primitive() {} + AABB_polyhedron_triangle_primitive(const AABB_polyhedron_triangle_primitive& + primitive) { + m_facet_handle = primitive.id(); + } + AABB_polyhedron_triangle_primitive(const Id& handle) + : m_facet_handle(handle) { }; + AABB_polyhedron_triangle_primitive(const Id* ptr) + : m_facet_handle(*ptr) { }; + template + AABB_polyhedron_triangle_primitive( Iterator it, + typename boost::enable_if< + boost::is_same + >::type* =0 + ) : m_facet_handle(*it) { } + + + // Default destructor, copy constructor and assignment operator are ok + + /// Returns by constructing on the fly the geometric datum wrapped by the primitive + Datum datum() const { + const Point& a = m_facet_handle->halfedge()->vertex()->point(); + const Point& b = m_facet_handle->halfedge()->next()->vertex()->point(); + const Point& c = m_facet_handle->halfedge()->next()->next()->vertex()->point(); + return Datum(a,b,c); + } + + /// Returns a point on the primitive + Point reference_point() const { + return m_facet_handle->halfedge()->vertex()->point(); + } + + /// Returns the identifier + const Id& id() const { + return m_facet_handle; + } + Id& id() { + return m_facet_handle; + } + +private: + /// The id, here a polyhedron facet handle + Id m_facet_handle; +}; // end class AABB_polyhedron_triangle_primitive + + +/** +* @class AABB_const_polyhedron_triangle_primitive +* +* +*/ +template +class AABB_const_polyhedron_triangle_primitive +{ +public: + /// AABBPrimitive types + typedef typename GeomTraits::Point_3 Point; + typedef typename GeomTraits::Triangle_3 Datum; + typedef typename Polyhedron::Facet_const_handle Id; + + /// Constructors + AABB_const_polyhedron_triangle_primitive(const Id& handle) + : m_facet_handle(handle) { }; + + // Default destructor, copy constructor and assignment operator are ok + + /// Returns by constructing on the fly the geometric datum wrapped by the primitive + Datum datum() const { + const Point& a = m_facet_handle->halfedge()->vertex()->point(); + const Point& b = m_facet_handle->halfedge()->next()->vertex()->point(); + const Point& c = m_facet_handle->halfedge()->next()->next()->vertex()->point(); + return Datum(a,b,c); + } + + /// Returns a point on the primitive + Point reference_point() const { + return m_facet_handle->halfedge()->vertex()->point(); + } + + /// Returns the identifier + Id id() const { + return m_facet_handle; + } + +private: + /// The id, here a polyhedron facet handle + Id m_facet_handle; +}; // end class AABB_polyhedron_triangle_primitive + +} // end namespace CGAL + + +#endif // CGAL_AABB_POLYHEDRON_TRIANGLE_PRIMITIVE_H_ diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h index 40527ff5161..80cb735ff60 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h @@ -1,5 +1,5 @@ -#ifndef CGAL_ALPHA_EXPANSION_GRAPH_CUT_H -#define CGAL_ALPHA_EXPANSION_GRAPH_CUT_H +#ifndef CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H +#define CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H #include #include @@ -82,12 +82,12 @@ public: int number_of_clusters = probability_matrix.size(); double min_cut = (std::numeric_limits::max)(); bool success; - CGAL::Timer gt; + Timer gt; gt.start(); do { success = false; for(int alpha = 0; alpha < number_of_clusters; ++alpha) { - CGAL::Timer t; + Timer t; t.start(); Graph graph; #if 0 @@ -174,13 +174,13 @@ public: int number_of_clusters = probability_matrix.size(); double min_cut = (std::numeric_limits::max)(); bool success; - CGAL::Timer gt; + Timer gt; gt.start(); double total_time = 0.0; do { success = false; for(int alpha = 0; alpha < number_of_clusters; ++alpha) { - CGAL::Timer t; + Timer t; t.start(); Graph graph(edges.size() + labels.size() + 2); // allocate using maximum possible size. @@ -272,7 +272,7 @@ public: int number_of_clusters = probability_matrix.size(); double min_cut = (std::numeric_limits::max)(); bool success; - CGAL::Timer gt; + Timer gt; gt.start(); int cluster_source = labels.size(); int cluster_sink = labels.size() + 1; @@ -285,7 +285,7 @@ public: do { success = false; for(int alpha = 0; alpha < number_of_clusters; ++alpha) { - CGAL::Timer t; + Timer t; t.start(); Graph graph(edges.begin(), edges.end(), labels.size()+2, edges.size()); #if 0 @@ -305,4 +305,4 @@ public: }; }//namespace internal }//namespace CGAL -#endif //CGAL_ALPHA_EXPANSION_GRAPH_CUT_H \ No newline at end of file +#endif //CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H \ No newline at end of file diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h index f28ce5d2043..f1aa599ada0 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h @@ -1,5 +1,5 @@ -#ifndef CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H -#define CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H +#ifndef CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H +#define CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H /* NEED TO BE DONE */ /* About implementation: * + There are a lot of parameters (with default values) and initialization types, @@ -243,8 +243,11 @@ protected: Gaussian_center new_center(initial_mean, initial_deviation, initial_mixing_coefficient); // if same point is choosen as a center twice, algorithm will not work - if ( is_already_center(new_center) ) --i; - else centers.push_back(new_center); + if(is_already_center(new_center)) { + --i; + } else { + centers.push_back(new_center); + } } calculate_initial_deviations(); } @@ -289,8 +292,11 @@ protected: Gaussian_center new_center(initial_mean, initial_deviation, initial_mixing_coefficient); // if same point is choosen as a center twice, algorithm will not work - if ( is_already_center(new_center) ) --i; - else centers.push_back(new_center); + if(is_already_center(new_center)) { + --i; + } else { + centers.push_back(new_center); + } } calculate_initial_deviations(); } @@ -501,4 +507,4 @@ protected: #ifdef SEG_DEBUG #undef SEG_DEBUG #endif -#endif //CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H \ No newline at end of file +#endif //CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H \ No newline at end of file diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h index e61e1d49254..62a473d2376 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h @@ -1,5 +1,5 @@ -#ifndef CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H -#define CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H +#ifndef CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H +#define CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H #include #include @@ -79,7 +79,7 @@ protected: /** * Finds closest center and adds itself to the closest center's mean. * @param centers available centers - * @return true if #center_id is changed + * @return true if #center_id is changed (i.e. previous center is changed) */ bool calculate_new_center(std::vector& centers) { int new_center_id = 0; @@ -156,8 +156,11 @@ protected: for(int i = 0; i < number_of_centers; ++i) { double initial_mean = points[rand() % points.size()].data; K_means_center new_center(initial_mean); - if ( is_already_center(new_center) ) --i; - else centers.push_back(new_center); + if(is_already_center(new_center)) { + --i; + } else { + centers.push_back(new_center); + } } } @@ -195,8 +198,11 @@ protected: - distance_square_cumulative.begin(); double initial_mean = points[selection_index].data; K_means_center new_center(initial_mean); - if ( is_already_center(new_center) ) --i; - else centers.push_back(new_center); + if(is_already_center(new_center)) { + --i; + } else { + centers.push_back(new_center); + } } } @@ -293,4 +299,4 @@ protected: #undef CGAL_DEFAULT_MAXIMUM_ITERATION #undef CGAL_DEFAULT_NUMBER_OF_RUN -#endif //CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H +#endif //CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H diff --git a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h index 212c44b7535..766a7663634 100644 --- a/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h +++ b/Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h @@ -1,43 +1,35 @@ -#ifndef CGAL_SDF_CALCULATION_H -#define CGAL_SDF_CALCULATION_H +#ifndef CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H +#define CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H #include #include #include #include +#include #include #include #include #include +#include #define CGAL_ANGLE_ST_DEV_DIVIDER 2.0 #define CGAL_ACCEPTANCE_RATE_THRESHOLD 0.5 #define CGAL_ST_DEV_MULTIPLIER 0.75 -#define CGAL_DEFAULT_CONE_ANGLE (2.0 / 3.0) * CGAL_PI -#define CGAL_DEFAULT_NUMBER_OF_RAYS 25 - namespace CGAL { namespace internal { -template +template < +class Polyhedron +> class SDF_calculation { public: - struct Parameters { - public: - double cone_angle; - int number_of_rays; - Parameters(double cone_angle = CGAL_DEFAULT_CONE_ANGLE, - int number_of_rays = CGAL_DEFAULT_NUMBER_OF_RAYS) - : cone_angle(cone_angle), number_of_rays(number_of_rays) { - } - }; //type definitions protected: typedef typename Polyhedron::Traits Kernel; @@ -45,17 +37,16 @@ protected: typedef typename Polyhedron::Facet Vertex; typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::Point_3 Point; - typedef typename Polyhedron::Facet_iterator Facet_iterator; - typedef typename Polyhedron::Facet_handle Facet_handle; - typedef typename Polyhedron::Vertex_handle Vertex_handle; + + typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator; typedef typename Kernel::Ray_3 Ray; typedef typename Kernel::Plane_3 Plane; typedef typename Kernel::Segment_3 Segment; - typedef typename CGAL::AABB_polyhedron_triangle_primitive + typedef typename AABB_const_polyhedron_triangle_primitive Primitive; - typedef typename CGAL::AABB_tree > + typedef typename CGAL::AABB_tree > Tree; typedef typename Tree::Object_and_primitive_id Object_and_primitive_id; @@ -68,7 +59,8 @@ protected: //Member variables protected: - Parameters parameters; + double cone_angle; + int number_of_rays; Disk_samples_list disk_samples_sparse; Disk_samples_list disk_samples_dense; @@ -77,31 +69,31 @@ protected: double multiplier_for_segment; public: - SDF_calculation(Parameters parameters = Parameters()) - : parameters(parameters), use_minimum_segment(false), - multiplier_for_segment(1) { + SDF_calculation(double cone_angle, int number_of_rays) + : cone_angle(cone_angle), number_of_rays(number_of_rays), + use_minimum_segment(false), multiplier_for_segment(1) { } - void calculate_sdf_values(Polyhedron& mesh, - std::map& sdf_values) { - int sparse_ray_count = parameters.number_of_rays; + template < class FacetValueMap > + void calculate_sdf_values(const Polyhedron& mesh, FacetValueMap sdf_values) { + int sparse_ray_count = number_of_rays; int dense_ray_count = sparse_ray_count * 2; disk_sampling_vogel_method(disk_samples_sparse, sparse_ray_count); disk_sampling_vogel_method(disk_samples_dense, dense_ray_count); Tree tree(mesh.facets_begin(), mesh.facets_end()); - for(Facet_iterator facet_it = mesh.facets_begin(); + for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); ++facet_it) { CGAL_precondition(facet_it->is_triangle()); //Mesh should contain triangles. double sdf = calculate_sdf_value_of_facet(facet_it, tree, disk_samples_sparse); - sdf_values.insert(std::pair(facet_it, sdf)); + boost::put(sdf_values, facet_it, sdf); } } protected: - double calculate_sdf_value_of_facet(const Facet_handle& facet, const Tree& tree, - const Disk_samples_list& samples) const { + double calculate_sdf_value_of_facet(Facet_const_iterator& facet, + const Tree& tree, const Disk_samples_list& samples) const { const Point& p1 = facet->halfedge()->vertex()->point(); const Point& p2 = facet->halfedge()->next()->vertex()->point(); const Point& p3 = facet->halfedge()->prev()->vertex()->point(); @@ -116,10 +108,10 @@ protected: //arrange_center_orientation(plane, normal, center); std::vector ray_distances, ray_weights; - ray_distances.reserve(parameters.number_of_rays); - ray_weights.reserve(parameters.number_of_rays); + ray_distances.reserve(number_of_rays); + ray_weights.reserve(number_of_rays); - const double length_of_normal = 1.0 / tan(parameters.cone_angle / 2.0); + const double length_of_normal = 1.0 / tan(cone_angle / 2.0); normal = normal * length_of_normal; // stores segment length, // making it too large might cause a non-filtered bboxes in traversal, @@ -158,7 +150,7 @@ protected: ray_direction = ray_direction / std::sqrt(ray_direction.squared_length()); ray_direction = ray_direction * (*segment_distance * multiplier_for_segment); - Segment segment(center, CGAL::operator+(center, ray_direction)); + Segment segment(center, operator+(center, ray_direction)); boost::tie(is_intersected, intersection_is_acute, min_distance) = cast_and_return_minimum(segment, tree, facet); if(!is_intersected) { //no intersection is found @@ -193,8 +185,8 @@ protected: ray_distances.push_back(min_distance); } double sdf_result, acceptance_rate; - boost::tie(sdf_result, acceptance_rate) = calculate_sdf_value_from_rays( - ray_distances, ray_weights); + boost::tie(sdf_result, acceptance_rate) = + remove_outliers_and_calculate_sdf_value(ray_distances, ray_weights); if(acceptance_rate > CGAL_ACCEPTANCE_RATE_THRESHOLD || samples.size() == disk_samples_dense.size()) { @@ -205,7 +197,7 @@ protected: template //Query can be templated for just Ray and Segment types. boost::tuple cast_and_return_minimum( - const Query& query, const Tree& tree, const Facet_handle& facet) const { + const Query& query, const Tree& tree, Facet_const_iterator& facet) const { // get<0> : if any intersection is found then true // get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true // get<2> : distance between ray/segment origin and intersection point. @@ -226,14 +218,14 @@ protected: for(typename std::list::iterator op_it = intersections.begin(); op_it != intersections.end() ; ++op_it) { - CGAL::Object object = op_it->first; + Object object = op_it->first; Primitive_id id = op_it->second; if(id == facet) { continue; //Since center is located on related facet, we should skip it if there is an intersection with it. } const Point* i_point; - if(!(i_point = CGAL::object_cast(&object))) { + if(!(i_point = object_cast(&object))) { continue; // Continue in case of segment. } @@ -266,7 +258,7 @@ protected: template //Query can be templated for just Ray and Segment types. boost::tuple cast_and_return_minimum_use_closest ( const Query& ray, const Tree& tree, - const Facet_handle& facet) const { + Facet_const_iterator& facet) const { // get<0> : if any intersection is found then true // get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true // get<2> : distance between ray/segment origin and intersection point. @@ -285,7 +277,7 @@ protected: } min_distance.get<0>() = true; // intersection is found - CGAL::Object object = intersection->first; + Object object = intersection->first; Primitive_id min_id = intersection->second; if(min_id == facet) { @@ -293,7 +285,7 @@ protected: } const Point* i_point; - if(!(i_point = CGAL::object_cast(&object))) { + if(!(i_point = object_cast(&object))) { return min_distance; } @@ -310,7 +302,8 @@ protected: min_distance.get<2>() = std::sqrt(min_i_ray.squared_length()); return min_distance; } - boost::tuple calculate_sdf_value_from_rays( + + boost::tuple remove_outliers_and_calculate_sdf_value( std::vector& ray_distances, std::vector& ray_weights) const { // get<0> : sdf value @@ -357,8 +350,8 @@ protected: for(std::vector::iterator dist_it = ray_distances.begin(); dist_it != ray_distances.end(); ++dist_it, ++w_it) { - // AF: replace fabs with CGAL::abs - if(CGAL::abs((*dist_it) - median_sdf) > (st_dev * CGAL_ST_DEV_MULTIPLIER)) { + + if(abs((*dist_it) - median_sdf) > (st_dev * CGAL_ST_DEV_MULTIPLIER)) { continue; } total_distance += (*dist_it) * (*w_it); @@ -385,17 +378,17 @@ protected: //if(plane.has_on_positive_side(center)) { return; } //Vector epsilon_normal = unit_normal * epsilon; //do { - // center = CGAL::operator+(center, epsilon_normal); + // center = operator+(center, epsilon_normal); //} while(!plane.has_on_positive_side(center)); //Option-2 //if(plane.has_on_positive_side(center)) { return; } - //double distance = sqrt(CGAL::squared_distance(plane, center)); + //double distance = sqrt(squared_distance(plane, center)); //distance = distance > epsilon ? (distance + epsilon) : epsilon; //Vector distance_normal = unit_normal * distance; // //do { - // center = CGAL::operator+(center, distance_normal); + // center = operator+(center, distance_normal); //} while(!plane.has_on_positive_side(center)); //Option-3 @@ -407,14 +400,14 @@ protected: Vector epsilon_normal = unit_normal * epsilon; do { - center = CGAL::operator+(center, epsilon_normal); + center = operator+(center, epsilon_normal); ray = Ray(center, unit_normal); } while(intersector(ray, plane)); } void disk_sampling_vogel_method(Disk_samples_list& samples, int ray_count) { - const double length_of_normal = 1.0 / tan(parameters.cone_angle / 2.0); - const double angle_st_dev = parameters.cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER; + const double length_of_normal = 1.0 / tan(cone_angle / 2.0); + const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER; const double golden_ratio = 3.0 - std::sqrt(5.0); #if 0 @@ -440,6 +433,5 @@ protected: #undef CGAL_ANGLE_ST_DEV_DIVIDER #undef CGAL_ST_DEV_MULTIPLIER #undef CGAL_ACCEPTANCE_RATE_THRESHOLD -#undef CGAL_DEFAULT_CONE_ANGLE -#undef CGAL_DEFAULT_NUMBER_OF_RAYS -#endif //CGAL_SDF_CALCULATION \ No newline at end of file + +#endif //CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H \ No newline at end of file diff --git a/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h new file mode 100644 index 00000000000..4e573d86f14 --- /dev/null +++ b/Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h @@ -0,0 +1,38 @@ +#ifndef CGAL_SURFACE_MESH_SEGMENTATION_MESH_SEGMENTATION_H +#define CGAL_SURFACE_MESH_SEGMENTATION_MESH_SEGMENTATION_H + +#include + +#define CGAL_DEFAULT_CONE_ANGLE (2.0 / 3.0) * CGAL_PI +#define CGAL_DEFAULT_NUMBER_OF_RAYS 25 +#define CGAL_DEFAULT_NUMBER_OF_CLUSTERS 5 +#define CGAL_DEFAULT_SMOOTHING_LAMBDA 23.0 + +namespace CGAL +{ + +template +void sdf_values_computation(const Polyhedron& polyhedron, + SDFPropertyMap sdf_values, + double cone_angle = CGAL_DEFAULT_CONE_ANGLE, + int number_of_rays = CGAL_DEFAULT_NUMBER_OF_RAYS) +{ + Surface_mesh_segmentation algorithm(polyhedron); + algorithm.calculate_sdf_values(sdf_values, cone_angle, number_of_rays); +} + +template +void surface_mesh_segmentation_from_sdf_values(const Polyhedron& polyhedron, + SDFPropertyMap sdf_values, SegmentPropertyMap segment_ids, + int number_of_centers = CGAL_DEFAULT_NUMBER_OF_CLUSTERS, + double smoothing_lambda = CGAL_DEFAULT_SMOOTHING_LAMBDA) +{ + Surface_mesh_segmentation algorithm(polyhedron); + algorithm.partition(sdf_values, segment_ids, number_of_centers, + smoothing_lambda); +} +}//namespace CGAL + +#undef CGAL_DEFAULT_CONE_ANGLE +#undef CGAL_DEFAULT_NUMBER_OF_RAYS +#endif // CGAL_SURFACE_MESH_SEGMENTATION_MESH_SEGMENTATION_H // \ No newline at end of file