aimsalgo  5.1.2
Neuroimaging image processing
meshmorphomat_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 
35 #ifndef AIMS_DISTANCEMAP_MESHMORPHOMAT_D_H
36 #define AIMS_DISTANCEMAP_MESHMORPHOMAT_D_H
37 
38 #include <aims/distancemap/meshmorphomat.h>
39 #include <aims/distancemap/stlsort.h>
40 #include <aims/distancemap/meshdistance.h>
41 #include <aims/connectivity/meshcc_d.h>
42 #include <stack>
43 #include <float.h>
44 
45 namespace aims
46 {
47  namespace meshdistance
48  {
49 
50 template<class T>
51 Texture<T> LineariseSulci( const AimsSurface<3,Void> & mesh,
52  const Texture<T> & inittex,
53  const T & label,
54  const T & Back, const T & For,
55  unsigned & min, unsigned & max, unsigned nbDil)
56 {
57 
58  Texture<float> dist;
59  Texture<T> tex = inittex;
60  bool allowUnreached = false;
61  const std::vector<Point3df> & vert = mesh.vertex();
62  unsigned n = vert.size(),imax,i,j,jmax,nbpt=0;
63  TimeTexture<float> distance(11,n);
64  float distMax = 0;
65  std::vector<unsigned> nodeMax;
66  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
67  ASSERT (n == tex.nItem() );
68 
69 
70  //init object map
71  for (i=0;i<n;++i)
72  {
73  if (inittex.item(i) == label)
74  {
75  tex.item(i) = Back;
76  ++nbpt;
77  }
78  else
79  tex.item(i) = For;
80  }
81 
82  if (nbpt == 0 || nbpt == 1)
83  {
84  std::cout << "The label " << label << " has only one node ! \n";
85  return(inittex);
86  }
87  //ncc = AimsMeshLabelNbConnectedComponent(mesh,tex,Back);
88  //ASSERT ( ncc == 1 );
89 
90  // neighbours map
91 
92  std::map<unsigned, std::set<unsigned> > neighbours;
93  unsigned v1, v2, v3;
94 
95  for( i=0; i<poly.size(); ++i )
96  {
97  v1 = poly[i][0];
98  v2 = poly[i][1];
99  v3 = poly[i][2];
100  if(tex.item(v1)!=For
101  && tex.item(v2)!=For)
102  {
103  neighbours[v1].insert( v2 );
104  neighbours[v2].insert( v1 );
105  }
106  if(tex.item(v1)!=For
107  && tex.item(v3)!=For)
108  {
109  neighbours[v1].insert( v3 );
110  neighbours[v3].insert( v1 );
111  }
112  if(tex.item(v2)!=For
113  && tex.item(v3)!=For)
114  {
115  neighbours[v2].insert( v3 );
116  neighbours[v3].insert( v2 );
117  }
118  }
119 
120  //init distance map
121  i=0;
122  while (tex.item(i) == For)
123  ++ i;
124  tex.item(i) = label;
125 
126  dist = MeshDistance<T>(mesh, tex, allowUnreached);
127  distance[0] = dist;
128 
129  for (i=0;i<n;++i)
130  if (distance[0].item(i) >= distMax)
131  {
132  distMax = distance[0].item(i);
133  imax = i;
134  }
135 
136  nodeMax.push_back(imax);
137 
138  //Distance map
139  for (j=1;j<5;++j)
140  {
141  distMax = 0;
142  for (i=0;i<n;++i)
143  {
144  if (inittex.item(i) > 0)
145  tex.item(i) = Back;
146  else
147  tex.item(i) = For;
148  }
149  tex.item(imax) = 1;
150  dist = MeshDistance<T>(mesh, tex, allowUnreached);
151  distance[j] = dist;
152  for (i=0;i<n;++i)
153  if (dist.item(i) > distMax )
154  {
155  distMax = dist.item(i);
156  imax = i;
157  }
158  nodeMax.push_back(imax);
159  }
160 
161  Texture<T> sulci(n);
162  std::set<unsigned>::iterator in,en;
163  float dmin;
164  min = imax = nodeMax[4];
165  max = jmax = nodeMax[3];
166  std::stack<unsigned> line;
167  i=imax;
168  while (i != jmax)
169  {
170  line.push(i);
171  dmin = FLT_MAX;
172  for (in = neighbours[i].begin(), en = neighbours[i].end(); in !=en ; ++in)
173  if (dist.item(*in) < dmin)
174  {
175  i = *in;
176  dmin = dist.item(*in);
177  }
178  sulci.item(i) = label;
179  }
180 
181  unsigned inc = 0, l = line.size();
182 
183  std::cout << "Size of the line: " << l << " nbDil: " << nbDil << std::endl;
184  if (nbDil != 0)
185  while (!line.empty())
186  {
187  ++inc;
188  i = line.top();
189  line.pop();
190  if (inc <= nbDil || inc >= l - nbDil)
191  sulci.item(i) = Back;
192  }
193 
194  for (i=0;i<n;++i)
195  {
196  if (inittex.item(i) != label)
197  sulci.item(i) = inittex.item(i);
198  }
199 
200  return (sulci);
201 }
202 
203 
204 
205 
206 
207 
208 
209 
210 template<class T>
211 Texture<T> SplitSulci( const AimsSurface<3,Void> & mesh,
212  const Texture<T> & inittex,
213  const T & Back, const T & For ,
214  const T label,
215  const unsigned min, const unsigned max)
216 {
217 
218 
219  const std::vector<Point3df> & vert = mesh.vertex();
220  unsigned ncc, n = vert.size(),i;
221  Texture<T> dilat(n),line(n),split(n);
222  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
223 
224 
225  ASSERT (n = inittex.nItem() );
226 
227  // Forbidden neighbours map
228 
229  std::set<unsigned> neighbours;
230  unsigned v1, v2, v3;
231  std::set<unsigned>::iterator is,es;
232 
233  i = min;
234  v1 = poly[i][0];
235  v2 = poly[i][1];
236  v3 = poly[i][2];
237  neighbours.insert(v1);
238  neighbours.insert(v2);
239  neighbours.insert(v3);
240  i = max;
241  v1 = poly[i][0];
242  v2 = poly[i][1];
243  v3 = poly[i][2];
244  neighbours.insert(v1);
245  neighbours.insert(v2);
246  neighbours.insert(v3);
247  neighbours.insert(min);
248  neighbours.insert(max);
249 
250  for (i=0;i<n;++i)
251  {
252  if (inittex.item(i) == label)
253  dilat.item(i) = label;
254  else
255  dilat.item(i) = Back;
256  }
257 
258 
259 
260  for (is=neighbours.begin(), es=neighbours.end();is != es; ++is)
261  //if (dilat.item(*is) == Back)
262  dilat.item(*is) = For;
263 
264 
265  dilat = MeshDilation<T>(mesh,dilat,Back,For,1,true);
266 
267  for (i=0;i<n;++i)
268  split.item(i) = dilat.item(i);
269 
270 
271  for (i=0;i<n;++i)
272  if (inittex.item(i) == label || dilat.item(i) == For )
273  split.item(i) = Back;
274 
275 
276  ncc = AimsMeshLabelNbConnectedComponent(mesh,split,label);
277  std::cout << "ncc: " << ncc << std::endl;
278 
279 
280  return(split);
281 
282 }
283 
284 
285 
286 //closing and thinning of the sulci
287 template<class T>
288 Texture<T> CloseSulci( const AimsSurface<3,Void> & mesh,
289  const Texture<T> & sulctex,
290  float nbDilMax,bool connexity,
291  const T & Back, const T & For,
292  const std::vector<std::list<unsigned> > & neigho)
293 {
294 
295  const std::vector<Point3df> & vert = mesh.vertex();
296  unsigned inc=0,nbLabel,i,/*nbDil,*/n = vert.size();
297  Texture<T> sulctemp(n),ccsulcus(n);
298  Texture<T> tex(n);
299  short ncc;
300  T label;
301  //std::set<short> labsulci;
302  std::set<T> cc_sulci_labels;//labels of the connex components
303  typename std::set<T>::iterator ilab,elab;
304  bool hole=true;
305 
306  // init
307  for (i=0; i<n; ++i)
308  {
309  cc_sulci_labels.insert(sulctex.item(i));
310  tex.item(i) = sulctex.item(i);
311  }
312 
313  cc_sulci_labels.erase(Back);
314  cc_sulci_labels.erase(For);
315  nbLabel = cc_sulci_labels.size();
316  std::cout << "\nNb of labels: " << nbLabel << std::endl;
317 
318 
319  for ( ilab = cc_sulci_labels.begin(), elab=cc_sulci_labels.end(); ilab!=elab; ++ilab)
320  {
321  ++inc;
322  //nbDil=0;
323  label = *ilab;
324  std::cout << inc << "/" << nbLabel << std::endl;
325  for (i=0; i<n; ++i)
326  {
327  if (sulctex.item(i) == label)
328  ccsulcus.item(i) = label;
329  else
330  ccsulcus.item(i) = Back;
331  }
332  ncc = AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label);
333  std::cout << "Nb initial of cc: " << ncc << std::endl;
334 
335  if (ncc == 1)
336  {
337  hole = HasHole(ccsulcus,mesh,label);
338 
339  if (hole)
340  std::cout << "With hole(s) \n";
341  else
342  std::cout << "Without hole \n";
343  }
344  float nbd=1;
345  unsigned NB_PT_MAX = 10,nbpt;
346 
347  //homotopic closing
348  while ( (ncc != 1 || hole) && nbd <= nbDilMax)
349  {
350  sulctemp = MeshDilation<T>( mesh, ccsulcus,Back,For,nbd,connexity);
351  ncc = AimsMeshLabelNbConnectedComponent(mesh,sulctemp,label);
352  std::cout << "After Dilation (" << nbd << ") : cc: " << ncc << std::endl ;
353  if (ncc == 1)
354  {
355  hole = HasHole(sulctemp,mesh,label);
356  if (hole)
357  std::cout << "With hole(s) \n";
358  else
359  std::cout << "Without hole \n";
360  }
361  nbd += 1;
362  }
363 
364  if (ncc != 1)
365  std::cerr <<"Not able to close the projected connnected component, but I continue.\n";
366  if (nbd == 1)
367  sulctemp = ccsulcus;
368  ccsulcus = MeshDilation<T>( mesh, ccsulcus,Back,For,1,connexity);
369  std::cout << "Skeletonization" << std::endl;
370  ccsulcus = MeshSkeletization( mesh,sulctemp,label,Back,neigho);
371  nbpt = NbOfPoint(sulctemp,label);
372  if ( nbpt > NB_PT_MAX)
373  ccsulcus = MeshSkeletization( mesh,ccsulcus,label,Back,neigho);//2 steps for the ebarbulage
374  else
375  std::cout << "Not enough point for the ebarbulage ( " << nbpt << " < " << NB_PT_MAX << ")" << std::endl;
376  ncc = AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label);
377  std::cout << "ncc: " << ncc << std::endl ;
378  hole = HasHole(ccsulcus,mesh,label);
379  if (hole)
380  std::cout << "With hole(s) \n";
381  else
382  std::cout << "Without hole \n";
383 
384  //old thinning
385  /*
386  if (ncc != 1 && nbDil > 0)
387  {
388  std::cout << "Impossible to close the sulci of label " << label << std::endl;
389  ccsulcus = AimsMeshFilterConnectedComponent( mesh, ccsulcus,label );
390  }
391  if (AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label) == 1)
392  {
393  std::cout << "Boundaries thinning \n";
394  ccsulcus = LineariseSulci<short>( mesh, ccsulcus ,label,Back,For,min,max,nbDil);
395  //std::cout << "Split sulci \n";
396  //ccsulcus = SplitSulci(mesh,ccsulcus,Back,For,label,min,max);
397 
398  }
399  */
400 
401  //check if the node is free
402  for (i=0;i<n;++i)
403  if (ccsulcus.item(i) == label)
404  {
405  if (sulctex.item(i)==label || sulctex.item(i) == Back )
406  tex.item(i)=label;
407  else
408  std::cout << "The point " << i << " cannot change its label \n";
409  }
410 
411 
412  for (i=0;i<n;++i)
413  if (tex.item(i) == label)
414  {
415  if (ccsulcus.item(i) == label)
416  tex.item(i) = label;
417  else
418  tex.item(i) = Back ;
419  }
420  }
421 
422  return(tex);
423 }
424 
425 
426 
427 //id but no thinning for label in the set label_forbidden
428 template<class T>
429 Texture<T> CloseSulci( const AimsSurface<3,Void> & mesh,
430  const Texture<T> & sulctex,
431  float nbDilMax,bool connexity,
432  const T & Back, const T & For,
433  const std::vector<std::list<unsigned> > & neigho,
434  const std::set<T> & label_forbidden)
435 {
436 
437  const std::vector<Point3df> & vert = mesh.vertex();
438  unsigned inc=0,nbLabel,i,/*nbDil,*/n = vert.size();
439  Texture<T> sulctemp(n),ccsulcus(n);
440  Texture<T> tex(n);
441  short ncc;
442  T label;
443  //std::set<short> labsulci;
444  std::set<T> cc_sulci_labels;//labels of the connex components
445  typename std::set<T>::iterator ilab,elab;
446  bool hole=true;
447 
448  // init
449  for (i=0; i<n; ++i)
450  {
451  cc_sulci_labels.insert(sulctex.item(i));
452  tex.item(i) = sulctex.item(i);
453  }
454 
455  cc_sulci_labels.erase(Back);
456  cc_sulci_labels.erase(For);
457  nbLabel = cc_sulci_labels.size();
458  std::cout << "\nNb of labels: " << nbLabel << std::endl;
459 
460  typename std::set<T>::iterator lfe = label_forbidden.end();
461 
462  for ( ilab = cc_sulci_labels.begin(), elab=cc_sulci_labels.end(); ilab!=elab; ++ilab)
463  {
464  ++inc;
465  //nbDil=0;
466  label = *ilab;
467  std::cout << inc << "/" << nbLabel << std::endl;
468  for (i=0; i<n; ++i)
469  {
470  if (sulctex.item(i) == label)
471  ccsulcus.item(i) = label;
472  else
473  ccsulcus.item(i) = Back;
474  }
475  ncc = AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label);
476  std::cout << "Nb initial of cc: " << ncc << std::endl;
477 
478  if (ncc == 1)
479  {
480  hole = HasHole(ccsulcus,mesh,label);
481 
482  if (hole)
483  std::cout << "With hole(s) \n";
484  else
485  std::cout << "Without hole \n";
486  }
487  float nbd=1;
488  unsigned NB_PT_MAX = 10,nbpt;
489 
490  //homotopic closing
491  while ( (ncc != 1 || hole) && nbd <= nbDilMax)
492  {
493  sulctemp = MeshDilation<T>( mesh, ccsulcus,Back,For,nbd,connexity);
494  ncc = AimsMeshLabelNbConnectedComponent(mesh,sulctemp,label);
495  std::cout << "After Dilation (" << nbd << ") : cc: " << ncc << std::endl ;
496  if (ncc == 1)
497  {
498  hole = HasHole(sulctemp,mesh,label);
499  if (hole)
500  std::cout << "With hole(s) \n";
501  else
502  std::cout << "Without hole \n";
503  }
504  nbd += 1;
505  }
506 
507  if (ncc != 1)
508  std::cerr <<"Not able to close the projected connnected component, but I continue.\n";
509  if (nbd == 1)
510  sulctemp = ccsulcus;
511  ccsulcus = MeshDilation<T>( mesh, ccsulcus,Back,For,1,connexity);
512  std::cout << "Skeletonization" << std::endl;
513 
514  if (label_forbidden.find(*ilab) == lfe)
515  {
516  ccsulcus = MeshSkeletization( mesh,sulctemp,label,Back,neigho);
517  nbpt = NbOfPoint(sulctemp,label);
518  if ( nbpt > NB_PT_MAX)
519  ccsulcus = MeshSkeletization( mesh,ccsulcus,label,Back,neigho);//2 steps for the ebarbulage
520  else
521  std::cout << "Not enough point for the ebarbulage ( " << nbpt << " < " << NB_PT_MAX << ")" << std::endl;
522  ncc = AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label);
523  std::cout << "ncc: " << ncc << std::endl ;
524  hole = HasHole(ccsulcus,mesh,label);
525  if (hole)
526  std::cout << "With hole(s) \n";
527  else
528  std::cout << "Without hole \n";
529  }
530  else
531  std::cout << "No thinning for the INSULA\n";
532 
533  //old thinning
534  /*
535  if (ncc != 1 && nbDil > 0)
536  {
537  std::cout << "Impossible to close the sulci of label " << label << std::endl;
538  ccsulcus = AimsMeshFilterConnectedComponent( mesh, ccsulcus,label );
539  }
540  if (AimsMeshLabelNbConnectedComponent(mesh,ccsulcus,label) == 1)
541  {
542  std::cout << "Boundaries thinning \n";
543  ccsulcus = LineariseSulci<short>( mesh, ccsulcus ,label,Back,For,min,max,nbDil);
544  //std::cout << "Split sulci \n";
545  //ccsulcus = SplitSulci(mesh,ccsulcus,Back,For,label,min,max);
546 
547  }
548  */
549 
550  //check if the node is free
551  for (i=0;i<n;++i)
552  if (ccsulcus.item(i) == label)
553  {
554  if (sulctex.item(i)==label || sulctex.item(i) == Back )
555  tex.item(i)=label;
556  else
557  std::cout << "The point " << i << " cannot change its label \n";
558  }
559 
560 
561  for (i=0;i<n;++i)
562  if (tex.item(i) == label)
563  {
564  if (ccsulcus.item(i) == label)
565  tex.item(i) = label;
566  else
567  tex.item(i) = Back ;
568  }
569  }
570 
571 
572  return(tex);
573 }
574 
575 
576 namespace
577 {
578 
579  template<class T> inline
580  bool SimplePoint( const Texture<T> & tex,
581  const std::list<unsigned> & neigho)
582  {
583 
584  std::list<unsigned>::const_iterator il=neigho.begin(), el=neigho.end();
585  short l=tex.item(*il);
586  int nb=0;//nb chgmt of value
587 
588  while(il != el)
589  {
590  if (tex.item(*il) != l)
591  {
592  ++nb;
593  l = tex.item(*il);
594  }
595  ++il;
596  }
597  return(nb == 1 || nb == 2);
598  }
599 
600 }
601 
602 template<class T> inline
603 bool HasHole( const Texture<T> &tex,
604  const AimsSurface<3,Void> & mesh, const T label )
605 {
606 
607  const std::vector<Point3df> & vert = mesh.vertex();
608  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
609  unsigned i, n = vert.size();
610  unsigned S=0,A=0,F=0;
611 
612  ASSERT( tex.nItem() == n );
613 
614  // neighbours map
615 
616  std::set<std::set<unsigned> > neighbours;
617  std::set<unsigned> temp;
618  unsigned v1, v2, v3;
619 
620  for( i=0; i<poly.size(); ++i )
621  {
622  v1 = poly[i][0];
623  v2 = poly[i][1];
624  v3 = poly[i][2];
625 
626  if(tex.item(v1)==label && tex.item(v2)==label)
627  {
628  temp.clear();
629  temp.insert(v1);
630  temp.insert(v2);
631  neighbours.insert(temp);
632  }
633 
634  if(tex.item(v1)==label && tex.item(v3)==label)
635  {
636  temp.clear();
637  temp.insert(v1);
638  temp.insert(v3);
639  neighbours.insert(temp);
640  }
641 
642  if(tex.item(v3)==label && tex.item(v2)==label)
643  {
644  temp.clear();
645  temp.insert(v2);
646  temp.insert(v3);
647  neighbours.insert(temp);
648  }
649 
650  if(tex.item(v3)==label && tex.item(v2)==label && tex.item(v1)==label)
651  ++F;
652  }
653 
654  //Nb de sommets S
655  for (i=0;i<n;++i)
656  if (tex.item(i) == label)
657  ++S;
658 
659  //Nb d'aretes A
660  A = neighbours.size();
661 
662 
663  return( (S - A + F != 1) );
664 
665 }
666 
667 
668 namespace
669 {
670 
671  template<class T> inline
672  bool ImmortalPoint( const Texture<T> & tex,
673  const std::list<unsigned> & neigho )
674  {
675 
676  std::list<unsigned>::const_iterator il=neigho.begin(), el=neigho.end();
677  short l=tex.item(*il);
678  int nb=0;//nb chgmt of value
679 
680  while(il != el)
681  {
682  if (tex.item(*il) != l)
683  {
684  ++nb;
685  l = tex.item(*il);
686  }
687  ++il;
688  }
689  il = neigho.begin();
690  if (tex.item(*il) != l)
691  ++nb;
692 
693  return(nb > 2 || nb ==1 );
694  }
695 
696 }
697 
698 
699 
700 template<class T> inline
701 unsigned NbOfPoint( const Texture<T> tex, const T label )
702 {
703  unsigned nb = 0;
704  for (unsigned i=0; i < tex.nItem(); ++i)
705  if (tex.item(i) == label)
706  ++nb;
707 
708  return(nb);
709 
710 }
711 
712 
713 template<class T>
714 Texture<T> MeshSkeletization( const AimsSurface<3,Void> & mesh,
715  const Texture<T> & inittex,
716  T label, T Back,
717  const std::vector<std::list<unsigned> > & neigho )
718 {
719 
720  unsigned i,n = mesh.vertex().size();
721  Texture<short> skelet(n),temp(n);
722  Texture<short> immortal(n); //0 : still alive , 1: immortal
723  Texture<T> tex(n);
724  std::list<unsigned>::const_iterator il,el;
725 
726  for (i=0; i<n; ++i)
727  if (inittex.item(i) == label)
728  skelet.item(i) = 1;
729 
730  //init
731 
732 
733  //Def initial immortal point map
734  for (i=0; i<n; ++i)
735  if ( skelet.item(i) == 1 && ImmortalPoint(skelet,neigho[i]) )
736  immortal.item(i) = 1;
737 
738  for (i=0; i<n; ++i)
739  if (skelet.item(i) == 1)
740  temp.item(i) = 0;
741  else
742  temp.item(i) = 1;
743 
744  Texture<float> dist;
745  std::multimap<float,unsigned> dfront;
746  std::multimap<float,unsigned>::iterator id,ed;
747  dist = MeshDistance(mesh,temp,true);
748 
749 
750  // node of the front are ordered according to their distance from the background
751  for (i=0; i<n; ++i)
752  if (skelet.item(i) == 1 && immortal.item(i) == 0)
753  for (il=neigho[i].begin(), el=neigho[i].end();il!=el;++il)
754  if (skelet.item(*il) == 0 )
755  dfront.insert( std::pair<float,unsigned>(dist.item(i),i) );
756 
757 
758  while (dfront.size()>0)
759  {
760 
761  for (id=dfront.begin(), ed=dfront.end();id != ed; ++id)
762  if ( immortal.item(id->second) == 0 )
763  {
764  skelet.item(id->second) = 0;
765  //upgrade immortal map
766  for (il=neigho[id->second].begin(), el=neigho[id->second].end();il!=el;++il)
767  if ( skelet.item(*il) == 1 && ImmortalPoint(skelet,neigho[*il]) )
768  immortal.item(*il) = 1;
769  }
770 
771  //upgrade front
772 
773  dfront.clear();
774  for (i=0; i<n; ++i)
775  if (skelet.item(i) == 1 && immortal.item(i) == 0)
776  for (il=neigho[i].begin(), el=neigho[i].end();il!=el;++il)
777  if (skelet.item(*il) == 0 )
778  dfront.insert( std::pair<float,unsigned>(dist.item(i),i) );
779  }
780 
781  for (i=0; i<n; ++i)
782  if ( skelet.item(i) == 1)
783  tex.item(i) =label;
784  else
785  tex.item(i) = Back;
786 
787  return(tex);
788 
789 }
790 
791 template <class T>
792 Texture<T> MeshDilation( const AimsSurface<3,Void> & mesh,
793  const Texture<T> & inittex,
794  const T & Back, const T & For,
795  const float dist,bool connexity )
796 {
797 
798  Texture<T> tex;
799  tex = MeshVoronoiT<T>( mesh, inittex, Back, For, dist, connexity,true );
800  return(tex);
801 
802 }
803 
804 template <class T>
805 Texture<T> MeshErosion( const AimsSurface<3,Void> & mesh,
806  const Texture<T> & inittex,
807  const T & Back, const T & For,
808  const float dist,bool connexity )
809 {
810 
811  Texture<T> tex;
812  tex = MeshVoronoiT<T>( mesh, inittex, Back,For,dist, connexity,false );
813  return(tex);
814 }
815 
816  }
817 }
818 
819 #endif
820 
ASSERT
#define ASSERT(EX)
AimsSurface::vertex
const std::vector< Point3df > & vertex() const
AimsSurface::polygon
const std::vector< AimsVector< uint, D > > & polygon() const
Texture::nItem
size_t nItem() const
Texture::item
const T & item(int n) const
TimeTexture::item
const T & item(int n) const
AimsMeshLabelNbConnectedComponent
unsigned AimsMeshLabelNbConnectedComponent(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T lab)
Definition: meshcc_d.h:160
meshcc_d.h
aims::meshdistance::SplitSulci
Texture< T > SplitSulci(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T &Back, const T &For, const T label, const unsigned min, const unsigned max)
Definition: meshmorphomat_d.h:211
aims::meshdistance::MeshErosion
Texture< T > MeshErosion(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T &Back, const T &For, const float dist, bool connexity)
Definition: meshmorphomat_d.h:805
aims::meshdistance::MeshDilation
Texture< T > MeshDilation(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T &Back, const T &For, const float dist, bool connexity)
Definition: meshmorphomat_d.h:792
aims::meshdistance::min
float min(float x, float y)
Definition: thickness.h:106
aims::meshdistance::CloseSulci
Texture< T > CloseSulci(const AimsSurface< 3, Void > &mesh, const Texture< T > &sulctex, float nbDilMax, bool connexity, const T &Back, const T &For, const std::vector< std::list< unsigned > > &neigho, const std::set< T > &label_forbidden)
Close the sulcus map (dilation and skeletization).
Definition: meshmorphomat_d.h:429
aims::meshdistance::max
float max(float x, float y)
Definition: thickness.h:98
aims::meshdistance::HasHole
bool HasHole(const Texture< T > &tex, const AimsSurface< 3, Void > &mesh, const T label)
Definition: meshmorphomat_d.h:603
aims::meshdistance::LineariseSulci
Texture< T > LineariseSulci(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const T &label, const T &Back, const T &For, unsigned &min, unsigned &max, unsigned nbDil)
Definition: meshmorphomat_d.h:51
aims::meshdistance::MeshSkeletization
Texture< T > MeshSkeletization(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, T label, T Back, const std::vector< std::list< unsigned > > &neigho)
Skeletization of the components of the innitex with label label.
Definition: meshmorphomat_d.h:714
aims::meshdistance::NbOfPoint
unsigned NbOfPoint(const Texture< T > tex, const T label)
Definition: meshmorphomat_d.h:701
aims::meshdistance::MeshDistance
Texture< float > MeshDistance(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, bool allowUnreached, float max_dist=FLT_MAX)
Computes a distance texture over a mesh.
Definition: meshdistance_d.h:42
split
std::vector< std::string > split(const std::string &text, const std::string &sep)