songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Using FacetIndexMap for in Surface_mesh_segmentation. 

(other structures which have a facet as key-value are converted to vectors)
Using const Polyhedron& as a parameter (previously I was using Polyhedron pointer).

mesh_segmentation.h added for API which includes free functions.

AABB_const_polyhedron_triangle_primitive.h is added for supporting Facet_const_iterator as a parameter to AABB Tree. (I am not sure that is the best solution.)
This commit is contained in:
Ílker Yaz 2012-07-31 01:01:18 +00:00
parent d2b5143ab7
commit 792ab1f525
8 changed files with 542 additions and 469 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitattributes vendored
View File

@ -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/segmentation_example_sdf.cpp -text
Surface_mesh_segmentation/example/Surface_mesh_segmentation/simple_example.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/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/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/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/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/K_means_clustering.h -text
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.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/copyright -text
Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/license.txt -text Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/license.txt -text
Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/maintainer -text Surface_mesh_segmentation/package_info/Surface_mesh_segmentation/maintainer -text

667
Surface_mesh_segmentation/include/CGAL/Surface_mesh_segmentation.h
View File

@ -24,6 +24,7 @@
#include <CGAL/Mesh_3/dihedral_angle_3.h> #include <CGAL/Mesh_3/dihedral_angle_3.h>
#include <boost/optional.hpp> #include <boost/optional.hpp>
#include <boost/property_map/property_map.hpp>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
@ -32,12 +33,15 @@
#include <algorithm> #include <algorithm>
#include <utility> #include <utility>
#include <queue> #include <queue>
#include <map>
//AF: macros must be prefixed with "CGAL_" //AF: macros must be prefixed with "CGAL_"
//IOY: done //IOY: done
#define CGAL_NORMALIZATION_ALPHA 5.0 #define CGAL_NORMALIZATION_ALPHA 5.0
#define CGAL_CONVEX_FACTOR 0.08 #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) //IOY: these are going to be removed at the end (no CGAL_ pref)
#define ACTIVATE_SEGMENTATION_DEBUG #define ACTIVATE_SEGMENTATION_DEBUG
#ifdef 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. * 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. * All preprocessing and postprocessing issues are handled here.
*/ */
template <class Polyhedron, class SDFCalculation = internal::SDF_calculation<Polyhedron> > template <
class Polyhedron,
class FacetIndexMap =
boost::associative_property_map<std::map<typename Polyhedron::Facet_const_iterator, int> >
>
class Surface_mesh_segmentation class Surface_mesh_segmentation
{ {
//type definitions //type definitions
@ -65,81 +73,110 @@ public:
typedef typename Polyhedron::Facet Vertex; typedef typename Polyhedron::Facet Vertex;
typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::Point_3 Point; 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<Polyhedron> SDF_calculation;
protected: protected:
typedef typename Kernel::Plane_3 Plane; typedef typename Kernel::Plane_3 Plane;
typedef std::map<Facet_handle, double> Face_value_map; // member variables
typedef std::map<Facet_handle, int> Face_center_map; public: // going to be protected !
typedef std::map<Facet_handle, int> Face_segment_map; const Polyhedron& mesh;
template <typename ValueTypeName> FacetIndexMap facet_index_map;
struct Compare_second_element {
bool operator()(const ValueTypeName& v1, const ValueTypeName& v2) const {
return v1.second < v2.second;
}
};
//member variables std::vector<double> sdf_values;
public: std::vector<int> centers;
Polyhedron* mesh; std::vector<int> segments;
Face_value_map sdf_values;
Face_center_map centers;
Face_segment_map segments;
int number_of_centers; int number_of_centers;
double smoothing_lambda; double smoothing_lambda;
std::map<Facet_const_iterator, int> facet_index_map_internal;
internal::Expectation_maximization fitter;/**< going to be removed */ // member functions
//member functions
public: public:
Surface_mesh_segmentation(Polyhedron* mesh): mesh(mesh) { Surface_mesh_segmentation(const Polyhedron& mesh)
#ifdef SEGMENTATION_PROFILE : mesh(mesh), facet_index_map(facet_index_map_internal) {
profile("profile.txt"); int facet_index = 0;
#endif 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()) { Surface_mesh_segmentation(Polyhedron* mesh, FacetIndexMap facet_index_map)
SEG_DEBUG(CGAL::Timer t) : 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<Facet_const_iterator, double> internal_map;
internal_map facet_value_map_internal;
boost::associative_property_map<internal_map> sdf_pmap(
facet_value_map_internal);
SDF_calculation(cone_angle, number_of_ray).calculate_sdf_values(mesh, sdf_pmap);
sdf_values = std::vector<double>(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 <class SDFPropertyMap>
void calculate_sdf_values(SDFPropertyMap sdf_pmap, double cone_angle,
int number_of_ray) {
SEG_DEBUG(Timer t)
SEG_DEBUG(t.start()) SEG_DEBUG(t.start())
sdf_values.clear(); SDF_calculation(cone_angle, number_of_rays).calculate_sdf_values(mesh,
SDFCalculation(parameters).calculate_sdf_values(*mesh, sdf_values); sdf_pmap);
sdf_values = std::vector<double>(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) SEG_DEBUG(std::cout << t.time() << std::endl)
check_zero_sdf_values(); check_zero_sdf_values();
smooth_sdf_values_with_bilateral(); 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) SEG_DEBUG(std::cout << t.time() << std::endl)
} }
void partition(int number_of_centers = 5, double smoothing_lambda = 23.0) { template <class FacetSegmentMap>
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->number_of_centers = number_of_centers;
this->smoothing_lambda = smoothing_lambda; this->smoothing_lambda = smoothing_lambda;
centers.clear(); centers.clear();
// log normalize sdf values
std::vector<double> sdf_vector; normalize_sdf_values();
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]);
}
// soft clustering using GMM-fitting initialized with k-means // soft clustering using GMM-fitting initialized with k-means
fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, internal::Expectation_maximization fitter(number_of_centers, sdf_values,
internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1); internal::Expectation_maximization::K_MEANS_INITIALIZATION, 1);
std::vector<int> labels; std::vector<int> labels;
fitter.fill_with_center_ids(labels); fitter.fill_with_center_ids(labels);
@ -156,21 +193,62 @@ public:
// apply graph cut // apply graph cut
internal::Alpha_expansion_graph_cut gc(edges, edge_weights, probability_matrix, internal::Alpha_expansion_graph_cut gc(edges, edge_weights, probability_matrix,
labels); labels);
centers = labels;
std::vector<int>::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_handle, int>(facet_it, (*center_it)));
}
// assign a segment id for each facet // assign a segment id for each facet
assign_segments(); 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: template <class FacetSegmentMap, class SDFPropertyMap>
double calculate_dihedral_angle_of_edge(const Halfedge_handle& edge) const { void partition(SDFPropertyMap sdf_pmap, FacetSegmentMap segment_pmap,
Facet_handle f1 = edge->facet(); int number_of_centers = CGAL_DEFAULT_NUMBER_OF_CLUSTERS,
Facet_handle f2 = edge->opposite()->facet(); double smoothing_lambda = CGAL_DEFAULT_SMOOTHING_LAMBDA) {
sdf_values = std::vector<double>(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<int> labels;
fitter.fill_with_center_ids(labels);
std::vector<std::vector<double> > probability_matrix;
fitter.fill_with_probabilities(probability_matrix);
log_normalize_probability_matrix(probability_matrix);
// calculating edge weights
std::vector<std::pair<int, int> > edges;
std::vector<double> 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_v1 = f2->halfedge()->vertex()->point();
const Point& f2_v2 = f2->halfedge()->next()->vertex()->point(); const Point& f2_v2 = f2->halfedge()->next()->vertex()->point();
@ -204,14 +282,14 @@ public:
return angle; 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& a = edge->vertex()->point();
const Point& b = edge->prev()->vertex()->point(); const Point& b = edge->prev()->vertex()->point();
const Point& c = edge->next()->vertex()->point(); const Point& c = edge->next()->vertex()->point();
const Point& d = edge->opposite()->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, // As far as I check: if, say, dihedral angle is 5, this returns 175,
// if 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; bool concave = n_angle > 0;
double angle = 1 + ((concave ? -1 : +1) * n_angle); double angle = 1 + ((concave ? -1 : +1) * n_angle);
@ -220,8 +298,8 @@ public:
} }
return angle; return angle;
//Facet_handle f1 = edge->facet(); //Facet_const_iterator f1 = edge->facet();
//Facet_handle f2 = edge->opposite()->facet(); //Facet_const_iterator f2 = edge->opposite()->facet();
// //
//const Point& f2_v1 = f2->halfedge()->vertex()->point(); //const Point& f2_v1 = f2->halfedge()->vertex()->point();
//const Point& f2_v2 = f2->halfedge()->next()->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_v1 = f1->halfedge()->vertex()->point();
//const Point& f1_v2 = f1->halfedge()->next()->vertex()->point(); //const Point& f1_v2 = f1->halfedge()->next()->vertex()->point();
//const Point& f1_v3 = f1->halfedge()->prev()->vertex()->point(); //const Point& f1_v3 = f1->halfedge()->prev()->vertex()->point();
//Vector f1_normal = CGAL::unit_normal(f1_v1, f1_v2, f1_v3); //Vector f1_normal = unit_normal(f1_v1, f1_v2, f1_v3);
//Vector f2_normal = CGAL::unit_normal(f2_v1, f2_v2, f2_v3); //Vector f2_normal = unit_normal(f2_v1, f2_v2, f2_v3);
// //
//double dot = f1_normal * f2_normal; //double dot = f1_normal * f2_normal;
//if(dot > 1.0) { dot = 1.0; } //if(dot > 1.0) { dot = 1.0; }
@ -257,37 +335,46 @@ public:
} }
void normalize_sdf_values() { void normalize_sdf_values() {
typedef typename Face_value_map::iterator fv_iterator; typedef std::vector<double>::iterator fv_iterator;
Compare_second_element<typename Face_value_map::value_type> comparator;
std::pair<fv_iterator, fv_iterator> min_max_pair = std::pair<fv_iterator, fv_iterator> min_max_pair =
CGAL::min_max_element(sdf_values.begin(), sdf_values.end(), comparator, min_max_element(sdf_values.begin(), sdf_values.end());
comparator);
double max_value = min_max_pair.second->second, double max_value = *min_max_pair.second, min_value = *min_max_pair.first;
min_value = min_max_pair.first->second;
double max_min_dif = max_value - min_value; double max_min_dif = max_value - min_value;
for(fv_iterator pair_it = sdf_values.begin(); pair_it != sdf_values.end(); for(fv_iterator it = sdf_values.begin(); it != sdf_values.end(); ++it) {
++pair_it) { double linear_normalized = (*it - min_value) / max_min_dif;
double linear_normalized = (pair_it->second - min_value) / max_min_dif;
double log_normalized = log(linear_normalized * CGAL_NORMALIZATION_ALPHA + 1) / double log_normalized = log(linear_normalized * CGAL_NORMALIZATION_ALPHA + 1) /
log(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<double>::iterator fv_iterator;
std::pair<fv_iterator, fv_iterator> 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() { void smooth_sdf_values() {
Face_value_map smoothed_sdf_values; std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
for(typename Face_value_map::iterator pair_it = sdf_values.begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
pair_it != sdf_values.end(); ++pair_it) { facet_it != mesh.facets_end(); ++facet_it) {
Facet_handle f = pair_it->first; typename Facet::Halfedge_around_facet_const_circulator facet_circulator =
typename Facet::Halfedge_around_facet_circulator facet_circulator = facet_it->facet_begin();
f->facet_begin();
double total_neighbor_sdf = 0.0; double total_neighbor_sdf = 0.0;
do { do {
total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; total_neighbor_sdf += get(sdf_values, facet_circulator->opposite()->facet());
} while( ++facet_circulator != f->facet_begin()); } while( ++facet_circulator != facet_it->facet_begin());
total_neighbor_sdf /= 3.0; 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; sdf_values = smoothed_sdf_values;
} }
@ -298,23 +385,22 @@ public:
const int iteration = 1; const int iteration = 1;
for(int i = 0; i < iteration; ++i) { for(int i = 0; i < iteration; ++i) {
Face_value_map smoothed_sdf_values; std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
for(typename Face_value_map::iterator pair_it = sdf_values.begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
pair_it != sdf_values.end(); ++pair_it) { facet_it != mesh.facets_end(); ++facet_it) {
Facet_handle facet = pair_it->first; std::map<Facet_const_iterator, int> neighbors;
std::map<Facet_handle, int> neighbors; get_neighbors_by_vertex(facet_it, neighbors, window_size);
get_neighbors_by_vertex(facet, neighbors, window_size);
double total_sdf_value = 0.0; double total_sdf_value = 0.0;
double total_weight = 0.0; double total_weight = 0.0;
for(typename std::map<Facet_handle, int>::iterator it = neighbors.begin(); for(typename std::map<Facet_const_iterator, int>::iterator it =
it != neighbors.end(); ++it) { neighbors.begin(); it != neighbors.end(); ++it) {
double weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0), double weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0),
2))); // window_size => 2*sigma 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; 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; sdf_values = smoothed_sdf_values;
} }
@ -326,17 +412,16 @@ public:
const int iteration = 1; const int iteration = 1;
for(int i = 0; i < iteration; ++i) { for(int i = 0; i < iteration; ++i) {
Face_value_map smoothed_sdf_values; std::vector<double> smoothed_sdf_values(mesh.size_of_facets());;
for(typename Face_value_map::iterator pair_it = sdf_values.begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
pair_it != sdf_values.end(); ++pair_it) { facet_it != mesh.facets_end(); ++facet_it) {
//Find neighbors and put their sdf values into a list //Find neighbors and put their sdf values into a list
Facet_handle facet = pair_it->first; std::map<Facet_const_iterator, int> neighbors;
std::map<Facet_handle, int> neighbors; get_neighbors_by_vertex(facet_it, neighbors, window_size);
get_neighbors_by_vertex(facet, neighbors, window_size);
std::vector<double> sdf_of_neighbors; std::vector<double> sdf_of_neighbors;
for(typename std::map<Facet_handle, int>::iterator it = neighbors.begin(); for(typename std::map<Facet_const_iterator, int>::iterator it =
it != neighbors.end(); ++it) { neighbors.begin(); it != neighbors.end(); ++it) {
sdf_of_neighbors.push_back(sdf_values[it->first]); sdf_of_neighbors.push_back(get(sdf_values, it->first));
} }
// Find median. // Find median.
double median_sdf = 0.0; double median_sdf = 0.0;
@ -351,7 +436,7 @@ public:
} else { } else {
median_sdf = sdf_of_neighbors[half_neighbor_count]; 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; sdf_values = smoothed_sdf_values;
} }
@ -366,51 +451,51 @@ public:
const int iteration = 1; const int iteration = 1;
for(int i = 0; i < iteration; ++i) { for(int i = 0; i < iteration; ++i) {
Face_value_map smoothed_sdf_values; std::vector<double> smoothed_sdf_values(mesh.size_of_facets());
for(typename Face_value_map::iterator pair_it = sdf_values.begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
pair_it != sdf_values.end(); ++pair_it) { facet_it != mesh.facets_end(); ++facet_it) {
Facet_handle facet = pair_it->first; //Facet_handle facet = facet_it;
std::map<Facet_handle, int> neighbors; std::map<Facet_const_iterator, int> neighbors;
get_neighbors_by_vertex(facet, neighbors, window_size); get_neighbors_by_vertex(facet_it, neighbors, window_size);
double total_sdf_value = 0.0, total_weight = 0.0; 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. // calculate deviation for range weighting.
double deviation = 0.0; double deviation = 0.0;
for(typename std::map<Facet_handle, int>::iterator it = neighbors.begin(); for(typename std::map<Facet_const_iterator, int>::iterator it =
it != neighbors.end(); ++it) { neighbors.begin(); it != neighbors.end(); ++it) {
deviation += std::pow(sdf_values[it->first] - current_sdf_value, 2); deviation += std::pow(get(sdf_values, it->first) - current_sdf_value, 2);
} }
deviation = std::sqrt(deviation / neighbors.size()); deviation = std::sqrt(deviation / neighbors.size());
if(deviation == 0.0) { if(deviation == 0.0) {
deviation = std::numeric_limits<double>::epsilon(); //this might happen deviation = std::numeric_limits<double>::epsilon(); //this might happen
} }
for(typename std::map<Facet_handle, int>::iterator it = neighbors.begin(); for(std::map<Facet_const_iterator, int>::iterator it = neighbors.begin();
it != neighbors.end(); ++it) { it != neighbors.end(); ++it) {
double spatial_weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0), double spatial_weight = exp(-0.5 * (std::pow(it->second / (window_size/2.0),
2))); // window_size => 2*sigma 2))); // window_size => 2*sigma
double domain_weight = exp(-0.5 * (std::pow((sdf_values[it->first] - double domain_weight = exp(-0.5 * (std::pow( (get(sdf_values,
current_sdf_value) / (std::sqrt(2.0)*deviation), 2))); it->first) - current_sdf_value) / (std::sqrt(2.0)*deviation), 2)));
double weight = spatial_weight * domain_weight; 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; 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; sdf_values = smoothed_sdf_values;
} }
} }
void get_neighbors_by_edge(const Facet_handle& facet, void get_neighbors_by_edge(Facet_const_iterator& facet,
std::map<Facet_handle, int>& neighbors, int max_level) { std::map<Facet_const_iterator, int>& neighbors, int max_level) {
typedef std::pair<Facet_handle, int> facet_level_pair; typedef std::pair<Facet_const_iterator, int> facet_level_pair;
std::queue<facet_level_pair> facet_queue; std::queue<facet_level_pair> facet_queue;
facet_queue.push(facet_level_pair(facet, 0)); facet_queue.push(facet_level_pair(facet, 0));
while(!facet_queue.empty()) { while(!facet_queue.empty()) {
const facet_level_pair& pair = facet_queue.front(); const facet_level_pair& pair = facet_queue.front();
bool inserted = neighbors.insert(pair).second; bool inserted = neighbors.insert(pair).second;
if(inserted && pair.second < max_level) { 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(); pair.first->facet_begin();
do { do {
facet_queue.push(facet_level_pair(facet_circulator->opposite()->facet(), facet_queue.push(facet_level_pair(facet_circulator->opposite()->facet(),
@ -421,20 +506,20 @@ public:
} }
} }
void get_neighbors_by_vertex(const Facet_handle& facet, void get_neighbors_by_vertex(Facet_const_iterator& facet,
std::map<Facet_handle, int>& neighbors, int max_level) { std::map<Facet_const_iterator, int>& neighbors, int max_level) {
typedef std::pair<Facet_handle, int> facet_level_pair; typedef std::pair<Facet_const_iterator, int> facet_level_pair;
std::queue<facet_level_pair> facet_queue; std::queue<facet_level_pair> facet_queue;
facet_queue.push(facet_level_pair(facet, 0)); facet_queue.push(facet_level_pair(facet, 0));
while(!facet_queue.empty()) { while(!facet_queue.empty()) {
const facet_level_pair& pair = facet_queue.front(); const facet_level_pair& pair = facet_queue.front();
bool inserted = neighbors.insert(pair).second; bool inserted = neighbors.insert(pair).second;
if(inserted && pair.second < max_level) { if(inserted && pair.second < max_level) {
Facet_handle facet = pair.first; Facet_const_iterator facet = pair.first;
Halfedge_handle edge = facet->halfedge(); Halfedge_const_iterator edge = facet->halfedge();
do { // loop on three vertices of the facet do { // loop on three vertices of the facet
Vertex_handle vertex = edge->vertex(); Vertex_const_iterator vertex = edge->vertex();
typename Facet::Halfedge_around_vertex_circulator vertex_circulator = typename Facet::Halfedge_around_vertex_const_circulator vertex_circulator =
vertex->vertex_begin(); vertex->vertex_begin();
do { // for each vertex loop on neighbor vertices do { // for each vertex loop on neighbor vertices
facet_queue.push(facet_level_pair(vertex_circulator->facet(), pair.second + 1)); facet_queue.push(facet_level_pair(vertex_circulator->facet(), pair.second + 1));
@ -444,34 +529,35 @@ public:
facet_queue.pop(); facet_queue.pop();
} }
} }
void check_zero_sdf_values() { void check_zero_sdf_values() {
// If there is any facet which has no sdf value, assign average sdf value of its neighbors // 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(); for(Facet_const_iterator facet_it = mesh.facets_begin();
pair_it != sdf_values.end(); ++pair_it) { facet_it != mesh.facets_end(); ++facet_it) {
if(pair_it->second == 0.0) { if(get(sdf_values, facet_it) == 0.0) {
typename Facet::Halfedge_around_facet_circulator facet_circulator = typename Facet::Halfedge_around_facet_const_circulator facet_circulator =
pair_it->first->facet_begin(); facet_it->facet_begin();
double total_neighbor_sdf = 0.0; double total_neighbor_sdf = 0.0;
do { do {
total_neighbor_sdf += sdf_values[facet_circulator->opposite()->facet()]; total_neighbor_sdf += get(sdf_values, facet_circulator->opposite()->facet());
} while( ++facet_circulator != pair_it->first->facet_begin()); } while( ++facet_circulator != facet_it->facet_begin());
pair_it->second = total_neighbor_sdf / 3.0; get(sdf_values, facet_it) = total_neighbor_sdf / 3.0;
} }
} }
// Find minimum sdf value other than 0 // Find minimum sdf value other than 0
double min_sdf = (std::numeric_limits<double>::max)(); double min_sdf = (std::numeric_limits<double>::max)();
for(typename Face_value_map::iterator pair_it = sdf_values.begin(); for(std::vector<double>::iterator it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it) { it != sdf_values.end(); ++it) {
if(pair_it->second < min_sdf && pair_it->second != 0.0) { if(*it < min_sdf && *it != 0.0) {
min_sdf = pair_it->second; min_sdf = *it;
} }
} }
// If still there is any facet which has no sdf value, assign minimum sdf value. // 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. // 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(); for(std::vector<double>::iterator it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it) { it != sdf_values.end(); ++it) {
if(pair_it->second == 0.0) { if(*it == 0.0) {
pair_it->second = min_sdf; *it = min_sdf;
} }
} }
} }
@ -495,19 +581,11 @@ public:
std::vector<std::pair<int, int> >& edges, std::vector<std::pair<int, int> >& edges,
std::vector<double>& edge_weights) const { std::vector<double>& edge_weights) const {
const double epsilon = 1e-5; const double epsilon = 1e-5;
//assign an id for every facet (facet-id) //edges and their weights. pair<int, int> stores facet-id pairs (see above) (may be using boost::tuple can be more suitable)
std::map<Facet_handle, int> facet_indices; for(Edge_const_iterator edge_it = mesh.edges_begin();
int index = 0; edge_it != mesh.edges_end(); ++edge_it) {
for(Facet_iterator facet_it = mesh->facets_begin(); int index_f1 = boost::get(facet_index_map, edge_it->facet());
facet_it != mesh->facets_end(); int index_f2 = boost::get(facet_index_map, edge_it->opposite()->facet());
++facet_it, ++index) {
facet_indices.insert(std::pair<Facet_handle, int>(facet_it, index));
}
//edges and their weights. pair<int, int> 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.push_back(std::pair<int, int>(index_f1, index_f2)); edges.push_back(std::pair<int, int>(index_f1, index_f2));
double angle = calculate_dihedral_angle_of_edge(edge_it); double angle = calculate_dihedral_angle_of_edge(edge_it);
@ -522,240 +600,47 @@ public:
} }
void assign_segments() { void assign_segments() {
segments.clear(); segments = std::vector<int>(mesh.size_of_facets(), -1);
for(Facet_iterator facet_it = mesh->facets_begin();
facet_it != mesh->facets_end(); ++facet_it) {
segments.insert(std::pair<Facet_handle, int>(facet_it, -1));
}
int segment_id = 0; int segment_id = 0;
for(Facet_iterator facet_it = mesh->facets_begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
facet_it != mesh->facets_end(); ++facet_it) { facet_it != mesh.facets_end(); ++facet_it) {
if(segments[facet_it] == -1) { if(get(segments, facet_it) == -1) {
depth_first_traversal(facet_it, segment_id); depth_first_traversal(facet_it, segment_id);
segment_id++; segment_id++;
} }
} }
} }
void depth_first_traversal(const Facet_handle& facet, int segment_id) { void depth_first_traversal(Facet_const_iterator& facet, int segment_id) {
segments[facet] = segment_id; get(segments, facet) = segment_id;
typename Facet::Halfedge_around_facet_circulator facet_circulator = typename Facet::Halfedge_around_facet_const_circulator facet_circulator =
facet->facet_begin(); facet->facet_begin();
double total_neighbor_sdf = 0.0; double total_neighbor_sdf = 0.0;
do { do {
Facet_handle neighbor = facet_circulator->opposite()->facet(); Facet_const_iterator neighbor = facet_circulator->opposite()->facet();
if(segments[neighbor] == -1 && centers[facet] == centers[neighbor]) { if(get(segments, neighbor) == -1
&& get(centers, facet) == get(centers, neighbor)) {
depth_first_traversal(neighbor, segment_id); depth_first_traversal(neighbor, segment_id);
} }
} while( ++facet_circulator != facet->facet_begin()); } while( ++facet_circulator != facet->facet_begin());
} }
/**
* Going to be removed
*/
void apply_GMM_fitting() {
centers.clear();
std::vector<double> 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<int> center_memberships;
fitter.fill_with_center_ids(center_memberships);
std::vector<int>::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_handle, int>(facet_it, (*center_it)));
}
}
/**
* Going to be removed
*/
void apply_K_means_clustering() {
centers.clear();
std::vector<double> 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<int> center_memberships;
clusterer.fill_with_center_ids(center_memberships);
std::vector<int>::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_handle, int>(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<double> 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<int> 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<int>::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_handle, int>(facet_it, (*center_it)));
}
}
/**
* Going to be removed
*/
void apply_GMM_fitting_and_K_means() {
centers.clear();
std::vector<double> 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<int> 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<int>::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_handle, int>(facet_it, (*center_it)));
}
}
/**
* Going to be removed
*/
void apply_graph_cut() {
std::vector<std::pair<int, int> > edges;
std::vector<double> edge_weights;
calculate_and_log_normalize_dihedral_angles(edges, edge_weights);
std::vector<std::vector<double> > probability_matrix; template<class T>
fitter.fill_with_probabilities(probability_matrix); T& get(std::vector<T>& data, const Facet_const_iterator& facet) {
return data[ boost::get(facet_index_map, facet) ];
std::vector<int> 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<std::vector<double> > probability_matrix(number_of_centers, std::vector<double>(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<double>::epsilon());
// }
//}
//
// std::vector<int> 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<int>::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_handle, int>(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<double> 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();
} }
public:
/** /**
* Going to be removed * Going to be removed
*/ */
void write_sdf_values(const char* file_name) { void write_sdf_values(const char* file_name) {
std::ofstream output(file_name); std::ofstream output(file_name);
for(Facet_iterator facet_it = mesh->facets_begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
facet_it != mesh->facets_end(); ++facet_it) { facet_it != mesh.facets_end(); ++facet_it) {
output << sdf_values[facet_it] << std::endl; output << get(sdf_values, facet_it) << std::endl;
} }
output.close(); output.close();
} }
@ -764,31 +649,29 @@ public:
*/ */
void read_sdf_values(const char* file_name) { void read_sdf_values(const char* file_name) {
std::ifstream input(file_name); std::ifstream input(file_name);
sdf_values.clear(); sdf_values = std::vector<double>(mesh.size_of_facets());
for(Facet_iterator facet_it = mesh->facets_begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
facet_it != mesh->facets_end(); ++facet_it) { facet_it != mesh.facets_end(); ++facet_it) {
double sdf_value; double sdf_value;
input >> sdf_value; input >> sdf_value;
sdf_values.insert(std::pair<Facet_handle, double>(facet_it, sdf_value)); get(sdf_values, facet_it) = sdf_value;
} }
} }
/** /**
* Going to be removed * Going to be removed
*/ */
void read_center_ids(const char* file_name) { void read_center_ids(const char* file_name) {
std::ifstream input(file_name); //std::ifstream input(file_name);
centers.clear(); //centers.clear();
int max_center = 0; //int max_center = 0;
for(Facet_iterator facet_it = mesh->facets_begin(); //for(Facet_const_iterator facet_it = mesh.facets_begin(); facet_it != mesh.facets_end(); ++facet_it)
facet_it != mesh->facets_end(); ++facet_it) { //{
int center_id; // int center_id;
input >> center_id; // input >> center_id;
centers.insert(std::pair<Facet_handle, int>(facet_it, center_id)); // centers.insert(std::pair<Facet_const_iterator, int>(facet_it, center_id));
if(center_id > max_center) { // if(center_id > max_center) { max_center = center_id; }
max_center = center_id; //}
} //number_of_centers = max_center + 1;
}
number_of_centers = max_center + 1;
} }
/** /**
@ -812,9 +695,9 @@ public:
void write_segment_ids(const char* file_name) { void write_segment_ids(const char* file_name) {
assign_segments(); assign_segments();
std::ofstream output(file_name); std::ofstream output(file_name);
for(Facet_iterator facet_it = mesh->facets_begin(); for(Facet_const_iterator facet_it = mesh.facets_begin();
facet_it != mesh->facets_end(); ++facet_it) { facet_it != mesh.facets_end(); ++facet_it) {
output << segments[facet_it] << std::endl; output << get(segments, facet_it) << std::endl;
} }
output.close(); output.close();
} }
@ -837,7 +720,7 @@ public:
use_minimum_segment = true; use_minimum_segment = true;
multiplier_for_segment = 1.0 + i; multiplier_for_segment = 1.0 + i;
CGAL::Timer t; Timer t;
t.start(); t.start();
calculate_sdf_values(); calculate_sdf_values();
@ -877,7 +760,7 @@ public:
use_minimum_segment = false; use_minimum_segment = false;
multiplier_for_segment = 1.0 + i * 0.2; multiplier_for_segment = 1.0 + i * 0.2;
CGAL::Timer t; Timer t;
t.start(); t.start();
calculate_sdf_values(); calculate_sdf_values();
@ -920,6 +803,10 @@ public:
#undef CGAL_NORMALIZATION_ALPHA #undef CGAL_NORMALIZATION_ALPHA
#undef CGAL_CONVEX_FACTOR #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 #ifdef SEG_DEBUG
#undef SEG_DEBUG #undef SEG_DEBUG

142
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h Normal file
View File

@ -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<typename GeomTraits, typename Polyhedron>
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<GeomTraits, Polyhedron> 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 <class Iterator>
AABB_polyhedron_triangle_primitive( Iterator it,
typename boost::enable_if<
boost::is_same<Id,typename Iterator::value_type>
>::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<typename GeomTraits, typename Polyhedron>
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_

18
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut.h
View File

@ -1,5 +1,5 @@
#ifndef CGAL_ALPHA_EXPANSION_GRAPH_CUT_H #ifndef CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H
#define CGAL_ALPHA_EXPANSION_GRAPH_CUT_H #define CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H
#include <CGAL/assertions.h> #include <CGAL/assertions.h>
#include <CGAL/Timer.h> #include <CGAL/Timer.h>
@ -82,12 +82,12 @@ public:
int number_of_clusters = probability_matrix.size(); int number_of_clusters = probability_matrix.size();
double min_cut = (std::numeric_limits<double>::max)(); double min_cut = (std::numeric_limits<double>::max)();
bool success; bool success;
CGAL::Timer gt; Timer gt;
gt.start(); gt.start();
do { do {
success = false; success = false;
for(int alpha = 0; alpha < number_of_clusters; ++alpha) { for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
CGAL::Timer t; Timer t;
t.start(); t.start();
Graph graph; Graph graph;
#if 0 #if 0
@ -174,13 +174,13 @@ public:
int number_of_clusters = probability_matrix.size(); int number_of_clusters = probability_matrix.size();
double min_cut = (std::numeric_limits<double>::max)(); double min_cut = (std::numeric_limits<double>::max)();
bool success; bool success;
CGAL::Timer gt; Timer gt;
gt.start(); gt.start();
double total_time = 0.0; double total_time = 0.0;
do { do {
success = false; success = false;
for(int alpha = 0; alpha < number_of_clusters; ++alpha) { for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
CGAL::Timer t; Timer t;
t.start(); t.start();
Graph graph(edges.size() + labels.size() + Graph graph(edges.size() + labels.size() +
2); // allocate using maximum possible size. 2); // allocate using maximum possible size.
@ -272,7 +272,7 @@ public:
int number_of_clusters = probability_matrix.size(); int number_of_clusters = probability_matrix.size();
double min_cut = (std::numeric_limits<double>::max)(); double min_cut = (std::numeric_limits<double>::max)();
bool success; bool success;
CGAL::Timer gt; Timer gt;
gt.start(); gt.start();
int cluster_source = labels.size(); int cluster_source = labels.size();
int cluster_sink = labels.size() + 1; int cluster_sink = labels.size() + 1;
@ -285,7 +285,7 @@ public:
do { do {
success = false; success = false;
for(int alpha = 0; alpha < number_of_clusters; ++alpha) { for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
CGAL::Timer t; Timer t;
t.start(); t.start();
Graph graph(edges.begin(), edges.end(), labels.size()+2, edges.size()); Graph graph(edges.begin(), edges.end(), labels.size()+2, edges.size());
#if 0 #if 0
@ -305,4 +305,4 @@ public:
}; };
}//namespace internal }//namespace internal
}//namespace CGAL }//namespace CGAL
#endif //CGAL_ALPHA_EXPANSION_GRAPH_CUT_H #endif //CGAL_SURFACE_MESH_SEGMENTATION_ALPHA_EXPANSION_GRAPH_CUT_H

20
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Expectation_maximization.h
View File

@ -1,5 +1,5 @@
#ifndef CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H #ifndef CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H
#define CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H #define CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H
/* NEED TO BE DONE */ /* NEED TO BE DONE */
/* About implementation: /* About implementation:
* + There are a lot of parameters (with default values) and initialization types, * + 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, Gaussian_center new_center(initial_mean, initial_deviation,
initial_mixing_coefficient); initial_mixing_coefficient);
// if same point is choosen as a center twice, algorithm will not work // if same point is choosen as a center twice, algorithm will not work
if ( is_already_center(new_center) ) --i; if(is_already_center(new_center)) {
else centers.push_back(new_center); --i;
} else {
centers.push_back(new_center);
}
} }
calculate_initial_deviations(); calculate_initial_deviations();
} }
@ -289,8 +292,11 @@ protected:
Gaussian_center new_center(initial_mean, initial_deviation, Gaussian_center new_center(initial_mean, initial_deviation,
initial_mixing_coefficient); initial_mixing_coefficient);
// if same point is choosen as a center twice, algorithm will not work // if same point is choosen as a center twice, algorithm will not work
if ( is_already_center(new_center) ) --i; if(is_already_center(new_center)) {
else centers.push_back(new_center); --i;
} else {
centers.push_back(new_center);
}
} }
calculate_initial_deviations(); calculate_initial_deviations();
} }
@ -501,4 +507,4 @@ protected:
#ifdef SEG_DEBUG #ifdef SEG_DEBUG
#undef SEG_DEBUG #undef SEG_DEBUG
#endif #endif
#endif //CGAL_SEGMENTATION_EXPECTATION_MAXIMIZATION_H #endif //CGAL_SURFACE_MESH_SEGMENTATION_EXPECTATION_MAXIMIZATION_H

22
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/K_means_clustering.h
View File

@ -1,5 +1,5 @@
#ifndef CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H #ifndef CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H
#define CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H #define CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H
#include <vector> #include <vector>
#include <cmath> #include <cmath>
@ -79,7 +79,7 @@ protected:
/** /**
* Finds closest center and adds itself to the closest center's mean. * Finds closest center and adds itself to the closest center's mean.
* @param centers available centers * @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<K_means_center>& centers) { bool calculate_new_center(std::vector<K_means_center>& centers) {
int new_center_id = 0; int new_center_id = 0;
@ -156,8 +156,11 @@ protected:
for(int i = 0; i < number_of_centers; ++i) { for(int i = 0; i < number_of_centers; ++i) {
double initial_mean = points[rand() % points.size()].data; double initial_mean = points[rand() % points.size()].data;
K_means_center new_center(initial_mean); K_means_center new_center(initial_mean);
if ( is_already_center(new_center) ) --i; if(is_already_center(new_center)) {
else centers.push_back(new_center); --i;
} else {
centers.push_back(new_center);
}
} }
} }
@ -195,8 +198,11 @@ protected:
- distance_square_cumulative.begin(); - distance_square_cumulative.begin();
double initial_mean = points[selection_index].data; double initial_mean = points[selection_index].data;
K_means_center new_center(initial_mean); K_means_center new_center(initial_mean);
if ( is_already_center(new_center) ) --i; if(is_already_center(new_center)) {
else centers.push_back(new_center); --i;
} else {
centers.push_back(new_center);
}
} }
} }
@ -293,4 +299,4 @@ protected:
#undef CGAL_DEFAULT_MAXIMUM_ITERATION #undef CGAL_DEFAULT_MAXIMUM_ITERATION
#undef CGAL_DEFAULT_NUMBER_OF_RUN #undef CGAL_DEFAULT_NUMBER_OF_RUN
#endif //CGAL_SEGMENTATION_K_MEANS_CLUSTERING_H #endif //CGAL_SURFACE_MESH_SEGMENTATION_K_MEANS_CLUSTERING_H

102
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/SDF_calculation.h
View File

@ -1,43 +1,35 @@
#ifndef CGAL_SDF_CALCULATION_H #ifndef CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H
#define CGAL_SDF_CALCULATION_H #define CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H
#include <CGAL/AABB_tree.h> #include <CGAL/AABB_tree.h>
#include <CGAL/AABB_traits.h> #include <CGAL/AABB_traits.h>
#include <CGAL/AABB_polyhedron_triangle_primitive.h> #include <CGAL/AABB_polyhedron_triangle_primitive.h>
#include <CGAL/internal/Surface_mesh_segmentation/AABB_traversal_traits.h> #include <CGAL/internal/Surface_mesh_segmentation/AABB_traversal_traits.h>
#include <CGAL/internal/Surface_mesh_segmentation/AABB_const_polyhedron_triangle_primitive.h>
#include <map> #include <map>
#include <vector> #include <vector>
#include <algorithm> #include <algorithm>
#include <boost/tuple/tuple.hpp> #include <boost/tuple/tuple.hpp>
#include <boost/property_map/property_map.hpp>
#define CGAL_ANGLE_ST_DEV_DIVIDER 2.0 #define CGAL_ANGLE_ST_DEV_DIVIDER 2.0
#define CGAL_ACCEPTANCE_RATE_THRESHOLD 0.5 #define CGAL_ACCEPTANCE_RATE_THRESHOLD 0.5
#define CGAL_ST_DEV_MULTIPLIER 0.75 #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 CGAL
{ {
namespace internal namespace internal
{ {
template <class Polyhedron> template <
class Polyhedron
>
class SDF_calculation class SDF_calculation
{ {
public: 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 //type definitions
protected: protected:
typedef typename Polyhedron::Traits Kernel; typedef typename Polyhedron::Traits Kernel;
@ -45,17 +37,16 @@ protected:
typedef typename Polyhedron::Facet Vertex; typedef typename Polyhedron::Facet Vertex;
typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::Point_3 Point; typedef typename Kernel::Point_3 Point;
typedef typename Polyhedron::Facet_iterator Facet_iterator;
typedef typename Polyhedron::Facet_handle Facet_handle; typedef typename Polyhedron::Facet_const_iterator Facet_const_iterator;
typedef typename Polyhedron::Vertex_handle Vertex_handle;
typedef typename Kernel::Ray_3 Ray; typedef typename Kernel::Ray_3 Ray;
typedef typename Kernel::Plane_3 Plane; typedef typename Kernel::Plane_3 Plane;
typedef typename Kernel::Segment_3 Segment; typedef typename Kernel::Segment_3 Segment;
typedef typename CGAL::AABB_polyhedron_triangle_primitive<Kernel, Polyhedron> typedef typename AABB_const_polyhedron_triangle_primitive<Kernel, Polyhedron>
Primitive; Primitive;
typedef typename CGAL::AABB_tree<CGAL::AABB_traits<Kernel, Primitive> > typedef typename CGAL::AABB_tree<AABB_traits<Kernel, Primitive> >
Tree; Tree;
typedef typename Tree::Object_and_primitive_id typedef typename Tree::Object_and_primitive_id
Object_and_primitive_id; Object_and_primitive_id;
@ -68,7 +59,8 @@ protected:
//Member variables //Member variables
protected: protected:
Parameters parameters; double cone_angle;
int number_of_rays;
Disk_samples_list disk_samples_sparse; Disk_samples_list disk_samples_sparse;
Disk_samples_list disk_samples_dense; Disk_samples_list disk_samples_dense;
@ -77,31 +69,31 @@ protected:
double multiplier_for_segment; double multiplier_for_segment;
public: public:
SDF_calculation(Parameters parameters = Parameters()) SDF_calculation(double cone_angle, int number_of_rays)
: parameters(parameters), use_minimum_segment(false), : cone_angle(cone_angle), number_of_rays(number_of_rays),
multiplier_for_segment(1) { use_minimum_segment(false), multiplier_for_segment(1) {
} }
void calculate_sdf_values(Polyhedron& mesh, template < class FacetValueMap >
std::map<Facet_handle, double>& sdf_values) { void calculate_sdf_values(const Polyhedron& mesh, FacetValueMap sdf_values) {
int sparse_ray_count = parameters.number_of_rays; int sparse_ray_count = number_of_rays;
int dense_ray_count = sparse_ray_count * 2; int dense_ray_count = sparse_ray_count * 2;
disk_sampling_vogel_method(disk_samples_sparse, sparse_ray_count); disk_sampling_vogel_method(disk_samples_sparse, sparse_ray_count);
disk_sampling_vogel_method(disk_samples_dense, dense_ray_count); disk_sampling_vogel_method(disk_samples_dense, dense_ray_count);
Tree tree(mesh.facets_begin(), mesh.facets_end()); 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) { facet_it != mesh.facets_end(); ++facet_it) {
CGAL_precondition(facet_it->is_triangle()); //Mesh should contain triangles. CGAL_precondition(facet_it->is_triangle()); //Mesh should contain triangles.
double sdf = calculate_sdf_value_of_facet(facet_it, tree, disk_samples_sparse); double sdf = calculate_sdf_value_of_facet(facet_it, tree, disk_samples_sparse);
sdf_values.insert(std::pair<Facet_handle, double>(facet_it, sdf)); boost::put(sdf_values, facet_it, sdf);
} }
} }
protected: protected:
double calculate_sdf_value_of_facet(const Facet_handle& facet, const Tree& tree, double calculate_sdf_value_of_facet(Facet_const_iterator& facet,
const Disk_samples_list& samples) const { const Tree& tree, const Disk_samples_list& samples) const {
const Point& p1 = facet->halfedge()->vertex()->point(); const Point& p1 = facet->halfedge()->vertex()->point();
const Point& p2 = facet->halfedge()->next()->vertex()->point(); const Point& p2 = facet->halfedge()->next()->vertex()->point();
const Point& p3 = facet->halfedge()->prev()->vertex()->point(); const Point& p3 = facet->halfedge()->prev()->vertex()->point();
@ -116,10 +108,10 @@ protected:
//arrange_center_orientation(plane, normal, center); //arrange_center_orientation(plane, normal, center);
std::vector<double> ray_distances, ray_weights; std::vector<double> ray_distances, ray_weights;
ray_distances.reserve(parameters.number_of_rays); ray_distances.reserve(number_of_rays);
ray_weights.reserve(parameters.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; normal = normal * length_of_normal;
// stores segment length, // stores segment length,
// making it too large might cause a non-filtered bboxes in traversal, // 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 / std::sqrt(ray_direction.squared_length());
ray_direction = ray_direction * (*segment_distance * multiplier_for_segment); 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, boost::tie(is_intersected, intersection_is_acute,
min_distance) = cast_and_return_minimum(segment, tree, facet); min_distance) = cast_and_return_minimum(segment, tree, facet);
if(!is_intersected) { //no intersection is found if(!is_intersected) { //no intersection is found
@ -193,8 +185,8 @@ protected:
ray_distances.push_back(min_distance); ray_distances.push_back(min_distance);
} }
double sdf_result, acceptance_rate; double sdf_result, acceptance_rate;
boost::tie(sdf_result, acceptance_rate) = calculate_sdf_value_from_rays( boost::tie(sdf_result, acceptance_rate) =
ray_distances, ray_weights); remove_outliers_and_calculate_sdf_value(ray_distances, ray_weights);
if(acceptance_rate > CGAL_ACCEPTANCE_RATE_THRESHOLD if(acceptance_rate > CGAL_ACCEPTANCE_RATE_THRESHOLD
|| samples.size() == disk_samples_dense.size()) { || samples.size() == disk_samples_dense.size()) {
@ -205,7 +197,7 @@ protected:
template <class Query> //Query can be templated for just Ray and Segment types. template <class Query> //Query can be templated for just Ray and Segment types.
boost::tuple<bool, bool, double> cast_and_return_minimum( boost::tuple<bool, bool, double> 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<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<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. // get<2> : distance between ray/segment origin and intersection point.
@ -226,14 +218,14 @@ protected:
for(typename std::list<Object_and_primitive_id>::iterator op_it = for(typename std::list<Object_and_primitive_id>::iterator op_it =
intersections.begin(); intersections.begin();
op_it != intersections.end() ; ++op_it) { op_it != intersections.end() ; ++op_it) {
CGAL::Object object = op_it->first; Object object = op_it->first;
Primitive_id id = op_it->second; Primitive_id id = op_it->second;
if(id == facet) { if(id == facet) {
continue; //Since center is located on related facet, we should skip it if there is an intersection with it. continue; //Since center is located on related facet, we should skip it if there is an intersection with it.
} }
const Point* i_point; const Point* i_point;
if(!(i_point = CGAL::object_cast<Point>(&object))) { if(!(i_point = object_cast<Point>(&object))) {
continue; // Continue in case of segment. continue; // Continue in case of segment.
} }
@ -266,7 +258,7 @@ protected:
template <class Query> //Query can be templated for just Ray and Segment types. template <class Query> //Query can be templated for just Ray and Segment types.
boost::tuple<bool, bool, double> cast_and_return_minimum_use_closest ( boost::tuple<bool, bool, double> cast_and_return_minimum_use_closest (
const Query& ray, const Tree& tree, 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<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<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. // get<2> : distance between ray/segment origin and intersection point.
@ -285,7 +277,7 @@ protected:
} }
min_distance.get<0>() = true; // intersection is found min_distance.get<0>() = true; // intersection is found
CGAL::Object object = intersection->first; Object object = intersection->first;
Primitive_id min_id = intersection->second; Primitive_id min_id = intersection->second;
if(min_id == facet) { if(min_id == facet) {
@ -293,7 +285,7 @@ protected:
} }
const Point* i_point; const Point* i_point;
if(!(i_point = CGAL::object_cast<Point>(&object))) { if(!(i_point = object_cast<Point>(&object))) {
return min_distance; return min_distance;
} }
@ -310,7 +302,8 @@ protected:
min_distance.get<2>() = std::sqrt(min_i_ray.squared_length()); min_distance.get<2>() = std::sqrt(min_i_ray.squared_length());
return min_distance; return min_distance;
} }
boost::tuple<double, double> calculate_sdf_value_from_rays(
boost::tuple<double, double> remove_outliers_and_calculate_sdf_value(
std::vector<double>& ray_distances, std::vector<double>& ray_distances,
std::vector<double>& ray_weights) const { std::vector<double>& ray_weights) const {
// get<0> : sdf value // get<0> : sdf value
@ -357,8 +350,8 @@ protected:
for(std::vector<double>::iterator dist_it = ray_distances.begin(); for(std::vector<double>::iterator dist_it = ray_distances.begin();
dist_it != ray_distances.end(); dist_it != ray_distances.end();
++dist_it, ++w_it) { ++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; continue;
} }
total_distance += (*dist_it) * (*w_it); total_distance += (*dist_it) * (*w_it);
@ -385,17 +378,17 @@ protected:
//if(plane.has_on_positive_side(center)) { return; } //if(plane.has_on_positive_side(center)) { return; }
//Vector epsilon_normal = unit_normal * epsilon; //Vector epsilon_normal = unit_normal * epsilon;
//do { //do {
// center = CGAL::operator+(center, epsilon_normal); // center = operator+(center, epsilon_normal);
//} while(!plane.has_on_positive_side(center)); //} while(!plane.has_on_positive_side(center));
//Option-2 //Option-2
//if(plane.has_on_positive_side(center)) { return; } //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; //distance = distance > epsilon ? (distance + epsilon) : epsilon;
//Vector distance_normal = unit_normal * distance; //Vector distance_normal = unit_normal * distance;
// //
//do { //do {
// center = CGAL::operator+(center, distance_normal); // center = operator+(center, distance_normal);
//} while(!plane.has_on_positive_side(center)); //} while(!plane.has_on_positive_side(center));
//Option-3 //Option-3
@ -407,14 +400,14 @@ protected:
Vector epsilon_normal = unit_normal * epsilon; Vector epsilon_normal = unit_normal * epsilon;
do { do {
center = CGAL::operator+(center, epsilon_normal); center = operator+(center, epsilon_normal);
ray = Ray(center, unit_normal); ray = Ray(center, unit_normal);
} while(intersector(ray, plane)); } while(intersector(ray, plane));
} }
void disk_sampling_vogel_method(Disk_samples_list& samples, int ray_count) { 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 length_of_normal = 1.0 / tan(cone_angle / 2.0);
const double angle_st_dev = parameters.cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER; const double angle_st_dev = cone_angle / CGAL_ANGLE_ST_DEV_DIVIDER;
const double golden_ratio = 3.0 - std::sqrt(5.0); const double golden_ratio = 3.0 - std::sqrt(5.0);
#if 0 #if 0
@ -440,6 +433,5 @@ protected:
#undef CGAL_ANGLE_ST_DEV_DIVIDER #undef CGAL_ANGLE_ST_DEV_DIVIDER
#undef CGAL_ST_DEV_MULTIPLIER #undef CGAL_ST_DEV_MULTIPLIER
#undef CGAL_ACCEPTANCE_RATE_THRESHOLD #undef CGAL_ACCEPTANCE_RATE_THRESHOLD
#undef CGAL_DEFAULT_CONE_ANGLE
#undef CGAL_DEFAULT_NUMBER_OF_RAYS #endif //CGAL_SURFACE_MESH_SEGMENTATION_SDF_CALCULATION_H
#endif //CGAL_SDF_CALCULATION

38
Surface_mesh_segmentation/include/CGAL/mesh_segmentation.h Normal file
View File

@ -0,0 +1,38 @@
#ifndef CGAL_SURFACE_MESH_SEGMENTATION_MESH_SEGMENTATION_H
#define CGAL_SURFACE_MESH_SEGMENTATION_MESH_SEGMENTATION_H
#include <CGAL/Surface_mesh_segmentation.h>
#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 <class Polyhedron, class SDFPropertyMap>
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<Polyhedron> algorithm(polyhedron);
algorithm.calculate_sdf_values(sdf_values, cone_angle, number_of_rays);
}
template <class Polyhedron, class SDFPropertyMap, class SegmentPropertyMap>
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<Polyhedron> 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 //