OpenWalnut  1.3.1
WGridRegular3D.h
1 //---------------------------------------------------------------------------
2 //
3 // Project: OpenWalnut ( http://www.openwalnut.org )
4 //
5 // Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
6 // For more information see http://www.openwalnut.org/copying
7 //
8 // This file is part of OpenWalnut.
9 //
10 // OpenWalnut is free software: you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // OpenWalnut is distributed in the hope that it will be useful,
16 // but WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public License
21 // along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
22 //
23 //---------------------------------------------------------------------------
24 
25 #ifndef WGRIDREGULAR3D_H
26 #define WGRIDREGULAR3D_H
27 
28 #include <cmath>
29 #include <string>
30 #include <utility>
31 #include <vector>
32 
33 #include <boost/array.hpp>
34 #include <boost/shared_ptr.hpp>
35 
36 #include <osg/Matrix>
37 #include <osg/Vec3>
38 
39 #include "../common/exceptions/WOutOfBounds.h"
40 #include "../common/exceptions/WPreconditionNotMet.h"
41 #include "../common/math/WLinearAlgebraFunctions.h"
42 #include "../common/math/WMatrix.h"
43 #include "../common/math/linearAlgebra/WLinearAlgebra.h"
44 #include "../common/WBoundingBox.h"
45 #include "../common/WCondition.h"
46 #include "../common/WDefines.h"
47 #include "../common/WProperties.h"
48 
49 #include "WGrid.h"
50 #include "WGridTransformOrtho.h"
51 
52 /**
53  * A grid that has parallelepiped cells which all have the same proportion. I.e.
54  * the samples along a single axis are equidistant. The distance of samples may
55  * vary between axes.
56  *
57  * \warning Positions on the upper bounddaries in x, y and z are considered outside the grid.
58  * \ingroup dataHandler
59  */
60 template< typename T >
61 class WGridRegular3DTemplate : public WGrid // NOLINT
62 {
63  // this (friend) is necessary to allow casting
64  template <class U>
65  friend class WGridRegular3DTemplate;
66  /**
67  * Only test are allowed as friends.
68  */
69  friend class WGridRegular3DTest;
70 public:
71  /**
72  * Convenience typedef for 3d vectors of the appropriate numerical type.
73  */
75 
76  /**
77  * Convenience typedef for a boost::shared_ptr< WGridRegular3DTemplate >.
78  */
79  typedef boost::shared_ptr< WGridRegular3DTemplate > SPtr;
80 
81  /**
82  * Convenience typedef for a boost::shared_ptr< const WGridRegular3DTemplate >.
83  */
84  typedef boost::shared_ptr< const WGridRegular3DTemplate > ConstSPtr;
85 
86  /**
87  * Convenience typedef for a boost::array< size_t, 8 >. Return type of getCellVertexIds.
88  */
89  typedef boost::array< size_t, 8 > CellVertexArray;
90 
91  /**
92  * Copy constructor.
93  * Copies the data from an WGridRegular3DTemplate object with arbitary numerical type.
94  *
95  * \param rhs A WGridRegular3DTemplate object, which mustn't have the same numerical type.
96  */
97  template< typename InputType >
98  WGridRegular3DTemplate( WGridRegular3DTemplate< InputType > const& rhs ); // NOLINT -- no explicit, this allows casts
99 
100  /**
101  * Defines the number of samples in each coordinate direction as ints,
102  * and the transformation of the grid via a grid transform.
103  *
104  * \param nbPosX number of positions along first axis
105  * \param nbPosY number of positions along second axis
106  * \param nbPosZ number of positions along third axis
107  * \param transform a grid transformation
108  */
109  WGridRegular3DTemplate( unsigned int nbPosX, unsigned int nbPosY, unsigned int nbPosZ,
111 
112  /**
113  * Returns the number of samples in x direction.
114  * \return The number of samples in x direction.
115  */
116  unsigned int getNbCoordsX() const;
117 
118  /**
119  * Returns the number of samples in y direction.
120  * \return The number of samples in y direction.
121  */
122  unsigned int getNbCoordsY() const;
123 
124  /**
125  * Returns the number of samples in z direction.
126  * \return The number of samples in z direction.
127  */
128  unsigned int getNbCoordsZ() const;
129 
130  /**
131  * Returns the distance between samples in x direction.
132  * \return The distance between samples in x direction.
133  */
134  T getOffsetX() const;
135 
136  /**
137  * Returns the distance between samples in y direction.
138  * \return The distance between samples in y direction.
139  */
140  T getOffsetY() const;
141 
142  /**
143  * Returns the distance between samples in z direction.
144  * \return The distance between samples in z direction.
145  */
146  T getOffsetZ() const;
147 
148  /**
149  * Returns the vector determining the direction of samples in x direction.
150  * Adding this vector to a grid position in world coordinates yields the position of the next sample
151  * along the grids (world coordinate) x-axis.
152  * \return The vector determining the direction of samples in x direction.
153  */
154  Vector3Type getDirectionX() const;
155 
156  /**
157  * Returns the vector determining the direction of samples in y direction.
158  * Adding this vector to a grid position in world coordinates yields the position of the next sample
159  * along the grids (world coordinate) y-axis.
160  * \return The vector determining the direction of samples in y direction.
161  */
162  Vector3Type getDirectionY() const;
163 
164  /**
165  * Returns the vector determining the direction of samples in z direction.
166  * Adding this vector to a grid position in world coordinates yields the position of the next sample
167  * along the grids (world coordinate) z-axis.
168  * \return The vector determining the direction of samples in z direction.
169  */
170  Vector3Type getDirectionZ() const;
171 
172  /**
173  * Returns the vector determining the unit (normalized) direction of samples in x direction.
174  * \return The vector determining the unit (normalized) direction of samples in x direction.
175  */
177 
178  /**
179  * Returns the vector determining the unit (normalized) direction of samples in y direction.
180  * \return The vector determining the unit (normalized) direction of samples in y direction.
181  */
183 
184  /**
185  * Returns the vector determining the unit (normalized) direction of samples in z direction.
186  * \return The vector determining the unit (normalized) direction of samples in z direction.
187  */
189 
190  /**
191  * Returns the position of the origin of the grid.
192  * \return The position of the origin of the grid.
193  */
194  Vector3Type getOrigin() const;
195 
196  /**
197  * Returns a 4x4 matrix that represents the grid's transformation.
198  * \return The grid's transformation.
199  */
201 
202  /**
203  * \copybrief WGrid::getBoundingBox()
204  * \return \copybrief WGrid::getBoundingBox()
205  */
207 
208  /**
209  * Returns the i-th position on the grid.
210  * \param i id of position to be obtained
211  * \return i-th position of the grid.
212  */
213  Vector3Type getPosition( unsigned int i ) const;
214 
215  /**
216  * Returns the position that is the iX-th in x direction, the iY-th in
217  * y direction and the iZ-th in z direction.
218  * \param iX id along first axis of position to be obtained
219  * \param iY id along second axis of position to be obtained
220  * \param iZ id along third axis of position to be obtained
221  * \return Position (iX,iY,iZ)
222  */
223  Vector3Type getPosition( unsigned int iX, unsigned int iY, unsigned int iZ ) const;
224 
225  /**
226  * Transforms world coordinates to texture coordinates.
227  * \param point The point with these coordinates will be transformed.
228  * \return point transformed into texture coordinate system
229  */
231 
232  /**
233  * Returns the i'th voxel where the given position belongs too.
234  *
235  * A voxel is a cuboid which surrounds a point on the grid.
236  *
237  * \verbatim
238  Voxel:
239  ______________ ____ (0.5, 0.5, 0.5)
240  /: /|
241  / : / |
242  / : / |
243  / : / |
244  _/____:_ ___ __/ |
245  | : | |
246  | : *<--|--------- grid point (0, 0, 0)
247  | :........|....|__
248  dz == 1| / | /
249  | / | / dy == 1
250  | / | /
251  _|/____________|/__
252  |<- dx == 1 ->|
253  -0.5,-0.5,-0.5
254  \endverbatim
255  *
256  * Please note the first voxel has only 1/8 of the size a normal voxel
257  * would have since all positions outside the grid do not belong
258  * to any voxel. Note: a cell is different to a voxel in terms of position.
259  * A voxel has a grid point as center whereas a cell has grid points as
260  * corners.
261  * \param pos Position for which we want to have the voxel number.
262  *
263  * \return Voxel number or -1 if the position refers to a point outside of
264  * the grid.
265  */
266  int getVoxelNum( const Vector3Type& pos ) const;
267 
268  /**
269  * returns the voxel index for a given discrete position in the grid
270  *
271  * \param x Position for which we want to have the voxel number.
272  * \param y Position for which we want to have the voxel number.
273  * \param z Position for which we want to have the voxel number.
274  *
275  * \return Voxel number or -1 if the position refers to a point outside of
276  * the grid.
277  */
278  int getVoxelNum( const size_t x, const size_t y, const size_t z ) const;
279  /**
280  * Computes the X coordinate of that voxel that contains the
281  * position pos.
282  *
283  * \param pos The position which selects the voxel for which the X
284  * coordinate is computed.
285  *
286  * \return The X coordinate or -1 if pos refers to point outside of the
287  * grid.
288  */
289  int getXVoxelCoord( const Vector3Type& pos ) const;
290 
291  /**
292  * Computes the Y coordinate of that voxel that contains the
293  * position pos.
294  *
295  * \param pos The position which selects the voxel for which the Y
296  * coordinate is computed.
297  *
298  * \return The Y coordinate or -1 if pos refers to point outside of the
299  * grid.
300  */
301  int getYVoxelCoord( const Vector3Type& pos ) const;
302 
303  /**
304  * Computes the Z coordinate of that voxel that contains the
305  * position pos.
306  *
307  * \param pos The position which selects the voxel for which the Z
308  * coordinate is computed.
309  *
310  * \return The Z coordinate or -1 if pos refers to point outside of the
311  * grid.
312  */
313  int getZVoxelCoord( const Vector3Type& pos ) const;
314 
315  /**
316  * Computes the voxel coordinates of that voxel which contains
317  * the position pos.
318  *
319  * \param pos The position selecting the voxel.
320  *
321  * \return A vector of ints where the first component is the X voxel
322  * coordinate, the second the Y component voxel coordinate and the last the
323  * Z component of the voxel coordinate. If the selecting position is
324  * outside of the grid then -1 -1 -1 is returned.
325  */
326  WVector3i getVoxelCoord( const Vector3Type& pos ) const;
327 
328  /**
329  * Computes the id of the cell containing the position pos. Note that the upper
330  * bound of the grid does not belong to any cell
331  *
332  * \param pos The position selecting the cell.
333  * \param success True if the position pos is inside the grid.
334  *
335  * \return id of the containing the position.
336  */
337  size_t getCellId( const Vector3Type& pos, bool* success ) const;
338 
339  /**
340  * Computes the ids of the vertices of a cell given by its id.
341  *
342  * \param cellId The id of the cell we want to know ther vertices of.
343  *
344  * \return Ids of vertices belonging to cell with given cellId.
345 
346  * \verbatim
347  z-axis y-axis
348  | /
349  | 6___/_7
350  |/: /|
351  4_:___5 |
352  | :...|.|
353  |.2 | 3
354  |_____|/ ____x-axis
355  0 1
356  \endverbatim
357  *
358  */
359  CellVertexArray getCellVertexIds( size_t cellId ) const;
360 
361  /**
362  * Computes the vertices for a voxel cuboid around the given point:
363  *
364  * \verbatim
365  z-axis y-axis
366  | /
367  | h___/_g
368  |/: /|
369  d_:___c |
370  | :...|.|
371  |.e | f
372  |_____|/ ____x-axis
373  a b
374  \endverbatim
375  *
376  * As you can see the order of the points is: a, b, c, d, e, f, g, h.
377  *
378  * \param point Center of the cuboid which must not necesarrily be a point
379  * of the grid.
380  * \param margin If you need to shrink the Voxel put here the delta > 0.
381  *
382  * \return Reference to a list of vertices which are the corner points of
383  * the cube. Note this must not be a voxel, but has the same size of the an
384  * voxel. If you need voxels at grid positions fill this function with
385  * voxel center positions aka grid points.
386  */
387  boost::shared_ptr< std::vector< Vector3Type > > getVoxelVertices( const Vector3Type& point,
388  const T margin = 0.0 ) const;
389 
390  /**
391  * Return the list of neighbour voxels.
392  *
393  * \throw WOutOfBounds If the voxel id is outside of the grid.
394  *
395  * \param id Number of the voxel for which the neighbours should be computed
396  *
397  * \return Vector of voxel ids which are all neighboured
398  */
399  std::vector< size_t > getNeighbours( size_t id ) const;
400 
401  /**
402  * Return the list of all neighbour voxels.
403  *
404  * \throw WOutOfBounds If the voxel id is outside of the grid.
405  *
406  * \param id Number of the voxel for which the neighbours should be computed
407  *
408  * \return Vector of voxel ids which are all neighboured
409  */
410  std::vector< size_t > getNeighbours27( size_t id ) const;
411 
412  /**
413  * Return the list of all neighbour voxels.
414  *
415  * \throw WOutOfBounds If the voxel id is outside of the grid.
416  *
417  * \param id Number of the voxel for which the neighbours should be computed
418  *
419  * \param range neighborhood range selected. It specifies the distance to count as neighbour in each direction.
420  *
421  * \return Vector of voxel ids which are all neighboured
422  */
423  std::vector< size_t > getNeighboursRange( size_t id, size_t range ) const;
424 
425  /**
426  * Return the list of all neighbour voxels.
427  *
428  * \throw WOutOfBounds If the voxel id is outside of the grid.
429  *
430  * \param id Number of the voxel for which the neighbours should be computed
431  *
432  * \return Vector of voxel ids which are all neighboured along the XY plane
433  */
434  std::vector< size_t > getNeighbours9XY( size_t id ) const;
435 
436  /**
437  * Return the list of all neighbour voxels.
438  *
439  * \throw WOutOfBounds If the voxel id is outside of the grid.
440  *
441  * \param id Number of the voxel for which the neighbours should be computed
442  *
443  * \return Vector of voxel ids which are all neighboured along the YZ plane
444  */
445  std::vector< size_t > getNeighbours9YZ( size_t id ) const;
446 
447  /**
448  * Return the list of all neighbour voxels.
449  *
450  * \throw WOutOfBounds If the voxel id is outside of the grid.
451  *
452  * \param id Number of the voxel for which the neighbours should be computed
453  *
454  * \return Vector of voxel ids which are all neighboured along the XZ plane
455  */
456  std::vector< size_t > getNeighbours9XZ( size_t id ) const;
457 
458  /**
459  * Decides whether a certain position is inside this grid or not.
460  *
461  * \param pos Position to test
462  *
463  * \return True if and only if the given point is inside or on boundary of this grid, otherwise false.
464  */
465  bool encloses( const Vector3Type& pos ) const;
466 
467  /**
468  * Return whether the transformations of the grid are only translation and/or scaling
469  * \return Transformation does not contain rotation?
470  */
471  bool isNotRotated() const;
472 
473  /**
474  * Returns the transformation used by this grid.
475  * \return The transformation.
476  */
478 
479 protected:
480 private:
481  /**
482  * Computes for the n'th component of the voxel coordinate where the voxel
483  * contains the position pos.
484  *
485  * \param pos The position for which the n'th component of the voxel
486  * coordinates should be computed.
487  * \param axis The number of the component. (0 == x-axis, 1 == y-axis, ...)
488  *
489  * \return The n'th component of the voxel coordinate
490  */
491  int getNVoxelCoord( const Vector3Type& pos, size_t axis ) const;
492 
493  /**
494  * Adds the specific information of this grid type to the
495  * informational properties.
496  */
498 
499  unsigned int m_nbPosX; //!< Number of positions in x direction
500  unsigned int m_nbPosY; //!< Number of positions in y direction
501  unsigned int m_nbPosZ; //!< Number of positions in z direction
502 
503  //! The grid's transformation.
505 };
506 
507 // Convenience typedefs
511 
512 template< typename T >
513 template< typename InputType >
515  WGrid( rhs.m_nbPosX * rhs.m_nbPosY * rhs.m_nbPosZ ),
516  m_nbPosX( rhs.m_nbPosX ),
517  m_nbPosY( rhs.m_nbPosY ),
518  m_nbPosZ( rhs.m_nbPosZ ),
519  m_transform( rhs.m_transform )
520 {
522 }
523 
524 template< typename T >
525 WGridRegular3DTemplate< T >::WGridRegular3DTemplate( unsigned int nbPosX, unsigned int nbPosY, unsigned int nbPosZ,
526  WGridTransformOrthoTemplate< T > const transform )
527  : WGrid( nbPosX * nbPosY * nbPosZ ),
528  m_nbPosX( nbPosX ),
529  m_nbPosY( nbPosY ),
530  m_nbPosZ( nbPosZ ),
531  m_transform( transform )
532 {
534 }
535 
536 template< typename T >
537 inline unsigned int WGridRegular3DTemplate< T >::getNbCoordsX() const
538 {
539  return m_nbPosX;
540 }
541 
542 template< typename T >
543 inline unsigned int WGridRegular3DTemplate< T >::getNbCoordsY() const
544 {
545  return m_nbPosY;
546 }
547 
548 template< typename T >
549 inline unsigned int WGridRegular3DTemplate< T >::getNbCoordsZ() const
550 {
551  return m_nbPosZ;
552 }
553 
554 template< typename T >
556 {
557  return m_transform.getOffsetX();
558 }
559 
560 template< typename T >
562 {
563  return m_transform.getOffsetY();
564 }
565 
566 template< typename T >
568 {
569  return m_transform.getOffsetZ();
570 }
571 
572 template< typename T >
574 {
575  return m_transform.getDirectionX();
576 }
577 
578 template< typename T >
580 {
581  return m_transform.getDirectionY();
582 }
583 
584 template< typename T >
586 {
587  return m_transform.getDirectionZ();
588 }
589 
590 template< typename T >
592 {
593  return m_transform.getUnitDirectionX();
594 }
595 
596 template< typename T >
598 {
599  return m_transform.getUnitDirectionY();
600 }
601 
602 template< typename T >
604 {
605  return m_transform.getUnitDirectionZ();
606 }
607 
608 template< typename T >
610 {
611  return m_transform.getOrigin();
612 }
613 
614 template< typename T >
616 {
617  return m_transform.getTransformationMatrix();
618 }
619 
620 template< typename T >
622 {
623  WBoundingBox result;
624  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( 0.0, 0.0, 0.0 ) ) );
625  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( getNbCoordsX() - 1, 0.0, 0.0 ) ) );
626  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( 0.0, getNbCoordsY() - 1, 0.0 ) ) );
627  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( getNbCoordsX() - 1, getNbCoordsY() - 1, 0.0 ) ) );
628  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( 0.0, 0.0, getNbCoordsZ() - 1 ) ) );
629  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( getNbCoordsX() - 1, 0.0, getNbCoordsZ() - 1 ) ) );
630  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( 0.0, getNbCoordsY() - 1, getNbCoordsZ() - 1 ) ) );
631  result.expandBy( m_transform.positionToWorldSpace( Vector3Type( getNbCoordsX() - 1, getNbCoordsY() - 1, getNbCoordsZ() - 1 ) ) );
632  return result;
633 }
634 
635 template< typename T >
637 {
638  return getPosition( i % m_nbPosX, ( i / m_nbPosX ) % m_nbPosY, i / ( m_nbPosX * m_nbPosY ) );
639 }
640 
641 template< typename T >
643  unsigned int iY,
644  unsigned int iZ ) const
645 {
646  Vector3Type i( iX, iY, iZ );
647  return m_transform.positionToWorldSpace( i );
648 }
649 
650 template< typename T >
652 {
653  Vector3Type r( m_transform.positionToGridSpace( point ) );
654 
655  // Scale to [0,1]
656  r[0] = r[0] / m_nbPosX;
657  r[1] = r[1] / m_nbPosY;
658  r[2] = r[2] / m_nbPosZ;
659 
660  // Correct the coordinates to have the position at the center of the texture voxel.
661  r[0] += 0.5 / m_nbPosX;
662  r[1] += 0.5 / m_nbPosY;
663  r[2] += 0.5 / m_nbPosZ;
664 
665  return r;
666 }
667 
668 template< typename T >
670 {
671  // Note: the reason for the +1 is that the first and last Voxel in a x-axis
672  // row are cut.
673  //
674  // y-axis
675  // _|_______ ___ this is the 3rd Voxel
676  // 1 | | | v
677  // |...............
678  // _|_:_|_:_|_:_|_:____ x-axis
679  // | : | : | : | :
680  // |.:...:...:...:.
681  // 0 1 2
682  int xVoxelCoord = getXVoxelCoord( pos );
683  int yVoxelCoord = getYVoxelCoord( pos );
684  int zVoxelCoord = getZVoxelCoord( pos );
685  if( xVoxelCoord == -1 || yVoxelCoord == -1 || zVoxelCoord == -1 )
686  {
687  return -1;
688  }
689  return xVoxelCoord
690  + yVoxelCoord * ( m_nbPosX )
691  + zVoxelCoord * ( m_nbPosX ) * ( m_nbPosY );
692 }
693 
694 template< typename T >
695 inline int WGridRegular3DTemplate< T >::getVoxelNum( const size_t x, const size_t y, const size_t z ) const
696 {
697  // since we use size_t here only a check for the upper bounds is needed
698  if( x > m_nbPosX || y > m_nbPosY || z > m_nbPosZ )
699  {
700  return -1;
701  }
702  return x + y * ( m_nbPosX ) + z * ( m_nbPosX ) * ( m_nbPosY );
703 }
704 
705 template< typename T >
707 {
708  // the current get*Voxel stuff is too complicated anyway
709  Vector3Type v = m_transform.positionToGridSpace( pos );
710 
711  // this part could be refactored into an inline function
712  T d;
713  v[ 2 ] = std::modf( v[ 0 ] + T( 0.5 ), &d );
714  int i = static_cast< int >( v[ 0 ] >= T( 0.0 ) && v[ 0 ] < m_nbPosX - T( 1.0 ) );
715  return -1 + i * static_cast< int >( T( 1.0 ) + d );
716 }
717 
718 template< typename T >
720 {
721  Vector3Type v = m_transform.positionToGridSpace( pos );
722 
723  T d;
724  v[ 0 ] = std::modf( v[ 1 ] + T( 0.5 ), &d );
725  int i = static_cast< int >( v[ 1 ] >= T( 0.0 ) && v[ 1 ] < m_nbPosY - T( 1.0 ) );
726  return -1 + i * static_cast< int >( T( 1.0 ) + d );
727 }
728 
729 template< typename T >
731 {
732  Vector3Type v = m_transform.positionToGridSpace( pos );
733 
734  T d;
735  v[ 0 ] = std::modf( v[ 2 ] + T( 0.5 ), &d );
736  int i = static_cast< int >( v[ 2 ] >= T( 0.0 ) && v[ 2 ] < m_nbPosZ - T( 1.0 ) );
737  return -1 + i * static_cast< int >( T( 1.0 ) + d );
738 }
739 
740 template< typename T >
742 {
743  WVector3i result;
744  result[0] = getXVoxelCoord( pos );
745  result[1] = getYVoxelCoord( pos );
746  result[2] = getZVoxelCoord( pos );
747  return result;
748 }
749 
750 template< typename T >
752 {
753  Vector3Type v = m_transform.positionToGridSpace( pos );
754 
755  T xCellId = floor( v[0] );
756  T yCellId = floor( v[1] );
757  T zCellId = floor( v[2] );
758 
759  *success = xCellId >= 0 && yCellId >=0 && zCellId >= 0 && xCellId < m_nbPosX - 1 && yCellId < m_nbPosY -1 && zCellId < m_nbPosZ -1;
760 
761  return xCellId + yCellId * ( m_nbPosX - 1 ) + zCellId * ( m_nbPosX - 1 ) * ( m_nbPosY - 1 );
762 }
763 
764 template< typename T >
766 {
768  size_t minVertexIdZ = cellId / ( ( m_nbPosX - 1 ) * ( m_nbPosY - 1 ) );
769  size_t remainderXY = cellId - minVertexIdZ * ( ( m_nbPosX - 1 ) * ( m_nbPosY - 1 ) );
770  size_t minVertexIdY = remainderXY / ( m_nbPosX - 1 );
771  size_t minVertexIdX = remainderXY % ( m_nbPosX - 1 );
772 
773  size_t minVertexId = minVertexIdX + minVertexIdY * m_nbPosX + minVertexIdZ * m_nbPosX * m_nbPosY;
774 
775  vertices[0] = minVertexId;
776  vertices[1] = vertices[0] + 1;
777  vertices[2] = minVertexId + m_nbPosX;
778  vertices[3] = vertices[2] + 1;
779  vertices[4] = minVertexId + m_nbPosX * m_nbPosY;
780  vertices[5] = vertices[4] + 1;
781  vertices[6] = vertices[4] + m_nbPosX;
782  vertices[7] = vertices[6] + 1;
783  return vertices;
784 }
785 
786 template< typename T >
787 boost::shared_ptr< std::vector< typename WGridRegular3DTemplate< T >::Vector3Type > > WGridRegular3DTemplate< T >::getVoxelVertices( const WGridRegular3DTemplate< T >::Vector3Type& point, const T margin ) const // NOLINT -- too long line
788 {
789  typedef boost::shared_ptr< std::vector< Vector3Type > > ReturnType;
790  ReturnType result = ReturnType( new std::vector< Vector3Type > );
791  result->reserve( 8 );
792  T halfMarginX = getOffsetX() / 2.0 - std::abs( margin );
793  T halfMarginY = getOffsetY() / 2.0 - std::abs( margin );
794  T halfMarginZ = getOffsetZ() / 2.0 - std::abs( margin );
795  result->push_back( Vector3Type( point[0] - halfMarginX, point[1] - halfMarginY, point[2] - halfMarginZ ) ); // a
796  result->push_back( Vector3Type( point[0] + halfMarginX, point[1] - halfMarginY, point[2] - halfMarginZ ) ); // b
797  result->push_back( Vector3Type( point[0] + halfMarginX, point[1] - halfMarginY, point[2] + halfMarginZ ) ); // c
798  result->push_back( Vector3Type( point[0] - halfMarginX, point[1] - halfMarginY, point[2] + halfMarginZ ) ); // d
799  result->push_back( Vector3Type( point[0] - halfMarginX, point[1] + halfMarginY, point[2] - halfMarginZ ) ); // e
800  result->push_back( Vector3Type( point[0] + halfMarginX, point[1] + halfMarginY, point[2] - halfMarginZ ) ); // f
801  result->push_back( Vector3Type( point[0] + halfMarginX, point[1] + halfMarginY, point[2] + halfMarginZ ) ); // g
802  result->push_back( Vector3Type( point[0] - halfMarginX, point[1] + halfMarginY, point[2] + halfMarginZ ) ); // h
803  return result;
804 }
805 
806 template< typename T >
807 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighbours( size_t id ) const
808 {
809  std::vector< size_t > neighbours;
810  size_t x = id % m_nbPosX;
811  size_t y = ( id / m_nbPosX ) % m_nbPosY;
812  size_t z = id / ( m_nbPosX * m_nbPosY );
813 
814  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
815  {
816  std::stringstream ss;
817  ss << "This point: " << id << " is not part of this grid: ";
818  ss << " nbPosX: " << m_nbPosX;
819  ss << " nbPosY: " << m_nbPosY;
820  ss << " nbPosZ: " << m_nbPosZ;
821  throw WOutOfBounds( ss.str() );
822  }
823  // for every neighbour we must check if its not on the boundary, it will be skipped otherwise
824  if( x > 0 )
825  {
826  neighbours.push_back( id - 1 );
827  }
828  if( x < m_nbPosX - 1 )
829  {
830  neighbours.push_back( id + 1 );
831  }
832  if( y > 0 )
833  {
834  neighbours.push_back( id - m_nbPosX );
835  }
836  if( y < m_nbPosY - 1 )
837  {
838  neighbours.push_back( id + m_nbPosX );
839  }
840  if( z > 0 )
841  {
842  neighbours.push_back( id - ( m_nbPosX * m_nbPosY ) );
843  }
844  if( z < m_nbPosZ - 1 )
845  {
846  neighbours.push_back( id + ( m_nbPosX * m_nbPosY ) );
847  }
848  return neighbours;
849 }
850 
851 template< typename T >
852 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighbours27( size_t id ) const
853 {
854  std::vector< size_t > neighbours;
855  size_t x = id % m_nbPosX;
856  size_t y = ( id / m_nbPosX ) % m_nbPosY;
857  size_t z = id / ( m_nbPosX * m_nbPosY );
858 
859  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
860  {
861  std::stringstream ss;
862  ss << "This point: " << id << " is not part of this grid: ";
863  ss << " nbPosX: " << m_nbPosX;
864  ss << " nbPosY: " << m_nbPosY;
865  ss << " nbPosZ: " << m_nbPosZ;
866  throw WOutOfBounds( ss.str() );
867  }
868  // for every neighbour we must check if its not on the boundary, it will be skipped otherwise
869  std::vector< int >tempResult;
870 
871  tempResult.push_back( getVoxelNum( x , y , z ) );
872  tempResult.push_back( getVoxelNum( x , y - 1, z ) );
873  tempResult.push_back( getVoxelNum( x , y + 1, z ) );
874  tempResult.push_back( getVoxelNum( x - 1, y , z ) );
875  tempResult.push_back( getVoxelNum( x - 1, y - 1, z ) );
876  tempResult.push_back( getVoxelNum( x - 1, y + 1, z ) );
877  tempResult.push_back( getVoxelNum( x + 1, y , z ) );
878  tempResult.push_back( getVoxelNum( x + 1, y - 1, z ) );
879  tempResult.push_back( getVoxelNum( x + 1, y + 1, z ) );
880 
881  tempResult.push_back( getVoxelNum( x , y , z - 1 ) );
882  tempResult.push_back( getVoxelNum( x , y - 1, z - 1 ) );
883  tempResult.push_back( getVoxelNum( x , y + 1, z - 1 ) );
884  tempResult.push_back( getVoxelNum( x - 1, y , z - 1 ) );
885  tempResult.push_back( getVoxelNum( x - 1, y - 1, z - 1 ) );
886  tempResult.push_back( getVoxelNum( x - 1, y + 1, z - 1 ) );
887  tempResult.push_back( getVoxelNum( x + 1, y , z - 1 ) );
888  tempResult.push_back( getVoxelNum( x + 1, y - 1, z - 1 ) );
889  tempResult.push_back( getVoxelNum( x + 1, y + 1, z - 1 ) );
890 
891  tempResult.push_back( getVoxelNum( x , y , z + 1 ) );
892  tempResult.push_back( getVoxelNum( x , y - 1, z + 1 ) );
893  tempResult.push_back( getVoxelNum( x , y + 1, z + 1 ) );
894  tempResult.push_back( getVoxelNum( x - 1, y , z + 1 ) );
895  tempResult.push_back( getVoxelNum( x - 1, y - 1, z + 1 ) );
896  tempResult.push_back( getVoxelNum( x - 1, y + 1, z + 1 ) );
897  tempResult.push_back( getVoxelNum( x + 1, y , z + 1 ) );
898  tempResult.push_back( getVoxelNum( x + 1, y - 1, z + 1 ) );
899  tempResult.push_back( getVoxelNum( x + 1, y + 1, z + 1 ) );
900 
901  for( size_t k = 0; k < tempResult.size(); ++k )
902  {
903  if( tempResult[k] != -1 )
904  {
905  neighbours.push_back( tempResult[k] );
906  }
907  }
908  return neighbours;
909 }
910 
911 template< typename T >
912 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighboursRange( size_t id, size_t range ) const
913 {
914  std::vector< size_t > neighbours;
915  size_t x = id % m_nbPosX;
916  size_t y = ( id / m_nbPosX ) % m_nbPosY;
917  size_t z = id / ( m_nbPosX * m_nbPosY );
918 
919  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
920  {
921  std::stringstream ss;
922  ss << "This point: " << id << " is not part of this grid: ";
923  ss << " nbPosX: " << m_nbPosX;
924  ss << " nbPosY: " << m_nbPosY;
925  ss << " nbPosZ: " << m_nbPosZ;
926  throw WOutOfBounds( ss.str() );
927  }
928  // for every neighbour we must check if its not on the boundary, it will be skipped otherwise
929  std::vector< int >tempResult;
930 
931  for( size_t xx = x - range; xx < x + range + 1; ++xx )
932  {
933  for( size_t yy = y - range; yy < y + range + 1; ++yy )
934  {
935  for( size_t zz = z - range; zz < z + range + 1; ++zz )
936  {
937  tempResult.push_back( getVoxelNum( xx, yy, zz ) );
938  }
939  }
940  }
941 
942  for( size_t k = 0; k < tempResult.size(); ++k )
943  {
944  if( tempResult[k] != -1 )
945  {
946  neighbours.push_back( tempResult[k] );
947  }
948  }
949  return neighbours;
950 }
951 
952 template< typename T >
953 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighbours9XY( size_t id ) const
954 {
955  std::vector< size_t > neighbours;
956  size_t x = id % m_nbPosX;
957  size_t y = ( id / m_nbPosX ) % m_nbPosY;
958  size_t z = id / ( m_nbPosX * m_nbPosY );
959 
960  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
961  {
962  std::stringstream ss;
963  ss << "This point: " << id << " is not part of this grid: ";
964  ss << " nbPosX: " << m_nbPosX;
965  ss << " nbPosY: " << m_nbPosY;
966  ss << " nbPosZ: " << m_nbPosZ;
967  throw WOutOfBounds( ss.str() );
968  }
969  // boundary check now happens in the getVoxelNum function
970  int vNum;
971 
972  vNum = getVoxelNum( x - 1, y, z );
973  if( vNum != -1 )
974  {
975  neighbours.push_back( vNum );
976  }
977  vNum = getVoxelNum( x - 1, y - 1, z );
978  if( vNum != -1 )
979  {
980  neighbours.push_back( vNum );
981  }
982  vNum = getVoxelNum( x, y - 1, z );
983  if( vNum != -1 )
984  {
985  neighbours.push_back( vNum );
986  }
987  vNum = getVoxelNum( x + 1, y - 1, z );
988  if( vNum != -1 )
989  {
990  neighbours.push_back( vNum );
991  }
992  vNum = getVoxelNum( x + 1, y, z );
993  if( vNum != -1 )
994  {
995  neighbours.push_back( vNum );
996  }
997  vNum = getVoxelNum( x + 1, y + 1, z );
998  if( vNum != -1 )
999  {
1000  neighbours.push_back( vNum );
1001  }
1002  vNum = getVoxelNum( x, y + 1, z );
1003  if( vNum != -1 )
1004  {
1005  neighbours.push_back( vNum );
1006  }
1007  vNum = getVoxelNum( x - 1, y + 1, z );
1008  if( vNum != -1 )
1009  {
1010  neighbours.push_back( vNum );
1011  }
1012  return neighbours;
1013 }
1014 
1015 template< typename T >
1016 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighbours9YZ( size_t id ) const
1017 {
1018  std::vector< size_t > neighbours;
1019  size_t x = id % m_nbPosX;
1020  size_t y = ( id / m_nbPosX ) % m_nbPosY;
1021  size_t z = id / ( m_nbPosX * m_nbPosY );
1022 
1023  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
1024  {
1025  std::stringstream ss;
1026  ss << "This point: " << id << " is not part of this grid: ";
1027  ss << " nbPosX: " << m_nbPosX;
1028  ss << " nbPosY: " << m_nbPosY;
1029  ss << " nbPosZ: " << m_nbPosZ;
1030  throw WOutOfBounds( ss.str() );
1031  }
1032  // boundary check now happens in the getVoxelNum function
1033  int vNum;
1034 
1035  vNum = getVoxelNum( x, y, z - 1 );
1036  if( vNum != -1 )
1037  {
1038  neighbours.push_back( vNum );
1039  }
1040  vNum = getVoxelNum( x, y - 1, z - 1 );
1041  if( vNum != -1 )
1042  {
1043  neighbours.push_back( vNum );
1044  }
1045  vNum = getVoxelNum( x, y - 1, z );
1046  if( vNum != -1 )
1047  {
1048  neighbours.push_back( vNum );
1049  }
1050  vNum = getVoxelNum( x, y - 1, z + 1 );
1051  if( vNum != -1 )
1052  {
1053  neighbours.push_back( vNum );
1054  }
1055  vNum = getVoxelNum( x, y, z + 1 );
1056  if( vNum != -1 )
1057  {
1058  neighbours.push_back( vNum );
1059  }
1060  vNum = getVoxelNum( x, y + 1, z + 1 );
1061  if( vNum != -1 )
1062  {
1063  neighbours.push_back( vNum );
1064  }
1065  vNum = getVoxelNum( x, y + 1, z );
1066  if( vNum != -1 )
1067  {
1068  neighbours.push_back( vNum );
1069  }
1070  vNum = getVoxelNum( x, y + 1, z - 1 );
1071  if( vNum != -1 )
1072  {
1073  neighbours.push_back( vNum );
1074  }
1075 
1076  return neighbours;
1077 }
1078 
1079 template< typename T >
1080 std::vector< size_t > WGridRegular3DTemplate< T >::getNeighbours9XZ( size_t id ) const
1081 {
1082  std::vector< size_t > neighbours;
1083  size_t x = id % m_nbPosX;
1084  size_t y = ( id / m_nbPosX ) % m_nbPosY;
1085  size_t z = id / ( m_nbPosX * m_nbPosY );
1086 
1087  if( x >= m_nbPosX || y >= m_nbPosY || z >= m_nbPosZ )
1088  {
1089  std::stringstream ss;
1090  ss << "This point: " << id << " is not part of this grid: ";
1091  ss << " nbPosX: " << m_nbPosX;
1092  ss << " nbPosY: " << m_nbPosY;
1093  ss << " nbPosZ: " << m_nbPosZ;
1094  throw WOutOfBounds( ss.str() );
1095  }
1096  // boundary check now happens in the getVoxelNum function
1097  int vNum;
1098 
1099  vNum = getVoxelNum( x, y, z - 1 );
1100  if( vNum != -1 )
1101  {
1102  neighbours.push_back( vNum );
1103  }
1104  vNum = getVoxelNum( x - 1, y, z - 1 );
1105  if( vNum != -1 )
1106  {
1107  neighbours.push_back( vNum );
1108  }
1109  vNum = getVoxelNum( x - 1, y, z );
1110  if( vNum != -1 )
1111  {
1112  neighbours.push_back( vNum );
1113  }
1114  vNum = getVoxelNum( x - 1, y, z + 1 );
1115  if( vNum != -1 )
1116  {
1117  neighbours.push_back( vNum );
1118  }
1119  vNum = getVoxelNum( x, y, z + 1 );
1120  if( vNum != -1 )
1121  {
1122  neighbours.push_back( vNum );
1123  }
1124  vNum = getVoxelNum( x + 1, y, z + 1 );
1125  if( vNum != -1 )
1126  {
1127  neighbours.push_back( vNum );
1128  }
1129  vNum = getVoxelNum( x + 1, y, z );
1130  if( vNum != -1 )
1131  {
1132  neighbours.push_back( vNum );
1133  }
1134  vNum = getVoxelNum( x + 1, y, z - 1 );
1135  if( vNum != -1 )
1136  {
1137  neighbours.push_back( vNum );
1138  }
1139 
1140  return neighbours;
1141 }
1142 
1143 template< typename T >
1145 {
1146  Vector3Type v = m_transform.positionToGridSpace( pos );
1147 
1148  if( v[ 0 ] < T( 0.0 ) || v[ 0 ] >= static_cast< T >( m_nbPosX - 1 ) )
1149  {
1150  return false;
1151  }
1152  if( v[ 1 ] < T( 0.0 ) || v[ 1 ] >= static_cast< T >( m_nbPosY - 1 ) )
1153  {
1154  return false;
1155  }
1156  if( v[ 2 ] < T( 0.0 ) || v[ 2 ] >= static_cast< T >( m_nbPosZ - 1 ) )
1157  {
1158  return false;
1159  }
1160  return true;
1161 }
1162 
1163 template< typename T >
1165 {
1166  return m_transform.isNotRotated();
1167 }
1168 
1169 template< typename T >
1171 {
1172  return m_transform;
1173 }
1174 
1175 template< typename T >
1177 {
1178  WPropInt xDim = m_infoProperties->addProperty( "X dimension: ",
1179  "The x dimension of the grid.",
1180  static_cast<int>( getNbCoordsX() ) );
1181  WPropInt yDim = m_infoProperties->addProperty( "Y dimension: ",
1182  "The y dimension of the grid.",
1183  static_cast<int>( getNbCoordsY() ) );
1184  WPropInt zDim = m_infoProperties->addProperty( "Z dimension: ",
1185  "The z dimension of the grid.",
1186  static_cast<int>( getNbCoordsZ() ) );
1187 
1188  WPropDouble xOffset = m_infoProperties->addProperty( "X offset: ",
1189  "The distance between samples in x direction",
1190  static_cast< double >( getOffsetX() ) );
1191  WPropDouble yOffset = m_infoProperties->addProperty( "Y offset: ",
1192  "The distance between samples in y direction",
1193  static_cast< double >( getOffsetY() ) );
1194  WPropDouble zOffset = m_infoProperties->addProperty( "Z offset: ",
1195  "The distance between samples in z direction",
1196  static_cast< double >( getOffsetZ() ) );
1197 
1198  WPropMatrix4X4 transformation = m_infoProperties->addProperty( "Transformation",
1199  "The transformation of this grid.",
1200  static_cast< WMatrix4d >( getTransform() ) );
1201 }
1202 
1203 // +----------------------+
1204 // | non-member functions |
1205 // +----------------------+
1206 
1207 /**
1208  * Convinience function returning all offsets per axis.
1209  * 0 : xAxis, 1 : yAxis, 2 : zAxis
1210  * \param grid The grid having the information.
1211  * \note Implementing this as NonMemberNonFriend was intentional.
1212  * \return Array of number of samples per axis.
1213  */
1214 template< typename T >
1215 inline boost::array< T, 3 > getOffsets( boost::shared_ptr< const WGridRegular3DTemplate< T > > grid )
1216 {
1217  boost::array< T, 3 > result = { { grid->getOffsetX(), grid->getOffsetY(), grid->getOffsetZ() } }; // NOLINT curly braces
1218  return result;
1219 }
1220 
1221 /**
1222  * Convinience function returning all number coords per axis.
1223  * 0 : xAxis, 1 : yAxis, 2 : zAxis
1224  * \param grid The grid having the information.
1225  * \note Implementing this as NonMemberNonFriend was intentional.
1226  * \return Array of number of samples per axis.
1227  */
1228 template< typename T >
1229 inline boost::array< unsigned int, 3 > getNbCoords( boost::shared_ptr< const WGridRegular3DTemplate< T > > grid )
1230 {
1231  boost::array< unsigned int, 3 > result = { { grid->getNbCoordsX(), grid->getNbCoordsY(), grid->getNbCoordsZ() } }; // NOLINT curly braces
1232  return result;
1233 }
1234 
1235 /**
1236  * Convinience function returning all axis directions.
1237  * 0 : xAxis, 1 : yAxis, 2 : zAxis
1238  * \param grid The grid having the information.
1239  * \note Implementing this as NonMemberNonFriend was intentional.
1240  * \return The direction of each axis as array
1241  */
1242 template< typename T >
1243 inline boost::array< typename WGridRegular3DTemplate< T >::Vector3Type, 3 > getDirections( boost::shared_ptr< const WGridRegular3DTemplate< T > > grid ) // NOLINT -- too long line
1244 {
1245  boost::array< typename WGridRegular3DTemplate< T >::Vector3Type, 3 > result = { { grid->getDirectionX(), grid->getDirectionY(), grid->getDirectionZ() } }; // NOLINT curly braces
1246  return result;
1247 }
1248 
1249 /**
1250  * Convinience function returning all axis unit directions.
1251  * 0 : xAxis, 1 : yAxis, 2 : zAxis
1252  * \param grid The grid having the information.
1253  * \note Implementing this as NonMemberNonFriend was intentional.
1254  * \return The direction of each axis as array
1255  */
1256 template< typename T >
1257 inline boost::array< typename WGridRegular3DTemplate< T >::Vector3Type, 3 > getUnitDirections( boost::shared_ptr< const WGridRegular3DTemplate< T > > grid ) // NOLINT -- too long line
1258 {
1259  boost::array< typename WGridRegular3DTemplate< T >::Vector3Type, 3 > result = { { grid->getUnitDirectionX(), grid->getUnitDirectionY(), grid->getUnitDirectionZ() } }; // NOLINT curly braces
1260  return result;
1261 }
1262 
1263 #endif // WGRIDREGULAR3D_H