Re: [eigen] Enhanced AlignedBox

[Manuel Yguel <manuel.yguel@xxxxxxxxx>]

Hi Gael,
thanks for the useful remarks.
I agree with all you said and I modify the code according to what you suggested.

I am not at all familiar with how the expressions work in Eigen, so
consider it a try.
It seems to me that there is probably a simplest way to write the
code, so correct me if there is one: I will learn.

I am thinking of the same kind of design that the one you suggested
for the random number generator.
I like the interface of the boost random number generators, so I will
try to make a proposition compatible with those ones.


- best regards,

Manuel

# HG changeset patch
# User Manuel Yguel <manuel.yguel@xxxxxxxxx>
# Date 1261419866 -3600
# Node ID 2d0813a2fb84155fa231076c992ae6370914f7f4
# Parent  8b3294611c0b3660a5e71293153aa2364df8c46f
Add convenience functions: sizes that get the lengths of the different sides of the bounding box, diagonal, volume that get the diagonal vector and the volume of the bounding box and sample that throw a point from the bounding box uniformly.

diff -r 8b3294611c0b -r 2d0813a2fb84 Eigen/src/Geometry/AlignedBox.h
--- a/Eigen/src/Geometry/AlignedBox.h	Mon Dec 21 11:35:08 2009 +0000
+++ b/Eigen/src/Geometry/AlignedBox.h	Mon Dec 21 19:24:26 2009 +0100
@@ -43,9 +43,28 @@
 public:
 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
   enum { AmbientDimAtCompileTime = _AmbientDim };
-  typedef _Scalar Scalar;
-  typedef typename NumTraits<Scalar>::Real RealScalar;
-  typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType;
+  typedef _Scalar                                   Scalar;
+  typedef typename NumTraits<Scalar>::Real          RealScalar;
+  typedef typename NumTraits<Scalar>::FloatingPoint FloatingPoint;
+  typedef Matrix<Scalar,AmbientDimAtCompileTime,1>  VectorType;
+
+  /** Define constants to name the corners of a 1D, 2D or 3D axis aligned bounding box */
+  typedef enum
+  {
+	/** 1D names */
+	Min=0, Max=1,
+
+	/** Added names for 2D */
+	BottomLeft=0, BottomRight=1,
+	TopLeft=2, TopRight=3,
+
+	/** Added names for 3D */
+	BottomLeftFloor=0, BottomRightFloor=1,
+	TopLeftFloor=2, TopRightFloor=3,
+	BottomLeftCeil=4, BottomRightCeil=5,
+	TopLeftCeil=6, TopRightCeil=7
+  } CornerType;
+
 
   /** Default constructor initializing a null box. */
   inline explicit AlignedBox()
@@ -53,13 +72,15 @@
 
   /** Constructs a null box with \a _dim the dimension of the ambient space. */
   inline explicit AlignedBox(int _dim) : m_min(_dim), m_max(_dim)
-  { setNull(); }
+  { setEmpty(); }
 
   /** Constructs a box with extremities \a _min and \a _max. */
-  inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_min(_min), m_max(_max) {}
+  template<typename OtherVectorType1, typename OtherVectorType2>
+  inline AlignedBox(const OtherVectorType1& _min, const OtherVectorType2& _max) : m_min(_min), m_max(_max) {}
 
   /** Constructs a box containing a single point \a p. */
-  inline explicit AlignedBox(const VectorType& p) : m_min(p), m_max(p) {}
+  template<typename OtherVectorType>
+  inline explicit AlignedBox(const OtherVectorType& p) : m_min(p), m_max(p) {}
 
   ~AlignedBox() {}
 
@@ -67,14 +88,21 @@
   inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : AmbientDimAtCompileTime; }
 
   /** \returns true if the box is null, i.e, empty. */
-  inline bool isNull() const { return (m_min.cwise() > m_max).any(); }
+  inline bool isNull() const { return isEmpty(); }
 
   /** Makes \c *this a null/empty box. */
-  inline void setNull()
+  inline void setNull() { setEmpty(); }
+
+  /** \returns true if the box is empty. */
+  inline bool isEmpty() const { return (m_min.cwise() > m_max).any(); }
+
+  /** Makes \c *this an empty box. */
+  inline void setEmpty()
   {
     m_min.setConstant( std::numeric_limits<Scalar>::max());
     m_max.setConstant(-std::numeric_limits<Scalar>::max());
   }
+
 
   /** \returns the minimal corner */
   inline const VectorType& min() const { return m_min; }
@@ -86,65 +114,158 @@
   inline VectorType& max() { return m_max; }
 
   /** \returns the center of the box */
