A.I.M.S algorithms


clusterArg_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_ROI_CLUSTERARG_D_H
36 #define AIMS_ROI_CLUSTERARG_D_H
37 
38 #include <aims/roi/clusterArg.h>
39 #include <aims/utility/threshold.h>
43 #include <aims/resampling/motion.h>
44 #include <aims/graph/graphmanip.h>
45 #include <aims/bucket/bucket.h>
46 #include <aims/mesh/surface.h>
48 #include <aims/mesh/mesher.h>
49 #include <graph/graph/graph.h>
50 #include <cartobase/smart/rcptr.h>
52 
53 namespace aims
54 {
55 
56  template<typename T>
57  void ClusterArgMaker::make( Graph & gr, const AimsData<T> & data )
58  {
59  // graph and global attributes
60 
61  std::string graphsyntax = gr.getSyntax();
62  std::string vertexsyntax = "cluster";
63  GraphElementCode bec, mec;
64 
65  enum GraphType
66  {
67  Cluster,
68  Fold
69  };
70 
71  GraphType gtype = Cluster;
72 
73  // bucket element code
74  bec.objectType = "Bucket";
75  bec.dataType = "VOID";
77 
78  // mesh element code
79  mec.objectType = "Mesh";
80  mec.dataType = "VOID";
82 
83  if( graphsyntax == "CorticalFoldArg" )
84  {
85  gtype = Fold;
86  vertexsyntax = "fold";
87 
88  bec.attribute = "aims_ss";
89  bec.global_index_attribute = "ss_label";
90  bec.local_file_attribute = "ss_filename";
91  bec.global_filename = "ss.bck";
92  bec.id = "ss";
93 
94  mec.attribute = "aims_Tmtktri";
95  mec.global_index_attribute = "Tmtktri_label";
96  mec.local_file_attribute = "Tmtktri_filename";
97  mec.global_filename = "tmtktri.mesh";
98  mec.id = "Tmtktri";
99  }
100  else
101  {
102  bec.attribute = "aims_bucket";
103  bec.global_index_attribute = "bucket_label";
104  bec.local_file_attribute = "bucket_filename";
105  bec.global_filename = "cluster.bck";
106  bec.id = "cluster";
107 
108  mec.attribute = "aims_Tmtktri";
109  mec.global_index_attribute = "Tmtktri_label";
110  mec.local_file_attribute = "Tmtktri_filename";
111  mec.global_filename = "Tmtktri.mesh";
112  mec.id = "Tmtktri";
113  }
114  bec.syntax = vertexsyntax;
115  mec.syntax = vertexsyntax;
116 
117  Vertex *v;
118 
119  gr.setProperty( graphsyntax + "_VERSION", std::string( "1.0" ) );
120  std::vector<float> vs;
121  vs.push_back( data.sizeX() );
122  vs.push_back( data.sizeY() );
123  vs.push_back( data.sizeZ() );
124  vs.push_back( data.sizeT() );
125  gr.setProperty( "voxel_size", vs );
126  std::vector<int> bb;
127  bb.push_back( 0 );
128  bb.push_back( 0 );
129  bb.push_back( 0 );
130  gr.setProperty( "boundingbox_min", bb );
131  bb[0] = data.dimX() - 1;
132  bb[1] = data.dimY() - 1;
133  bb[2] = data.dimZ() - 1;
134  gr.setProperty( "boundingbox_max", bb );
135  std::string outdir = _fileout;
136  std::string::size_type pos = outdir.rfind( '.' );
137  if( pos != std::string::npos )
138  outdir.erase( pos, outdir.length() - pos );
139  pos = outdir.rfind( '/' );
140  if( pos != std::string::npos )
141  outdir.erase( 0, pos + 1 );
142  outdir += ".data";
143  gr.setProperty( "filename_base", outdir );
144  std::vector<int> colors;
145  colors.push_back( 0 );
146  colors.push_back( 255 );
147  colors.push_back( 0 );
148  gr.setProperty( bec.local_file_attribute, colors );
149  gr.setProperty( bec.global_index_attribute, colors );
150 
152  objtable( new GraphElementTable );
153  std::map<std::string,GraphElementCode>
154  & cltab = (*objtable)[ vertexsyntax ];
155  cltab[ bec.id ] = bec;
156  if( _domesh )
157  {
158  cltab[ mec.id ] = mec;
159  gr.setProperty( mec.local_file_attribute, colors );
160  gr.setProperty( mec.global_index_attribute, colors );
161  }
162  gr.setProperty( "aims_objects_table", objtable );
163 
164  // thresholding
165 
166  bool thresh = false;
168 
169  if( _upth != FLT_MAX )
170  {
171  thresh = true;
172  if( _lowth != -FLT_MAX )
173  thtype = AIMS_BETWEEN;
174  else
175  {
176  thtype = AIMS_LOWER_OR_EQUAL_TO;
177  _lowth = _upth;
178  _upth = 0;
179  }
180  }
181  else if( _lowth != -FLT_MAX )
182  {
183  thresh = true;
184  thtype = AIMS_GREATER_OR_EQUAL_TO;
185  _upth = 0;
186  }
187  else if( _binarize )
188  {
189  thresh = true;
190  if( _lowth == -FLT_MAX )
191  {
192  _lowth = 0;
193  thtype = AIMS_GREATER_THAN;
194  }
195  // _upth = 0;
196  }
197 
198  AimsData<T> datat;
199 
200  if( thresh )
201  {
202  AimsThreshold<T, T> th( thtype, (T) _lowth,
203  (T) _upth );
204  std::cout << "threshold: " << _lowth << ", bin: " << _binarize
205  << std::endl;
206  if( _binarize )
207  datat = th.bin( data );
208  else
209  datat = th( data );
210  }
211  else
212  datat = data;
213 
214  // connected components
215 
216  std::cout << "extracting connected components...\n";
217 
219  datac( datat.dimX(), datat.dimY(), datat.dimZ(), datat.dimT(), 1 );
221  conv2.convert( datat, datac );
222  datac.fillBorder( -1 );
223  /*
224  const AttributedHeader
225  *ah = dynamic_cast<const AttributedHeader *>( data.header() );
226  if( ah )
227  datac.setHeader( new AttributedHeader( *ah ) );
228  datac.setSizeXYZT( data.sizeX(), data.sizeY(), data.sizeZ(),
229  data.sizeT() );*/
230 
231  // find background level
232  short bg = 0;
233  bool bgfound = false;
234  int x, y, z;
235 
236  for( z=0; z<datat.dimZ() && !bgfound; ++z )
237  for( y=0; y<datat.dimY() && !bgfound; ++y )
238  for( x=0; x<datat.dimX(); ++x )
239  if( datat( x, y, z ) == 0 )
240  {
241  bg = datac( x, y, z );
242  bgfound = true;
243  break;
244  }
245 
246  AimsConnectedComponent( datac, _connectivity, bg, false, _minsize, 0 );
247  std::cout << "done\n";
248 
249  // transform
250 
251  Motion motion;
252  if( !_matrix.empty() )
253  {
254  std::cout << "reading transformation...\n";
255  std::ifstream ms( _matrix.c_str() );
256  if( !ms )
257  {
258  std::cerr << "Can't open transformation file " << _matrix
259  << std::endl;
260  throw carto::file_error( _matrix );
261  }
262  std::vector<float> transl( 3 ), rot( 9 );
263  ms >> transl[0] >> transl[1] >> transl[2];
264  ms >> rot[0] >> rot[1] >> rot[2];
265  ms >> rot[3] >> rot[4] >> rot[5];
266  ms >> rot[6] >> rot[7] >> rot[8];
267  if( !ms )
268  {
269  std::cerr << "transformation file " << _matrix << " has wrong format\n";
270  throw carto::file_error( _matrix );
271  }
272  std::cout << "done\n";
273  // translation is in other form: t' = r^-1 * t
274  AimsData<float> & rm = motion.rotation();
275  rm(0,0) = rot[0];
276  rm(0,1) = rot[1];
277  rm(0,2) = rot[2];
278  rm(1,0) = rot[3];
279  rm(1,1) = rot[4];
280  rm(1,2) = rot[5];
281  rm(2,0) = rot[6];
282  rm(2,1) = rot[7];
283  rm(2,2) = rot[8];
284  Motion mot2 = motion.inverse();
285  Point3df & mtr = motion.translation();
286  mtr[0] = transl[0];
287  mtr[1] = transl[1];
288  mtr[2] = transl[2];
289 
290  Point3df t2 = mot2.transform( transl[0], transl[1], transl[2] );
291  transl[0] = t2[0];
292  transl[1] = t2[1];
293  transl[2] = t2[2];
294  gr.setProperty( "Talairach_translation", transl );
295  gr.setProperty( "Talairach_rotation", rot );
296  // Talairach_scale ?
297  }
298 
299  // resampling
300  // well, we shouldn't resample anymore!
301 
302  // buckets
303 
304  std::cout << "converting to buckets...\n";
306  BucketMap<Void> bck;
307  conv.convert( datac, bck );
308  std::cout << "done, " << bck.size() << " buckets\n";
309 
310  // meshes
311 
312  std::map<size_t, AimsSurfaceTriangle> mesh;
313 
314  if( _domesh )
315  {
316  std::cout << "meshing components...\n";
317  std::map<size_t, std::list<AimsSurfaceTriangle> > meshl;
318  Mesher mesher;
319  mesher.setDecimation( 100.0, _deciMaxClearance, _deciMaxError, 120.0 );
320  mesher.setMinFacetNumber( _minFacetNumber );
321  mesher.setSmoothing( Mesher::LAPLACIAN, 5, 0.4 );
322  mesher.setSmoothingLaplacian( 180. );
323  mesher.doit( datac, meshl );
324  std::cout << "done\n";
325  //if( onenode )
326  //{
327  std::cout << "regrouping meshes...\n";
328  std::map<size_t, std::list<AimsSurfaceTriangle> >::iterator
329  im, em = meshl.end();
330  for( im=meshl.begin(); im!=em; ++im )
331  {
332  //std::cout << "* " << im->first << ", " << im->second.size();
333  AimsSurfaceTriangle & mesh2 = mesh[im->first];
334  //im->second.erase( im->second.begin() );
335  SurfaceManip::meshMerge( mesh2, im->second );
336  //std::cout << " -> " << mesh2.vertex().size() << " vertices\n";
337  }
338  std::cout << "done\n";
339  //}
340  }
341 
342  // vertices
343 
344  BucketMap<Void> *b;
345  BucketMap<Void>::iterator ib, eb = bck.end();
346  char num[50];
347  float voxvol = bck.sizeX() * bck.sizeY() * bck.sizeZ();
348  int index = 0;
350  Point3d pt;
352  std::vector<float> pf(3);
353  Point3df pm;
354 
355  std::cout << "making cluster vertices...\n";
356  for( ib=bck.begin(); ib!=eb; ++ib )
357  {
358  std::cout << "cluster " << ib->first << "..." << std::endl;
359  v = gr.addVertex( vertexsyntax );
360  b = new BucketMap<Void>;
361  b->setSizeXYZT( bck.sizeX(), bck.sizeY(), bck.sizeZ(), bck.sizeT() );
362  (*b)[0] = ib->second;
363  v->setProperty( bec.attribute, carto::rc_ptr<BucketMap<Void> >( b ) );
364  v->setProperty( "index", (int) ib->first );
365  pt = ib->second.begin()->first;
366  sprintf( num, "%d", (int) datat( pt[0], pt[1], pt[2] ) ); // label
367  v->setProperty( "name", std::string( num ) );
368  v->setProperty( "point_number", (int) ib->second.size() );
369  v->setProperty( "size", voxvol * ib->second.size() );
370  if( gtype == Fold )
371  {
372  v->setProperty( "skeleton_label", index++ );
373  v->setProperty( "ss_point_number", (int) ib->second.size() );
374  v->setProperty( "bottom_point_number", (int) 0 );
375  }
376  if( _domesh )
377  {
378  s = new AimsSurfaceTriangle( mesh[ ib->first ] );
379  v->setProperty( mec.attribute,
381  }
382  // cluster center
383  pf[0] = pf[1] = pf[2] = 0;
384  BucketMap<Void>::Bucket & bk = b->begin()->second;
385  for( ibb=bk.begin(), ebb=bk.end(); ibb!=ebb; ++ibb )
386  {
387  pf[0] += ibb->first[0];
388  pf[1] += ibb->first[1];
389  pf[2] += ibb->first[2];
390  }
391  pf[0] = vs[0] * pf[0] / bk.size();
392  pf[1] = vs[1] * pf[1] / bk.size();
393  pf[2] = vs[2] * pf[2] / bk.size();
394  v->setProperty( "gravity_center", pf );
395  // in talairach space
396  pm = motion.transform( pf[0], pf[1], pf[2] );
397  pf[0] = pm[0];
398  pf[1] = pm[1];
399  pf[2] = pm[2];
400  if( gtype == Fold )
401  v->setProperty( "refgravity_center", pf );
402  else
403  v->setProperty( "Talairach(mm)", pf );
404  sprintf( num, "%f %f %f", pf[0], pf[1], pf[2] );
405  v->setProperty( "label", std::string( num ) );
406  }
407  std::cout << "clusters created\n";
408 
409  // Relations ?
410  }
411 
412 }
413 
414 #endif
Vertex * addVertex(const std::string &s="")
std::string local_file_attribute
int dimX() const
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
float sizeY() const
void setSmoothing(SmoothingType smoothType, int smoothIt, float smoothRate)
std::map< std::string, std::map< std::string, GraphElementCode > > GraphElementTable
int dimT() const
void setMinFacetNumber(uint m)
Definition: mesher.h:169
std::string global_index_attribute
void make(Graph &gr, const AimsData< T > &data)
Definition: clusterArg_d.h:57
float sizeX() const
float sizeT() const
float sizeZ() const
void fillBorder(const T &val)
Point3dd transform(double x, double y, double z) const
void setSizeXYZT(float sizex, float sizey, float sizez, float sizet)
AimsData< U > bin(const AimsData< T > &sqv)
Definition: mesher.h:51
float sizeT() const
AimsData< float > & rotation()
void doit(const AimsData< short > &thing, const std::string &name, const std::string &mode="binar")
threshold_t
void AimsConnectedComponent(AimsData< T > &data, aims::Connectivity::Type connectivity, std::map< T, size_t > &valids, const T &backgrnd=0, bool bin=true, size_t minSize=0, size_t numMax=0, bool verbose=true)
float sizeZ() const
void setSmoothingLaplacian(float smoothFeatureAngle)
float sizeX() const
float sizeY() const
static void meshMerge(AimsTimeSurface< D, T > &dst, const AimsTimeSurface< D, T > &add)
int dimZ() const
int dimY() const
AffineTransformation3d inverse() const
std::string global_filename
void setDecimation(float deciReductionRate, float deciMaxClearance, float deciMaxError, float deciFeatureAngle)