moab
moab::OrientedBox Class Reference

Oriented bounding box. More...

#include <OrientedBox.hpp>

List of all members.

Classes

struct  CovarienceData

Public Member Functions

 OrientedBox ()
 OrientedBox (const CartVect axis[3], const CartVect &center)
double inner_radius () const
 radius of inscribed sphere
double outer_radius () const
 radius of circumscribed sphere
double outer_radius_squared (const double reps) const
 square of (radius+at least epsilon) of circumsphere
double inner_radius_squared (const double reps) const
 square of (radius-epsilon) of inscribed sphere
double volume () const
 volume of box
CartVect dimensions () const
 number of dimensions for which box is not flat
double area () const
 largest side area
CartVect scaled_axis (int index) const
 get vector in direction of axis, from box center to face
bool contained (const CartVect &point, double tolerance) const
bool intersect_ray (const CartVect &ray_start_point, const CartVect &ray_unit_direction, const double distance_tolerance, const double *nonnegatve_ray_len=0, const double *negative_ray_len=0) const
void closest_location_in_box (const CartVect &input_position, CartVect &output_position) const
 Find closest position on/within box to input position.
ErrorCode make_hex (EntityHandle &hex, Interface *instance)
 Construct a hexahedral element with the same shape as this box.

Static Public Member Functions

static ErrorCode tag_handle (Tag &handle_out, Interface *instance, const char *name)
 get tag handle for storing oriented box
static ErrorCode compute_from_vertices (OrientedBox &result, Interface *instance, const Range &vertices)
 Calculate an oriented box from a set of vertices.
static ErrorCode compute_from_2d_cells (OrientedBox &result, Interface *instance, const Range &elements)
 Calculate an oriented box from a set of 2D elements.
static ErrorCode covariance_data_from_tris (CovarienceData &result, Interface *moab_instance, const Range &elements)
static ErrorCode compute_from_covariance_data (OrientedBox &result, Interface *moab_instance, const CovarienceData *orient_array, unsigned orient_array_length, const Range &vertices)
static ErrorCode compute_from_covariance_data (OrientedBox &result, Interface *moab_instance, CovarienceData &orientation_data, const Range &vertices)

Public Attributes

CartVect center
 Box center.
CartVect axis [3]
 Box axes, unit vectors sorted by extent of box along axis.
CartVect length
 distance from center to plane along each axis
double radius
 outer radius (1/2 diagonal length) of box

Detailed Description

Oriented bounding box.

Definition at line 40 of file OrientedBox.hpp.


Constructor & Destructor Documentation

Definition at line 52 of file OrientedBox.hpp.

{}
moab::OrientedBox::OrientedBox ( const CartVect  axis[3],
const CartVect center 
)

Definition at line 94 of file OrientedBox.cpp.

  : center(mid)
{
    // re-order axes by length
  CartVect len( axes[0].length(), axes[1].length(), axes[2].length() );
  axis[0] = axes[0];
  axis[1] = axes[1];
  axis[2] = axes[2];
  
  if (len[2] < len[1])
  {
    if (len[2] < len[0]) {
      std::swap( len[0], len[2] );
      std::swap( axis[0], axis[2] );
    }
  }
  else if (len[1] < len[0]) {
    std::swap( len[0], len[1] );
    std::swap( axis[0], axis[1] );
  }
  if (len[1] > len[2]) {
    std::swap( len[1], len[2] );
    std::swap( axis[1], axis[2] );
  }
  
#if MB_ORIENTED_BOX_UNIT_VECTORS
  this->length = len;
  if (len[0] > 0.0)
    axis[0] /= len[0];
  if (len[1] > 0.0)
    axis[1] /= len[1];
  if (len[2] > 0.0)
    axis[2] /= len[2];
#endif

#if MB_ORIENTED_BOX_OUTER_RADIUS
  radius = len.length();
#endif
}

Member Function Documentation

double moab::OrientedBox::area ( ) const [inline]

largest side area

Definition at line 228 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return 4 * length[1] * length[2];
#else
  return 4 * (axis[1] * axis[2]).length();
#endif
}
void moab::OrientedBox::closest_location_in_box ( const CartVect input_position,
CartVect output_position 
) const

Find closest position on/within box to input position.

Find the closest position in the solid box to the input position. If the input position is on or within the box, then the output position will be the same as the input position. If the input position is outside the box, the outside position will be the closest point on the box boundary to the input position.

Definition at line 688 of file OrientedBox.cpp.

