aimsdata 6.0.0
Neuroimaging data handling
surfacemanip_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#ifndef AIMS_SURFACE_SURFACEOPERATION_D_H
35#define AIMS_SURFACE_SURFACEOPERATION_D_H
36
37#include <aims/mesh/surfacemanip.h>
38#include <aims/mesh/texture.h>
39#include <aims/resampling/motion.h>
40
41namespace aims
42{
43
44 template<int D, class T> std::vector<std::set<uint> >
45 SurfaceManip::surfaceNeighbours( const AimsSurface<D,T> & surf )
46 {
47 const std::vector< AimsVector<uint,D> > & poly = surf.polygon();
48 uint n = poly.size();
49 std::vector<std::set<uint> > neigh( surf.vertex().size() );
50
51 for ( uint i=0; i<n; ++i )
52 for ( uint j=0; j<D; ++j )
53 for ( uint k=0; k<D; ++k )
54 if ( j != k )
55 neigh[poly[i][j]].insert( poly[i][k] );
56
57 return neigh;
58 }
59
60 template<int D, class T> std::vector<std::set<uint> >
61 SurfaceManip::surfaceNeighbours( const AimsTimeSurface<D,T> & surf )
62 {
63 const std::vector< AimsVector<uint,D> > & poly = surf.polygon();
64 uint n = poly.size();
65 std::vector<std::set<uint> > neigh( surf.vertex().size() );
66
67 for ( uint i=0; i<n; ++i )
68 for ( uint j=0; j<D; ++j )
69 for ( uint k=0; k<D; ++k )
70 if ( j != k )
71 neigh[poly[i][j]].insert( poly[i][k] );
72
73 return neigh;
74 }
75
76 template<int D, class T> std::vector<std::set<uint> >
77 SurfaceManip::surfaceNeighbours2ndOrder( const AimsTimeSurface<D,T> & surf )
78 {
79 const std::vector< AimsVector<uint,D> > & poly = surf.polygon();
80 uint n = poly.size();
81 std::vector<std::set<uint> > neigh( surf.vertex().size() );
82 std::vector<std::set<uint> > neigh_final( surf.vertex().size() );
83
84 for ( uint i=0; i<n; ++i )
85 for ( uint j=0; j<D; ++j )
86 for ( uint k=0; k<D; ++k )
87 if ( j != k )
88 neigh[poly[i][j]].insert( poly[i][k] );
89
90 for ( uint i=0; i<neigh.size() ; ++i)
91 {
92 std::set<uint> voisin=neigh[i];
93 std::set<uint> voisin_final=neigh[i];
94 std::set<uint>::iterator voisin_pt=voisin.begin();
95 for ( ; voisin_pt!=voisin.end(); ++voisin_pt)
96 {
97 uint j=*voisin_pt;
98 std::set<uint> vois2=neigh[j];
99 std::set<uint>::iterator v_pt=vois2.begin();
100 for ( ; v_pt!=vois2.end(); ++v_pt)
101 voisin_final.insert(*v_pt);
102 }
103 neigh_final[i]=voisin_final;
104 }
105
106 return neigh_final;
107 }
108
109
110 template<int D, class T> std::vector<std::map<uint, float> >
111 SurfaceManip::surfaceNeighbourDistance( const AimsTimeSurface<D,T> & surf )
112 {
113 const std::vector< AimsVector<uint,D> > & poly = surf.polygon();
114 uint n = poly.size();
115 const std::vector<Point3df> & vert = surf.vertex();
116 std::vector<std::map<uint, float> > neighDist(surf.vertex().size());
117 Point3df diff;
118 float dist;
119
120 for ( uint i=0; i<n; ++i )
121 for ( uint j=0; j<D; ++j )
122 for ( uint k=0; k<D; ++k )
123 if ( j != k )
124 {
125 diff=vert[poly[i][j]] - vert[poly[i][k]]; // COURANT - VOISIN
126 dist = diff.norm(); //DISTANCE ENTRE LES DEUX;
127 neighDist[poly[i][j]][poly[i][k]] = dist;
128 }
129
130 return neighDist;
131 }
132
133 template<int D, class T>
134 void SurfaceManip::invertSurfacePolygons( AimsTimeSurface<D,T> & surface )
135 {
136 typename AimsTimeSurface<D,T>::iterator is, fs=surface.end();
137 typename std::vector< AimsVector<uint,D> >::iterator ip, fp;
138 uint n = D/2;
139 for ( is=surface.begin(); is!=fs; ++is )
140 {
141 AimsSurface<D, T> & surf = is->second;
142 std::vector< AimsVector<uint,D> > & poly = surf.polygon();
143
144 for ( ip=poly.begin(), fp=poly.end(); ip!=fp; ++ip )
145 {
146 AimsVector<uint,D> & p = *ip;
147 for ( uint i=0; i<n; ++i )
148 {
149 uint vert = p[i];
150 p[i] = p[D-i-1];
151 p[D-i-1] = vert;
152 }
153 }
154 }
155
156 carto::Object mat;
157 if( surface.header().hasProperty( "material" ) )
158 mat = surface.header().getProperty( "material" );
159 else
160 mat = carto::Object::value( carto::Dictionary() );
161 std::string fface = "counterclockwise";
162 if( mat->hasProperty( "front_face" ) )
163 {
164 std::string oldface = mat->getProperty( "front_face" )->getString();
165 if( oldface == "counterclockwise" )
166 fface = "clockwise";
167 }
168 mat->setProperty( "front_face", fface );
169 surface.header().setProperty( "material", mat );
170
171 }
172
173
174 template<int D, class T>
175 void SurfaceManip::meshMerge( AimsTimeSurface<D,T> & dst,
176 const AimsTimeSurface<D,T> & add )
177 {
178 typename AimsTimeSurface<D, T>::const_iterator is, es = add.end();
179 dst.header().copyProperties( carto::Object::reference( add.header() ) );
180
181 for( is=add.begin(); is!=es; ++is )
182 {
183 int t = is->first;
184 const AimsSurface<D, T> & imesh = is->second;
185 AimsSurface<D, T> & omesh = dst[t];
186 std::vector<Point3df> & vert = omesh.vertex();
187 std::vector<Point3df> & norm = omesh.normal();
188 std::vector<AimsVector<uint,D> > & poly = omesh.polygon();
189 const std::vector<Point3df> & vert2 = imesh.vertex();
190 const std::vector<Point3df> & norm2 = imesh.normal();
191 const std::vector<AimsVector<uint,D> > & poly2 = imesh.polygon();
192 unsigned i, j, n = vert.size(), m = vert2.size(), p = poly2.size();
193
194 // copy vertices & normals
195 for ( i=0; i<m; ++i )
196 vert.push_back( vert2[i] );
197 for ( i=0, m=norm2.size(); i<m; ++i )
198 norm.push_back( norm2[i] );
199
200 // translate & copy polygons
201 for ( i=0; i<p; ++i )
202 {
203 AimsVector<uint,D> pol;
204 for ( j=0; j<D; ++j )
205 pol[j] = poly2[i][j] + n;
206 poly.push_back( pol );
207 }
208 if( t == 0 )
209 dst.header().setProperty( "vertex_number", int( n + m ) );
210 }
211 }
212
213 template<int D, class T>
214 void SurfaceManip::meshMerge( AimsTimeSurface<D,T> & dst,
215 const std::list<AimsTimeSurface<D,T> > & src )
216 {
217 dst.erase();
218 typename std::list<AimsTimeSurface<D,T> >::const_iterator
219 im, em = src.end();
220 for ( im=src.begin(); im!=em; ++im )
221 meshMerge( dst, *im );
222 }
223
224
225 template<int D, class T>
226 void SurfaceManip::meshTransform( AimsTimeSurface<D,T> & mesh,
227 const Motion & motion )
228 {
229 typename AimsTimeSurface<D,T>::iterator itime, tend = mesh.end();
230 float norm;
231
232 for( itime=mesh.begin(); itime!=tend; ++itime )
233 {
234 typedef std::vector< Point3df > Points;
235 Points& vertex = itime->second.vertex();
236 Points& normal = itime->second.normal();
237 unsigned nnorm = normal.size();
238 for( Points::size_type k = 0; k < vertex.size(); ++k )
239 {
240 vertex[ k ] = motion.transform( vertex[ k ] );
241 /* Denis: applying this transformation to normals is WRONG if
242 there is any scale factor ! */
243 if( k < nnorm )
244 {
245 normal[ k ] = motion.transform( normal[ k ] )
246 - motion.translation();
247 // normalize normals
248 norm = normal[k][0] * normal[k][0] + normal[k][1] * normal[k][1]
249 + normal[k][2] * normal[k][2];
250 norm = 1. / sqrt( norm );
251 normal[k][0] *= norm;
252 normal[k][1] *= norm;
253 normal[k][2] *= norm;
254 }
255 }
256 }
257 if( !motion.isDirect() )
258 invertSurfacePolygons( mesh );
259 }
260
261
262 template <int D, typename T>
263 AimsTimeSurface<D,T> *
264 SurfaceManip::meshExtract( const AimsTimeSurface<D,T> & mesh,
265 const TimeTexture<int16_t> & tex,
266 int16_t value,
267 std::vector<size_t> ** overtIndex )
268 {
269 AimsTimeSurface<D,T> *omesh = new AimsTimeSurface<D,T>;
270 typename AimsTimeSurface<D,T>::const_iterator imt, jmt, emt = mesh.end();
271 typename TimeTexture<int16_t>::const_iterator itt, jtt, ett = tex.end();
272 if( overtIndex )
273 * overtIndex = new std::vector<size_t>;
274 jmt = mesh.begin();
275 imt = jmt;
276 jtt = tex.begin();
277 itt = jtt;
278 int ti = 0;
279 while( jmt != emt || jtt != ett )
280 {
281 if( jmt != emt )
282 imt = jmt;
283 if( jtt != ett )
284 itt = jtt;
285 AimsSurface<D,T> & mh = (*omesh)[ ti ];
286 std::vector<Point3df> & overt = mh.vertex();
287 std::vector<Point3df> & onorm = mh.normal();
288 std::vector<AimsVector<uint, D> > & opoly = mh.polygon();
289 const std::vector<Point3df> & vert = imt->second.vertex();
290 const std::vector<Point3df> & norm = imt->second.normal();
291 const std::vector<AimsVector<uint, D> > & poly = imt->second.polygon();
292 typename std::vector<AimsVector<uint, D> >::const_iterator
293 im, em = poly.end();
294 int i = 0, k;
295 const Texture<int16_t> & tx = itt->second;
296 std::map<uint, uint> vmap;
297 typename std::map<uint, uint>::iterator iv, ev = vmap.end();
298 unsigned nt;
299 bool absmaj, hasdist;
300 int16_t y, mval = 0;
301 Point3df g;
302 float dist, dist2;
303 for( im=poly.begin(); im != em; ++im )
304 {
305 std::map<int16_t, unsigned> vals;
306 nt = 0;
307 absmaj = false;
308 g = Point3df( 0, 0, 0 );
309 // let's have a vote
310 for( k=0; k<D; ++k )
311 {
312 g += vert[ (*im)[k] ];
313 y = tx[ (*im)[k] ];
314 unsigned & x = vals[y];
315 ++x;
316 if( x > nt )
317 {
318 nt = x;
319 mval = y;
320 absmaj = true;
321 }
322 else if( x == nt )
323 absmaj = false;
324 }
325 if( vals[ value ] == nt )
326 {
327 // 'value' is amongst the most taken value
328 if( !absmaj )
329 {
330 // take nearest vertex from barycenter with value in majority
331 g /= D;
332 dist = 0;
333 hasdist = false;
334 for( k=0; k<D; ++k )
335 {
336 y = tx[ (*im)[k] ];
337 if( vals[ y ] == nt )
338 {
339 dist2 = ( vert[ (*im)[k] ] - g ).norm2();
340 if( !hasdist || dist2 < dist )
341 {
342 dist = dist2;
343 mval = y;
344 }
345 }
346 }
347 if( mval == value )
348 absmaj = true;
349 }
350 // else absmaj == true: value has absolute majority
351 }
352 else
353 // value is in minority
354 absmaj = false;
355 if( absmaj )
356 { // copy polygon
357 AimsVector<uint, D> p;
358 for( k=0; k<D; ++k )
359 {
360 iv = vmap.find( (*im)[k] );
361 if( iv == ev )
362 {
363 vmap[ (*im)[k] ] = i;
364 p[k] = i;
365 overt.push_back( vert[ (*im)[k] ] );
366 onorm.push_back( norm[ (*im)[k] ] );
367 if( overtIndex )
368 (** overtIndex).push_back((*im)[k]);
369 ++i;
370 }
371 else
372 p[k] = iv->second;
373 }
374 opoly.push_back( p );
375 }
376 }
377 ++ti;
378 if( jmt != emt )
379 ++jmt;
380 if( jtt != ett )
381 ++jtt;
382 }
383 return omesh;
384 }
385
386
387 template <int D, typename T>
388 TimeTexture<float>*
389 SurfaceManip::meshDensity( const AimsTimeSurface<D,T> & mesh,
390 bool asDistance )
391 {
392 TimeTexture<float> *tex = new TimeTexture<float>;
393 typename AimsTimeSurface<D,T>::const_iterator is, es = mesh.end();
394 uint i;
395
396 for( is=mesh.begin(); is!=es; ++is )
397 {
398 const AimsSurface<D,T> & surf = is->second;
399 const std::vector<Point3df> & vert = surf.vertex();
400 const std::vector<AimsVector<uint,D> > & poly = surf.polygon();
401
402 // create corresponding timestep texture
403 std::vector<float> & tx = (*tex)[ is->first ].data();
404 std::vector<unsigned> counts;
405 tx.reserve( vert.size() );
406 tx.insert( tx.end(), vert.size(), 0. );
407 counts.reserve( vert.size() );
408 counts.insert( counts.end(), vert.size(), 0 );
409
410 // record edges sizes on all polygons
411 typename std::vector<AimsVector<uint,D> >::const_iterator
412 ip, ep = poly.end(), jp;
413 float n;
414
415 for( ip=poly.begin(); ip!=ep; ++ip )
416 {
417 for( i=0; i<D; ++i )
418 {
419 // edge distance, accounted on both vertices of the edge
420 uint vi = (*ip)[i];
421 uint vj = (*ip)[ (i+1) % D ];
422 n = ( vert[ vi ] - vert[ vj ] ).norm();
423 tx[ vi ] += n;
424 tx[ vj ] += n;
425 // count
426 ++counts[ vi ];
427 ++counts[ vj ];
428 }
429 }
430
431 // now average distances, and invert them to get a density
432 std::vector<float>::iterator it, et = tx.end();
433 std::vector<unsigned>::const_iterator itc;
434 if( asDistance )
435 for( it=tx.begin(), itc=counts.begin(); it!=et; ++it, ++itc )
436 if( *itc != 0. )
437 *it /= float( *itc );
438 else
439 *it = 0.; // to avoid NaN
440 else
441 for( it=tx.begin(), itc=counts.begin(); it!=et; ++it, ++itc )
442 *it = float( *itc ) / *it;
443 }
444
445 return tex;
446 }
447
448
449 template <int D, typename T>
450 TimeTexture<float>*
451 SurfaceManip::meshEdgeLengthRatioTexture(
452 const AimsTimeSurface<D,T> & nummesh,
453 const AimsTimeSurface<D,T> & denommesh )
454 {
455 TimeTexture<float> *tex = new TimeTexture<float>;
456 typename AimsTimeSurface<D,T>::const_iterator is, es = nummesh.end(),
457 ids, eds = denommesh.end();
458 uint i;
459
460 for( is=nummesh.begin(), ids=denommesh.begin(); is!=es && ids!=eds;
461 ++is, ++ids )
462 {
463 const AimsSurface<D,T> & surf = is->second;
464 const std::vector<Point3df> & vert = surf.vertex();
465 const std::vector<AimsVector<uint,D> > & poly = surf.polygon();
466 const std::vector<Point3df> & dvert = ids->second.vertex();
467
468 // create corresponding timestep texture
469 std::vector<float> & tx = (*tex)[ is->first ].data();
470 std::vector<unsigned> counts;
471 tx.reserve( vert.size() );
472 tx.insert( tx.end(), vert.size(), 0. );
473 counts.reserve( vert.size() );
474 counts.insert( counts.end(), vert.size(), 0 );
475
476 // record edges sizes on all polygons
477 typename std::vector<AimsVector<uint,D> >::const_iterator
478 ip, ep = poly.end(), jp;
479 float n, dn;
480
481 for( ip=poly.begin(); ip!=ep; ++ip )
482 {
483 for( i=0; i<D; ++i )
484 {
485 // edge distance, accounted on both vertices of the edge
486 uint vi = (*ip)[i];
487 uint vj = (*ip)[ (i+1) % D ];
488 n = ( vert[ vi ] - vert[ vj ] ).norm();
489 dn = ( dvert[ vi ] - dvert[ vj ] ).norm();
490 if( dn == 0. )
491 n = 1.; // what to do else ? (leave inf ?)
492 else
493 n /= dn;
494 tx[ vi ] = std::max( tx[ vi ], n );
495 tx[ vj ] = std::max( tx[ vj ], n );
496 // count
497 ++counts[ vi ];
498 ++counts[ vj ];
499 }
500 }
501
502 /*
503 // now average distances
504 std::vector<float>::iterator it, et = tx.end();
505 std::vector<unsigned>::const_iterator itc;
506 for( it=tx.begin(), itc=counts.begin(); it!=et; ++it, ++itc )
507 if( *itc != 0. )
508 *it /= float( *itc );
509 else
510 *it = 0.; // to avoid NaN
511 */
512 }
513
514 return tex;
515 }
516
517
518 namespace internal
519 {
520 inline uint take_mid_point_and_insert(
521 std::map<std::pair<uint,uint>, uint> & segments,
522 uint i, uint j,
523 const std::vector<Point3df> & vert,
524 std::vector<Point3df> & overt )
525 {
526 std::pair<uint,uint> seg = std::make_pair( i, j );
527 std::map<std::pair<uint,uint>, uint>::const_iterator
528 iseg = segments.find( seg );
529 uint iout;
530 if( iseg == segments.end() )
531 {
532 Point3df p = ( vert[i] + vert[j] ) / 2;
533 iout = overt.size();
534 overt.push_back( p );
535 segments[ seg ] = iout;
536 }
537 else
538 iout = iseg->second;
539 return iout;
540 }
541
542 }
543
544
545 template <typename T>
546 AimsSurface<2,Void> *
547 SurfaceManip::meshTextureBoundary( const AimsSurface<3,Void> & mesh,
548 const Texture<T> & tex, T region )
549 {
550 const std::vector<Point3df> & vert = mesh.vertex();
551 const std::vector<AimsVector<uint,3> > & poly = mesh.polygon();
552 AimsSurface<2,Void> *outmesh = new AimsSurface<2,Void>;
553 std::map<std::pair<uint,uint>, uint> segments;
554 std::vector<AimsVector<uint,3> >::const_iterator ip, ep = poly.end();
555 uint i, j, k;
556 bool ii, jj, kk;
557 std::vector<Point3df> & overt = outmesh->vertex();
558 std::vector<AimsVector<uint,2> > & opoly = outmesh->polygon();
559
560 for( ip=poly.begin(); ip!=ep; ++ip )
561 {
562 bool done = false;
563 i = (*ip)[0];
564 j = (*ip)[1];
565 k = (*ip)[2];
566 // order indices
567 if( i > j )
568 {
569 i = j;
570 j = (*ip)[0];
571 }
572 if( i > k )
573 {
574 k = i;
575 i = (*ip)[2];
576 }
577 if( j > k )
578 {
579 uint w = k;
580 k = j;
581 j = w;
582 }
583 // find boundaries in triangle
584 if( region >= 0 )
585 {
586 if( tex[i] == region )
587 ii = true;
588 else
589 ii = false;
590 if( tex[j] == region )
591 jj = true;
592 else
593 jj = false;
594 if( tex[k] == region )
595 kk = true;
596 else
597 kk = false;
598 }
599 else
600 {
601 if( tex[i] != tex[j] && tex[i] != tex[k] && tex[j] != tex[k] )
602 {
603 // all different
604 Point3df center = ( vert[i] + vert[j] + vert[k] ) / 3.;
605 uint iop = overt.size();
606 overt.push_back( center );
607 uint iout = internal::take_mid_point_and_insert( segments, i, j, vert,
608 overt );
609 opoly.push_back( AimsVector<uint,2>( iout, iop ) );
610 iout = internal::take_mid_point_and_insert( segments, i, k, vert,
611 overt );
612 opoly.push_back( AimsVector<uint,2>( iout, iop ) );
613 iout = internal::take_mid_point_and_insert( segments, j, k, vert,
614 overt );
615 opoly.push_back( AimsVector<uint,2>( iout, iop ) );
616 done = true;
617 }
618 else
619 {
620 ii = false;
621 if( tex[i] == tex[j] )
622 jj = false;
623 else
624 jj = true;
625 if( tex[i] == tex[k] )
626 kk = false;
627 else
628 kk = true;
629 }
630 }
631
632 if( !done && ( ii != jj || ii != kk ) )
633 {
634 std::vector<uint> ipts;
635 if( ii != jj )
636 ipts.push_back( internal::take_mid_point_and_insert( segments, i, j,
637 vert, overt ) );
638 if( ii != kk )
639 ipts.push_back( internal::take_mid_point_and_insert( segments, i, k,
640 vert, overt ) );
641 if( jj != kk )
642 ipts.push_back( internal::take_mid_point_and_insert( segments, j, k,
643 vert, overt ) );
644 opoly.push_back( AimsVector<uint,2>( ipts[0], ipts[1] ) );
645 }
646 }
647 return outmesh;
648 }
649
650
651 template <typename T>
652 AimsTimeSurface<2,Void> *
653 SurfaceManip::meshTextureBoundary( const AimsSurfaceTriangle & mesh,
654 const TimeTexture<T> & tex, T region )
655 {
656 AimsSurfaceTriangle::const_iterator imesh, jmesh, emesh = mesh.end();
657 imesh = mesh.begin();
658 jmesh = imesh;
659 if( jmesh != emesh )
660 ++jmesh;
661 typename TimeTexture<T>::const_iterator
662 itex = tex.begin(), jtex=itex, etex = tex.end();
663 if( jtex != etex )
664 ++jtex;
665 AimsTimeSurface<2,Void> *outmesh = new AimsTimeSurface<2,Void>;
666 while( imesh != emesh && itex != etex )
667 {
668 unsigned mtimestep = imesh->first;
669 unsigned ttimestep = itex->first;
670 AimsSurface<2, Void> *tmesh = meshTextureBoundary( imesh->second,
671 itex->second,
672 region );
673 (*outmesh)[ std::max( mtimestep, ttimestep ) ] = *tmesh;
674 delete tmesh;
675 if( jmesh == emesh && jtex == etex )
676 break;
677 if( jmesh != emesh && ( jtex == etex || jmesh->first < jtex->first ) )
678 {
679 imesh = jmesh;
680 ++jmesh;
681 }
682 else if( jtex != etex
683 && ( jmesh == emesh || jmesh->first > jtex->first ) )
684 {
685 itex = jtex;
686 ++jtex;
687 }
688 else
689 {
690 if( jmesh != emesh )
691 {
692 imesh = jmesh;
693 ++jmesh;
694 }
695 if( jtex != etex )
696 {
697 itex = jtex;
698 ++jtex;
699 }
700 }
701 }
702 return outmesh;
703 }
704
705
706 template <int D, typename T>
707 void SurfaceManip::sortPolygonsAlongDirection( AimsTimeSurface<D, T> & mesh,
708 int timestep,
709 const Point3df & direction )
710 {
711 const std::vector<Point3df> & vert = mesh[timestep].vertex();
712 std::vector<AimsVector<uint, D> > & poly = mesh[timestep].polygon();
713 std::multimap<float, AimsVector<uint, D> > sorted;
714 typename std::vector<AimsVector<uint, D> >::iterator
715 ip, ep = poly.end();
716 typename std::multimap<float, AimsVector<uint, D> >::iterator
717 im, em = sorted.end();
718 int d;
719 Point3df center;
720
721 // sort using a map
722 for( ip=poly.begin(); ip!=ep; ++ip )
723 {
724 center = vert[(*ip)[0]];
725 for( d=1; d<D; ++d )
726 center += vert[(*ip)[d]];
727 sorted.insert( std::make_pair( center.dot( direction ), *ip ) );
728 }
729
730 // rewrite polygons
731 for( im=sorted.begin(), ip=poly.begin(); im!=em; ++im, ++ip )
732 *ip = im->second;
733 }
734
735
736 template <int D, typename T>
737 std::pair<carto::rc_ptr< AimsTimeSurface<D, Void> >,
738 carto::rc_ptr<TimeTexture<T> > >
739 SurfaceManip::splitTexturedMesh( const AimsTimeSurface<D, T> & texmesh )
740 {
741 carto::rc_ptr< AimsTimeSurface<D, Void> > mesh(
742 new AimsTimeSurface<D, Void> );
743 carto::rc_ptr< TimeTexture<T> > tex( new TimeTexture<T> );
744
745 typename AimsTimeSurface<D, T>::const_iterator im, em = texmesh.end();
746 for( im=texmesh.begin(); im!=em; ++im )
747 {
748 AimsSurface<D> & m = (*mesh)[im->first];
749 Texture<T> & t = (*tex)[im->first];
750 m.vertex() = im->second.vertex();
751 m.normal() = im->second.normal();
752 m.polygon() = im->second.polygon();
753 t.data() = im->second.texture();
754 }
755
756 mesh->header().copyProperties(
757 carto::Object::reference( texmesh.header() ) );
758 tex->header().copyProperties(
759 carto::Object::reference( texmesh.header() ) );
760
761 return std::make_pair( mesh, tex );
762 }
763
764
765 template <int D, typename T>
766 carto::rc_ptr<AimsTimeSurface<D, T> >
767 SurfaceManip::joinTexturedMesh( const AimsTimeSurface<D, Void> & mesh,
768 const TimeTexture<T> & texture )
769 {
770 carto::rc_ptr< AimsTimeSurface<D, T> > texmesh(
771 new AimsTimeSurface<D, T> );
772
773 typename AimsTimeSurface<D, Void>::const_iterator im, em = mesh.end();
774 typename TimeTexture<T>::const_iterator it, et = texture.end();
775
776 std::set<int> timesteps;
777 for( im=mesh.begin(); im!=em; ++im )
778 timesteps.insert( im->first );
779 for( it=texture.begin(); it!=et; ++it )
780 timesteps.insert( it->first );
781 std::set<int>::iterator is, es = timesteps.end();
782
783 for( is=timesteps.begin(); is!=es; ++is )
784 {
785 im = mesh.lower_bound( *is );
786 if( im == em )
787 break;
788 it = texture.lower_bound( *is );
789
790 AimsSurface<D, T> & m = (*texmesh)[*is];
791 m.vertex() = im->second.vertex();
792 m.normal() = im->second.normal();
793 m.polygon() = im->second.polygon();
794 if( it != et )
795 m.texture() = it->second.data();
796 }
797
798 texmesh->header().copyProperties(
799 carto::Object::reference( mesh.header() ) );
800 texmesh->header().copyProperties(
801 carto::Object::reference( texture.header() ) );
802
803 return texmesh;
804 }
805
806
807}
808
809#endif
AimsSurface
The template class to manage a mesh.
Definition surface.h:69
AimsSurface::normal
const std::vector< Point3df > & normal() const
Get a const reference to the vector of normals.
Definition surface.h:98
AimsSurface::polygon
const std::vector< AimsVector< uint, D > > & polygon() const
Get a const reference to the vector of polygons.
Definition surface.h:108
AimsSurface::texture
const std::vector< T > & texture() const
Get a const reference to the vector of textures.
Definition surface.h:103
AimsSurface::vertex
const std::vector< Point3df > & vertex() const
Get a const reference to the vector of vertices.
Definition surface.h:93
AimsTimeSurface
The template class to manage a mesh with time if needed.
Definition surface.h:317
AimsTimeSurface::erase
void erase()
Clear all the meshes.
Definition surface.h:519
AimsTimeSurface::header
const aims::PythonHeader & header() const
Get the header.
Definition surface.h:332
AimsTimeSurface::const_iterator
std::map< int, AimsSurface< D, T > >::const_iterator const_iterator
Definition surface.h:321
AimsTimeSurface::polygon
const std::vector< AimsVector< uint, D > > & polygon() const
Get a const reference to the vector of polygons of the 0 surface.
Definition surface.h:366
AimsTimeSurface::iterator
std::map< int, AimsSurface< D, T > >::iterator iterator
Definition surface.h:319
AimsTimeSurface::vertex
const std::vector< Point3df > & vertex() const
Get a const reference to the vector of verteces of the surface of index 0.
Definition surface.h:348
AimsVector::dot
T dot(const AimsVector< T, D > &other) const
AimsVector::norm
float norm() const
Texture::data
const std::vector< T > & data() const
Definition texture.h:81
TimeTexture::const_iterator
std::map< int, Texture< T > >::const_iterator const_iterator
Definition texture.h:125
TimeTexture::header
const aims::PythonHeader & header() const
Get the header.
Definition texture.h:143
aims::SurfaceManip::sortPolygonsAlongDirection
static void sortPolygonsAlongDirection(AimsTimeSurface< D, T > &mesh, int timestep, const Point3df &direction)
Sort polygons along a given direction.
Definition surfacemanip_d.h:707
aims::SurfaceManip::meshExtract
static AimsTimeSurface< D, T > * meshExtract(const AimsTimeSurface< D, T > &mesh, const TimeTexture< int16_t > &tex, int16_t value, std::vector< size_t > **overtIndex=0)
Extracts a sub-mesh region from a label texture.
Definition surfacemanip_d.h:264
aims::SurfaceManip::invertSurfacePolygons
static void invertSurfacePolygons(AimsTimeSurface< D, T > &surface)
Definition surfacemanip_d.h:134
aims::SurfaceManip::meshEdgeLengthRatioTexture
static TimeTexture< float > * meshEdgeLengthRatioTexture(const AimsTimeSurface< D, T > &nummesh, const AimsTimeSurface< D, T > &denommesh)
Calculate an edge length ratio in edges of two meshes with the same topology.
Definition surfacemanip_d.h:451
aims::SurfaceManip::surfaceNeighbourDistance
static std::vector< std::map< uint, float > > surfaceNeighbourDistance(const AimsTimeSurface< D, T > &surf)
Definition surfacemanip_d.h:111
aims::SurfaceManip::meshTransform
static void meshTransform(AimsTimeSurface< D, T > &mesh, const Motion &trans)
Applies a transformation to a mesh.
Definition surfacemanip_d.h:226
aims::SurfaceManip::meshMerge
static void meshMerge(AimsTimeSurface< D, T > &dst, const AimsTimeSurface< D, T > &add)
concatenates 2 meshes into one (adds the second to the first one)
Definition surfacemanip_d.h:175
aims::SurfaceManip::splitTexturedMesh
static std::pair< carto::rc_ptr< AimsTimeSurface< D, Void > >, carto::rc_ptr< TimeTexture< T > > > splitTexturedMesh(const AimsTimeSurface< D, T > &texmesh)
Split a textured mesh into a mesh and a texture.
Definition surfacemanip_d.h:739
aims::SurfaceManip::joinTexturedMesh
static carto::rc_ptr< AimsTimeSurface< D, T > > joinTexturedMesh(const AimsTimeSurface< D, Void > &mesh, const TimeTexture< T > &texture)
Join a mesh and a texture into a textured mesh.
Definition surfacemanip_d.h:767
aims::SurfaceManip::meshDensity
static TimeTexture< float > * meshDensity(const AimsTimeSurface< D, T > &mesh, bool asDistance=false)
Calculate a mesh density: inverse of the average edges distance If asDistance is true,...
Definition surfacemanip_d.h:389
aims::SurfaceManip::surfaceNeighbours
static std::vector< std::set< uint > > surfaceNeighbours(const AimsSurface< D, T > &surf)
Definition surfacemanip_d.h:45
aims::SurfaceManip::surfaceNeighbours2ndOrder
static std::vector< std::set< uint > > surfaceNeighbours2ndOrder(const AimsTimeSurface< D, T > &surf)
Definition surfacemanip_d.h:77
aims::SurfaceManip::meshTextureBoundary
static AimsSurface< 2, Void > * meshTextureBoundary(const AimsSurface< 3, Void > &mesh, const Texture< T > &tex, T region)
Extracts the boundary of region of value <region> of the input texture, on the mesh.
Definition surfacemanip_d.h:547
carto::DictionaryInterface::copyProperties
virtual void copyProperties(Object source)
carto::Object::reference
Object reference(Object &value)
carto::Object::value
static Object value()
carto::TypedObject::getProperty
virtual bool getProperty(const std::string &, Object &) const
carto::TypedObject::setProperty
virtual void setProperty(const std::string &, Object)
carto::TypedObject::hasProperty
virtual bool hasProperty(const std::string &) const
carto::rc_ptr
soma::AffineTransformation3dBase::transform
Point3dd transform(double x, double y, double z) const
soma::AffineTransformationBase::isDirect
bool isDirect() const override
aims::internal::take_mid_point_and_insert
uint take_mid_point_and_insert(std::map< std::pair< uint, uint >, uint > &segments, uint i, uint j, const std::vector< Point3df > &vert, std::vector< Point3df > &overt)
Definition surfacemanip_d.h:520
aims
The class for EcatSino data write operation.
Definition borderfiller.h:13
carto::Dictionary
std::map< std::string, Object > Dictionary
Motion
aims::AffineTransformation3d Motion
Definition spmR.h:52
AimsSurfaceTriangle
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
Definition surface.h:595
uint
unsigned int uint
Point3df
AimsVector< float, 3 > Point3df
norm2
float norm2(const AimsVector< T, D > &v1)
norm
float norm(const AimsVector< T, D > &v1)