SMIL  1.0.4
Dendrogram.hpp
1 /*
2  * Copyright (c) 2011-2015, Matthieu FAESSEL and ARMINES
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
14  * names of its contributors may be used to endorse or promote products
15  * derived from this software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
18  * AND ANY
19  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
20  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE
22  * LIABLE FOR ANY
23  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
24  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
25  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
26  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29  *
30  * last modification by Amin Fehri
31  *
32  */
33 
34 #ifndef _DENDROGRAM_HPP
35 #define _DENDROGRAM_HPP
36 
37 #include "DendroNode.hpp"
38 
39 using namespace std;
40 
41 namespace smil
42 {
46  template <class MarkerLabelT = size_t, class NodeT = size_t,
47  class ValGraphT = size_t>
48  class Dendrogram
49  {
50  protected:
51  size_t nbrNodes;
52  size_t nbrMarkers;
53  std::vector<DendroNode<MarkerLabelT> *> dendroNodes;
54 
55  public:
57 
59  typedef std::map<NodeT, ValGraphT> NodeValuesType;
60  typedef std::vector<Edge<NodeT, ValGraphT>> EdgeListType;
61 
63  Dendrogram() : nbrNodes(0), nbrMarkers(0), dendroNodes(0){};
65  Dendrogram(const Dendrogram &dendrogramToCopy)
66  {
67  nbrNodes = dendrogramToCopy
68  .nbrNodes; // N leafs and (N-1) internal nodes = (2N-1)
69  nbrMarkers = dendrogramToCopy.nbrMarkers;
70  std::vector<DendroNodeType *> dendroNodesToCopy =
71  dendrogramToCopy.dendroNodes;
72  for (size_t i = 0; i < dendrogramToCopy.dendroNodes.size(); i++) {
73  dendroNodes.push_back(new DendroNodeType);
74  }
75  for (size_t i = 0; i < nbrNodes; i++) {
76  DendroNodeType &curNode = *dendroNodesToCopy[i];
77  DendroNodeType &newNode = *dendroNodes[i];
78 
79  newNode.setLabel(curNode.getLabel());
80  newNode.setMarker(curNode.getMarker());
81  newNode.setNbMarkersUnder(curNode.getNbMarkersUnder());
82  newNode.setMarkersCount(curNode.getMarkersCount());
83  newNode.setValuation(curNode.getValuation());
84  newNode.setLeafValuation(curNode.getLeafValuation());
85  newNode.setInternalNodeValuationInitial(
87  newNode.setInternalNodeValuationFinal(
88  curNode.getInternalNodeValuationFinal());
89  newNode.setIsInternalNode(curNode.getIsInternalNode());
90  newNode.setLookupProgeny(curNode.getLookupProgeny());
91  newNode.setMoments(curNode.getMoments());
92  newNode.setContoursCount(curNode.getContoursCount());
93  newNode.setContoursSize(curNode.getContoursSize());
94  }
95  for (size_t i = 0; i < nbrNodes; i++) {
96  DendroNodeType &newNode = *dendroNodes[i];
97  MarkerLabelT researchedLabel = newNode.getLabel();
98 
99  DendroNodeType &correspondingNode =
100  *dendrogramToCopy.researchLabel(researchedLabel);
101  MarkerLabelT fatherLabel = correspondingNode.getFather()->getLabel();
102  newNode.setFather(this->researchLabel(fatherLabel));
103 
104  if (correspondingNode.getIsInternalNode() == true) {
105  MarkerLabelT childLeftLabel =
106  correspondingNode.getChildLeft()->getLabel();
107  MarkerLabelT childRightLabel =
108  correspondingNode.getChildRight()->getLabel();
109  MarkerLabelT neighborLeftLabel =
110  correspondingNode.getNeighborLeft()->getLabel();
111  MarkerLabelT neighborRightLabel =
112  correspondingNode.getNeighborRight()->getLabel();
113  newNode.setChildLeft(this->researchLabel(childLeftLabel));
114  newNode.setChildRight(this->researchLabel(childRightLabel));
115  newNode.setNeighborLeft(this->researchLabel(neighborLeftLabel));
116  newNode.setNeighborRight(this->researchLabel(neighborRightLabel));
117  }
118  }
119  };
120 
127  : nbrNodes(0), nbrMarkers(0), dendroNodes(0)
128  {
129  mst.sortEdges(
130  true); // sort Edges of the MST by increasing weights of edges
131  // Extract informations from the MST
132  size_t leavesNbr = mst.getNodeNbr();
133  size_t internalNodesNbr = mst.getEdgeNbr();
134  NodeValuesType &mstNodes = mst.getNodeValues();
135  EdgeListType & mstEdges = mst.getEdges();
136 
137  // Set the number of required nodes and creates them in dendroNodes
138  nbrNodes = leavesNbr + internalNodesNbr;
139  // nbrNodes = 2*leavesNbr - 1 ;
140  for (size_t i = 0; i < nbrNodes; i++) {
141  DendroNodeType *nNode = new DendroNodeType;
142  dendroNodes.push_back(nNode);
143  }
144  // Filling the leaves
145  for (size_t i = 0; i < leavesNbr; i++) {
146  // dendroNodes is filled with leavesNbr
147  // leaves and then with internalNodesNbr
148  // internal nodes
149  DendroNodeType &curNode = *dendroNodes[i];
150  // so Leaves have labels from 0 to leavesNbr-1
151  curNode.setLabel(mstNodes.find(i)->first);
152  curNode.setLeafValuation(mstNodes.find(i)->second);
153  curNode.setValuation(mstNodes.find(i)->second);
154  // initialize the lookupProgeny for each leaf with 1 for the label of
155  // the leaf, and 0 elsewhere
156  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
157  nLookupProgeny.at(mstNodes.find(i)->first) = 1;
158  curNode.setLookupProgeny(nLookupProgeny);
159  // initialize the contoursCount
160  std::vector<MarkerLabelT> initContoursCount(leavesNbr, 0);
161  curNode.setContoursCount(initContoursCount);
162  }
163 
164  // Filling the internal nodes
165  for (size_t i = 0; i < internalNodesNbr; i++) {
166  // dendroFromMSTNodes is filled with leavesNbr leaves and then
167  // with internalNodesNbr internal nodes
168  DendroNodeType &curNode = *dendroNodes[leavesNbr + i];
169  curNode.setIsInternalNode(1);
170  curNode.setLabel(leavesNbr + i);
171  curNode.setInternalNodeValuationInitial(mstEdges[i].weight);
172  MarkerLabelT NeighborLeftLabel =
173  min(mstEdges[i].source, mstEdges[i].target);
174  MarkerLabelT NeighborRightLabel =
175  max(mstEdges[i].source, mstEdges[i].target);
176 
177  // dendroNodesLabels[NeighborLeftLabel] to modify
178  curNode.setNeighborLeft(dendroNodes[NeighborLeftLabel]);
179  // dendroNodesLabels[NeighborRightLabel] to modify
180  curNode.setNeighborRight(dendroNodes[NeighborRightLabel]);
181 
182  curNode.setChildLeft(curNode.getNeighborLeft()->getAncestor());
183  curNode.setChildRight(curNode.getNeighborRight()->getAncestor());
184 
185  curNode.getChildLeft()->setFather(&curNode);
186  curNode.getChildRight()->setFather(&curNode);
187 
188  std::vector<MarkerLabelT> lookupChildLeft =
189  curNode.getChildLeft()->getLookupProgeny();
190  std::vector<MarkerLabelT> lookupChildRight =
191  curNode.getChildRight()->getLookupProgeny();
192  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
193  for (MarkerLabelT i = 0; i < lookupChildLeft.size(); i++) {
194  nLookupProgeny.at(i) =
195  max(lookupChildRight.at(i), lookupChildLeft.at(i));
196  }
197  // new lookupProgeny for each leaf, with 1 for the label of
198  // the leaf, and 0 elsewhere
199  curNode.setLookupProgeny(nLookupProgeny);
200  }
201  // last node parameters
202  DendroNodeType &lastNode = *dendroNodes[leavesNbr + internalNodesNbr - 1];
203  lastNode.setFather(&lastNode); // the last internal node is its own father
204 
205  }; // end Dendrogram(mst)
206 
208 
218  smil::Image<MarkerLabelT> &imLabels,
219  smil::Image<MarkerLabelT> &imIn, const size_t nbOfMoments = 5)
220  : nbrNodes(0), nbrMarkers(0), dendroNodes(0)
221  {
222  // sort Edges of the MST by increasing weights of edges
223  mst.sortEdges(true);
224  // Extract informations from the MST
225  size_t leavesNbr = mst.getNodeNbr();
226  size_t internalNodesNbr = mst.getEdgeNbr();
227  NodeValuesType &mstNodes = mst.getNodeValues();
228  EdgeListType & mstEdges = mst.getEdges();
229 
230  // Set the number of required nodes and creates them in dendroNodes
231  nbrNodes = leavesNbr + internalNodesNbr;
232  // nbrNodes = 2*leavesNbr - 1 ;
233  for (size_t i = 0; i < nbrNodes; i++) {
234  DendroNodeType *nNode = new DendroNodeType;
235  dendroNodes.push_back(nNode);
236  }
237  // Filling the leaves
238  for (size_t i = 0; i < leavesNbr; i++) {
239  // dendroNodes is filled with leavesNbr leaves
240  // and then with internalNodesNbr internal nodes
241  DendroNodeType &curNode = *dendroNodes[i];
242  // Leaves have labels from 0 to leavesNbr-1
243  curNode.setLabel(mstNodes.find(i)->first);
244  curNode.setLeafValuation(mstNodes.find(i)->second);
245  curNode.setValuation(mstNodes.find(i)->second);
246  // initialize the moments
247  std::vector<double> initMoments(nbOfMoments);
248  curNode.setMoments(initMoments);
249  // initialize the lookupProgeny for each leaf
250  // with 1 for the label of the leaf, and 0 elsewhere
251  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
252  nLookupProgeny.at(mstNodes.find(i)->first) = 1;
253  curNode.setLookupProgeny(nLookupProgeny);
254  // initialize the contoursCount
255  std::vector<MarkerLabelT> initContoursCount(leavesNbr, 0);
256  curNode.setContoursCount(initContoursCount);
257  }
258 
259  // Filling the internal nodes
260  for (size_t i = 0; i < internalNodesNbr; i++) {
261  // dendroFromMSTNodes is filled with leavesNbr leaves and then
262  // with internalNodesNbr internal nodes
263  DendroNodeType &curNode = *dendroNodes[leavesNbr + i];
264  curNode.setIsInternalNode(1);
265  curNode.setLabel(leavesNbr + i);
266  curNode.setInternalNodeValuationInitial(mstEdges[i].weight);
267  MarkerLabelT NeighborLeftLabel =
268  min(mstEdges[i].source, mstEdges[i].target);
269  MarkerLabelT NeighborRightLabel =
270  max(mstEdges[i].source, mstEdges[i].target);
271 
272  // dendroNodesLabels[NeighborLeftLabel] to modify
273  curNode.setNeighborLeft(dendroNodes[NeighborLeftLabel]);
274  // dendroNodesLabels[NeighborRightLabel] to modify
275  curNode.setNeighborRight(dendroNodes[NeighborRightLabel]);
276 
277  curNode.setChildLeft(curNode.getNeighborLeft()->getAncestor());
278  curNode.setChildRight(curNode.getNeighborRight()->getAncestor());
279 
280  curNode.getChildLeft()->setFather(&curNode);
281  curNode.getChildRight()->setFather(&curNode);
282 
283  std::vector<MarkerLabelT> lookupChildLeft =
284  curNode.getChildLeft()->getLookupProgeny();
285  std::vector<MarkerLabelT> lookupChildRight =
286  curNode.getChildRight()->getLookupProgeny();
287  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
288  for (MarkerLabelT i = 0; i < lookupChildLeft.size(); i++) {
289  nLookupProgeny.at(i) =
290  max(lookupChildRight.at(i), lookupChildLeft.at(i));
291  }
292  // new lookupProgeny for each leaf, with 1 for the label
293  // of the leaf, and 0 elsewhere
294  curNode.setLookupProgeny(nLookupProgeny);
295  }
296  // last node parameters
297  DendroNodeType &lastNode = *dendroNodes[leavesNbr + internalNodesNbr - 1];
298  // the last internal node is its own father
299  lastNode.setFather(&lastNode);
300 
301  // Complete the moments values
302  this->setMomentsContours(imLabels, imIn, nbOfMoments);
303  }; // end Dendrogram(mst,imLabels,imMarkers)
304 
306  virtual ~Dendrogram()
307  {
308  for (size_t i = 0; i < nbrNodes; i++) {
309  delete dendroNodes[i];
310  }
311  };
314  {
315  return Dendrogram(*this);
316  };
318  void sortNodes(bool reverse = false)
319  {
320  if (!reverse) // ... by growing valuationInitial
321  std::stable_sort(dendroNodes.begin(), dendroNodes.end(),
322  DendroNodeType::isInferior);
323  else // ... by decreasing valuationInitial
324  std::stable_sort(dendroNodes.begin(), dendroNodes.end(),
325  DendroNodeType::isSuperior);
326  };
329  {
330  std::stable_sort(dendroNodes.begin(), dendroNodes.end(),
331  DendroNodeType::isSuperior);
332  };
333  void addDendroNodes(DendroNodeType *dendroNode)
334  {
335  dendroNodes.push_back(dendroNode);
336  };
337  DendroNodeType *researchLabel(MarkerLabelT researchedLabel) const
338  {
339  DendroNodeType *dendronodeToReturn(0);
340  for (size_t i = 0; i < nbrNodes; i++) {
341  DendroNodeType &curNode = *dendroNodes[i];
342  if (curNode.getLabel() == researchedLabel) {
343  dendronodeToReturn = dendroNodes[i];
344  }
345  }
346  return dendronodeToReturn;
347  };
348  // TODO : check -1
349  // void setMarkers(Image<MarkerLabelT> &imLabels,
350  // Image<MarkerLabelT> &imMarkers) {
351  // // get the Image makers regions, with the marker value instead of the
352  // label Image<MarkerLabelT> imMarkersbis =
353  // Image<MarkerLabelT>(imMarkers); std::map<MarkerLabelT, MarkerLabelT>
354  // lookupMapNewMarkers; // to get the markers 1,2,3,...
355  //
356  // int sizeImMarkers[3];
357  // imMarkers.getSize(sizeImMarkers);
358  // // imMarkers traversal to fill the lookupMapNewMarkers
359  // MarkerLabelT nMarker = 1;
360  // lookupMapNewMarkers[0] = 0;
361  // for (int i = 0; i < sizeImMarkers[0]; i++) {
362  // for (int j = 0; j < sizeImMarkers[1]; j++) {
363  // MarkerLabelT markerValue = imMarkers.getPixel(i, j);
364  // if (markerValue != 0 && (lookupMapNewMarkers.find(markerValue) ==
365  // lookupMapNewMarkers.end())) {
366  // lookupMapNewMarkers[markerValue] = nMarker;
367  // nMarker++;
368  // }
369  // }
370  // }
371  // applyLookup(
372  // imMarkers, lookupMapNewMarkers,
373  // imMarkersbis); // imMarkersbis = new imMarkers with values
374  // 1,2,3,...
375  // int nbMarkers = maxVal(imMarkersbis);
376  // // cout << "minVal : " << minVal(imMarkersbis,true) << "
377  // maxVal "
378  // // << maxVal(imMarkersbis) << endl;
379  // // imMarkersbis.show();
380  // // Gui::execLoop();
381  // this->setNbrMarkers(nbMarkers);
382  //
383  // // initialize the markersCount
384  // for (size_t i = 0; i < nbrNodes; i++) {
385  // DendroNodeType &curNode = *dendroNodes[i];
386  // std::vector<MarkerLabelT> nMarkersCount(nbrMarkers, 0);
387  // curNode.setMarkersCount(nMarkersCount);
388  // }
389  //
390  // std::map<MarkerLabelT, MarkerLabelT> lookupMapMarkers; // to get the
391  // marker
392  // // value if
393  // it's the
394  // // only one
395  // for a
396  // // given label
397  // //
398  // (label,marker)
399  // std::map<MarkerLabelT, std::vector<MarkerLabelT> >
400  // lookupMapMarkersCount; //(label,markerCount)
401  // int sizeIm[3];
402  // imLabels.getSize(sizeIm);
403  //
404  // // imLabels traversal to initialize the lookupMapMarkers
405  // for (int i = 0; i < sizeIm[0]; i++) {
406  // for (int j = 0; j < sizeIm[1]; j++) {
407  // MarkerLabelT labelValue = imLabels.getPixel(i, j);
408  // lookupMapMarkers[labelValue] = 0;
409  // lookupMapMarkersCount[labelValue] =
410  // std::vector<MarkerLabelT>(nbrMarkers, 0);
411  // }
412  // }
413  //
414  // // imMarkersbis traversal
415  // for (int i = 0; i < sizeIm[0]; i++) {
416  // for (int j = 0; j < sizeIm[1]; j++) {
417  // MarkerLabelT markerValue = imMarkersbis.getPixel(i, j);
418  // if (markerValue != 0) {
419  // MarkerLabelT labelValue = imLabels.getPixel(i, j);
420  // lookupMapMarkersCount[labelValue].at(markerValue - 1) =
421  // lookupMapMarkersCount[labelValue].at(markerValue - 1) + 1;
422  // if (lookupMapMarkers[labelValue] != 0 &&
423  // lookupMapMarkers[labelValue] != -1 &&
424  // lookupMapMarkers[labelValue] != markerValue) {
425  // lookupMapMarkers[labelValue] = -1;
426  // } else if (lookupMapMarkers[labelValue] == 0) {
427  // lookupMapMarkers[labelValue] = markerValue;
428  // }
429  // }
430  // }
431  // } // end imMarkers traversal
432  // // lookupMapMarkers traversal
433  // for (typename std::map<MarkerLabelT, MarkerLabelT>::iterator iter =
434  // lookupMapMarkers.begin();
435  // iter != lookupMapMarkers.end(); ++iter) {
436  // if (iter->second == -1) {
437  // iter->second = 0;
438  // }
439  // }
440  // for (typename std::map<MarkerLabelT, MarkerLabelT>::iterator iter =
441  // lookupMapMarkers.begin();
442  // iter != lookupMapMarkers.end(); ++iter) {
443  // if (iter->second != 0) {
444  // MarkerLabelT researchedLabel = iter->first;
445  // MarkerLabelT markerValue = iter->second;
446  // // cout << "researchedLabel : " << researchedLabel << " ,
447  // // markerValue : " << markerValue << endl;
448  // this->researchLabel(researchedLabel)->setMarker(markerValue);
449  // }
450  // }
451  // // lookupMapMarkersCount traversal
452  // for (typename std::map<MarkerLabelT, std::vector<MarkerLabelT>
453  // >::iterator
454  // iter = lookupMapMarkersCount.begin();
455  // iter != lookupMapMarkersCount.end(); ++iter) {
456  // // cout << "researchedLabel : " << iter->first << " ,
457  // // markersValues : " << iter->second[25] << endl;
458  // // cout << "label : " << iter->first << endl;
459  // // for (int i = 0;i<iter->second.size();i++){
460  // // if (iter->second[i]!=0){
461  // // cout << " marker = " << i << " ; nombre = " <<
462  // // iter->second[i] << endl;
463  // // }
464  // // }
465  // // cout << std::accumulate((iter->second).begin(),
466  // // (iter->second).end(), 0) << endl;
467  // this->researchLabel(iter->first)->setMarkersCount(iter->second);
468  // this->researchLabel(iter->first)
469  // ->setNbMarkersUnder(
470  // std::accumulate((iter->second).begin(),
471  // (iter->second).end(), 0));
472  // }
473  //
474  // this->reorganize(); // to get the new nbMarkersUnder and markersCount
475  // // applyLookup(imLabels,lookupMapMarkers,imMarkersbis);
476  // // imMarkersbis.show();
477  // // Gui::execLoop();
478  // }
479 
480  void setMomentsContours(Image<MarkerLabelT> &imLabels,
481  Image<MarkerLabelT> &imIn, const size_t nbOfMoments)
482  {
483  size_t leavesNbr = nbrNodes / 2 + 1;
484 
485  std::map<MarkerLabelT, std::vector<double>> lookupMapMoments;
486  std::map<MarkerLabelT, std::vector<MarkerLabelT>> lookupMapContours;
487  size_t sizeIm[3];
488  imLabels.getSize(sizeIm);
489 
490  // imLabels traversal to initialize the lookupMapMoments and
491  // lookupMapContours
492  for (size_t i = 0; i < sizeIm[0]; i++) {
493  for (size_t j = 0; j < sizeIm[1]; j++) {
494  MarkerLabelT labelValue = imLabels.getPixel(i, j);
495  if (lookupMapMoments.count(labelValue) == 0) {
496  std::vector<double> initMoments(nbOfMoments);
497  lookupMapMoments[labelValue] = initMoments;
498  std::vector<MarkerLabelT> initContours(leavesNbr, 0);
499  lookupMapContours[labelValue] = initContours;
500  }
501  }
502  }
503 
504  // imIn/imLabels traversal to complete lookupMapMoments
505  for (size_t i = 0; i < sizeIm[0]; i++) {
506  for (size_t j = 0; j < sizeIm[1]; j++) {
507  MarkerLabelT imValue = imIn.getPixel(i, j);
508  MarkerLabelT labelValue = imLabels.getPixel(i, j);
509  for (size_t k = 0; k < nbOfMoments; k++) {
510  lookupMapMoments[labelValue][k] =
511  lookupMapMoments[labelValue][k] + std::pow(imValue, k);
512  }
513  }
514  } // end imIn traversal
515 
516  // imLabels traversal to complete lookupMapContours
517  // we avoid borders
518  for (size_t i = 1; i < sizeIm[0] - 1; i++) {
519  for (size_t j = 1; j < sizeIm[1] - 1; j++) {
520  MarkerLabelT labelValue = imLabels.getPixel(i, j);
521  if (imLabels.getPixel(i - 1, j - 1) != labelValue) {
522  MarkerLabelT newLabel = imLabels.getPixel(i - 1, j - 1);
523  lookupMapContours[labelValue].at(newLabel) =
524  lookupMapContours[labelValue].at(newLabel) + 1;
525  } else if (imLabels.getPixel(i - 1, j) != labelValue) {
526  MarkerLabelT newLabel = imLabels.getPixel(i - 1, j);
527  lookupMapContours[labelValue].at(newLabel) =
528  lookupMapContours[labelValue].at(newLabel) + 1;
529  } else if (imLabels.getPixel(i - 1, j + 1) != labelValue) {
530  MarkerLabelT newLabel = imLabels.getPixel(i - 1, j + 1);
531  lookupMapContours[labelValue].at(newLabel) =
532  lookupMapContours[labelValue].at(newLabel) + 1;
533  } else if (imLabels.getPixel(i, j - 1) != labelValue) {
534  MarkerLabelT newLabel = imLabels.getPixel(i, j - 1);
535  lookupMapContours[labelValue].at(newLabel) =
536  lookupMapContours[labelValue].at(newLabel) + 1;
537  } else if (imLabels.getPixel(i, j) != labelValue) {
538  MarkerLabelT newLabel = imLabels.getPixel(i, j);
539  lookupMapContours[labelValue].at(newLabel) =
540  lookupMapContours[labelValue].at(newLabel) + 1;
541  } else if (imLabels.getPixel(i, j + 1) != labelValue) {
542  MarkerLabelT newLabel = imLabels.getPixel(i, j + 1);
543  lookupMapContours[labelValue].at(newLabel) =
544  lookupMapContours[labelValue].at(newLabel) + 1;
545  } else if (imLabels.getPixel(i + 1, j - 1) != labelValue) {
546  MarkerLabelT newLabel = imLabels.getPixel(i + 1, j - 1);
547  lookupMapContours[labelValue].at(newLabel) =
548  lookupMapContours[labelValue].at(newLabel) + 1;
549  } else if (imLabels.getPixel(i + 1, j) != labelValue) {
550  MarkerLabelT newLabel = imLabels.getPixel(i + 1, j);
551  lookupMapContours[labelValue].at(newLabel) =
552  lookupMapContours[labelValue].at(newLabel) + 1;
553  } else if (imLabels.getPixel(i + 1, j + 1) != labelValue) {
554  MarkerLabelT newLabel = imLabels.getPixel(i + 1, j + 1);
555  lookupMapContours[labelValue].at(newLabel) =
556  lookupMapContours[labelValue].at(newLabel) + 1;
557  }
558  }
559  } // end imLabels traversal
560 
561  // lookupMapMoments traversal
562  for (typename std::map<MarkerLabelT, std::vector<double>>::iterator iter =
563  lookupMapMoments.begin();
564  iter != lookupMapMoments.end(); ++iter) {
565  MarkerLabelT researchedLabel = iter->first;
566  std::vector<double> momentsValues = iter->second;
567  // cout << "label = " << researchedLabel << endl;
568  // cout << "moments " << endl;
569  // cout << "0 : " << momentsValues[0] << endl;
570  // cout << "1 : " << momentsValues[1] << endl;
571  // cout << "2 : " << momentsValues[2] << endl;
572  // cout << "3 : " << momentsValues[3] << endl;
573  // cout << "4 : " << momentsValues[4] << endl;
574  researchLabel(researchedLabel)->setMoments(momentsValues);
575  }
576 
577  // lookupMapContours traversal
578  for (typename std::map<MarkerLabelT, std::vector<MarkerLabelT>>::iterator
579  iter = lookupMapContours.begin();
580  iter != lookupMapContours.end(); ++iter) {
581  MarkerLabelT researchedLabel = iter->first;
582  std::vector<MarkerLabelT> contoursCountValues = iter->second;
583  // cout << "label = " << researchedLabel << endl;
584  // cout << "contoursCountValues " << endl;
585  // cout << "0 : " << contoursCountValues[0] << endl;
586  // cout << "1 : " << contoursCountValues[1] << endl;
587  // cout << "2 : " << contoursCountValues[2] << endl;
588  // cout << "3 : " << contoursCountValues[3] << endl;
589  // cout << "4 : " << contoursCountValues[4] << endl;
590  researchLabel(researchedLabel)->setContoursCount(contoursCountValues);
591  }
592  this->reorganize();
593  };
594 
595  void reorganize()
596  {
597  size_t leavesNbr = nbrNodes / 2 + 1;
598  size_t internalNodesNbr = nbrNodes / 2;
599  // Clean the kinship between nodes
600  for (size_t i = 0; i < nbrNodes; i++) {
601  DendroNodeType &curNode = *dendroNodes[i];
602  curNode.setChildLeft(0);
603  curNode.setChildRight(0);
604  curNode.setFather(0);
605  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
606  curNode.setLookupProgeny(nLookupProgeny);
607 
608  if (curNode.getIsInternalNode()) {
609  std::vector<MarkerLabelT> nMarkersCount(nbrMarkers, 0);
610  curNode.setMarkersCount(nMarkersCount);
611  }
612  }
613  // Filling the leaves
614  for (size_t i = 0; i < nbrNodes; i++) {
615  DendroNodeType &curNode = *dendroNodes[i];
616  // initialize the lookupProgeny
617  if (!curNode.getIsInternalNode()) {
618  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
619  nLookupProgeny.at(curNode.getLabel()) = 1;
620  // new lookupProgeny for each leaf, with 1 for the label
621  // of the leaf, and 0 elsewhere
622  curNode.setLookupProgeny(nLookupProgeny);
623  }
624  }
625  // we sort nodes by growing internalNodeValuationInitial
626  this->sortNodes();
627  // Filling the internal nodes
628  for (size_t i = 0; i < nbrNodes; i++) {
629  DendroNodeType &curNode = *dendroNodes[i];
630  if (curNode.getIsInternalNode()) {
631  curNode.setChildLeft(curNode.getNeighborLeft()->getAncestor());
632  curNode.setChildRight(curNode.getNeighborRight()->getAncestor());
633 
634  curNode.getChildLeft()->setFather(&curNode);
635  curNode.getChildRight()->setFather(&curNode);
636 
637  // Computing the new lookupProgeny
638  std::vector<MarkerLabelT> lookupChildLeft =
639  curNode.getChildLeft()->getLookupProgeny();
640  std::vector<MarkerLabelT> lookupChildRight =
641  curNode.getChildRight()->getLookupProgeny();
642  std::vector<MarkerLabelT> nLookupProgeny(leavesNbr, 0);
643  for (size_t i = 0; i < lookupChildLeft.size(); i++) {
644  nLookupProgeny.at(i) =
645  max(lookupChildRight.at(i), lookupChildLeft.at(i));
646  }
647  // new lookupProgeny for each leaf, with 1 for the label
648  // of the leaf, and 0 elsewhere
649  curNode.setLookupProgeny(nLookupProgeny);
650 
651  // Computing the new nbMarkersUnder
652  MarkerLabelT nbMarkersUnderLeft =
653  curNode.getChildLeft()->getNbMarkersUnder();
654  MarkerLabelT nbMarkersUnderRight =
655  curNode.getChildRight()->getNbMarkersUnder();
656  curNode.setNbMarkersUnder(nbMarkersUnderLeft + nbMarkersUnderRight);
657 
658  // Computing the new markersCount
659  std::vector<MarkerLabelT> markersCountLeft =
660  curNode.getChildLeft()->getMarkersCount();
661  std::vector<MarkerLabelT> markersCountRight =
662  curNode.getChildRight()->getMarkersCount();
663  std::vector<MarkerLabelT> nMarkersCount(nbrMarkers, 0);
664  for (size_t i = 0; i < markersCountLeft.size(); i++) {
665  nMarkersCount.at(i) =
666  markersCountLeft.at(i) + markersCountRight.at(i);
667  }
668  curNode.setMarkersCount(nMarkersCount);
669 
670  // Computing the new moments vector
671  std::vector<double> momentsChildLeft =
672  curNode.getChildLeft()->getMoments();
673  std::vector<double> momentsChildRight =
674  curNode.getChildRight()->getMoments();
675  std::vector<double> nMoments(momentsChildLeft.size(), 0);
676  for (size_t i = 0; i < momentsChildLeft.size(); i++) {
677  nMoments.at(i) = momentsChildLeft.at(i) + momentsChildRight.at(i);
678  }
679  // new moments vector for each node
680  curNode.setMoments(nMoments);
681 
682  // Computing the new fidelity term
683 
684  // Computing the new contoursCount
685  std::vector<MarkerLabelT> contoursCountLeft =
686  curNode.getChildLeft()->getContoursCount();
687  std::vector<MarkerLabelT> contoursCountRight =
688  curNode.getChildRight()->getContoursCount();
689 
690  std::vector<MarkerLabelT> nContoursCount(leavesNbr, 0);
691  for (size_t i = 0; i < contoursCountLeft.size(); i++) {
692  nContoursCount.at(i) =
693  contoursCountLeft.at(i) + contoursCountRight.at(i);
694  }
695  // new lookupProgeny for each leaf, with 1 for the label
696  // of the leaf, and 0 elsewhere
697  curNode.setContoursCount(nContoursCount);
698  // Computing the new contoursSize
699  double nContoursSize = 0;
700  for (size_t i = 0; i < contoursCountLeft.size(); i++) {
701  nContoursSize = nContoursSize +
702  nContoursCount.at(i) * (1 - nLookupProgeny.at(i));
703  }
704  curNode.setContoursSize(nContoursSize);
705  }
706  // set the contours size for leaves too
707  else if (!curNode.getIsInternalNode()) {
708  std::vector<MarkerLabelT> lookupProgeny = curNode.getLookupProgeny();
709  std::vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
710 
711  double nContoursSize = 0;
712  for (size_t i = 0; i < contoursCount.size(); i++) {
713  nContoursSize =
714  nContoursSize + contoursCount.at(i) * (1 - lookupProgeny.at(i));
715  }
716  curNode.setContoursSize(nContoursSize);
717  }
718  }
719  // last node parameters
720  DendroNodeType &lastNode = *dendroNodes[leavesNbr + internalNodesNbr - 1];
721  lastNode.setFather(&lastNode); // the last internal node is its own father
722  };
723 
731  void normalize(const std::string typeOfNormalization = "reg")
732  {
733  // size_t leavesNbr = nbrNodes/2+1;
734  size_t internalNodesNbr = nbrNodes / 2;
735 
736  if (typeOfNormalization == "reg") {
737  // we sort nodes by growing internalNodeValuationInitial
738  this->sortNodes();
739  double counter = 1;
740  // Filling the internal nodes
741  for (size_t i = 0; i < nbrNodes; i++) {
742  DendroNodeType &curNode = *dendroNodes.at(i);
743  if (curNode.getIsInternalNode()) {
744  double newVal = counter / internalNodesNbr;
745  curNode.setInternalNodeValuationFinal(newVal);
746  counter = counter + 1;
747  }
748  }
749  this->putValuationsFinalInInitial();
750  this->reorganize();
751  } else if (typeOfNormalization == "maxnorm") {
752  // we sort nodes by decreasing internalNodeValuationInitial
753  this->sortNodes(true);
754  double maxVal;
755 
756  // get the biggest finite value
757  for (size_t i = 0; i < nbrNodes; i++) {
758  DendroNodeType &curNodeInit = *dendroNodes.at(i);
759  maxVal = curNodeInit.getInternalNodeValuationInitial();
760  // we avoid dividing by infinite value
761  if (!isinf(maxVal)) {
762  break;
763  }
764  }
765  // Filling the internal nodes
766  for (size_t i = 0; i < nbrNodes; i++) {
767  DendroNodeType &curNode = *dendroNodes.at(i);
768  if (curNode.getIsInternalNode()) {
769  double newVal =
770  double(curNode.getInternalNodeValuationInitial() / maxVal);
771  curNode.setInternalNodeValuationFinal(newVal);
772  }
773  }
774  this->putValuationsFinalInInitial();
775  this->reorganize();
776  } else {
777  std::cerr << "normalize(const std::string typeOfNormalization) -> "
778  "typeOfNormalization must be chosen within: "
779  << "reg, maxnorm" << endl;
780  }
781  }
782 #ifndef SWIG
786  {
787  for (size_t i = 0; i < nbrNodes; i++) {
788  DendroNodeType &curNode = *dendroNodes[i];
789  if (curNode.getIsInternalNode() == 1) {
790  double temp = curNode.getInternalNodeValuationFinal();
791  curNode.setInternalNodeValuationInitial(temp);
792  curNode.setInternalNodeValuationFinal(0);
793  }
794  }
795  };
796  void setValuationsToZero()
797  {
798  for (size_t i = 0; i < nbrNodes; i++) {
799  DendroNodeType &curNode = *dendroNodes[i];
800  if (curNode.getIsInternalNode() == 1) {
801  curNode.setValuation(0);
802  }
803  if (curNode.getIsInternalNode() == 0) {
804  curNode.setValuation(curNode.getLeafValuation());
805  }
806  }
807  }
808 #endif
812  double lbd)
813  { // Dendrogram<MarkerLabelT,NodeT,ValGraphT>&
814  // dendrogram
815  std::vector<DendroNodeType *> &dendroNodes = this->getDendroNodes();
816  size_t nodeNbr = this->getNbrNodes();
817  // size_t internalNodesNbr = associated_mst.getEdgeNbr();
818  // size_t leavesNbr = associated_mst.getNodeNbr();
819  for (size_t i = 0; i < nodeNbr; i++) {
820  DendroNodeType &curNode = *dendroNodes.at(i);
821  if (curNode.getInternalNodeValuationInitial() >= lbd &&
822  curNode.getIsInternalNode() == true) {
823  MarkerLabelT srcToRemove = curNode.getNeighborLeft()->getLabel();
824  MarkerLabelT targetToRemove = curNode.getNeighborRight()->getLabel();
825  associated_mst.removeEdge(srcToRemove, targetToRemove);
826  associated_mst.removeEdge(targetToRemove, srcToRemove);
827  }
828  }
829  };
830 
845  const std::string typeOfHierarchy, const int nParam = 50,
847  const std::string typeOfTransform = "erode",
848  const StrElt & se = DEFAULT_SE)
849  {
850  // sort by increasing values of internalNodeValuationInitial
851  this->sortNodes();
852  std::vector<DendroNodeType *> &dendroNodes = this->getDendroNodes();
853  if (typeOfHierarchy == "none") {
854  } else if (typeOfHierarchy == "surfacic") {
855  // only one dendroNodes traversal, and only on internal nodes
856  // because leaves already have surface as valuations
857  for (size_t i = 0; i < dendroNodes.size(); i++) {
858  DendroNodeType &curNode = *dendroNodes[i];
859  // we verify that it's an internal node
860  if (curNode.getIsInternalNode() == 1) {
861  if (curNode.getValuation() == 0) {
862  double childLeftValuation =
863  curNode.getChildLeft()->getValuation();
864  double childRightValuation =
865  curNode.getChildRight()->getValuation();
866  curNode.setValuation(childLeftValuation + childRightValuation);
867  curNode.setInternalNodeValuationFinal(
868  fmin(childLeftValuation, childRightValuation));
869  } else {
870  curNode.setInternalNodeValuationFinal(
872  }
873  }
874  }
875  this->putValuationsFinalInInitial();
876  this->reorganize();
877  this->setValuationsToZero();
878  } else if (typeOfHierarchy == "surfacicImageReturn") {
879  // only one dendroNodes traversal, and only on internal nodes,
880  // because leaves already have surface as valuations
881  for (size_t i = 0; i < dendroNodes.size(); i++) {
882  DendroNodeType &curNode = *dendroNodes[i];
883  // we verify that it's an internal node
884  if (curNode.getIsInternalNode() == 1) {
885  std::vector<MarkerLabelT> lookupChildLeft =
886  curNode.getChildLeft()->getLookupProgeny();
887  std::vector<MarkerLabelT> lookupChildRight =
888  curNode.getChildRight()->getLookupProgeny();
889 
890  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
891  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
892  lookupMapChildLeft[j] = lookupChildLeft.at(j);
893  }
894 
895  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
896  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
897  lookupMapChildRight[j] = lookupChildRight.at(j);
898  }
899 
900  smil::Image<MarkerLabelT> imTmpLeft(imMosa);
901  smil::Image<MarkerLabelT> imTmpRight(imMosa);
902 
903  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
904  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
905 
906  if (typeOfTransform == "erode") {
907  erode(imTmpLeft, imTmpLeft, se);
908  erode(imTmpRight, imTmpRight, se);
909  } else if (typeOfTransform == "dilate") {
910  dilate(imTmpLeft, imTmpLeft, se);
911  dilate(imTmpRight, imTmpRight, se);
912  } else if (typeOfTransform == "open") {
913  open(imTmpLeft, imTmpLeft, se);
914  open(imTmpRight, imTmpRight, se);
915  } else if (typeOfTransform == "close") {
916  close(imTmpLeft, imTmpLeft, se);
917  close(imTmpRight, imTmpRight, se);
918  } else {
919  cout << "Please choose typeOfTransform in the following: erode, "
920  "dilate, open, close"
921  << endl;
922  }
923  double childLeftSurf = area(imTmpLeft);
924  double childRightSurf = area(imTmpRight);
925 
926  double childLeftValuation = childLeftSurf;
927  double childRightValuation = childRightSurf;
928  curNode.setValuation(childLeftValuation + childRightValuation);
929  curNode.setInternalNodeValuationFinal(
930  fmin(childLeftValuation, childRightValuation));
931  }
932  }
933  this->putValuationsFinalInInitial();
934  this->reorganize();
935  this->setValuationsToZero();
936  } else if (typeOfHierarchy == "stochasticSurfacic") {
937  // First dendroNodes traversal, to get the "surfacic" dendrogram
938  for (size_t i = 0; i < dendroNodes.size(); i++) {
939  DendroNodeType &curNode = *dendroNodes[i];
940  // we verify that it's an internal node
941  if (curNode.getIsInternalNode() == 1) {
942  if (curNode.getValuation() == 0) {
943  double childLeftValuation =
944  curNode.getChildLeft()->getValuation();
945  double childRightValuation =
946  curNode.getChildRight()->getValuation();
947  curNode.setValuation(childLeftValuation + childRightValuation);
948  } else {
949  curNode.setInternalNodeValuationFinal(
951  }
952  }
953  }
954 
955  DendroNodeType &curNodeTmp = *dendroNodes[1];
956  // get the total surface of the domain
957  double totalSurface = curNodeTmp.getAncestor()->getValuation();
958  // Second dendroNodes traversal, to get the stochasticSurfacic
959  // dendrogram
960  for (size_t i = 0; i < dendroNodes.size(); i++) {
961  DendroNodeType &curNode = *dendroNodes[i];
962  // we verify that it's an internal node
963  if (curNode.getIsInternalNode() == 1) {
964  double childLeftSurf = curNode.getChildLeft()->getValuation();
965  double childRightSurf = curNode.getChildRight()->getValuation();
966  double newVal =
967  1 - pow(1 - (childLeftSurf / totalSurface), nParam) -
968  pow(1 - (childRightSurf / totalSurface), nParam) +
969  pow(1 - ((childLeftSurf + childRightSurf) / totalSurface),
970  nParam);
971  if (newVal <= 0) {
972  newVal = 0.000000001;
973  }
974  curNode.setInternalNodeValuationFinal(newVal);
975  }
976  }
977  this->putValuationsFinalInInitial();
978  this->reorganize();
979  this->setValuationsToZero();
980  } else if (typeOfHierarchy == "stochasticExtinctionSurfacic") {
981  // First dendroNodes traversal, to get the "surfacic" dendrogram
982  for (size_t i = 0; i < dendroNodes.size(); i++) {
983  DendroNodeType &curNode = *dendroNodes[i];
984  // we verify that it's an internal node
985  if (curNode.getIsInternalNode() == 1) {
986  if (curNode.getValuation() == 0) {
987  double childLeftValuation =
988  curNode.getChildLeft()->getValuation();
989  double childRightValuation =
990  curNode.getChildRight()->getValuation();
991  curNode.setValuation(childLeftValuation + childRightValuation);
992  } else {
993  curNode.setInternalNodeValuationFinal(
995  }
996  }
997  }
998 
999  DendroNodeType &curNodeTmp = *dendroNodes[1];
1000  // get the total surface of the domain
1001  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1002  // Second dendroNodes traversal, to get the stochasticSurfacic
1003  // dendrogram
1004  for (size_t i = 0; i < dendroNodes.size(); i++) {
1005  DendroNodeType &curNode = *dendroNodes[i];
1006  // we verify that it's an internal node
1007  if (curNode.getIsInternalNode() == 1) {
1008  double childLeftSurf = curNode.getChildLeft()->getValuation();
1009  double childRightSurf = curNode.getChildRight()->getValuation();
1010  double newVal =
1011  1 -
1012  pow(1 - std::min(childLeftSurf, childRightSurf) / totalSurface,
1013  nParam);
1014  curNode.setInternalNodeValuationFinal(newVal);
1015  }
1016  }
1017  this->putValuationsFinalInInitial();
1018  this->reorganize();
1019  this->setValuationsToZero();
1020  } else if (typeOfHierarchy == "stochasticSurfacicImageReturn") {
1021  // First dendroNodes traversal, to get the "surfacic" dendrogram
1022  for (size_t i = 0; i < dendroNodes.size(); i++) {
1023  DendroNodeType &curNode = *dendroNodes[i];
1024  // we verify that it's an internal node
1025  if (curNode.getIsInternalNode() == 1) {
1026  if (curNode.getValuation() == 0) {
1027  double childLeftValuation =
1028  curNode.getChildLeft()->getValuation();
1029  double childRightValuation =
1030  curNode.getChildRight()->getValuation();
1031  curNode.setValuation(childLeftValuation + childRightValuation);
1032  }
1033  }
1034  }
1035 
1036  DendroNodeType &curNodeTmp = *dendroNodes[1];
1037  // get the total surface of the domain
1038  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1039  // Second dendroNodes traversal, to get the stochasticSurfacic
1040  // dendrogram
1041  for (size_t i = 0; i < dendroNodes.size(); i++) {
1042  DendroNodeType &curNode = *dendroNodes[i];
1043  // we verify that it's an internal node
1044  if (curNode.getIsInternalNode() == 1) {
1045  std::vector<MarkerLabelT> lookupChildLeft =
1046  curNode.getChildLeft()->getLookupProgeny();
1047  std::vector<MarkerLabelT> lookupChildRight =
1048  curNode.getChildRight()->getLookupProgeny();
1049 
1050  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
1051  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
1052  lookupMapChildLeft[j] = lookupChildLeft.at(j);
1053  }
1054 
1055  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
1056  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
1057  lookupMapChildRight[j] = lookupChildRight.at(j);
1058  }
1059 
1060  smil::Image<MarkerLabelT> imTmpLeft(imMosa);
1061  smil::Image<MarkerLabelT> imTmpRight(imMosa);
1062 
1063  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
1064  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
1065 
1066  if (typeOfTransform == "erode") {
1067  erode(imTmpLeft, imTmpLeft, se);
1068  erode(imTmpRight, imTmpRight, se);
1069  } else if (typeOfTransform == "dilate") {
1070  dilate(imTmpLeft, imTmpLeft, se);
1071  dilate(imTmpRight, imTmpRight, se);
1072  } else if (typeOfTransform == "open") {
1073  open(imTmpLeft, imTmpLeft, se);
1074  open(imTmpRight, imTmpRight, se);
1075  } else if (typeOfTransform == "close") {
1076  close(imTmpLeft, imTmpLeft, se);
1077  close(imTmpRight, imTmpRight, se);
1078  } else {
1079  cout << "Please choose typeOfTransform in the following: erode, "
1080  "dilate, open, close"
1081  << endl;
1082  }
1083  double childLeftSurf = area(imTmpLeft);
1084  double childRightSurf = area(imTmpRight);
1085 
1086  double newVal =
1087  1 - pow(1 - (childLeftSurf / totalSurface), nParam) -
1088  pow(1 - (childRightSurf / totalSurface), nParam) +
1089  pow(1 - ((childLeftSurf + childRightSurf) / totalSurface),
1090  nParam);
1091  if (newVal <= 0) {
1092  newVal = 0.000000001;
1093  }
1094  curNode.setInternalNodeValuationFinal(newVal);
1095  }
1096  }
1097  this->putValuationsFinalInInitial();
1098  this->reorganize();
1099  this->setValuationsToZero();
1100  } else if (typeOfHierarchy == "stochasticExtinctionSurfacicImageReturn") {
1101  for (size_t i = 0; i < dendroNodes.size(); i++) {
1102  DendroNodeType &curNode = *dendroNodes[i];
1103  // we verify that it's an internal node
1104  if (curNode.getIsInternalNode() == 1) {
1105  if (curNode.getValuation() == 0) {
1106  double childLeftValuation =
1107  curNode.getChildLeft()->getValuation();
1108  double childRightValuation =
1109  curNode.getChildRight()->getValuation();
1110  curNode.setValuation(childLeftValuation + childRightValuation);
1111  } else {
1112  curNode.setInternalNodeValuationFinal(
1114  }
1115  }
1116  }
1117 
1118  DendroNodeType &curNodeTmp = *dendroNodes[1];
1119  // get the total surface of the domain
1120  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1121  // Second dendroNodes traversal, to get the stochasticSurfacic
1122  // dendrogram
1123  for (size_t i = 0; i < dendroNodes.size(); i++) {
1124  DendroNodeType &curNode = *dendroNodes[i];
1125  // we verify that it's an internal node
1126  if (curNode.getIsInternalNode() == 1) {
1127  std::vector<MarkerLabelT> lookupChildLeft =
1128  curNode.getChildLeft()->getLookupProgeny();
1129  std::vector<MarkerLabelT> lookupChildRight =
1130  curNode.getChildRight()->getLookupProgeny();
1131 
1132  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
1133  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
1134  lookupMapChildLeft[j] = lookupChildLeft.at(j);
1135  }
1136 
1137  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
1138  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
1139  lookupMapChildRight[j] = lookupChildRight.at(j);
1140  }
1141 
1142  smil::Image<MarkerLabelT> imTmpLeft(imMosa);
1143  smil::Image<MarkerLabelT> imTmpRight(imMosa);
1144 
1145  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
1146  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
1147 
1148  if (typeOfTransform == "erode") {
1149  erode(imTmpLeft, imTmpLeft, se);
1150  erode(imTmpRight, imTmpRight, se);
1151  } else if (typeOfTransform == "dilate") {
1152  dilate(imTmpLeft, imTmpLeft, se);
1153  dilate(imTmpRight, imTmpRight, se);
1154  } else if (typeOfTransform == "open") {
1155  open(imTmpLeft, imTmpLeft, se);
1156  open(imTmpRight, imTmpRight, se);
1157  } else if (typeOfTransform == "close") {
1158  close(imTmpLeft, imTmpLeft, se);
1159  close(imTmpRight, imTmpRight, se);
1160  } else {
1161  cout << "Please choose typeOfTransform in the following: erode, "
1162  "dilate, open, close"
1163  << endl;
1164  }
1165  double childLeftSurf = area(imTmpLeft);
1166  double childRightSurf = area(imTmpRight);
1167 
1168  double newVal =
1169  1 -
1170  pow(1 - std::min(childLeftSurf, childRightSurf) / totalSurface,
1171  nParam);
1172 
1173  if (newVal <= 0) {
1174  newVal = 0.000000001;
1175  }
1176  curNode.setInternalNodeValuationFinal(newVal);
1177  }
1178  }
1179  this->putValuationsFinalInInitial();
1180  this->reorganize();
1181  this->setValuationsToZero();
1182  } else if (typeOfHierarchy == "stochasticSurfacicImageReturnNonPoint") {
1183  // First dendroNodes traversal, to get the "surfacic" dendrogram
1184  for (size_t i = 0; i < dendroNodes.size(); i++) {
1185  DendroNodeType &curNode = *dendroNodes[i];
1186  // we verify that it's an internal node
1187  if (curNode.getIsInternalNode() == 1) {
1188  if (curNode.getValuation() == 0) {
1189  double childLeftValuation =
1190  curNode.getChildLeft()->getValuation();
1191  double childRightValuation =
1192  curNode.getChildRight()->getValuation();
1193  curNode.setValuation(childLeftValuation + childRightValuation);
1194  }
1195  }
1196  }
1197 
1198  DendroNodeType &curNodeTmp = *dendroNodes[1];
1199  // get the total surface of the domain
1200  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1201  smil::Image<MarkerLabelT> imTmpLeft(imMosa);
1202  smil::Image<MarkerLabelT> imTmpRight(imMosa);
1203  smil::Image<MarkerLabelT> imTmpInter(imMosa);
1204  smil::Image<MarkerLabelT> imTmpInterbis(imMosa);
1205 
1206  // Second dendroNodes traversal, to get the stochasticSurfacic
1207  // dendrogram
1208  for (size_t i = 0; i < dendroNodes.size(); i++) {
1209  DendroNodeType &curNode = *dendroNodes[i];
1210  // we verify that it's an internal node
1211  if (curNode.getIsInternalNode() == 1) {
1212  std::vector<MarkerLabelT> lookupChildLeft =
1213  curNode.getChildLeft()->getLookupProgeny();
1214  std::vector<MarkerLabelT> lookupChildRight =
1215  curNode.getChildRight()->getLookupProgeny();
1216 
1217  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
1218  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
1219  lookupMapChildLeft[j] = lookupChildLeft.at(j);
1220  }
1221 
1222  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
1223  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
1224  lookupMapChildRight[j] = lookupChildRight.at(j);
1225  }
1226 
1227  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
1228  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
1229  dilate(imTmpLeft, imTmpLeft, se);
1230  dilate(imTmpRight, imTmpRight, se);
1231 
1232  mul(imTmpLeft, imTmpRight, imTmpInter);
1233  double LRSurf = area(imTmpInter);
1234 
1235  MarkerLabelT one = 1;
1236  MarkerLabelT zero = 0;
1237 
1238  sub(imTmpLeft, imTmpRight, imTmpInter);
1239 
1240  compare(imTmpInter, "==", one, one, zero, imTmpInterbis);
1241  double L_RSurf = area(imTmpInterbis);
1242 
1243  sub(imTmpRight, imTmpLeft, imTmpInter);
1244  compare(imTmpInter, "==", one, one, zero, imTmpInterbis);
1245  double R_LSurf = area(imTmpInterbis);
1246 
1247  double newVal =
1248  pow(1 - (LRSurf / totalSurface), nParam) *
1249  (1 - pow(1 - (L_RSurf / (totalSurface - LRSurf)), nParam) -
1250  pow(1 - (R_LSurf / (totalSurface - LRSurf)), nParam) +
1251  pow(1 - (L_RSurf + R_LSurf) / (totalSurface - LRSurf),
1252  nParam));
1253 
1254  if (newVal <= 0) {
1255  newVal = 0.000000001;
1256  }
1257  curNode.setInternalNodeValuationFinal(newVal);
1258  }
1259  }
1260  this->putValuationsFinalInInitial();
1261  this->reorganize();
1262  this->setValuationsToZero();
1263  } else if (typeOfHierarchy == "volumic") {
1264  // First dendroNodes traversal, to get the "surfacic" dendrogram
1265  for (size_t i = 0; i < dendroNodes.size(); i++) {
1266  DendroNodeType &curNode = *dendroNodes[i];
1267  // we verify that it's an internal node
1268  if (curNode.getIsInternalNode() == 1) {
1269  if (curNode.getValuation() == 0) {
1270  double childLeftValuation =
1271  curNode.getChildLeft()->getValuation();
1272  double childRightValuation =
1273  curNode.getChildRight()->getValuation();
1274  curNode.setValuation(childLeftValuation + childRightValuation);
1275  }
1276  }
1277  }
1278  // Second dendroNodes traversal, to get the "volumic" dendrogram
1279  for (size_t i = 0; i < dendroNodes.size(); i++) {
1280  DendroNodeType &curNode = *dendroNodes[i];
1281  // we verify that the node is not the ancestor
1282  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1283  double height =
1284  curNode.getFather()->getInternalNodeValuationInitial();
1285  double surface = curNode.getValuation();
1286  curNode.setValuation(height * surface);
1287  }
1288  }
1289  // Third dendroNodes traversal, to get the final
1290  // valuation of internal nodes
1291  for (size_t i = 0; i < dendroNodes.size(); i++) {
1292  DendroNodeType &curNode = *dendroNodes[i];
1293  // we verify that it's an internal node
1294  if (curNode.getIsInternalNode() == 1) {
1295  double childLeftValuation = curNode.getChildLeft()->getValuation();
1296  double childRightValuation =
1297  curNode.getChildRight()->getValuation();
1298  double height = curNode.getInternalNodeValuationInitial();
1299  curNode.setInternalNodeValuationFinal(
1300  height * fmin(childLeftValuation, childRightValuation));
1301  }
1302  }
1303  this->putValuationsFinalInInitial();
1304  this->reorganize();
1305  this->setValuationsToZero();
1306  } else if (typeOfHierarchy == "volumicImageReturn") {
1307  // First dendroNodes traversal, to get the "surfacic" dendrogram
1308  for (size_t i = 0; i < dendroNodes.size(); i++) {
1309  DendroNodeType &curNode = *dendroNodes[i];
1310  // we verify that it's an internal node
1311  if (curNode.getIsInternalNode() == 1) {
1312  std::vector<MarkerLabelT> lookupChildLeft =
1313  curNode.getChildLeft()->getLookupProgeny();
1314  std::vector<MarkerLabelT> lookupChildRight =
1315  curNode.getChildRight()->getLookupProgeny();
1316 
1317  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
1318  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
1319  lookupMapChildLeft[j] = lookupChildLeft.at(j);
1320  }
1321 
1322  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
1323  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
1324  lookupMapChildRight[j] = lookupChildRight.at(j);
1325  }
1326 
1327  smil::Image<MarkerLabelT> imTmpLeft(imMosa);
1328  smil::Image<MarkerLabelT> imTmpRight(imMosa);
1329 
1330  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
1331  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
1332 
1333  if (typeOfTransform == "erode") {
1334  erode(imTmpLeft, imTmpLeft, se);
1335  erode(imTmpRight, imTmpRight, se);
1336  } else if (typeOfTransform == "dilate") {
1337  dilate(imTmpLeft, imTmpLeft, se);
1338  dilate(imTmpRight, imTmpRight, se);
1339  } else if (typeOfTransform == "open") {
1340  open(imTmpLeft, imTmpLeft, se);
1341  open(imTmpRight, imTmpRight, se);
1342  } else if (typeOfTransform == "close") {
1343  close(imTmpLeft, imTmpLeft, se);
1344  close(imTmpRight, imTmpRight, se);
1345  } else {
1346  cout << "Please choose typeOfTransform in the following: erode, "
1347  "dilate, open, close"
1348  << endl;
1349  }
1350  double childLeftSurf = area(imTmpLeft);
1351  double childRightSurf = area(imTmpRight);
1352 
1353  curNode.setValuation(childLeftSurf + childRightSurf);
1354  }
1355  }
1356  // Second dendroNodes traversal, to get the "volumic" dendrogram
1357  for (size_t i = 0; i < dendroNodes.size(); i++) {
1358  DendroNodeType &curNode = *dendroNodes[i];
1359  // we verify that the node is not the ancestor
1360  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1361  double height =
1362  curNode.getFather()->getInternalNodeValuationInitial();
1363  double surface = curNode.getValuation();
1364  curNode.setValuation(height * surface);
1365  }
1366  }
1367  // Third dendroNodes traversal, to get the final
1368  // valuation of internal nodes
1369  for (size_t i = 0; i < dendroNodes.size(); i++) {
1370  DendroNodeType &curNode = *dendroNodes[i];
1371  // we verify that it's an internal node
1372  if (curNode.getIsInternalNode() == 1) {
1373  double childLeftValuation = curNode.getChildLeft()->getValuation();
1374  double childRightValuation =
1375  curNode.getChildRight()->getValuation();
1376  double height = curNode.getInternalNodeValuationInitial();
1377  double newVal =
1378  height * fmin(childLeftValuation, childRightValuation);
1379 
1380  if (newVal <= 0) {
1381  newVal = 0.000000001;
1382  };
1383  curNode.setInternalNodeValuationFinal(newVal);
1384  }
1385  }
1386  this->putValuationsFinalInInitial();
1387  this->reorganize();
1388  this->setValuationsToZero();
1389  } else if (typeOfHierarchy == "stochasticVolumic") {
1390  double totalSurface = 0;
1391  double totalDepth = 0;
1392 
1393  // First dendroNodes traversal, to get the "surfacic" dendrogram
1394  for (size_t i = 0; i < dendroNodes.size(); i++) {
1395  DendroNodeType &curNode = *dendroNodes[i];
1396  // we verify that it's an internal node
1397  if (curNode.getIsInternalNode() == 1) {
1398  if (curNode.getValuation() == 0) {
1399  double childLeftValuation =
1400  curNode.getChildLeft()->getValuation();
1401  double childRightValuation =
1402  curNode.getChildRight()->getValuation();
1403  curNode.setValuation(childLeftValuation + childRightValuation);
1404  }
1405  }
1406  }
1407 
1408  DendroNodeType &curNodeTmp = *dendroNodes[1];
1409  // get the total surface of the domain
1410  totalSurface = curNodeTmp.getAncestor()->getValuation();
1411  // get the total depth of the domain
1412  totalDepth =
1413  curNodeTmp.getAncestor()->getInternalNodeValuationInitial();
1414  double totalVolume = totalSurface * totalDepth;
1415  // Second dendroNodes traversal, to get the "volumic" dendrogram
1416  for (size_t i = 0; i < dendroNodes.size(); i++) {
1417  DendroNodeType &curNode = *dendroNodes[i];
1418  // we verify that the node is not the ancestor
1419  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1420  double height =
1421  curNode.getFather()->getInternalNodeValuationInitial();
1422  double surface = curNode.getValuation();
1423  curNode.setValuation(height * surface);
1424  }
1425  }
1426 
1427  // Third dendroNodes traversal, to get the
1428  // stochasticVolumic dendrogram
1429  for (size_t i = 0; i < dendroNodes.size(); i++) {
1430  DendroNodeType &curNode = *dendroNodes[i];
1431  // we verify that it's an internal node
1432  if (curNode.getIsInternalNode() == 1) {
1433  double childLeftVol = curNode.getChildLeft()->getValuation();
1434  double childRightVol = curNode.getChildRight()->getValuation();
1435 
1436  double newVal =
1437  1 - pow(1 - (childLeftVol / totalVolume), nParam) -
1438  pow(1 - (childRightVol / totalVolume), nParam) +
1439  pow(1 - (childLeftVol + childRightVol) / totalVolume, nParam);
1440  curNode.setInternalNodeValuationFinal(newVal);
1441  }
1442  }
1443  this->putValuationsFinalInInitial();
1444  this->reorganize();
1445  this->setValuationsToZero();
1446  } else if (typeOfHierarchy == "stochasticExtinctionVolumic") {
1447  double totalSurface = 0;
1448  double totalDepth = 0;
1449 
1450  // First dendroNodes traversal, to get the "surfacic" dendrogram
1451  for (size_t i = 0; i < dendroNodes.size(); i++) {
1452  DendroNodeType &curNode = *dendroNodes[i];
1453  // we verify that it's an internal node
1454  if (curNode.getIsInternalNode() == 1) {
1455  if (curNode.getValuation() == 0) {
1456  double childLeftValuation =
1457  curNode.getChildLeft()->getValuation();
1458  double childRightValuation =
1459  curNode.getChildRight()->getValuation();
1460  curNode.setValuation(childLeftValuation + childRightValuation);
1461  }
1462  }
1463  }
1464 
1465  DendroNodeType &curNodeTmp = *dendroNodes[1];
1466  // get the total surface of the domain
1467  totalSurface = curNodeTmp.getAncestor()->getValuation();
1468  // get the total depth of the domain
1469  totalDepth =
1470  curNodeTmp.getAncestor()->getInternalNodeValuationInitial();
1471 
1472  double totalVolume = totalSurface * totalDepth;
1473  // Second dendroNodes traversal, to get the "volumic" dendrogram
1474  for (size_t i = 0; i < dendroNodes.size(); i++) {
1475  DendroNodeType &curNode = *dendroNodes[i];
1476  // we verify that the node is not the ancestor
1477  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1478  // double height =
1479  double height =
1480  curNode.getFather()->getInternalNodeValuationInitial();
1481  double surface = curNode.getValuation();
1482  curNode.setValuation(height * surface);
1483  }
1484  }
1485  // Third dendroNodes traversal, to get the stochasticVolumic dendrogram
1486  for (size_t i = 0; i < dendroNodes.size(); i++) {
1487  DendroNodeType &curNode = *dendroNodes[i];
1488  // we verify that it's an internal node
1489  if (curNode.getIsInternalNode() == 1) {
1490  double childLeftVol = curNode.getChildLeft()->getValuation();
1491  double childRightVol = curNode.getChildRight()->getValuation();
1492 
1493  double newVal =
1494  1 - pow(1 - std::min(childLeftVol, childRightVol) / totalVolume,
1495  nParam);
1496 
1497  if (newVal <= 0) {
1498  newVal = 0.000000001;
1499  }
1500  curNode.setInternalNodeValuationFinal(newVal);
1501  }
1502  }
1503  this->putValuationsFinalInInitial();
1504  this->reorganize();
1505  this->setValuationsToZero();
1506  } else if (typeOfHierarchy == "stochasticVolumicImageReturn") {
1507  double totalSurface = 0;
1508  double totalDepth = 0;
1509 
1510  // First dendroNodes traversal, to get the "surfacic" dendrogram
1511  for (size_t i = 0; i < dendroNodes.size(); i++) {
1512  DendroNodeType &curNode = *dendroNodes[i];
1513  // we verify that it's an internal node
1514  if (curNode.getIsInternalNode() == 1) {
1515  std::vector<MarkerLabelT> lookupChildLeft =
1516  curNode.getChildLeft()->getLookupProgeny();
1517  std::vector<MarkerLabelT> lookupChildRight =
1518  curNode.getChildRight()->getLookupProgeny();
1519 
1520  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildLeft;
1521  for (MarkerLabelT j = 0; j < lookupChildLeft.size(); j++) {
1522  lookupMapChildLeft[j] = lookupChildLeft.at(j);
1523  }
1524 
1525  std::map<MarkerLabelT, MarkerLabelT> lookupMapChildRight;
1526  for (MarkerLabelT j = 0; j < lookupChildRight.size(); j++) {
1527  lookupMapChildRight[j] = lookupChildRight.at(j);
1528  }
1529 
1530  smil::Image<UINT8> imTmpLeft(imMosa);
1531  smil::Image<UINT8> imTmpRight(imMosa);
1532 
1533  applyLookup(imMosa, lookupMapChildLeft, imTmpLeft);
1534  applyLookup(imMosa, lookupMapChildRight, imTmpRight);
1535 
1536  if (typeOfTransform == "erode") {
1537  erode(imTmpLeft, imTmpLeft, se);
1538  erode(imTmpRight, imTmpRight, se);
1539  } else if (typeOfTransform == "dilate") {
1540  dilate(imTmpLeft, imTmpLeft, se);
1541  dilate(imTmpRight, imTmpRight, se);
1542  } else if (typeOfTransform == "open") {
1543  open(imTmpLeft, imTmpLeft, se);
1544  open(imTmpRight, imTmpRight, se);
1545  } else if (typeOfTransform == "close") {
1546  close(imTmpLeft, imTmpLeft, se);
1547  close(imTmpRight, imTmpRight, se);
1548  } else {
1549  cout << "Please choose typeOfTransform in the following: erode, "
1550  "dilate, open, close"
1551  << endl;
1552  }
1553  double childLeftSurf = area(imTmpLeft);
1554  double childRightSurf = area(imTmpRight);
1555  curNode.setValuation(childLeftSurf + childRightSurf);
1556  }
1557  }
1558 
1559  DendroNodeType &curNodeTmp = *dendroNodes[1];
1560  // get the total surface of the domain
1561  totalSurface = curNodeTmp.getAncestor()->getValuation();
1562  // get the total depth of the domain
1563  totalDepth =
1564  curNodeTmp.getAncestor()->getInternalNodeValuationInitial();
1565  double totalVolume = totalSurface * totalDepth;
1566  // Second dendroNodes traversal, to get the "volumic"
1567  for (size_t i = 0; i < dendroNodes.size(); i++) {
1568  DendroNodeType &curNode = *dendroNodes[i];
1569  // we verify that the node is not the ancestor
1570  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1571  double height =
1572  curNode.getFather()->getInternalNodeValuationInitial();
1573  double surface = curNode.getValuation();
1574  curNode.setValuation(height * surface);
1575  }
1576  }
1577  // Third dendroNodes traversal, to get the stochasticVolumic dendrogram
1578  for (size_t i = 0; i < dendroNodes.size(); i++) {
1579  DendroNodeType &curNode = *dendroNodes[i];
1580  // we verify that it's an internal node
1581  if (curNode.getIsInternalNode() == 1) {
1582  double childLeftVol = curNode.getChildLeft()->getValuation();
1583  double childRightVol = curNode.getChildRight()->getValuation();
1584 
1585  double newVal =
1586  1 - pow(1 - (childLeftVol / totalVolume), nParam) -
1587  pow(1 - (childRightVol / totalVolume), nParam) +
1588  pow(1 - (childLeftVol + childRightVol) / totalVolume, nParam);
1589  if (newVal <= 0) {
1590  newVal = 0.000000001;
1591  };
1592  curNode.setInternalNodeValuationFinal(newVal);
1593  }
1594  }
1595  this->putValuationsFinalInInitial();
1596  this->reorganize();
1597  this->setValuationsToZero();
1598  } else if (typeOfHierarchy == "marker") {
1599  // First dendroNodes traversal, to set the markers values
1600  for (size_t i = 0; i < dendroNodes.size(); i++) {
1601  DendroNodeType &curNode = *dendroNodes[i];
1602  // we verify that it's a leaf
1603  if (curNode.getIsInternalNode() == 0) {
1604  if (curNode.getMarker() == 0) {
1605  curNode.setMarker(curNode.getValuation());
1606  }
1607  }
1608  }
1609  // Second dendroNodes traversal, to get the "marker" dendrogram
1610  for (size_t i = 0; i < dendroNodes.size(); i++) {
1611  DendroNodeType &curNode = *dendroNodes[i];
1612  // we verify that it's an internal node
1613  if (curNode.getIsInternalNode() == 1) {
1614  double markerLeft = curNode.getChildLeft()->getMarker();
1615  double markerRight = curNode.getChildRight()->getMarker();
1616  curNode.setMarker(fmax(markerLeft, markerRight));
1617  }
1618  }
1619  // Third dendroNodes traversal, to get the final
1620  // valuation of internal nodes
1621  for (size_t i = 0; i < dendroNodes.size(); i++) {
1622  DendroNodeType &curNode = *dendroNodes[i];
1623  if (curNode.getIsInternalNode() == 1) {
1624  double markerLeft = curNode.getChildLeft()->getMarker();
1625  double markerRight = curNode.getChildRight()->getMarker();
1626  curNode.setInternalNodeValuationFinal(
1627  fmin(markerLeft, markerRight));
1628  }
1629  }
1630  this->putValuationsFinalInInitial();
1631  this->reorganize();
1632  this->setValuationsToZero();
1633  } else if (typeOfHierarchy == "stochasticSurfacicCountMarkers") {
1634  // First dendroNodes traversal, to get the "surfacic" dendrogram
1635  for (size_t i = 0; i < dendroNodes.size(); i++) {
1636  DendroNodeType &curNode = *dendroNodes[i];
1637  curNode.setValuation(curNode.getNbMarkersUnder());
1638  }
1639 
1640  DendroNodeType &curNodeTmp = *dendroNodes[1];
1641  // get the total surface of the domain
1642  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1643  // cout << "totalSurface = " << totalSurface << endl;
1644  // Second dendroNodes traversal, to get the stochasticSurfacic
1645  // dendrogram
1646  for (size_t i = 0; i < dendroNodes.size(); i++) {
1647  DendroNodeType &curNode = *dendroNodes[i];
1648  // we verify that it's an internal node
1649  if (curNode.getIsInternalNode() == 1) {
1650  double childLeftSurf = curNode.getChildLeft()->getValuation();
1651  double childRightSurf = curNode.getChildRight()->getValuation();
1652  double newVal =
1653  1 - pow(1 - (childLeftSurf / totalSurface), nParam) -
1654  pow(1 - (childRightSurf / totalSurface), nParam) +
1655  pow(1 - ((childLeftSurf + childRightSurf) / totalSurface),
1656  nParam);
1657  curNode.setInternalNodeValuationFinal(newVal);
1658  }
1659  }
1660  this->putValuationsFinalInInitial();
1661  this->reorganize();
1662  this->setValuationsToZero();
1663  } else if (typeOfHierarchy == "stochasticSurfacicSumMarkers") {
1664  // First dendroNodes traversal, to get the "surfacic" dendrogram
1665  for (size_t i = 0; i < dendroNodes.size(); i++) {
1666  DendroNodeType &curNode = *dendroNodes[i];
1667  curNode.setValuation(curNode.getNbMarkersUnder());
1668  }
1669 
1670  DendroNodeType &curNodeTmp = *dendroNodes[1];
1671  // get the total surface of the domain
1672  double totalSurface = curNodeTmp.getAncestor()->getValuation();
1673  // Second dendroNodes traversal, to get the stochasticSurfacic
1674  // dendrogram
1675  for (size_t i = 0; i < dendroNodes.size(); i++) {
1676  DendroNodeType &curNode = *dendroNodes[i];
1677  // we verify that it's an internal node
1678  if (curNode.getIsInternalNode() == 1) {
1679  double childLeftSurf = curNode.getChildLeft()->getValuation();
1680  double childRightSurf = curNode.getChildRight()->getValuation();
1681  double newVal =
1682  1 - pow(1 - (childLeftSurf / totalSurface), nParam) -
1683  pow(1 - (childRightSurf / totalSurface), nParam) +
1684  pow(1 - ((childLeftSurf + childRightSurf) / totalSurface),
1685  nParam);
1686  curNode.setInternalNodeValuationFinal(newVal);
1687  }
1688  }
1689  this->putValuationsFinalInInitial();
1690  this->reorganize();
1691  this->setValuationsToZero();
1692  } else if (typeOfHierarchy == "stochasticVolumicCountMarkers") {
1693  // First dendroNodes traversal, to get the "surfacic" dendrogram
1694  for (size_t i = 0; i < dendroNodes.size(); i++) {
1695  DendroNodeType &curNode = *dendroNodes[i];
1696  curNode.setValuation(curNode.getNbMarkersUnder());
1697  }
1698  double totalSurface = 0;
1699  DendroNodeType &curNodeTmp = *dendroNodes[1];
1700  // get the total surface of the domain
1701  totalSurface = curNodeTmp.getAncestor()->getValuation();
1702  // Second dendroNodes traversal, to get the "volumic" dendrogram
1703  for (size_t i = 0; i < dendroNodes.size(); i++) {
1704  DendroNodeType &curNode = *dendroNodes[i];
1705  // we verify that the node is not the ancestor
1706  if (curNode.getLabel() != curNode.getFather()->getLabel()) {
1707  // double height =
1708  // curNode.getFather()->getInternalNodeValuationInitial();
1709  double height =
1710  curNode.getFather()->getInternalNodeValuationInitial();
1711  double surface = curNode.getValuation();
1712  // curNode.setValuation(surface);
1713  curNode.setValuation(height * surface);
1714  }
1715  }
1716 
1717  // Second dendroNodes traversal, to get the stochasticSurfacic
1718  // dendrogram
1719  for (size_t i = 0; i < dendroNodes.size(); i++) {
1720  DendroNodeType &curNode = *dendroNodes[i];
1721  // we verify that it's an internal node
1722  if (curNode.getIsInternalNode() == 1) {
1723  double childLeftSurf = curNode.getChildLeft()->getValuation();
1724  double childRightSurf = curNode.getChildRight()->getValuation();
1725  double newVal =
1726  1 - pow(1 - (childLeftSurf / totalSurface), nParam) -
1727  pow(1 - (childRightSurf / totalSurface), nParam) +
1728  pow(1 - ((childLeftSurf + childRightSurf) / totalSurface),
1729  nParam);
1730  curNode.setInternalNodeValuationFinal(newVal);
1731  }
1732  }
1733  this->putValuationsFinalInInitial();
1734  this->reorganize();
1735  this->setValuationsToZero();
1736  } else if (typeOfHierarchy == "waterfall") {
1737  // dendroNodes traversal, to set the waterfall and diameters values
1738  for (size_t i = 0; i < dendroNodes.size(); i++) {
1739  DendroNodeType &curNode = *dendroNodes[i];
1740  // we verify that it's an internal node
1741  if (curNode.getIsInternalNode() == 1) {
1742  if (curNode.getValuation() == 0) {
1743  double childLeftValuation =
1744  curNode.getChildLeft()->getValuation();
1745  double childRightValuation =
1746  curNode.getChildRight()->getValuation();
1747  double childLeftDiameter =
1748  curNode.getChildLeft()->getInternalNodeValuationFinal();
1749  double childRightDiameter =
1750  curNode.getChildRight()->getInternalNodeValuationFinal();
1751 
1752  double newVal = 1 + min(childLeftValuation, childRightValuation);
1753  curNode.setValuation(newVal);
1754  curNode.setInternalNodeValuationFinal(
1755  max(newVal, max(childLeftDiameter, childRightDiameter)));
1756  }
1757  }
1758  }
1759  this->putValuationsFinalInInitial();
1760  this->reorganize();
1761  this->setValuationsToZero();
1762  } else if (typeOfHierarchy == "persistence") {
1763  for (size_t i = 0; i < dendroNodes.size(); i++) {
1764  DendroNodeType &curNode = *dendroNodes[i];
1765  // we verify that it's an internal node
1766  if (curNode.getIsInternalNode() == 1) {
1767  double childLeftINValuation =
1768  curNode.getChildLeft()->getInternalNodeValuationInitial();
1769  double childRightINValuation =
1770  curNode.getChildRight()->getInternalNodeValuationInitial();
1771  double INValuation = curNode.getInternalNodeValuationInitial();
1772 
1773  double newVal = INValuation - std::max(childLeftINValuation,
1774  childRightINValuation);
1775  curNode.setInternalNodeValuationFinal(newVal);
1776  }
1777  }
1778  this->putValuationsFinalInInitial();
1779  this->reorganize();
1780  this->setValuationsToZero();
1781  } else {
1782  cout
1783  << "void "
1784  "Dendrogram::HierarchicalConstruction(Dendrogram& "
1785  "dendrogram,std::string typeOfHierarchy) \n"
1786  << "Please choose one of the following hierarchies: surfacic, "
1787  "volumic, stochasticSurfacic, stochasticVolumic, "
1788  << "surfacicImageReturn, volumicImageReturn, "
1789  "stochasticSurfacicImageReturn, stochasticVolumicImageReturn "
1790  << ",stochasticSurfacicCountMarkers,stochasticVolumicCountMarkers, "
1791  << "stochasticExtinctionSurfacic, stochasticExtinctionVolumic, "
1792  << "stochasticExtinctionSurfacicImageReturn, "
1793  << "waterfall,marker,persistence" << endl;
1794  }
1795  };
1796 
1798  void EnergyConstruction(const std::string typeOfHierarchy)
1799  {
1800  // ,double lamb
1801  // sort by increasing values of internalNodeValuationInitial
1802  this->sortNodes();
1803 
1804  std::vector<DendroNodeType *> &dendroNodes = this->getDendroNodes();
1805 
1806  if (typeOfHierarchy == "none") {
1807  } else if (typeOfHierarchy == "simplifiedMumfordShah") {
1808  for (size_t i = 0; i < dendroNodes.size();
1809  i++) { // dendroNodes traversal
1810  DendroNodeType &curNode = *dendroNodes[i];
1811  if (curNode.getIsInternalNode() == 1) {
1812  std::vector<double> childLeftMoments =
1813  curNode.getChildLeft()->getMoments();
1814  std::vector<double> childRightMoments =
1815  curNode.getChildRight()->getMoments();
1816  std::vector<double> moments = curNode.getMoments();
1817 
1818  double contoursValue = curNode.getContoursSize();
1819  double contoursValueLeft =
1820  curNode.getChildLeft()->getContoursSize();
1821  double contoursValueRight =
1822  curNode.getChildRight()->getContoursSize();
1823 
1824  double fidelityTerm =
1825  moments.at(2) - (moments.at(1) * moments.at(1)) / moments.at(0);
1826  double fidelityTermLeft =
1827  childLeftMoments.at(2) -
1828  (childLeftMoments.at(1) * childLeftMoments.at(1)) /
1829  childLeftMoments.at(0);
1830  double fidelityTermRight =
1831  childRightMoments.at(2) -
1832  (childRightMoments.at(1) * childRightMoments.at(1)) /
1833  childRightMoments.at(0);
1834 
1835  double newLamb =
1836  (fidelityTerm - fidelityTermLeft - fidelityTermRight) /
1837  (contoursValueLeft + contoursValueRight - contoursValue);
1838 
1839  curNode.setInternalNodeValuationFinal(newLamb);
1840  curNode.setValuation(newLamb);
1841  }
1842  }
1843  this->putValuationsFinalInInitial();
1844  this->reorganize();
1845 
1846  // We compute the energies with the new organization of nodes obtained
1847  // with the reorganize function
1848 
1849  // sort nodes by decreasing values of
1850  // internalNodeValuationInitial
1851  this->sortNodes(true);
1852  // initialization of contoursValue and fidelityValue
1853  double contoursValue = 0.0;
1854  DendroNodeType & highestNode = *dendroNodes[0];
1855  std::vector<double> momentsHighestNode = highestNode.getMoments();
1856  double fidelityValue =
1857  momentsHighestNode.at(2) -
1858  (momentsHighestNode.at(1) * momentsHighestNode.at(1)) /
1859  momentsHighestNode.at(0);
1860 
1861  // dendroNodes traversal
1862  for (size_t i = 0; i < dendroNodes.size(); i++) {
1863  DendroNodeType &curNode = *dendroNodes[i];
1864  // we verify that it's an internal node
1865  if (curNode.getIsInternalNode() == 1) {
1866  std::vector<double> moments = curNode.getMoments();
1867  std::vector<double> momentsLeft =
1868  curNode.getChildLeft()->getMoments();
1869  std::vector<double> momentsRight =
1870  curNode.getChildRight()->getMoments();
1871 
1872  double fidelityTerm =
1873  moments.at(2) - (moments.at(1) * moments.at(1)) / moments.at(0);
1874  double fidelityTermLeft =
1875  momentsLeft.at(2) -
1876  (momentsLeft.at(1) * momentsLeft.at(1)) / momentsLeft.at(0);
1877  double fidelityTermRight =
1878  momentsRight.at(2) -
1879  (momentsRight.at(1) * momentsRight.at(1)) / momentsRight.at(0);
1880 
1881  double Lamb = curNode.getInternalNodeValuationInitial();
1882 
1883  contoursValue = contoursValue + curNode.getContoursSize();
1884  fidelityValue = fidelityValue - fidelityTerm + fidelityTermLeft +
1885  fidelityTermRight;
1886 
1887  double newEnergy = fidelityValue + Lamb * contoursValue;
1888  curNode.setEnergy(newEnergy);
1889  }
1890  }
1891 
1892  this->setValuationsToZero();
1893  }
1894 
1895  else if (typeOfHierarchy == "varianceEnergy") {
1896  for (size_t i = 0; i < dendroNodes.size();
1897  i++) { // dendroNodes traversal
1898  DendroNodeType &curNode = *dendroNodes[i];
1899  if (curNode.getIsInternalNode() == 1) {
1900  std::vector<double> childLeftMoments =
1901  curNode.getChildLeft()->getMoments();
1902  std::vector<double> childRightMoments =
1903  curNode.getChildRight()->getMoments();
1904  std::vector<double> moments = curNode.getMoments();
1905 
1906  double fidelityTerm =
1907  moments.at(2) - (moments.at(1) * moments.at(1)) / moments.at(0);
1908  double fidelityTermLeft =
1909  childLeftMoments.at(2) -
1910  (childLeftMoments.at(1) * childLeftMoments.at(1)) /
1911  childLeftMoments.at(0);
1912  double fidelityTermRight =
1913  childRightMoments.at(2) -
1914  (childRightMoments.at(1) * childRightMoments.at(1)) /
1915  childRightMoments.at(0);
1916 
1917  double newLamb =
1918  (fidelityTerm - fidelityTermLeft - fidelityTermRight);
1919 
1920  curNode.setInternalNodeValuationFinal(newLamb);
1921  curNode.setValuation(newLamb);
1922  }
1923  }
1924  this->putValuationsFinalInInitial();
1925  this->reorganize();
1926 
1927  // We compute the energies with the new organization of nodes obtained
1928  // with the reorganize function
1929 
1930  // sort nodes by decreasing values of
1931  // internalNodeValuationInitial
1932  this->sortNodes(true);
1933  // initialization of contoursValue and fidelityValue
1934  double contoursValue = 0.0;
1935  DendroNodeType & highestNode = *dendroNodes[0];
1936  std::vector<double> momentsHighestNode = highestNode.getMoments();
1937  double fidelityValue =
1938  momentsHighestNode.at(2) -
1939  (momentsHighestNode.at(1) * momentsHighestNode.at(1)) /
1940  momentsHighestNode.at(0);
1941 
1942  // dendroNodes traversal
1943  for (size_t i = 0; i < dendroNodes.size(); i++) {
1944  DendroNodeType &curNode = *dendroNodes[i];
1945  // we verify that it's an internal node
1946  if (curNode.getIsInternalNode() == 1) {
1947  std::vector<double> moments = curNode.getMoments();
1948  std::vector<double> momentsLeft =
1949  curNode.getChildLeft()->getMoments();
1950  std::vector<double> momentsRight =
1951  curNode.getChildRight()->getMoments();
1952 
1953  double fidelityTerm =
1954  moments.at(2) - (moments.at(1) * moments.at(1)) / moments.at(0);
1955  double fidelityTermLeft =
1956  momentsLeft.at(2) -
1957  (momentsLeft.at(1) * momentsLeft.at(1)) / momentsLeft.at(0);
1958  double fidelityTermRight =
1959  momentsRight.at(2) -
1960  (momentsRight.at(1) * momentsRight.at(1)) / momentsRight.at(0);
1961 
1962  double Lamb = curNode.getInternalNodeValuationInitial();
1963 
1964  contoursValue = contoursValue + curNode.getContoursSize();
1965  fidelityValue = fidelityValue - fidelityTerm + fidelityTermLeft +
1966  fidelityTermRight;
1967 
1968  double newEnergy = fidelityValue + Lamb * contoursValue;
1969  curNode.setEnergy(newEnergy);
1970  }
1971  }
1972 
1973  this->setValuationsToZero();
1974  }
1975 
1976  else {
1977  cout << "void Dendrogram::EnergyConstruction(Dendrogram& "
1978  "dendrogram,std::string typeOfHierarchy) \n"
1979  << "Please choose one of the following hierarchies: "
1980  "simplifiedMumfordShah, varianceEnergy "
1981  << endl;
1982  }
1983  }
1984 
1986  std::vector<DendroNodeType *> &getDendroNodes()
1987  {
1988  return dendroNodes;
1989  };
1990 
1991  // DendroNodeType*
1992  // &getDendroNodeLabel(std::vector<DendroNodeType*>
1993  // dendroNodes, MarkerLabelT label){
1994  // DendroNodeType &curNode = *dendroNodes[0];
1995  // for (size_t i = 0 ; i<dendroNodes.size() ; i++){
1996  // &curNode = *dendroNodes[i];
1997  // if (curNode.getLabel() == label){
1998  // break;
1999  // }
2000  // }
2001  // return &curNode;
2002  // }
2003  void setNbrNodes(size_t nNbrNodes)
2004  {
2005  nbrNodes = nNbrNodes;
2006  }
2007  size_t getNbrNodes()
2008  {
2009  return nbrNodes;
2010  }
2011  void setNbrMarkers(size_t nNbrMarkers)
2012  {
2013  nbrMarkers = nNbrMarkers;
2014  }
2015  size_t getNbrMarkers()
2016  {
2017  return nbrMarkers;
2018  }
2019 
2027  double getNodeValue(size_t nodeIndex, string nameOfValueWanted)
2028  {
2029  DendroNodeType &curNode = *dendroNodes.at(nodeIndex);
2030  if (nameOfValueWanted == "valuation") {
2031  return curNode.getValuation();
2032  } else if (nameOfValueWanted == "internalNodeValuationInitial") {
2033  return curNode.getInternalNodeValuationInitial();
2034  } else if (nameOfValueWanted == "internalNodeValuationFinal") {
2035  return curNode.getInternalNodeValuationFinal();
2036  } else if (nameOfValueWanted == "label") {
2037  return curNode.getLabel();
2038  } else if (nameOfValueWanted == "labelChildRight") {
2039  return curNode.getChildRight()->getLabel();
2040  } else if (nameOfValueWanted == "labelChildLeft") {
2041  return curNode.getChildLeft()->getLabel();
2042  } else if (nameOfValueWanted == "labelNeighborRight") {
2043  return curNode.getNeighborRight()->getLabel();
2044  } else if (nameOfValueWanted == "labelNeighborLeft") {
2045  return curNode.getNeighborLeft()->getLabel();
2046  } else if (nameOfValueWanted == "labelFather") {
2047  return curNode.getFather()->getLabel();
2048  } else if (nameOfValueWanted == "isInternalNode") {
2049  return curNode.getIsInternalNode();
2050  } else if (nameOfValueWanted == "marker") {
2051  return curNode.getMarker();
2052  } else if (nameOfValueWanted == "energy") {
2053  return curNode.getEnergy();
2054  } else if (nameOfValueWanted == "contours") {
2055  return curNode.getContoursSize();
2056  } else {
2057  cout
2058  << "int Dendrogram::getNodeValue(int nodeIndex,string "
2059  "nameOfValueWanted) -> nameOfValueWanted must be chosen in this "
2060  "list: valuation, internalNodeValuationInitial "
2061  << ", internalNodeValuationInitial, label, marker." << endl;
2062  return 0.0;
2063  }
2064  };
2065 
2074  void setNodeValue(size_t nodeIndex, string nameOfValueWanted, double value)
2075  {
2076  DendroNodeType &curNode = *dendroNodes.at(nodeIndex);
2077  if (nameOfValueWanted == "valuation") {
2078  curNode.setValuation(value);
2079  } else if (nameOfValueWanted == "internalNodeValuationInitial") {
2080  curNode.setInternalNodeValuationInitial(value);
2081  } else if (nameOfValueWanted == "internalNodeValuationFinal") {
2082  curNode.setInternalNodeValuationFinal(value);
2083  } else if (nameOfValueWanted == "label") {
2084  curNode.setLabel(value);
2085  } else if (nameOfValueWanted == "labelChildRight") {
2086  curNode.getChildRight()->setLabel(value);
2087  } else if (nameOfValueWanted == "labelChildLeft") {
2088  curNode.getChildLeft()->setLabel(value);
2089  } else if (nameOfValueWanted == "labelNeighborRight") {
2090  curNode.getNeighborRight()->setLabel(value);
2091  } else if (nameOfValueWanted == "labelNeighborLeft") {
2092  curNode.getNeighborLeft()->setLabel(value);
2093  } else if (nameOfValueWanted == "labelFather") {
2094  curNode.getFather()->setLabel(value);
2095  } else if (nameOfValueWanted == "isInternalNode") {
2096  curNode.setIsInternalNode(value);
2097  } else if (nameOfValueWanted == "marker") {
2098  curNode.setMarker(value);
2099  } else {
2100  cout << "int Dendrogram::setNodeValue(int nodeIndex,string "
2101  "nameOfValueWanted, double value) -> nameOfValueWanted must be "
2102  "chosen in this list: valuation, internalNodeValuationInitial "
2103  << ", internalNodeValuationInitial, label, marker." << endl;
2104  }
2105  };
2106 
2114  std::vector<MarkerLabelT> getLookupProgeny(size_t nodeIndex,
2115  string nameOfValueWanted)
2116  {
2117  DendroNodeType & curNode = *dendroNodes.at(nodeIndex);
2118  std::vector<MarkerLabelT> lookupToReturn;
2119  if (nameOfValueWanted == "current") {
2120  lookupToReturn = curNode.getLookupProgeny();
2121  } else if (nameOfValueWanted == "childLeft") {
2122  lookupToReturn = curNode.getChildLeft()->getLookupProgeny();
2123  } else if (nameOfValueWanted == "childRight") {
2124  lookupToReturn = curNode.getChildRight()->getLookupProgeny();
2125  } else {
2126  cout << "std::vector<MarkerLabelT> getLookupProgeny(size_t "
2127  "nodeIndex,string nameOfValueWanted) -> nameOfValueWanted must "
2128  "be chosen in this list: current, "
2129  << ", childLeft, childRight." << endl;
2130  }
2131  return lookupToReturn;
2132  }
2133 
2141  std::vector<MarkerLabelT> getThreshLookupProgeny(double thresh,
2142  string nameOfValueWanted)
2143  {
2144  std::vector<MarkerLabelT> lookupToReturn;
2145  this->sortNodes(
2146  "true"); // sort nodes by decreasing internalNodeValuationInitial
2147  if (thresh >= 0) {
2148  if (nameOfValueWanted == "current") {
2149  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2150  DendroNodeType &curNode = *dendroNodes.at(i);
2151  if (curNode.getInternalNodeValuationInitial() <= thresh) {
2152  lookupToReturn = curNode.getLookupProgeny();
2153  break;
2154  }
2155  }
2156  }
2157 
2158  else if (nameOfValueWanted == "childLeft") {
2159  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2160  DendroNodeType &curNode = *dendroNodes.at(i);
2161  if (curNode.getInternalNodeValuationInitial() <= thresh) {
2162  lookupToReturn = curNode.getChildLeft()->getLookupProgeny();
2163  break;
2164  }
2165  }
2166  }
2167 
2168  else if (nameOfValueWanted == "childRight") {
2169  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2170  DendroNodeType &curNode = *dendroNodes.at(i);
2171  if (curNode.getInternalNodeValuationInitial() <= thresh) {
2172  lookupToReturn = curNode.getChildRight()->getLookupProgeny();
2173  break;
2174  }
2175  }
2176  } else {
2177  cout << "std::vector<MarkerLabelT> "
2178  "getThreshLookupProgeny(thresh,nameOfValueWanted) -> "
2179  "nameOfValueWanted must be chosen in this list: current, "
2180  << ", childLeft, childRight." << endl;
2181  }
2182  } else {
2183  cout << "getThreshLookupProgeny(thresh,nameOfValueWanted) "
2184  << "-> the thresh specified must be higher than 0" << endl;
2185  }
2186 
2187  return lookupToReturn;
2188  }
2189 
2199  std::vector<std::vector<double>>
2200  getPersistences(string typeOfMatrix = "sparse",
2201  string momentOfLife = "death")
2202  {
2203  size_t nbrRegions = nbrNodes / 2 + 1;
2204  // persistences matrix
2205  std::vector<std::vector<double>> persistences(
2206  nbrRegions, vector<double>(nbrRegions, -1.0));
2207 
2208  if (typeOfMatrix == "sparse" && momentOfLife == "death") {
2209  // First dendrogram traversal to initialize the links between nodes in
2210  // the persistences matrix (we put the value to -0.5 when there is a
2211  // link)
2212  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2213  DendroNodeType &curNode = *dendroNodes.at(i);
2214  // we check if it is a leaf
2215  if (curNode.getIsInternalNode() == 0) {
2216  MarkerLabelT label1 = curNode.getLabel();
2217  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2218  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2219  MarkerLabelT label2 = t;
2220  if (contoursCount.at(t) != 0) {
2221  if (persistences[label1][label2] == -1) {
2222  persistences[label1][label2] = -0.5;
2223  persistences[label2][label1] = -0.5;
2224  }
2225  }
2226  }
2227  }
2228  }
2229  // sort nodes by decreasing internalNodeValuationInitial
2230  this->sortNodes(true);
2231  // Second dendrogram traversal by decreasing
2232  // internalNodeValuationInitial to complete the persistences matrix
2233  MarkerLabelT count = 0;
2234  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2235  DendroNodeType &curNode = *dendroNodes.at(i);
2236  if (curNode.getIsInternalNode() == 1) {
2237  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2238  curNode.getChildLeft()->getLookupProgeny();
2239  std::vector<MarkerLabelT> lookupProgenyChildRight =
2240  curNode.getChildRight()->getLookupProgeny();
2241  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2242  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2243  q++) {
2244  if (lookupProgenyChildLeft.at(p) == 1 &&
2245  lookupProgenyChildRight.at(q) == 1 &&
2246  persistences[p][q] == -0.5) {
2247  persistences[p][q] = nbrRegions - count;
2248  persistences[q][p] = nbrRegions - count;
2249  }
2250  }
2251  }
2252  count++;
2253  }
2254  }
2255  // we set negative values of the matrix to zero
2256  for (size_t i = 0; i < persistences.size(); i++) {
2257  for (size_t j = 0; j < persistences[i].size(); j++)
2258  if (persistences[i][j] <= 0.0) {
2259  persistences[i][j] = 0.0;
2260  persistences[j][i] = 0.0;
2261  }
2262  }
2263  } // end if "sparse" and "death"
2264 
2265  else if (typeOfMatrix == "full" && momentOfLife == "death") {
2266  persistences = std::vector<std::vector<double>>(
2267  nbrRegions, vector<double>(nbrRegions, -0.5));
2268  // sort nodes by decreasing internalNodeValuationInitial
2269  this->sortNodes(true);
2270  // Second dendrogram traversal by decreasing
2271  // internalNodeValuationInitial to complete the persistences matrix
2272  MarkerLabelT count = 0;
2273  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2274  DendroNodeType &curNode = *dendroNodes.at(i);
2275  if (curNode.getIsInternalNode() == 1) {
2276  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2277  curNode.getChildLeft()->getLookupProgeny();
2278  std::vector<MarkerLabelT> lookupProgenyChildRight =
2279  curNode.getChildRight()->getLookupProgeny();
2280  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2281  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2282  q++) {
2283  if (lookupProgenyChildLeft.at(p) == 1 &&
2284  lookupProgenyChildRight.at(q) == 1 &&
2285  persistences[p][q] == -0.5) {
2286  persistences[p][q] = nbrRegions - count;
2287  persistences[q][p] = nbrRegions - count;
2288  }
2289  }
2290  }
2291  count++;
2292  }
2293  }
2294  // we set negative values of the matrix to zero
2295  for (size_t i = 0; i < persistences.size(); i++) {
2296  for (size_t j = 0; j < persistences[i].size(); j++)
2297  if (persistences[i][j] <= 0.0) {
2298  persistences[i][j] = 0.0;
2299  persistences[j][i] = 0.0;
2300  }
2301  }
2302  } // end if "full" and "death"
2303  else if (typeOfMatrix == "sparse" && momentOfLife == "life") {
2304  // First dendrogram traversal initialize the links between nodes in the
2305  // persistences matrix
2306  // (we put the value to -0.5 when there is a link)
2307  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2308  DendroNodeType &curNode = *dendroNodes.at(i);
2309  if (curNode.getIsInternalNode() == 0) { // we check if it is a leaf
2310  MarkerLabelT label1 = curNode.getLabel();
2311  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2312  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2313  MarkerLabelT label2 = t;
2314  if (contoursCount.at(t) != 0) {
2315  if (persistences[label1][label2] == -1) {
2316  persistences[label1][label2] = 1;
2317  persistences[label2][label1] = 1;
2318  }
2319  }
2320  }
2321  }
2322  }
2323  // we set negative values of the matrix to zero
2324  for (size_t i = 0; i < persistences.size(); i++) {
2325  for (size_t j = 0; j < persistences[i].size(); j++)
2326  if (persistences[i][j] <= 0.0) {
2327  persistences[i][j] = 0.0;
2328  persistences[j][i] = 0.0;
2329  }
2330  }
2331  } // end if "sparse" and "life"
2332  else if (typeOfMatrix == "full" && momentOfLife == "life") {
2333  // we begin by generating the same matrix as in the "sparse"/"death"
2334  // version
2335  // First dendrogram traversal to initialize the links between nodes in
2336  // the persistences matrix (we put the value to -0.5 when there is a
2337  // link)
2338  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2339  DendroNodeType &curNode = *dendroNodes.at(i);
2340  if (curNode.getIsInternalNode() == 0) { // we check if it is a leaf
2341  MarkerLabelT label1 = curNode.getLabel();
2342  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2343  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2344  MarkerLabelT label2 = t;
2345  if (contoursCount.at(t) != 0) {
2346  if (persistences[label1][label2] == -1) {
2347  persistences[label1][label2] = -0.5;
2348  persistences[label2][label1] = -0.5;
2349  }
2350  }
2351  }
2352  }
2353  }
2354  // sort nodes by decreasing internalNodeValuationInitial
2355  this->sortNodes(true);
2356  // Second dendrogram traversal by decreasing
2357  // internalNodeValuationInitial to complete the persistences matrix
2358  MarkerLabelT count = 0;
2359  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2360  DendroNodeType &curNode = *dendroNodes.at(i);
2361  if (curNode.getIsInternalNode() == 1) {
2362  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2363  curNode.getChildLeft()->getLookupProgeny();
2364  std::vector<MarkerLabelT> lookupProgenyChildRight =
2365  curNode.getChildRight()->getLookupProgeny();
2366  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2367  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2368  q++) {
2369  if (lookupProgenyChildLeft.at(p) == 1 &&
2370  lookupProgenyChildRight.at(q) == 1 &&
2371  persistences[p][q] == -0.5) {
2372  persistences[p][q] = nbrRegions - count;
2373  persistences[q][p] = nbrRegions - count;
2374  }
2375  }
2376  }
2377  count++;
2378  }
2379  }
2380  // we set negative values of the matrix to zero
2381  for (size_t i = 0; i < persistences.size(); i++) {
2382  for (size_t j = 0; j < persistences[i].size(); j++)
2383  if (persistences[i][j] <= 0.0) {
2384  persistences[i][j] = 0.0;
2385  persistences[j][i] = 0.0;
2386  }
2387  }
2388 
2389  // Now that we have the "sparse"/"death" matrix, we apply a recursive
2390  // lexicographic-distance-based
2391  // algorithm to infer the desired matrix
2392  std::vector<std::vector<double>> persistencesInit(persistences);
2393  std::vector<std::vector<double>> persistencesTmp(
2394  nbrRegions, vector<double>(nbrRegions, 0.0));
2395 
2396  // recursive max-min term-by-term product of the persistences matrix
2397  // until equilibrium state
2398  while (persistences != persistencesTmp) {
2399  persistencesTmp = persistences;
2400  for (size_t i = 0; i < nbrRegions; i++) {
2401  for (size_t j = i; j < nbrRegions; j++) {
2402  std::vector<double> minVec(nbrRegions, 0.0);
2403  for (size_t k = 0; k < nbrRegions; k++) {
2404  minVec[k] = min(persistences[i][k], persistencesInit[k][j]);
2405  }
2406  persistences[i][j] = *max_element(minVec.begin(), minVec.end());
2407  persistences[j][i] = persistences[i][j];
2408  }
2409  }
2410  }
2411  } // end if "full" and "life"
2412  else {
2413  cout << "getPersistences(string typeOfMatrix = ''sparse'',string "
2414  "momentOfLife = ''death'') -> please choose "
2415  << "typeOfMatrix in the following : ''sparse'', ''full''; and "
2416  "momentOfLife in the following : ''death'',''life'' "
2417  << endl;
2418  }
2419 
2420  return persistences;
2421  } // end getPersistences
2422 
2431  std::vector<std::vector<double>>
2432  getSaliences(string typeOfMatrix = "sparse", string momentOfLife = "death")
2433  {
2434  size_t nbrRegions = nbrNodes / 2 + 1;
2435  // saliences matrix
2436  std::vector<std::vector<double>> saliences(
2437  nbrRegions, vector<double>(nbrRegions, -1.0));
2438 
2439  if (typeOfMatrix == "sparse" && momentOfLife == "death") {
2440  // First dendrogram traversal initialize the links between nodes in the
2441  // saliences matrix (we put the value to -0.5 when there is a link)
2442  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2443  DendroNodeType &curNode = *dendroNodes.at(i);
2444  if (curNode.getIsInternalNode() == 0) { // we check if it is a leaf
2445  MarkerLabelT label1 = curNode.getLabel();
2446  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2447  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2448  MarkerLabelT label2 = t;
2449  if (contoursCount.at(t) != 0) {
2450  if (saliences[label1][label2] == -1) {
2451  saliences[label1][label2] = -0.5;
2452  saliences[label2][label1] = -0.5;
2453  }
2454  }
2455  }
2456  }
2457  }
2458  // sort nodes by decreasing internalNodeValuationInitial
2459  this->sortNodes(true);
2460  // Second dendrogram traversal by decreasing
2461  // internalNodeValuationInitial to complete the saliences matrix
2462  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2463  DendroNodeType &curNode = *dendroNodes.at(i);
2464  if (curNode.getIsInternalNode() == 1) {
2465  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2466  curNode.getChildLeft()->getLookupProgeny();
2467  std::vector<MarkerLabelT> lookupProgenyChildRight =
2468  curNode.getChildRight()->getLookupProgeny();
2469  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2470  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2471  q++) {
2472  if (lookupProgenyChildLeft.at(p) == 1 &&
2473  lookupProgenyChildRight.at(q) == 1 &&
2474  saliences[p][q] == -0.5) {
2475  saliences[p][q] = curNode.getInternalNodeValuationInitial();
2476  saliences[q][p] = curNode.getInternalNodeValuationInitial();
2477  }
2478  }
2479  }
2480  }
2481  }
2482 
2483  // we set the non positive values of the matrix to zero
2484  for (size_t i = 0; i < saliences.size(); i++) {
2485  for (size_t j = 0; j < saliences[i].size(); j++)
2486  if (saliences[i][j] <= 0.0) {
2487  saliences[i][j] = 0.0;
2488  saliences[j][i] = 0.0;
2489  }
2490  }
2491  } // end if "sparse" && "death"
2492  else if (typeOfMatrix == "full" && momentOfLife == "death") {
2493  saliences = std::vector<std::vector<double>>(
2494  nbrRegions, vector<double>(nbrRegions, -0.5));
2495 
2496  // sort nodes by decreasing internalNodeValuationInitial
2497  this->sortNodes(true);
2498  // Second dendrogram traversal by decreasing
2499  // internalNodeValuationInitial to complete the saliences matrix
2500  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2501  DendroNodeType &curNode = *dendroNodes.at(i);
2502  if (curNode.getIsInternalNode() == 1) {
2503  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2504  curNode.getChildLeft()->getLookupProgeny();
2505  std::vector<MarkerLabelT> lookupProgenyChildRight =
2506  curNode.getChildRight()->getLookupProgeny();
2507  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2508  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2509  q++) {
2510  if (lookupProgenyChildLeft.at(p) == 1 &&
2511  lookupProgenyChildRight.at(q) == 1 &&
2512  saliences[p][q] == -0.5) {
2513  saliences[p][q] = curNode.getInternalNodeValuationInitial();
2514  saliences[q][p] = curNode.getInternalNodeValuationInitial();
2515  }
2516  }
2517  }
2518  }
2519  }
2520 
2521  // we set the non positive values of the matrix to zero
2522  for (size_t i = 0; i < saliences.size(); i++) {
2523  for (size_t j = 0; j < saliences[i].size(); j++)
2524  if (saliences[i][j] <= 0.0) {
2525  saliences[i][j] = 0.0;
2526  saliences[j][i] = 0.0;
2527  }
2528  }
2529  } // end if "full" && "death"
2530  else if (typeOfMatrix == "sparse" && momentOfLife == "life") {
2531  // First dendrogram traversal initialize the links between nodes in the
2532  // saliences matrix (we put the value to -0.5 when there is a link)
2533  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2534  DendroNodeType &curNode = *dendroNodes.at(i);
2535  if (curNode.getIsInternalNode() == 0) { // we check if it is a leaf
2536  MarkerLabelT label1 = curNode.getLabel();
2537  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2538  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2539  MarkerLabelT label2 = t;
2540  if (contoursCount.at(t) != 0) {
2541  if (saliences[label1][label2] == -1) {
2542  saliences[label1][label2] = 0.01;
2543  saliences[label2][label1] = 0.01;
2544  }
2545  }
2546  }
2547  }
2548  }
2549  // we set the non positive values of the matrix to zero
2550  for (size_t i = 0; i < saliences.size(); i++) {
2551  for (size_t j = 0; j < saliences[i].size(); j++)
2552  if (saliences[i][j] <= 0.0) {
2553  saliences[i][j] = 0.0;
2554  saliences[j][i] = 0.0;
2555  }
2556  }
2557 
2558  } // end if "sparse" && "life"
2559  else if (typeOfMatrix == "full" && momentOfLife == "life") {
2560  saliences = std::vector<std::vector<double>>(
2561  nbrRegions, vector<double>(nbrRegions, -0.5));
2562  // we begin by generating the same matrix as in the "sparse"/"death"
2563  // version
2564  // First dendrogram traversal initialize the links between nodes in the
2565  // saliences matrix
2566  //(we put the value to -0.5 when there is a link)
2567  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2568  DendroNodeType &curNode = *dendroNodes.at(i);
2569  if (curNode.getIsInternalNode() == 0) { // we check if it is a leaf
2570  MarkerLabelT label1 = curNode.getLabel();
2571  vector<MarkerLabelT> contoursCount = curNode.getContoursCount();
2572  for (MarkerLabelT t = 0; t < contoursCount.size(); t++) {
2573  MarkerLabelT label2 = t;
2574  if (contoursCount.at(t) != 0) {
2575  if (saliences[label1][label2] == -1) {
2576  saliences[label1][label2] = -0.5;
2577  saliences[label2][label1] = -0.5;
2578  }
2579  }
2580  }
2581  }
2582  }
2583 
2584  // sort nodes by decreasing internalNodeValuationInitial
2585  this->sortNodes(true);
2586  // Second dendrogram traversal by decreasing
2587  // internalNodeValuationInitial to complete the saliences matrix
2588  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2589  DendroNodeType &curNode = *dendroNodes.at(i);
2590  if (curNode.getIsInternalNode() == 1) {
2591  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2592  curNode.getChildLeft()->getLookupProgeny();
2593  std::vector<MarkerLabelT> lookupProgenyChildRight =
2594  curNode.getChildRight()->getLookupProgeny();
2595  for (MarkerLabelT p = 0; p < lookupProgenyChildLeft.size(); p++) {
2596  for (MarkerLabelT q = 0; q < lookupProgenyChildRight.size();
2597  q++) {
2598  if (lookupProgenyChildLeft.at(p) == 1 &&
2599  lookupProgenyChildRight.at(q) == 1 &&
2600  saliences[p][q] == -0.5) {
2601  saliences[p][q] = curNode.getInternalNodeValuationInitial();
2602  saliences[q][p] = curNode.getInternalNodeValuationInitial();
2603  }
2604  }
2605  }
2606  }
2607  }
2608  // we set the non positive values of the matrix to zero
2609  for (size_t i = 0; i < saliences.size(); i++) {
2610  for (size_t j = 0; j < saliences[i].size(); j++)
2611  if (saliences[i][j] <= 0.0) {
2612  saliences[i][j] = 0.0;
2613  saliences[j][i] = 0.0;
2614  }
2615  }
2616  // Now that we have the "sparse"/"death" matrix, we apply a recursive
2617  // lexicographic-distance-based
2618  // algorithm to infer the desired matrix
2619  std::vector<std::vector<double>> saliencesInit(saliences);
2620  std::vector<std::vector<double>> saliencesTmp(
2621  nbrRegions, vector<double>(nbrRegions, 0.0));
2622 
2623  // recursive max-min term-by-term product of the persistences matrix
2624  // until equilibrium state
2625  while (saliences != saliencesTmp) {
2626  saliencesTmp = saliences;
2627  for (size_t i = 0; i < nbrRegions; i++) {
2628  for (size_t j = i; j < nbrRegions; j++) {
2629  std::vector<double> minVec(nbrRegions, 0.0);
2630  for (size_t k = 0; k < nbrRegions; k++) {
2631  minVec[k] = min(saliences[i][k], saliencesInit[k][j]);
2632  }
2633  saliences[i][j] = *max_element(minVec.begin(), minVec.end());
2634  saliences[j][i] = saliences[i][j];
2635  }
2636  }
2637  }
2638  } // end if "full" && "life"
2639  else {
2640  cerr << "getSaliences(string typeOfMatrix = ''sparse'',momentOfLife = "
2641  "''death'') -> please choose typeOfMatrix"
2642  << " in the following : sparse, full; and please choose "
2643  "momentOfLife in the following : ''life'', ''death'' "
2644  << endl;
2645  }
2646 
2647  return saliences;
2648  } // end getSaliences
2649 
2658  {
2659  double dasguptaCF = 0;
2660 
2661  // sort nodes by decreasing internalNodeValuationInitial
2662  this->sortNodes(true);
2663 
2664  // We extract pairs of labels in the form of labelsin and labelsout
2665  // vectors
2666  std::vector<Edge<NodeT, ValGraphT>> edges = completeGraph.getEdges();
2667 
2668  std::vector<NodeT> labelsin = std::vector<NodeT>(edges.size(), 0);
2669  std::vector<NodeT> labelsout = std::vector<NodeT>(edges.size(), 0);
2670  std::vector<ValGraphT> weights = std::vector<ValGraphT>(edges.size(), 0);
2671 
2672  for (MarkerLabelT i = 0; i < edges.size(); i++) {
2673  labelsin.at(i) = edges.at(i).source;
2674  labelsout.at(i) = edges.at(i).source;
2675  weights.at(i) = edges.at(i).source;
2676  }
2677 
2678  // dendrogram traversal to compute Dasgupta cost function
2679  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2680  DendroNodeType &curNode = *dendroNodes.at(i);
2681  // we check if it is an internal node
2682  if (curNode.getIsInternalNode() == 1) {
2683  std::vector<MarkerLabelT> lookupProgeny = curNode.getLookupProgeny();
2684  std::vector<MarkerLabelT> lookupProgenyChildLeft =
2685  curNode.getChildLeft()->getLookupProgeny();
2686  std::vector<MarkerLabelT> lookupProgenyChildRight =
2687  curNode.getChildRight()->getLookupProgeny();
2688 
2689  double nbLeavesUnder =
2690  std::accumulate(lookupProgeny.begin(), lookupProgeny.end(), 0.);
2691 
2692  for (MarkerLabelT s = 0; s < edges.size(); s++) {
2693  MarkerLabelT labelIn = labelsin.at(s);
2694  MarkerLabelT labelOut = labelsout.at(s);
2695  double currentWeight = weights.at(s);
2696 
2697  if (lookupProgeny.at(labelIn) == lookupProgeny.at(labelOut) &&
2698  lookupProgenyChildLeft.at(labelIn) !=
2699  lookupProgenyChildRight.at(labelOut)) {
2700  dasguptaCF += nbLeavesUnder * currentWeight;
2701  }
2702  }
2703 
2704  // for (MarkerLabelT s = 0; s < lookupProgeny.size(); s++) {
2705  // for (MarkerLabelT t = s+1; t < lookupProgeny.size(); t++) {
2706  // if (lookupProgeny.at(s)==lookupProgeny.at(t) &&
2707  // lookupProgenyChildLeft.at(s)!=lookupProgenyChildRight.at(t)){
2708  // dasguptaCF += nbLeavesUnder*curVal;
2709  // }
2710  // }
2711  // }
2712  }
2713  }
2714  return dasguptaCF;
2715  }
2716 
2717  double computeScoreInteractiveSeg()
2718  {
2719  double score = 0;
2720 
2721  std::map<MarkerLabelT, std::vector<MarkerLabelT>> lookupMapMarkersLabels;
2722  std::map<std::pair<MarkerLabelT, MarkerLabelT>, std::vector<double>>
2723  lookupMapMarkersPairsValuations;
2724  std::map<std::pair<MarkerLabelT, MarkerLabelT>, std::vector<double>>
2725  lookupMapSameMarkersPairsValuations;
2726 
2727  // first dendrogram traversal to get marked nodes by labels in
2728  // lookupMapMarkers
2729  for (MarkerLabelT i = 0; i < nbrNodes; i++) {
2730  DendroNodeType &curNode = *dendroNodes.at(i);
2731  if (curNode.getMarker() != 0) {
2732  MarkerLabelT markerValue = curNode.getMarker();
2733  MarkerLabelT labelValue = curNode.getLabel();
2734  lookupMapMarkersLabels[markerValue].push_back(labelValue);
2735  }
2736  }
2737  // lookupMapMarkersLabels traversal to initialize
2738  // lookupMapMarkersPairsValuations
2739  for (typename std::map<MarkerLabelT, std::vector<MarkerLabelT>>::iterator
2740  iter = lookupMapMarkersLabels.begin();
2741  iter != lookupMapMarkersLabels.end(); ++iter) {
2742  typename std::map<MarkerLabelT, std::vector<MarkerLabelT>>::iterator
2743  iter2 = iter;
2744  iter2++;
2745  for (; iter2 != lookupMapMarkersLabels.end(); ++iter2) {
2746  if (iter->first != iter2->first) {
2747  std::pair<MarkerLabelT, MarkerLabelT> pairMarkers =
2748  make_pair(iter->first, iter2->first);
2749  lookupMapMarkersPairsValuations[pairMarkers] =
2750  std::vector<double>();
2751  }
2752  }
2753  std::pair<MarkerLabelT, MarkerLabelT> pairSameMarkers =
2754  make_pair(iter->first, iter->first);
2755  lookupMapSameMarkersPairsValuations[pairSameMarkers] =
2756  std::vector<double>();
2757  }
2758 
2759  // lookupMapMarkersPairsValuations traversal to complete the valuations
2760  // field
2761  for (typename std::map<std::pair<MarkerLabelT, MarkerLabelT>,
2762  std::vector<double>>::iterator iter =
2763  lookupMapMarkersPairsValuations.begin();
2764  iter != lookupMapMarkersPairsValuations.end(); iter++) {
2765  MarkerLabelT marker1 = iter->first.first;
2766  MarkerLabelT marker2 = iter->first.second;
2767  std::vector<MarkerLabelT> labelsMarker1 =
2768  lookupMapMarkersLabels[marker1];
2769  std::vector<MarkerLabelT> labelsMarker2 =
2770  lookupMapMarkersLabels[marker2];
2771 
2772  for (size_t i = 0; i < labelsMarker1.size(); i++) {
2773  for (size_t j = 0; j < labelsMarker2.size(); j++) {
2774  double scoreTmp = 0;
2775  MarkerLabelT label1 = labelsMarker1[i];
2776  MarkerLabelT label2 = labelsMarker2[j];
2777 
2778  for (size_t d = 0; d < nbrNodes; d++) {
2779  DendroNodeType & curNode = *dendroNodes.at(d);
2780  std::vector<MarkerLabelT> lookupProgenyd =
2781  curNode.getLookupProgeny();
2782  if (lookupProgenyd[label1] != 0 && lookupProgenyd[label2] != 0) {
2783  if (scoreTmp == 0) {
2784  scoreTmp = curNode.getInternalNodeValuationInitial();
2785  } else {
2786  if (curNode.getInternalNodeValuationInitial() < scoreTmp) {
2787  scoreTmp = curNode.getInternalNodeValuationInitial();
2788  }
2789  }
2790  }
2791  }
2792  if (scoreTmp != 0) {
2793  (iter->second).push_back(scoreTmp);
2794  }
2795  } // end for (int j =0; j<labelsMarker2.size();j++){
2796  } // end for (int i = 0;i<labelsMarker1.size();i++){
2797  }
2798 
2799  // lookupMapSameMarkersPairsValuations traversal to complete the
2800  // valuations field
2801  for (typename std::map<std::pair<MarkerLabelT, MarkerLabelT>,
2802  std::vector<double>>::iterator iter =
2803  lookupMapSameMarkersPairsValuations.begin();
2804  iter != lookupMapSameMarkersPairsValuations.end(); iter++) {
2805  MarkerLabelT marker = iter->first.first;
2806  std::vector<MarkerLabelT> labelsMarker = lookupMapMarkersLabels[marker];
2807  for (size_t i = 0; i < labelsMarker.size(); i++) {
2808  for (size_t j = 0; j < labelsMarker.size(); j++) {
2809  double scoreTmp = 0;
2810  MarkerLabelT label1 = labelsMarker[i];
2811  MarkerLabelT label2 = labelsMarker[j];
2812  for (size_t d = 0; d < nbrNodes; d++) {
2813  DendroNodeType & curNode = *dendroNodes.at(d);
2814  std::vector<MarkerLabelT> lookupProgenyd =
2815  curNode.getLookupProgeny();
2816  if (lookupProgenyd[label1] != 0 && lookupProgenyd[label2] != 0) {
2817  if (scoreTmp == 0) {
2818  scoreTmp = curNode.getInternalNodeValuationInitial();
2819  } else {
2820  if (curNode.getInternalNodeValuationInitial() < scoreTmp) {
2821  scoreTmp = curNode.getInternalNodeValuationInitial();
2822  }
2823  }
2824  }
2825  }
2826  (iter->second).push_back(scoreTmp);
2827  } // end for (int j =0; j<labelsMarker2.size();j++){
2828  } // end for (int i = 0;i<labelsMarker1.size();i++){
2829  }
2830 
2831  // check if lookupMapMarkersPairsValuations is correctly filled
2832  // for (typename
2833  // std::map<std::pair<MarkerLabelT,MarkerLabelT>,std::vector<double>
2834  // >::iterator iter = lookupMapMarkersPairsValuations.begin();
2835  // iter!=lookupMapMarkersPairsValuations.end();iter++){
2836  // cout << "m1 = " << iter->first.first << " m2 = " <<
2837  // iter->first.second << endl;
2838  // for (int i = 0; i<(iter->second).size();i++){
2839  // cout << (iter->second)[i]<<endl;
2840  // }
2841  // }
2842 
2843  // computation of the score
2844  // Mean and variances for pairs of different markers
2845  double meanDiffMarkers = 0;
2846  double varDiffMarkers = 0;
2847  double sizeDiffMarkers = 0;
2848  for (typename std::map<std::pair<MarkerLabelT, MarkerLabelT>,
2849  std::vector<double>>::iterator iter =
2850  lookupMapMarkersPairsValuations.begin();
2851  iter != lookupMapMarkersPairsValuations.end(); iter++) {
2852  sizeDiffMarkers = sizeDiffMarkers + 1;
2853  double mean_ind = 0;
2854  double var_ind = 0;
2855  for (MarkerLabelT i = 0; i < (iter->second).size(); i++) {
2856  mean_ind = mean_ind + (iter->second)[i];
2857  }
2858  mean_ind = mean_ind / (iter->second).size();
2859  for (MarkerLabelT i = 0; i < (iter->second).size(); i++) {
2860  var_ind = var_ind + (mean_ind - (iter->second)[i]) *
2861  (mean_ind - (iter->second)[i]);
2862  }
2863  var_ind = var_ind / (iter->second).size();
2864  meanDiffMarkers = meanDiffMarkers + mean_ind;
2865  varDiffMarkers = varDiffMarkers + var_ind;
2866  }
2867  meanDiffMarkers = meanDiffMarkers / sizeDiffMarkers;
2868  varDiffMarkers = varDiffMarkers / sizeDiffMarkers;
2869 
2870  // mean and variances for pairs of same markers
2871  double meanSameMarkers = 0;
2872  double varSameMarkers = 0;
2873  double sizeSameMarkers = 0;
2874  for (typename std::map<std::pair<MarkerLabelT, MarkerLabelT>,
2875  std::vector<double>>::iterator iter =
2876  lookupMapSameMarkersPairsValuations.begin();
2877  iter != lookupMapSameMarkersPairsValuations.end(); iter++) {
2878  sizeSameMarkers = sizeSameMarkers + 1;
2879  double mean_ind = 0;
2880  double var_ind = 0;
2881  for (MarkerLabelT i = 0; i < (iter->second).size(); i++) {
2882  mean_ind = mean_ind + (iter->second)[i];
2883  }
2884  mean_ind = mean_ind / (iter->second).size();
2885  for (MarkerLabelT i = 0; i < (iter->second).size(); i++) {
2886  varSameMarkers = varSameMarkers + (mean_ind - (iter->second)[i]) *
2887  (mean_ind - (iter->second)[i]);
2888  }
2889  var_ind = var_ind / (iter->second).size();
2890  meanSameMarkers = meanSameMarkers + mean_ind;
2891  varSameMarkers = varSameMarkers + var_ind;
2892  }
2893  meanSameMarkers = meanSameMarkers / sizeSameMarkers;
2894  varSameMarkers = varSameMarkers / sizeSameMarkers;
2895 
2896  score = (varSameMarkers * varSameMarkers) /
2897  (std::abs(meanSameMarkers - meanDiffMarkers));
2898  return score;
2899  }; // end computeScoreInteractiveSeg
2900 
2901  }; // end Dendrogram
2902 
2903  // TODO: debug distDendro (problem with wrapping in Python using SWIG -
2904  // segfault)
2905  // template <class MarkerLabelT, class NodeT,class ValGraphT>
2906  // double distDendro(Dendrogram<MarkerLabelT,NodeT,ValGraphT> &
2907  // d1,Dendrogram<MarkerLabelT,NodeT,ValGraphT> & d2,int order){
2908  // double distToReturn = 0;
2909  // if (d1.getNbrNodes()==d2.getNbrNodes()){
2910  // int nbrNodes = d1.getNbrNodes();
2911  // std::vector<DendroNode<MarkerLabelT>*>&dendroNodes1 =
2912  // d1.getDendroNodes();
2913  // for (int i = 0; i<nbrNodes; i++){
2914  // DendroNode<MarkerLabelT> &node1 = *dendroNodes1[i];
2915  // MarkerLabelT researchedLabel = node1.getLabel();
2916  // DendroNode<MarkerLabelT> &node2 =
2917  // *d2.researchLabel(researchedLabel);
2918  // distToReturn = distToReturn
2919  // +
2920  // std::abs(std::pow(node1.getInternalNodeValuationInitial()-node2.getInternalNodeValuationInitial(),order));
2921  // }
2922  // }
2923  // else{
2924  // cout << "The two dendrograms do not have the same number of nodes" <<
2925  // endl;
2926  // }
2927  // return distToReturn;
2928  // };
2929 
2930  // template <class MarkerLabelT, class NodeT,class ValGraphT>
2931  // Dendrogram<MarkerLabelT,NodeT,ValGraphT>
2932  // infDendro(Dendrogram<MarkerLabelT,NodeT,ValGraphT> &
2933  // d1,Dendrogram<MarkerLabelT,NodeT,ValGraphT> & d2){
2934  // Dendrogram<MarkerLabelT,NodeT,ValGraphT> dendroToReturn(d1);
2935  // if (d1.getNbrNodes()==d2.getNbrNodes()){
2936  // int nbrNodes = d1.getNbrNodes();
2937  // std::vector<DendroNode<MarkerLabelT>*>&dendroNodes1 =
2938  // d1.getDendroNodes();
2939  // for (int i = 0; i<nbrNodes; i++){
2940  // DendroNode<MarkerLabelT> &node1 = *dendroNodes1[i];
2941  // MarkerLabelT researchedLabel = node1.getLabel();
2942  // DendroNode<MarkerLabelT> &node2 =
2943  // *d2.researchLabel(researchedLabel);
2944  // DendroNode<MarkerLabelT> &nodeToModify =
2945  // *dendroToReturn.researchLabel(researchedLabel);
2946  // nodeToModify.setInternalNodeValuationInitial(std::min(node1.getInternalNodeValuationInitial(),
2947  // node2.getInternalNodeValuationInitial()));
2948  // }
2949  // }
2950  // return dendroToReturn;
2951  // }
2952 
2953  // template <class MarkerLabelT, class NodeT,class ValGraphT>
2954  // void supDendro(Dendrogram<MarkerLabelT,NodeT,ValGraphT> &
2955  // d1,Dendrogram<MarkerLabelT,NodeT,ValGraphT> & d2,
2956  // Dendrogram<MarkerLabelT,NodeT,ValGraphT> & dOut){
2957  // dOut = d1.clone();
2958  // if (dOut.getNbrNodes()==d2.getNbrNodes()){
2959  // int nbrNodes = dOut.getNbrNodes();
2960  // std::vector<DendroNode<MarkerLabelT>*>&dendroNodesOut =
2961  // dOut.getDendroNodes();
2962  // for (int i = 0; i<nbrNodes; i++){
2963  // DendroNode<MarkerLabelT> &node1 = *dendroNodesOut[i];
2964  // MarkerLabelT researchedLabel = node1.getLabel();
2965  // DendroNode<MarkerLabelT> &node2 =
2966  // *d2.researchLabel(researchedLabel);
2967  // DendroNode<MarkerLabelT> &nodeToModify =
2968  // *dOut.researchLabel(researchedLabel);
2969  // nodeToModify.setInternalNodeValuationInitial(std::max(node1.getInternalNodeValuationInitial(),
2970  // node2.getInternalNodeValuationInitial()));
2971  // };
2972  // }
2973  // // dOut.reorganize();
2974  // }
2975  //
2976  // template <class MarkerLabelT, class NodeT,class ValGraphT>
2977  // void infDendro(Dendrogram<MarkerLabelT,NodeT,ValGraphT> &
2978  // d1,Dendrogram<MarkerLabelT,NodeT,ValGraphT> & d2,
2979  // Dendrogram<MarkerLabelT,NodeT,ValGraphT> & dOut){
2980  // dOut = d1.clone();
2981  // if (dOut.getNbrNodes()==d2.getNbrNodes()){
2982  // int nbrNodes = dOut.getNbrNodes();
2983  // std::vector<DendroNode<MarkerLabelT>*>&dendroNodesOut =
2984  // dOut.getDendroNodes();
2985  // for (int i = 0; i<nbrNodes; i++){
2986  // DendroNode<MarkerLabelT> &node1 = *dendroNodesOut[i];
2987  // MarkerLabelT researchedLabel = node1.getLabel();
2988  // DendroNode<MarkerLabelT> &node2 =
2989  // *d2.researchLabel(researchedLabel);
2990  // DendroNode<MarkerLabelT> &nodeToModify =
2991  // *dOut.researchLabel(researchedLabel);
2992  // nodeToModify.setInternalNodeValuationInitial(std::min(node1.getInternalNodeValuationInitial(),
2993  // node2.getInternalNodeValuationInitial()));
2994  // };
2995  // }
2996  // // dOut.reorganize();
2997  // }
2998 
3008  template <class MarkerLabelT, class NodeT, class ValGraphT>
3009  void
3010  getLevelHierarchySeg(const smil::Image<MarkerLabelT> & imLabels,
3011  Dendrogram<MarkerLabelT, NodeT, ValGraphT> &dendrogram,
3012  Graph<NodeT, ValGraphT> & associated_mst,
3013  double desiredNbrRegions,
3015  {
3016  Dendrogram<MarkerLabelT, NodeT, ValGraphT> dendrobis = dendrogram.clone();
3017  Graph<NodeT, ValGraphT> mstbis = associated_mst.clone();
3018  dendrobis.sortNodes(true);
3019 
3020  MarkerLabelT desiredNbrRegionsbis = desiredNbrRegions;
3021  if (desiredNbrRegions <= 1) {
3022  desiredNbrRegionsbis =
3023  (dendrobis.getNbrNodes() / 2 + 1) * desiredNbrRegions;
3024  }
3025  if ((desiredNbrRegionsbis > 1) &&
3026  (desiredNbrRegionsbis <= (dendrobis.getNbrNodes() / 2 + 1))) {
3027  double lamb = dendrobis.getNodeValue(desiredNbrRegionsbis - 2,
3028  "internalNodeValuationInitial");
3029  dendrobis.removeMSTEdgesDendrogram(mstbis, lamb);
3030  graphToMosaic(imLabels, mstbis, imOut);
3031  } else {
3032  cerr << "getLevelHierarchySeg(imFinePartition,dendro,associated_mst,"
3033  "desiredNbrRegions,imOut) "
3034  << "-> the desiredNbrRegions specified is higher than the number of "
3035  "regions"
3036  << endl;
3037  }
3038  }
3048  template <class MarkerLabelT, class NodeT, class ValGraphT>
3049  void
3050  getThreshHierarchySeg(const smil::Image<MarkerLabelT> & imLabels,
3051  Dendrogram<MarkerLabelT, NodeT, ValGraphT> &dendrogram,
3052  Graph<NodeT, ValGraphT> &associated_mst, double thresh,
3054  {
3055  if (thresh >= 0) {
3056  Dendrogram<MarkerLabelT, NodeT, ValGraphT> dendrobis = dendrogram.clone();
3057  Graph<NodeT, ValGraphT> mstbis = associated_mst.clone();
3058  dendrobis.sortNodes(true);
3059  dendrobis.removeMSTEdgesDendrogram(mstbis, thresh);
3060  graphToMosaic(imLabels, mstbis, imOut);
3061  } else {
3062  cout << "getLevelHierarchySeg(imFinePartition,dendro,associated_mst,"
3063  "thresh,"
3064  "imOut) "
3065  << "-> the thresh specified must be higher than 0" << endl;
3066  }
3067  }
3068 
3069 } // namespace smil
3070 
3071 #endif // _DENDROGRAM_HPP
DendroNode : node of a Dendrogram.
Definition: DendroNode.hpp:57
double getInternalNodeValuationInitial()
Setters and getters.
Definition: DendroNode.hpp:196
Dendrogram.
Definition: Dendrogram.hpp:49
Dendrogram(const Dendrogram &dendrogramToCopy)
Copy constructor.
Definition: Dendrogram.hpp:65
Dendrogram(Graph< NodeT, ValGraphT > &mst)
Dendrogram constructor from a MST graph.
Definition: Dendrogram.hpp:126
std::vector< MarkerLabelT > getLookupProgeny(size_t nodeIndex, string nameOfValueWanted)
Get the lookup of the progeny of a node.
Definition: Dendrogram.hpp:2114
Dendrogram()
Default constructor.
Definition: Dendrogram.hpp:63
void sortReverseNodes()
Reorganizes dendroNodes by decreasing valuationInitial.
Definition: Dendrogram.hpp:328
void normalize(const std::string typeOfNormalization="reg")
Dendrogram ultrametric values normalization.
Definition: Dendrogram.hpp:731
void removeMSTEdgesDendrogram(Graph< NodeT, ValGraphT > &associated_mst, double lbd)
Given an internal node index lambda of the dendrogram, remove corresponding edge in the associated MS...
Definition: Dendrogram.hpp:811
double getNodeValue(size_t nodeIndex, string nameOfValueWanted)
Access a value of a node in the dendrogram.
Definition: Dendrogram.hpp:2027
std::vector< std::vector< double > > getSaliences(string typeOfMatrix="sparse", string momentOfLife="death")
Compute saliences matrix (salience expresses ultrametric values)
Definition: Dendrogram.hpp:2432
std::vector< MarkerLabelT > getThreshLookupProgeny(double thresh, string nameOfValueWanted)
Get the lookup of the progeny below a certain ultrametric threshold.
Definition: Dendrogram.hpp:2141
void HierarchicalConstruction(const std::string typeOfHierarchy, const int nParam=50, const smil::Image< MarkerLabelT > &imMosa=smil::Image< MarkerLabelT >(), const std::string typeOfTransform="erode", const StrElt &se=DEFAULT_SE)
Computes a new hierarchy from a given dendrogram hierarchy.
Definition: Dendrogram.hpp:844
virtual ~Dendrogram()
Destructor.
Definition: Dendrogram.hpp:306
void setNodeValue(size_t nodeIndex, string nameOfValueWanted, double value)
Manually modify a value of a node in the dendrogram.
Definition: Dendrogram.hpp:2074
std::vector< std::vector< double > > getPersistences(string typeOfMatrix="sparse", string momentOfLife="death")
Compute persistences matrix (persistence measures only the order of fusion)
Definition: Dendrogram.hpp:2200
void EnergyConstruction(const std::string typeOfHierarchy)
Computes a new hierarchy from a given dendrogram hierarchy.
Definition: Dendrogram.hpp:1798
void putValuationsFinalInInitial()
Put the internalNodeValuationFinal in internalNodeValuationInitial and then set internalNodeValuation...
Definition: Dendrogram.hpp:785
double computeDasguptaCF(Graph< NodeT, ValGraphT > completeGraph)
Compute dasgupta score (paper "A cost function for similarity-based hierarchical clustering")
Definition: Dendrogram.hpp:2657
void sortNodes(bool reverse=false)
Reorganizes dendroNodes...
Definition: Dendrogram.hpp:318
Dendrogram(Graph< NodeT, ValGraphT > &mst, smil::Image< MarkerLabelT > &imLabels, smil::Image< MarkerLabelT > &imIn, const size_t nbOfMoments=5)
Constructor from a MST graph and labels/image to compute the moments.
Definition: Dendrogram.hpp:217
Dendrogram clone()
Clone Dendrogram.
Definition: Dendrogram.hpp:313
std::vector< DendroNodeType * > & getDendroNodes()
Setters and Getters.
Definition: Dendrogram.hpp:1986
Non-oriented graph.
Definition: DGraph.hpp:197
size_t getEdgeNbr()
getEdgeNbr() -
Definition: DGraph.hpp:371
size_t getNodeNbr()
getNodeNbr() -
Definition: DGraph.hpp:363
void removeEdge(const size_t index)
Remove an edge.
Definition: DGraph.hpp:425
EdgeListType & getEdges()
getEdges() - Get a vector containing the graph edges
Definition: DGraph.hpp:534
void sortEdges(bool reverse=false)
Sort edges (by weight as defined by the operator < of class Edge)
Definition: DGraph.hpp:334
NodeValuesType & getNodeValues()
getNodeValues() -
Definition: DGraph.hpp:542
Main Image class.
Definition: DImage.hpp:57
Base structuring element.
Definition: DStructuringElement.h:70
RES_T mul(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
mul() - Multiply two images
Definition: DImageArith.hpp:749
RES_T sub(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
sub() - Subtraction between two images
Definition: DImageArith.hpp:160
RES_T pow(const Image< T1 > &imIn, Image< T2 > &imOut, double exponent=2)
pow() - power of an image
Definition: DImageArith.hpp:905
RES_T compare(const Image< T1 > &imIn1, const char *compareType, const Image< T1 > &imIn2, const Image< T2 > &trueIm, const Image< T2 > &falseIm, Image< T2 > &imOut)
compare() - Compare each pixel of two images (or an image and a value)
Definition: DImageArith.hpp:1244
RES_T applyLookup(const Image< T1 > &imIn, const mapT &_map, Image< T2 > &imOut, T2 defaultValue=T2(0))
applyLookup() - Apply a lookup map to a labeled image
Definition: DImageArith.hpp:1990
RES_T open(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
open() - Morphological grayscale opening
Definition: DMorphoFilter.hpp:123
RES_T close(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
close() - Morphological grayscale closing
Definition: DMorphoFilter.hpp:70
size_t area(const Image< T > &imIn)
area() - Area of an image
Definition: DMeasures.hpp:91
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
maxVal() - Max value of an image
Definition: DMeasures.hpp:348
RES_T erode(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE, const T borderVal=ImDtTypes< T >::max())
erode() - Morphological grayscale erosion
Definition: DMorphoBase.hpp:112
RES_T dilate(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE, const T borderVal=ImDtTypes< T >::min())
dilate() - Morphological grayscale dilation
Definition: DMorphoBase.hpp:65