{
    // get coordinates on box axes
  const CartVect from_center = input_position - center;

#if MB_ORIENTED_BOX_UNIT_VECTORS
  CartVect local( from_center % axis[0],
                    from_center % axis[1],
                    from_center % axis[2] );

  for (int i = 0; i < 3; ++i) {
    if (local[i] < -length[i])
      local[i] = -length[i];
    else if (local[i] > length[i])
      local[i] =  length[i];
  }
#else
  CartVect local( (from_center % axis[0]) / (axis[0] % axis[0]),
                    (from_center % axis[1]) / (axis[1] % axis[1]),
                    (from_center % axis[2]) / (axis[2] % axis[2]) );

  for (int i = 0; i < 3; ++i) {
    if (local[i] < -1.0)
      local[i] = -1.0;
    else if (local[i] > 1.0)
      local[i] = 1.0;
  }
#endif

  output_position = center
                  + local[0] * axis[0] 
                  + local[1] * axis[1]
                  + local[2] * axis[2];
}
ErrorCode moab::OrientedBox::compute_from_2d_cells ( OrientedBox result,
Interface instance,
const Range elements 
) [static]

Calculate an oriented box from a set of 2D elements.

Definition at line 331 of file OrientedBox.cpp.

{
    // Get orientation data from elements
  CovarienceData data;
  ErrorCode rval = covariance_data_from_tris( data, instance, elements );
  if (MB_SUCCESS != rval)
    return rval;
  
    // get vertices from elements
  Range points;
  rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
  if (MB_SUCCESS != rval)
    return rval;
    
    // Calculate box given points and orientation data
  return compute_from_covariance_data( result, instance, data, points );
}
ErrorCode moab::OrientedBox::compute_from_covariance_data ( OrientedBox result,
Interface moab_instance,
const CovarienceData orient_array,
unsigned  orient_array_length,
const Range vertices 
) [static]

Calculate an OrientedBox given an arrray of CovarienceData and the list of vertices the box is to bound.

Definition at line 396 of file OrientedBox.cpp.

{
    // Sum input CovarienceData structures
  CovarienceData data_sum( Matrix3(0.0), CartVect(0.0), 0.0 );
  for (const CovarienceData* const end = data+data_length; data != end; ++data) {
    data_sum.matrix += data->matrix;
    data_sum.center += data->center;
    data_sum.area += data->area;
  }
    // Compute box from sum of structs
  return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
}
ErrorCode moab::OrientedBox::compute_from_covariance_data ( OrientedBox result,
Interface moab_instance,
CovarienceData orientation_data,
const Range vertices 
) [static]

Calculate an OrientedBox given a CovarienceData struct and the list of points the box is to bound.

Definition at line 351 of file OrientedBox.cpp.

{
  if (data.area <= 0.0) {
    CartVect axis[3] = { CartVect(0.), CartVect(0.), CartVect(0.) };
    result = OrientedBox( axis, CartVect(0.) );
    return MB_SUCCESS;
  }

    // get center from sum
  result.center = data.center / data.area;

    // get covariance matrix from moments
  data.matrix /= 12 * data.area;
  data.matrix -= outer_product( result.center, result.center );

    // get axes (Eigenvectors) from covariance matrix
  double lamda[3];
  moab::Matrix::EigenDecomp( data.matrix, lamda, result.axis );

    // We now have only the axes.  Calculate proper center
    // and extents for enclosed points.
  return box_from_axes( result, instance, vertices );
}      
ErrorCode moab::OrientedBox::compute_from_vertices ( OrientedBox result,
Interface instance,
const Range vertices 
) [static]

Calculate an oriented box from a set of vertices.

Definition at line 241 of file OrientedBox.cpp.

{
  const Range::iterator begin = vertices.lower_bound( MBVERTEX );
  const Range::iterator end = vertices.upper_bound( MBVERTEX );
  size_t count = 0;
  
    // compute mean
  CartVect v;
  result.center = CartVect( 0, 0, 0 );
  for (Range::iterator i = begin; i != end; ++i)
  {
    ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
    if (MB_SUCCESS != rval)
      return rval;
    result.center += v;
    ++count;
  }
  result.center /= count;
  
    // compute covariance matrix
  Matrix3 a( 0.0 );
  for (Range::iterator i = begin; i != end; ++i)
  {
    ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
    if (MB_SUCCESS != rval)
      return rval;
  
    v -= result.center;
    a += outer_product( v, v );
  }
  a /= count;

    // Get axes (Eigenvectors) from covariance matrix
  double lambda[3];
  moab::Matrix::EigenDecomp( a, lambda, result.axis );
  
    // Calculate center and extents of box given orientation defined by axes
  return box_from_axes( result, instance, vertices );
}
bool moab::OrientedBox::contained ( const CartVect point,
double  tolerance 
) const

