A.I.M.S algorithms


projection.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_PROJECTION_H
36 #define AIMS_DISTANCEMAP_PROJECTION_H
37 
38 #include <cstdlib>
39 #include <aims/mesh/texture.h>
40 #include <aims/mesh/surface.h>
41 #include <aims/data/data.h>
44 #include <aims/math/gausslu.h>
45 #include <aims/math/eigen.h>
46 #include <set>
47 #include <list>
48 #include <map>
49 #include <float.h>
50 
51 
52 class Graph;
53 namespace aims
54 {
55  namespace meshdistance
56  {
57 
58 
66  const Texture<float> & curvtex,
67  const AimsData <short> & bottom_vol,
68  const AimsData <short> & surface_vol,
69  float alpha, float dmin,
70  int MINCC,
71  const std::map <short,std::string> & trans,
72  const std::set<std::string> & labels,
73  float maxdil,float maxdil_mesh,float alpha_reg,
74  bool connexity,
75  const std::vector<std::list<unsigned> > & neigho);
79  const AimsData <short> & bottom_vol,
80  const AimsData <short> & surface_vol,
81  const Point3df & CA,
82  float demin,float dpmin,
83  int MINCC,
84  const std::map <short,std::string> & trans,
85  const std::set<std::string> & labels,
86  float maxdil,float maxdil_mesh,float alpha_reg,
87  bool connexity,
88  const std::vector<std::list<unsigned> > & neigho);
89 
90 
93  const AimsData<short> &sulcvol,
94  const AimsSurface<3,Void> & mesh,
95  std::set<short> &cc_sulci_labels,
96  std::map<short,short > &tri_sulci,
97  std::map<unsigned, Point3dfSet > &cc_sulci_coord,
98  std::map<Point3df, Point3df,
99  Point3dfCompare > &initend,
100  const Texture<float> &curvtex,
101  float alpha,float dmin);
102 
105  ( const AimsData<short> &ccvol, const AimsData<short> &sulcvol,
106  const AimsSurface<3,Void> & mesh,
107  std::set<short> &cc_sulci_labels,
108  std::map<short,short > &tri_sulci,
109  const std::map<Point3df, Point3dfSet,
110  Point3dfCompare > &neigh,
111  std::map<unsigned, Point3dfSet > &cc_sulci_coord,
112  std::map<Point3df, Point3df, Point3dfCompare > &initend,
113  float demin, float dpmin,
114  short label_insula_left,short label_insula_right,
115  const Point3df & CA );
116 
117  void
119  const AimsSurface<3,Void> & mesh,
120  const std::map<Point3df, Point3df, Point3dfCompare >
121  &initend,
122  const Point3dfSet &setpointi,
123  short label,
124  float dmin,float alpha_reg );
125 
127  const Texture<short> &tex,
128  const AimsSurface<3,Void> & mesh,
129  short val_domain,
130  short back);
131 
134  const AimsSurface<3,Void> & mesh,
135  short back);
136 
139  void
140  NeighbourInCC(std::map<Point3df, Point3dfSet, Point3dfCompare > &neigh,
141  const AimsData <short> & bvol,const AimsData <short> & svol,
142  Connectivity::Type connectivity,
143  short label, const unsigned size_neigh,
144  unsigned max_points = 50);
145 
146 
147 
148  //inline functions
149 
150  //first : barycenter of the input pt
151  //second : normal vector
152 
153  inline std::pair<Point3df,Point3df> NormalFromPoints( const Point3dfSet & init_point )
154  {
155  std::pair<Point3df,Point3df> output;
156  unsigned inc = 0, i=0,j=0,n = init_point.size();
157  Point3df mean = Point3df(0,0,0),temp;
158  Point3dfSet norm_point;
159  AimsData<float> cov(3,3),w(3) ;
160  Point3dfSet::const_iterator is,es;
161 
162  for (is = init_point.begin(), es = init_point.end(); is != es; ++is)
163  mean += *is;
164 
165  mean = mean / (float)n;
166 
167  for (is = init_point.begin(), es = init_point.end(); is != es; ++is,
168  ++inc)
169  norm_point.insert(*is - mean);
170 
171  cov = 0;
172  for (is = norm_point.begin(), es = norm_point.end(); is != es; ++is,
173  ++inc)
174  {
175  temp = *is;
176  for (i=0; i<3; ++i)
177  for (j=0; j<3; ++j)
178  cov(i,j) += temp[i] * temp[j];
179  }
180  AimsEigen< float > eigen;
181  AimsData< float > value = eigen.doit( cov );
182 
183  eigen.sort(cov,value);
184 
185  temp[0] = cov(0,2);
186  temp[1] = cov(1,2);
187  temp[2] = cov(2,2);
188 
189  temp.normalize();
190  output.second = temp;
191  output.first = mean;
192 
193  return(output);
194  }
195 
196  //Point3df : normal vector to the plane passing through pt
197  inline Point3df NormalFromPoints( const Point3df & pt,
198  const Point3dfSet & init_point )
199  {
200  unsigned inc = 0, i=0,j=0;
201  Point3df temp;
202  Point3dfSet norm_point;
203  AimsData<float> cov(3,3),w(3) ;
204  Point3dfSet::const_iterator is,es;
205 
206 
207 
208  for (is = init_point.begin(), es = init_point.end(); is != es; ++is, ++inc)
209  norm_point.insert(*is - pt);
210 
211  cov = 0;
212  for (is = norm_point.begin(), es = norm_point.end(); is != es; ++is, ++inc)
213  {
214  temp = *is;
215  for (i=0; i<3; ++i)
216  for (j=0; j<3; ++j)
217  cov(i,j) += temp[i] * temp[j];
218  }
219  AimsEigen< float > eigen;
220  AimsData< float > value = eigen.doit( cov );
221 
222  eigen.sort(cov,value);
223 
224  temp[0] = cov(0,2);
225  temp[1] = cov(1,2);
226  temp[2] = cov(2,2);
227 
228  temp.normalize();
229 
230 
231  return(temp);
232  }
233 
234  inline Point3df cross(const Point3df &A, const Point3df &B)
235  {
236  Point3df C;
237  C[0] = A[1] * B[2] - A[2] * B[1];
238  C[1] = A[2] * B[0] - A[0] * B[2];
239  C[2] = A[0] * B[1] - A[1] * B[0];
240  return (C);
241 
242  }
243 
244 
245  //Give the point of intersection
246  //between a line
247  //(defined by a point 'pt'
248  //and a direction 'direction_line' )
249  // and a plane
250  //(defined by 3 points s1, s2, s3)
251  inline Point3df IntersectionPointOfPlanAndLine(const Point3df &s1,const Point3df &s2,const Point3df &s3,
252  const Point3df & direction_line, const Point3df & pt)
253  {
254  float a0, b0, c0;
255  float xb,yb,zb;
256  float a,b,c;
257  float xm,ym,zm;
258  Point3df pt_intersection, normal_tri;
259 
260 
261  ASSERT(direction_line != Point3df(0,0,0));
262  normal_tri = cross(s1 - s2,s1 - s3);
263  normal_tri.normalize();
264 
265  //mesh triangle normal
266  a0 = normal_tri[0];
267  b0 = normal_tri[1];
268  c0 = normal_tri[2];
269 
270  //one node of the triangle
271  xb = s1[0];
272  yb = s1[1];
273  zb = s1[2];
274 
275  //direction of the line
276  a = direction_line[0];
277  b = direction_line[1];
278  c = direction_line[2];
279 
280  //coord of the pt defining the line
281  xm = pt[0];
282  ym = pt[1];
283  zm = pt[2];
284 
285  pt_intersection[0] = (c * c0 * xm +
286  a * c0 * zb -
287  a * b0 * ym +
288  a * a0 * xb +
289  a * b0 * yb +
290  b0 * b * xm -
291  c0 * a * zm) /
292  ( c * c0 + a * a0 + b0 * b);
293 
294 
295  pt_intersection[1] = ( b * c0 * zb +
296  b * a0 * xb +
297  b * b0 * yb -
298  b * xm * a0 +
299  c * ym * c0 +
300  ym * a * a0 -
301  b * zm * c0 ) /
302  (c * c0 + a * a0 + b0 * b);
303 
304  pt_intersection[2] = ( c * c0 * zb
305  - b0 * c * ym
306  + c * a0 * xb
307  + c * b0 * yb
308  + a * zm * a0
309  + b0 * b * zm
310  - c * xm * a0 ) /
311  (c * c0 + a * a0 + b0 * b);
312 
313  return (pt_intersection);
314  }
315 
316 
317  /* Distance from a point (A) to the plane defined by
318  one of its point : bipoint.first
319  its normal (norm = 1): bipoint.second */
320  inline float DistancePointToPlane(const Point3df &A, const std::pair<Point3df,Point3df> &plane )
321  {
322 
323  float d;
324  Point3df temp = plane.first - A;
325  d = fabs( temp.dot(plane.second) );
326 
327  return(d);
328  }
329 
330 
331  /* Distance from a point (A) to a line defined by
332  one of its point : bipoint.first
333  its direction vector (norm = 1): bipoint.second */
334  inline float DistancePointToLine(const Point3df &A, const std::pair<Point3df,Point3df> &line )
335  {
336 
337  float d;
338  Point3df temp = A - line.first;
339  d = cross(temp,line.second).norm();
340 
341  return(d);
342  }
343 
344  //Compute the mean normal n0 of the cc
345  // n1 fix the sens of each normal
346  // Only pt with a least MIN_NEIGH neighbours are taken into account
347  // It avoid considering extremal points whose normal estimation is bad
348  // The normal of these pt can be seen as outliers
349  inline Point3df meanNormalofCC(const std::vector<Point3df> &cc,
350  const Point3dfNeigh & neigh,
351  unsigned MIN_NEIGH)
352  {
353  Point3df PtR,n0 = Point3df(0,0,0),n1= Point3df(0,0,0),normal_plan;
354  Point3dfNeigh::const_iterator is,es = neigh.end();
355  std::vector<Point3df>::const_iterator icc,ecc = cc.end();
356  bool ok1 = false;
357  //Fix an arbitrary orientation (n1) for the normals of each pt of the cc
358  icc = cc.begin();
359  while (!ok1 && icc != ecc)
360  {
361  PtR = *icc;
362  is = neigh.find(PtR);
363  if (is != es)
364  {
365  if (is->second.size() > MIN_NEIGH)
366  {
367  ok1 = false;
368  normal_plan = NormalFromPoints(is->first,is->second);
369  }
370  }
371  else
372  std::cerr << "Pb1: cannot find the neighbours of the point "
373  << PtR <<"\n";
374  ++icc;
375  }
376  n1 = normal_plan;
377 
378  for(icc = cc.begin(); icc != ecc; ++icc)
379  {
380  PtR = *icc;
381  is = neigh.find(PtR);
382  if (is != es)
383  {
384  if (is->second.size() > MIN_NEIGH)
385  normal_plan = NormalFromPoints(is->first,is->second);
386  }
387  else
388  std::cerr << "Pb2: cannot find the neighbours of the point "
389  << PtR <<"\n";
390  normal_plan *= normal_plan.dot(n1);
391  normal_plan.normalize();
392  n0 += normal_plan ;
393  }
394 
395  n0.normalize();
396  std::cout << "CC mean normal : "<< n0 << std::endl;
397 
398  return (n0);
399  }
400 
401 
402 
404  extern const short MESHDISTANCE_FORBIDDEN;
405  extern const short MESHDISTANCE_UNREACHED;
406 
407  }
408 }
409 
410 
411 #endif
const short MESHDISTANCE_FORBIDDEN
global variable...
Definition: meshvoronoi.h:85
Texture< short > VolumeParcellation2MeshParcellation(const AimsData< short > &initVol, const AimsSurface< 3, Void > &mesh, short back)
TimeTexture< short > SulcusVolume2Texture(const AimsSurface< 3, Void > &mesh, const Texture< float > &curvtex, const AimsData< short > &bottom_vol, const AimsData< short > &surface_vol, float alpha, float dmin, int MINCC, const std::map< short, std::string > &trans, const std::set< std::string > &labels, float maxdil, float maxdil_mesh, float alpha_reg, bool connexity, const std::vector< std::list< unsigned > > &neigho)
Projection of the voxels relative to the bottom of the sulci to the edges of a triangulation.
Point3df meanNormalofCC(const std::vector< Point3df > &cc, const Point3dfNeigh &neigh, unsigned MIN_NEIGH)
Definition: projection.h:349
Point3df cross(const Point3df &A, const Point3df &B)
Definition: projection.h:234
AimsVector< T, D > & normalize()
float DistancePointToPlane(const Point3df &A, const std::pair< Point3df, Point3df > &plane)
Definition: projection.h:320
Texture< short > FirstSulciProjectionWithCurvatureMap(const AimsData< short > &ccvol, const AimsData< short > &sulcvol, const AimsSurface< 3, Void > &mesh, std::set< short > &cc_sulci_labels, std::map< short, short > &tri_sulci, std::map< unsigned, Point3dfSet > &cc_sulci_coord, std::map< Point3df, Point3df, Point3dfCompare > &initend, const Texture< float > &curvtex, float alpha, float dmin)
void NeighbourInCC(std::map< Point3df, Point3dfSet, Point3dfCompare > &neigh, const AimsData< short > &bvol, const AimsData< short > &svol, Connectivity::Type connectivity, short label, const unsigned size_neigh, unsigned max_points=50)
Give the neighbours of a point The neighbours have the same label as the point.
void AffinePointProjection(Texture< short > &tex, const AimsSurface< 3, Void > &mesh, const std::map< Point3df, Point3df, Point3dfCompare > &initend, const Point3dfSet &setpointi, short label, float dmin, float alpha_reg)
AimsData< T > doit(AimsData< T > &, AimsData< T > *wi=NULL)
Eigen system resolution function.
std::pair< Point3df, Point3df > NormalFromPoints(const Point3dfSet &init_point)
Definition: projection.h:153
float DistancePointToLine(const Point3df &A, const std::pair< Point3df, Point3df > &line)
Definition: projection.h:334
iterator end()
void sort(AimsData< T > &eigenvectors, AimsData< T > &eigenvalues, AimsData< T > *wi=NULL)
Sort the eigenvectors and eigenvalues in decreasing order.
T dot(const AimsVector< T, D > &other) const
float norm() const
#define ASSERT(EX)
Texture< short > FirstSulciProjectionWithInterpolationPlane(const AimsData< short > &ccvol, const AimsData< short > &sulcvol, const AimsSurface< 3, Void > &mesh, std::set< short > &cc_sulci_labels, std::map< short, short > &tri_sulci, const std::map< Point3df, Point3dfSet, Point3dfCompare > &neigh, std::map< unsigned, Point3dfSet > &cc_sulci_coord, std::map< Point3df, Point3df, Point3dfCompare > &initend, float demin, float dpmin, short label_insula_left, short label_insula_right, const Point3df &CA)
Point3df IntersectionPointOfPlanAndLine(const Point3df &s1, const Point3df &s2, const Point3df &s3, const Point3df &direction_line, const Point3df &pt)
Definition: projection.h:251
const short MESHDISTANCE_UNREACHED
Definition: meshvoronoi.h:86
AimsData< short > MeshParcellation2Volume(const AimsData< short > &initVol, const Texture< short > &tex, const AimsSurface< 3, Void > &mesh, short val_domain, short back)