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

First version of alpha expansion graph cut.

This commit is contained in:
Ílker Yaz 2012-06-19 18:17:19 +00:00
parent 9d28b10663
commit 8e0538d2c7
3 changed files with 190 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -126,6 +126,8 @@ public:
bool use_minimum_segment;
double multiplier_for_segment;
static long miss_counter;
internal::Expectation_maximization fitter;
//member functions
public:
Surface_mesh_segmentation(Polyhedron* mesh,
@ -196,8 +198,8 @@ inline Surface_mesh_segmentation<Polyhedron>::Surface_mesh_segmentation(
#endif
//
//write_sdf_values("sdf_values_sample_dino_ws.txt");
//read_sdf_values("sdf_values_sample_elephant.txt");
//apply_GMM_fitting();
//read_sdf_values("sdf_values_sample_camel.txt");
apply_GMM_fitting();
apply_graph_cut();
}
@ -292,15 +294,13 @@ Surface_mesh_segmentation<Polyhedron>::calculate_sdf_value_of_facet(
segment_distance = min_distance; //first assignment of the segment_distance
} else { // use segment_distance to limit rays as segments
ray_direction = ray_direction / sqrt(ray_direction.squared_length());
ray_direction = ray_direction * (*segment_distance * multiplier_for_segment);
Segment segment(center, CGAL::operator+(center, ray_direction));
boost::tie(is_intersected, intersection_is_acute,
min_distance) = cast_and_return_minimum(segment, tree, facet);
if(!is_intersected) {
//continue; // for utopia case - just continue on miss
if(!is_intersected) { //no intersection is found
++miss_counter;
//Ray ray(segment.target(), ray_direction);
Ray ray(segment.target(), ray_direction);
boost::tie(is_intersected, intersection_is_acute,
min_distance) = cast_and_return_minimum(ray, tree, facet);
@ -308,8 +308,9 @@ Surface_mesh_segmentation<Polyhedron>::calculate_sdf_value_of_facet(
continue;
}
min_distance += CGAL::sqrt(
segment.squared_length()); // since we our source is segment.target()
} else if(!intersection_is_acute) {
segment.squared_length()); // since our source is segment.target()
} else if(
!intersection_is_acute) { // intersection is found, but it is not acceptable (so, do not continue ray-segment casting)
continue;
}
@ -897,8 +898,8 @@ inline void Surface_mesh_segmentation<Polyhedron>::apply_GMM_fitting()
}
SEG_DEBUG(CGAL::Timer t)
SEG_DEBUG(t.start())
// apply em with 5 runs, number of runs might become a parameter.
internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10);
//internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 10);
fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 10);
SEG_DEBUG(std::cout << t.time() << std::endl)
std::vector<int> center_memberships;
fitter.fill_with_center_ids(center_memberships);
@ -927,7 +928,6 @@ inline void Surface_mesh_segmentation<Polyhedron>::apply_K_means_clustering()
pair_it != sdf_values.end(); ++pair_it, ++center_it) {
centers.insert(std::pair<Facet_handle, int>(pair_it->first, (*center_it)));
}
//center_memberships_temp = center_memberships; //remove
}
template <class Polyhedron>
inline void
@ -944,8 +944,9 @@ Surface_mesh_segmentation<Polyhedron>::apply_GMM_fitting_with_K_means_init()
std::vector<int> center_memberships;
clusterer.fill_with_center_ids(center_memberships);
//std::vector<int> center_memberships = center_memberships_temp;
internal::Expectation_maximization fitter(number_of_centers, sdf_vector,
center_memberships);
//internal::Expectation_maximization fitter(number_of_centers, sdf_vector, center_memberships);
fitter = internal::Expectation_maximization(number_of_centers, sdf_vector,
center_memberships);
center_memberships.clear();
fitter.fill_with_center_ids(center_memberships);
std::vector<int>::iterator center_it = center_memberships.begin();
@ -976,7 +977,7 @@ void Surface_mesh_segmentation<Polyhedron>::apply_graph_cut()
int index_f2 = facet_indices[edge_it->opposite()->facet()];
edges.push_back(std::pair<int, int>(index_f1, index_f2));
angle = -log(angle);
angle *= 15; //lambda, will be variable.
angle *= 10; //lambda, will be variable.
// we may also want to consider edge lengths, also penalize convex angles.
edge_weights.push_back(angle);
}
@ -987,7 +988,7 @@ void Surface_mesh_segmentation<Polyhedron>::apply_graph_cut()
pair_it != sdf_values.end(); ++pair_it) {
sdf_vector.push_back(pair_it->second);
}
internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 50);
//internal::Expectation_maximization fitter(number_of_centers, sdf_vector, 3);
//fill probability matrix.
std::vector<std::vector<double> > probability_matrix;
fitter.fill_with_minus_log_probabilities(probability_matrix);

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

@ -43,10 +43,173 @@ public:
typedef Traits::vertex_iterator Vertex_iterator;
typedef Traits::edge_iterator Edge_iterator;
Alpha_expansion_graph_cut(std::vector<std::pair<int, int> >& edges,
std::vector<double>& edge_weights, std::vector<int>& labels,
std::vector<std::vector<double> >& probability_matrix,
Alpha_expansion_graph_cut(const std::vector<std::pair<int, int> >& edges,
const std::vector<double>& edge_weights, std::vector<int>& labels,
const std::vector<std::vector<double> >& probability_matrix,
std::vector<int>& center_ids) {
apply_alpha_expansion_2(edges, edge_weights, probability_matrix, labels);
center_ids = labels;
}
boost::tuple<Edge_descriptor, Edge_descriptor>
add_edge_and_reverse(Vertex_descriptor& v1, Vertex_descriptor& v2, double w1,
double w2, Graph& graph) {
Edge_descriptor v1_v2, v2_v1;
bool v1_v2_added, v2_v1_added;
tie(v1_v2, v1_v2_added) = boost::add_edge(v1, v2, graph);
tie(v2_v1, v2_v1_added) = boost::add_edge(v2, v1, graph);
CGAL_assertion(v1_v2_added && v2_v1_added);
//put edge capacities
boost::put(boost::edge_reverse, graph, v1_v2, v2_v1);
boost::put(boost::edge_reverse, graph, v2_v1, v1_v2);
//map reverse edges
boost::put(boost::edge_capacity, graph, v1_v2, w1);
boost::put(boost::edge_capacity, graph, v2_v1, w2);
return boost::make_tuple(v1_v2, v2_v1);
}
void apply_alpha_expansion(const std::vector<std::pair<int, int> >& edges,
const std::vector<double>& edge_weights,
const std::vector<std::vector<double> >& probability_matrix,
std::vector<int>& labels) {
int number_of_clusters = probability_matrix.size();
double min_cut = (std::numeric_limits<double>::max)();
bool success;
do {
success = false;
for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
Graph graph;
Vertex_descriptor cluster_source = boost::add_vertex(graph);
Vertex_descriptor cluster_sink = boost::add_vertex(graph);
std::vector<Vertex_descriptor> inserted_vertices;
// For E-Data
// add every facet as a vertex to graph, put edges to source-sink vertices
for(std::size_t vertex_i = 0; vertex_i < probability_matrix[0].size();
++vertex_i) {
Vertex_descriptor new_vertex = boost::add_vertex(graph);
inserted_vertices.push_back(new_vertex);
double source_weight = probability_matrix[alpha][vertex_i];
double sink_weight = (labels[vertex_i] == alpha) ?
(std::numeric_limits<double>::max)() :
probability_matrix[labels[vertex_i]][vertex_i];
add_edge_and_reverse(cluster_source, new_vertex, source_weight, 0.0, graph);
add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, 0.0, graph);
}
// For E-Smooth
// add edge between every vertex,
std::vector<double>::const_iterator weight_it = edge_weights.begin();
for(std::vector<std::pair<int, int> >::const_iterator edge_it = edges.begin();
edge_it != edges.end();
++edge_it, ++weight_it) {
Vertex_descriptor v1 = inserted_vertices[edge_it->first];
Vertex_descriptor v2 = inserted_vertices[edge_it->second];
add_edge_and_reverse(v1, v2, *weight_it, *weight_it, graph);
}
double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source,
cluster_sink);
if(min_cut - flow < flow * 1e-3) {
continue;
}
std::cout << "prev flow: " << min_cut << "new flow: " << flow << std::endl;
min_cut = flow;
success = true;
//update labeling
for(std::size_t vertex_i = 0; vertex_i < inserted_vertices.size(); ++vertex_i) {
boost::default_color_type color = boost::get(boost::vertex_color, graph,
inserted_vertices[vertex_i]);
if(labels[vertex_i] != alpha && color == ColorTraits::white()) { //new comers
labels[vertex_i] = alpha;
}
}
}
} while(success);
}
void apply_alpha_expansion_2(const std::vector<std::pair<int, int> >& edges,
const std::vector<double>& edge_weights,
const std::vector<std::vector<double> >& probability_matrix,
std::vector<int>& labels) {
int number_of_clusters = probability_matrix.size();
double min_cut = (std::numeric_limits<double>::max)();
bool success;
do {
success = false;
for(int alpha = 0; alpha < number_of_clusters; ++alpha) {
Graph graph;
Vertex_descriptor cluster_source = boost::add_vertex(graph);
Vertex_descriptor cluster_sink = boost::add_vertex(graph);
std::vector<Vertex_descriptor> inserted_vertices;
// For E-Data
// add every facet as a vertex to graph, put edges to source-sink vertices
for(std::size_t vertex_i = 0; vertex_i < probability_matrix[0].size();
++vertex_i) {
Vertex_descriptor new_vertex = boost::add_vertex(graph);
inserted_vertices.push_back(new_vertex);
double source_weight = probability_matrix[alpha][vertex_i];
double sink_weight = (labels[vertex_i] == alpha) ?
(std::numeric_limits<double>::max)() :
probability_matrix[labels[vertex_i]][vertex_i];
add_edge_and_reverse(cluster_source, new_vertex, source_weight, 0.0, graph);
add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, 0.0, graph);
}
// For E-Smooth
// add edge between every vertex,
std::vector<double>::const_iterator weight_it = edge_weights.begin();
for(std::vector<std::pair<int, int> >::const_iterator edge_it = edges.begin();
edge_it != edges.end();
++edge_it, ++weight_it) {
Vertex_descriptor v1 = inserted_vertices[edge_it->first],
v2 = inserted_vertices[edge_it->second];
int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
if(label_1 == label_2) {
double w1 = label_1 == alpha ? 0 : *weight_it;
add_edge_and_reverse(v1, v2, w1, w1, graph);
} else {
Vertex_descriptor inbetween = boost::add_vertex(graph);
add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
double w1 = label_1 == alpha ? 0 : *weight_it;
double w2 = label_2 == alpha ? 0 : *weight_it;
add_edge_and_reverse(inbetween, v1, w1, w1, graph);
add_edge_and_reverse(inbetween, v2, w2, w2, graph);
}
}
double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source,
cluster_sink);
if(min_cut - flow < flow * 1e-3) {
continue;
}
std::cout << "prev flow: " << min_cut << "new flow: " << flow << std::endl;
min_cut = flow;
success = true;
//update labeling
for(std::size_t vertex_i = 0; vertex_i < inserted_vertices.size(); ++vertex_i) {
boost::default_color_type color = boost::get(boost::vertex_color, graph,
inserted_vertices[vertex_i]);
if(labels[vertex_i] != alpha && color == ColorTraits::white()) { //new comers
labels[vertex_i] = alpha;
}
}
}
} while(success);
}
void apply_simple_cut(const std::vector<std::pair<int, int> >& edges,
const std::vector<double>& edge_weights, std::vector<int>& labels,
const std::vector<std::vector<double> >& probability_matrix,
std::vector<int>& center_ids) {
//Graph graph(edges.begin(), edges.end(), vertex_weights.size(), edge_weights.size());
Graph graph;
Vertex_descriptor cluster_source = boost::add_vertex(graph);
@ -86,8 +249,8 @@ public:
boost::put(boost::edge_reverse, graph, sink_e, sink_e_rev);
boost::put(boost::edge_reverse, graph, sink_e_rev, sink_e);
}
std::vector<double>::iterator weight_it = edge_weights.begin();
for(std::vector<std::pair<int, int> >::iterator edge_it = edges.begin();
std::vector<double>::const_iterator weight_it = edge_weights.begin();
for(std::vector<std::pair<int, int> >::const_iterator edge_it = edges.begin();
edge_it != edges.end();
++edge_it, ++weight_it) {
Vertex_descriptor v1 = inserted_vertices[edge_it->first];
@ -101,18 +264,14 @@ public:
CGAL_assertion(edge_added && edge_rev_added);
//put edge capacities
if(labels[edge_it->first] == labels[edge_it->second]) {
boost::put(boost::edge_capacity, graph, edge, *weight_it);
boost::put(boost::edge_capacity, graph, edge_rev, *weight_it);
} else {
boost::put(boost::edge_capacity, graph, edge, *weight_it);
boost::put(boost::edge_capacity, graph, edge_rev, *weight_it);
}
boost::put(boost::edge_capacity, graph, edge, *weight_it);
boost::put(boost::edge_capacity, graph, edge_rev, *weight_it);
//map reverse edges
boost::put(boost::edge_reverse, graph, edge, edge_rev);
boost::put(boost::edge_reverse, graph, edge_rev, edge);
}
double flow = boost::boykov_kolmogorov_max_flow(graph, cluster_source,
cluster_sink);
@ -121,19 +280,12 @@ public:
vertex_it != inserted_vertices.end(); ++vertex_it) {
boost::default_color_type color = boost::get(boost::vertex_color, graph,
*vertex_it);
if(color == ColorTraits::black() ) { // in sink
if(color == ColorTraits::black()) { // in sink
center_ids.push_back(1);
} else {
center_ids.push_back(0);
}
}
std::cout << flow << std::endl;
//for(std::pair<Vertex_iterator, Vertex_iterator> pair_v = boost::vertices(graph);
// pair_v.first != pair_v.second; ++(pair_v.first))
//{
// if(cluster_source == (*pair_v.first) || cluster_sink == (*pair_v.first)) { continue; }
//
//}
}
};
}//namespace internal

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

@ -70,6 +70,7 @@ protected:
std::vector<std::vector<double> > membership_matrix;
public:
Expectation_maximization() { }
/* For uniform initialization, with one run */
Expectation_maximization(int number_of_centers, const std::vector<double>& data,
double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER)
@ -136,7 +137,7 @@ public:
sum += probability;
probabilities[center_i][point_i] = probability;
}
#if 0
#if 1
// pdf values scaled so that their sum will equal to 1.
for(std::size_t point_i = 0; point_i < points.size(); ++point_i) {
double probability = probabilities[center_i][point_i] / sum;