Test if point is contained in box

Definition at line 379 of file OrientedBox.cpp.

{
  CartVect from_center = point - center;
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return fabs(from_center % axis[0]) - length[0] <= tol &&
         fabs(from_center % axis[1]) - length[1] <= tol &&
         fabs(from_center % axis[2]) - length[2] <= tol ;
#else
  for (int i = 0; i < 3; ++i) {
    double length = axis[i].length();
    if (fabs(from_center % axis[i]) - length*length > length*tol)
      return false;
  }
  return true;
#endif
}
ErrorCode moab::OrientedBox::covariance_data_from_tris ( CovarienceData result,
Interface moab_instance,
const Range elements 
) [static]

Calculate a CovarienceData struct from a list of triangles

Definition at line 283 of file OrientedBox.cpp.

{
  ErrorCode rval;
  const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
  const Range::iterator end = elements.lower_bound( CN::TypeDimensionMap[3].first );
  
    // compute mean and moments
  result.matrix = Matrix3(0.0);
  result.center = CartVect(0.0);
  result.area = 0.0;
  for (Range::iterator i = begin; i != end; ++i)
  {
    const EntityHandle* conn;
    int conn_len;
    rval = instance->get_connectivity( *i, conn, conn_len );
    if (MB_SUCCESS != rval)
      return rval;
    
      // for each triangle in the 2-D cell
    for (int j = 2; j < conn_len; ++j)
    {
      EntityHandle vertices[3] = { conn[0], conn[j-1], conn[j] };
      CartVect coords[3];
      rval = instance->get_coords( vertices, 3, coords[0].array() );
      if (MB_SUCCESS != rval)
        return rval;
      
        // edge vectors
      const CartVect edge0 = coords[1] - coords[0];
      const CartVect edge1 = coords[2] - coords[0];
      const CartVect centroid = (coords[0] + coords[1] + coords[2]) / 3;
      const double tri_area2 = (edge0 * edge1).length();
      result.area += tri_area2;
      result.center += tri_area2 * centroid;
      
      result.matrix += tri_area2 * (9 * outer_product( centroid,  centroid  ) +
                                    outer_product( coords[0], coords[0] ) +
                                    outer_product( coords[1], coords[1] ) +
                                    outer_product( coords[2], coords[2] ));
    } // for each triangle
  } // for each element

  return MB_SUCCESS;
}

number of dimensions for which box is not flat

Definition at line 219 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return 2.0 * length;
#else
  return 2.0 * CartVect( axis[0].length(), axis[1].length(), axis[2].length() );
#endif
}
double moab::OrientedBox::inner_radius ( ) const [inline]

radius of inscribed sphere

Definition at line 163 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return length[0];
#else
  return axis[0].length();
#endif
}
double moab::OrientedBox::inner_radius_squared ( const double  reps) const [inline]

square of (radius-epsilon) of inscribed sphere

Definition at line 199 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return (length[0]-reps) * (length[0]-reps);
#else
  CartVect tmp = axis[0];
  tmp -= CartVect(reps,reps,reps);
  return (tmp % tmp);
#endif
}
bool moab::OrientedBox::intersect_ray ( const CartVect ray_start_point,
const CartVect ray_unit_direction,
const double  distance_tolerance,
const double *  nonnegatve_ray_len = 0,
const double *  negative_ray_len = 0 
) const

Test for intersection of a ray (or line segment) with this box. Ray length limits are used to optimize Monte Carlo particle tracking.

Parameters:
ray_start_pointThe base point of the ray
ray_unit_directionThe direction of the ray (must be unit length)
distance_toleranceTolerance to use in intersection checks
nonnegative_ray_lenOptional length of ray in forward direction
negative_ray_lenOptional length of ray in reverse direction

Definition at line 510 of file OrientedBox.cpp.