-  inline VectorType center() const { return (m_min + m_max) / 2; }
+  inline
+  const CwiseUnaryOp< ei_scalar_quotient1_op<Scalar>, VectorType > center() const
+  {
+	  typedef ei_scalar_quotient1_op<Scalar>                 Div2_t;
+	  typedef CwiseUnaryOp< Div2_t, VectorType >             DivOp_t;
+	  typedef ei_scalar_sum_op<Scalar>                       Sum_t;
+	  typedef CwiseBinaryOp< Sum_t, VectorType, VectorType > SumOp_t;
+
+	  return DivOp_t( SumOp_t( m_min, m_max ), Div2_t(2) );
+  }
+
+  /** \returns the lengths of the sides of the bounding box.
+   * Note that this function does not get the same
+   * result for integral or floating scalar types: see
+   */
+  inline
+  const CwiseBinaryOp< ei_scalar_difference_op<Scalar>, VectorType, VectorType> sizes() const
+  {
+	  typedef ei_scalar_difference_op<Scalar>                       Diff_t;
+	  typedef CwiseBinaryOp< Diff_t, VectorType, VectorType >       DiffOp_t;
+	  return DiffOp_t( m_max, m_min );
+  }
+
+  /** \returns the volume of the bounding box */
+  inline Scalar volume() const {
+	  return sizes().prod(); }
+
+  /** \returns an expression for the bounding box diagonal vector
+   * if the length of the diagonal is needed: diagonal().norm()
+   * will provide it.
+   * */
+  inline CwiseBinaryOp< ei_scalar_difference_op<Scalar>, VectorType, VectorType> diagonal() const {
+	  return sizes(); }
+
+  /** \returns the vertex of the bounding box at the corner defined by
+   * the corner-id corner. It works only for a 1D, 2D or 3D bounding box.
+   * For 1D bounding boxes corners are named by 2 enum constants:
+   * BottomLeft and BottomRight.
+   * For 2D bounding boxes, corners are named by 4 enum constants:
+   * BottomLeft, BottomRight, TopLeft, TopRight.
+   * For 3D bounding boxes, the following names are added:
+   * BottomLeftCeil, BottomRightCeil, TopLeftCeil, TopRightCeil
+   * */
+  inline
+  VectorType corner( CornerType corner ) const
+  {
+	  EIGEN_STATIC_ASSERT( _AmbientDim <= 3,
+			  THIS_METHOD_IS_ONLY_FOR_VECTORS_OF_A_SPECIFIC_SIZE );
+
+	  VectorType res;
+
+	  int mult=1;
+	  for( int d=0; d<dim(); ++d )
+	  {
+		  if( mult & corner ){
+			  res[d] = m_max[d]; }
+		  else{
+			  res[d] = m_min[d]; }
+		  mult*=2;
+	  }
+	  return res;
+  }
+
+  /** \returns a random point inside the bounding box sampled with
+   * a uniform distribution */
+  inline
+  VectorType sample() const
+  {
+	  VectorType r;
+	  for( int d=0; d<dim(); ++d )
+	  {
+		  if( NumTraits<Scalar>::HasFloatingPoint )
+		  {
+			  r[d] = m_min[d] + (m_max[d]-m_min[d])
+					  * ( ei_random<Scalar>() + ei_random_amplitude<Scalar>() )
+					  / (2*ei_random_amplitude<Scalar>() );
+		  }
+		  else{
+			  r[d] = ei_random( m_min[d], m_max[d] ); }
+	  }
+	  return r;
+  }
 
   /** \returns true if the point \a p is inside the box \c *this. */
-  inline bool contains(const VectorType& p) const
+  template<typename Derived>
+  inline bool contains(const MatrixBase<Derived>& p) const
   { return (m_min.cwise()<=p).all() && (p.cwise()<=m_max).all(); }
 
   /** \returns true if the box \a b is entirely inside the box \c *this. */
-  inline bool contains(const AlignedBox& b) const
+  template<typename Os, int Od>
+  inline bool contains(const AlignedBox<Os,Od>& b) const
   { return (m_min.cwise()<=b.min()).all() && (b.max().cwise()<=m_max).all(); }
 
   /** Extends \c *this such that it contains the point \a p and returns a reference to \c *this. */
-  inline AlignedBox& extend(const VectorType& p)
+  template<typename Derived>
+  inline AlignedBox& extend(const MatrixBase<Derived>& p)
   { m_min = m_min.cwise().min(p); m_max = m_max.cwise().max(p); return *this; }
 
   /** Extends \c *this such that it contains the box \a b and returns a reference to \c *this. */
