SMIL  0.9.1
DMorphoMaxTree.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'' AND ANY
18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE LIABLE FOR ANY
21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
28 
29 
30 #ifndef _D_MAX_TREE
31 #define _D_MAX_TREE
32 
33 #include "DMorphImageOperations.hpp"//BMI
34 #include "DMorphoHierarQ.hpp"//BMI
35 #include "Core/include/DImage.h"
36 #include "Base/include/private/DImageHistogram.hpp"
37 #include "Morpho/include/private/DMorphoMaxTreeCriteria.hpp"
38 
39 #include <complex>
40 #include <math.h>
41 namespace smil
42 {
50 #ifndef SWIG
51 
52 
53 #define COLUMN_NBR 20// 16
54 #define COLUMN_SIZE 2097152//131072
55 
56 #define COLUMN_SHIFT (COLUMN_NBR+1)
57 #define COLUMN_MOD (COLUMN_SIZE-1)
58 
59 
60 #define GET_TREE_OBJ(type, node) type[node >> COLUMN_SHIFT][node & COLUMN_MOD]
61 
62 
63  // BEGIN BMI
64 
65 
66 template <class T, class CriterionT, class OffsetT=size_t, class LabelT = UINT32>
67 class MaxTree2
68 {
69 private:
70  size_t GRAY_LEVEL_NBR;
71  Image<T> const *img;
72 
73  // PriorityQueue<T, OffsetT> pq;
74  // HierarchicalQueue<T,size_t, STD_Queue<size_t> > &hq;
76 
77  T **levels;
78 
79  LabelT **children;
80  LabelT **brothers;
81  CriterionT **criteria;
82 
83  LabelT *labels, curLabel;
84 
85  bool initialized;
86 
87  size_t imWidth ;
88  size_t imHeight ;
89 
90 
91  void reset()
92  {
93  if (!initialized)
94  return;
95 
96  for (OffsetT i=0;i<COLUMN_NBR;i++)
97  {
98  if (!levels[i])
99  break;
100  delete[] levels[i]; levels[i] = NULL;
101  delete[] children[i]; children[i] = NULL;
102  delete[] brothers[i]; brothers[i] = NULL;
103  delete[] criteria[i]; criteria[i] = NULL;
104  }
105  }
106 
107  void allocatePage(UINT page)
108  {
109  levels[page] = new T[COLUMN_SIZE]();
110  children[page] = new LabelT[COLUMN_SIZE]();
111  brothers[page] = new LabelT[COLUMN_SIZE]();
112  criteria[page] = new CriterionT[COLUMN_SIZE]();
113  }
114 
115  T initialize(const Image<T> &imIn, LabelT *img_eti, const StrElt &se)
116  {
117  imIn.getSize(imSize);
118 
119 
120  // BMI BEGIN
121  sePtsNbr = sePts.size();
122  dOffsets.clear();
123  sePts.clear();
124  oddSE = se.odd;
125 
126  // set an offset distance for each se point (!=0,0,0)
127  for(vector<IntPoint>::const_iterator it = se.points.begin() ; it!=se.points.end() ; it++)
128  if (it->x!=0 || it->y!=0 || it->z!=0)
129  {
130  sePts.push_back(*it);
131  dOffsets.push_back(it->x + it->y*imSize[0] + it->z*imSize[0]*imSize[1]);
132  }
133 
134  sePtsNbr = sePts.size();
135  // BMI END
136 
137  if (initialized)
138  reset();
139 
140  this->img = &imIn;
141  typename ImDtTypes<T>::lineType pix = img->getPixels();
142 
143  T minValue = ImDtTypes<T>::max();
144  T tMinV = ImDtTypes<T>::min();
145  OffsetT minOff = 0;
146  for (size_t i=0;i<img->getPixelCount();i++)
147  if (pix[i]<minValue)
148  {
149  minValue = pix[i];
150  minOff = i;
151  if (minValue==tMinV)
152  break;
153  }
154 
155  allocatePage(0);
156 
157  curLabel = 1;
158  levels[0][curLabel] = minValue;
159 
160  memset(labels, 0, GRAY_LEVEL_NBR*sizeof(LabelT));// BMI labels? quel rapport avec img_eti?
161 
162  img_eti[minOff] = curLabel;
163  labels[minValue] = curLabel;
164 
165  hq.initialize(*img);
166  hq.push(minValue, minOff);
167 
168 
169  size_t x, y, z;
170  img->getCoordsFromOffset(minOff, x, y, z);
171  // std::cout<<"PUSH:offset="<<minOff<<"(x,y)="<<x<<", "<<y<<", val="<<minValue<<"\n";
172 
173 
174  // getCriterion(curLabel).xmin = getCriterion(curLabel).xmax = x; // A voir comment on peut melanger des criteres... BMI
175  // getCriterion(curLabel).ymin = getCriterion(curLabel).ymax = y;
176  // getCriterion(curLabel).zmin = getCriterion(curLabel).zmax = z;
177 
178  getCriterion(curLabel).initialize();
179  curLabel++;
180 
181  initialized = true;
182 
183  imWidth = img->getWidth();
184  imHeight = img->getHeight();
185 
186  return minValue;
187  }
188 
189  int nextLowerLabel(T valeur)
190  {
191  if ((curLabel & COLUMN_MOD) == 0)
192  allocatePage(curLabel >> COLUMN_SHIFT);
193 
194  getLevel(curLabel) = valeur;
195  int i;
196  for(i=valeur-1;labels[i]==0;i--);
197 
198  getChild(curLabel) = getChild(labels[i]);
199  getChild(labels[i]) = curLabel;
200  getBrother(curLabel) = getBrother(getChild(curLabel));
201  getBrother(getChild(curLabel)) = 0;
202  getCriterion(curLabel).reset();
203  return curLabel++;
204  }
205 
206  int nextHigherLabel(T parent_valeur, T valeur)
207  {
208  if ((curLabel & COLUMN_MOD) == 0)
209  allocatePage(curLabel >> COLUMN_SHIFT);
210 
211  getLevel(curLabel) = valeur;
212  getBrother(curLabel) = getChild(labels[parent_valeur]);
213  getChild(labels[parent_valeur]) = curLabel;
214  getCriterion(curLabel).reset();
215  return curLabel++;
216  }
217 
218  bool subFlood(typename ImDtTypes<T>::lineType imgPix, LabelT *img_eti, OffsetT p, OffsetT p_suiv)
219  {
220  LabelT indice;
221 
222  if (imgPix[p_suiv]>imgPix[p])
223  {
224  LabelT j;
225  for(j=imgPix[p]+1;j<imgPix[p_suiv];j++)
226  labels[j]=0;
227  indice = img_eti[p_suiv] = labels[j] = nextHigherLabel(imgPix[p], imgPix[p_suiv]);
228 
229  }
230  else if (labels[imgPix[p_suiv]]==0)
231  indice = img_eti[p_suiv] = labels[imgPix[p_suiv]] = nextLowerLabel(imgPix[p_suiv]);
232  else
233  indice = img_eti[p_suiv] = labels[imgPix[p_suiv]];
234 
235  size_t x, y, z;
236  img->getCoordsFromOffset(p_suiv, x, y, z);
237 
238  getCriterion(indice).update(x,y,z);
239  hq.push(imgPix[p_suiv], p_suiv);
240 
241  // std::cout<<"PUSH:offset="<<p_suiv<<"("<<x<<","<<y<<"), val="<<int(imgPix[p_suiv])<<"\n";
242 
243  if (imgPix[p_suiv]>imgPix[p])
244  {
245  hq.push(imgPix[p], p);
246  // std::cout<<"PUSH_P:offset="<<p<<", val="<<int(imgPix[p])<<"\n";
247  return true;
248  }
249  return false;
250  }
251 
252  void flood(const Image<T> &img, UINT *img_eti, unsigned int level)
253  {
254  OffsetT p;
255 
256  typename ImDtTypes<T>::lineType imgPix = img.getPixels();
257  size_t x0, y0, z0;
258 
259 
260  while( (hq.getHigherLevel()>=level) && !hq.isEmpty())
261  {
262  p = hq.pop();
263 
264  img.getCoordsFromOffset(p, x0, y0, z0);
265  // std::cout<<"POP:offset="<<p<<"(x,y)="<<x0<<", "<<y0<<"\n";
266  // std::cout<<"-------------------------\n";
267 
268 
269  bool oddLine = oddSE && ((y0)%2);
270 
271  int x, y, z; // not size_t in order to (possibly be negative!)
272  OffsetT p_suiv;
273 
274  for(UINT i=0;i<sePtsNbr;i++)
275  {
276  IntPoint &pt = sePts[i];
277  x = x0 + pt.x;
278  y = y0 + pt.y;
279  z = z0 + pt.z;
280 
281  if (oddLine)
282  x += (((y+1)%2)!=0);
283 
284  if (x>=0 && x<(int)imSize[0] && y>=0 && y<(int)imSize[1] && z>=0 && z<(int)imSize[2])
285  {
286  p_suiv = p + dOffsets[i];
287  if (oddLine)
288  p_suiv += (((y+1)%2)!=0);
289 
290  if(img_eti[p_suiv]==0){
291  // std::cout<<"("<<x0<<","<<y0<<"),"<<"("<<x<<","<<y<<"),"<<p<<","<<p_suiv<<"\n";
292  if(subFlood(imgPix, img_eti, p, p_suiv))
293  break;
294  }
295  }
296  }// for each ngb
297 
298  }// while hq.notEmpty
299  }//void flood
300 
301 public:
302  MaxTree2():hq(true)// TOTO BMI
303  {
304  GRAY_LEVEL_NBR = ImDtTypes<T>::max()-ImDtTypes<T>::min()+1;
305  // hq.reverse();// change priority order (max first)
306  // hq = HierarchicalQueue<T,OffsetT> (true);
307 
308  children = new LabelT*[COLUMN_NBR]();
309  brothers = new LabelT*[COLUMN_NBR]();
310  levels = new T*[COLUMN_NBR]();
311  criteria = new CriterionT*[COLUMN_NBR]();
312  labels = new LabelT[GRAY_LEVEL_NBR];
313 
314  initialized = false;
315 
316  }
317  ~MaxTree2()
318  {
319  reset();
320 
321  delete[] children;
322  delete[] brothers;
323  delete[] levels;
324  delete[] criteria;
325  delete[] labels;
326  }
327 protected:
328  size_t imSize[3];
329 
330  vector<IntPoint> sePts;
331  UINT sePtsNbr;
332  bool oddSE;
333  vector<int> dOffsets;
334 
335 public:
336 
337  inline CriterionT &getCriterion(const LabelT node)
338  {
339  return GET_TREE_OBJ(criteria, node);
340  }
341 
342  inline T &getLevel(const LabelT node)
343  {
344  return GET_TREE_OBJ(levels, node);
345  }
346 
347  inline LabelT &getChild(const LabelT node)
348  {
349  return GET_TREE_OBJ(children, node);
350  }
351 
352  inline LabelT &getBrother(const LabelT node)
353  {
354  return GET_TREE_OBJ(brothers, node);
355  }
356 
357  inline LabelT getLabelMax()
358  {
359  return curLabel;
360  }
361 
362  inline int getImWidth()
363  {
364  return imWidth;
365  }
366 
367  inline int getImHeight()
368  {
369  return imHeight;
370  }
371 
372  int build(const Image<T> &img, LabelT *img_eti, const StrElt &se)
373  {
374  T minValue = initialize(img, img_eti,se);
375 
376  flood(img, img_eti, minValue);// BMI: dOffset already contains se information
377  updateCriteria(labels[minValue]);
378  return labels[minValue];
379  }
380 
381  // BMI (from Andres)
383  void updateCriteria(const int node);
384 
385 
386  CriterionT PREVIOUS_OLD_FROMJONATHAN_updateCriteria(LabelT node)
387  {
388  LabelT child = getChild(node);
389  while (child!=0)
390  {
391  CriterionT c = updateCriteria(child);
392  getCriterion(node).ymin = MIN(getCriterion(node).ymin, c.ymin);
393  getCriterion(node).ymax = MAX(getCriterion(node).ymax, c.ymax);
394  getCriterion(node).xmin = MIN(getCriterion(node).xmin, c.xmin);
395  getCriterion(node).xmax = MAX(getCriterion(node).xmax, c.xmax);
396  getCriterion(node).merge(getCriterion(child));
397  child = getBrother(child);
398  }
399  return getCriterion(node);
400  }
401 
402 };
403 
404  // END BMI
405 
406 // Update criteria on a given max-tree node.// From Andres
407 template <class T, class CriterionT, class OffsetT, class LabelT>
409 {
410  LabelT child = getChild(node);
411  while ( child != 0 )
412  {
413  updateCriteria( child );
414  getCriterion( node ).merge(&getCriterion( child ) );
415  child = getBrother( child );
416  }
417 }
418 
419 
420 // NEW BMI # ##################################################
421 //(tree, transformee_node, indicatrice_node, child, stopSize, (T)0, 0, hauteur, tree.getLevel(root), tree.getLevel(root));
422 template <class T, class CriterionT, class OffsetT, class LabelT, class tAttType>
423 void ComputeDeltaUO(MaxTree2<T,CriterionT,OffsetT,LabelT> &tree, T* transformee_node, tAttType* indicatrice_node, int node, int nParent, T prev_residue, tAttType stop, UINT delta, int isPrevMaxT){
424 
425  //self,node = 1, nParent =0, stop=0, delta = 0, isPrevMaxT = 0):
426  int child; // index node
427  T current_residue;
428  UINT cNode, cParent; // attributes
429  T lNode, lParent; // node levels, the same type than input image
430 
431  cNode = tree.getCriterion(node).getAttributeValue();//ymax-tree.getCriterion(node).ymin+1;// #current criterion
432  lNode = tree.getLevel(node);// #current level
433 
434 
435  cParent = tree.getCriterion(nParent).getAttributeValue();//ymax-tree.getCriterion(nParent).ymin+1;// #current criterion
436  lParent = tree.getLevel(nParent);// #current level
437 
438  int flag;
439 
440  if ((cParent - cNode) <= delta){
441  flag = 1;
442  }
443  else{
444  flag = 0;
445  }
446  if (flag){
447  current_residue = prev_residue + lNode - lParent ;
448  }
449  else{
450  current_residue = lNode - lParent;
451  }
452 
453  transformee_node[node] = transformee_node[nParent];
454  indicatrice_node[node] = indicatrice_node[nParent];
455 
456  int isMaxT = 0;
457 
458  if(cNode < stop){
459  if (current_residue > transformee_node[node]){
460  // std::cout<<"UPDATE RES\n";
461  isMaxT = 1;
462  transformee_node[node] = current_residue;
463  if(! (isPrevMaxT && flag)){
464  indicatrice_node[node] = cNode + 1;
465  }
466 
467  }
468  }
469  else{
470  indicatrice_node[node] = 0;
471  }
472  child=tree.getChild(node);
473  while (child!=0){
474  ComputeDeltaUO(tree, transformee_node, indicatrice_node, child, node, current_residue, stop, delta,isMaxT);
475  child = tree.getBrother(child);
476  }
477 
478 }
479  template <class T1, class T2>
480  void compute_max(MaxTree2<T1,HeightCriterion,size_t,UINT32> &tree, T1* transformee_node, T2* indicatrice_node, UINT32 node, T2 stop, T1 max_tr, T2 max_in, T2 hauteur_parent, T1 valeur_parent, T1 previous_value)
481 {
482  T1 m;
483  T1 max_node;
484  T2 max_criterion;
485  UINT32 child;
486  T2 hauteur = tree.getCriterion(node).getAttributeValue();//ymax-tree.getCriterion(node).ymin+1;
487  // std::cout<<"IN COMPUTE_MAX:"<<"; node="<<node<<"\n";
488  m = (hauteur==hauteur_parent) ? tree.getLevel(node)-previous_value : tree.getLevel(node)-valeur_parent;
489  if (hauteur>=stop)
490  {
491  max_node = max_tr;
492  max_criterion = 0;//max_in;
493  transformee_node[node] = max_node;
494 
495  indicatrice_node[node]=0;
496  child=tree.getChild(node);
497  }
498  else
499  {
500  if (m>max_tr)
501  {
502  max_node=m;
503  max_criterion=hauteur;
504  } else
505  {
506  max_node=max_tr;
507  max_criterion=max_in;
508  }
509  transformee_node[node]=max_node;
510  indicatrice_node[node]=max_criterion+1;
511  child=tree.getChild(node);
512  }
513  if (hauteur==hauteur_parent)
514  {
515 
516  while (child!=0)
517  {
518  // if(child > tree.getLabelMax()){
519  // std::cout<<"ERROR call child:"<<child<<"\n";
520  // }
521  if (hauteur_parent>stop)
522  compute_max(tree, transformee_node, indicatrice_node, child, stop, max_node, max_criterion, hauteur, tree.getLevel(node), previous_value);
523  else
524  compute_max(tree, transformee_node, indicatrice_node, child, stop, max_node, max_criterion, hauteur, tree.getLevel(node)/*valeur_parent*/, previous_value);
525  child = tree.getBrother(child);
526  }
527  }
528  else
529  {
530  while (child!=0)
531  {
532  // if(child > tree.getLabelMax()){
533  // std::cout<<"ERROR call child:"<<child<<"\n";
534  // }
535 
536  compute_max(tree, transformee_node, indicatrice_node, child, stop, max_node, max_criterion, hauteur, tree.getLevel(node), valeur_parent);
537  child = tree.getBrother(child);
538  }
539  }
540 }
541 
542  template <class T1, class T2>
543 void compute_contrast(MaxTree2<T1,HeightCriterion,size_t,UINT32> &tree, T1* transformee_node, T2* indicatrice_node, UINT32 root, T2 stopSize,UINT delta = 0)
544 {
545 
546  UINT32 child;
547  T2 hauteur;
548  // std::cout<<"IN compute_contrast\n";
549 
550  transformee_node[root]=0;
551  indicatrice_node[root]=0;
552  // std::cout<<"IN compute_contrast 2\n";
553  // tree.updateCriteria(root);
554 
555  hauteur = tree.getCriterion(root).getAttributeValue();//ymax - tree.getCriterion(root).ymin+1;
556  child = tree.getChild(root);
557  if(delta == 0){
558  while (child!=0)
559  {
560  // std::cout<<"child"<<child<<"; level"<<tree.getLevel(child);
561 
562  compute_max(tree, transformee_node, indicatrice_node, child, stopSize, (T1)0, (T2)0, hauteur, tree.getLevel(root), tree.getLevel(root));
563  child = tree.getBrother(child);
564  }
565  }// END delta == 0
566  else{
567  while (child!=0)
568  {
569  // std::cout<<"child"<<child<<"; level"<<tree.getLevel(child);
570  ComputeDeltaUO(tree, transformee_node, indicatrice_node, child, root/*parent*/, (T1)0/* prev_residue*/, stopSize /*stop*/, delta, 0 /*isPrevMaxT*/);
571 
572  child = tree.getBrother(child);
573  }
574  }//END dela != 0
575 }
576 
577 
578 #endif // SWIG
579 
593  template <class T1, class T2>
594  RES_T ultimateOpen(const Image<T1> &imIn, Image<T1> &imTrans, Image<T2> &imIndic, const StrElt &se = DEFAULT_SE, T2 stopSize=0, UINT delta = 0)
595  {
596 
597  ASSERT_ALLOCATED(&imIn, &imTrans, &imIndic);
598  ASSERT_SAME_SIZE(&imIn, &imTrans, &imIndic);
599 
600  if (stopSize==0)
601  stopSize = imIn.getHeight()-1;
602 
603  int imSize = imIn.getPixelCount();
604  UINT *img_eti = new UINT[imSize]();
605 
606 
608  UINT32 root = tree.build(imIn, img_eti,se);
609 
610  std::cout<<"ULTIMATE OPEN, after tree.build"<<"NB VERTEX="<<tree.getLabelMax()<<"\n";
611  T1 *transformee_node = new T1[tree.getLabelMax()]();
612  T2 *indicatrice_node = new T2[tree.getLabelMax()]();
613  // std::cout<<"ULTIMATE OPEN, after memory allocation\n";
614  compute_contrast(tree, transformee_node, indicatrice_node, root, stopSize,delta);
615  // std::cout<<"ULTIMATE OPEN, compute_contrast\n";
616 
617  typename ImDtTypes<T1>::lineType transformeePix = imTrans.getPixels();
618  typename ImDtTypes<T2>::lineType indicatricePix = imIndic.getPixels();
619 
620  for(int i=0;i<imSize;i++)
621  {
622  transformeePix[i]=transformee_node[img_eti[i]];
623  indicatricePix[i]=indicatrice_node[img_eti[i]];
624  }
625 
626  delete[] img_eti;
627  delete[] transformee_node;
628  delete[] indicatrice_node;
629 
630  imTrans.modified();
631  imIndic.modified();
632 
633  return RES_OK;
634  }
635 
636 
637 
638 #ifndef SWIG
639 // NEW BMI # ##################################################
640 //(tree, transformee_node, indicatrice_node, child, stopSize, (T)0, 0, hauteur, tree.getLevel(root), tree.getLevel(root));
641 template <class T, class CriterionT, class OffsetT, class LabelT, class tAttType>
642 void ComputeDeltaUOMSER(MaxTree2<T,CriterionT,OffsetT,LabelT> &tree, T* transformee_node, tAttType* indicatrice_node, int node, int nParent, int first_ancestor, tAttType stop, UINT delta, UINT method, int isPrevMaxT,UINT minArea=0,T threshold=0,T mymax=0){
643  // method: 1 (MSER), 2 (RGR)
644 
645  // "node": the current node; "nParent": its direct parent (allows
646  // attribute comparison for Delta versions); "first_ancestor": the
647  // first ancestor with a different attribute (used for residue
648  // computation level(first_ancestor) - level(node) and for area stability computation
649 
650  //self,node = 1, nParent =0, stop=0, delta = 0, isPrevMaxT = 0):
651  int child; // index node
652  T current_residue, stab_residue;
653  UINT hNode, hParent; // attributes
654  size_t aNode, aParent,aAncestor;
655  T lNode, lParent, lAncestor; // node levels, the same type than input image
656  float stability;
657 
658 
659  hNode = tree.getCriterion(node).getAttributeValue().H;// #current criterion
660  aNode = tree.getCriterion(node).getAttributeValue().A;
661  lNode = tree.getLevel(node);// #current level
662 
663  hParent = tree.getCriterion(nParent).getAttributeValue().H;// #current criterion
664  aParent = tree.getCriterion(nParent).getAttributeValue().A;
665  aAncestor = tree.getCriterion(first_ancestor).getAttributeValue().A;
666 
667  lParent = tree.getLevel(nParent);// #current level
668  lAncestor = tree.getLevel(first_ancestor);// #current level
669  int flag;
670 
671 
672  if ((hParent - hNode) <= delta){
673  flag = 1;
674  }
675  else{
676  flag = 0;
677  }
678  float relativeGrowthRate;//BMI RGR
679  double lgARoot = std::log((double)tree.getImWidth()*tree.getImHeight());
680  double factor ;
681 
682  if(method == 2){//RGR
683  factor = ImDtTypes<T>::max()/((double)100*lgARoot);//BMI RGR La valeur max sera max/10 pour un changement d'aire de 1 a taille_im en pixels
684  }
685  else if (method == 3){//MSER_sub
686  factor = ImDtTypes<T>::max()/((double)100.0);//BMI
687  factor = mymax/((double)25.0);//BMI
688  // std::cout<<"mymax="<<int(mymax)<<"; factor3="<<factor<<"\n";
689  }
690  if (flag){
691  current_residue = lNode - lAncestor ;
692  if(method == 1){ // mser stability
693  stability = 1 - ((aAncestor - aNode)*1.0/(aAncestor));
694  stab_residue = round(current_residue * stability);
695  }
696  else if (method == 2){// relative growth rate
697  relativeGrowthRate = factor * (std::log((double)aAncestor ) - std::log((double) aNode ));
698  // std::cout<<"aAncestor:"<<aAncestor.value<<", aNode="<<aNode.value<<", rgr="<<relativeGrowthRate<<"\n";
699  if(current_residue > relativeGrowthRate){
700  stab_residue = round(current_residue - relativeGrowthRate);
701  }
702  else{stab_residue = 0;}
703 
704  // std::cout<<"Hg="<<cParent<<"; Hp="<<cNode<<"; Levels="<<int(lNode)<<", "<<int(lAncestor)<<";;;";
705  // std::cout<<"current="<<int(current_residue)<<"; rel="<<relativeGrowthRate<<"; stb_res="<<int(stab_residue)<<"\n";
706  }// end RGR
707  else if (method == 3){// MSER sustraire
708  stability = factor*( ((aAncestor - aNode)*1.0/(aAncestor)) );
709  // std::cout<<"factor="<<factor<<"current="<<current_residue<<"stab="<<stability<<"\n";
710  if(current_residue > stability){
711  stab_residue = round(current_residue - stability);
712  }
713  else{stab_residue = 0;}
714  }// end MSER_sub
715 
716  // std::cout<<"Ag="<<aAncestor.value<<"; Ap="<<aNode.value<<"; diff_logs="<<(std::log((double)aAncestor.value ) - std::log((double) aNode.value ))<<"; ";
717  }// if flag
718  else{
719  current_residue = (lNode - lParent);
720  if(method == 1){ // mser stability
721  stability = 1 - ((aParent - aNode)*1.0/(aParent));// relative growth rate
722  stab_residue = round(current_residue * stability);
723  }
724  else if (method == 2){// relative growth rate
725  relativeGrowthRate = factor * (std::log((double)aParent ) - std::log((double) aNode ));
726  // std::cout<<"aAncestor:"<<aAncestor.value<<", aNode="<<aNode.value<<", rgr="<<relativeGrowthRate<<"\n";
727  if(current_residue > relativeGrowthRate){
728  stab_residue = round(current_residue - relativeGrowthRate);
729  }
730  else{stab_residue = 0;}
731  // std::cout<<"NFHg="<<cParent<<"; Hp="<<cNode<<"; Levels="<<int(lNode)<<", "<<int(lParent)<<";;;";
732  // std::cout<<"current="<<int(current_residue)<<"; rel="<<int(relativeGrowthRate)<<"; stb_res="<<int(stab_residue)<<"\n";
733  }// end RGR
734  else if (method == 3){
735  stability = factor*( ((aParent - aNode)*1.0/(aParent)));// msersub
736  // std::cout<<"factor="<<factor<<"current="<<current_residue<<"stab="<<stability<<"\n";
737  if(current_residue > stability){
738  stab_residue = round(current_residue - stability);
739  }
740  else{stab_residue = 0;}
741 
742  }
743  else{
744  std::cout<<"ERROR UNKOWN METHOD\n";
745  }
746  }// if not flag
747 
748 
749 
750  transformee_node[node] = transformee_node[nParent];
751  indicatrice_node[node] = indicatrice_node[nParent];
752  int isMaxT = 0;
753  if(hNode < stop){
754  if (stab_residue > transformee_node[node] && stab_residue > threshold && aNode > minArea){
755  isMaxT = 1;
756  transformee_node[node] = stab_residue;
757 
758  if(! (isPrevMaxT && flag)){
759  indicatrice_node[node] = hNode + 1;
760  }
761 
762  }
763  }
764  else{
765  indicatrice_node[node] = 0;
766  }
767  child=tree.getChild(node);
768  while (child!=0){
769  if(flag && (hNode < stop)){
770  ComputeDeltaUOMSER(tree, transformee_node, indicatrice_node, child, node,first_ancestor, stop, delta,method,isMaxT,minArea,threshold,mymax);
771  }
772  else{
773 
774  ComputeDeltaUOMSER(tree, transformee_node, indicatrice_node, child, node, nParent,stop, delta,method,isMaxT,minArea,threshold,mymax);
775  }
776  child = tree.getBrother(child);
777  }
778 }
779 
780 
781 inline void computeFillAspectRatioFactor(UINT wNode,UINT cNode,UINT area,UINT width,UINT height,float &fillRatio,float &AspectRatio)
782 {
783  //compute fillRatio and AspectRatio
784  UINT minHW, maxHW;
785  minHW = MIN(wNode,cNode);
786  maxHW = MAX(wNode,cNode);
787  fillRatio = area*1.0/(wNode * cNode * 1.0);
788  AspectRatio = minHW*1.0/maxHW ;
789 
790  if(AspectRatio <= 0.4)// && fillRatio < 0.4)
791  AspectRatio = 0.0;
792  else
793  AspectRatio = pow(AspectRatio,1.0/3);
794 
795  if(fillRatio <= 0.2 || fillRatio >= 0.9 )
796  fillRatio = 0.0;
797  else{
798  fillRatio = abs(fillRatio - 0.55);
799  fillRatio = pow((1.0- (fillRatio*1.0/0.35)),1.0/3);
800  }
801 
802  if(AspectRatio <= 0.4 && fillRatio >= 0.4){//a verifier
803  AspectRatio = 1;//0.9
804  fillRatio = 1;
805  }//0.9
806 
807  if(area < 20 || area > (0.9* width* height) || cNode < 5 ){//|| cNode > (tree.getImHeight()/3)){
808  AspectRatio = 0.0;
809  fillRatio = 0.0;
810  }
811 }
812 
813 
814 // NEW BMI # ##################################################
815 //(tree, transformee_node, indicatrice_node, child, stopSize, (T)0, 0, hauteur, tree.getLevel(root), tree.getLevel(root));
816  template <class T, class CriterionT, class OffsetT, class LabelT, class tAttType>
817  void ComputeDeltaUOMSERSC(MaxTree2<T,CriterionT,OffsetT,LabelT> &tree, T* transformee_node, tAttType* indicatrice_node,int node, int nParent, int first_ancestor, tAttType stop, UINT delta, int isPrevMaxT){
818 
819  // "node": the current node; "nParent": its direct parent (allows
820  // attribute comparison for Delta versions); "first_ancestor": the
821  // first ancestor with a different attribute (used for residue
822  // computation level(first_ancestor) - level(node) and for area stability computation
823 
824  //self,node = 1, nParent =0, stop=0, delta = 0, isPrevMaxT = 0):
825  int child; // index node
826  T current_residue, stab_residue;
827  UINT hNode, hParent, wNode; // attributes
828  size_t aNode, aParent,aAncestor;
829  T lNode, lParent,lAncestor; /*wParent set but not used*/ // node levels, the same type than input image
830  float stability, fillRatio, AspectRatio, fac;
831 
832 
833  hNode = tree.getCriterion(node).getAttributeValue().H;// #current criterion
834  wNode = tree.getCriterion(node).getAttributeValue().W;// #width
835  aNode = tree.getCriterion(node).getAttributeValue().A;
836  lNode = tree.getLevel(node);// #current level
837 
838 
839  hParent = tree.getCriterion(nParent).getAttributeValue().H;// #current criterion
840 // wParent = tree.getCriterion(nParent).xmax-tree.getCriterion(nParent).xmin+1;// #width
841  aParent = tree.getCriterion(nParent).getAttributeValue().A;
842  aAncestor = tree.getCriterion(first_ancestor).getAttributeValue().A;
843  lParent = tree.getLevel(nParent);// #current level
844  lAncestor = tree.getLevel(first_ancestor);// #current level
845  int flag;
846 
847 
848  if ((hParent - hNode) <= delta){
849  flag = 1;
850  }
851  else{
852  flag = 0;
853  }
854 
855  computeFillAspectRatioFactor(wNode,hNode,aNode,tree.getImWidth(),tree.getImHeight(),fillRatio,AspectRatio);
856 
857  if (flag){//no significant attribute change
858  stability = pow((1.0 - ((aAncestor - aNode)*1.0/aAncestor)),1.0/3);
859  fac = (stability * AspectRatio * fillRatio);
860  current_residue = lNode - lAncestor ;
861  stab_residue = round(current_residue * fac);
862  }
863  else{
864  stability = pow((1.0 - ((aParent - aNode)*1.0/aParent)),1.0/3);
865  fac = (stability * AspectRatio * fillRatio);
866  current_residue = (lNode - lParent);
867  stab_residue = round(current_residue * fac);
868  }
869 
870 
871  transformee_node[node] = transformee_node[nParent];
872  indicatrice_node[node] = indicatrice_node[nParent];
873 
874 
875  int isMaxT = 0;
876  if(hNode < stop){
877 
878  if (stab_residue > transformee_node[node]){
879  isMaxT = 1;
880  transformee_node[node] = stab_residue;
881 
882  if(! (isPrevMaxT && flag)) {
883  indicatrice_node[node] = hNode + 1;
884  }
885 
886  }
887  }
888  else
889  indicatrice_node[node] = 0;
890 
891  child=tree.getChild(node);
892  while (child!=0){
893  if(flag && (hNode < stop)){
894  ComputeDeltaUOMSERSC(tree, transformee_node, indicatrice_node, child, node,first_ancestor, stop, delta,isMaxT);
895  }
896  else{
897  ComputeDeltaUOMSERSC(tree, transformee_node, indicatrice_node, child, node, nParent,stop, delta,isMaxT);
898  }
899  child = tree.getBrother(child);
900  }
901 
902 }
903 
904  template <class T1, class T2>
905  void compute_contrast_MSER(MaxTree2<T1,HWACriterion,size_t,UINT32> &tree, T1* transformee_node, T2* indicatrice_node, UINT32 root, T2 stopSize,UINT delta = 0, UINT method = 2, UINT minArea=0,T1 threshold=0, bool use_textShape =0)
906 {
907 
908  int child;
909 // UINT hauteur;
910 
911  transformee_node[root]=0;
912  indicatrice_node[root]=0;
913 
914  // tree.updateCriteria(root); // already in tree building function
915  // hauteur = tree.getCriterion(root).ymax - tree.getCriterion(root).ymin+1;
916  child = tree.getChild(root);
917 
918  // BEGIN COMPUTE DYNAMIC, BMI
919 
920  UINT32 mynode;
921  T1 mylevel, mymax;
922  mymax = 0;
923  for(mynode = 0; mynode < tree.getLabelMax(); mynode ++){
924  mylevel = tree.getLevel(mynode);
925  if(mylevel > mymax){
926  mymax = mylevel;
927  }
928  }
929  std::cout<<"mymax="<<mymax<<"\n";
930  // END COMPUTE DYNAMIC, BMI
931  while (child!=0)
932  {
933  if(!use_textShape)
934  ComputeDeltaUOMSER(tree, transformee_node, indicatrice_node, child,
935  root/*parent*/,root /*first_ancestor*/,
936  stopSize /*stop*/, delta, method,0 /*isPrevMaxT*/,
937  minArea,threshold,mymax);
938 
939  else{
940  ComputeDeltaUOMSERSC(tree, transformee_node, indicatrice_node, child,
941  root/*parent*/,root /*first_ancestor*/,
942  stopSize /*stop*/, delta, 0 /*isPrevMaxT*/);
943  }
944  child = tree.getBrother(child);
945 
946  }
947 
948 }
949 
950 
951 
952 
953 #endif // SWIG
954 
955  template <class T1, class T2>
956  RES_T ultimateOpenMSER(const Image<T1> &imIn, Image<T1> &imTrans, Image<T2> &imIndic, const StrElt &se = DEFAULT_SE, T2 stopSize=0, UINT delta = 0, UINT method = 2, UINT minArea=0,T1 threshold=0, bool use_textShape =0)
957  {
958  ASSERT_ALLOCATED(&imIn, &imTrans, &imIndic);
959  ASSERT_SAME_SIZE(&imIn, &imTrans, &imIndic);
960 
961  int imSize = imIn.getPixelCount();
962  UINT *img_eti = new UINT[imSize]();
963 
965  int root = tree.build(imIn, img_eti,se);
966 
967  if(stopSize == 0){
968  stopSize= imIn.getHeight()-1;
969  }
970 
971 
972  T1 *transformee_node = new T1[tree.getLabelMax()]();
973  T2 *indicatrice_node = new T2[tree.getLabelMax()]();
974 
975  compute_contrast_MSER(tree, transformee_node, indicatrice_node, root, stopSize,delta,method,minArea,threshold,use_textShape);
976 
977 
978  typename ImDtTypes<T1>::lineType transformeePix = imTrans.getPixels();
979  typename ImDtTypes<T2>::lineType indicatricePix = imIndic.getPixels();
980 
981  for(int i=0;i<imSize;i++)
982  {
983  transformeePix[i]=transformee_node[img_eti[i]];
984  indicatricePix[i]=indicatrice_node[img_eti[i]];
985  }
986 
987  delete[] img_eti;
988  delete[] transformee_node;
989  delete[] indicatrice_node;
990 
991  imTrans.modified();
992  imIndic.modified();
993 
994  return RES_OK;
995  }
996 
997 
998 
999 
1000 
1001 #ifndef SWIG
1002 
1003  template <class T, class CriterionT, class OffsetT, class LabelT, class tAttType>
1004  void processMaxTree(MaxTree2<T,CriterionT,OffsetT,LabelT> &tree, LabelT node, T* lut_out, T previousLevel, tAttType stop)
1005 {
1006  //MORPHEE_ENTER_FUNCTION("s_OpeningTree::computeMaxTree");
1007  LabelT child;
1008  T nodeLevel = tree.getLevel(node);
1009 
1010  T currentLevel = ( tree.getCriterion( node ).getAttributeValue() < stop ) ? previousLevel : nodeLevel;
1011 
1012  lut_out[node] = currentLevel;
1013  // std::cout<<"node:"<<node<<"->"<<int(currentLevel)<<"\n";
1014 
1015  child = tree.getChild(node);
1016  while (child!=0)
1017  {
1018  processMaxTree(tree, child, lut_out, currentLevel, stop);
1019  child = tree.getBrother(child);
1020  }
1021 
1022 }
1023 
1024  template <class T, class CriterionT, class OffsetT, class LabelT, class tAttType>
1025  void compute_AttributeOpening(MaxTree2<T,CriterionT,OffsetT,LabelT> &tree, T* lut_node, LabelT root, tAttType stopSize)
1026 {
1027  LabelT child;
1028 
1029  lut_node[root] = tree.getLevel(root);
1030 
1031  // tree.updateCriteria(root);
1032 
1033  //DEBUG
1034  // LabelT lab;
1035 
1036  // for(lab = 1; lab < tree.getLabelMax(); lab ++){
1037  // std::cout<<"lab="<<lab<<", att="<<tree.getCriterion(lab).getAttributeValue()<<"\n";
1038  // }
1039 
1040  // END DEBUG
1041  child = tree.getChild(root);
1042 
1043  while (child!=0)
1044  {
1045  processMaxTree(tree, child, lut_node, tree.getLevel(root), stopSize);
1046  child = tree.getBrother(child);
1047  }
1048 
1049 
1050 }// compute_AttributeOpening
1051 
1052 
1053  template <class T, class CriterionT, class OffsetT, class LabelT>
1054  RES_T attributeOpen(const Image<T> &imIn, Image<T> &imOut, size_t stopSize, const StrElt &se)
1055  {
1056  ASSERT_ALLOCATED(&imIn, &imOut);
1057  ASSERT_SAME_SIZE(&imIn, &imOut);
1058 
1059  size_t imSize = imIn.getPixelCount();
1060  LabelT *img_eti = new LabelT[imSize]();
1061 
1063  LabelT root = tree.build(imIn, img_eti,se);
1064 
1065 
1066  // BEGIN BMI, DEBUG
1067  // for(size_t i=0;i<imSize;i++) {
1068  // if(i%5 == 0){
1069  // std::cout<<"\n";
1070  // }
1071  // std::cout<<img_eti[i]<<",";
1072  // }
1073  // END BEGIN BMI, DEBUG
1074  T *out_node = new T[tree.getLabelMax()]();
1075 
1076  compute_AttributeOpening(tree, out_node, root, stopSize);
1077 
1078 
1079  typename ImDtTypes<T>::lineType outPix = imOut.getPixels();
1080 
1081  for(size_t i=0;i<imSize;i++)
1082  outPix[i] = out_node[img_eti[i]];
1083 
1084  delete[] img_eti;
1085  delete[] out_node;
1086 
1087  imOut.modified();
1088 
1089  return RES_OK;
1090  }
1091 
1092 #endif // SWIG
1093 
1104  template <class T>
1105  RES_T heightOpen(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1106  {
1107  return attributeOpen<T, HeightCriterion,size_t,UINT32>(imIn, imOut, stopSize,se);
1108  }// END heightOpen
1109 
1110  template <class T>
1111  RES_T heightClose(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1112  {
1113  ASSERT_ALLOCATED(&imIn, &imOut);
1114  ASSERT_SAME_SIZE(&imIn, &imOut);
1115 
1116  ImageFreezer freeze(imOut);
1117 
1118  Image<T> tmpIm(imIn);
1119  inv(imIn, tmpIm);
1120  RES_T res = attributeOpen<T, HeightCriterion,size_t,UINT32>(tmpIm, imOut, stopSize, se);
1121  inv(imOut, imOut);
1122 
1123  return res;
1124  }// END heightClose
1125 
1136  template <class T>
1137  RES_T widthOpen(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1138  {
1139  return attributeOpen<T, WidthCriterion,size_t,UINT32>(imIn, imOut, stopSize,se);
1140  }// END widthOpen
1141 
1142  template <class T>
1143  RES_T widthClose(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1144  {
1145  ASSERT_ALLOCATED(&imIn, &imOut);
1146  ASSERT_SAME_SIZE(&imIn, &imOut);
1147 
1148  ImageFreezer freeze(imOut);
1149 
1150  Image<T> tmpIm(imIn);
1151  inv(imIn, tmpIm);
1152  RES_T res = attributeOpen<T, WidthCriterion,size_t,UINT32>(tmpIm, imOut, stopSize, se);
1153  inv(imOut, imOut);
1154 
1155  return res;
1156  }// END widthClose
1157 
1168  template <class T>
1169  RES_T areaOpen(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1170  {
1171  return attributeOpen<T, AreaCriterion,size_t,UINT32>(imIn, imOut, stopSize,se);
1172  }
1173 
1184  template <class T>
1185  RES_T areaClose(const Image<T> &imIn, size_t stopSize, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
1186  {
1187  ASSERT_ALLOCATED(&imIn, &imOut);
1188  ASSERT_SAME_SIZE(&imIn, &imOut);
1189 
1190  ImageFreezer freeze(imOut);
1191 
1192  Image<T> tmpIm(imIn);
1193  inv(imIn, tmpIm);
1194  RES_T res = attributeOpen<T, AreaCriterion,size_t,UINT32>(tmpIm, imOut, stopSize, se);
1195  inv(imOut, imOut);
1196 
1197  return res;
1198  }
1199 
1202 } // namespace smil
1203 
1204 
1205 #endif // _D_MAX_TREE_HPP
1206 
RES_T areaOpen(const Image< T > &imIn, size_t stopSize, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Area opening.
Definition: DMorphoMaxTree.hpp:1169
RES_T widthOpen(const Image< T > &imIn, size_t stopSize, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Width opening.
Definition: DMorphoMaxTree.hpp:1137
RES_T threshold(const Image< T > &imIn, T minVal, T maxVal, T_out trueVal, T_out falseVal, Image< T_out > &imOut)
Image threshold.
Definition: DImageHistogram.hpp:170
Definition: DColorConvert.h:38
void updateCriteria(const int node)
Update criteria of a given max-tree node.
Definition: DMorphoMaxTree.hpp:408
Definition: DBaseImage.h:235
vector< IntPoint > points
List of neighbor points.
Definition: DStructuringElement.h:125
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
RES_T inv(const Image< T > &imIn, Image< T > &imOut)
Invert an image.
Definition: DImageArith.hpp:375
RES_T areaClose(const Image< T > &imIn, size_t stopSize, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Area closing.
Definition: DMorphoMaxTree.hpp:1185
RES_T ultimateOpen(const Image< T1 > &imIn, Image< T1 > &imTrans, Image< T2 > &imIndic, const StrElt &se=DEFAULT_SE, T2 stopSize=0, UINT delta=0)
Ultimate Opening using the max-trees.
Definition: DMorphoMaxTree.hpp:594
RES_T log(const Image< T > &imIn, Image< T > &imOut, int base=0)
Logarithm.
Definition: DImageArith.hpp:833
size_t getHeight() const
Get image height.
Definition: DBaseImage.h:91
RES_T heightOpen(const Image< T > &imIn, size_t stopSize, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Height opening.
Definition: DMorphoMaxTree.hpp:1105
size_t area(const Image< T > &imIn)
Area of an image.
Definition: DMeasures.hpp:70
Definition: DMorphoMaxTree.hpp:67
Base structuring element.
Definition: DStructuringElement.h:51
virtual void modified()
Trigger modified event (allows to force display update)
Definition: DImage.hxx:223
size_t getWidth() const
Get image width.
Definition: DBaseImage.h:87
size_t getPixelCount() const
Get the number of pixels.
Definition: DBaseImage.h:150
Definition: DTypes.hpp:78
Main Image class.
Definition: DQVtkViewer.hpp:44
void getCoordsFromOffset(size_t off, size_t &x, size_t &y, size_t &z) const
Get x,y(,z) coordinates for a given offset.
Definition: DBaseImage.h:182
Definition: DCommon.h:104