aimsalgo  5.1.2
Neuroimaging image processing
meshparcellation_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_MESHPARCELLATION_D_H
36 #define AIMS_DISTANCEMAP_MESHPARCELLATION_D_H
37 
38 #include <aims/distancemap/meshparcellation.h>
39 #include <aims/distancemap/meshdistance.h>
40 #include <aims/distancemap/meshmorphomat.h>
41 #include <aims/distancemap/stlsort.h>
42 #include <aims/mesh/surfaceOperation.h>
43 #include <aims/mesh/surface.h>
44 #include <cartobase/containers/nditerator.h>
45 #include <stack>
46 #include <set>
47 #include <float.h>
48 
49 namespace aims
50 {
51  namespace meshdistance
52  {
53 
54 // Select the gyrus seeds
55 template <class T>
56 Texture<std::set<T> >
57 gyrusSeedDefinition( const AimsSurface<3,Void> & mesh,
58  const Texture<T> & inittex,
59  const std::set<T> & setBack, const std::set<T> & setFor,
60  const std::set<std::set<T> > & labelAllowed)
61 {
62 
63  const std::vector<Point3df> & vert = mesh.vertex();
64  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
65  unsigned i, n = vert.size();
66  Texture<std::set<T> > tex(n);
67  std::map <unsigned, std::set<T> > labelNeigh;//vertex -> labels of neighbours
68  typename std::map <unsigned, std::set<T> >::iterator il,el;
69 
70  ASSERT( inittex.nItem() == n );
71 
72  // neighbours map
73  std::map<unsigned, std::set<unsigned> > neighbours;
74  unsigned v1, v2, v3;
75 
76  for( i=0; i<poly.size(); ++i )
77  {
78  v1 = poly[i][0];
79  v2 = poly[i][1];
80  v3 = poly[i][2];
81  labelNeigh[v1].insert( inittex.item(v2) );
82  labelNeigh[v2].insert( inittex.item(v1) );
83  labelNeigh[v1].insert( inittex.item(v3) );
84  labelNeigh[v3].insert( inittex.item(v1) );
85  labelNeigh[v3].insert( inittex.item(v2) );
86  labelNeigh[v2].insert( inittex.item(v3) );
87  }
88 
89  std::set<std::set<short> >::const_iterator ia,dl;
90  std::set<short> setLabel,setTemp;
91 
92  // bool ok;
93  // for ( il = labelNeigh.begin(), el = labelNeigh.end(); il!= el; ++il)
94  // {
95  // ok = false;
96  // if ( (il->second).size() >= 2 )
97  // {
98  // setLabel.clear();
99  // setLabel = il->second;
100  // for (ia = labelAllowed.begin(),dl = labelAllowed.end(); ia != dl; ++ia)
101  // {
102  // setTemp.clear();
103  // std::insert_iterator<std::set<short> > ii(setTemp, setTemp.begin() );
104  // set_intersection(setLabel.begin(),setLabel.end(),ia->begin(),ia->end(),ii );
105  // if (setTemp.size() >= 2)
106  // {
107  // tex.item(il->first) = *ia;
108  // ok = true;
109  // break;
110  // }
111  // }
112  // if (!ok)
113  // tex.item(il->first) = setBack;
114  // }
115  // else
116  // {
117  // if ( il->second != setFor )
118  // tex.item(il->first) = setBack;
119  // else
120  // tex.item(il->first) = setFor;
121  // }
122  // }
123 
124  Texture<std::set<T> > tex2(n);
125  unsigned interSize,sizeNeigh;
126  typename std::map<unsigned,std::set<T> >::reverse_iterator im,em;
127  std::multimap<unsigned,std::set<T> > matchedLabel;
128 
129  //bool ok;
130  for ( il = labelNeigh.begin(), el = labelNeigh.end(); il!= el; ++il)
131  {
132  //ok = false;
133  sizeNeigh = (il->second).size();
134  //2 neighbours
135 
136  /*
137  if ( sizeNeigh == 2)
138  {
139  setLabel = il->second;
140  for (ia = labelAllowed.begin(),dl = labelAllowed.end(); ia != dl; ++ia)
141  if( (*ia).size() == 2 )
142  {
143  setTemp.clear();
144  std::insert_iterator<std::set<T> > ii(setTemp, setTemp.begin() );
145  set_intersection(setLabel.begin(),setLabel.end(),ia->begin(),ia->end(),ii );
146  if (setTemp.size() == 2 )
147  {
148  tex.item(il->first) = *ia;
149  ok = true;
150  break;
151  }
152  }
153  if (!ok)
154  tex.item(il->first) = setBack;
155  }
156  else
157  */
158  //3 neighbours
159  //if ( sizeNeigh == 3)
160  if ( sizeNeigh > 1 )
161  {
162  setLabel = il->second;
163  matchedLabel.clear();
164 
165  typename std::set<T>::iterator i1,i2;
166 
167  for (ia = labelAllowed.begin(),dl = labelAllowed.end(); ia != dl; ++ia)
168  {
169  setTemp.clear();
170  std::insert_iterator<std::set<T> > ii(setTemp, setTemp.begin() );
171  set_intersection(setLabel.begin(),setLabel.end(),ia->begin(),ia->end(),ii );
172  interSize = setTemp.size();
173  if ( interSize >= 2 )
174  matchedLabel.insert(std::pair<unsigned, std::set<T> >(interSize,*ia) ) ;
175  }
176 
177  if (matchedLabel.size() == 1)
178  {
179  im = matchedLabel.rbegin(); //best matched label
180  tex.item(il->first) = im->second;
181  }
182  else
183  {
184  tex.item(il->first) = setBack;
185  //if (matchedLabel.size() > 1)
186  //std::cerr << "Cannot define the gyrus label" << std::endl;
187  }
188  }
189  else
190  {
191  if ( il->second != setFor )
192  tex.item(il->first) = setBack;
193  else
194  tex.item(il->first) = setFor;
195  }
196  }
197 
198 
199  return (tex);
200 }
201 
202 // Dilate the gyrus s
203 // in the SKIZ
204 template <class T>
205 Texture<T>
206 gyrusSeedDilationInSKIZ( const AimsSurface<3,Void> & mesh,
207  const Texture<T> & seed,
208  const Texture<T> & skiz,
209  const T & Back,
210  const T & For)
211 {
212 
213  const std::vector<Point3df> & vert = mesh.vertex();
214  unsigned i, n = vert.size();
215  Texture<T> tex(n);
216 
217  ASSERT( seed.nItem() == n || skiz.nItem() == n);
218 
219  for(i=0; i<n; ++i)
220  if( skiz.item( i ) == Back )
221  tex.item( i ) = seed.item(i);
222  else
223  tex.item( i ) = For;
224 
225  tex = MeshDilation( mesh,tex, Back, For,FLT_MAX,false);
226 
227  for(i=0; i<n; ++i)
228  if( skiz.item( i ) != Back )
229  tex.item ( i ) = seed.item( i );
230 
231 
232  return (tex);
233 }
234 
235 
236 // Extract the point of the border of a parcellisation
237 template <class T>
238 Texture<std::set<T> >
239 MeshBorderVoronoi( const AimsSurface<3,Void> & mesh,
240  const Texture<T> & inittex,
241  const std::set<T> & setBack, const std::set<T> & setFor )
242 {
243 
244  const std::vector<Point3df> & vert = mesh.vertex();
245  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
246  unsigned i, n = vert.size();
247  Texture<std::set<T> > tex(n);
248  std::map <unsigned, std::set<T> > labelNeigh;//vertex -> labels of neighbours
249  typename std::map <unsigned, std::set<T> >::iterator il,el;
250 
251  ASSERT( inittex.nItem() == n );
252 
253  // neighbours map
254  std::map<unsigned, std::set<unsigned> > neighbours;
255  unsigned v1, v2, v3;
256 
257  for( i=0; i<poly.size(); ++i )
258  {
259  v1 = poly[i][0];
260  v2 = poly[i][1];
261  v3 = poly[i][2];
262  labelNeigh[v1].insert( inittex.item(v2) );
263  labelNeigh[v2].insert( inittex.item(v1) );
264  labelNeigh[v1].insert( inittex.item(v3) );
265  labelNeigh[v3].insert( inittex.item(v1) );
266  labelNeigh[v3].insert( inittex.item(v2) );
267  labelNeigh[v2].insert( inittex.item(v3) );
268  }
269 
270  unsigned sizeNeigh;
271 
272  for ( il = labelNeigh.begin(), el = labelNeigh.end(); il!= el; ++il)
273  {
274  sizeNeigh = (il->second).size();
275  if ( sizeNeigh >= 2)
276  tex.item(il->first) = setBack;
277  else
278  tex.item(il->first) = setFor;
279  }
280 
281 
282  return (tex);
283 }
284 
285 
286 
287 
288 
289 
290 
291 // Convert the texture of set to a texture of short labels
292 template <class T>
293 Texture<short>
294 border2Texture( const Texture<std::set<T> > &tex,
295  const AimsSurface<3,Void> & mesh,
296  const std::set<T> & setBack, const std::set<T> & setFor)
297 {
298  const std::vector<Point3df> & vert = mesh.vertex();
299  unsigned i, n = vert.size();;
300  Texture<short> outtex(n);
301  std::set<std::set<T>,SetCompare<T> > labSet;
302  typename std::set<std::set<T>,SetCompare<T> >::iterator il,el;
303  std::map<std::set<T>,short,SetCompare<T> > newLab;
304  unsigned inc = 0;
305 
306  for (i=0;i<n;++i)
307  labSet.insert( tex.item(i) );
308 
309  //typename std::set<T>::iterator is,es;
310  /*
311  cout << "Label of the frontier \n";
312  for (il=labSet.begin(), el = labSet.end(); il!=el;++il)
313  {
314  for (is=il->begin(),es=il->end();is!=es;++is)
315  cout << *is << " ";
316  cout << std::endl;
317  }
318  */
319  il=labSet.begin();
320  if (*il == setFor)
321  {
322  newLab[*il] = -1;
323  ++il;
324  }
325  if (*il == setBack)
326  {
327  newLab[*il] = 0 ;
328  ++il;
329  }
330  el=labSet.end();
331 
332  while(il!=el)
333  {
334  ++inc;
335  newLab[*il]=inc;
336  ++il;
337  }
338 
339  for (i=0;i<n;++i)
340  if ( tex.item(i).size() > 1 )
341  outtex.item(i) = newLab[tex.item(i)];
342 
343  return(outtex);
344 }
345 
346 
347 
348 template<class T, class U>
349 Texture<short>
350 Voronoi2toTexture( const Texture<std::set<T> > & vor,
351  const AimsSurface<3,Void> & mesh,
352  const std::set<T> & setBack, const std::set<T> & setFor )
353 {
354  std::set<std::set<U>,SetCompare<U> > lab_front_sulci;
355  typename std::set<std::set<U>,SetCompare<U> >::iterator il,el;
356  std::set<U> setTemp;
357  typename std::set<T>::iterator is,es;
358  const std::vector<Point3df> & vert = mesh.vertex();
359  unsigned i, n = vert.size();
360  Texture<std::set<U> > sTex(n);
361  Texture<short> Tex(n);
362  std::map<std::set<U>, short> lab_set_short;
363  //map<T,U> cc = cc_sulci;
364  ASSERT( n == vor.nItem() );
365 
366  for (i=0;i<n;++i)
367  {
368  setTemp.clear();
369  for (is = vor.item(i).begin(), es = vor.item(i).end();is != es; ++is)
370  setTemp.insert( *is );
371  //setTemp.insert( cc[*is] );
372  sTex.item(i) = setTemp;
373  lab_front_sulci.insert(setTemp);
374  }
375 
376  unsigned inc = 1;
377 
378  lab_set_short[setFor] = MESHDISTANCE_FORBIDDEN;
379  lab_set_short[setBack] = 0;
380 
381  for ( il=lab_front_sulci.begin(), el= lab_front_sulci.end(); il!=el;++il )
382  if (*il != setFor && *il != setBack)
383  {
384  lab_set_short[*il] = inc;
385  ++inc;
386  }
387  /*
388  cout << "Correspondance label frontier -> sulci\n";
389  for ( il=lab_front_sulci.begin(), el= lab_front_sulci.end(); il!=el;++il )
390  {
391  for (is = (*il).begin(), es = (*il).end(); is !=es; ++is)
392  cout << *is << " " ;
393  cout << " -> " << lab_set_short[*il] << std::endl;
394  }
395  */
396  for (i=0;i<n;++i)
397  Tex.item(i) = lab_set_short[sTex.item(i)];
398 
399  return(Tex);
400 
401 }
402 
403 
404 
405 template<class T>
406 std::map<T,float>
407 SurfaceParcel( const Texture<T> & tex, const AimsSurface<3,Void> & mesh )
408 {
409 
410  std::map<T,float> stat;
411  const std::vector<Point3df> & vertex = mesh.vertex();
412  unsigned i,n = vertex.size();
413  const std::vector< AimsVector<uint,3> > & poly = mesh.polygon();
414  ASSERT( n == tex.nItem() );
415  std::map<T,std::set<unsigned> > labels;
416  typename std::map<T,std::set<unsigned> >::iterator il,el=labels.end();
417  T lab;
418  for (i=0;i<n;++i)
419  labels[tex.item(i)].insert(i);
420 
421  // blob surface
422 
423  for (il=labels.begin();il != el; ++il)
424  {
425  std::set<unsigned> vset = il->second;
426  lab = il->first;
427  std::set<unsigned>::iterator evset = vset.end();
428  float surface=0;
429  Point3df zero( 0.0, 0.0, 0.0 ),AB,AC,H,CH;
430 
431  for( unsigned ii=0, n=poly.size(); ii!=n; ++ii )
432  if( vset.find( poly[ii][0] ) != evset
433  && vset.find( poly[ii][1] ) != evset
434  && vset.find( poly[ii][2] ) != evset )
435  {
436  const AimsVector<uint,3> & pol = poly[ii];
437  AB = vertex[ pol[ 1 ] ] - vertex[ pol[ 0 ] ];
438  AC = vertex[ pol[ 2 ] ] - vertex[ pol[ 0 ] ];
439  H = vertex[ pol[ 0 ] ] + ( AB.dot( AC ) / norm2( AB ) ) * AB;
440  CH = H - vertex[ pol[ 2 ] ];
441  if ( AB != zero )
442  surface += norm( CH ) * norm( AB ) / 2;
443  }
444  stat[lab] = surface;
445  //cout << lab << " -> " << surface << " mm2" << std::endl;
446  }
447 
448  return(stat);
449 }
450 
451 
452 
453 
454 
455 template<class T>
456 std::map<T,float>
457 VolumeParcel( const carto::rc_ptr<carto::Volume<T> > & vol )
458 {
459 
460  std::map<T,float> stat;
461  std::set<T> labels;
462  typename std::set<T>::iterator il,el;
463  int x,y,z;
464  //T lab;
465  std::vector<float> vs = vol->getVoxelSize();
466  float voxelVol = vs[0] * vs[1] * vs[2];
467 
468  const T *p, *pp;
469  carto::const_line_NDIterator<T> it( &vol->at( 0 ), vol->getSize(),
470  vol->getStrides(), true );
471  for( ; !it.ended(); ++it )
472  {
473  p = &*it;
474  for( pp=p + it.line_length(); p!=pp; it.inc_line_ptr( p ) )
475  labels.insert( *p );
476  }
477 
478  for( il = labels.begin(), el = labels.end(); il != el; ++il )
479  stat[*il] = 0;
480 
481  for( it.reset(); !it.ended(); ++it )
482  {
483  p = &*it;
484  for( pp=p + it.line_length(); p!=pp; it.inc_line_ptr( p ) )
485  stat[*p] += voxelVol;
486  }
487 
488  return(stat);
489 }
490 
491  }
492 }
493 
494 #endif
ASSERT
#define ASSERT(EX)
AimsSurface::vertex
const std::vector< Point3df > & vertex() const
AimsSurface::polygon
const std::vector< AimsVector< uint, D > > & polygon() const
AimsVector::dot
T dot(const AimsVector< T, D > &other) const
Texture::nItem
size_t nItem() const
Texture::item
const T & item(int n) const
iterator
T * iterator
const_iterator
const T * const_iterator
carto::NDIterator_base::ended
bool ended() const
carto::NDIterator_base::reset
void reset()
carto::const_line_NDIterator
carto::const_line_NDIterator::line_length
long line_length() const
carto::const_line_NDIterator::inc_line_ptr
void inc_line_ptr(const T *&p) const
carto::rc_ptr
aims::meshdistance::MeshBorderVoronoi
Texture< std::set< T > > MeshBorderVoronoi(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const std::set< T > &setBack, const std::set< T > &setFor)
Extract the boundarie of a voronoi diagram i.e.
Definition: meshparcellation_d.h:239
aims::meshdistance::border2Texture
Texture< short > border2Texture(const Texture< std::set< T > > &tex, const AimsSurface< 3, Void > &mesh, const std::set< T > &setBack, const std::set< T > &setFor)
Convert a multidimensional boundary map to a texture of short.
Definition: meshparcellation_d.h:294
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::VolumeParcel
std::map< T, float > VolumeParcel(const carto::rc_ptr< carto::Volume< T > > &vol)
Give the surface of labelled parcels defined from volume.
Definition: meshparcellation_d.h:457
aims::meshdistance::Voronoi2toTexture
Texture< short > Voronoi2toTexture(const Texture< std::set< T > > &vor, const AimsSurface< 3, Void > &mesh, const std::set< T > &setBack, const std::set< T > &setFor)
Definition: meshparcellation_d.h:350
aims::meshdistance::gyrusSeedDefinition
Texture< std::set< T > > gyrusSeedDefinition(const AimsSurface< 3, Void > &mesh, const Texture< T > &inittex, const std::set< T > &setBack, const std::set< T > &setFor, const std::set< std::set< T > > &labelAllowed)
Define the gyrus seed from the SKIZ.
Definition: meshparcellation_d.h:57
aims::meshdistance::gyrusSeedDilationInSKIZ
Texture< T > gyrusSeedDilationInSKIZ(const AimsSurface< 3, Void > &mesh, const Texture< T > &seed, const Texture< T > &skiz, const T &Back, const T &For)
Definition: meshparcellation_d.h:206
aims::meshdistance::MESHDISTANCE_FORBIDDEN
const short MESHDISTANCE_FORBIDDEN
global variable...
Definition: meshvoronoi.h:84
aims::meshdistance::SurfaceParcel
std::map< T, float > SurfaceParcel(const Texture< T > &tex, const AimsSurface< 3, Void > &mesh)
Give the surface of labelled parcels defined from textures.
Definition: meshparcellation_d.h:407
norm
AIMSDATA_API float norm(const Tensor &thing)
norm2
float norm2(const AimsVector< T, D > &v1)