-  inline AlignedBox& extend(const AlignedBox& b)
+  template<typename Os, int Od>
+  inline AlignedBox& extend(const AlignedBox<Os,Od>& b)
   { m_min = m_min.cwise().min(b.m_min); m_max = m_max.cwise().max(b.m_max); return *this; }
 
   /** Clamps \c *this by the box \a b and returns a reference to \c *this. */
-  inline AlignedBox& clamp(const AlignedBox& b)
+  template<typename Os, int Od>
+  inline AlignedBox& clamp(const AlignedBox<Os,Od>& b)
   { m_min = m_min.cwise().max(b.m_min); m_max = m_max.cwise().min(b.m_max); return *this; }
 
   /** Returns an AlignedBox that is the intersection of \a b and \c *this */
-  inline AlignedBox intersection(const AlignedBox &b) const
+  template<typename Os, int Od>
+  inline AlignedBox intersection(const AlignedBox<Os,Od> &b) const
   { return AlignedBox(m_min.cwise().max(b.m_min), m_max.cwise().min(b.m_max)); }
 
   /** Returns an AlignedBox that is the union of \a b and \c *this */
-  inline AlignedBox merged(const AlignedBox &b) const
+  template<typename Os, int Od>
+  inline AlignedBox merged(const AlignedBox<Os,Od> &b) const
   { return AlignedBox(m_min.cwise().min(b.m_min), m_max.cwise().max(b.m_max)); }
 
   /** Translate \c *this by the vector \a t and returns a reference to \c *this. */
-  inline AlignedBox& translate(const VectorType& t)
+  template<typename Derived>
+  inline AlignedBox& translate(const MatrixBase<Derived>& t)
   { m_min += t; m_max += t; return *this; }
 
   /** \returns the squared distance between the point \a p and the box \c *this,
     * and zero if \a p is inside the box.
     * \sa exteriorDistance()
     */
-  inline Scalar squaredExteriorDistance(const VectorType& p) const;
+  template<typename Derived>
+  inline Scalar squaredExteriorDistance(const MatrixBase<Derived>& p) const;
 
   /** \returns the squared distance between the boxes \a b and \c *this,
     * and zero if the boxes intersect.
     * \sa exteriorDistance()
     */
-  inline Scalar squaredExteriorDistance(const AlignedBox& b) const;
+  template<typename Os,int Od>
+  inline Scalar squaredExteriorDistance(const AlignedBox<Os,Od>& b) const;
 
   /** \returns the distance between the point \a p and the box \c *this,
     * and zero if \a p is inside the box.
     * \sa squaredExteriorDistance()
     */
-  inline Scalar exteriorDistance(const VectorType& p) const
-  { return ei_sqrt(squaredExteriorDistance(p)); }
+  template<typename Derived>
+  inline FloatingPoint exteriorDistance(const MatrixBase<Derived>& p) const
+  { return ei_sqrt(FloatingPoint(squaredExteriorDistance(p))); }
 
   /** \returns the distance between the boxes \a b and \c *this,
     * and zero if the boxes intersect.
     * \sa squaredExteriorDistance()
     */
-  inline Scalar exteriorDistance(const AlignedBox& b) const
-  { return ei_sqrt(squaredExteriorDistance(b)); }
+  template<typename Os, int Od>
+  inline FloatingPoint exteriorDistance(const AlignedBox<Os,Od>& b) const
+  { return ei_sqrt(FloatingPoint(squaredExteriorDistance(b))); }
 
   /** \returns \c *this with scalar type casted to \a NewScalarType
     *
@@ -171,7 +292,7 @@
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const
+  bool isApprox(const AlignedBox& other, RealScalar prec = dummy_precision<Scalar>()) const
   { return m_min.isApprox(other.m_min, prec) && m_max.isApprox(other.m_max, prec); }
 
 protected:
@@ -179,32 +300,48 @@
   VectorType m_min, m_max;
 };
 
+
+
 template<typename Scalar,int AmbiantDim>
-inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const VectorType& p) const
+template<typename Derived>
+inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const MatrixBase<Derived>& p) const
 {
   Scalar dist2 = 0.;
   Scalar aux;
   for (int k=0; k<dim(); ++k)
   {
-    if ((aux = (p[k]-m_min[k]))<0.)
-      dist2 += aux*aux;
-    else if ( (aux = (m_max[k]-p[k]))<0. )
-      dist2 += aux*aux;
+    if( m_min[k] > p[k] )
+	{
+		aux = m_min[k] - p[k];
+		dist2 += aux*aux;
+	}
+    else if( p[k] > m_max[k] )
+	{
+		aux = p[k] - m_max[k];
+		dist2 += aux*aux;
+	}
   }
   return dist2;
 }
 
 template<typename Scalar,int AmbiantDim>
-inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const AlignedBox& b) const
+template<typename Os, int Od>
+inline Scalar AlignedBox<Scalar,AmbiantDim>::squaredExteriorDistance(const AlignedBox<Os,Od>& b) const
 {
   Scalar dist2 = 0.;
   Scalar aux;
   for (int k=0; k<dim(); ++k)
   {
-    if ((aux = (b.m_min[k]-m_max[k]))>0.)
-      dist2 += aux*aux;
-    else if ( (aux = (m_min[k]-b.m_max[k]))>0. )
-      dist2 += aux*aux;
+	  if( m_min[k] > b.m_max[k] )
+	  {
+		  aux = m_min[k] - b.m_max[k];
+		  dist2 += aux*aux;
+	  }
+	  else if( b.m_min[k] > m_max[k] )
+	  {
+		  aux = b.m_min[k] - m_max[k];
+		  dist2 += aux*aux;
+	  }
   }
   return dist2;
 }

# HG changeset patch
# User Manuel Yguel <manuel.yguel@xxxxxxxxx>
# Date 1261425818 -3600
# Node ID 4e4fe8c87fbbebe6431c3ab99ad339d71bf2b3c6
# Parent  2d0813a2fb84155fa231076c992ae6370914f7f4
Tests corresponding to the adding functions: sizes, diagonal, volume, sample

diff -r 2d0813a2fb84 -r 4e4fe8c87fbb test/geo_alignedbox.cpp
--- a/test/geo_alignedbox.cpp	Mon Dec 21 19:24:26 2009 +0100
+++ b/test/geo_alignedbox.cpp	Mon Dec 21 21:03:38 2009 +0100
@@ -27,6 +27,9 @@
 #include <Eigen/LU>
 #include <Eigen/QR>
 
+#include<iostream>
+using namespace std;
+
 template<typename BoxType> void alignedbox(const BoxType& _box)
 {
   /* this test covers the following files:
@@ -40,6 +43,8 @@
 
   VectorType p0 = VectorType::Random(dim);
   VectorType p1 = VectorType::Random(dim);
+  while( p1 == p0 ){
+      p1 =  VectorType::Random(dim); }
   RealScalar s1 = ei_random<RealScalar>(0,1);
 
   BoxType b0(dim);
@@ -49,20 +54,13 @@
   b0.extend(p0);
   b0.extend(p1);
   VERIFY(b0.contains(p0*s1+(Scalar(1)-s1)*p1));
-  VERIFY(!b0.contains(p0 + (1+s1)*(p1-p0)));
+  VERIFY(!b0.contains(p0 + (2+s1)*(p1-p0)));
 
   (b2 = b0).extend(b1);
   VERIFY(b2.contains(b0));
   VERIFY(b2.contains(b1));
   VERIFY_IS_APPROX(b2.clamp(b0), b0);
 
-  // casting
-  const int Dim = BoxType::AmbientDimAtCompileTime;
-  typedef typename GetDifferentType<Scalar>::type OtherScalar;
-  AlignedBox<OtherScalar,Dim> hp1f = b0.template cast<OtherScalar>();
-  VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),b0);
-  AlignedBox<Scalar,Dim> hp1d = b0.template cast<Scalar>();
-  VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),b0);
 
   // alignment -- make sure there is no memory alignment assertion
   BoxType *bp0 = new BoxType(dim);
@@ -70,13 +68,117 @@
   bp0->extend(*bp1);
   delete bp0;
   delete bp1;
+
+  // sampling
+  for( int i=0; i<10; ++i )
+  {
+      VectorType r = b0.sample();
+      VERIFY(b0.contains(r));
+  }
+
 }
+
+
+
+template<typename BoxType>
+void alignedboxCastTests(const BoxType& _box)
+{
+  // casting
+  const int dim = _box.dim();
+  typedef typename BoxType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef Matrix<Scalar, BoxType::AmbientDimAtCompileTime, 1> VectorType;
+
+  VectorType p0 = VectorType::Random(dim);
+  VectorType p1 = VectorType::Random(dim);
+
+  BoxType b0(dim);
+
+  b0.extend(p0);
+  b0.extend(p1);
+
+  const int Dim = BoxType::AmbientDimAtCompileTime;
+  typedef typename GetDifferentType<Scalar>::type OtherScalar;
+  AlignedBox<OtherScalar,Dim> hp1f = b0.template cast<OtherScalar>();
+  VERIFY_IS_APPROX(hp1f.template cast<Scalar>(),b0);
+  AlignedBox<Scalar,Dim> hp1d = b0.template cast<Scalar>();
+  VERIFY_IS_APPROX(hp1d.template cast<Scalar>(),b0);
+}
+
+
+void specificTest1()
+{
+    Vector2f m; m << -1.0f, -2.0f;
+    Vector2f M; M <<  1.0f,  5.0f;
+
+    typedef AlignedBox<float,2>  BoxType;
+    BoxType box( m, M );
+
+    Vector2f sides = M-m;
+    VERIFY_IS_APPROX(sides, box.sizes() );
+    VERIFY_IS_APPROX(sides[1], box.sizes()[1] );
+    VERIFY_IS_APPROX(sides[1], box.sizes().maxCoeff() );
+    VERIFY_IS_APPROX(sides[0], box.sizes().minCoeff() );
+
+    VERIFY_IS_APPROX( 14.0f, box.volume() );
+    VERIFY_IS_APPROX( 53.0f, box.diagonal().squaredNorm() );
+    VERIFY_IS_APPROX( ei_sqrt( 53.0f ), box.diagonal().norm() );
+
+    VERIFY_IS_APPROX( m, box.corner( BoxType::BottomLeft ) );
+    VERIFY_IS_APPROX( M, box.corner( BoxType::TopRight ) );
+    Vector2f bottomRight; bottomRight << M[0], m[1];
+    Vector2f topLeft; topLeft << m[0], M[1];
+    VERIFY_IS_APPROX( bottomRight, box.corner( BoxType::BottomRight ) );
+    VERIFY_IS_APPROX( topLeft, box.corner( BoxType::TopLeft ) );
+}
+
+
+void specificTest2()
+{
+    Vector3i m; m << -1, -2, 0;
+    Vector3i M; M <<  1,  5, 3;
+
+    typedef AlignedBox<int,3>  BoxType;
+    BoxType box( m, M );
+
+    Vector3i sides = M-m;
+    VERIFY_IS_APPROX(sides, box.sizes() );
+    VERIFY_IS_APPROX(sides[1], box.sizes()[1] );
+    VERIFY_IS_APPROX(sides[1], box.sizes().maxCoeff() );
+    VERIFY_IS_APPROX(sides[0], box.sizes().minCoeff() );
+
+    VERIFY_IS_APPROX( 42, box.volume() );
+    VERIFY_IS_APPROX( 62, box.diagonal().squaredNorm() );
+
+    VERIFY_IS_APPROX( m, box.corner( BoxType::BottomLeftFloor ) );
+    VERIFY_IS_APPROX( M, box.corner( BoxType::TopRightCeil ) );
+    Vector3i bottomRightFloor; bottomRightFloor << M[0], m[1], m[2];
+    Vector3i topLeftFloor; topLeftFloor << m[0], M[1], m[2];
+    VERIFY_IS_APPROX( bottomRightFloor, box.corner( BoxType::BottomRightFloor ) );
+    VERIFY_IS_APPROX( topLeftFloor, box.corner( BoxType::TopLeftFloor ) );
+}
+
 
 void test_geo_alignedbox()
 {
-  for(int i = 0; i < g_repeat; i++) {
+  for(int i = 0; i < g_repeat; i++)
+  {
     CALL_SUBTEST_1( alignedbox(AlignedBox<float,2>()) );
-    CALL_SUBTEST_2( alignedbox(AlignedBox<float,3>()) );
-    CALL_SUBTEST_3( alignedbox(AlignedBox<double,4>()) );
+    CALL_SUBTEST_2( alignedboxCastTests(AlignedBox<float,2>()) );
+
+    CALL_SUBTEST_3( alignedbox(AlignedBox<float,3>()) );
+    CALL_SUBTEST_4( alignedboxCastTests(AlignedBox<float,3>()) );
+
+    CALL_SUBTEST_5( alignedbox(AlignedBox<double,4>()) );
+    CALL_SUBTEST_6( alignedboxCastTests(AlignedBox<double,4>()) );
+
+    CALL_SUBTEST_7( alignedbox(AlignedBox<double,1>()) );
+    CALL_SUBTEST_8( alignedboxCastTests(AlignedBox<double,1>()) );
+
+    CALL_SUBTEST_9( alignedbox(AlignedBox<int,1>()) );
+    CALL_SUBTEST_10( alignedbox(AlignedBox<int,2>()) );
+    CALL_SUBTEST_11( alignedbox(AlignedBox<int,3>()) );
   }
+  CALL_SUBTEST_12( specificTest1() );
+  CALL_SUBTEST_13( specificTest2() );
 }