mirror of
https://github.com/truebelief/cc-treeiso-plugin.git
synced 2026-08-31 01:20:29 +08:00
Daniel Girardeau-Montaut
707 lines
34 KiB
C++
707 lines
34 KiB
C++
#pragma once
|
|
|
|
//Local
|
|
#include "Graph.h"
|
|
|
|
//System
|
|
#include <math.h>
|
|
#include <queue>
|
|
#include <iostream>
|
|
#include <fstream>
|
|
|
|
//Boost
|
|
#include <boost/graph/boykov_kolmogorov_max_flow.hpp>
|
|
|
|
namespace CP
|
|
{
|
|
template <typename T>
|
|
struct CPparameter
|
|
{
|
|
T reg_strenth; //regularization strength, multiply the edge weight
|
|
uint32_t cutoff; //minimal component size
|
|
uint32_t flow_steps; //number of steps in the optimal binary cut computation
|
|
uint32_t kmeans_ite; //number of iteration in the kmeans sampling
|
|
uint32_t kmeans_resampling; //number of kmeans re-intilialization
|
|
uint32_t verbose; //verbosity
|
|
uint32_t max_ite_main; //max number of iterations in the main loop
|
|
bool backward_step; //indicates if a backward step should be performed
|
|
double stopping_ratio; //when (E(t-1) - E(t) / (E(0) - E(t)) is too small, the algorithm stops
|
|
fidelityType fidelity; //the fidelity function
|
|
double smoothing; //smoothing term (for Kl divergence only)
|
|
bool parallel; //enable/disable parrallelism
|
|
T weight_decay; //for continued optimization of the flow steps
|
|
};
|
|
|
|
template <typename T>
|
|
struct CutPursuit
|
|
{
|
|
Graph<T> main_graph; //the Graph structure containing the main structure
|
|
Graph<T> reduced_graph; //the reduced graph whose vertices are the connected component
|
|
std::vector<std::vector<VertexDescriptor<T>>> components; //contains the list of the vertices in each component
|
|
std::vector<VertexDescriptor<T>> root_vertex; //the root vertex for each connected components
|
|
std::vector<bool> saturated_components; //is the component saturated (uncuttable)
|
|
std::vector<std::vector<EdgeDescriptor>> borders; //the list of edges forming the borders between the connected components
|
|
VertexDescriptor<T> source; //source vertex for graph cut
|
|
VertexDescriptor<T> sink; //sink vertex
|
|
uint32_t dim; // dimension of the data
|
|
uint32_t nVertex; // number of data point
|
|
uint32_t nEdge; // number of edges between vertices (not counting the edge to source/sink)
|
|
CP::VertexIterator<T> lastIterator; //iterator pointing to the last vertex which is neither sink nor source
|
|
CPparameter<T> parameter;
|
|
|
|
CutPursuit(uint32_t nbVertex = 1)
|
|
{
|
|
this->main_graph = Graph<T>(nbVertex);
|
|
this->reduced_graph = Graph<T>(1);
|
|
this->components = std::vector<std::vector<VertexDescriptor<T>>>(1);
|
|
this->root_vertex = std::vector<VertexDescriptor<T>>(1, 0);
|
|
this->saturated_components = std::vector<bool>(1, false);
|
|
this->source = VertexDescriptor<T>();
|
|
this->sink = VertexDescriptor<T>();
|
|
this->dim = 1;
|
|
this->nVertex = 1;
|
|
this->nEdge = 0;
|
|
this->parameter.flow_steps = 3;
|
|
this->parameter.kmeans_ite = 5;
|
|
this->parameter.kmeans_resampling = 3;
|
|
this->parameter.verbose = 2;
|
|
this->parameter.max_ite_main = 6;
|
|
this->parameter.backward_step = true;
|
|
this->parameter.stopping_ratio = 0.0001;
|
|
this->parameter.fidelity = L2;
|
|
this->parameter.smoothing = 0.1;
|
|
this->parameter.parallel = true;
|
|
this->parameter.weight_decay = static_cast<T>(0.7);
|
|
}
|
|
|
|
//=============================================================================================
|
|
std::pair<std::vector<T>, std::vector<T>> run()
|
|
{
|
|
//first initilialize the structure
|
|
this->initialize();
|
|
if (this->parameter.verbose > 0)
|
|
{
|
|
std::cout << "Graph " << boost::num_vertices(this->main_graph) << " vertices and "
|
|
<< boost::num_edges(this->main_graph) << " edges and observation of dimension "
|
|
<< this->dim << '\n';
|
|
}
|
|
T energy_zero = this->compute_energy().first; //energy with 1 component
|
|
T old_energy = energy_zero; //energy at the previous iteration
|
|
//vector with time and energy, useful for benchmarking
|
|
std::vector<T> energy_out(this->parameter.max_ite_main), time_out(this->parameter.max_ite_main);
|
|
TimeStack ts; ts.tic();
|
|
//the main loop
|
|
for (uint32_t ite_main = 1; ite_main <= this->parameter.max_ite_main; ite_main++)
|
|
{
|
|
//--------those two lines are the whole iteration-------------------------
|
|
size_t saturation = this->split(); //compute optimal binary partition
|
|
this->reduce(); //compute the new reduced graph
|
|
//-------end of the iteration - rest is stopping check and display------
|
|
std::pair<T, T> energy = this->compute_energy();
|
|
energy_out.push_back((energy.first + energy.second));
|
|
time_out.push_back(ts.tocDouble());
|
|
if (this->parameter.verbose > 1)
|
|
{
|
|
printf("Iteration %3i - %4i components - ", ite_main, static_cast<int>(this->components.size()));
|
|
printf("Saturation %5.1f %% - ", (100.0 * saturation) / this->nVertex);
|
|
switch (this->parameter.fidelity)
|
|
{
|
|
case L2:
|
|
{
|
|
printf("Quadratic Energy %4.3f %% - ", 100 * (energy.first + energy.second) / energy_zero);
|
|
break;
|
|
}
|
|
case linear:
|
|
{
|
|
printf("Linear Energy %10.1f - ", energy.first + energy.second);
|
|
break;
|
|
}
|
|
case KL:
|
|
{
|
|
printf("KL Energy %4.3f %% - ", 100 * (energy.first + energy.second) / energy_zero);
|
|
break;
|
|
}
|
|
case SPG:
|
|
{
|
|
printf("Quadratic Energy %4.3f %% - ", 100 * (energy.first + energy.second) / energy_zero);
|
|
break;
|
|
}
|
|
}
|
|
std::cout << "Timer " << ts.toc() << std::endl;
|
|
}
|
|
//----stopping checks-----
|
|
if (saturation == this->nVertex)
|
|
{ //all components are saturated
|
|
if (this->parameter.verbose > 1)
|
|
{
|
|
std::cout << "All components are saturated" << std::endl;
|
|
}
|
|
break;
|
|
}
|
|
if ((old_energy - energy.first - energy.second) / (old_energy)
|
|
< this->parameter.stopping_ratio)
|
|
{ //relative energy progress stopping criterion
|
|
if (this->parameter.verbose > 1)
|
|
{
|
|
std::cout << "Stopping criterion reached" << std::endl;
|
|
}
|
|
break;
|
|
}
|
|
if (ite_main >= this->parameter.max_ite_main)
|
|
{ //max number of iteration
|
|
if (this->parameter.verbose > 1)
|
|
{
|
|
std::cout << "Max number of iteration reached" << std::endl;
|
|
}
|
|
break;
|
|
}
|
|
old_energy = energy.first + energy.second;
|
|
}
|
|
if (this->parameter.cutoff > 0)
|
|
{
|
|
this->cutoff();
|
|
}
|
|
return std::pair<std::vector<T>, std::vector<T>>(energy_out, time_out);
|
|
}
|
|
|
|
//=============================================================================================
|
|
//=========== VIRTUAL METHODS DEPENDING ON THE CHOICE OF FIDELITY FUNCTION =====================
|
|
//=============================================================================================
|
|
//
|
|
//=============================================================================================
|
|
//============================= SPLIT ===========================================
|
|
//=============================================================================================
|
|
virtual size_t split()
|
|
{
|
|
//compute the optimal binary partition
|
|
return 0;
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================ compute_energy_L2 ====================================
|
|
//=============================================================================================
|
|
virtual std::pair<T, T> compute_energy()
|
|
{
|
|
//compute the current energy
|
|
return std::pair<T, T>(0, 0);
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================= COMPUTE_VALUE =========================================
|
|
//=============================================================================================
|
|
virtual std::pair<std::vector<T>, T> compute_value(const uint32_t & ind_com)
|
|
{
|
|
//compute the optimal the values associated with the current partition
|
|
return std::pair<std::vector<T>, T>(std::vector<T>(0), 0);
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================= COMPUTE_MERGE_GAIN =========================================
|
|
//=============================================================================================
|
|
virtual std::pair<std::vector<T>, T> compute_merge_gain(const VertexDescriptor<T> & comp1
|
|
, const VertexDescriptor<T> & comp2)
|
|
{
|
|
//compute the gain of mergeing two connected components
|
|
return std::pair<std::vector<T>, T>(std::vector<T>(0), 0);
|
|
}
|
|
|
|
//=============================================================================================
|
|
//========================== END OF VIRTUAL METHODS ===========================================
|
|
//=============================================================================================
|
|
|
|
//=============================================================================================
|
|
//============================= INITIALIZE ===========================================
|
|
//=============================================================================================
|
|
void initialize()
|
|
{
|
|
//build the reduced graph with one component, fill the first vector of components
|
|
//and add the sink and source nodes
|
|
VertexIterator<T> ite_ver, ite_ver_end;
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
this->components[0] = std::vector<VertexDescriptor<T>>(0);//(this->nVertex);
|
|
this->root_vertex[0] = *boost::vertices(this->main_graph).first;
|
|
this->nVertex = static_cast<uint32_t>(boost::num_vertices(this->main_graph));
|
|
this->nEdge = static_cast<uint32_t>(boost::num_edges(this->main_graph));
|
|
//--------compute the first reduced graph----------------------------------------------------------
|
|
for (boost::tie(ite_ver, ite_ver_end) = boost::vertices(this->main_graph);
|
|
ite_ver != ite_ver_end; ++ite_ver)
|
|
{
|
|
this->components[0].push_back(*ite_ver);
|
|
}
|
|
this->lastIterator = ite_ver;
|
|
this->compute_value(0);
|
|
//--------build the link to source and sink--------------------------------------------------------
|
|
this->source = boost::add_vertex(this->main_graph);
|
|
this->sink = boost::add_vertex(this->main_graph);
|
|
uint32_t eIndex = static_cast<uint32_t>(boost::num_edges(this->main_graph));
|
|
ite_ver = boost::vertices(this->main_graph).first;
|
|
for (uint32_t ind_ver = 0; ind_ver < this->nVertex; ind_ver++)
|
|
{
|
|
// note that source and edge will have many nieghbors, and hence boost::edge should never be called to get
|
|
// the in_edge. use the out_edge and then reverse_Edge
|
|
addDoubledge<T>(this->main_graph, this->source, boost::vertex(ind_ver, this->main_graph), 0.,
|
|
eIndex, edge_attribute_map, false);
|
|
eIndex += 2;
|
|
addDoubledge<T>(this->main_graph, boost::vertex(ind_ver, this->main_graph), this->sink, 0.,
|
|
eIndex, edge_attribute_map, false);
|
|
eIndex += 2;
|
|
++ite_ver;
|
|
}
|
|
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================ COMPUTE_REDUCE_VALUE ====================================
|
|
//=============================================================================================
|
|
void compute_reduced_value()
|
|
{
|
|
for (uint32_t ind_com = 0; ind_com < this->components.size(); ++ind_com)
|
|
{ //compute the reduced value of each component
|
|
compute_value(ind_com);
|
|
}
|
|
}
|
|
|
|
//=============================================================================================
|
|
//============================= ACTIVATE_EDGES ==========================================
|
|
//=============================================================================================
|
|
size_t activate_edges(bool allows_saturation = true)
|
|
{ //this function analyzes the optimal binary partition to detect:
|
|
//- saturated components (i.e. uncuttable)
|
|
//- new activated edges
|
|
VertexAttributeMap<T> vertex_attribute_map = boost::get(boost::vertex_bundle, this->main_graph);
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
//saturation is the proportion of nodes in saturated components
|
|
size_t saturation = 0;
|
|
uint32_t nb_comp = static_cast<uint32_t>(this->components.size());
|
|
//---- first check if the component are saturated-------------------------
|
|
//#pragma omp parallel for if (this->parameter.parallel) schedule(dynamic)
|
|
for (uint32_t ind_com = 0; ind_com < nb_comp; ind_com++)
|
|
{
|
|
if (this->saturated_components[ind_com])
|
|
{ //ind_com is saturated, we increement saturation by ind_com size
|
|
saturation += this->components[ind_com].size();
|
|
continue;
|
|
}
|
|
std::vector<T> totalWeight(2, 0);
|
|
for (uint32_t ind_ver = 0; ind_ver < this->components[ind_com].size(); ind_ver++)
|
|
{
|
|
bool isSink
|
|
= (vertex_attribute_map(this->components[ind_com][ind_ver]).color
|
|
== vertex_attribute_map(this->sink).color);
|
|
if (isSink)
|
|
{
|
|
totalWeight[0] += vertex_attribute_map(this->components[ind_com][ind_ver]).weight;
|
|
}
|
|
else
|
|
{
|
|
totalWeight[1] += vertex_attribute_map(this->components[ind_com][ind_ver]).weight;
|
|
}
|
|
}
|
|
if (allows_saturation && ((totalWeight[0] == 0) || (totalWeight[1] == 0)))
|
|
{
|
|
//the component is saturated
|
|
this->saturateComponent(ind_com);
|
|
saturation += this->components[ind_com].size();
|
|
}
|
|
}
|
|
//----check which edges have been activated----
|
|
EdgeIterator<T> ite_edg, ite_edg_end;
|
|
uint32_t color_v1, color_v2, color_combination;
|
|
for (boost::tie(ite_edg, ite_edg_end) = boost::edges(this->main_graph);
|
|
ite_edg != ite_edg_end; ++ite_edg)
|
|
{
|
|
if (!edge_attribute_map(*ite_edg).realEdge)
|
|
{
|
|
continue;
|
|
}
|
|
color_v1 = vertex_attribute_map(boost::source(*ite_edg, this->main_graph)).color;
|
|
color_v2 = vertex_attribute_map(boost::target(*ite_edg, this->main_graph)).color;
|
|
//color_source = 0, color_sink = 4, uncolored = 1
|
|
//we want an edge when a an interface source/sink
|
|
//this corresponds to a sum of 4
|
|
//for the case of uncolored nodes we arbitrarily chose source-uncolored
|
|
color_combination = color_v1 + color_v2;
|
|
if ((color_combination == 0) || (color_combination == 2) || (color_combination == 2)
|
|
|| (color_combination == 8))
|
|
{ //edge between two vertices of the same color
|
|
continue;
|
|
}
|
|
//the edge is active!
|
|
edge_attribute_map(*ite_edg).isActive = true;
|
|
edge_attribute_map(*ite_edg).capacity = 0;
|
|
vertex_attribute_map(boost::source(*ite_edg, this->main_graph)).isBorder = true;
|
|
vertex_attribute_map(boost::target(*ite_edg, this->main_graph)).isBorder = true;
|
|
}
|
|
return saturation;
|
|
}
|
|
|
|
//=============================================================================================
|
|
//============================= REDUCE ===========================================
|
|
//=============================================================================================
|
|
void reduce()
|
|
{ //compute the reduced graph, and if need be performed a backward check
|
|
this->compute_connected_components();
|
|
if (this->parameter.backward_step)
|
|
{ //compute the structure of the reduced graph
|
|
this->compute_reduced_graph();
|
|
//check for beneficial merges
|
|
this->merge(false);
|
|
}
|
|
else
|
|
{ //compute only the value associated to each connected components
|
|
this->compute_reduced_value();
|
|
}
|
|
}
|
|
|
|
//=============================================================================================
|
|
//============================== compute_connected_components=========================================
|
|
//=============================================================================================
|
|
void compute_connected_components()
|
|
{ //this function compute the connected components of the graph with active edges removed
|
|
//the boolean vector indicating wether or not the edges and vertices have been seen already
|
|
//the root is the first vertex of a component
|
|
//this function is written such that the new components are appended at the end of components
|
|
//this allows not to recompute saturated component
|
|
VertexAttributeMap<T> vertex_attribute_map = boost::get(boost::vertex_bundle, this->main_graph);
|
|
VertexIndexMap<T> vertex_index_map = get(boost::vertex_index, this->main_graph);
|
|
//indicate which edges and nodes have been seen already by the dpsearch
|
|
std::vector<bool> edges_seen(this->nEdge, false);
|
|
std::vector<bool> vertices_seen(this->nVertex + 2, false);
|
|
vertices_seen[vertex_index_map(this->source)] = true;
|
|
vertices_seen[vertex_index_map(this->sink)] = true;
|
|
//-------- start with the known roots------------------------------------------------------
|
|
//#pragma omp parallel for if (this->parameter.parallel) schedule(dynamic)
|
|
for (uint32_t ind_com = 0; ind_com < this->root_vertex.size(); ind_com++)
|
|
{
|
|
VertexDescriptor<T> root = this->root_vertex[ind_com]; //the first vertex of the component
|
|
if (this->saturated_components[ind_com])
|
|
{ //this component is saturated, we don't need to recompute it
|
|
for (uint32_t ind_ver = 0; ind_ver < this->components[ind_com].size(); ++ind_ver)
|
|
{
|
|
vertices_seen[vertex_index_map(this->components[ind_com][ind_ver])] = true;
|
|
}
|
|
}
|
|
else
|
|
{ //compute the new content of this component
|
|
this->components.at(ind_com) = connected_comp_from_root(root, this->components.at(ind_com).size(), vertices_seen, edges_seen);
|
|
}
|
|
}
|
|
//----now look for components that did not already exists----
|
|
VertexIterator<T> ite_ver;
|
|
for (ite_ver = boost::vertices(this->main_graph).first;
|
|
ite_ver != this->lastIterator; ++ite_ver)
|
|
{
|
|
if (vertices_seen[vertex_index_map(*ite_ver)])
|
|
{
|
|
continue;
|
|
} //this vertex is not currently in a connected component
|
|
VertexDescriptor<T> root = *ite_ver; //we define it as the root of a new component
|
|
size_t current_component_size = this->components[vertex_attribute_map(root).in_component].size();
|
|
this->components.push_back(connected_comp_from_root(root, current_component_size, vertices_seen, edges_seen));
|
|
this->root_vertex.push_back(root);
|
|
this->saturated_components.push_back(false);
|
|
}
|
|
this->components.shrink_to_fit();
|
|
}
|
|
|
|
//=============================================================================================
|
|
//============================== CONNECTED_COMP_FROM_ROOT=========================================
|
|
//=============================================================================================
|
|
inline std::vector<VertexDescriptor<T>> connected_comp_from_root(const VertexDescriptor<T>& root
|
|
, const size_t& size_comp, std::vector<bool>& vertices_seen, std::vector<bool>& edges_seen)
|
|
{
|
|
//this function compute the connected component of the graph with active edges removed
|
|
// associated with the root ROOT by performing a depth search first
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
VertexIndexMap<T> vertex_index_map = get(boost::vertex_index, this->main_graph);
|
|
EdgeIndexMap<T> edge_index_map = get(&EdgeAttribute<T>::index, this->main_graph);
|
|
std::vector<VertexDescriptor<T>> vertices_added; //the vertices in the current connected component
|
|
// vertices_added contains the vertices that have been added to the current coomponent
|
|
vertices_added.reserve(size_comp);
|
|
//heap_explore contains the vertices to be added to the current component
|
|
std::vector<VertexDescriptor<T>> vertices_to_add;
|
|
vertices_to_add.reserve(size_comp);
|
|
VertexDescriptor<T> vertex_current; //the node being consideed
|
|
EdgeDescriptor edge_current, edge_reverse; //the edge being considered
|
|
//fill the heap with the root node
|
|
vertices_to_add.push_back(root);
|
|
while (vertices_to_add.size() > 0)
|
|
{ //as long as there are vertices left to add
|
|
vertex_current = vertices_to_add.back(); //the current node is the last node to add
|
|
vertices_to_add.pop_back(); //remove the current node from the vertices to add
|
|
if (vertices_seen[vertex_index_map(vertex_current)])
|
|
{ //this vertex has already been treated
|
|
continue;
|
|
}
|
|
vertices_added.push_back(vertex_current); //we add the current node
|
|
vertices_seen[vertex_index_map(vertex_current)] = true; //and flag it as seen
|
|
//----we now explore the neighbors of current_node
|
|
typename boost::graph_traits<Graph<T>>::out_edge_iterator ite_edg, ite_edg_end;
|
|
for (boost::tie(ite_edg, ite_edg_end) = boost::out_edges(vertex_current, this->main_graph);
|
|
ite_edg != ite_edg_end; ++ite_edg)
|
|
{ //explore edges leaving current_node
|
|
edge_current = *ite_edg;
|
|
if (edge_attribute_map(*ite_edg).isActive || (edges_seen[edge_index_map(edge_current)]))
|
|
{ //edge is either active or treated, we skip it
|
|
continue;
|
|
}
|
|
//the target of this edge is a node to add
|
|
edge_reverse = edge_attribute_map(edge_current).edge_reverse;
|
|
edges_seen[edge_index_map(edge_current)] = true;
|
|
edges_seen[edge_index_map(edge_reverse)] = true;
|
|
vertices_to_add.push_back(boost::target(edge_current, this->main_graph));
|
|
}
|
|
}
|
|
vertices_added.shrink_to_fit();
|
|
return vertices_added;
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================ COMPUTE_REDUCE_GRAPH ====================================
|
|
//=============================================================================================
|
|
void compute_reduced_graph()
|
|
{ //compute the adjacency structure between components as well as weight and value of each component
|
|
//this is stored in the reduced graph structure
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
VertexAttributeMap<T> vertex_attribute_map = boost::get(boost::vertex_bundle, this->main_graph);
|
|
this->reduced_graph = Graph<T>(this->components.size());
|
|
VertexAttributeMap<T> component_attribute_map = boost::get(boost::vertex_bundle, this->reduced_graph);
|
|
//----fill the value sof the reduced graph----
|
|
#ifdef OPENMP
|
|
#pragma omp parallel for schedule(dynamic)
|
|
#endif
|
|
for (uint32_t ind_com = 0; ind_com < this->components.size(); ind_com++)
|
|
{
|
|
std::pair<std::vector<T>, T> component_values_and_weight = this->compute_value(ind_com);
|
|
//----fill the value and weight field of the reduced graph-----------------------------
|
|
VertexDescriptor<T> reduced_vertex = boost::vertex(ind_com, this->reduced_graph);
|
|
component_attribute_map[reduced_vertex] = VertexAttribute<T>(this->dim);
|
|
component_attribute_map(reduced_vertex).weight
|
|
= component_values_and_weight.second;
|
|
for (uint32_t i_dim = 0; i_dim < this->dim; i_dim++)
|
|
{
|
|
component_attribute_map(reduced_vertex).value[i_dim] = component_values_and_weight.first[i_dim];
|
|
}
|
|
}
|
|
//------compute the edges of the reduced graph
|
|
EdgeAttributeMap<T> border_edge_attribute_map = boost::get(boost::edge_bundle, this->reduced_graph);
|
|
this->borders.clear();
|
|
EdgeDescriptor edge_current, border_edge_current;
|
|
uint32_t ind_border_edge = 0, comp1, comp2, component_source, component_target;
|
|
VertexDescriptor<T> source_component, target_component;
|
|
bool reducedEdgeExists;
|
|
typename boost::graph_traits<Graph<T>>::edge_iterator ite_edg, ite_edg_end;
|
|
for (boost::tie(ite_edg, ite_edg_end) = boost::edges(this->main_graph); ite_edg != ite_edg_end; ++ite_edg)
|
|
{
|
|
if (!edge_attribute_map(*ite_edg).realEdge)
|
|
{ //edges linking the source or edge node do not take part
|
|
continue;
|
|
}
|
|
edge_current = *ite_edg;
|
|
//compute the connected components of the source and target of current_edge
|
|
comp1 = vertex_attribute_map(boost::source(edge_current, this->main_graph)).in_component;
|
|
comp2 = vertex_attribute_map(boost::target(edge_current, this->main_graph)).in_component;
|
|
if (comp1 == comp2)
|
|
{ //this edge links two nodes in the same connected component
|
|
continue;
|
|
}
|
|
//by convention we note component_source the smallest index and
|
|
//component_target the largest
|
|
component_source = std::min(comp1, comp2);
|
|
component_target = std::max(comp1, comp2);
|
|
//retrieve the corresponding vertex in the reduced graph
|
|
source_component = boost::vertex(component_source, this->reduced_graph);
|
|
target_component = boost::vertex(component_target, this->reduced_graph);
|
|
//try to add the border-edge linking those components in the reduced graph
|
|
boost::tie(border_edge_current, reducedEdgeExists)
|
|
= boost::edge(source_component, target_component, this->reduced_graph);
|
|
if (!reducedEdgeExists)
|
|
{ //this border-edge did not already existed in the reduced graph
|
|
//border_edge_current = boost::add_edge(source_component, target_component, this->reduced_graph).first;
|
|
border_edge_current = boost::add_edge(source_component, target_component, this->reduced_graph).first;
|
|
border_edge_attribute_map(border_edge_current).index = ind_border_edge;
|
|
border_edge_attribute_map(border_edge_current).weight = 0;
|
|
ind_border_edge++;
|
|
//create a new entry for the borders list containing this border
|
|
this->borders.push_back(std::vector<EdgeDescriptor>(0));
|
|
}
|
|
//add the weight of the current edge to the weight of the border-edge
|
|
border_edge_attribute_map(border_edge_current).weight += 0.5*edge_attribute_map(edge_current).weight;
|
|
this->borders[border_edge_attribute_map(border_edge_current).index].push_back(edge_current);
|
|
}
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================ MERGE ====================================
|
|
//=============================================================================================
|
|
uint32_t merge(bool is_cutoff)
|
|
{
|
|
// TODO: right now we only do one loop through the heap of potential mergeing, and only
|
|
//authorize one mergeing per component. We could update the gain and merge until it is no longer
|
|
//beneficial
|
|
//check wether the energy can be decreased by removing edges from the reduced graph
|
|
//----load graph structure---
|
|
VertexAttributeMap<T> vertex_attribute_map = boost::get(boost::vertex_bundle, this->main_graph);
|
|
VertexAttributeMap<T> component_attribute_map = boost::get(boost::vertex_bundle, this->reduced_graph);
|
|
EdgeAttributeMap<T> border_edge_attribute_map = boost::get(boost::edge_bundle, this->reduced_graph);
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
VertexIndexMap<T> component_index_map = boost::get(boost::vertex_index, this->reduced_graph);
|
|
//-----------------------------------
|
|
EdgeDescriptor border_edge_current;
|
|
typename boost::graph_traits<Graph<T>>::edge_iterator ite_border, ite_border_end;
|
|
typename std::vector<EdgeDescriptor>::iterator ite_border_edge;
|
|
VertexDescriptor<T> source_component, target_component;
|
|
uint32_t ind_source_component, ind_target_component, border_edge_currentIndex;
|
|
//gain_current is the vector of gains associated with each mergeing move
|
|
//std::vector<T> gain_current(boost::num_edges(this->reduced_graph));
|
|
//we store in merge_queue the potential mergeing with a priority on the potential gain
|
|
std::priority_queue<ComponentsFusion<T>, std::vector<ComponentsFusion<T>>, lessComponentsFusion<T>> merge_queue;
|
|
T gain; // the gain obtained by removing the border corresponding to the edge in the reduced graph
|
|
for (boost::tie(ite_border, ite_border_end) = boost::edges(this->reduced_graph); ite_border != ite_border_end; ++ite_border)
|
|
{
|
|
//a first pass go through all the edges in the reduced graph and compute the gain obtained by
|
|
//mergeing the corresponding vertices
|
|
border_edge_current = *ite_border;
|
|
border_edge_currentIndex = border_edge_attribute_map(border_edge_current).index;
|
|
//retrieve the two components corresponding to this border
|
|
source_component = boost::source(border_edge_current, this->reduced_graph);
|
|
target_component = boost::target(border_edge_current, this->reduced_graph);
|
|
if (is_cutoff && component_attribute_map(source_component).weight >= this->parameter.cutoff
|
|
&&component_attribute_map(target_component).weight >= this->parameter.cutoff)
|
|
{
|
|
continue;
|
|
}
|
|
ind_source_component = static_cast<uint32_t>(component_index_map(source_component));
|
|
ind_target_component = static_cast<uint32_t>(component_index_map(target_component));
|
|
//----now compute the gain of mergeing those two components-----
|
|
// compute the fidelity lost by mergeing the two connected components
|
|
std::pair<std::vector<T>, T> merge_gain = compute_merge_gain(source_component, target_component);
|
|
// the second part is due to the removing of the border
|
|
gain = merge_gain.second
|
|
+ border_edge_attribute_map(border_edge_current).weight * this->parameter.reg_strenth;
|
|
//mergeing_information store the indexes of the components as well as the edge index and the gain
|
|
//in a structure ordered by the gain
|
|
ComponentsFusion<T> mergeing_information(ind_source_component, ind_target_component, border_edge_currentIndex, gain);
|
|
mergeing_information.merged_value = merge_gain.first;
|
|
if (is_cutoff || gain > 0)
|
|
{ //it is beneficial to merge those two components
|
|
//we add them to the merge_queue
|
|
merge_queue.push(mergeing_information);
|
|
//gain_current.at(border_edge_currentIndex) = gain;
|
|
}
|
|
}
|
|
uint32_t n_merged = 0;
|
|
//----go through the priority queue of merges and perform them as long as it is beneficial---
|
|
//is_merged indicate which components no longer exists because they have been merged with a neighboring component
|
|
std::vector<bool> is_merged(this->components.size(), false);
|
|
//to_destroy indicates the components that are needed to be removed
|
|
std::vector<bool> to_destroy(this->components.size(), false);
|
|
while (merge_queue.size() > 0)
|
|
{ //loop through the potential mergeing and accept the ones that decrease the energy
|
|
ComponentsFusion<T> mergeing_information = merge_queue.top();
|
|
if (!is_cutoff && mergeing_information.merge_gain <= 0)
|
|
{ //no more mergeing provide a gain in energy
|
|
break;
|
|
}
|
|
merge_queue.pop();
|
|
if (is_merged.at(mergeing_information.comp1) || (is_merged.at(mergeing_information.comp2)))
|
|
{
|
|
//at least one of the components have already been merged
|
|
continue;
|
|
}
|
|
n_merged++;
|
|
//---proceed with the fusion of comp1 and comp2----
|
|
//add the vertices of comp2 to comp1
|
|
this->components[mergeing_information.comp1].insert(this->components[mergeing_information.comp1].end()
|
|
, components[mergeing_information.comp2].begin(), this->components[mergeing_information.comp2].end());
|
|
//if comp1 was saturated it might not be anymore
|
|
this->saturated_components[mergeing_information.comp1] = false;
|
|
//the new weight is the sum of both weights
|
|
component_attribute_map(mergeing_information.comp1).weight
|
|
+= component_attribute_map(mergeing_information.comp2).weight;
|
|
//the new value is already computed in mergeing_information
|
|
component_attribute_map(mergeing_information.comp1).value = mergeing_information.merged_value;
|
|
//we deactivate the border between comp1 and comp2
|
|
for (ite_border_edge = this->borders.at(mergeing_information.border_index).begin();
|
|
ite_border_edge != this->borders.at(mergeing_information.border_index).end(); ++ite_border_edge)
|
|
{
|
|
edge_attribute_map(*ite_border_edge).isActive = false;
|
|
}
|
|
is_merged.at(mergeing_information.comp1) = true;
|
|
is_merged.at(mergeing_information.comp2) = true;
|
|
to_destroy.at(mergeing_information.comp2) = true;
|
|
}
|
|
//we now rebuild the vectors components, rootComponents and saturated_components
|
|
std::vector<std::vector<VertexDescriptor<T>>> new_components;
|
|
std::vector<VertexDescriptor<T>> new_root_vertex;
|
|
std::vector<bool> new_saturated_components;
|
|
uint32_t ind_new_component = 0;
|
|
for (uint32_t ind_com = 0; ind_com < this->components.size(); ind_com++)
|
|
{
|
|
if (to_destroy.at(ind_com))
|
|
{ //this component has been removed
|
|
continue;
|
|
}//this components is kept
|
|
new_components.push_back(this->components.at(ind_com));
|
|
new_root_vertex.push_back(this->root_vertex.at(ind_com));
|
|
new_saturated_components.push_back(saturated_components.at(ind_com));
|
|
//if (is_merged.at(ind_com))
|
|
//{ //we need to update the value of the vertex in this component
|
|
for (uint32_t ind_ver = 0; ind_ver < this->components[ind_com].size(); ++ind_ver)
|
|
{
|
|
vertex_attribute_map(this->components[ind_com][ind_ver]).value
|
|
= component_attribute_map(boost::vertex(ind_com, this->reduced_graph)).value;
|
|
vertex_attribute_map(this->components[ind_com][ind_ver]).in_component
|
|
= ind_new_component;//ind_com;
|
|
}
|
|
//}
|
|
ind_new_component++;
|
|
}
|
|
this->components = new_components;
|
|
this->root_vertex = new_root_vertex;
|
|
this->saturated_components = new_saturated_components;
|
|
return n_merged;
|
|
}
|
|
|
|
//=============================================================================================
|
|
//================================ CUTOFF ====================================
|
|
//=============================================================================================
|
|
void cutoff()
|
|
{
|
|
int i = 0;
|
|
uint32_t n_merged;
|
|
while (true)
|
|
{
|
|
//this->compute_connected_components();
|
|
this->compute_reduced_graph();
|
|
n_merged = merge(true);
|
|
i++;
|
|
if (n_merged == 0 || i > 50)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
//===============================================================================================
|
|
//========================= saturateComponent ===================================================
|
|
//===============================================================================================
|
|
inline void saturateComponent(const uint32_t & ind_com)
|
|
{ //this component is uncuttable and needs to be removed from further graph-cuts
|
|
EdgeAttributeMap<T> edge_attribute_map = boost::get(boost::edge_bundle, this->main_graph);
|
|
this->saturated_components[ind_com] = true;
|
|
for (uint32_t i_ver = 0; i_ver < this->components[ind_com].size(); i_ver++)
|
|
{
|
|
VertexDescriptor<T> desc_v = this->components[ind_com][i_ver];
|
|
// because of the adjacency structure NEVER access edge (source,v) directly!
|
|
EdgeDescriptor edg_ver2source = boost::edge(desc_v, this->source, this->main_graph).first;
|
|
EdgeDescriptor edg_source2ver = edge_attribute_map(edg_ver2source).edge_reverse; //use edge_reverse instead
|
|
EdgeDescriptor edg_sink2ver = boost::edge(desc_v, this->sink, this->main_graph).first;
|
|
// we set the capacities of edges to source and sink to zero
|
|
edge_attribute_map(edg_source2ver).capacity = 0.;
|
|
edge_attribute_map(edg_sink2ver).capacity = 0.;
|
|
}
|
|
}
|
|
};
|
|
}
|