A.I.M.S algorithms


primalSketchUtil.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_PRIMALSKETCH_PRIMALSKETCHUTIL_H
36 #define AIMS_PRIMALSKETCH_PRIMALSKETCHUTIL_H
37 
38 #include <cstdlib>
41 #include <aims/bucket/bucket.h>
42 #include <aims/graph/graphmanip.h>
43 #include <aims/mesh/surfacegen.h>
45 #include <unistd.h>
46 #include <graph/graph/graph.h>
47 #include <set>
48 #include <map>
49 #include <list>
50 #include <utility>
51 #include <math.h>
52 
53 namespace aims
54 {
55  // Functions that deal with primal sketches
56 
57  //------------------------------------- Declarations ----------------------------------------------------------------
58 
59  template<typename T> AimsData<short> *GetSSBlobImage(PrimalSketch<AimsData<T>, AimsData<T> > *sketch);
60  template<int D, class T> TimeTexture<float> GetSSBlobTexture(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch);
61  template<int D, class T> AimsSurfaceTriangle *GetSSBlobMesh(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch);
62 
63  template<typename T> AimsBucket<Void> *GetSSBlobBucket(PrimalSketch<AimsData<T>, AimsData<T> > *sketch);
64  template<int D, class T> AimsSurfaceTriangle *GetSSBlobTore(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch);
65 
66  template<typename Geom, typename Text> void AddBlobsToPSGraph(PrimalSketch<Geom,Text> *sketch, Graph *graph);
67  template<typename T> void AddBlobsToPSGraph(PrimalSketch<AimsData<T>, AimsData<T> > *sketch, Graph *graph);
68  template<int D, class T> void AddBlobsToPSGraph(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch, Graph *graph);
69 
70  //------------------------------------ Definitions ------------------------------------------------------------------
71 
72  template<typename T> AimsData<short> *GetSSBlobImage(PrimalSketch<AimsData<T>, AimsData<T> > *sketch)
73  {
74  typedef Point3d Site;
75 
76  AimsData<short> *image;
77  int x,y,z,sx,sy,sz, i;
78  float dx, dy, dz;
79  std::list<ScaleSpaceBlob<Site>*> ssBlobList=sketch->BlobSet();
80  std::list<ScaleSpaceBlob<Site>*>::iterator itBlob=ssBlobList.begin();
81  GreyLevelBlob<Site> *glBlob;
82  std::set<Site,ltstr_p3d<Site> > points;
83  std::set<Site,ltstr_p3d<Site> >::iterator itPoints;
84 
85  sx=sketch->scaleSpace()->GetOriginalImage().dimX();
86  sy=sketch->scaleSpace()->GetOriginalImage().dimY();
87  sz=sketch->scaleSpace()->GetOriginalImage().dimZ();
88  dx=sketch->scaleSpace()->GetOriginalImage().sizeX();
89  dy=sketch->scaleSpace()->GetOriginalImage().sizeY();
90  dz=sketch->scaleSpace()->GetOriginalImage().sizeZ();
91 
92  image=new AimsData<short>(sx, sy, sz);
93  image->setSizeXYZT(dx,dy,dz);
94 
95  for (z=0; z<sz; z++)
96  for (y=0; y<sy; y++)
97  for (x=0; x<sx; x++)
98  (*image)(x,y,z)=0;
99 
100  for (; itBlob!=ssBlobList.end(); ++itBlob)
101  {
102  glBlob=(*itBlob)->GlBlobRep();
103  points=glBlob->GetListePoints();
104  itPoints=points.begin();
105  for (; itPoints!=points.end(); ++itPoints)
106  {
107  x=(*itPoints)[0]; y=(*itPoints)[1]; z=(*itPoints)[2];
108  if ((i=(*image)(x,y,z))==0)
109  (*image)(x,y,z)=100 + ((*itBlob)->Label() * 10);
110  else
111  (*image)(x,y,z)=105 + ((*itBlob)->Label() * 10);
112  }
113  }
114 
115  return image;
116  }
117 
118 
119  //------------------------------------------------------------------------------------------
120 
121  template<int D, class T> TimeTexture<float> GetSSBlobTexture(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch)
122  {
123  typedef std::pair<Point3df,uint> Site;
124 
126  std::list<ScaleSpaceBlob<Site>*> ssBlobList=sketch->BlobSet();
127  std::list<ScaleSpaceBlob<Site>*>::iterator itBlob=ssBlobList.begin();
128  GreyLevelBlob<Site> *glBlob;
129  std::set<Site,ltstr_p3d<Site> > points;
130  std::set<Site,ltstr_p3d<Site> >::iterator itPoints;
131  int label;
132 
133  textBlob=new TexturedData<AimsSurface<D, Void>, Texture<T> >(sketch->scaleSpace()->Mesh(), &(sketch->scaleSpace()->GetOriginalImage()));
134 
135  SiteIterator<AimsSurface<D,Void> > itSite=textBlob->siteBegin();
136 
137 
138  std::set<float> scale= sketch->scaleSpace()->GetScaleList();
139 
140 
141  TimeTexture<float> tex(scale.size(),sketch->scaleSpace()->Mesh()->vertex().size());
142  for (uint i=0;i<scale.size();i++)
143  for (uint j=0;j<sketch->scaleSpace()->Mesh()->vertex().size();j++)
144  tex[i].item(j)=0;
145 
146  float sc;
147  std::set<float>::iterator it;
148  uint test=0;
149  uint i=0;
150  for (; itBlob!=ssBlobList.end(); ++itBlob)
151  {
152  label=(*itBlob)->Label();
153  glBlob=(*itBlob)->GlBlobRep();
154  points=glBlob->GetListePoints();
155  itPoints=points.begin();
156 // if (points.size() == 1){
157  std::list<GreyLevelBlob<Site>*> glBlobs=(*itBlob)->glBlobs;
158  std::list<GreyLevelBlob<Site>*>::iterator glit;
159  for (glit = glBlobs.begin() ; glit != glBlobs.end() ; glit++){
160  points=(*glit)->GetListePoints();
161  itPoints=points.begin();
162  sc = (*glit)->GetScale();
163  for (i=0,it=scale.begin();
164  *it != sc && it!=scale.end();it++,i++) {}
165 
166  for (; itPoints!=points.end(); ++itPoints)
167  {
168  tex[i].item((*itPoints).second)= (*itBlob)->GetMeasurements().t; //100+test*10;
169  }
170  }
171  test++;
172 // }
173 // uint i=0;
174 // float delta=1000.0;
175 // for (it=scale.begin();it!=scale.end();it++){
176 // if (fabs(*it-(*itBlob)->ScaleRep()) < delta)
177 // {
178 // delta=fabs(*it-(*itBlob)->ScaleRep());
179 // sc = *it;
180 // }
181 // }
182 // cout << "sc:" << sc << " ";
183  for (i=0,it=scale.begin();*it != sc && it!=scale.end();
184  it++,i++) {}
185 
186  }
187 
188  return tex;
189 // return textBlob->GetTexture();
190  }
191 
192  //------------------------------------------------------------------------------------------
193 
194  template<typename T> AimsBucket<Void> *GetSSBlobBucket(PrimalSketch<AimsData<T>, AimsData<T> > *sketch)
195  {
196  typedef Point3d Site;
197  AimsBucket<Void> *blobBck;
198  float dx, dy, dz;
199  std::list<ScaleSpaceBlob<Site>*> ssBlobList=sketch->BlobSet();
200  std::list<ScaleSpaceBlob<Site>*>::iterator itBlob=ssBlobList.begin();
201  GreyLevelBlob<Site> *glBlob;
202  std::set<Site,ltstr_p3d<Site> > points;
203  std::set<Site,ltstr_p3d<Site> >::iterator itPoints;
204  AimsBucketItem<Void> bckItem;
205 
206  dx=sketch->scaleSpace()->GetOriginalImage().sizeX();
207  dy=sketch->scaleSpace()->GetOriginalImage().sizeY();
208  dz=sketch->scaleSpace()->GetOriginalImage().sizeZ();
209 
210  blobBck=new AimsBucket<Void>;
211  blobBck->setSizeX(dx);
212  blobBck->setSizeY(dy);
213  blobBck->setSizeZ(dz);
214 
215  for (; itBlob!=ssBlobList.end(); ++itBlob)
216  {
217  glBlob=(*itBlob)->GlBlobRep();
218  points=glBlob->GetListePoints();
219  itPoints=points.begin();
220  for (; itPoints!=points.end(); ++itPoints)
221  {
222  bckItem.location()=(*itPoints);
223  (*blobBck)[(*itBlob)->Label()].push_back(bckItem);
224  }
225  }
226  return blobBck;
227  }
228 
229  //------------------------------------------------------------------------------------------
230  template<int D, class T> AimsSurfaceTriangle *GetSSBlobMesh(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch, AimsSurfaceTriangle *mesh)
231  {
232  typedef std::pair<Point3df,uint> Site; // ATTENTION : cette fonction ne marche que pour les surfaces triangulees
233  // et pas pour D quelconque. C'est un peu verole mais bon...
234  // Peut-?tre le primalSketch de surfaces devrait-il ?tre limit? aux triangulations
235  AimsSurfaceTriangle *meshBlob;
236  std::list<ScaleSpaceBlob<Site>*> ssBlobList=sketch->BlobSet();
237  std::list<ScaleSpaceBlob<Site>*>::iterator ssblobit=ssBlobList.begin();
238  std::set<Site,ltstr_p3d<Site> > points;
239 
240  std::map<int,int> tableNodes;
241  int label, index,countNode=0,countPoly=0;
242  uint count=0, count2=0;
243  Point3df coord;
244  AimsTimeSurface<3,Void> oneMesh, tempMesh;
245  std::vector<Point3df> nodes;
246  std::vector<Point3dd> pts;
247  std::vector< AimsVector<uint,3> > poly, allPoly;
248  std::vector< AimsVector<uint,3> >::iterator itPoly;
249  AimsVector<uint,3> tri;
250 
251  meshBlob=new AimsSurfaceTriangle();
252  (*mesh)[0].updateNormals();
253 
254  std::cout << "==============================" << std::endl;
255  for (; ssblobit!=ssBlobList.end(); ++ssblobit ) { // ON PARCOURT LES SSBLOBS
256  oneMesh=AimsSurfaceTriangle();
257  label=(*ssblobit)->Label();
258 
259  std::cout << "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\bssb n°" << label << "("<< count2++ <<"/" << ssBlobList.size()<<") " << std::flush;
260 
261  std::list<GreyLevelBlob<Site>*>::iterator glblobit=(*ssblobit)->glBlobs.begin();
262  std::map<float,GreyLevelBlob<Site>*> glblobmap;
263  for (;glblobit != (*ssblobit)->glBlobs.end();glblobit++){
264  glblobmap[(*glblobit)->GetScale()] = *glblobit;
265  }
266  std::map<float,GreyLevelBlob<Site>*>::iterator mapit=glblobmap.begin();
267 
268  for (count = 0;count<1 && count <glblobmap.size();count++,mapit++){ // ON PARCOURT LES GLBLOBS DU SSB
269  tempMesh=AimsSurfaceTriangle();
270  countNode=0;
271  countPoly=0;
272  tableNodes.clear();
273  points=(*mapit).second->GetListePoints();
274  points=(*ssblobit)->GlBlobRep()->GetListePoints();
275  std::set<Site,ltstr_p3d<Site> >::iterator itPoints=points.begin();
276  std::set<uint> auxpts;
277  for (; itPoints!=points.end(); ++itPoints){
278  index=(*itPoints).second;
279  auxpts.insert(index);
280  }
281  for (uint i=0; i<(*mesh)[0].polygon().size(); i++){
282  if ((auxpts.find((*mesh)[0].polygon()[i][0]) != auxpts.end()) // si le polytgone ne contient que des
283  && (auxpts.find((*mesh)[0].polygon()[i][1]) != auxpts.end()) // points du blob
284  && (auxpts.find((*mesh)[0].polygon()[i][2]) != auxpts.end())){ // on le garde
285  for (uint j=0;j<3;j++){
286  Point3df node((*mesh)[0].vertex()[(*mesh)[0].polygon()[i][j]]);
287  Point3df norm((*mesh)[0].normal()[(*mesh)[0].polygon()[i][j]]);
288  norm *= 0.01;
289  node += norm;
290  tempMesh[0].vertex().push_back(Point3df(node[0],node[1],node[2]));
291  }
292  tempMesh[0].polygon().push_back(AimsVector<uint,3>( tempMesh[0].vertex().size()-3,tempMesh[0].vertex().size()-2,tempMesh[0].vertex().size()-1));
293  }
294  }
295 
296  // JUSTE LES CYLINDRES DE LIAISON
297 // if ((*ssblobit)->TopBifurcation()->Type() != APPEAR){
298 // AimsSurfaceTriangle *msh;
299 // std::list<ScaleSpaceBlob<Site>*> topblobs = (*ssblobit)->TopBifurcation()->TopBlobs();
300 // Point3df v1 = (*mesh).vertex()[(*ssblobit)->GlBlobRep()->GetMaximum()->_node.second];
301 // Point3df n1 = (*mesh).normal()[(*ssblobit)->GlBlobRep()->GetMaximum()->_node.second];
302 // n1 *= (*mapit).second->GetScale()*3;
303 // v1 += n1;
304 // std::list<ScaleSpaceBlob<Site>*>::iterator ittb=topblobs.begin();
305 // for (; ittb !=topblobs.end();++ittb){
306 // Point3df v2 = (*mesh).vertex()[(*ittb)->GlBlobRep()->GetMaximum()->_node.second];
307 // Point3df n2 = (*mesh).normal()[(*ittb)->GlBlobRep()->GetMaximum()->_node.second];
308 // n2 *= (*ittb)->ScaleMin();
309 // v2 += n2;
310 // msh = SurfaceGenerator::cylinder(v1 ,v2 , 0.1, 0.1, 6,true,true );
311 // SurfaceManip::meshMerge( tempMesh, *msh );
312 // }
313 // }
314 //
315  SurfaceManip::meshMerge( oneMesh, tempMesh );
316  }
317 
318 
319  if (oneMesh[0].polygon().size()!=0){
320  (*meshBlob)[label]=AimsSurface<3,Void>();
321  (*meshBlob)[label].polygon()=oneMesh[0].polygon();
322  (*meshBlob)[label].vertex()=oneMesh[0].vertex();
323  (*meshBlob)[label].updateNormals();
324  }
325  else {
326 // if (true){
327  AimsSurfaceTriangle *msh;
328  Point3df auxnode = (*mesh)[0].vertex()[(*ssblobit)->GlBlobRep()->GetMaximum()->_node.second];
329  msh = SurfaceGenerator::sphere(auxnode , 0.5,6 );
330  (*meshBlob)[label]=AimsSurface<3,Void>();
331  (*meshBlob)[label].polygon()=(*msh).polygon();
332  (*meshBlob)[label].vertex()=(*msh).vertex();
333  (*meshBlob)[label].updateNormals();
334  }
335 
336  }
337  std::cout << "END" << std::endl;
338 
339  return meshBlob;
340  }
341 
342  template<typename T> AimsBucket<Void> *GetSSBlobMeshBucket(AimsSurfaceTriangle *mesh)
343  {
344  AimsBucket<Void> *blobBck;
345  AimsBucketItem<Void> bckItem;
346 
347  blobBck=new AimsBucket<Void>;
348  blobBck->setSizeX(1.0);
349  blobBck->setSizeY(1.0);
350  blobBck->setSizeZ(1.0);
351 
352  for (uint i=0; i<(*mesh).size(); i++) {
353  for (uint j=0; j<(*mesh)[i].vertex().size(); j++)
354  {
355  bckItem.location()=static_cast<Point3d>((*mesh)[i].vertex()[j]);
356  (*blobBck)[i].push_back(bckItem);
357  }
358  }
359  return blobBck;
360  }
361 
362 
363  //------------------------------------------------------------------------------------------
364  template<typename Geom, typename Text> void AddBlobsToPSGraph(PrimalSketch<Geom,Text> *sketch, Graph *graph)
365  {
366  }
367 
368 
369  //------------------------------------------------------------------------------------------
370  template<typename T> void AddBlobsToPSGraph(PrimalSketch<AimsData<T>, AimsData<T> > *sketch, Graph *graph) {
371  std::set<Vertex * > vertSet;
372  std::set<Vertex *>::iterator itVert;
373  Vertex *node;
374  BucketMap<Void> *bckMap;
375  int label;
376  aims::GraphManip manip;
378  std::vector<int> bounding_min, bounding_max;
379  std::vector<float> resolution;
380 
381  bounding_min.push_back(0);
382  bounding_min.push_back(0);
383  bounding_min.push_back(0);
384  bounding_max.push_back(sketch->scaleSpace()->GetOriginalImage().dimX());
385  bounding_max.push_back(sketch->scaleSpace()->GetOriginalImage().dimY());
386  bounding_max.push_back(sketch->scaleSpace()->GetOriginalImage().dimZ());
387  resolution.push_back(sketch->scaleSpace()->GetOriginalImage().sizeX());
388  resolution.push_back(sketch->scaleSpace()->GetOriginalImage().sizeY());
389  resolution.push_back(sketch->scaleSpace()->GetOriginalImage().sizeZ());
390 
391  graph->setProperty( "boundingbox_min", bounding_min );
392  graph->setProperty( "boundingbox_max", bounding_max );
393  graph->setProperty( "voxel_size", resolution );
394 
395  bckMap = new BucketMap<Void>( *GetSSBlobBucket(sketch) );
396  bckMap->setSizeX( resolution[0] );
397  bckMap->setSizeY( resolution[1] );
398  bckMap->setSizeZ( resolution[2] );
399  vertSet = graph->vertices();
400  for ( itVert = vertSet.begin() ; itVert != vertSet.end() ; ++itVert ) {
401  node =*itVert;
402  node->getProperty( "index", label );
403  ptrBck=carto::rc_ptr<BucketMap<Void> >( new BucketMap<Void> );
404  (*ptrBck)[0] = (*bckMap)[label];
405  ptrBck->setSizeX( resolution[0] );
406  ptrBck->setSizeY( resolution[1] );
407  ptrBck->setSizeZ( resolution[2] );
408  manip.storeAims( *graph, node, "ssblob", ptrBck );
409  node->setProperty( "ssblob_label", label );
410  }
411  }
412 
413  //------------------------------------------------------------------------------------------
414  template<int D, class T> void AddBlobsToPSGraph(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch, Graph *graph){
415  std::set<Vertex * > vertSet;
416  std::set<Vertex *>::iterator itVert;
417  Vertex *node;
418  AimsSurfaceTriangle *tore;
419  BucketMap<Void> *bckMap;
420 
421  int label;
422  aims::GraphManip manip;
425 
426  std::vector<int> bounding_min, bounding_max;
427  std::vector<float> bounding_minf, bounding_maxf;
428  std::vector<float> resolution;
429  Point3df mini, maxi;
430 
431  AimsSurfaceTriangle *mesh;
432  mesh = new AimsSurfaceTriangle();
433 
434  if (sketch->scaleSpace()->AuxMesh() != NULL){
435  mesh = sketch->scaleSpace()->AuxMesh();
436  }
437  else {
438  AimsSurface<3,Void> auxsurf(*(sketch->scaleSpace()->Mesh()));
439  (*mesh)[0] = auxsurf;
440  }
441 
442  for (uint z=0;z<3;z++){ mini[z] = 100000000.0; maxi[z] = -100000000.0; }
443 
444  for (uint i=0;i<(*mesh)[0].vertex().size();i++)
445  for (uint z=0;z<3;z++){
446  mini[z] = std::min((*mesh)[0].vertex()[i][z], mini[z]);
447  maxi[z] = std::max((*mesh)[0].vertex()[i][z], maxi[z]);
448  }
449 
450  resolution.push_back(1);
451  resolution.push_back(1);
452  resolution.push_back(1);
453  bounding_min.push_back(-1);
454  bounding_min.push_back(-1);
455  bounding_min.push_back(-1);
456  bounding_max.push_back(1);
457  bounding_max.push_back(1);
458  bounding_max.push_back(1);
459  graph->setProperty("boundingbox_min", bounding_min);
460  graph->setProperty("boundingbox_max", bounding_max);
461  graph->setProperty("voxel_size", resolution);
462 
463  tore=GetSSBlobMesh(sketch, mesh);
464  bckMap=new BucketMap<Void>();
465  bckMap->setSizeX(resolution[0]);
466  bckMap->setSizeY(resolution[1]);
467  bckMap->setSizeZ(resolution[2]);
468  vertSet=graph->vertices();
469 
470  for (itVert=vertSet.begin(); itVert!=vertSet.end(); ++itVert){
471  node=*itVert;
472  node->getProperty("index", label);
473  // Rajouter le bckMap[label] dans les attributs
475  (*ptrTore)[0]=(*tore)[label];
476 // for (uint j=0; j< (*tore)[label].vertex().size(); j++){
477 // pair<Point3d,Void> aux;
478 // aux.first = Point3d((int)(((*tore)[label].vertex()[j][0]+resolution[0]/2.0)/resolution[0]),(int)(((*tore)[label].vertex()[j][1]+resolution[1]/2.0)/resolution[1]),(int)(((*tore)[label].vertex()[j][2]+resolution[2]/2.0)/resolution[2]));
479 // (*bckMap)[label].insert(aux);
480 // }
481 
482  manip.storeAims( *graph, node, "ssblob", ptrTore );
483  ptrBck=carto::rc_ptr<BucketMap<Void> >(new BucketMap<Void>);
484 
485  (*ptrBck)[0]=(*bckMap)[label];
486  ptrBck->setSizeX(resolution[0]);
487  ptrBck->setSizeY(resolution[1]);
488  ptrBck->setSizeZ(resolution[2]);
489 // manip.storeAims( *graph, node, "ssblobbucket", ptrBck );
490 // node->setProperty("ssblobbucket_label", label);
491  node->setProperty("ssblob_label", label);
492  }
493 
494  }
495 
496  //------------------------------------------------------------------------------------------
497  template<int D, class T> AimsSurfaceTriangle *GetSSBlobTore(PrimalSketch<AimsSurface<D, Void>, Texture<T> > *sketch)
498  {
499  AimsSurfaceTriangle *blobMesh, *blobTore;
500 
501 
502  blobMesh=GetSSBlobMesh(sketch);
503 
504  // Code from here is a simple extension from the ScaleSpace::surface2Tore function
505  // by Arnaud Cachia
506 
507  unsigned i, t=0, p, t1, t2, t3;
508  //AimsSurfaceTriangle::iterator itMesh;
509  unsigned nedge = 0;
510  AimsSurfaceTriangle *msh, *tmpMesh;
511 
512  blobTore=new AimsSurfaceTriangle;
513 
514  //for (itMesh==blobMesh->begin();itMesh!=blobMesh->end();++itMesh)
515  for (t=0; t<blobMesh->size(); t++)
516  {
517  const std::vector<Point3df> vert = (*blobMesh)[t].vertex();
518  const std::vector<AimsVector<uint,3> > poly = (*blobMesh)[t].polygon();
519  std::map<std::pair<unsigned, unsigned>, unsigned> edges;
520  std::map<std::pair<unsigned, unsigned>, unsigned>::iterator ie, eie = edges.end();
521  p = poly.size();
522 
523  for( i=0; i<p; ++i )
524  {
525  t1 = poly[i][0];
526  t2 = poly[i][1];
527  t3 = poly[i][2];
528 
529  if( t1 < t2 )
530  ++edges[ std::pair<unsigned, unsigned>( t1, t2 ) ];
531  else
532  ++edges[ std::pair<unsigned, unsigned>( t2, t1 ) ];
533  if( t1 < t3 )
534  ++edges[ std::pair<unsigned, unsigned>( t1, t3 ) ];
535  else
536  ++edges[ std::pair<unsigned, unsigned>( t3, t1 ) ];
537  if( t2 < t3 )
538  ++edges[ std::pair<unsigned, unsigned>( t2, t3 ) ];
539  else
540  ++edges[ std::pair<unsigned, unsigned>( t3, t2 ) ];
541  }
542 
543  tmpMesh=new AimsSurfaceTriangle;
544  for( ie=edges.begin(); ie!=eie; ++ie)
545  if( (*ie).second == 1 )
546  {
547  ++nedge;
548  msh = SurfaceGenerator::cylinder( vert[(*ie).first.first], vert[(*ie).first.second], 0.3, 0.3, 3, false );
549  SurfaceManip::meshMerge( *tmpMesh, *msh );
550  delete msh;
551  }
552  (*blobTore)[t]=(*tmpMesh)[0];
553  delete tmpMesh;
554  }
555 
556  return( blobTore );
557  }
558 
559 }
560 
561 #endif
float min(float x, float y)
Definition: thickness.h:105
void setSizeY(float sizey)
float max(float x, float y)
Definition: thickness.h:97
static AimsSurfaceTriangle * cylinder(const carto::GenericObject &params)
AIMSDATA_API AimsTimeSurface< 3, Void > AimsSurfaceTriangle
void setSizeX(float sizex)
const std::set< Vertex * > & vertices() const
SiteIterator< Geom > siteBegin()
const std::vector< AimsVector< uint, D > > & polygon() const
const std::vector< Point3df > & vertex() const
const AimsVector< short, 3 > & location() const
unsigned int uint
static void meshMerge(AimsTimeSurface< D, T > &dst, const AimsTimeSurface< D, T > &add)
void setSizeXYZT(float sizex=1.0f, float sizey=1.0f, float sizez=1.0f, float sizet=1.0f)
void setSizeZ(float sizez)
AIMSDATA_API float norm(const Tensor &thing)
const std::vector< AimsVector< uint, D > > & polygon() const
static AimsSurfaceTriangle * sphere(const carto::GenericObject &params)
static void storeAims(Graph &graph, GraphObject *vertex, const std::string &attribute, carto::rc_ptr< T > obj)