MOAB: Mesh Oriented datABase  (version 5.2.1)
OrientedBox.cpp
Go to the documentation of this file.
00001 /*
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 /*
00017  * The algorithms for the calculation of the oriented box from a
00018  * set of points or a set of cells was copied from the implementation
00019  " in the "Visualization Toolkit".  J.K. - 2006-07-19
00020  *
00021  * Program:   Visualization Toolkit
00022  * Module:    $RCSfile$
00023  *
00024  * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
00025  * All rights reserved.
00026  * See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
00027  */
00028 
00029 /**\file OrientedBox.cpp
00030  *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
00031  *\date 2006-07-18
00032  */
00033 
00034 #include "moab/Interface.hpp"
00035 #include "moab/CN.hpp"
00036 #include "moab/OrientedBox.hpp"
00037 #include "moab/Range.hpp"
00038 #include <ostream>
00039 #include <assert.h>
00040 #include <limits>
00041 
00042 namespace moab
00043 {
00044 
00045 std::ostream& operator<<( std::ostream& s, const OrientedBox& b )
00046 {
00047     return s << b.center << " + " << b.axes.col( 0 )
00048 #if MB_ORIENTED_BOX_UNIT_VECTORS
00049              << ":" << b.length[0]
00050 #endif
00051              << " x " << b.axes.col( 1 )
00052 #if MB_ORIENTED_BOX_UNIT_VECTORS
00053              << ":" << b.length[1]
00054 #endif
00055              << " x " << b.axes.col( 2 )
00056 #if MB_ORIENTED_BOX_UNIT_VECTORS
00057              << ":" << b.length[2]
00058 #endif
00059         ;
00060 }
00061 
00062 /**\brief Find closest point on line
00063  *
00064  * Find the point on the line for which a line trough the
00065  * input point \a p and the result position is orthogonal to
00066  * the input line.
00067  * \param p  The point for which to find the perpendicular
00068  * \param b  A point on the line
00069  * \param m  The direction of the line
00070  * \return   The location on the line specified as 't' in the
00071  *           formula t * m + b
00072  */
00073 static double point_perp( const CartVect& p,   // closest to this point
00074                           const CartVect& b,   // point on line
00075                           const CartVect& m )  // line direction
00076 {
00077 #if MB_ORIENTED_BOX_UNIT_VECTORS
00078     double t = ( m % ( p - b ) );
00079 #else
00080     double t           = ( m % ( p - b ) ) / ( m % m );
00081 #endif
00082     return Util::is_finite( t ) ? t : 0.0;
00083 }
00084 
00085 void OrientedBox::order_axes_by_length( double ax1_len, double ax2_len, double ax3_len )
00086 {
00087 
00088     CartVect len( ax1_len, ax2_len, ax3_len );
00089 
00090     if( len[2] < len[1] )
00091     {
00092         if( len[2] < len[0] )
00093         {
00094             std::swap( len[0], len[2] );
00095             axes.swapcol( 0, 2 );
00096         }
00097     }
00098     else if( len[1] < len[0] )
00099     {
00100         std::swap( len[0], len[1] );
00101         axes.swapcol( 0, 1 );
00102     }
00103     if( len[1] > len[2] )
00104     {
00105         std::swap( len[1], len[2] );
00106         axes.swapcol( 1, 2 );
00107     }
00108 
00109 #if MB_ORIENTED_BOX_UNIT_VECTORS
00110     length = len;
00111     if( len[0] > 0.0 ) axes.colscale( 0, 1.0 / len[0] );
00112     if( len[1] > 0.0 ) axes.colscale( 1, 1.0 / len[1] );
00113     if( len[2] > 0.0 ) axes.colscale( 2, 1.0 / len[2] );
00114 #endif
00115 
00116 #if MB_ORIENTED_BOX_OUTER_RADIUS
00117     radius = len.length();
00118 #endif
00119 }
00120 
00121 OrientedBox::OrientedBox( const CartVect axes_in[3], const CartVect& mid ) : center( mid )
00122 {
00123 
00124     axes = Matrix3( axes_in[0], axes_in[1], axes_in[2], false );
00125 
00126     order_axes_by_length( axes_in[0].length(), axes_in[1].length(), axes_in[2].length() );
00127 }
00128 
00129 OrientedBox::OrientedBox( const Matrix3& axes_mat, const CartVect& mid ) : center( mid ), axes( axes_mat )
00130 {
00131     order_axes_by_length( axes.col( 0 ).length(), axes.col( 1 ).length(), axes.col( 2 ).length() );
00132 }
00133 
00134 ErrorCode OrientedBox::tag_handle( Tag& handle_out, Interface* instance, const char* name )
00135 {
00136     // We're going to assume this when mapping the OrientedBox
00137     // to tag data, so assert it.
00138 #if MB_ORIENTED_BOX_OUTER_RADIUS
00139     const int rad_size = 1;
00140 #else
00141     const int rad_size = 0;
00142 #endif
00143 #if MB_ORIENTED_BOX_UNIT_VECTORS
00144     const int SIZE = rad_size + 15;
00145 #else
00146     const int SIZE     = rad_size + 12;
00147 #endif
00148     assert( sizeof( OrientedBox ) == SIZE * sizeof( double ) );
00149 
00150     return instance->tag_get_handle( name, SIZE, MB_TYPE_DOUBLE, handle_out, MB_TAG_DENSE | MB_TAG_CREAT );
00151 }
00152 
00153 /**\brief Common code for box calculation
00154  *
00155  * Given the orientation of the box and an approximate center,
00156  * calculate the exact center and extents of the box.
00157  *
00158  *\param result.center  As input, the approximate center of the box.
00159  *                      As output, the exact center of the box.
00160  *\param result.axes    As input, directions of principal axes corresponding
00161  *                      to the orientation of the box.  Axes are assumed to
00162  *                      be unit-length on input.  Output will include extents
00163  *                      of box.
00164  *\param points  The set of points the box should contain.
00165  */
00166 static ErrorCode box_from_axes( OrientedBox& result, Interface* instance, const Range& points )
00167 {
00168     ErrorCode rval;
00169 
00170     // project points onto axes to get box extents
00171     CartVect min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() );
00172     for( Range::iterator i = points.begin(); i != points.end(); ++i )
00173     {
00174         CartVect coords;
00175         rval = instance->get_coords( &*i, 1, coords.array() );MB_CHK_ERR( rval );
00176 
00177         for( int d = 0; d < 3; ++d )
00178         {
00179             const double t = point_perp( coords, result.center, result.axes.col( d ) );
00180             if( t < min[d] ) min[d] = t;
00181             if( t > max[d] ) max[d] = t;
00182         }
00183     }
00184 
00185     // We now have a box defined by three orthogonal line segments
00186     // that intersect at the center of the box.  Each line segment
00187     // is defined as result.center + t * result.axes[i], where the
00188     // range of t is [min[i], max[i]].
00189 
00190     // Calculate new center
00191     const CartVect mid = 0.5 * ( min + max );
00192     result.center += mid[0] * result.axes.col( 0 ) + mid[1] * result.axes.col( 1 ) + mid[2] * result.axes.col( 2 );
00193 
00194     // reorder axes by length
00195     CartVect range = 0.5 * ( max - min );
00196     if( range[2] < range[1] )
00197     {
00198         if( range[2] < range[0] )
00199         {
00200             std::swap( range[0], range[2] );
00201             result.axes.swapcol( 0, 2 );
00202         }
00203     }
00204     else if( range[1] < range[0] )
00205     {
00206         std::swap( range[0], range[1] );
00207         result.axes.swapcol( 0, 1 );
00208     }
00209     if( range[1] > range[2] )
00210     {
00211         std::swap( range[1], range[2] );
00212         result.axes.swapcol( 1, 2 );
00213     }
00214 
00215     // scale axis to encompass all points, divide in half
00216 #if MB_ORIENTED_BOX_UNIT_VECTORS
00217     result.length = range;
00218 #else
00219     result.axes.colscale( 0, range[0] );
00220     result.axes.colscale( 1, range[1] );
00221     result.axes.colscale( 2, range[2] );
00222 #endif
00223 
00224 #if MB_ORIENTED_BOX_OUTER_RADIUS
00225     result.radius = range.length();
00226 #endif
00227 
00228     return MB_SUCCESS;
00229 }
00230 
00231 ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result, Interface* instance, const Range& vertices )
00232 {
00233     const Range::iterator begin = vertices.lower_bound( MBVERTEX );
00234     const Range::iterator end   = vertices.upper_bound( MBVERTEX );
00235     size_t count                = 0;
00236 
00237     // compute mean
00238     CartVect v;
00239     result.center = CartVect( 0, 0, 0 );
00240     for( Range::iterator i = begin; i != end; ++i )
00241     {
00242         ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
00243         if( MB_SUCCESS != rval ) return rval;
00244         result.center += v;
00245         ++count;
00246     }
00247     result.center /= count;
00248 
00249     // compute covariance matrix
00250     Matrix3 a( 0.0 );
00251     for( Range::iterator i = begin; i != end; ++i )
00252     {
00253         ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
00254         if( MB_SUCCESS != rval ) return rval;
00255 
00256         v -= result.center;
00257         a += outer_product( v, v );
00258     }
00259     a /= count;
00260 
00261     // Get axes (Eigenvectors) from covariance matrix
00262     CartVect lambda;
00263     a.eigen_decomposition( lambda, result.axes );
00264 
00265     // Calculate center and extents of box given orientation defined by axes
00266     return box_from_axes( result, instance, vertices );
00267 }
00268 
00269 ErrorCode OrientedBox::covariance_data_from_tris( CovarienceData& result, Interface* instance, const Range& elements )
00270 {
00271     ErrorCode rval;
00272     const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
00273     const Range::iterator end   = elements.lower_bound( CN::TypeDimensionMap[3].first );
00274 
00275     // compute mean and moments
00276     result.matrix = Matrix3( 0.0 );
00277     result.center = CartVect( 0.0 );
00278     result.area   = 0.0;
00279     for( Range::iterator i = begin; i != end; ++i )
00280     {
00281         const EntityHandle* conn = NULL;
00282         int conn_len             = 0;
00283         rval                     = instance->get_connectivity( *i, conn, conn_len );
00284         if( MB_SUCCESS != rval ) return rval;
00285 
00286         // for each triangle in the 2-D cell
00287         for( int j = 2; j < conn_len; ++j )
00288         {
00289             EntityHandle vertices[3] = { conn[0], conn[j - 1], conn[j] };
00290             CartVect coords[3];
00291             rval = instance->get_coords( vertices, 3, coords[0].array() );
00292             if( MB_SUCCESS != rval ) return rval;
00293 
00294             // edge vectors
00295             const CartVect edge0    = coords[1] - coords[0];
00296             const CartVect edge1    = coords[2] - coords[0];
00297             const CartVect centroid = ( coords[0] + coords[1] + coords[2] ) / 3;
00298             const double tri_area2  = ( edge0 * edge1 ).length();
00299             result.area += tri_area2;
00300             result.center += tri_area2 * centroid;
00301 
00302             result.matrix +=
00303                 tri_area2 * ( 9 * outer_product( centroid, centroid ) + outer_product( coords[0], coords[0] ) +
00304                               outer_product( coords[1], coords[1] ) + outer_product( coords[2], coords[2] ) );
00305         }  // for each triangle
00306     }      // for each element
00307 
00308     return MB_SUCCESS;
00309 }
00310 
00311 ErrorCode OrientedBox::compute_from_2d_cells( OrientedBox& result, Interface* instance, const Range& elements )
00312 {
00313     // Get orientation data from elements
00314     CovarienceData data;
00315     ErrorCode rval = covariance_data_from_tris( data, instance, elements );
00316     if( MB_SUCCESS != rval ) return rval;
00317 
00318     // get vertices from elements
00319     Range points;
00320     rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
00321     if( MB_SUCCESS != rval ) return rval;
00322 
00323     // Calculate box given points and orientation data
00324     return compute_from_covariance_data( result, instance, data, points );
00325 }
00326 
00327 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result, Interface* instance, CovarienceData& data,
00328                                                      const Range& vertices )
00329 {
00330     if( data.area <= 0.0 )
00331     {
00332         Matrix3 empty_axes( 0.0 );
00333         result = OrientedBox( empty_axes, CartVect( 0. ) );
00334         return MB_SUCCESS;
00335     }
00336 
00337     // get center from sum
00338     result.center = data.center / data.area;
00339 
00340     // get covariance matrix from moments
00341     data.matrix /= 12 * data.area;
00342     data.matrix -= outer_product( result.center, result.center );
00343 
00344     // get axes (Eigenvectors) from covariance matrix
00345     CartVect lambda;
00346     data.matrix.eigen_decomposition( lambda, result.axes );
00347 
00348     // We now have only the axes.  Calculate proper center
00349     // and extents for enclosed points.
00350     return box_from_axes( result, instance, vertices );
00351 }
00352 
00353 bool OrientedBox::contained( const CartVect& point, double tol ) const
00354 {
00355     CartVect from_center = point - center;
00356 #if MB_ORIENTED_BOX_UNIT_VECTORS
00357     return fabs( from_center % axes.col( 0 ) ) - length[0] <= tol &&
00358            fabs( from_center % axes.col( 1 ) ) - length[1] <= tol &&
00359            fabs( from_center % axes.col( 2 ) ) - length[2] <= tol;
00360 #else
00361     for( int i = 0; i < 3; ++i )
00362     {
00363         double length = axes.col( i ).length();
00364         if( fabs( from_center % axes.col( i ) ) - length * length > length * tol ) return false;
00365     }
00366     return true;
00367 #endif
00368 }
00369 
00370 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result, Interface* moab_instance,
00371                                                      const CovarienceData* data, unsigned data_length,
00372                                                      const Range& vertices )
00373 {
00374     // Sum input CovarienceData structures
00375     CovarienceData data_sum( Matrix3( 0.0 ), CartVect( 0.0 ), 0.0 );
00376     for( const CovarienceData* const end = data + data_length; data != end; ++data )
00377     {
00378         data_sum.matrix += data->matrix;
00379         data_sum.center += data->center;
00380         data_sum.area += data->area;
00381     }
00382     // Compute box from sum of structs
00383     return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
00384 }
00385 
00386 // bool OrientedBox::contained( const OrientedBox& box, double tol ) const
00387 //{
00388 //  for (int i = -1; i < 2; i += 2)
00389 //  {
00390 //    for (int j = -1; j < 2; j += 2)
00391 //    {
00392 //      for (int k = -1; k < 2; k += 2)
00393 //      {
00394 //        CartVect corner( center );
00395 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
00396 //        corner += i * box.length[0] * box.axis.col(0);
00397 //        corner += j * box.length[1] * box.axis.col(1);
00398 //        corner += k * box.length[2] * box.axis.col(2);
00399 //#else
00400 //        corner += i * box.axis.col(0);
00401 //        corner += j * box.axis.col(1);
00402 //        corner += k * box.axis.col(2);
00403 //#endif
00404 //        if (!contained( corner, tol ))
00405 //          return false;
00406 //      }
00407 //    }
00408 //  }
00409 //  return true;
00410 //}
00411 
00412 /* This is a helper function to check limits on ray length, turning the box-ray
00413  * intersection test into a box-segment intersection test. Use this to test the
00414  * limits against one side (plane) of the box. The side of the box (plane) is
00415  * normal to an axis.
00416  *
00417  *   normal_par_pos  Coordinate of particle's position along axis normal to side of box
00418  *   normal_par_dir  Coordinate of particle's direction along axis normal to side of box
00419  *   half_extent     Distance between center of box and side of box
00420  *   nonneg_ray_len  Maximum ray length in positive direction (in front of origin)
00421  *   neg_ray_len     Maximum ray length in negative direction (behind origin)
00422  *   return          true if intersection with plane occurs within distance limits
00423  *
00424  * ray equation:   intersection = origin + dist*direction
00425  * plane equation: intersection.plane_normal = half_extent
00426  *
00427  * Assume plane_normal and direction are unit vectors. Combine equations.
00428  *
00429  *     (origin + dist*direction).plane_normal = half_extent
00430  *     origin.plane_normal + dist*direction.plane_normal = half_extent
00431  *     dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
00432  *
00433  * Although this solves for distance, avoid floating point division.
00434  *
00435  *     dist*direction.plane_normal = half_extent - origin.plane_normal
00436  *
00437  * Use inequalities to test dist against ray length limits. Be aware that
00438  * inequalities change due to sign of direction.plane_normal.
00439  */
00440 inline bool check_ray_limits( const double normal_par_pos, const double normal_par_dir, const double half_extent,
00441                               const double* nonneg_ray_len, const double* neg_ray_len )
00442 {
00443 
00444     const double extent_pos_diff = half_extent - normal_par_pos;
00445 
00446     // limit in positive direction
00447     if( nonneg_ray_len )
00448     {  // should be 0 <= t <= nonneg_ray_len
00449         assert( 0 <= *nonneg_ray_len );
00450         if( normal_par_dir > 0 )
00451         {  // if/else if needed for pos/neg divisor
00452             if( *nonneg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff >= 0 ) return true;
00453         }
00454         else if( normal_par_dir < 0 )
00455         {
00456             if( *nonneg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff <= 0 ) return true;
00457         }
00458     }
00459     else
00460     {  // should be 0 <= t
00461         if( normal_par_dir > 0 )
00462         {  // if/else if needed for pos/neg divisor
00463             if( extent_pos_diff >= 0 ) return true;
00464         }
00465         else if( normal_par_dir < 0 )
00466         {
00467             if( extent_pos_diff <= 0 ) return true;
00468         }
00469     }
00470 
00471     // limit in negative direction
00472     if( neg_ray_len )
00473     {  // should be neg_ray_len <= t < 0
00474         assert( 0 >= *neg_ray_len );
00475         if( normal_par_dir > 0 )
00476         {  // if/else if needed for pos/neg divisor
00477             if( *neg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff < 0 ) return true;
00478         }
00479         else if( normal_par_dir < 0 )
00480         {
00481             if( *neg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff > 0 ) return true;
00482         }
00483     }
00484 
00485     return false;
00486 }
00487 
00488 /* This implementation copied from cgmMC (overlap.C).
00489  * Original author:  Tim Tautges?
00490  */
00491 bool OrientedBox::intersect_ray( const CartVect& ray_origin, const CartVect& ray_direction, const double reps,
00492                                  const double* nonneg_ray_len, const double* neg_ray_len ) const
00493 {
00494     // test distance from box center to line
00495     const CartVect cx       = center - ray_origin;
00496     const double dist_s     = cx % ray_direction;
00497     const double dist_sq    = cx % cx - ( dist_s * dist_s );
00498     const double max_diagsq = outer_radius_squared( reps );
00499 
00500     // For the largest sphere, no intersections exist if discriminant is negative.
00501     // Geometrically, if distance from box center to line is greater than the
00502     // longest diagonal, there is no intersection.
00503     // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
00504     if( dist_sq > max_diagsq ) return false;
00505 
00506     // If the closest possible intersection must be closer than nonneg_ray_len. Be
00507     // careful with absolute value, squaring distances, and subtracting squared
00508     // distances.
00509     if( nonneg_ray_len )
00510     {
00511         assert( 0 <= *nonneg_ray_len );
00512         double max_len;
00513         if( neg_ray_len )
00514         {
00515             assert( 0 >= *neg_ray_len );
00516             max_len = std::max( *nonneg_ray_len, -( *neg_ray_len ) );
00517         }
00518         else
00519         {
00520             max_len = *nonneg_ray_len;
00521         }
00522         const double temp = fabs( dist_s ) - max_len;
00523         if( 0.0 < temp && temp * temp > max_diagsq ) return false;
00524     }
00525 
00526     // if smaller than shortest diagonal, we do hit
00527     if( dist_sq < inner_radius_squared( reps ) )
00528     {
00529         // nonnegative direction
00530         if( dist_s >= 0.0 )
00531         {
00532             if( nonneg_ray_len )
00533             {
00534                 if( *nonneg_ray_len > dist_s ) return true;
00535             }
00536             else
00537             {
00538                 return true;
00539             }
00540             // negative direction
00541         }
00542         else
00543         {
00544             if( neg_ray_len && *neg_ray_len < dist_s ) return true;
00545         }
00546     }
00547 
00548     // get transpose of axes
00549     Matrix3 B = axes.transpose();
00550 
00551     // transform ray to box coordintae system
00552     CartVect par_pos = B * ( ray_origin - center );
00553     CartVect par_dir = B * ray_direction;
00554 
00555     // Fast Rejection Test: Ray will not intersect if it is going away from the box.
00556     // This will not work for rays with neg_ray_len. length[0] is half of box width
00557     // along axes.col(0).
00558     const double half_x = length[0] + reps;
00559     const double half_y = length[1] + reps;
00560     const double half_z = length[2] + reps;
00561     if( !neg_ray_len )
00562     {
00563         if( ( par_pos[0] > half_x && par_dir[0] >= 0 ) || ( par_pos[0] < -half_x && par_dir[0] <= 0 ) ) return false;
00564 
00565         if( ( par_pos[1] > half_y && par_dir[1] >= 0 ) || ( par_pos[1] < -half_y && par_dir[1] <= 0 ) ) return false;
00566 
00567         if( ( par_pos[2] > half_z && par_dir[2] >= 0 ) || ( par_pos[2] < -half_z && par_dir[2] <= 0 ) ) return false;
00568     }
00569 
00570     // test if ray_origin is inside box
00571     if( par_pos[0] <= half_x && par_pos[0] >= -half_x && par_pos[1] <= half_y && par_pos[1] >= -half_y &&
00572         par_pos[2] <= half_z && par_pos[2] >= -half_z )
00573         return true;
00574 
00575     // test two xy plane
00576     if( fabs( par_dir[0] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <=
00577             fabs( par_dir[2] * half_x ) &&  // test against x extents using z
00578         fabs( par_dir[1] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <=
00579             fabs( par_dir[2] * half_y ) &&  // test against y extents using z
00580         check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
00581         return true;
00582     if( fabs( par_dir[0] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <= fabs( par_dir[2] * half_x ) &&
00583         fabs( par_dir[1] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <= fabs( par_dir[2] * half_y ) &&
00584         check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
00585         return true;
00586 
00587     // test two xz plane
00588     if( fabs( par_dir[0] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
00589         fabs( par_dir[2] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
00590         check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
00591         return true;
00592     if( fabs( par_dir[0] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
00593         fabs( par_dir[2] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
00594         check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
00595         return true;
00596 
00597     // test two yz plane
00598     if( fabs( par_dir[1] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
00599         fabs( par_dir[2] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
00600         check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
00601         return true;
00602     if( fabs( par_dir[1] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
00603         fabs( par_dir[2] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
00604         check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
00605         return true;
00606 
00607     return false;
00608 }
00609 
00610 ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
00611 {
00612     ErrorCode rval;
00613     int signs[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
00614                         { -1, -1, 1 },  { 1, -1, 1 },  { 1, 1, 1 },  { -1, 1, 1 } };
00615 
00616     std::vector< EntityHandle > vertices;
00617     for( int i = 0; i < 8; ++i )
00618     {
00619         CartVect coords( center );
00620         for( int j = 0; j < 3; ++j )
00621         {
00622 #if MB_ORIENTED_BOX_UNIT_VECTORS
00623             coords += signs[i][j] * ( axes.col( j ) * length[j] );
00624 #else
00625             coords += signs[i][j] * axes.col( j );
00626 #endif
00627         }
00628         EntityHandle handle;
00629         rval = instance->create_vertex( coords.array(), handle );
00630         if( MB_SUCCESS != rval )
00631         {
00632             instance->delete_entities( &vertices[0], vertices.size() );
00633             return rval;
00634         }
00635         vertices.push_back( handle );
00636     }
00637 
00638     rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
00639     if( MB_SUCCESS != rval )
00640     {
00641         instance->delete_entities( &vertices[0], vertices.size() );
00642         return rval;
00643     }
00644 
00645     return MB_SUCCESS;
00646 }
00647 
00648 void OrientedBox::closest_location_in_box( const CartVect& input_position, CartVect& output_position ) const
00649 {
00650     // get coordinates on box axes
00651     const CartVect from_center = input_position - center;
00652 
00653 #if MB_ORIENTED_BOX_UNIT_VECTORS
00654     CartVect local( from_center % axes.col( 0 ), from_center % axes.col( 1 ), from_center % axes.col( 2 ) );
00655 
00656     for( int i = 0; i < 3; ++i )
00657     {
00658         if( local[i] < -length[i] )
00659             local[i] = -length[i];
00660         else if( local[i] > length[i] )
00661             local[i] = length[i];
00662     }
00663 #else
00664     CartVect local( ( from_center % axes.col( 0 ) ) / ( axes.col( 0 ) % axes.col( 0 ) ),
00665                     ( from_center % axes.col( 1 ) ) / ( axes.col( 1 ) % axes.col( 1 ) ),
00666                     ( from_center % axes.col( 2 ) ) / ( axes.col( 2 ) % axes.col( 2 ) ) );
00667 
00668     for( int i = 0; i < 3; ++i )
00669     {
00670         if( local[i] < -1.0 )
00671             local[i] = -1.0;
00672         else if( local[i] > 1.0 )
00673             local[i] = 1.0;
00674     }
00675 #endif
00676 
00677     output_position = center + local[0] * axes.col( 0 ) + local[1] * axes.col( 1 ) + local[2] * axes.col( 2 );
00678 }
00679 
00680 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines