diff --git a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h index f21087e3ad2..f9039f4ffa7 100644 --- a/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h +++ b/Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h @@ -14,6 +14,20 @@ * +) Deciding how to generate rays in cone: for now using "polar angle" and "accept-reject (square)" and "concentric mapping" techniques */ + +//#include "Expectation_maximization.h" +//#include "K_means_clustering.h" +//#include "Timer.h" + +#include +#include + +#include +#include +#include +#include +#include + #include #include #include @@ -21,22 +35,17 @@ #include #include -//#include "Timer.h" -//#include "Expectation_maximization.h" -//#include "K_means_clustering.h" -#include -#include -#include - -#include -#include -#include -#include - #define LOG_5 1.60943791 #define NORMALIZATION_ALPHA 4.0 #define ANGLE_ST_DEV_DIVIDER 3.0 +//#define ACTIVATE_SEGMENTATION_DEBUG +#ifdef ACTIVATE_SEGMENTATION_DEBUG +#define SEG_DEBUG(x) x; +#else +#define SEG_DEBUG(x) +#endif + namespace CGAL { @@ -52,7 +61,9 @@ public: typedef typename Polyhedron::Facet_iterator Facet_iterator; typedef typename Polyhedron::Facet_handle Facet_handle; protected: - typedef typename Kernel::Ray_3 Ray; + typedef typename Kernel::Ray_3 Ray; + typedef typename Kernel::Plane_3 Plane; + typedef typename CGAL::AABB_polyhedron_triangle_primitive Primitive; typedef typename CGAL::AABB_tree > @@ -105,8 +116,7 @@ public: void calculate_sdf_values(); double calculate_sdf_value_of_facet (const Facet_handle& facet, - const Point& center, - const Vector& normal_const, const Tree& tree) const; + const Tree& tree) const; void cast_and_return_minimum (const Ray& ray, const Tree& tree, const Facet_handle& facet, bool& is_found, double& min_distance) const; @@ -118,6 +128,9 @@ public: double calculate_sdf_value_from_rays_with_trimmed_mean ( std::vector& ray_distances, std::vector& ray_weights) const; + void arrange_center_orientation(const Plane& plane, const Vector& unit_normal, + Point& center) const; + void disk_sampling_rejection(); void disk_sampling_polar_mapping(); void disk_sampling_concentric_mapping(); @@ -141,13 +154,14 @@ inline Surface_mesh_segmentation::Surface_mesh_segmentation( : mesh(mesh), cone_angle(cone_angle), number_of_rays_sqrt(number_of_rays_sqrt), number_of_centers(number_of_centers), log_file("log_file.txt") { - //Timer t; + SEG_DEBUG(Timer t) + disk_sampling_concentric_mapping(); calculate_sdf_values(); - //std::cout << t; - apply_GMM_fitting_with_K_means_init(); + SEG_DEBUG(std::cout << t) + apply_GMM_fitting(); //write_sdf_values("sdf_values_sample_dino_2.txt"); - //read_sdf_values("sdf_values_sample_camel.txt"); + //read_sdf_values("sdf_values_sample_dino.txt"); } template @@ -159,14 +173,8 @@ inline void Surface_mesh_segmentation::calculate_sdf_values() facet_it != mesh->facets_end(); ++facet_it) { CGAL_precondition(facet_it->is_triangle()); //Mesh should contain triangles. - Point v1 = facet_it->halfedge()->vertex()->point(); - Point v2 = facet_it->halfedge()->next()->vertex()->point(); - Point v3 = facet_it->halfedge()->next()->next()->vertex()->point(); - Point center = CGAL::centroid(v1, v2, v3); - Vector normal = CGAL::unit_normal(v1, v2, - v3) * -1.0; //Assuming triangles are CCW oriented. //SL: cone angle and number of rays should be parameters. - double sdf = calculate_sdf_value_of_facet(facet_it, center, normal, tree); + double sdf = calculate_sdf_value_of_facet(facet_it, tree); sdf_values.insert(std::pair(facet_it, sdf)); } normalize_sdf_values(); @@ -176,14 +184,22 @@ inline void Surface_mesh_segmentation::calculate_sdf_values() template inline double Surface_mesh_segmentation::calculate_sdf_value_of_facet( - const Facet_handle& facet, const Point& center, - const Vector& normal_const, const Tree& tree) const + const Facet_handle& facet, const Tree& tree) const { - typename Kernel::Plane_3 plane(center, normal_const); + Point p1 = facet->halfedge()->vertex()->point(); + Point p2 = facet->halfedge()->next()->vertex()->point(); + Point p3 = facet->halfedge()->next()->next()->vertex()->point(); + Point center = CGAL::centroid(p1, p2, p3); + Vector normal = CGAL::unit_normal(p1, p2, + p3) * -1.0; //Assuming triangles are CCW oriented. + + Plane plane(center, normal); Vector v1 = plane.base1(); Vector v2 = plane.base2(); - v1 = v1 / CGAL::sqrt(v1.squared_length()); - v2 = v2 / CGAL::sqrt(v2.squared_length()); + v1 = v1 / sqrt(v1.squared_length()); + v2 = v2 / sqrt(v2.squared_length()); + + arrange_center_orientation(plane, normal, center); int ray_count = number_of_rays_sqrt * number_of_rays_sqrt; @@ -192,7 +208,7 @@ Surface_mesh_segmentation::calculate_sdf_value_of_facet( ray_weights.reserve(ray_count); double length_of_normal = 1.0 / tan(cone_angle / 2); - Vector normal = normal_const * length_of_normal; + normal = normal * length_of_normal; for(Disk_samples_list::const_iterator sample_it = disk_samples.begin(); sample_it != disk_samples.end(); ++sample_it) { @@ -256,6 +272,8 @@ void Surface_mesh_segmentation::cast_and_return_minimum( CGAL::ORIGIN + min_normal) != CGAL::ACUTE) { is_found = false; } + //total_intersection += intersections.size(); + //total_cast++; } template @@ -403,6 +421,49 @@ Surface_mesh_segmentation::calculate_sdf_value_from_rays_with_trimme return total_distance / total_weights; } +/* Slightly move center towards inverse normal of facet. + * Parameter plane is constructed by inverse normal. + */ +template +void Surface_mesh_segmentation::arrange_center_orientation( + const Plane& plane, const Vector& unit_normal, Point& center) const +{ + /* + * Not sure how to decide on how much I should move center ? + */ + double epsilon = 1e-8; + // Option-1 + //if(plane.has_on_positive_side(center)) { return; } + //Vector epsilon_normal = unit_normal * epsilon; + //do { + // center = CGAL::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)); + //distance = distance > epsilon ? (distance + epsilon) : epsilon; + //Vector distance_normal = unit_normal * distance; + // + //do { + // center = CGAL::operator+(center, distance_normal); + //} while(!plane.has_on_positive_side(center)); + + //Option-3 + Ray ray(center, unit_normal); + Kernel::Do_intersect_3 intersector = Kernel().do_intersect_3_object(); + if(!intersector(ray, plane)) { + return; + } + + Vector epsilon_normal = unit_normal * epsilon; + do { + center = CGAL::operator+(center, epsilon_normal); + ray = Ray(center, unit_normal); + } while(intersector(ray, plane)); + +} + template inline void Surface_mesh_segmentation::disk_sampling_rejection() { @@ -522,7 +583,7 @@ inline void Surface_mesh_segmentation::smooth_sdf_values() total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; } while( ++facet_circulator != f->facet_begin()); total_neighbor_sdf /= 3; - smoothed_sdf_values[f] = (sdf_values[f] + total_neighbor_sdf) / 2; + smoothed_sdf_values[f] = (pair_it->second + total_neighbor_sdf) / 2; } sdf_values = smoothed_sdf_values; } @@ -610,6 +671,7 @@ inline void Surface_mesh_segmentation::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) { double sdf_value; @@ -640,4 +702,8 @@ inline void Surface_mesh_segmentation::read_center_ids( #undef ANGLE_ST_DEV_DIVIDER #undef LOG_5 #undef NORMALIZATION_ALPHA + +#ifdef SEG_DEBUG +#undef SEG_DEBUG +#endif #endif //CGAL_SURFACE_MESH_SEGMENTATION_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 dc0c5537a6c..34b1ff63dc9 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 @@ -2,7 +2,7 @@ #define CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H /* NEED TO BE DONE */ /* About implementation: - * Calculating probability multiple times + * Calculating probability multiple times, solution: storing matrix(cluster, point) for probability values. * About safe division: where centers have no members, in graph cut it may happen */ #include @@ -12,6 +12,15 @@ #include #include +#include +#include +//#define ACTIVATE_SEGMENTATION_EM_DEBUG +#ifdef ACTIVATE_SEGMENTATION_EM_DEBUG +#define SEG_DEBUG(x) x; +#else +#define SEG_DEBUG(x) +#endif + namespace CGAL { @@ -102,13 +111,16 @@ protected: public: Expectation_maximization(int number_of_centers, const std::vector& data, const std::vector& initial_centers = std::vector(), - int maximum_iteration = 100) - : points(data.begin(), data.end()), threshold(1e-3), + int number_of_runs = 50, int maximum_iteration = 100) + : points(data.begin(), data.end()), threshold(1e-4), maximum_iteration(maximum_iteration), is_converged(false), seed(time(NULL)) { srand(seed); - //calculate_fitting_with_multiple_run(number_of_centers, 50); - initiate_centers(number_of_centers, initial_centers); - calculate_fitting(); + if(initial_centers.empty()) { + calculate_fitting_with_multiple_run(number_of_centers, number_of_runs); + } else { + initiate_centers(number_of_centers, initial_centers); + calculate_fitting(); + } } void fill_with_center_ids(std::vector& data_centers) { @@ -133,21 +145,26 @@ protected: void initiate_centers_randomly(int number_of_centers) { centers.clear(); /* Randomly generate means of centers */ + //vector center_indexes; + double initial_mixing_coefficient = 1.0; double initial_deviation = 1.0 / (2.0 * number_of_centers); - double initial_mixing_coefficient = 1.0 / number_of_centers; for(int i = 0; i < number_of_centers; ++i) { - double initial_mean = points[rand() % points.size()].data; + double random_index = rand() % points.size(); + //center_indexes.push_back(random_index); + Gaussian_point mean_point = points[random_index]; + double initial_mean = mean_point.data; centers.push_back(Gaussian_center(initial_mean, initial_deviation, initial_mixing_coefficient)); } sort(centers.begin(), centers.end()); + //write_random_centers("center_indexes.txt", center_indexes); } void initiate_centers_uniformly(int number_of_centers) { centers.clear(); /* Uniformly generate centers */ double initial_deviation = 1.0 / (2.0 * number_of_centers); - double initial_mixing_coefficient = 1.0 / number_of_centers; + double initial_mixing_coefficient = 1.0; for(int i = 0; i < number_of_centers; ++i) { double initial_mean = (i + 1.0) / (number_of_centers + 1.0); centers.push_back(Gaussian_center(initial_mean, initial_deviation, @@ -203,6 +220,14 @@ protected: point_it->calculate_total_membership(centers); } } + double calculate_deviation(const Gaussian_point& point) const { + double deviation = 0.0; + for(std::vector::const_iterator point_it = points.begin(); + point_it != points.end(); ++point_it) { + deviation += pow(point_it->data - point.data, 2); + } + return sqrt(deviation / points.size()); + } /*Calculates new parameter values for each cluster */ void calculate_parameters() { for(std::vector::iterator center_it = centers.begin(); @@ -239,9 +264,9 @@ protected: while(!is_converged && iteration_count++ < maximum_iteration) { prev_likelihood = likelihood; likelihood = iterate(); - is_converged = likelihood - prev_likelihood < threshold * likelihood; + is_converged = likelihood - prev_likelihood < threshold * fabs(likelihood); } - //std::cout << likelihood << " " << iteration_count << std::endl; + SEG_DEBUG(std::cout << likelihood << " " << iteration_count << std::endl) return likelihood; } @@ -254,6 +279,7 @@ protected: initiate_centers_randomly(number_of_centers); double likelihood = calculate_fitting(); + write_center_parameters("center_paramters.txt"); if(likelihood > max_likelihood) { max_centers = centers; max_likelihood = likelihood; @@ -261,6 +287,28 @@ protected: } centers = max_centers; } + + void write_random_centers(const char* file_name, + const std::vector& center_indexes) { + std::ofstream output(file_name, std::ios_base::app); + for(std::vector::const_iterator it = center_indexes.begin(); + it != center_indexes.end(); ++it) { + output << points[(*it)].data << std::endl; + } + output.close(); + } + + void write_center_parameters(const char* file_name) { + std::ofstream output(file_name, std::ios_base::app); + for(std::vector::iterator center_it = centers.begin(); + center_it != centers.end(); ++center_it) { + output << "mean: " << center_it->mean << " dev: " << center_it->deviation + << " mix: " << center_it->mixing_coefficient << std::endl; + } + } }; }//namespace CGAL +#ifdef SEG_DEBUG +#undef SEG_DEBUG +#endif #endif //CGAL_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 f927b285c3c..77e3f8de740 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 @@ -7,6 +7,13 @@ #include #include +//#define ACTIVATE_SEGMENTATION_K_MEANS_DEBUG +#ifdef ACTIVATE_SEGMENTATION_K_MEANS_DEBUG +#define SEG_DEBUG(x) x; +#else +#define SEG_DEBUG(x) +#endif + namespace CGAL { @@ -84,8 +91,9 @@ public: srand(seed); calculate_clustering_with_multiple_run(number_of_centers, number_of_run); } - + /* For each data point, data_center is filled by its center's id. */ void fill_with_center_ids(std::vector& data_centers) { + data_centers.reserve(points.size()); for(std::vector::iterator point_it = points.begin(); point_it != points.end(); ++point_it) { data_centers.push_back(point_it->center_id); @@ -115,18 +123,21 @@ public: bool any_center_changed = true; while(any_center_changed && iteration_count++ < maximum_iteration) { any_center_changed = false; + /* For each point, calculate its new center */ for(std::vector::iterator point_it = points.begin(); point_it != points.end(); ++point_it) { bool center_changed = point_it->calculate_new_center(centers); any_center_changed |= center_changed; } + /* For each center, calculate its new mean */ for(std::vector::iterator center_it = centers.begin(); center_it != centers.end(); ++center_it) { center_it->calculate_mean(); } } is_converged = !any_center_changed; - //std::cout << iteration_count << " " << is_converged << std::endl; + SEG_DEBUG(std::cout << iteration_count << " " << (is_converged ? "converged" : + "not converged") << std::endl) } void calculate_clustering_with_multiple_run(int number_of_centers, @@ -170,4 +181,7 @@ public: } }; }//namespace CGAL +#ifdef SEG_DEBUG +#undef SEG_DEBUG +#endif #endif //CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H