ACloudViewer  3.9.4
A Modern Library for 3D Data Processing
API.h
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // - CloudViewer: www.cloudViewer.org -
3 // ----------------------------------------------------------------------------
4 // Copyright (c) 2018-2024 www.cloudViewer.org
5 // SPDX-License-Identifier: MIT
6 // ----------------------------------------------------------------------------
7 
8 #pragma once
9 #ifdef _OPENMP
10 #include <omp.h>
11 #endif
12 #include "CutPursuit_KL.h"
13 #include "CutPursuit_L2.h"
14 #include "CutPursuit_Linear.h"
15 #include "CutPursuit_SPG.h"
16 //**********************************************************************************
17 //*******************************L0-CUT
18 // PURSUIT*************************************
19 //**********************************************************************************
20 // Greedy graph cut based algorithm to solve the generalized minimal
21 // partition problem
22 //
23 // Cut Pursuit: fast algorithms to learn piecewise constant functions on
24 // general weighted graphs, Loic Landrieu and Guillaume Obozinski,2016.
25 //
26 // Produce a piecewise constant approximation of signal $y$ structured
27 // by the graph G=(V,e,mu,w) with mu the node weight and w the edgeweight:
28 // argmin \sum_{i \IN V}{mu_i * phi(x_I, y_I)}
29 //+ \sum_{(i,j) \IN E}{w_{i,j} 1(x_I != x_J)}
30 //
31 // phi(X,Y) the fidelity function (3 are implemented)
32 //(x != y) the function equal to 1 if x!=y and 0 else
33 //
34 // LOIC LANDRIEU 2017
35 //
36 //=======================SYNTAX===================================================
37 //---------------REGULARIZATION---------------------------------------------------
38 // C style inputs
39 // void cut_pursuit(const uint32_t n_nodes, const uint32_t n_edges, const
40 // uint32_t nObs
41 // ,const T * observation, const uint32_t * Eu, const uint32_t * Ev
42 // ,const T * edgeWeight, const T * nodeWeight
43 // ,T * solution, const T lambda, const uint32_t cutoff, const T mode,
44 // const T speed, const T weight_decay , const float verbose)
45 // C++ style input
46 // void cut_pursuit(const uint32_t n_nodes, const uint32_t n_edges, const
47 // uint32_t nObs
48 // , std::vector< std::vector<T> > & observation
49 // , const std::vector<uint32_t> & Eu, const std::vector<uint32_t> & Ev
50 // ,const std::vector<T> & edgeWeight, const std::vector<T> &
51 // nodeWeight ,std::vector< std::vector<T> > & solution, const T
52 // lambda, const uint32_t cutoff, const T mode, const T speed, const T
53 // weight_decay , const float verbose)
54 // when D = 1
55 // void cut_pursuit(const uint32_t n_nodes, const uint32_t n_edges, const
56 // uint32_t nObs
57 // , std::vector<T> & observation
58 // , const std::vector<uint32_t> & Eu, const std::vector<uint32_t> & Ev
59 // ,const std::vector<T> & edgeWeight, const std::vector<T> &
60 // nodeWeight ,std::vector<T> & solution, const T lambda, const
61 // uint32_t cutoff, const T mode, const T speed
62 // , const float verbose)
63 
64 //-----INPUT-----
65 // 1x1 uint32_t n_nodes = number of nodes
66 // 1x1 uint32_t n_edges = number of edges
67 // 1x1 uint32_t nObs = dimension of data on each node
68 // NxD float observation : the observed signal
69 // Ex1 uint32_t Eu, Ev: the origin and destination of each node
70 // Ex1 float edgeWeight: the edge weight
71 // Nx1 float nodeWeight: the node weight
72 // 1x1 float lambda : the regularization strength
73 // 1x1 uint32_t cutoff : minimal component size
74 // 1x1 float mode : the fidelity function
75 // 0 : linear (for simplex bound data)
76 // 1 : quadratic (default)
77 // 0<a<1: KL with a smoothing (for simplex bound data)
78 // 1x1 float speed : parametrization impacting performance
79 // 0 : slow but precise
80 // 1 : recommended (default)
81 // 2 : fast but approximated (no backward step)
82 // 3 : ludicrous - for prototyping (no backward step)
83 // 1x1 float weight_decay : weight decay to compute the optimal binary partition
84 // 1x1 bool verose : verbosity
85 // 0 : silent
86 // 1 : recommended (default)
87 // 2 : chatty
88 //-----OUTPUT-----
89 // Nx1 float solution: piecewise constant approximation
90 // Nx1 uint32_t inComponent: for each node, in which component it belongs
91 // n_node_redx1 cell components : for each component, list of the nodes
92 // 1x1 n_node_red : number of components
93 // 1x1 uint32_t n_edges_red : number of edges in reduced graph
94 // n_edges_redx1 uint32_t Eu_red, Ev_red : source and target of reduced edges
95 // n_edges_redx1 float edgeWeight_red: weights of reduced edges
96 // n_node_redx1 float nodeWeight_red: weights of reduced nodes
97 //---------------SEGMENTATION--------------------------------------------------
98 // for the segmentation, the functions has a few extra argumens allowing to
99 // record the structrue of the reduced graph
100 // C++ style input
101 // void cut_pursuit(const uint32_t n_nodes, const uint32_t n_edges, const
102 // uint32_t nObs
103 // , std::vector< std::vector<T> > & observation
104 // , const std::vector<uint32_t> & Eu, const std::vector<uint32_t> & Ev
105 // ,const std::vector<T> & edgeWeight, const std::vector<T> &
106 // nodeWeight ,std::vector< std::vector<T> > & solution, , const
107 // std::vector<uint32_t> & in_component
108 // , std::vector< std::vector<uint32_t> > & components
109 // , uint32_t & n_nodes_red, uint32_t & n_edges_red
110 // , std::vector<uint32_t> & Eu_red, std::vector<uint32_t> & Ev_red
111 // , std::vector<T> & edgeWeight_red, std::vector<T> & nodeWeight_red
112 // , const T lambda, const T mode, const T speed, const T weight_decay
113 // , const float verbose)
114 //-----EXTRA INPUT-----
115 // Nx1 uint32_t inComponent: for each node, in which component it belongs
116 // 1x1 n_node_red : number of components
117 // 1x1 uint32_t n_edges_red : number of edges in reduced graph
118 // n_node_redx1 cell components : for each component, list of the nodes
119 // n_edges_redx1 uint32_t Eu_red, Ev_red : source and target of reduced edges
120 // n_edges_redx1 float edgeWeight_red: weights of reduced edges
121 // n_node_redx1 float nodeWeight_red: weights of reduced nodes
122 
123 namespace CP {
124 //===========================================================================
125 //===================== CREATE_CP ===================================
126 //===========================================================================
127 
128 template <typename T>
129 CP::CutPursuit<T>* create_CP(const T mode, const float verbose) {
130  CP::CutPursuit<float>* cp = NULL;
131  fidelityType fidelity = L2;
132  if (mode == 0) {
133  if (verbose > 0) {
134  std::cout << " WITH LINEAR FIDELITY" << std::endl;
135  }
136  fidelity = linear;
137  cp = new CP::CutPursuit_Linear<float>();
138  } else if (mode == 1) {
139  if (verbose > 0) {
140  std::cout << " WITH L2 FIDELITY" << std::endl;
141  }
142  fidelity = L2;
143  cp = new CP::CutPursuit_L2<float>();
144  } else if (mode > 0 && mode < 1) {
145  if (verbose > 0) {
146  std::cout << " WITH KULLBACK-LEIBLER FIDELITY SMOOTHING : " << mode
147  << std::endl;
148  }
149  fidelity = KL;
150  cp = new CP::CutPursuit_KL<float>();
151  cp->parameter.smoothing = mode;
152  } else if (mode == -1) {
153  if (verbose > 0) {
154  std::cout << " WITH ALTERNATE L2 NORM : " << mode << std::endl;
155  }
156  fidelity = SPG;
157  cp = new CP::CutPursuit_KL<float>();
158  cp->parameter.smoothing = mode;
159  } else if (mode == 2) {
160  if (verbose > 0) {
161  std::cout << " PARTITION MODE WITH SPATIAL INFORMATION : " << mode
162  << std::endl;
163  }
164  fidelity = SPG;
165  cp = new CP::CutPursuit_SPG<float>();
166  cp->parameter.smoothing = mode;
167  } else {
168  std::cout << " UNKNOWN MODE, SWICTHING TO L2 FIDELITY" << std::endl;
169  fidelity = L2;
170  cp = new CP::CutPursuit_L2<float>();
171  }
172  cp->parameter.fidelity = fidelity;
173  cp->parameter.verbose = verbose;
174  return cp;
175 }
176 
177 //===========================================================================
178 //===================== cut_pursuit C++-style ============================
179 //===========================================================================
180 template <typename T>
181 void cut_pursuit(const uint32_t n_nodes,
182  const uint32_t n_edges,
183  const uint32_t nObs,
184  std::vector<std::vector<T>>& observation,
185  const std::vector<uint32_t>& Eu,
186  const std::vector<uint32_t>& Ev,
187  const std::vector<T>& edgeWeight,
188  const std::vector<T>& nodeWeight,
189  std::vector<std::vector<T>>& solution,
190  const T lambda,
191  const uint32_t cutoff,
192  const T mode,
193  const T speed,
194  const T weight_decay,
195  const float verbose) { // C-style ++ interface
196  std::srand(1);
197  if (verbose > 0) {
198  std::cout << "L0-CUT PURSUIT";
199  }
200  //--------parameterization---------------------------------------------
201  CP::CutPursuit<T>* cp = create_CP(mode, verbose);
202  set_speed(cp, speed, weight_decay, verbose);
203  set_up_CP(cp, n_nodes, n_edges, nObs, observation, Eu, Ev, edgeWeight,
204  nodeWeight);
205  cp->parameter.reg_strenth = lambda;
206  cp->parameter.cutoff = cutoff;
207  //-------run the optimization------------------------------------------
208  cp->run();
209  //------------write the solution-----------------------------
210  VertexAttributeMap<T> vertex_attribute_map =
211  boost::get(boost::vertex_bundle, cp->main_graph);
212  VertexIterator<T> ite_nod = boost::vertices(cp->main_graph).first;
213  for (uint32_t ind_nod = 0; ind_nod < n_nodes; ind_nod++) {
214  for (uint32_t ind_dim = 0; ind_dim < nObs; ind_dim++) {
215  solution[ind_nod][ind_dim] =
216  vertex_attribute_map[*ite_nod].value[ind_dim];
217  }
218  ite_nod++;
219  }
220  delete cp;
221 }
222 
223 //===========================================================================
224 //===================== cut_pursuit segmentation light C++-style
225 //================
226 //===========================================================================
227 template <typename T>
228 void cut_pursuit(const uint32_t n_nodes,
229  const uint32_t n_edges,
230  const uint32_t nObs,
231  std::vector<std::vector<T>>& observation,
232  const std::vector<uint32_t>& Eu,
233  const std::vector<uint32_t>& Ev,
234  const std::vector<T>& edgeWeight,
235  const std::vector<T>& nodeWeight,
236  std::vector<std::vector<T>>& solution,
237  std::vector<uint32_t>& in_component,
238  std::vector<std::vector<uint32_t>>& components,
239  const T lambda,
240  const uint32_t cutoff,
241  const T mode,
242  const T speed,
243  const T weight_decay,
244  const float verbose) { // C-style ++ interface
245  std::srand(1);
246 
247  if (verbose > 0) {
248  std::cout << "L0-CUT PURSUIT";
249  }
250  //--------parameterization---------------------------------------------
251  CP::CutPursuit<T>* cp = create_CP(mode, verbose);
252 
253  set_speed(cp, speed, weight_decay, verbose);
254  set_up_CP(cp, n_nodes, n_edges, nObs, observation, Eu, Ev, edgeWeight,
255  nodeWeight);
256  cp->parameter.reg_strenth = lambda;
257  cp->parameter.cutoff = cutoff;
258  //-------run the optimization------------------------------------------
259  cp->run();
260  cp->compute_reduced_graph();
261  //------------resize the vectors-----------------------------
262  uint32_t n_nodes_red =
263  static_cast<uint32_t>(boost::num_vertices(cp->reduced_graph));
264  in_component.resize(n_nodes);
265  components.resize(n_nodes_red);
266  //------------write the solution-----------------------------
267  VertexAttributeMap<T> vertex_attribute_map =
268  boost::get(boost::vertex_bundle, cp->main_graph);
269  VertexIterator<T> ite_nod = boost::vertices(cp->main_graph).first;
270  for (uint32_t ind_nod = 0; ind_nod < n_nodes; ind_nod++) {
271  for (uint32_t ind_dim = 0; ind_dim < nObs; ind_dim++) {
272  solution[ind_nod][ind_dim] =
273  vertex_attribute_map[*ite_nod].value[ind_dim];
274  }
275  ite_nod++;
276  }
277  //------------fill the components-----------------------------
278  VertexIndexMap<T> vertex_index_map =
279  get(boost::vertex_index, cp->main_graph);
280  for (uint32_t ind_nod_red = 0; ind_nod_red < n_nodes_red; ind_nod_red++) {
281  size_t component_size = cp->components[ind_nod_red].size();
282  components[ind_nod_red] = std::vector<uint32_t>(component_size, 0);
283  for (size_t ind_nod = 0; ind_nod < component_size; ind_nod++) {
284  components[ind_nod_red][ind_nod] = static_cast<uint32_t>(
285  vertex_index_map(cp->components[ind_nod_red][ind_nod]));
286  }
287  }
288  ite_nod = boost::vertices(cp->main_graph).first;
289  for (uint32_t ind_nod = 0; ind_nod < n_nodes; ind_nod++) {
290  in_component[ind_nod] = vertex_attribute_map[*ite_nod].in_component;
291  ite_nod++;
292  }
293  delete cp;
294 }
295 
296 //===========================================================================
297 //===================== SET_UP_CP C++ style ============================
298 //===========================================================================
299 template <typename T>
300 void set_up_CP(CP::CutPursuit<T>* cp,
301  const uint32_t n_nodes,
302  const uint32_t n_edges,
303  const uint32_t nObs,
304  const std::vector<std::vector<T>> observation,
305  const std::vector<uint32_t> Eu,
306  const std::vector<uint32_t> Ev,
307  const std::vector<T> edgeWeight,
308  const std::vector<T> nodeWeight) {
309  cp->main_graph = Graph<T>(n_nodes);
310  cp->dim = nObs;
311  //--------fill the vertices--------------------------------------------
312  VertexAttributeMap<T> vertex_attribute_map =
313  boost::get(boost::vertex_bundle, cp->main_graph);
314  VertexIterator<T> ite_nod = boost::vertices(cp->main_graph).first;
315  // the node attributes used to fill each node
316  for (uint32_t ind_nod = 0; ind_nod < n_nodes; ind_nod++) {
317  VertexAttribute<T> v_attribute(nObs);
318  for (uint32_t i_dim = 0; i_dim < nObs;
319  i_dim++) { // fill the observation of v_attribute
320  v_attribute.observation[i_dim] = observation[ind_nod][i_dim];
321  } // and its weight
322  v_attribute.weight = nodeWeight[ind_nod];
323  // set the attributes of the current node
324  vertex_attribute_map[*ite_nod++] = v_attribute;
325  }
326  //--------build the edges-----------------------------------------------
327  EdgeAttributeMap<T> edge_attribute_map =
328  boost::get(boost::edge_bundle, cp->main_graph);
329  uint32_t true_ind_edg =
330  0; // this index count the number of edges ACTUALLY added
331  for (uint32_t ind_edg = 0; ind_edg < n_edges;
332  ind_edg++) { // add edges in each direction
333  if (addDoubledge(
334  cp->main_graph, boost::vertex(Eu[ind_edg], cp->main_graph),
335  boost::vertex(Ev[ind_edg], cp->main_graph),
336  edgeWeight[ind_edg], true_ind_edg, edge_attribute_map)) {
337  true_ind_edg += 2;
338  }
339  }
340 }
341 
342 //===========================================================================
343 //===================== SET SPEED ===================================
344 //===========================================================================
345 template <typename T>
346 void set_speed(CP::CutPursuit<T>* cp,
347  const T speed,
348  const T weight_decay,
349  const float verbose) {
350  if (speed == 4) {
351  if (verbose > 0) {
352  std::cout << "PARAMETERIZATION = SPECIAL SUPERPOINTGRAPH"
353  << std::endl;
354  }
355  cp->parameter.flow_steps = 3;
356  cp->parameter.weight_decay = weight_decay;
357  cp->parameter.kmeans_ite = 5;
358  cp->parameter.kmeans_resampling = 10;
359  cp->parameter.max_ite_main = 15;
360  cp->parameter.backward_step = true;
361  cp->parameter.stopping_ratio = 0.05;
362  }
363  if (speed == 3) {
364  if (verbose > 0) {
365  std::cout << "PARAMETERIZATION = LUDICROUS SPEED" << std::endl;
366  }
367  cp->parameter.flow_steps = 1;
368  cp->parameter.weight_decay = weight_decay;
369  cp->parameter.kmeans_ite = 3;
370  cp->parameter.kmeans_resampling = 1;
371  cp->parameter.max_ite_main = 5;
372  cp->parameter.backward_step = false;
373  cp->parameter.stopping_ratio = 0.1;
374  }
375  if (speed == 2) {
376  if (verbose > 0) {
377  std::cout << "PARAMETERIZATION = FAST" << std::endl;
378  }
379  cp->parameter.flow_steps = 2;
380  cp->parameter.weight_decay = weight_decay;
381  cp->parameter.kmeans_ite = 5;
382  cp->parameter.kmeans_resampling = 2;
383  cp->parameter.max_ite_main = 5;
384  cp->parameter.backward_step = true;
385  cp->parameter.stopping_ratio = 0.05;
386  } else if (speed == 0) {
387  if (verbose > 0) {
388  std::cout << "PARAMETERIZATION = SLOW" << std::endl;
389  }
390  cp->parameter.flow_steps = 4;
391  cp->parameter.weight_decay = weight_decay;
392  cp->parameter.kmeans_ite = 8;
393  cp->parameter.kmeans_resampling = 5;
394  cp->parameter.max_ite_main = 20;
395  cp->parameter.backward_step = true;
396  cp->parameter.stopping_ratio = 0.001;
397  } else if (speed == 1) {
398  if (verbose > 0) {
399  std::cout << "PARAMETERIZATION = STANDARD" << std::endl;
400  }
401  cp->parameter.flow_steps = 3;
402  cp->parameter.weight_decay = weight_decay;
403  cp->parameter.kmeans_ite = 5;
404  cp->parameter.kmeans_resampling = 2;
405  cp->parameter.max_ite_main = 10;
406  cp->parameter.backward_step = true;
407  cp->parameter.stopping_ratio = 0.01;
408  }
409 }
410 } // namespace CP
CutPursuit_KL.h
CutPursuit_L2.h
CutPursuit_Linear.h
CutPursuit_SPG.h
vertex_index
long vertex_index
Definition: FileIOFactory.cpp:47
fidelityType
fidelityType
Definition: Common.h:27
L2
@ L2
Definition: Common.h:27
linear
@ linear
Definition: Common.h:27
KL
@ KL
Definition: Common.h:27
SPG
@ SPG
Definition: Common.h:27
NULL
#define NULL
Definition: UVUnwrapping.cpp:40
CP
Definition: API.h:123
CP::addDoubledge
bool addDoubledge(Graph< T > &g, const VertexDescriptor< T > &source, const VertexDescriptor< T > &target, const T weight, uint32_t eIndex, EdgeAttributeMap< T > &edge_attribute_map, bool real=true)
Definition: Graph.h:110
CP::set_speed
void set_speed(CP::CutPursuit< T > *cp, const T speed, const T weight_decay, const float verbose)
Definition: API.h:346
CP::VertexAttributeMap
boost::vec_adj_list_vertex_property_map< Graph< T >, Graph< T > *, VertexAttribute< T >, VertexAttribute< T > &, boost::vertex_bundle_t > VertexAttributeMap
Definition: Graph.h:92
CP::VertexIterator
typename boost::graph_traits< Graph< T > >::vertex_iterator VertexIterator
Definition: Graph.h:82
CP::cut_pursuit
void cut_pursuit(const uint32_t n_nodes, const uint32_t n_edges, const uint32_t nObs, std::vector< std::vector< T >> &observation, const std::vector< uint32_t > &Eu, const std::vector< uint32_t > &Ev, const std::vector< T > &edgeWeight, const std::vector< T > &nodeWeight, std::vector< std::vector< T >> &solution, const T lambda, const uint32_t cutoff, const T mode, const T speed, const T weight_decay, const float verbose)
Definition: API.h:181
CP::set_up_CP
void set_up_CP(CP::CutPursuit< T > *cp, const uint32_t n_nodes, const uint32_t n_edges, const uint32_t nObs, const std::vector< std::vector< T >> observation, const std::vector< uint32_t > Eu, const std::vector< uint32_t > Ev, const std::vector< T > edgeWeight, const std::vector< T > nodeWeight)
Definition: API.h:300
CP::Graph
typename boost::adjacency_list< boost::vecS, boost::vecS, boost::directedS, VertexAttribute< T >, EdgeAttribute< T > > Graph
Definition: Graph.h:71
CP::VertexIndexMap
typename boost::property_map< Graph< T >, boost::vertex_index_t >::type VertexIndexMap
Definition: Graph.h:103
CP::EdgeAttributeMap
boost::adj_list_edge_property_map< boost::directed_tag, EdgeAttribute< T >, EdgeAttribute< T > &, uint64_t, CP::EdgeAttribute< T >, boost::edge_bundle_t > EdgeAttributeMap
Definition: Graph.h:100
CP::create_CP
CP::CutPursuit< T > * create_CP(const T mode, const float verbose)
Definition: API.h:129
QtCompat::endl
QTextStream & endl(QTextStream &stream)
Definition: QtCompat.h:718
CP::CutPursuit_KL
Definition: CutPursuit_KL.h:14
CP::CutPursuit_L2
Definition: CutPursuit_L2.h:14
CP::CutPursuit_Linear
Definition: CutPursuit_Linear.h:14
CP::CutPursuit_SPG
Definition: CutPursuit_SPG.h:14
CP::CutPursuit::parameter
CPparameter< T > parameter
Definition: CutPursuit.h:65
CP::CutPursuit::main_graph
Graph< T > main_graph
Definition: CutPursuit.h:45
CP::CutPursuit::components
std::vector< std::vector< VertexDescriptor< T > > > components
Definition: CutPursuit.h:49
CP::CutPursuit::compute_reduced_graph
void compute_reduced_graph()
Definition: CutPursuit.h:541
CP::CutPursuit::dim
uint32_t dim
Definition: CutPursuit.h:59
CP::CutPursuit::run
std::pair< std::vector< T >, std::vector< T > > run()
Definition: CutPursuit.h:94
CP::CutPursuit::reduced_graph
Graph< T > reduced_graph
Definition: CutPursuit.h:46
CP::VertexAttribute::observation
std::vector< T > observation
Definition: Graph.h:37