MOAB: Mesh Oriented datABase  (version 5.4.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 ([email protected])
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 <cassert>
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,
00328                                                      Interface* instance,
00329                                                      CovarienceData& data,
00330                                                      const Range& vertices )
00331 {
00332     if( data.area <= 0.0 )
00333     {
00334         Matrix3 empty_axes( 0.0 );
00335         result = OrientedBox( empty_axes, CartVect( 0. ) );
00336         return MB_SUCCESS;
00337     }
00338 
00339     // get center from sum
00340     result.center = data.center / data.area;
00341 
00342     // get covariance matrix from moments
00343     data.matrix /= 12 * data.area;
00344     data.matrix -= outer_product( result.center, result.center );
00345 
00346     // get axes (Eigenvectors) from covariance matrix
00347     CartVect lambda;
00348     data.matrix.eigen_decomposition( lambda, result.axes );
00349 
00350     // We now have only the axes.  Calculate proper center
00351     // and extents for enclosed points.
00352     return box_from_axes( result, instance, vertices );
00353 }
00354 
00355 bool OrientedBox::contained( const CartVect& point, double tol ) const
00356 {
00357     CartVect from_center = point - center;
00358 #if MB_ORIENTED_BOX_UNIT_VECTORS
00359     return fabs( from_center % axes.col( 0 ) ) - length[0] <= tol &&
00360            fabs( from_center % axes.col( 1 ) ) - length[1] <= tol &&
00361            fabs( from_center % axes.col( 2 ) ) - length[2] <= tol;
00362 #else
00363     for( int i = 0; i < 3; ++i )
00364     {
00365         double length = axes.col( i ).length();
00366         if( fabs( from_center % axes.col( i ) ) - length * length > length * tol ) return false;
00367     }
00368     return true;
00369 #endif
00370 }
00371 
00372 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
00373                                                      Interface* moab_instance,
00374                                                      const CovarienceData* data,
00375                                                      unsigned data_length,
00376                                                      const Range& vertices )
00377 {
00378     // Sum input CovarienceData structures
00379     CovarienceData data_sum( Matrix3( 0.0 ), CartVect( 0.0 ), 0.0 );
00380     for( const CovarienceData* const end = data + data_length; data != end; ++data )
00381     {
00382         data_sum.matrix += data->matrix;
00383         data_sum.center += data->center;
00384         data_sum.area += data->area;
00385     }
00386     // Compute box from sum of structs
00387     return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
00388 }
00389 
00390 // bool OrientedBox::contained( const OrientedBox& box, double tol ) const
00391 //{
00392 //  for (int i = -1; i < 2; i += 2)
00393 //  {
00394 //    for (int j = -1; j < 2; j += 2)
00395 //    {
00396 //      for (int k = -1; k < 2; k += 2)
00397 //      {
00398 //        CartVect corner( center );
00399 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
00400 //        corner += i * box.length[0] * box.axis.col(0);
00401 //        corner += j * box.length[1] * box.axis.col(1);
00402 //        corner += k * box.length[2] * box.axis.col(2);
00403 //#else
00404 //        corner += i * box.axis.col(0);
00405 //        corner += j * box.axis.col(1);
00406 //        corner += k * box.axis.col(2);
00407 //#endif
00408 //        if (!contained( corner, tol ))
00409 //          return false;
00410 //      }
00411 //    }
00412 //  }
00413 //  return true;
00414 //}
00415 
00416 /* This is a helper function to check limits on ray length, turning the box-ray
00417  * intersection test into a box-segment intersection test. Use this to test the
00418  * limits against one side (plane) of the box. The side of the box (plane) is
00419  * normal to an axis.
00420  *
00421  *   normal_par_pos  Coordinate of particle's position along axis normal to side of box
00422  *   normal_par_dir  Coordinate of particle's direction along axis normal to side of box
00423  *   half_extent     Distance between center of box and side of box
00424  *   nonneg_ray_len  Maximum ray length in positive direction (in front of origin)
00425  *   neg_ray_len     Maximum ray length in negative direction (behind origin)
00426  *   return          true if intersection with plane occurs within distance limits
00427  *
00428  * ray equation:   intersection = origin + dist*direction
00429  * plane equation: intersection.plane_normal = half_extent
00430  *
00431  * Assume plane_normal and direction are unit vectors. Combine equations.
00432  *
00433  *     (origin + dist*direction).plane_normal = half_extent
00434  *     origin.plane_normal + dist*direction.plane_normal = half_extent
00435  *     dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
00436  *
00437  * Although this solves for distance, avoid floating point division.
00438  *
00439  *     dist*direction.plane_normal = half_extent - origin.plane_normal
00440  *
00441  * Use inequalities to test dist against ray length limits. Be aware that
00442  * inequalities change due to sign of direction.plane_normal.
00443  */
00444 inline bool check_ray_limits( const double normal_par_pos,
00445                               const double normal_par_dir,
00446                               const double half_extent,
00447                               const double* nonneg_ray_len,
00448                               const double* neg_ray_len )
00449 {
00450 
00451     const double extent_pos_diff = half_extent - normal_par_pos;
00452 
00453     // limit in positive direction
00454     if( nonneg_ray_len )
00455     {  // should be 0 <= t <= nonneg_ray_len
00456         assert( 0 <= *nonneg_ray_len );
00457         if( normal_par_dir > 0 )
00458         {  // if/else if needed for pos/neg divisor
00459             if( *nonneg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff >= 0 ) return true;
00460         }
00461         else if( normal_par_dir < 0 )
00462         {
00463             if( *nonneg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff <= 0 ) return true;
00464         }
00465     }
00466     else
00467     {  // should be 0 <= t
00468         if( normal_par_dir > 0 )
00469         {  // if/else if needed for pos/neg divisor
00470             if( extent_pos_diff >= 0 ) return true;
00471         }
00472         else if( normal_par_dir < 0 )
00473         {
00474             if( extent_pos_diff <= 0 ) return true;
00475         }
00476     }
00477 
00478     // limit in negative direction
00479     if( neg_ray_len )
00480     {  // should be neg_ray_len <= t < 0
00481         assert( 0 >= *neg_ray_len );
00482         if( normal_par_dir > 0 )
00483         {  // if/else if needed for pos/neg divisor
00484             if( *neg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff < 0 ) return true;
00485         }
00486         else if( normal_par_dir < 0 )
00487         {
00488             if( *neg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff > 0 ) return true;
00489         }
00490     }
00491 
00492     return false;
00493 }
00494 
00495 /* This implementation copied from cgmMC (overlap.C).
00496  * Original author:  Tim Tautges?
00497  */
00498 bool OrientedBox::intersect_ray( const CartVect& ray_origin,
00499                                  const CartVect& ray_direction,
00500                                  const double reps,
00501                                  const double* nonneg_ray_len,
00502                                  const double* neg_ray_len ) const
00503 {
00504     // test distance from box center to line
00505     const CartVect cx       = center - ray_origin;
00506     const double dist_s     = cx % ray_direction;
00507     const double dist_sq    = cx % cx - ( dist_s * dist_s );
00508     const double max_diagsq = outer_radius_squared( reps );
00509 
00510     // For the largest sphere, no intersections exist if discriminant is negative.
00511     // Geometrically, if distance from box center to line is greater than the
00512     // longest diagonal, there is no intersection.
00513     // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
00514     if( dist_sq > max_diagsq ) return false;
00515 
00516     // If the closest possible intersection must be closer than nonneg_ray_len. Be
00517     // careful with absolute value, squaring distances, and subtracting squared
00518     // distances.
00519     if( nonneg_ray_len )
00520     {
00521         assert( 0 <= *nonneg_ray_len );
00522         double max_len;
00523         if( neg_ray_len )
00524         {
00525             assert( 0 >= *neg_ray_len );
00526             max_len = std::max( *nonneg_ray_len, -( *neg_ray_len ) );
00527         }
00528         else
00529         {
00530             max_len = *nonneg_ray_len;
00531         }
00532         const double temp = fabs( dist_s ) - max_len;
00533         if( 0.0 < temp && temp * temp > max_diagsq ) return false;
00534     }
00535 
00536     // if smaller than shortest diagonal, we do hit
00537     if( dist_sq < inner_radius_squared( reps ) )
00538     {
00539         // nonnegative direction
00540         if( dist_s >= 0.0 )
00541         {
00542             if( nonneg_ray_len )
00543             {
00544                 if( *nonneg_ray_len > dist_s ) return true;
00545             }
00546             else
00547             {
00548                 return true;
00549             }
00550             // negative direction
00551         }
00552         else
00553         {
00554             if( neg_ray_len && *neg_ray_len < dist_s ) return true;
00555         }
00556     }
00557 
00558     // get transpose of axes
00559     Matrix3 B = axes.transpose();
00560 
00561     // transform ray to box coordintae system
00562     CartVect par_pos = B * ( ray_origin - center );
00563     CartVect par_dir = B * ray_direction;
00564 
00565     // Fast Rejection Test: Ray will not intersect if it is going away from the box.
00566     // This will not work for rays with neg_ray_len. length[0] is half of box width
00567     // along axes.col(0).
00568     const double half_x = length[0] + reps;
00569     const double half_y = length[1] + reps;
00570     const double half_z = length[2] + reps;
00571     if( !neg_ray_len )
00572     {
00573         if( ( par_pos[0] > half_x && par_dir[0] >= 0 ) || ( par_pos[0] < -half_x && par_dir[0] <= 0 ) ) return false;
00574 
00575         if( ( par_pos[1] > half_y && par_dir[1] >= 0 ) || ( par_pos[1] < -half_y && par_dir[1] <= 0 ) ) return false;
00576 
00577         if( ( par_pos[2] > half_z && par_dir[2] >= 0 ) || ( par_pos[2] < -half_z && par_dir[2] <= 0 ) ) return false;
00578     }
00579 
00580     // test if ray_origin is inside box
00581     if( par_pos[0] <= half_x && par_pos[0] >= -half_x && par_pos[1] <= half_y && par_pos[1] >= -half_y &&
00582         par_pos[2] <= half_z && par_pos[2] >= -half_z )
00583         return true;
00584 
00585     // test two xy plane
00586     if( fabs( par_dir[0] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <=
00587             fabs( par_dir[2] * half_x ) &&  // test against x extents using z
00588         fabs( par_dir[1] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <=
00589             fabs( par_dir[2] * half_y ) &&  // test against y extents using z
00590         check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
00591         return true;
00592     if( fabs( par_dir[0] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <= fabs( par_dir[2] * half_x ) &&
00593         fabs( par_dir[1] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <= fabs( par_dir[2] * half_y ) &&
00594         check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
00595         return true;
00596 
00597     // test two xz plane
00598     if( fabs( par_dir[0] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
00599         fabs( par_dir[2] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
00600         check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
00601         return true;
00602     if( fabs( par_dir[0] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
00603         fabs( par_dir[2] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
00604         check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
00605         return true;
00606 
00607     // test two yz plane
00608     if( fabs( par_dir[1] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
00609         fabs( par_dir[2] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
00610         check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
00611         return true;
00612     if( fabs( par_dir[1] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
00613         fabs( par_dir[2] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
00614         check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
00615         return true;
00616 
00617     return false;
00618 }
00619 
00620 ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
00621 {
00622     ErrorCode rval;
00623     int signs[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
00624                         { -1, -1, 1 },  { 1, -1, 1 },  { 1, 1, 1 },  { -1, 1, 1 } };
00625 
00626     std::vector< EntityHandle > vertices;
00627     for( int i = 0; i < 8; ++i )
00628     {
00629         CartVect coords( center );
00630         for( int j = 0; j < 3; ++j )
00631         {
00632 #if MB_ORIENTED_BOX_UNIT_VECTORS
00633             coords += signs[i][j] * ( axes.col( j ) * length[j] );
00634 #else
00635             coords += signs[i][j] * axes.col( j );
00636 #endif
00637         }
00638         EntityHandle handle;
00639         rval = instance->create_vertex( coords.array(), handle );
00640         if( MB_SUCCESS != rval )
00641         {
00642             instance->delete_entities( &vertices[0], vertices.size() );
00643             return rval;
00644         }
00645         vertices.push_back( handle );
00646     }
00647 
00648     rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
00649     if( MB_SUCCESS != rval )
00650     {
00651         instance->delete_entities( &vertices[0], vertices.size() );
00652         return rval;
00653     }
00654 
00655     return MB_SUCCESS;
00656 }
00657 
00658 void OrientedBox::closest_location_in_box( const CartVect& input_position, CartVect& output_position ) const
00659 {
00660     // get coordinates on box axes
00661     const CartVect from_center = input_position - center;
00662 
00663 #if MB_ORIENTED_BOX_UNIT_VECTORS
00664     CartVect local( from_center % axes.col( 0 ), from_center % axes.col( 1 ), from_center % axes.col( 2 ) );
00665 
00666     for( int i = 0; i < 3; ++i )
00667     {
00668         if( local[i] < -length[i] )
00669             local[i] = -length[i];
00670         else if( local[i] > length[i] )
00671             local[i] = length[i];
00672     }
00673 #else
00674     CartVect local( ( from_center % axes.col( 0 ) ) / ( axes.col( 0 ) % axes.col( 0 ) ),
00675                     ( from_center % axes.col( 1 ) ) / ( axes.col( 1 ) % axes.col( 1 ) ),
00676                     ( from_center % axes.col( 2 ) ) / ( axes.col( 2 ) % axes.col( 2 ) ) );
00677 
00678     for( int i = 0; i < 3; ++i )
00679     {
00680         if( local[i] < -1.0 )
00681             local[i] = -1.0;
00682         else if( local[i] > 1.0 )
00683             local[i] = 1.0;
00684     }
00685 #endif
00686 
00687     output_position = center + local[0] * axes.col( 0 ) + local[1] * axes.col( 1 ) + local[2] * axes.col( 2 );
00688 }
00689 
00690 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines