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

Demo update 

K-means: implemented k-means++ for initialization (just for try)
This commit is contained in:
Ílker Yaz 2012-06-28 17:31:50 +00:00
parent e7e45a3587
commit 23465d67a2
7 changed files with 358 additions and 321 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
Polyhedron/demo/Polyhedron/Mesh_segmentation_dialog.ui
View File

@ -5,8 +5,8 @@
<rect>
<x>0</x>
<y>0</y>
<width>328</width>
<height>149</height>
<width>352</width>
<height>151</height>
</rect>
</property>
<property name="windowTitle" >
@ -15,30 +15,6 @@
<layout class="QVBoxLayout" >
<item>
<layout class="QGridLayout" >
<item row="0" column="1" >
<widget class="QSpinBox" name="Number_of_rays_spin_box" >
<property name="value" >
<number>7</number>
</property>
</widget>
</item>
<item row="2" column="1" >
<widget class="QSpinBox" name="Number_of_clusters_spin_box" >
<property name="value" >
<number>2</number>
</property>
</widget>
</item>
<item row="2" column="0" >
<widget class="QLabel" name="Number_of_clusters_label" >
<property name="text" >
<string>Number of clusters:</string>
</property>
<property name="alignment" >
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
</widget>
</item>
<item row="0" column="0" >
<widget class="QLabel" name="Number_of_rays_label" >
<property name="text" >
@ -49,6 +25,13 @@
</property>
</widget>
</item>
<item row="0" column="1" >
<widget class="QSpinBox" name="Number_of_rays_spin_box" >
<property name="value" >
<number>5</number>
</property>
</widget>
</item>
<item row="1" column="0" >
<widget class="QLabel" name="Cone_angle_label" >
<property name="text" >
@ -70,19 +53,53 @@
</widget>
</item>
<item row="3" column="0" >
<widget class="QCheckBox" name="Show_sdf_values_check_box" >
<widget class="QLabel" name="Number_of_clusters_label" >
<property name="text" >
<string>Show SDF values</string>
<string>Number of clusters:</string>
</property>
<property name="checked" >
<bool>true</bool>
<property name="alignment" >
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
</widget>
</item>
<item row="3" column="1" >
<widget class="QCheckBox" name="Show_clusters_check_box" >
<widget class="QSpinBox" name="Number_of_clusters_spin_box" >
<property name="value" >
<number>5</number>
</property>
</widget>
</item>
<item row="4" column="0" >
<widget class="QLabel" name="Smoothness_label" >
<property name="text" >
<string>Show clusters</string>
<string>Smoothness:</string>
</property>
<property name="alignment" >
<set>Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter</set>
</property>
</widget>
</item>
<item row="4" column="1" >
<widget class="QDoubleSpinBox" name="Smoothness_spin_box" >
<property name="singleStep" >
<double>0.500000000000000</double>
</property>
<property name="value" >
<double>15.500000000000000</double>
</property>
</widget>
</item>
<item row="2" column="0" >
<widget class="Line" name="line" >
<property name="orientation" >
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="2" column="1" >
<widget class="Line" name="line_2" >
<property name="orientation" >
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
@ -104,9 +121,16 @@
</spacer>
</item>
<item>
<widget class="QPushButton" name="Apply_button" >
<widget class="QPushButton" name="SDF_button" >
<property name="text" >
<string>Apply</string>
<string>Caluclate SDF</string>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="Partition_button" >
<property name="text" >
<string>Partition</string>
</property>
</widget>
</item>

73
Polyhedron/demo/Polyhedron/Polyhedron_demo_mesh_segmentation_plugin.cpp
View File

@ -37,10 +37,11 @@ public:
}
}
void colorize(CGAL::Surface_mesh_segmentation<Polyhedron>& segmentation, std::vector<QColor>& color_vector, bool sdf, bool cluster);
void colorize(CGAL::Surface_mesh_segmentation<Polyhedron>& segmentation, std::vector<QColor>& color_vector, bool sdf);
public slots:
void on_actionSegmentation_triggered();
void on_Apply_button_clicked();
void on_Partition_button_clicked();
void on_SDF_button_clicked();
private:
QAction* actionSegmentation;
@ -63,17 +64,16 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_actionSegmentation_triggered()
ui_dialog = new Ui::Mesh_segmentation_dialog();
ui_dialog->setupUi(dialog);
connect(ui_dialog->Apply_button, SIGNAL(clicked()), this, SLOT(on_Apply_button_clicked()));
connect(ui_dialog->Partition_button, SIGNAL(clicked()), this, SLOT(on_Partition_button_clicked()));
connect(ui_dialog->SDF_button, SIGNAL(clicked()), this, SLOT(on_SDF_button_clicked()));
dialog->show();
}
//void calculata_sdf_values(CGAL::Surface_mesh_segmentation<Polyhedron>* segmentation)
//{
// segmentation->calculate_sdf_values();
//}
void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked()
void Polyhedron_demo_mesh_segmentation_plugin::on_SDF_button_clicked()
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_polyhedron_item* item = qobject_cast<Scene_polyhedron_item*>(scene->item(index));
Scene_polyhedron_with_color_item* item_colored = qobject_cast<Scene_polyhedron_with_color_item*>(scene->item(index));
@ -84,8 +84,6 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked()
int number_of_rays_sqrt = ui_dialog->Number_of_rays_spin_box->value();
double cone_angle = (ui_dialog->Cone_angle_spin_box->value() / 180.0) * CGAL_PI;
int number_of_clusters = ui_dialog->Number_of_clusters_spin_box->value();
bool show_sdf_values = ui_dialog->Show_sdf_values_check_box->isChecked();
bool show_clusters = ui_dialog->Show_clusters_check_box->isChecked();
if(!item_colored)
{
@ -99,7 +97,7 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked()
CGAL::Surface_mesh_segmentation<Polyhedron>* segmentation
= new CGAL::Surface_mesh_segmentation<Polyhedron>(new_item->polyhedron(), number_of_rays_sqrt, cone_angle, number_of_clusters);
colorize(*segmentation, color_vector, show_sdf_values, show_clusters);
colorize(*segmentation, color_vector, true);
new_item->set_color_vector(color_vector);
new_item->segmentation = segmentation;
@ -122,40 +120,73 @@ void Polyhedron_demo_mesh_segmentation_plugin::on_Apply_button_clicked()
// ui_dialog->Sdf_value_calculation_bar->setValue(static_cast<int>(segmentation->get_process_of_sdf_calculation() * 100));
segmentation->calculate_sdf_values();
//segmentation->apply_GMM_fitting_with_K_means_init();
}
segmentation->number_of_centers = number_of_clusters;
segmentation->apply_GMM_fitting_with_K_means_init();
}
colorize(*segmentation, color_vector, show_sdf_values, show_clusters);
colorize(*segmentation, color_vector, true);
item_colored->set_color_vector(color_vector);
scene->itemChanged(index);
scene->setSelectedItem(index);
}
QApplication::restoreOverrideCursor();
}
void Polyhedron_demo_mesh_segmentation_plugin::on_Partition_button_clicked()
{
const Scene_interface::Item_id index = scene->mainSelectionIndex();
Scene_polyhedron_with_color_item* item_colored = qobject_cast<Scene_polyhedron_with_color_item*>(scene->item(index));
if(!item_colored) { return; }
QApplication::setOverrideCursor(Qt::WaitCursor);
int number_of_rays_sqrt = ui_dialog->Number_of_rays_spin_box->value();
double cone_angle = (ui_dialog->Cone_angle_spin_box->value() / 180.0) * CGAL_PI;
int number_of_clusters = ui_dialog->Number_of_clusters_spin_box->value();
double smoothness = ui_dialog->Smoothness_spin_box->value();
std::vector<QColor> color_vector;
CGAL::Surface_mesh_segmentation<Polyhedron>* segmentation = item_colored->segmentation;
segmentation->number_of_centers = number_of_clusters;
segmentation->apply_GMM_fitting_and_K_means();
segmentation->smoothing_lambda = smoothness;
segmentation->apply_graph_cut();
segmentation->assign_segments();
colorize(*segmentation, color_vector, false);
item_colored->set_color_vector(color_vector);
scene->itemChanged(index);
scene->setSelectedItem(index);
QApplication::restoreOverrideCursor();
//qApp->processEvents();
}
void Polyhedron_demo_mesh_segmentation_plugin::colorize(CGAL::Surface_mesh_segmentation<Polyhedron>& segmentation,
std::vector<QColor>& color_vector, bool sdf, bool cluster)
std::vector<QColor>& color_vector, bool sdf)
{
QColor predefined_list[] = { Qt::red, Qt::darkRed, Qt::green, Qt::darkGreen, Qt::blue, Qt::darkBlue
, Qt::cyan, Qt::darkCyan, Qt::magenta, Qt::darkMagenta, Qt::yellow, Qt::darkYellow, Qt::gray, Qt::darkGray, Qt::lightGray, Qt::white, Qt::black};
color_vector.reserve(segmentation.sdf_values.size());
int findex = 0;
for(CGAL::Surface_mesh_segmentation<Polyhedron>::Facet_iterator facet_it = segmentation.mesh->facets_begin();
facet_it != segmentation.mesh->facets_end(); ++facet_it, ++findex)
{
facet_it->set_patch_id(findex);
int color = 0;
int sdf_color = (int) (255 * segmentation.sdf_values[facet_it]);
int cluster_color = (int) (255.0 / segmentation.number_of_centers) * segmentation.centers[facet_it];
if(sdf && cluster)
//int cluster_color = (int) (255.0 / segmentation.number_of_centers) * segmentation.centers[facet_it];
if(sdf)
{
color = (sdf_color + cluster_color) / 2;
int sdf_color = (int) (255 * segmentation.sdf_values[facet_it]);
color_vector.push_back(QColor(sdf_color,sdf_color,sdf_color));
}
else
{
color = sdf ? sdf_color : cluster_color;
color_vector.push_back(predefined_list[segmentation.segments[facet_it] % 17]);
}
color_vector.push_back(QColor(color,color,color));
}
}

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

@ -73,8 +73,8 @@ public:
typedef typename Polyhedron::Edge_iterator Edge_iterator;
typedef typename Polyhedron::Vertex_handle Vertex_handle;
protected:
typedef typename Kernel::Ray_3 Ray;
typedef typename Kernel::Plane_3 Plane;
typedef typename Kernel::Ray_3 Ray;
typedef typename Kernel::Plane_3 Plane;
typedef typename Kernel::Segment_3 Segment;
typedef typename CGAL::AABB_polyhedron_triangle_primitive<Kernel, Polyhedron>
@ -90,9 +90,10 @@ protected:
typedef std::map<Facet_handle, int> Face_center_map;
typedef std::map<Halfedge_handle, double> Edge_angle_map;
typedef std::map<Facet_handle, int> Face_segment_map;
/*Sampled points from disk, t1 = coordinate-x, t2 = coordinate-y, t3 = angle with cone-normal (weight). */
typedef CGAL::Triple<double, double, double> Disk_sample;
typedef std::vector<CGAL::Triple<double, double, double> > Disk_samples_list;
typedef CGAL::Triple<double, double, double> Disk_sample;
typedef std::vector<Disk_sample> Disk_samples_list;
template <typename ValueTypeName>
struct Compare_second_element {
@ -151,12 +152,12 @@ public:
boost::tuple<bool, bool, double> cast_and_return_minimum_use_closest(
const Query& ray, const Tree& tree, const Facet_handle& facet) const;
double calculate_sdf_value_from_rays (std::vector<double>& ray_distances,
std::vector<double>& ray_weights) const;
double calculate_sdf_value_from_rays_with_mean (std::vector<double>&
double calculate_sdf_value_from_rays(std::vector<double>& ray_distances,
std::vector<double>& ray_weights) const;
double calculate_sdf_value_from_rays_with_mean(std::vector<double>&
ray_distances, std::vector<double>& ray_weights) const;
double calculate_sdf_value_from_rays_with_trimmed_mean(std::vector<double>&
ray_distances, std::vector<double>& ray_weights) const;
double calculate_sdf_value_from_rays_with_trimmed_mean (
std::vector<double>& ray_distances, std::vector<double>& ray_weights) const;
void arrange_center_orientation(const Plane& plane, const Vector& unit_normal,
Point& center) const;
@ -174,12 +175,13 @@ public:
void apply_GMM_fitting();
void apply_K_means_clustering();
void apply_GMM_fitting_with_K_means_init();
void apply_GMM_fitting_and_K_means();
void apply_graph_cut();
void apply_graph_cut_with_EM();
void apply_graph_cut_multiple_run(int number_of_run = 5);
void assign_segments();
void dfs(Facet_handle facet, int segment_id);
void write_sdf_values(const char* file_name);
@ -212,10 +214,10 @@ inline Surface_mesh_segmentation<Polyhedron>::Surface_mesh_segmentation(
//write_sdf_values("sdf_values_sample_teddy.txt");
//write_sdf_values("18_3.txt");
//read_sdf_values("sdf_values_sample_camel.txt");
apply_GMM_fitting_with_K_means_init();
apply_graph_cut();
//read_sdf_values("C:/Users/Ilker/Documents/Visual Studio 2008/Projects/SurfaceMeshSegmentation/SurfaceMeshSegmentation/models/sdf_values_sample_elephant.txt");
//read_sdf_values("D:/GSoC/MeshsegBenchmark-1.0-full/MeshsegBenchmark-1.0/data/off/173.txt");
//apply_GMM_fitting_with_K_means_init();
//apply_graph_cut();
}
template <class Polyhedron>
@ -230,13 +232,6 @@ inline void Surface_mesh_segmentation<Polyhedron>::calculate_sdf_values()
double sdf = calculate_sdf_value_of_facet(facet_it, tree);
sdf_values.insert(std::pair<Facet_handle, double>(facet_it, sdf));
}
//std::cout << Listing_intersection_traits_ray_or_segment_triangle
// < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list<Object_and_primitive_id> > >::inter_counter << std::endl;
//std::cout << Listing_intersection_traits_ray_or_segment_triangle
// < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list<Object_and_primitive_id> > >::true_inter_counter << std::endl;
//std::cout << Listing_intersection_traits_ray_or_segment_triangle
// < typename Tree::AABB_traits, Ray, std::back_insert_iterator< std::list<Object_and_primitive_id> > >::do_inter_counter << std::endl;
normalize_sdf_values();
smooth_sdf_values();
}
@ -348,16 +343,17 @@ Surface_mesh_segmentation<Polyhedron>::calculate_sdf_value_of_facet(
ray_weights);
}
// just for Ray and Segment
template <class Polyhedron>
template <class Query>
template <class Query> //Query can be templated for just Ray and Segment types.
boost::tuple<bool, bool, double>
Surface_mesh_segmentation<Polyhedron>::cast_and_return_minimum(
const Query& query, const Tree& tree, const Facet_handle& facet) const
{
boost::tuple<bool, bool, double> min_distance = boost::make_tuple(false, false,
0);
// get<0> : if any intersection is found then true
// get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true
// get<2> : distance between ray/segment origin and intersection point.
boost::tuple<bool, bool, double> min_distance(false, false, 0);
std::list<Object_and_primitive_id> intersections;
#if 1
//SL: the difference with all_intersections is that in the traversal traits, we do do_intersect before calling intersection.
@ -381,13 +377,9 @@ Surface_mesh_segmentation<Polyhedron>::cast_and_return_minimum(
}
const Point* i_point;
//Point i_point;
//AF: Use object_cast as it is faster than assign
//IOY: ok, also using pointer-returning version to get rid of copying
if(!(i_point = CGAL::object_cast<Point>(&object))) {
continue;
continue; // Continue in case of segment.
}
//if(!CGAL::assign(i_point, object)) { continue; } //What to do here (in case of intersection object is a segment), I am not sure ???
Vector i_ray = (query.source() - *i_point);
double new_distance = i_ray.squared_length();
@ -412,7 +404,6 @@ Surface_mesh_segmentation<Polyhedron>::cast_and_return_minimum(
return min_distance;
}
min_distance.get<1>() = true; // founded intersection is acceptable.
CGAL_assertion(min_distance.get<2>() > 0);
min_distance.get<2>() = sqrt(min_distance.get<2>());
return min_distance;
}
@ -424,12 +415,10 @@ Surface_mesh_segmentation<Polyhedron>::cast_and_return_minimum_use_closest (
const Query& ray, const Tree& tree,
const Facet_handle& facet) const
{
//static double dist = 0.1;
//boost::optional<double> min_distance_2 = dist++;
//return min_distance_2;
boost::tuple<bool, bool, double> min_distance = boost::make_tuple(false, false,
0);
// get<0> : if any intersection is found then true
// get<1> : if found intersection is acceptable (i.e. accute angle with surface normal) then true
// get<2> : distance between ray/segment origin and intersection point.
boost::tuple<bool, bool, double> min_distance(false, false, 0);
#if 1
Closest_intersection_traits<typename Tree::AABB_traits, Query> traversal_traits;
tree.traversal(ray, traversal_traits);
@ -733,7 +722,7 @@ double Surface_mesh_segmentation<Polyhedron>::calculate_dihedral_angle_of_edge_2
double n_angle = CGAL::Mesh_3::dihedral_angle(a, b, c, d) / 180.0;
bool n_concave = n_angle > 0;
double folded_angle = 1 + ((n_concave ? -1 : +1) * n_angle);
folded_angle = (CGAL::max)(folded_angle, epsilon);
folded_angle = (std::max)(folded_angle, epsilon);
if(!n_concave) {
return epsilon; // we may want to also penalize convex angles as well...
@ -905,28 +894,65 @@ inline void Surface_mesh_segmentation<Polyhedron>::smooth_sdf_values()
}
sdf_values = smoothed_sdf_values;
}
template <class Polyhedron>
inline void Surface_mesh_segmentation<Polyhedron>::apply_GMM_fitting()
{
centers.clear();
std::vector<double> sdf_vector;
sdf_vector.reserve(sdf_values.size());
for(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it) {
sdf_vector.push_back(pair_it->second);
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, 15);
fitter = internal::Expectation_maximization(number_of_centers, sdf_vector, 20);
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(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it, ++center_it) {
centers.insert(std::pair<Facet_handle, int>(pair_it->first, (*center_it)));
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)));
}
}
template <class Polyhedron>
inline void
Surface_mesh_segmentation<Polyhedron>::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, 20);
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, center_memberships);
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)));
}
}
@ -936,17 +962,19 @@ inline void Surface_mesh_segmentation<Polyhedron>::apply_K_means_clustering()
centers.clear();
std::vector<double> sdf_vector;
sdf_vector.reserve(sdf_values.size());
for(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it) {
sdf_vector.push_back(pair_it->second);
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(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it, ++center_it) {
centers.insert(std::pair<Facet_handle, int>(pair_it->first, (*center_it)));
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
}
@ -957,9 +985,10 @@ Surface_mesh_segmentation<Polyhedron>::apply_GMM_fitting_with_K_means_init()
centers.clear();
std::vector<double> sdf_vector;
sdf_vector.reserve(sdf_values.size());
for(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it) {
sdf_vector.push_back(pair_it->second);
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;
@ -971,9 +1000,10 @@ Surface_mesh_segmentation<Polyhedron>::apply_GMM_fitting_with_K_means_init()
center_memberships.clear();
fitter.fill_with_center_ids(center_memberships);
std::vector<int>::iterator center_it = center_memberships.begin();
for(typename Face_value_map::iterator pair_it = sdf_values.begin();
pair_it != sdf_values.end(); ++pair_it, ++center_it) {
centers.insert(std::pair<Facet_handle, int>(pair_it->first, (*center_it)));
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)));
}
}
@ -998,11 +1028,12 @@ 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 = (CGAL::max)(angle, 1e-5);
angle = (std::max)(angle, 1e-5);
angle *= smoothing_lambda; //lambda, will be variable.
// we may also want to consider edge lengths, also penalize convex angles.
edge_weights.push_back(angle);
}
//apply gmm fitting
std::vector<double> sdf_vector;
sdf_vector.reserve(sdf_values.size());

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