{
  // test distance from box center to line
  const CartVect cx       = center - ray_origin;
  const double dist_s     = cx % ray_direction;
  const double dist_sq    = cx % cx - (dist_s*dist_s);
  const double max_diagsq = outer_radius_squared(reps);
  
  // For the largest sphere, no intersections exist if discriminant is negative.
  // Geometrically, if distance from box center to line is greater than the 
  // longest diagonal, there is no intersection.
  // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
  if(dist_sq > max_diagsq) return false;

  // If the closest possible intersection must be closer than nonneg_ray_len. Be 
  // careful with absolute value, squaring distances, and subtracting squared
  // distances.
  if (nonneg_ray_len) {
    assert(0<=*nonneg_ray_len);
    double max_len;
    if(neg_ray_len) {
      assert(0>=*neg_ray_len);
      max_len = std::max(*nonneg_ray_len,-(*neg_ray_len));
    } else {
      max_len = *nonneg_ray_len;
    }
    const double temp = fabs(dist_s) - max_len;
    if(0.0<temp && temp*temp>max_diagsq) return false;
  } 

  // if smaller than shortest diagonal, we do hit
  if (dist_sq < inner_radius_squared(reps)) {
    // nonnegative direction
    if(dist_s>=0.0 ) {
      if(nonneg_ray_len) {
        if(*nonneg_ray_len>dist_s) return true;
      } else {
        return true;
      }
    // negative direction 
    } else {
      if(neg_ray_len && *neg_ray_len<dist_s) return true;
    }
  }
  
    // get transpose of axes
    // Note: if axes were stored as a matrix, could skip
    // transpose and just switch order of operands in
    // matrix-vector multiplies below. - J.K.
  //Matrix3 B( axis[0][0], axis[1][0], axis[2][0],
  //             axis[0][1], axis[1][1], axis[2][1],
  //             axis[0][2], axis[1][2], axis[2][2] );
  Matrix3 B( axis[0][0], axis[0][1], axis[0][2],
               axis[1][0], axis[1][1], axis[1][2],
               axis[2][0], axis[2][1], axis[2][2] );
  //CartVect T = B * -center;
  
    // transform ray to box coordintae system
  //CartVect par_pos = T + B * b;
  CartVect par_pos = B * (ray_origin - center);
  CartVect par_dir = B * ray_direction;

  // Fast Rejection Test: Ray will not intersect if it is going away from the box.
  // This will not work for rays with neg_ray_len. length[0] is half of box width 
  // along axis[0].
  const double half_x = length[0] + reps;
  const double half_y = length[1] + reps;
  const double half_z = length[2] + reps;
  if(!neg_ray_len) {
    if ((par_pos[0] >  half_x && par_dir[0] >= 0) ||
    (par_pos[0] < -half_x && par_dir[0] <= 0))
      return false;
  
    if ((par_pos[1] >  half_y && par_dir[1] >= 0) ||
    (par_pos[1] < -half_y && par_dir[1] <= 0))
      return false;
    
    if ((par_pos[2] >  half_z && par_dir[2] >= 0) ||
    (par_pos[2] < -half_z && par_dir[2] <= 0))
      return false;
  }

  // test if ray_origin is inside box
  if (par_pos[0] <= half_x && par_pos[0] >= -half_x &&
      par_pos[1] <= half_y && par_pos[1] >= -half_y &&
      par_pos[2] <= half_z && par_pos[2] >= -half_z)
    return true;

    //test two xy plane
  if (fabs(par_dir[0] * (half_z - par_pos[2]) + par_dir[2] * par_pos[0]) 
        <= fabs(par_dir[2] * half_x) &&  // test against x extents using z
      fabs(par_dir[1] * (half_z - par_pos[2]) + par_dir[2] * par_pos[1]) 
        <= fabs(par_dir[2] * half_y) &&  // test against y extents using z
      check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[0] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[0]) 
        <= fabs(par_dir[2] * half_x) && 
      fabs(par_dir[1] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[1]) 
        <= fabs(par_dir[2] * half_y) &&
      check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
    return true;

    //test two xz plane
  if (fabs(par_dir[0] * (half_y - par_pos[1]) + par_dir[1] * par_pos[0]) 
        <= fabs(par_dir[1] * half_x) && 
      fabs(par_dir[2] * (half_y - par_pos[1]) + par_dir[1] * par_pos[2]) 
        <= fabs(par_dir[1] * half_z) &&
      check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[0] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[0]) 
        <= fabs(par_dir[1] * half_x) && 
      fabs(par_dir[2] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[2])
        <= fabs(par_dir[1] * half_z) &&
      check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
    return true;

    //test two yz plane
  if (fabs(par_dir[1] * (half_x - par_pos[0]) + par_dir[0] * par_pos[1]) 
        <= fabs(par_dir[0] * half_y) &&
      fabs(par_dir[2] * (half_x - par_pos[0]) + par_dir[0] * par_pos[2]) 
        <= fabs(par_dir[0] * half_z) &&
      check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[1] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[1])
        <= fabs(par_dir[0] * half_y) &&
      fabs(par_dir[2] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[2]) 
        <= fabs(par_dir[0] * half_z) &&
      check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
    return true;

  return false;
}

