aimsalgo 6.0.0
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
49namespace aims
50{
51 namespace meshdistance
52 {
53
54// Select the gyrus seeds
55template <class T>
56Texture<std::set<T> >
57gyrusSeedDefinition( 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
204template <class T>
205Texture<T>
206gyrusSeedDilationInSKIZ( 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
237template <class T>
238Texture<std::set<T> >
239MeshBorderVoronoi( 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
292template <class T>
293Texture<short>
294border2Texture( 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
348template<class T, class U>
349Texture<short>
350Voronoi2toTexture( 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
405template<class T>
406std::map<T,float>
407SurfaceParcel( 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
455template<class T>
456std::map<T,float>
457VolumeParcel( 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::polygon
const std::vector< AimsVector< uint, D > > & polygon() const
AimsSurface::vertex
const std::vector< Point3df > & vertex() 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
carto::NDIterator_base::ended
bool ended() const
carto::NDIterator_base::reset
void reset()
carto::const_line_NDIterator
carto::const_line_NDIterator::line_length
std::ptrdiff_t 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)
Point3df
AimsVector< float, 3 > Point3df
norm2
float norm2(const AimsVector< T, D > &v1)