@ -23,7 +23,6 @@ public:
typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS,
// 4 vertex properties (nested)
//boost::no_property
boost::property<boost::vertex_index_t, int,
boost::property<boost::vertex_color_t, boost::default_color_type,
boost::property<boost::vertex_distance_t, double,
@ -47,8 +46,8 @@ public:
const std::vector<double>& edge_weights,
const std::vector<std::vector<double> >& probability_matrix,
std::vector<int>& labels, double* result = NULL) {
double min_cut = apply_alpha_expansion_2(edges, edge_weights,
probability_matrix, labels);
double min_cut = apply_alpha_expansion(edges, edge_weights, probability_matrix,
labels);
if(result != NULL) {
*result = min_cut;
}
@ -75,10 +74,10 @@ public:
return boost::make_tuple(v1_v2, v2_v1);
}
double 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) {
double 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;
@ -97,13 +96,14 @@ public:
inserted_vertices.push_back(new_vertex);
double source_weight = probability_matrix[alpha][vertex_i];
// since it is expansion move, current alpha labeled vertices will be assigned to alpha again,
// making sink_weight 'infinity' guarantee this.
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, source_weight,
graph);
add_edge_and_reverse(new_vertex, cluster_sink, sink_weight, sink_weight, graph);
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,
@ -115,17 +115,15 @@ public:
v2 = inserted_vertices[edge_it->second];
int label_1 = labels[edge_it->first], label_2 = labels[edge_it->second];
if(label_1 == label_2) {
// actually no need for this, since two alpha labeled vertices will not be seperated
// (their edges between sink is infitinity)
double w1 = (label_1 == alpha) ? 0 : *weight_it;
//double w1 = *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;
//double w1 = *weight_it;
//double w2 = *weight_it;
add_edge_and_reverse(inbetween, v1, w1, w1, graph);
add_edge_and_reverse(inbetween, v2, w2, w2, graph);
add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
@ -144,7 +142,8 @@ public:
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
if(labels[vertex_i] != alpha
&& color == ColorTraits::white()) { //new comers (expansion occurs)
labels[vertex_i] = alpha;
}
}
@ -153,171 +152,6 @@ public:
} while(success);
return min_cut;
}
#if 0
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) {
#if 1
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);
#else
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;
double w1 = *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;
double w1 = *weight_it;
double w2 = *weight_it;
add_edge_and_reverse(inbetween, v1, w1, w1, graph);
add_edge_and_reverse(inbetween, v2, w2, w2, graph);
//add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
}
#endif
}
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);
Vertex_descriptor cluster_sink = boost::add_vertex(graph);
std::vector<Vertex_descriptor> inserted_vertices;
//inserted_vertices.reserve(vertex_weights.size());
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);
Edge_descriptor source_e, source_e_rev, sink_e, sink_e_rev;
bool source_e_added, source_e_rev_added, sink_e_added, sink_e_rev_added;
tie(source_e, source_e_added) = boost::add_edge(cluster_source, new_vertex,
graph);
tie(source_e_rev, source_e_rev_added) = boost::add_edge(new_vertex,
cluster_source, graph);
tie(sink_e, sink_e_added) = boost::add_edge(new_vertex, cluster_sink, graph);
tie(sink_e_rev, sink_e_rev_added) = boost::add_edge(cluster_sink, new_vertex,
graph);
CGAL_assertion(source_e_added && source_e_rev_added && sink_e_added
&& sink_e_rev_added);
//put vertex capacities (to edges between itself and source & sink)
boost::put(boost::edge_capacity, graph, source_e,
probability_matrix[0][vertex_i]);
boost::put(boost::edge_capacity, graph, source_e_rev, 0);
boost::put(boost::edge_capacity, graph, sink_e,
probability_matrix[1][vertex_i]);
boost::put(boost::edge_capacity, graph, sink_e_rev, 0);
//map reverse edges
boost::put(boost::edge_reverse, graph, source_e, source_e_rev);
boost::put(boost::edge_reverse, graph, source_e_rev, source_e);
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>::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];
bool edge_added, edge_rev_added;
Edge_descriptor edge, edge_rev;
tie(edge, edge_added) = boost::add_edge(v1, v2, graph);
tie(edge_rev, edge_rev_added) = boost::add_edge(v2, v1, graph);
CGAL_assertion(edge_added && edge_rev_added);
//put edge capacities
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);
for(std::vector<Vertex_descriptor>::iterator vertex_it =
inserted_vertices.begin();
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
center_ids.push_back(1);
} else {
center_ids.push_back(0);
}
}
}
#endif
};
}//namespace internal
}//namespace CGAL

16
Surface_mesh_segmentation/include/CGAL/internal/Surface_mesh_segmentation/Alpha_expansion_graph_cut_with_EM.h
View File

@ -122,16 +122,12 @@ public:
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;
//double w1 = *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;
//double w1 = *weight_it;
//double w2 = *weight_it;
add_edge_and_reverse(inbetween, v1, w1, w1, graph);
add_edge_and_reverse(inbetween, v2, w2, w2, graph);
add_edge_and_reverse(inbetween, cluster_sink, *weight_it, 0.0, graph);
@ -156,13 +152,13 @@ public:
}
}
labels_result = labels;
}
Expectation_maximization em(probability_matrix.size(), data, labels);
labels.clear();
probability_matrix.clear();
em.fill_with_center_ids(labels);
em.fill_with_minus_log_probabilities(probability_matrix);
}
int matrix_size = probability_matrix.size();
probability_matrix.clear();
Expectation_maximization(data).refresh_parameters_and_probabilities(matrix_size,
labels, probability_matrix);
} while(success);
// recalculate clustering
labels = labels_result;

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

@ -19,6 +19,7 @@
#define DEF_MAX_ITER 10
#define DEF_THRESHOLD 1e-4
#define USE_MATRIX true
#define DEF_SEED 1340818006
#define ACTIVATE_SEGMENTATION_EM_DEBUG
#ifdef ACTIVATE_SEGMENTATION_EM_DEBUG
@ -65,19 +66,22 @@ public:
double threshold;
int maximum_iteration;
bool is_converged;
double final_likelihood;
protected:
unsigned int seed;
std::vector<std::vector<double> > membership_matrix;
public:
Expectation_maximization() { }
Expectation_maximization(const std::vector<double>& data): points(data) { }
/* 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)
: points(data), threshold(threshold), maximum_iteration(maximum_iteration),
is_converged(false), seed(0),
membership_matrix(std::vector<std::vector<double> >(number_of_centers,
std::vector<double>(points.size()))) {
std::vector<double>(points.size()))),
final_likelihood(-(std::numeric_limits<double>::max)()) {
initiate_centers_uniformly(number_of_centers);
fit();
}
@ -86,7 +90,8 @@ public:
std::vector<int>& initial_center_ids,
double threshold = DEF_THRESHOLD, int maximum_iteration = DEF_MAX_ITER)
: points(data), threshold(threshold), maximum_iteration(maximum_iteration),
is_converged(false), seed(0) {
is_converged(false), seed(0),
final_likelihood(-(std::numeric_limits<double>::max)()) {
number_of_centers = initiate_centers_from_memberships(number_of_centers,
initial_center_ids);
membership_matrix = std::vector<std::vector<double> >(number_of_centers,
@ -101,8 +106,13 @@ public:
is_converged(false),
seed(static_cast<unsigned int>(time(NULL))),
membership_matrix(std::vector<std::vector<double> >(number_of_centers,
std::vector<double>(points.size()))) {
std::vector<double>(points.size()))),
final_likelihood(-(std::numeric_limits<double>::max)()) {
#if 0
srand(seed);
#else
srand(DEF_SEED);
#endif
fit_with_multiple_run(number_of_centers, number_of_runs);
}
@ -139,9 +149,9 @@ public:
}
for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) {
double probability = probabilities[center_i][point_i] / total_probability;
probability = (CGAL::max)(probability, epsilon);
probability = (std::max)(probability, epsilon);
probability = -log(probability);
probabilities[center_i][point_i] = (CGAL::max)(probability,
probabilities[center_i][point_i] = (std::max)(probability,
1e-5); // this is another epsilon, edge can not hold 0 to source in graph-cut.
}
}
@ -173,12 +183,82 @@ public:
// {
// double probability = probabilities[center_i][point_i];
// probability = (probability - min_value) / max_min_dif;
// probability = (CGAL::max)(probability, epsilon);
// probability = (std::max)(probability, epsilon);
// probabilities[center_i][point_i] = -log(probability);
// }
// #endif
//}
}
// experimental, going to be removed most prob !
void refresh_parameters_and_probabilities(int number_of_centers,
std::vector<int>& initial_center_ids,
std::vector<std::vector<double> >& probabilities) {
std::vector<bool> cluster_exist(number_of_centers, false);
std::vector<int> cluster_shift(number_of_centers, 0);
for(std::vector<int>::iterator id_it = initial_center_ids.begin();
id_it != initial_center_ids.end(); ++id_it) {
cluster_exist[*id_it] = true;
}
/* Calculate mean */
int number_of_points = initial_center_ids.size();
centers = std::vector<Gaussian_center>(number_of_centers);
std::vector<int> member_count(number_of_centers, 0);
for(int i = 0; i < number_of_points; ++i) {
int center_id = initial_center_ids[i];
double data = points[i];
centers[center_id].mean += data;
member_count[center_id] += 1;
}
/* Assign mean, and mixing coef */
for(int i = 0; i < number_of_centers; ++i) {
if(!cluster_exist[i]) {
continue;
}
centers[i].mean /= member_count[i];
centers[i].mixing_coefficient = member_count[i] / static_cast<double>
(number_of_points);
}
/* Calculate deviation */
for(int i = 0; i < number_of_points; ++i) {
int center_id = initial_center_ids[i];
double data = points[i];
centers[center_id].deviation += pow(data - centers[center_id].mean, 2);
}
for(int i = 0; i < number_of_centers; ++i) {
if(!cluster_exist[i]) {
continue;
}
centers[i].deviation = sqrt(centers[i].deviation / member_count[i]);
}
double epsilon = 1e
-5; // this epsilon should be consistent with epsilon in calculate_dihedral_angle_of_edge!
probabilities = std::vector<std::vector<double> >
(number_of_centers, std::vector<double>(points.size()));
for(std::size_t point_i = 0; point_i < points.size(); ++point_i) {
double total_probability = 0.0;
for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) {
if(!cluster_exist[center_i]) {
probabilities[center_i][point_i] = 0;
continue;
}
double probability = centers[center_i].probability(points[point_i]);
total_probability += probability;
probabilities[center_i][point_i] = probability;
}
for(std::size_t center_i = 0; center_i < centers.size(); ++center_i) {
double probability = probabilities[center_i][point_i] / total_probability;
probability = (std::max)(probability, epsilon);
probability = -log(probability);
probabilities[center_i][point_i] = (std::max)(probability,
1e-5); // this is another epsilon, edge can not hold 0 to source in graph-cut.
}
}
//sort(centers.begin(), centers.end());
}
protected:
// Initialization functions for centers.
@ -213,6 +293,7 @@ protected:
sort(centers.begin(), centers.end());
}
int initiate_centers_from_memberships(int number_of_centers,
std::vector<int>& initial_center_ids) {
/* For handling clusters that have 0 members (in initial_center_ids) */
@ -336,29 +417,28 @@ protected:
while(!is_converged && iteration_count++ < maximum_iteration) {
prev_likelihood = likelihood;
likelihood = iterate(iteration_count == 1);
is_converged = likelihood - prev_likelihood < threshold * fabs(likelihood);
is_converged = likelihood - prev_likelihood < threshold * std::fabs(likelihood);
}
SEG_DEBUG(std::cout << "likelihood: " << likelihood << "iteration: " <<
iteration_count << std::endl)
if(final_likelihood < likelihood) {
final_likelihood = likelihood;
}
return likelihood;
}
// Calls fit() number_of_run times, and stores best found centers
void fit_with_multiple_run(int number_of_centers, int number_of_run) {
double max_likelihood = -(std::numeric_limits<double>::max)();
std::vector<Gaussian_center> max_centers;
while(number_of_run-- > 0) {
initiate_centers_randomly(number_of_centers);
double likelihood = fit();
//write_center_parameters("center_paramters.txt");
if(likelihood > max_likelihood) {
if(likelihood == final_likelihood) {
max_centers = centers;
max_likelihood = likelihood;
}
}
SEG_DEBUG(std::cout << "max likelihood: " << max_likelihood << std::endl)
SEG_DEBUG(std::cout << "max likelihood: " << final_likelihood << std::endl)
centers = max_centers;
}
@ -383,6 +463,7 @@ protected:
};
}//namespace internal
}//namespace CGAL
#undef DEF_SEED
#undef DEF_MAX_ITER
#undef DEF_THRESHOLD
#undef USE_MATRIX

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

@ -6,7 +6,9 @@
#include <ctime>
#include <cstdlib>
#include <limits>
#include <algorithm>
#define DEF_SEED 1340818006
//#define ACTIVATE_SEGMENTATION_K_MEANS_DEBUG
#ifdef ACTIVATE_SEGMENTATION_K_MEANS_DEBUG
#define SEG_DEBUG(x) x;
@ -59,9 +61,9 @@ public:
* Returns true if center is changed.*/
bool calculate_new_center(std::vector<K_means_center>& centers) {
int new_center_id = 0;
double min_distance = fabs(centers[0].mean - data);
double min_distance = std::fabs(centers[0].mean - data);
for(std::size_t i = 1; i < centers.size(); ++i) {
double new_distance = fabs(centers[i].mean - data);
double new_distance = std::fabs(centers[i].mean - data);
if(new_distance < min_distance) {
new_center_id = i;
min_distance = new_distance;
@ -87,11 +89,15 @@ protected:
public:
K_means_clustering(int number_of_centers, const std::vector<double>& data,
int number_of_run = 100, int maximum_iteration = 100)
int number_of_run = 20, int maximum_iteration = 100)
: points(data.begin(), data.end()), maximum_iteration(maximum_iteration),
is_converged(false),
seed(static_cast<unsigned int>(time(NULL))) {
#if 0
srand(seed);
#else
srand(DEF_SEED);
#endif
calculate_clustering_with_multiple_run(number_of_centers, number_of_run);
}
/* For each data point, data_center is filled by its center's id. */
@ -121,6 +127,35 @@ public:
sort(centers.begin(), centers.end());
}
void initiate_centers_k_means_plus_plus(int number_of_centers) {
centers.clear();
std::vector<double> distance_square_cumulative(points.size());
std::vector<double> distance_square(points.size(),
(std::numeric_limits<double>::max)());
double initial_mean = points[rand() % points.size()].data;
centers.push_back(K_means_center(initial_mean));
for(int i = 1; i < number_of_centers; ++i) {
double cumulative_distance_square = 0.0;
for(std::size_t j = 0; j < points.size(); ++j) {
double new_distance = pow(centers.back().mean - points[j].data, 2);
if(new_distance < distance_square[j]) {
distance_square[j] = new_distance;
}
cumulative_distance_square += distance_square[j];
distance_square_cumulative[j] = cumulative_distance_square;
}
double random_ds = (rand() / static_cast<double>(RAND_MAX)) *
(distance_square_cumulative.back());
int selection_index = lower_bound(distance_square_cumulative.begin(),
distance_square_cumulative.end(), random_ds)
- distance_square_cumulative.begin();
double initial_mean = points[selection_index].data;
centers.push_back(K_means_center(initial_mean));
}
}
void calculate_clustering() {
int iteration_count = 0;
bool any_center_changed = true;
@ -149,8 +184,12 @@ public:
double error = (std::numeric_limits<double>::max)();
while(number_of_run-- > 0) {
clear_center_ids();
#if 0
initiate_centers_randomly(number_of_centers);
#else
initiate_centers_k_means_plus_plus(number_of_centers);
#endif
calculate_clustering();
#if 0
double new_error = sum_of_squares();
@ -198,6 +237,7 @@ public:
};
}//namespace internal
}//namespace CGAL
#undef DEF_SEED
#ifdef SEG_DEBUG
#undef SEG_DEBUG
#endif