Construct a hexahedral element with the same shape as this box.

Definition at line 647 of file OrientedBox.cpp.

{
  ErrorCode rval;
  int signs[8][3] = { { -1, -1, -1 },
                      {  1, -1, -1 },
                      {  1,  1, -1 },
                      { -1,  1, -1 },
                      { -1, -1,  1 },
                      {  1, -1,  1 },
                      {  1,  1,  1 },
                      { -1,  1,  1 } };
                      
  std::vector<EntityHandle> vertices;
  for (int i = 0; i < 8; ++i)
  {
    CartVect coords(center);
    for (int j = 0; j < 3; ++j){
#if MB_ORIENTED_BOX_UNIT_VECTORS
      coords += signs[i][j] * (axis[j]*length[j]);
#else
      coords += signs[i][j] * axis[j];
#endif
    }
    EntityHandle handle;
    rval = instance->create_vertex( coords.array(), handle );
    if (MB_SUCCESS != rval) {
      instance->delete_entities( &vertices[0], vertices.size() );
      return rval;
    }
    vertices.push_back( handle );
  }
  
  rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
  if (MB_SUCCESS != rval) {
    instance->delete_entities( &vertices[0], vertices.size() );
    return rval;
  }
  
  return MB_SUCCESS;
}
double moab::OrientedBox::outer_radius ( ) const [inline]

radius of circumscribed sphere

Definition at line 172 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_OUTER_RADIUS
  return radius;
#elif MB_ORIENTED_BOX_UNIT_VECTORS
  return length.length();
#else
  return (axis[0] + axis[1] + axis[2]).length();
#endif
}
double moab::OrientedBox::outer_radius_squared ( const double  reps) const [inline]

square of (radius+at least epsilon) of circumsphere

Definition at line 184 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_OUTER_RADIUS
  return (radius+reps)*(radius+reps);
#elif MB_ORIENTED_BOX_UNIT_VECTORS
  CartVect tmp(length[0]+reps,length[1]+reps,length[2]+reps);
  return tmp % tmp;
#else
  CartVect half_diag = axis[0] + axis[1] + axis[2];
  half_diag += CartVect(reps,reps,reps);
  return half_diag % half_diag;
#endif
}
CartVect moab::OrientedBox::scaled_axis ( int  index) const [inline]

get vector in direction of axis, from box center to face

Definition at line 237 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return length[index] * axis[index];
#else
  return axis[index];
#endif
}
ErrorCode moab::OrientedBox::tag_handle ( Tag handle_out,
Interface instance,
const char *  name 
) [static]

get tag handle for storing oriented box

Get the handle for the tag with the specified name and check that the tag is appropriate for storing instances of OrientedBox. The resulting tag may be used to store instances of OrientedBox directly.

Parameters:
handle_outThe TagHandle, passed back to caller
nameThe tag name
createIf true, tag will be created if it does not exist

Definition at line 134 of file OrientedBox.cpp.

{
    // We're going to assume this when mapping the OrientedBox
    // to tag data, so assert it.  
#if MB_ORIENTED_BOX_OUTER_RADIUS
  const int rad_size = 1;
#else
  const int rad_size = 0;
#endif
#if MB_ORIENTED_BOX_UNIT_VECTORS
  const int SIZE = rad_size + 15;
#else
  const int SIZE = rad_size + 12;
#endif
  assert( sizeof(OrientedBox) == SIZE*sizeof(double) );
  
  return instance->tag_get_handle( name, SIZE, MB_TYPE_DOUBLE,
                                   handle_out, MB_TAG_DENSE|MB_TAG_CREAT );
}
double moab::OrientedBox::volume ( ) const [inline]

volume of box

Definition at line 210 of file OrientedBox.hpp.

{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return 8 * length[0] * length[1] * length[2];
#else
  return fabs(8 * axis[0] % (axis[1] * axis[2]));
#endif
}

Member Data Documentation

Box axes, unit vectors sorted by extent of box along axis.

Definition at line 44 of file OrientedBox.hpp.

Box center.

Definition at line 43 of file OrientedBox.hpp.

distance from center to plane along each axis

Definition at line 46 of file OrientedBox.hpp.

outer radius (1/2 diagonal length) of box

Definition at line 49 of file OrientedBox.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines