A.I.M.S algorithms


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