DGtal  1.5.beta
Display3D.ih
1 /**
2  * This program is free software: you can redistribute it and/or
3  * modify it under the terms of the GNU Lesser General Public License
4  * as published by the Free Software Foundation, either version 3 of
5  * the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file Display3D.ih
19  * @author Bertrand Kerautret (\c kerautre@loria.fr )
20  * LORIA (CNRS, UMR 7503), University of Nancy, France
21  *
22  * @date 2011/08/08
23  *
24  * Implementation of inline methods defined in Display3D.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 ///////////////////////////////////////////////////////////////////////////////
30 // IMPLEMENTATION of inline methods.
31 ///////////////////////////////////////////////////////////////////////////////
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 #include <cmath>
36 #include "DGtal/io/CDrawableWithDisplay3D.h"
37 #include "DGtal/io/Display3DFactory.h"
38 #include "DGtal/io/writers/MeshWriter.h"
39 
40 
41 //tmp
42 #include <iostream>
43 
44 //////////////////////////////////////////////////////////////////////////////
45 
46 ///////////////////////////////////////////////////////////////////////////////
47 // Implementation of inline methods //
48 
49 
50 // //
51 ///////////////////////////////////////////////////////////////////////////////
52 
53 ///////////////////////////////////////////////////////////////////////////////
54 // Implementation of inline functions and external operators //
55 
56 
57 template < typename Space ,typename KSpace >
58 inline
59 void
60 DGtal::Display3D< Space ,KSpace >::setFillColor(DGtal::Color aColor)
61 {
62  myCurrentFillColor=aColor;
63 }
64 
65 template < typename Space ,typename KSpace >
66 inline
67 void
68 DGtal::Display3D< Space ,KSpace >::setFillTransparency(unsigned char alpha)
69 {
70  myCurrentFillColor.alpha(alpha);
71 }
72 
73 template < typename Space ,typename KSpace >
74 inline
75 void
76 DGtal::Display3D< Space ,KSpace >::setLineColor(DGtal::Color aColor)
77 {
78  myCurrentLineColor=aColor;
79 }
80 
81 
82 template < typename Space ,typename KSpace >
83 inline
84 void
85 DGtal::Display3D< Space ,KSpace >::setKSpace( const KSpace & aKSpace )
86 {
87  myKSpace = aKSpace;
88  *myCellEmbedder = CellEmbedder( myKSpace );
89  *mySCellEmbedder = SCellEmbedder( myKSpace );
90 }
91 
92 
93 template < typename Space ,typename KSpace >
94 inline
95 DGtal::Color
96 DGtal::Display3D< Space ,KSpace >::getLineColor()
97 {
98  return myCurrentLineColor;
99 }
100 
101 template < typename Space ,typename KSpace >
102 inline
103 DGtal::Color
104 DGtal::Display3D< Space ,KSpace >::getFillColor()
105 {
106  return myCurrentFillColor;
107 }
108 
109 //-----------------------------------------------------------------------------
110 template < typename Space ,typename KSpace >
111 inline
112 void
113 DGtal::Display3D< Space ,KSpace >::
114 setName3d( DGtal::int32_t name )
115 {
116  myName3d = name;
117 }
118 
119 //-----------------------------------------------------------------------------
120 template < typename Space ,typename KSpace >
121 inline
122 DGtal::int32_t
123 DGtal::Display3D< Space ,KSpace >::
124 name3d() const
125 {
126  return myName3d;
127 }
128 //-----------------------------------------------------------------------------
129 template < typename Space ,typename KSpace >
130 inline
131 void
132 DGtal::Display3D< Space ,KSpace >::
133 setSelectCallback3D( SelectCallbackFct fct, void* data,
134  DGtal::int32_t min_name, DGtal::int32_t max_name )
135 {
136  mySelectCallBackFcts.insert( SelectCallbackFctStore( fct, data, min_name, max_name ) );
137 }
138 //-----------------------------------------------------------------------------
139 template < typename Space ,typename KSpace >
140 inline
141 typename DGtal::Display3D< Space ,KSpace >::SelectCallbackFct
142 DGtal::Display3D< Space ,KSpace >::
143 getSelectCallback3D( DGtal::int32_t aName, void*& data ) const
144 {
145  typedef typename std::set<SelectCallbackFctStore>::const_iterator SetConstIterator;
146  SelectCallbackFctStore seek( 0, 0, aName, aName );
147  SetConstIterator it_up = mySelectCallBackFcts.upper_bound( seek );
148  if ( ( it_up != mySelectCallBackFcts.end() ) && it_up->isSelected( aName ) )
149  {
150  data = it_up->data;
151  return it_up->fct;
152  }
153  if (it_up == mySelectCallBackFcts.begin()){
154  return 0;
155  }
156  --it_up;
157 
158  if ( ( it_up != mySelectCallBackFcts.end() ) && it_up->isSelected( aName ) )
159  {
160  data = it_up->data;
161  return it_up->fct;
162  }
163  return 0;
164 }
165 
166 
167 //-----------------------------------------------------------------------------
168 template < typename Space ,typename KSpace >
169 inline
170 std::string
171 DGtal::Display3D< Space ,KSpace >::getMode( const std::string & objectName ) const
172 {
173  ModeMapping::const_iterator itm = myModes.find( objectName );
174  return itm == myModes.end() ? "" : itm->second;
175 }
176 
177 
178 template < typename Space ,typename KSpace >
179 inline
180 DGtal::int32_t
181 DGtal::Display3D< Space ,KSpace >::createNewCubeList()
182 {
183  // looking for an empty key
184  DGtal::int32_t aKey=0;
185  bool found = myCubesMap.count(aKey) == 0;
186  do{
187  aKey++;
188  found = (myCubesMap.count(aKey) == 0) &&
189  (myQuadsMap.count(aKey) == 0);
190  }while (!found && aKey < std::numeric_limits<DGtal::int32_t>::max());
191  if (found){
192  myName3d = aKey;
193  return aKey;
194  }
195  return -1;
196 }
197 
198 template < typename Space ,typename KSpace >
199 inline
200 bool
201 DGtal::Display3D< Space ,KSpace >::deleteCubeList(const DGtal::int32_t idList)
202 {
203  return myCubesMap.erase(idList);
204 }
205 
206 
207 template < typename Space ,typename KSpace >
208 inline
209 void
210 DGtal::Display3D< Space ,KSpace >::createNewBallList(std::string str)
211 {
212  std::vector< BallD3D > v;
213  myBallSetList.push_back(v);
214  myBallSetNameList.push_back(str);
215 }
216 
217 template < typename Space ,typename KSpace >
218 inline
219 void
220 DGtal::Display3D< Space ,KSpace >::createNewLineList(std::string str)
221 {
222  std::vector< LineD3D > v;
223  myLineSetList.push_back(v);
224  myLineSetNameList.push_back(str);
225 }
226 
227  template < typename Space ,typename KSpace >
228  inline
229  DGtal::int32_t
230  DGtal::Display3D< Space ,KSpace >::createNewQuadList()
231  {
232  // looking for an empty key
233  DGtal::int32_t aKey=0;
234  bool found = myQuadsMap.count(aKey) == 0;
235  do{
236  aKey++;
237  found = (myCubesMap.count(aKey) == 0) &&
238  (myQuadsMap.count(aKey) == 0);
239  }while (!found && aKey < std::numeric_limits<DGtal::int32_t>::max());
240  if (found){
241  myName3d = aKey;
242  return aKey;
243  }
244 
245  return -1;
246  }
247 
248 
249 
250 template < typename Space ,typename KSpace >
251 inline
252 bool
253 DGtal::Display3D< Space ,KSpace >::deleteQuadList(const DGtal::int32_t idList)
254 {
255  return myQuadsMap.erase(idList);
256 }
257 
258 
259 
260 template < typename Space ,typename KSpace >
261 inline
262 void
263 DGtal::Display3D< Space ,KSpace >::createNewTriangleList(std::string str)
264 {
265  std::vector< TriangleD3D > v;
266  myTriangleSetList.push_back(v);
267  myTriangleSetNameList.push_back(str);
268 }
269 
270 template < typename Space ,typename KSpace >
271 inline
272 void
273 DGtal::Display3D< Space ,KSpace >::createNewPolygonList(std::string str)
274 {
275  std::vector< PolygonD3D > v;
276  myPolygonSetList.push_back(v);
277  myPolygonSetNameList.push_back(str);
278 }
279 
280 
281 template < typename Space ,typename KSpace >
282 inline
283 void
284 DGtal::Display3D< Space ,KSpace >::addCube(const RealPoint &center, double width)
285 {
286  //because the width in the cube parameters is actually the distance between the center and the faces
287  width = width/2;
288  updateBoundingBox(center);
289  CubeD3D v;
290  v.center = center;
291  v.color = getFillColor();
292  v.width = width;
293  v.name = name3d();
294  myCubesMap[ v.name ].push_back( v );
295 
296 }
297 
298 
299 template < typename Space ,typename KSpace >
300 inline
301 void
302 DGtal::Display3D< Space ,KSpace >::addBall(const RealPoint &center,
303  const double radius,
304  const unsigned int resolution)
305 {
306  updateBoundingBox(center);
307  BallD3D p;
308  p.center = center;
309  p.color = getFillColor();
310  p.isSigned = false;
311  p.signPos = false;
312  p.radius = radius;
313  p.resolution = resolution;
314  p.name = name3d();
315  if (myBallSetList.size()== 0)
316  createNewBallList("Ball Root");
317  (myBallSetList.at(myBallSetList.size()-1)).push_back(p);
318 }
319 
320 
321 
322 template < typename Space ,typename KSpace >
323 inline
324 void
325 DGtal::Display3D< Space ,KSpace >::addLine(const RealPoint &p1, const RealPoint &p2, const double width)
326 {
327  updateBoundingBox(p1);
328  updateBoundingBox(p2);
329  LineD3D l;
330  l.point1 = p1;
331  l.point2 = p2;
332  l.color = getLineColor();
333  l.width = width;
334  l.isSigned = false;
335  l.signPos = false;
336  l.name = name3d();
337  if (myLineSetList.size()== 0)
338  createNewLineList("Line Root");
339 
340  (myLineSetList.at(myLineSetList.size()-1)).push_back(l);
341 }
342 
343 template < typename Space ,typename KSpace >
344 inline
345 void
346 DGtal::Display3D< Space ,KSpace >::addQuad(const RealPoint &p1, const RealPoint &p2,
347  const RealPoint &p3, const RealPoint &p4)
348 {
349 
350  double u[3]; double v [3]; double n [3];
351  u[0]=p2[0]-p1[0]; u[1]=p2[1]-p1[1]; u[2]=p2[2]-p1[2];
352  v[0]=p3[0]-p1[0]; v[1]=p3[1]-p1[1]; v[2]=p3[2]-p1[2];
353  cross(n, u, v );
354  normalize(n);
355  addQuadWithNormal(p1, p2, p3, p4, RealPoint(n[0], n[1], n[2]), false);
356 }
357 template < typename Space ,typename KSpace >
358 inline
359 void
360 DGtal::Display3D< Space ,KSpace >::addQuadWithNormal(const RealPoint &p1, const RealPoint &p2,
361  const RealPoint &p3, const RealPoint &p4,
362  const RealPoint &n,
363  const bool enableReorientation,
364  const bool enableDoubleFace)
365 {
366  QuadD3D aQuad,aQuad2;
367  updateBoundingBox(p1);
368  updateBoundingBox(p2);
369  updateBoundingBox(p3);
370  updateBoundingBox(p4);
371  double u[3];
372  u[0]=p2[0]-p1[0]; u[1]=p2[1]-p1[1]; u[2]=p2[2]-p1[2];
373  double v[3];
374  v[0]=p3[0]-p2[0]; v[1]=p3[1]-p2[1]; v[2]=p3[2]-p2[2];
375  double crossprod[3];
376  cross(crossprod, u,v);
377  double epsilon = 0.0001;
378 
379  if(u[0]==0.0 && u[1]==0.0 && u[2]==0.0)
380  {
381  trace.warning()<< "Warning quad not added due to identical vertex... "<<std::endl;
382  return;
383  }
384 
385  if ((enableReorientation)
386  && ((crossprod[0]*n[0] + crossprod[1]*n[1] +crossprod[2]*n[2]) < 0.0))
387  {
388  aQuad.point1 = p1;
389  aQuad.point2 = p4;
390  aQuad.point3 = p3;
391  aQuad.point4 = p2;
392  aQuad.nx = n[0];
393  aQuad.ny = n[1];
394  aQuad.nz = n[2];
395  aQuad.color = getFillColor();
396 
397  if (enableDoubleFace)
398  {
399  aQuad2.point1 = (p1 - epsilon*n);
400  aQuad2.point2 = (p2 - epsilon*n);
401  aQuad2.point3 = (p3 - epsilon*n);
402  aQuad2.point4 = (p4 - epsilon*n);
403  aQuad2.nx = -n[0];
404  aQuad2.ny = -n[1];
405  aQuad2.nz = -n[2];
406  aQuad2.color = getFillColor();
407  }
408  }
409  else
410  { // Note: order of points is changed
411  aQuad.point1 = p1;
412  aQuad.point2 = p2;
413  aQuad.point3 = p3;
414  aQuad.point4 = p4;
415  aQuad.nx = n[0];
416  aQuad.ny = n[1];
417  aQuad.nz = n[2];
418  aQuad.color = getFillColor();
419 
420  if (enableDoubleFace)
421  {
422  aQuad2.point1 = (p1 - epsilon*n);
423  aQuad2.point2 = (p4 - epsilon*n);
424  aQuad2.point3 = (p3 - epsilon*n);
425  aQuad2.point4 = (p2 - epsilon*n);
426  aQuad2.nx = -n[0];
427  aQuad2.ny = -n[1];
428  aQuad2.nz = -n[2];
429  aQuad2.color = getFillColor();
430  }
431  }
432  aQuad.name = name3d();
433  aQuad2.name = name3d();
434 
435  myQuadsMap[ aQuad.name ] .push_back( aQuad );
436  if (enableDoubleFace)
437  myQuadsMap[ aQuad2.name ].push_back( aQuad2 );
438 }
439 
440 
441 template < typename Space ,typename KSpace >
442 inline
443 void
444 DGtal::Display3D< Space ,KSpace >::addTriangle(const RealPoint &p1, const RealPoint &p2, const RealPoint &p3)
445 {
446  updateBoundingBox(p1);
447  updateBoundingBox(p2);
448  updateBoundingBox(p3);
449 
450  TriangleD3D aTriangle;
451  double u[3]; double v [3]; double n [3];
452  u[0]=p2[0]-p1[0]; u[1]=p2[1]-p1[1]; u[2]=p2[2]-p1[2];
453  v[0]=p3[0]-p1[0]; v[1]=p3[1]-p1[1]; v[2]=p3[2]-p1[2];
454  cross(n, u, v );
455  cross(n, u, v );
456  normalize(n);
457 
458  if(u[0]==0.0 && u[1]==0.0 && u[2]==0.0)
459  {
460  trace.error()<< "Warning triangle not added due to identical vertex... "<<std::endl;
461  return;
462  }
463  aTriangle.point1 = p1;
464  aTriangle.point2 = p2;
465  aTriangle.point3 = p3;
466  aTriangle.nx = n[0];
467  aTriangle.ny = n[1];
468  aTriangle.nz = n[2];
469  aTriangle.color = getFillColor();
470  aTriangle.name = name3d();
471  if (myTriangleSetList.size()== 0)
472  createNewTriangleList("Triangle Root");
473 
474  (myTriangleSetList.at(myTriangleSetList.size()-1)).push_back(aTriangle);
475 }
476 
477 
478 template < typename Space ,typename KSpace >
479 inline
480 void
481 DGtal::Display3D< Space ,KSpace >::addPolygon(const std::vector<RealPoint> &vertices)
482 {
483 
484  ASSERT_MSG(vertices.size()>2, "Polygon must have at least two vertices");
485 
486  PolygonD3D polygon;
487  for(unsigned int i=0; i< vertices.size();i++)
488  {
489  polygon.vertices.push_back(vertices[i]);
490  updateBoundingBox(vertices[i]);
491  }
492 
493  double x1 = vertices[0][0];
494  double y1 = vertices[0][1];
495  double z1 = vertices[0][2];
496 
497  double x2 = vertices[1][0];
498  double y2 = vertices[1][1];
499  double z2 = vertices[1][2];
500 
501  double x3 = vertices[2][0];
502  double y3 = vertices[2][1];
503  double z3 = vertices[2][2];
504 
505  double u[3]; double v[3]; double n[3];
506  u[0]=x2-x1; u[1]=y2-y1; u[2]=z2-z1;
507  v[0]=x3-x1; v[1]=y3-y1; v[2]=z3-z1;
508  cross(n, u, v );
509  normalize(n);
510 
511  if(u[0]==0.0 && u[1]==0.0 && u[2]==0.0)
512  {
513  trace.error()<< "Warning polygon not added due to identical vertices... "<<std::endl;
514  return;
515  }
516 
517  polygon.nx = n[0];
518  polygon.ny = n[1];
519  polygon.nz = n[2];
520  polygon.color = getFillColor();
521  polygon.name = name3d();
522  if (myPolygonSetList.size()== 0)
523  createNewPolygonList("Polygon Root");
524 
525  (myPolygonSetList.at(myPolygonSetList.size()-1)).push_back(polygon);
526 }
527 
528 
529 
530 template < typename Space ,typename KSpace >
531 inline
532 void
533 DGtal::Display3D< Space ,KSpace >::addPrism(const RealPoint &baseQuadCenter,
534  bool xSurfel, bool ySurfel, bool zSurfel,
535  double sizeShiftFactor,
536  double sizeFactor, bool isSigned, bool aSign)
537 {
538  updateBoundingBox(baseQuadCenter);
539  double retract = 0.05*(sizeShiftFactor+myCurrentfShiftVisuPrisms);
540  double width = 0.03*(sizeShiftFactor+myCurrentfShiftVisuPrisms);
541 
542 
543  double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
544  double x5, y5, z5, x6, y6, z6, x7, y7, z7, x8, y8, z8;
545  double dx, dy, dz;
546 
547 
548  ASSERT( xSurfel || ySurfel || zSurfel );
549  boost::ignore_unused_variable_warning( xSurfel );
550 
551  if(zSurfel)
552  {
553  x1= baseQuadCenter[0]-(0.5*sizeFactor); y1= baseQuadCenter[1]-(0.5*sizeFactor); z1= baseQuadCenter[2]-0.5;
554  x2= baseQuadCenter[0]+(0.5*sizeFactor); y2= baseQuadCenter[1]-(0.5*sizeFactor); z2= baseQuadCenter[2]-0.5;
555  x3= baseQuadCenter[0]+(0.5*sizeFactor); y3= baseQuadCenter[1]+(0.5*sizeFactor); z3= baseQuadCenter[2]-0.5;
556  x4= baseQuadCenter[0]-(0.5*sizeFactor); y4= baseQuadCenter[1]+(0.5*sizeFactor); z4= baseQuadCenter[2]-0.5;
557  y1+=retract; y2+=retract; y3-=retract; y4-=retract;
558  x1+=retract; x2-=retract; x3-=retract; x4+=retract;
559  dx=0.0; dy=0.0; dz=width;
560  }
561  else if(ySurfel)
562  {
563  x1= baseQuadCenter[0]-(0.5*sizeFactor); y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-(0.5*sizeFactor);
564  x2= baseQuadCenter[0]-(0.5*sizeFactor); y2= baseQuadCenter[1]-0.5; z2= baseQuadCenter[2]+(0.5*sizeFactor);
565  x3= baseQuadCenter[0]+(0.5*sizeFactor); y3= baseQuadCenter[1]-0.5; z3= baseQuadCenter[2]+(0.5*sizeFactor);
566  x4= baseQuadCenter[0]+(0.5*sizeFactor); y4= baseQuadCenter[1]-0.5; z4= baseQuadCenter[2]-(0.5*sizeFactor);
567  x1+=retract; x2+=retract; x3-=retract; x4-=retract;
568  z1+=retract; z2-=retract; z3-=retract; z4+=retract;
569  dx=0.0; dy=width; dz=0.0;
570  }
571  else
572  {
573  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-(0.5*sizeFactor); z1= baseQuadCenter[2]-(0.5*sizeFactor);
574  x2= baseQuadCenter[0]-0.5; y2= baseQuadCenter[1]+(0.5*sizeFactor); z2= baseQuadCenter[2]-(0.5*sizeFactor);
575  x3= baseQuadCenter[0]-0.5; y3= baseQuadCenter[1]+(0.5*sizeFactor); z3= baseQuadCenter[2]+(0.5*sizeFactor);
576  x4= baseQuadCenter[0]-0.5; y4= baseQuadCenter[1]-(0.5*sizeFactor); z4= baseQuadCenter[2]+(0.5*sizeFactor);
577  y1+=retract; y2-=retract; y3-=retract; y4+=retract;
578  z1+=retract; z2+=retract; z3-=retract; z4-=retract;
579  dx=width; dy=0.0; dz=0.0;
580  }
581 
582  double xcenter= (x1+x2+x3+x4)/4.0;
583  double ycenter= (y1+y2+y3+y4)/4.0;
584  double zcenter= (z1+z2+z3+z4)/4.0;
585 
586  x5=x1-dx; y5=y1-dy; z5= z1-dz;
587  x6=x2-dx; y6=y2-dy; z6= z2-dz;
588  x7=x3-dx; y7=y3-dy; z7= z3-dz;
589  x8=x4-dx; y8=y4-dy; z8= z4-dz;
590 
591  x1=x1+dx; y1=y1+dy; z1= z1+dz;
592  x2=x2+dx; y2=y2+dy; z2= z2+dz;
593  x3=x3+dx; y3=y3+dy; z3= z3+dz;
594  x4=x4+dx; y4=y4+dy; z4= z4+dz;
595 
596  //if cell is oriented positively we retrac the upper face:
597  if(isSigned && aSign)
598  {
599  x1= 0.6*x1+0.4*xcenter; x2= 0.6*x2+0.4*xcenter; x3= 0.6*x3+0.4*xcenter; x4= 0.6*x4+0.4*xcenter;
600  y1= 0.6*y1+0.4*ycenter; y2= 0.6*y2+0.4*ycenter; y3= 0.6*y3+0.4*ycenter; y4= 0.6*y4+0.4*ycenter;
601  z1= 0.6*z1+0.4*zcenter; z2= 0.6*z2+0.4*zcenter; z3= 0.6*z3+0.4*zcenter; z4= 0.6*z4+0.4*zcenter;
602  }else if (isSigned)
603  {
604  x5= 0.6*x5+0.4*xcenter; x6= 0.6*x6+0.4*xcenter; x7= 0.6*x7+0.4*xcenter; x8= 0.6*x8+0.4*xcenter;
605  y5= 0.6*y5+0.4*ycenter; y6= 0.6*y6+0.4*ycenter; y7= 0.6*y7+0.4*ycenter; y8= 0.6*y8+0.4*ycenter;
606  z5= 0.6*z5+0.4*zcenter; z6= 0.6*z6+0.4*zcenter; z7= 0.6*z7+0.4*zcenter; z8= 0.6*z8+0.4*zcenter;
607  }
608 
609  //main up face
610  QuadD3D qFaceUp;
611  double normaleUp [3];
612  normaleUp[0] = dx!=0.0? 1.0:0.0;
613  normaleUp[1] = dy!=0.0 ? 1.0:0.0;
614  normaleUp[2] = dz!=0.0? 1.0:0.0;
615  qFaceUp.nx= normaleUp[0];
616  qFaceUp.ny= normaleUp[1];
617  qFaceUp.nz= normaleUp[2];
618  qFaceUp.point1[0]=x1; qFaceUp.point1[1]=y1; qFaceUp.point1[2]= z1;
619  qFaceUp.point2[0]=x2; qFaceUp.point2[1]=y2; qFaceUp.point2[2]= z2;
620  qFaceUp.point3[0]=x3; qFaceUp.point3[1]=y3; qFaceUp.point3[2]= z3;
621  qFaceUp.point4[0]=x4; qFaceUp.point4[1]=y4; qFaceUp.point4[2]= z4;
622  qFaceUp.color = myCurrentFillColor;
623  qFaceUp.name = name3d();
624  myPrismList.push_back(qFaceUp);
625 
626  //main down face
627  QuadD3D qFaceDown;
628  qFaceDown.nx= -normaleUp[0];
629  qFaceDown.ny= -normaleUp[1];
630  qFaceDown.nz= -normaleUp[2];
631  qFaceDown.point1[0]=x5; qFaceDown.point1[1]=y5; qFaceDown.point1[2]= z5;
632  qFaceDown.point2[0]=x8; qFaceDown.point2[1]=y8; qFaceDown.point2[2]= z8;
633  qFaceDown.point3[0]=x7; qFaceDown.point3[1]=y7; qFaceDown.point3[2]= z7;
634  qFaceDown.point4[0]=x6; qFaceDown.point4[1]=y6; qFaceDown.point4[2]= z6;
635  qFaceDown.color = myCurrentFillColor;
636  qFaceDown.name = name3d();
637  myPrismList.push_back(qFaceDown);
638 
639  //small face 1
640  QuadD3D qFace1;
641  double vF1[3]; double v1 [3]; double n1 [3];
642  vF1[0] = x2-x1; vF1[1] = y2-y1; vF1[2] = z2-z1;
643  v1[0] = x5-x1; v1[1] = y5-y1; v1[2] = z5-z1;
644  cross(n1, v1,vF1);
645  normalize(n1);
646  qFace1.nx=n1[0]; qFace1.ny=n1[1]; qFace1.nz=n1[2];
647  qFace1.point1[0]= x1; qFace1.point1[1] =y1; qFace1.point1[2]=z1;
648  qFace1.point2[0]= x5; qFace1.point2[1] =y5; qFace1.point2[2]=z5;
649  qFace1.point3[0]= x6; qFace1.point3[1] =y6; qFace1.point3[2]=z6;
650  qFace1.point4[0]= x2; qFace1.point4[1] =y2; qFace1.point4[2]=z2;
651  qFace1.color = myCurrentFillColor;
652  qFace1.name = name3d();
653  myPrismList.push_back(qFace1);
654 
655  //small face 2
656  QuadD3D qFace2;
657  double vF2[3]; double v2 [3]; double n2[3];
658  vF2[0]= x3-x2; vF2[1]= y3-y2; vF2[2]= z3-z2;
659  v2[0] = x6-x2; v2[1] = y6-y2; v2[2] = z6-z2;
660  cross(n2, v2, vF2);
661  normalize(n2);
662  qFace2.nx=n2[0]; qFace2.ny=n2[1]; qFace2.nz=n2[2];
663  qFace2.point1[0]= x2; qFace2.point1[1] =y2; qFace2.point1[2]=z2;
664  qFace2.point2[0]= x6; qFace2.point2[1] =y6; qFace2.point2[2]=z6;
665  qFace2.point3[0]= x7; qFace2.point3[1] =y7; qFace2.point3[2]=z7;
666  qFace2.point4[0]= x3; qFace2.point4[1] =y3; qFace2.point4[2]=z3;
667  qFace2.color = myCurrentFillColor;
668  qFace2.name = name3d();
669  myPrismList.push_back(qFace2);
670 
671  //small face 3
672  QuadD3D qFace3;
673  double vF3[3]; double v3 [3]; double n3[3];
674  vF3[0] = x4-x3; vF3[1] = y4-y3; vF3[2] = z4-z3;
675  v3[0] = x7-x3; v3[1] = y7-y3; v3[2] = z7-z3;
676  cross(n3, v3, vF3);
677  normalize(n3);
678  qFace3.nx=n3[0]; qFace3.ny=n3[1]; qFace3.nz=n3[2];
679  qFace3.point1[0]= x3; qFace3.point1[1] =y3; qFace3.point1[2]=z3;
680  qFace3.point2[0]= x7; qFace3.point2[1] =y7; qFace3.point2[2]=z7;
681  qFace3.point3[0]= x8; qFace3.point3[1] =y8; qFace3.point3[2]=z8;
682  qFace3.point4[0]= x4; qFace3.point4[1] =y4; qFace3.point4[2]=z4;
683  qFace3.color = myCurrentFillColor;
684  qFace3.name = name3d();
685  myPrismList.push_back(qFace3);
686 
687  //small face 4
688  QuadD3D qFace4;
689  double vF4[3]; double v4 [3]; double n4[3];
690  vF4[0] = x1-x4; vF4[1] = y1-y4; vF4[2] = z1-z4;
691  v4[0] = x8-x4; v4[1] = y8-y4; v4[2] = z8-z4;
692  cross(n4, v4, vF4);
693  normalize(n4);
694  qFace4.nx=n4[0]; qFace4.ny=n4[1]; qFace4.nz=n4[2];
695  qFace4.point1[0]= x4; qFace4.point1[1] =y4; qFace4.point1[2]=z4;
696  qFace4.point2[0]= x8; qFace4.point2[1] =y8; qFace4.point2[2]=z8;
697  qFace4.point3[0]= x5; qFace4.point3[1] =y5; qFace4.point3[2]=z5;
698  qFace4.point4[0]= x1; qFace4.point4[1] =y1; qFace4.point4[2]=z1;
699  qFace4.color = myCurrentFillColor;
700  qFace4.name = name3d();
701  myPrismList.push_back(qFace4);
702 
703 }
704 
705 template < typename Space ,typename KSpace >
706 inline
707 void
708 DGtal::Display3D< Space ,KSpace >::addQuadFromSurfelCenter(const RealPoint &baseQuadCenter,
709  bool xSurfel, bool ySurfel, bool zSurfel)
710 {
711  updateBoundingBox(baseQuadCenter);
712  ASSERT( xSurfel || ySurfel || zSurfel );
713  boost::ignore_unused_variable_warning( xSurfel );
714 
715  double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
716 
717  if(zSurfel)
718  {
719  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
720  x2= baseQuadCenter[0]+0.5; y2= baseQuadCenter[1]-0.5; z2= baseQuadCenter[2]-0.5;
721  x3= baseQuadCenter[0]+0.5; y3= baseQuadCenter[1]+0.5; z3= baseQuadCenter[2]-0.5;
722  x4= baseQuadCenter[0]-0.5; y4= baseQuadCenter[1]+0.5; z4= baseQuadCenter[2]-0.5;
723  }
724  else if(ySurfel)
725  {
726  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
727  x2= baseQuadCenter[0]-0.5; y2= baseQuadCenter[1]-0.5; z2= baseQuadCenter[2]+0.5;
728  x3= baseQuadCenter[0]+0.5; y3= baseQuadCenter[1]-0.5; z3= baseQuadCenter[2]+0.5;
729  x4= baseQuadCenter[0]+0.5; y4= baseQuadCenter[1]-0.5; z4= baseQuadCenter[2]-0.5;
730  }
731  else
732  {
733  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
734  x2= baseQuadCenter[0]-0.5; y2= baseQuadCenter[1]+0.5; z2= baseQuadCenter[2]-0.5;
735  x3= baseQuadCenter[0]-0.5; y3= baseQuadCenter[1]+0.5; z3= baseQuadCenter[2]+0.5;
736  x4= baseQuadCenter[0]-0.5; y4= baseQuadCenter[1]-0.5; z4= baseQuadCenter[2]+0.5;
737  }
738 
739  addQuad(RealPoint(x1, y1, z1), RealPoint(x2 ,y2, z2),
740  RealPoint(x3, y3, z3), RealPoint(x4, y4, z4));
741 }
742 
743 
744 template < typename Space ,typename KSpace >
745 inline
746 void
747 DGtal::Display3D< Space ,KSpace >::addQuadFromSurfelCenterWithNormal(const RealPoint &baseQuadCenter,
748  bool xSurfel, bool ySurfel, bool zSurfel,
749  const RealVector &aNormal,
750  const bool enableReorientation,
751  const bool sign,
752  const bool enableDoubleFace )
753 {
754  updateBoundingBox(baseQuadCenter);
755  ASSERT( xSurfel || ySurfel || zSurfel );
756  boost::ignore_unused_variable_warning( xSurfel );
757 
758  double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
759 
760  if(zSurfel)
761  {
762  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
763  x2= baseQuadCenter[0]+0.5; y2= baseQuadCenter[1]-0.5; z2= baseQuadCenter[2]-0.5;
764  x3= baseQuadCenter[0]+0.5; y3= baseQuadCenter[1]+0.5; z3= baseQuadCenter[2]-0.5;
765  x4= baseQuadCenter[0]-0.5; y4= baseQuadCenter[1]+0.5; z4= baseQuadCenter[2]-0.5;
766  }
767  else
768  if(ySurfel)
769  {
770  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
771  x2= baseQuadCenter[0]-0.5; y2= baseQuadCenter[1]-0.5; z2= baseQuadCenter[2]+0.5;
772  x3= baseQuadCenter[0]+0.5; y3= baseQuadCenter[1]-0.5; z3= baseQuadCenter[2]+0.5;
773  x4= baseQuadCenter[0]+0.5; y4= baseQuadCenter[1]-0.5; z4= baseQuadCenter[2]-0.5;
774  }
775  else
776  {
777  x1= baseQuadCenter[0]-0.5; y1= baseQuadCenter[1]-0.5; z1= baseQuadCenter[2]-0.5;
778  x2= baseQuadCenter[0]-0.5; y2= baseQuadCenter[1]+0.5; z2= baseQuadCenter[2]-0.5;
779  x3= baseQuadCenter[0]-0.5; y3= baseQuadCenter[1]+0.5; z3= baseQuadCenter[2]+0.5;
780  x4= baseQuadCenter[0]-0.5; y4= baseQuadCenter[1]-0.5; z4= baseQuadCenter[2]+0.5;
781  }
782 
783  //If we need to use the surfel sign and if it is negative, we
784  //reorder the points (see drawOriented... vs drawUnoriented methods
785  //in Display3DFactory).
786  if ((!enableReorientation) && (!sign))
787  addQuadWithNormal(RealPoint(x1, y1, z1),RealPoint(x4, y4, z4),
788  RealPoint(x3, y3, z3),
789  RealPoint(x2 ,y2, z2), aNormal,
790  enableReorientation,
791  enableDoubleFace);
792  else
793  addQuadWithNormal(RealPoint(x1, y1, z1), RealPoint(x2 ,y2, z2),
794  RealPoint(x3, y3, z3), RealPoint(x4, y4, z4), aNormal,
795  enableReorientation,
796  enableDoubleFace);
797 }
798 
799 
800 
801 
802 // add multiple triangular faces which form a hexagonal-based pyramid
803 template < typename Space ,typename KSpace >
804 inline
805 void
806 DGtal::Display3D< Space ,KSpace >::addCone(const RealPoint &p1, const RealPoint &p2,
807  double width)
808 {
809  updateBoundingBox(p1);
810  updateBoundingBox(p2);
811 
812  int nbPoints = 6; //the number of point on the base of the cone
813  double degree = 360 / nbPoints; // the angle between two points
814  double radian = degree * M_PI/180.0;
815  double norm; // norm of the vectors
816 
817  // A(x1,y1,z1) and B(x2,y2,z2) vector AB
818  double xab, yab, zab;
819 
820  xab = p2[0]-p1[0];
821  yab = p2[1]-p1[1];
822  zab = p2[2]-p1[2];
823 
824  norm = sqrt( xab*xab + yab*yab + zab*zab);
825  if (norm == 0) return;
826  xab /= norm;
827  yab /= norm;
828  zab /= norm;
829 
830  //take a third point M different from A and B
831  double xm = 0;
832  double ym = 0;
833  double zm = 0;
834 
835  while (p1[0] == xm || p2[0] ==xm) xm++;
836  while (p1[1] == ym || p2[1] ==ym) ym++;
837  while (p1[2] == zm || p2[2] ==zm) zm++;
838 
839  //vector AM
840  double xam, yam, zam;
841  //normal vector n = AB^AM
842  double xn, yn, zn;
843 
844  // ensure that M is not collinear to AB
845  while (true)
846  {
847  xam = xm-p1[0];
848  yam = ym-p1[1];
849  zam = zm-p1[2];
850 
851  norm = sqrt( xam*xam + yam*yam + zam*zam);
852  ASSERT( norm > 0 );
853  xam /= norm;
854  yam /= norm;
855  zam /= norm;
856 
857  xn = yab*zam - yam*zab;
858  yn = xam*zab - xab*zam;
859  zn = xab*yam - xam*yab;
860 
861  //divide n by its norm
862  norm = sqrt( xn*xn + yn*yn + zn*zn);
863  if (norm > 0) break;
864  xm += .1;
865  ym -= M_PI;
866  }
867 
868  ASSERT( norm > 0 );
869  xn /= norm;
870  yn /= norm;
871  zn /= norm;
872 
873  //the first point of the base
874  double xf, yf, zf;
875 
876  xf = p2[0] + width * xn;
877  yf = p2[1] + width * yn;
878  zf = p2[2] + width * zn;
879 
880  //two following points of the base
881  double xb1 = xf;
882  double yb1 = yf;
883  double zb1 = zf;
884  double xb2,yb2, zb2;
885 
886  //r = AB^n for the future rotation
887  double xr, yr, zr;
888 
889  createNewTriangleList("Cone");
890  for (int i =0; i < nbPoints-1; i ++)
891  {
892  // r = AB^n
893  xr = yab*zn - yn*zab;
894  yr = xn*zab - xab*zn;
895  zr = xab*yn - xn*yab;
896 
897  //rotate n by degree
898  xn = cos(radian)*xn + sin(radian)*xr;
899  yn = cos(radian)*yn + sin(radian)*yr;
900  zn = cos(radian)*zn + sin(radian)*zr;
901 
902  //divide n by its norm
903  norm = sqrt( xn*xn + yn*yn + zn*zn);
904  if (norm > 0)
905  {
906  xn = xn / norm;
907  yn = yn / norm;
908  zn = zn / norm;
909  }
910 
911  // calculate one point with the normal vector at a distance width
912  xb2 = p2[0] + width * xn;
913  yb2 = p2[1] + width * yn;
914  zb2 = p2[2] + width * zn;
915 
916  //adding the triangles associated with this point
917  addTriangle(p1,
918  RealPoint(xb1, yb1, zb1),
919  RealPoint(xb2, yb2, zb2));
920  addTriangle(p2,
921  RealPoint(xb1, yb1, zb1),
922  RealPoint(xb2, yb2, zb2));
923 
924  xb1 = xb2;
925  yb1 = yb2;
926  zb1 = zb2;
927  }
928 
929  //Last part to close the cone
930  addTriangle(p1,
931  RealPoint(xb1, yb1,zb1),
932  RealPoint(xf, yf, zf));
933  addTriangle(p2,
934  RealPoint(xb1, yb1, zb1),
935  RealPoint(xf, yf, zf));
936 }
937 
938 // add multiple triangular faces which form a hexagonal-based cylinder
939 template < typename Space ,typename KSpace >
940 inline
941 void
942 DGtal::Display3D< Space ,KSpace >::addCylinder(const RealPoint &p1, const RealPoint &p2,
943  const double width)
944 {
945  updateBoundingBox(p1);
946  updateBoundingBox(p2);
947 
948  int nbPoints = 6; //the number of point on the base of the cone
949  double degree = 360 / nbPoints; // the angle between two points
950  double radian = degree * M_PI/180.0;
951  double norm; // norm of the vectors
952 
953  // A(p1[0],p1[1],p1[2]) and B(p2[0],p2[1],p2[2]) vector AB
954  double xab, yab, zab;
955 
956  xab = p2[0]-p1[0];
957  yab = p2[1]-p1[1];
958  zab = p2[2]-p1[2];
959 
960  norm = sqrt( xab*xab + yab*yab + zab*zab);
961  assert (norm > 0);
962  xab = xab / norm;
963  yab = yab / norm;
964  zab = zab / norm;
965 
966  //take a third point M different from A and B
967  double xm = 0;
968  double ym = 0;
969  double zm = 0;
970 
971  while (p1[0] == xm || p2[0] ==xm) xm++;
972  while (p1[1] == ym || p2[1] ==ym) ym++;
973  while (p1[2] == zm || p2[2] ==zm) zm++;
974 
975  //vector AM
976  double xam, yam, zam;
977 
978  xam = xm-p1[0];
979  yam = ym-p1[1];
980  zam = zm-p1[2];
981 
982  norm = sqrt( xam*xam + yam*yam + zam*zam);
983  assert (norm > 0);
984  xam = xam / norm;
985  yam = yam / norm;
986  zam = zam / norm;
987 
988  //normal vector n = AB^AM
989  double xn, yn, zn;
990 
991  xn = yab*zam - yam*zab;
992  yn = xam*zab - xab*zam;
993  zn = xab*yam - xam*yab;
994 
995  //divide n by its norm
996  norm = sqrt( xn*xn + yn*yn + zn*zn);
997  assert (norm > 0);
998  xn = xn / norm;
999  yn = yn / norm;
1000  zn = zn / norm;
1001 
1002 
1003  //the first point of the bases
1004  double xbf, ybf, zbf;
1005  double xaf, yaf, zaf;
1006 
1007  xbf = p2[0] + width * xn;
1008  ybf = p2[1] + width * yn;
1009  zbf = p2[2] + width * zn;
1010 
1011  xaf = p1[0] + width * xn;
1012  yaf = p1[1] + width * yn;
1013  zaf = p1[2] + width * zn;
1014 
1015  //two following points of the bases
1016  double xb1 = xbf;
1017  double yb1 = ybf;
1018  double zb1 = zbf;
1019  double xb2,yb2, zb2;
1020 
1021  double xa1 = xaf;
1022  double ya1 = yaf;
1023  double za1 = zaf;
1024  double xa2,ya2, za2;
1025 
1026  //r = AB^n for the future rotation
1027  double xr, yr, zr;
1028 
1029  createNewTriangleList("Cylinder");
1030  //createNewQuadList("Cylinder");
1031  for (int i =0; i < nbPoints-1; i ++)
1032  {
1033  // r = AB^n
1034  xr = yab*zn - yn*zab;
1035  yr = xn*zab - xab*zn;
1036  zr = xab*yn - xn*yab;
1037 
1038  //rotate n by degree
1039  xn = cos(radian)*xn + sin(radian)*xr;
1040  yn = cos(radian)*yn + sin(radian)*yr;
1041  zn = cos(radian)*zn + sin(radian)*zr;
1042 
1043  //divide n by its norm (used when a bug appear)
1044  norm = sqrt( xn*xn + yn*yn + zn*zn);
1045  assert (norm > 0);
1046  xn = xn / norm;
1047  yn = yn / norm;
1048  zn = zn / norm;
1049 
1050  // calculate one point with the normal vector at a distance width
1051  xb2 = p2[0] + width * xn;
1052  yb2 = p2[1] + width * yn;
1053  zb2 = p2[2] + width * zn;
1054 
1055  xa2 = p1[0] + width * xn;
1056  ya2 = p1[1] + width * yn;
1057  za2 = p1[2] + width * zn;
1058 
1059  //adding the triangles and the quad associated with those points
1060  addTriangle(p1,
1061  RealPoint(xa1, ya1, za1),
1062  RealPoint(xa2, ya2, za2));
1063  addTriangle(p2,
1064  RealPoint(xb1, yb1, zb1),
1065  RealPoint(xb2, yb2, zb2));
1066  addQuad(RealPoint(xb1, yb1, zb1),
1067  RealPoint(xb2, yb2, zb2),
1068  RealPoint(xa2, ya2, za2),
1069  RealPoint(xa1, ya1, za1));
1070 
1071  xb1 = xb2;
1072  yb1 = yb2;
1073  zb1 = zb2;
1074 
1075  xa1 = xa2;
1076  ya1 = ya2;
1077  za1 = za2;
1078 
1079  }
1080 
1081  //Last part to close the cylinder
1082  addTriangle(p1,
1083  RealPoint(xa1, ya1, za1),
1084  RealPoint(xaf, yaf, zaf));
1085  addTriangle(p2,
1086  RealPoint(xb1, yb1, zb1),
1087  RealPoint(xbf, ybf, zbf));
1088  addQuad(RealPoint(xb1, yb1, zb1),
1089  RealPoint(xbf, ybf, zbf),
1090  RealPoint(xaf, yaf, zaf),
1091  RealPoint(xa1, ya1, za1));
1092 }
1093 
1094 
1095 
1096 template < typename Space ,typename KSpace >
1097 inline
1098 void
1099 DGtal::Display3D< Space ,KSpace >::addClippingPlane(double a, double b, double c, double d, bool drawPlane)
1100 {
1101  ClippingPlaneD3D cp;
1102  cp.a=a; cp.b=b; cp.c=c; cp.d=d;
1103  myClippingPlaneList.push_back(cp);
1104  if(drawPlane)
1105  {
1106  double x1,x2,x3,x4;
1107  double y1,y2,y3,y4;
1108  double z1,z2,z3,z4;
1109  double norm = sqrt(a*a+b*b+c*c);
1110  double dec=0.1;
1111  // Z dominant projection of the upper face
1112  if(std::abs(c)>=std::abs(b) && std::abs(c) >= std::abs(a))
1113  {
1114  x1= myBoundingPtUp[0]+a*dec/norm; y1= myBoundingPtUp[1]+b*dec/norm;
1115  z1 = c*dec/norm +(-d-a*myBoundingPtUp[0]-b*myBoundingPtUp[1])/c;
1116  x2= myBoundingPtLow[0]+a*dec/norm; y2= myBoundingPtUp[1]+b*dec/norm;
1117  z2= c*dec/norm+(-d-a*myBoundingPtLow[0]-b*myBoundingPtUp[1])/c;
1118  x3= myBoundingPtLow[0]+a*dec/norm; y3= myBoundingPtLow[1]+b*dec/norm;
1119  z3= c*dec/norm+(-d-a*myBoundingPtLow[0]-b*myBoundingPtLow[1])/c;
1120  x4= myBoundingPtUp[0]+a*dec/norm; y4= myBoundingPtLow[1]+b*dec/norm;
1121  z4= c*dec/norm+(-d-a*myBoundingPtUp[0]-b*myBoundingPtLow[1])/c;
1122  if(c>0)
1123  {
1124  addQuad(RealPoint(x4,y4,z4),RealPoint(x3,y3,z3),
1125  RealPoint(x2,y2,z2),RealPoint(x1,y1,z1) );
1126  }else
1127  {
1128  addQuad(RealPoint(x1,y1,z1), RealPoint( x2,y2,z2),
1129  RealPoint(x3,y3,z3), RealPoint(x4,y4,z4));
1130  }
1131  }// Y dominant projection of the upper face
1132  else if(std::abs(b)>=std::abs(c) && std::abs(b) >= std::abs(a))
1133  {
1134  x1= myBoundingPtUp[0]+a*dec/norm; z1= myBoundingPtUp[2]+c*dec/norm;
1135  y1= b*dec/norm +(-d-a*myBoundingPtUp[0]-c*myBoundingPtUp[2])/b;
1136  x2= myBoundingPtUp[0]+a*dec/norm; z2= myBoundingPtLow[2]+c*dec/norm;
1137  y2= b*dec/norm +(-d-a*myBoundingPtUp[0]-c*myBoundingPtLow[2])/b;
1138  x3= myBoundingPtLow[0]+a*dec/norm; z3= myBoundingPtLow[2]+c*dec/norm;
1139  y3= b*dec/norm +(-d-a*myBoundingPtLow[0]-c*myBoundingPtLow[2])/b;
1140  x4= myBoundingPtLow[0]+a*dec/norm; z4= myBoundingPtUp[2]+c*dec/norm;
1141  y4= b*dec/norm +(-d-a*myBoundingPtLow[0]-c*myBoundingPtUp[2])/b;
1142  if(b>0)
1143  {
1144  addQuad(RealPoint(x4,y4,z4),RealPoint(x3, y3,z3),
1145  RealPoint(x2,y2,z2), RealPoint(x1,y1,z1));
1146  }else
1147  {
1148  addQuad(RealPoint(x1,y1,z1), RealPoint(x2,y2,z2),
1149  RealPoint(x3,y3,z3), RealPoint(x4,y4,z4));
1150  }
1151  }// X dominant projection of the upper face
1152  else if(std::abs(a)>=std::abs(c) && std::abs(a) >= std::abs(b))
1153  {
1154  y1= myBoundingPtUp[1]+b*dec/norm; z1= myBoundingPtUp[2]+c*dec/norm;
1155  x1= a*dec/norm +(-d-b*myBoundingPtUp[1]-c*myBoundingPtUp[2])/a;
1156  y2= myBoundingPtLow[1]+b*dec/norm; z2= myBoundingPtUp[2]+c*dec/norm;
1157  x2= a*dec/norm +(-d-b*myBoundingPtLow[1]-c*myBoundingPtUp[2])/a;
1158  y3= myBoundingPtLow[1]+b*dec/norm; z3= myBoundingPtLow[2]+c*dec/norm;
1159  x3= a*dec/norm +(-d-b*myBoundingPtLow[1]-c*myBoundingPtLow[2])/a;
1160  y4= myBoundingPtUp[1]+b*dec/norm; z4= myBoundingPtLow[2]+c*dec/norm;
1161  x4= a*dec/norm +(-d-b*myBoundingPtUp[1]-c*myBoundingPtLow[2])/a;
1162 
1163  if(a>0)
1164  {
1165 
1166  addQuad(RealPoint(x4,y4,z4),
1167  RealPoint(x3,y3,z3),
1168  RealPoint(x2,y2,z2),
1169  RealPoint(x1,y1,z1));
1170  }else
1171  {
1172  addQuad( RealPoint(x1,y1,z1), RealPoint(x2,y2,z2),
1173  RealPoint(x3,y3,z3),RealPoint(x4,y4,z4));
1174  }
1175  }
1176  }
1177 }
1178 
1179 
1180 template < typename Space ,typename KSpace >
1181 inline
1182 void
1183 DGtal::Display3D< Space ,KSpace >::updateBoundingBox(const RealPoint &p)
1184 {
1185  if (myBoundingPtEmptyTag)
1186  {
1187  myBoundingPtLow[0]= p[0];
1188  myBoundingPtLow[1]= p[1];
1189  myBoundingPtLow[2]= p[2];
1190  myBoundingPtUp[0]= p[0];
1191  myBoundingPtUp[1]= p[1];
1192  myBoundingPtUp[2]= p[2];
1193  myBoundingPtEmptyTag = false;
1194  myBoundingPtChangedTag = true;
1195  }
1196  else
1197  {
1198  if ( p[ 0 ] < myBoundingPtLow[ 0 ] )
1199  {
1200  myBoundingPtLow[ 0 ] = p[ 0 ];
1201  myBoundingPtChangedTag = true;
1202  }
1203  if ( p[ 1 ] < myBoundingPtLow[ 1 ] )
1204  {
1205  myBoundingPtLow[ 1 ] = p[ 1 ];
1206  myBoundingPtChangedTag = true;
1207  }
1208 
1209  if ( p[ 2 ] < myBoundingPtLow[ 2 ] )
1210  {
1211  myBoundingPtLow[ 2 ] = p[ 2 ];
1212  myBoundingPtChangedTag = true;
1213  }
1214 
1215  if ( p[ 0 ] > myBoundingPtUp[ 0 ] )
1216  {
1217  myBoundingPtUp[ 0 ] = p[ 0 ];
1218  myBoundingPtChangedTag = true;
1219  }
1220  if ( p[ 1 ] > myBoundingPtUp[ 1 ] )
1221  {
1222  myBoundingPtUp[ 1 ] = p[ 1 ];
1223  myBoundingPtChangedTag = true;
1224  }
1225  if ( p[ 2 ] > myBoundingPtUp[ 2 ] )
1226  {
1227  myBoundingPtUp[ 2 ] = p[ 2 ];
1228  myBoundingPtChangedTag = true;
1229  }
1230  }
1231 }
1232 
1233 
1234 template < typename Space ,typename KSpace >
1235 inline
1236 typename DGtal::Display3D< Space ,KSpace >::RealPoint
1237 DGtal::Display3D< Space ,KSpace >::embed( const typename Space::Point & dp ) const
1238 {
1239  ASSERT( myEmbedder->isValid());
1240  return myEmbedder->embed(dp);
1241 }
1242 
1243 template < typename Space ,typename KSpace >
1244 inline
1245 typename DGtal::Display3D< Space ,KSpace >::RealPoint
1246 DGtal::Display3D< Space ,KSpace >::embedKS( const typename KSpace::SCell & scell ) const
1247 {
1248  ASSERT( mySCellEmbedder->isValid());
1249  return mySCellEmbedder->embed(scell);
1250 }
1251 
1252 template < typename Space ,typename KSpace >
1253 inline
1254 typename DGtal::Display3D< Space ,KSpace >::RealPoint
1255 DGtal::Display3D< Space ,KSpace >::embedKS( const TransformedPrism & aTrans ) const
1256 {
1257  ASSERT( mySCellEmbedder->isValid());
1258  return mySCellEmbedder->embed(aTrans.mySurfel);
1259 }
1260 
1261 
1262 template < typename Space ,typename KSpace >
1263 inline
1264 typename DGtal::Display3D< Space ,KSpace >::RealPoint
1265 DGtal::Display3D< Space ,KSpace >::embedK( const typename KSpace::Cell & cell ) const
1266 {
1267  ASSERT( myCellEmbedder->isValid());
1268  return myCellEmbedder->embed(cell);
1269 }
1270 
1271 template < typename Space ,typename KSpace >
1272 inline
1273 void
1274 DGtal::Display3D< Space ,KSpace >::exportToMesh(DGtal::Mesh<RealPoint> &aMesh) const
1275 {
1276  unsigned int vertexIndex=0;
1277 
1278  // Export of SurfelPrism lists (generated from KhalimskyCell or SignedKhalimskyCell (through addPrism))
1279  for (unsigned int i=0; i< myPrismList.size(); i++)
1280  {
1281  QuadD3D quad = myPrismList.at(i);
1282  RealPoint p1, p2, p3, p4;
1283  p1=quad.point1;
1284  p2=quad.point2;
1285  p3=quad.point3;
1286  p4=quad.point4;
1287  aMesh.addVertex(p1);
1288  aMesh.addVertex(p2);
1289  aMesh.addVertex(p3);
1290  aMesh.addVertex(p4);
1291  aMesh.addQuadFace(vertexIndex, vertexIndex+1, vertexIndex+2,vertexIndex+3,
1292  quad.color);
1293  vertexIndex+=4;
1294  }
1295 
1296 
1297  // Export QuadList
1298  for (typename QuadsMap::const_iterator it = myQuadsMap.begin(); it != myQuadsMap.end(); it++)
1299  {
1300  for (typename std::vector<QuadD3D>::const_iterator aQuad = it->second.begin(); aQuad!=it->second.end();aQuad ++)
1301  {
1302  RealPoint p1, p2, p3, p4;
1303  p1=aQuad->point1;
1304  p2=aQuad->point2;
1305  p3=aQuad->point3;
1306  p4=aQuad->point4;
1307  aMesh.addVertex(p4);
1308  aMesh.addVertex(p3);
1309  aMesh.addVertex(p2);
1310  aMesh.addVertex(p1);
1311  aMesh.addQuadFace(vertexIndex,vertexIndex+1, vertexIndex+2, vertexIndex+3,
1312  aQuad->color);
1313  vertexIndex+=4;
1314  }
1315  }
1316 
1317  // Export the mesh from TriangleList
1318  for ( typename std::vector<std::vector<TriangleD3D> >::const_iterator it =myTriangleSetList.begin(); it != myTriangleSetList.end(); it++)
1319  {
1320  for (typename std::vector<TriangleD3D>::const_iterator aTriangle = it->begin(); aTriangle!=it->end();aTriangle ++)
1321  {
1322  RealPoint p1, p2, p3;
1323  p1=aTriangle->point1;
1324  p2=aTriangle->point2;
1325  p3=aTriangle->point3;
1326  aMesh.addVertex(p1);
1327  aMesh.addVertex(p2);
1328  aMesh.addVertex(p3);
1329  aMesh.addTriangularFace(vertexIndex, vertexIndex+1, vertexIndex+2,
1330  aTriangle->color);
1331  vertexIndex+=3;
1332  }
1333  }
1334 
1335  // Export of cubeSet (generated from addCube)
1336  // Export CubesList
1337  for (typename CubesMap::const_iterator it = myCubesMap.begin(); it != myCubesMap.end(); it++)
1338  {
1339  for (typename std::vector<CubeD3D>::const_iterator itCube = it->second.begin(); itCube!=it->second.end(); itCube++)
1340  {
1341 
1342  CubeD3D cube = *itCube;
1343  RealPoint p1, p2, p3, p4, p5, p6, p7, p8;
1344  double width= cube.width;
1345 
1346  p1 = RealPoint(cube.center[0]-width, cube.center[1]+width, cube.center[2]+width);
1347  p2 = RealPoint(cube.center[0]+width, cube.center[1]+width, cube.center[2]+width);
1348  p3 = RealPoint(cube.center[0]+width, cube.center[1]-width, cube.center[2]+width);
1349  p4 = RealPoint(cube.center[0]-width, cube.center[1]-width, cube.center[2]+width);
1350  p5 = RealPoint(cube.center[0]-width, cube.center[1]+width, cube.center[2]-width);
1351  p6 = RealPoint(cube.center[0]+width, cube.center[1]+width, cube.center[2]-width);
1352  p7 = RealPoint(cube.center[0]+width, cube.center[1]-width, cube.center[2]-width);
1353  p8 = RealPoint(cube.center[0]-width, cube.center[1]-width, cube.center[2]-width);
1354 
1355  aMesh.addVertex(p1);
1356  aMesh.addVertex(p2);
1357  aMesh.addVertex(p3);
1358  aMesh.addVertex(p4);
1359  aMesh.addVertex(p5);
1360  aMesh.addVertex(p6);
1361  aMesh.addVertex(p7);
1362  aMesh.addVertex(p8);
1363 
1364  //z+
1365  aMesh.addQuadFace(vertexIndex, vertexIndex+3, vertexIndex+2, vertexIndex+1,
1366  cube.color);
1367  //z-
1368  aMesh.addQuadFace(vertexIndex+4, vertexIndex+5, vertexIndex+6, vertexIndex+7,
1369  cube.color);
1370  //y+
1371  aMesh.addQuadFace(vertexIndex+1, vertexIndex+2, vertexIndex+6, vertexIndex+5,
1372  cube.color);
1373  //y-
1374  aMesh.addQuadFace(vertexIndex, vertexIndex+4, vertexIndex+7, vertexIndex+3,
1375  cube.color);
1376  //x+
1377  aMesh.addQuadFace(vertexIndex, vertexIndex+1, vertexIndex+5, vertexIndex+4,
1378  cube.color);
1379  //x-
1380  aMesh.addQuadFace(vertexIndex+3, vertexIndex+7, vertexIndex+6, vertexIndex+2,
1381  cube.color);
1382 
1383  vertexIndex+=8;
1384  }
1385  }
1386 }
1387 
1388 
1389 template < typename Space ,typename KSpace >
1390 template<typename TDrawableWithDisplay3D>
1391 inline
1392 DGtal::Display3D< Space ,KSpace >&
1393 DGtal::Display3D< Space ,KSpace >::operator <<( const TDrawableWithDisplay3D & object )
1394 {
1395  // BOOST_CONCEPT_ASSERT((CDrawableWithDisplay3D< TDrawableWithDisplay3D >));
1396 
1397  DGtal::Display3DFactory<Space,KSpace>::draw(*this, object);
1398  return *this;
1399 }
1400 
1401 
1402 template < typename Space ,typename KSpace >
1403 inline
1404 std::ostream&
1405 DGtal::operator<< ( std::ostream & out,
1406  const Display3D< Space ,KSpace >& object )
1407 {
1408  object.selfDisplay ( out );
1409  return out;
1410 }
1411 
1412 
1413 template < typename Space ,typename KSpace >
1414 inline
1415 void
1416 DGtal::operator>> ( const Display3D< Space ,KSpace >&aDisplay3D, DGtal::Mesh< typename Display3D<Space , KSpace >::RealPoint > &aMesh)
1417 {
1418  aDisplay3D.exportToMesh(aMesh);
1419 }
1420 
1421 
1422 template < typename Space ,typename KSpace >
1423 inline
1424 void
1425 DGtal::operator>> ( const Display3D< Space ,KSpace >&aDisplay3D, std::string aFilename)
1426 {
1427  // exporting with a mesh containing color (parameter constructor to true):
1428  DGtal::Mesh<typename Display3D<Space , KSpace >::RealPoint> mesh(true);
1429  aDisplay3D >> mesh;
1430  trace.info() << "generating faces done." << std::endl;
1431  mesh >> aFilename;
1432  trace.info() << "file exported in file: " << aFilename << std::endl;
1433 }
1434 
1435 
1436 template < typename Space ,typename KSpace >
1437 inline
1438 void
1439 DGtal::Display3D< Space ,KSpace >::cross (double dst[3], double srcA[3], double srcB[3])
1440 {
1441  dst[0] = srcA[1]*srcB[2] - srcA[2]*srcB[1];
1442  dst[1] = srcA[2]*srcB[0] - srcA[0]*srcB[2];
1443  dst[2] = srcA[0]*srcB[1] - srcA[1]*srcB[0];
1444 }
1445 
1446 
1447 template < typename Space ,typename KSpace >
1448 inline
1449 void
1450 DGtal::Display3D< Space ,KSpace >::normalize (double vec[3])
1451 {
1452  const double squaredLen = vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2];
1453  vec[0] /= sqrt (squaredLen);
1454  vec[1] /= sqrt (squaredLen);
1455  vec[2] /= sqrt (squaredLen);
1456 }
1457 
1458 template < typename Space ,typename KSpace >
1459 inline
1460 bool
1461 DGtal::Display3D< Space ,KSpace >::isValid() const
1462 {
1463  return myEmbedder->isValid()
1464  && myCellEmbedder->isValid()
1465  && mySCellEmbedder->isValid();
1466 }
1467 
1468 template < typename Space ,typename KSpace >
1469 inline
1470 void
1471 DGtal::Display3D< Space ,KSpace >::clear()
1472 {
1473  myCubesMap.clear();
1474  myLineSetList.clear();
1475  myBallSetList.clear();
1476  myClippingPlaneList.clear();
1477  myPrismList.clear();
1478  myQuadsMap.clear();
1479  myTriangleSetList.clear();
1480  myPolygonSetList.clear();
1481  myCubeSetNameList.clear();
1482  myLineSetNameList.clear();
1483  myBallSetNameList.clear();
1484  myClippingPlaneNameList.clear();
1485  myPrismNameList.clear();
1486 
1487 myTriangleSetNameList.clear();
1488  myPolygonSetNameList.clear();
1489 
1490  //Bounding box reset
1491  myBoundingPtEmptyTag = true;
1492  for (unsigned int i=0; i< 3; i++)
1493  {
1494  myBoundingPtUp[i] = 0.0;
1495  myBoundingPtLow[i] = 0.0;
1496  }
1497 }
1498 
1499 
1500 
1501 // //
1502 ///////////////////////////////////////////////////////////////////////////////