moab
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 Coroporation, 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 implemenation
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 
00034 #include "moab/Interface.hpp"
00035 #include "moab/CN.hpp"
00036 #include "OrientedBox.hpp"
00037 #include "moab/Range.hpp"
00038 #include "moab/Matrix3.hpp"
00039 #include <ostream>
00040 #include <assert.h>
00041 #include <limits>
00042 
00043 namespace moab {
00044 
00045 #if defined(_MSC_VER) || defined(__MINGW32__)
00046 #  include <float.h>
00047 #  define finite(A) _finite(A)
00048 #endif
00049  
00050 std::ostream& operator<<( std::ostream& s, const OrientedBox& b )
00051 {
00052   return s << b.center 
00053            << " + " 
00054            << b.axis[0] 
00055 #if MB_ORIENTED_BOX_UNIT_VECTORS
00056            << ":" << b.length[0] 
00057 #endif
00058            << " x " 
00059            << b.axis[1] 
00060 #if MB_ORIENTED_BOX_UNIT_VECTORS
00061            << ":" << b.length[1] 
00062 #endif
00063            << " x " 
00064            << b.axis[2]
00065 #if MB_ORIENTED_BOX_UNIT_VECTORS
00066            << ":" << b.length[2] 
00067 #endif
00068             ;
00069 }
00070 
00082 static double point_perp( const CartVect& p,   // closest to this point
00083                           const CartVect& b,   // point on line
00084                           const CartVect& m )  // line direction
00085 {
00086 #if MB_ORIENTED_BOX_UNIT_VECTORS
00087   double t = (m % (p - b));
00088 #else
00089   double t = (m % (p - b)) / (m % m);
00090 #endif
00091   return finite(t) ? t : 0.0;
00092 }
00093 
00094 OrientedBox::OrientedBox( const CartVect axes[3], const CartVect& mid )
00095   : center(mid)
00096 {
00097     // re-order axes by length
00098   CartVect len( axes[0].length(), axes[1].length(), axes[2].length() );
00099   axis[0] = axes[0];
00100   axis[1] = axes[1];
00101   axis[2] = axes[2];
00102   
00103   if (len[2] < len[1])
00104   {
00105     if (len[2] < len[0]) {
00106       std::swap( len[0], len[2] );
00107       std::swap( axis[0], axis[2] );
00108     }
00109   }
00110   else if (len[1] < len[0]) {
00111     std::swap( len[0], len[1] );
00112     std::swap( axis[0], axis[1] );
00113   }
00114   if (len[1] > len[2]) {
00115     std::swap( len[1], len[2] );
00116     std::swap( axis[1], axis[2] );
00117   }
00118   
00119 #if MB_ORIENTED_BOX_UNIT_VECTORS
00120   this->length = len;
00121   if (len[0] > 0.0)
00122     axis[0] /= len[0];
00123   if (len[1] > 0.0)
00124     axis[1] /= len[1];
00125   if (len[2] > 0.0)
00126     axis[2] /= len[2];
00127 #endif
00128 
00129 #if MB_ORIENTED_BOX_OUTER_RADIUS
00130   radius = len.length();
00131 #endif
00132 }
00133 
00134 ErrorCode OrientedBox::tag_handle( Tag& handle_out,
00135                                        Interface* instance,
00136                                        const char* name)
00137 {
00138     // We're going to assume this when mapping the OrientedBox
00139     // to tag data, so assert it.  
00140 #if MB_ORIENTED_BOX_OUTER_RADIUS
00141   const int rad_size = 1;
00142 #else
00143   const int rad_size = 0;
00144 #endif
00145 #if MB_ORIENTED_BOX_UNIT_VECTORS
00146   const int SIZE = rad_size + 15;
00147 #else
00148   const int SIZE = rad_size + 12;
00149 #endif
00150   assert( sizeof(OrientedBox) == SIZE*sizeof(double) );
00151   
00152   return instance->tag_get_handle( name, SIZE, MB_TYPE_DOUBLE,
00153                                    handle_out, MB_TAG_DENSE|MB_TAG_CREAT );
00154 }
00155 
00169 static ErrorCode box_from_axes( OrientedBox& result,
00170                                   Interface* instance,
00171                                   const Range& points )
00172 { 
00173   ErrorCode rval;
00174   
00175     // project points onto axes to get box extents
00176   CartVect min(std::numeric_limits<double>::max()), 
00177              max(-std::numeric_limits<double>::max());
00178   for (Range::iterator i = points.begin(); i != points.end(); ++i)
00179   {
00180     CartVect coords;
00181     rval = instance->get_coords( &*i, 1, coords.array() );
00182     if (MB_SUCCESS != rval)
00183       return rval;
00184     
00185     for (int d = 0; d < 3; ++d)
00186     {
00187       double t = point_perp( coords, result.center, result.axis[d] );
00188       if (t < min[d])
00189         min[d] = t;
00190       if (t > max[d])
00191         max[d] = t;
00192     }
00193   }
00194   
00195     // We now have a box defined by three orthogonal line segments
00196     // that intersect at the center of the box.  Each line segment
00197     // is defined as result.center + t * result.axis[i], where the
00198     // range of t is [min[i], max[i]].
00199   
00200     // Calculate new center
00201   CartVect mid = 0.5 * (min + max);
00202   result.center += mid[0] * result.axis[0] +
00203                    mid[1] * result.axis[1] +
00204                    mid[2] * result.axis[2];
00205   
00206     // reorder axes by length
00207   CartVect range = 0.5 * (max - min);
00208   if (range[2] < range[1])
00209   {
00210     if (range[2] < range[0]) {
00211       std::swap( range[0], range[2] );
00212       std::swap( result.axis[0], result.axis[2] );
00213     }
00214   }
00215   else if (range[1] < range[0]) {
00216     std::swap( range[0], range[1] );
00217     std::swap( result.axis[0], result.axis[1] );
00218   }
00219   if (range[1] > range[2]) {
00220     std::swap( range[1], range[2] );
00221     std::swap( result.axis[1], result.axis[2] );
00222   }
00223 
00224     // scale axis to encompass all points, divide in half
00225 #if MB_ORIENTED_BOX_UNIT_VECTORS
00226   result.length = range;
00227 #else
00228   result.axis[0] *= range[0];
00229   result.axis[1] *= range[1];
00230   result.axis[2] *= range[2];
00231 #endif
00232 
00233 #if MB_ORIENTED_BOX_OUTER_RADIUS
00234   result.radius = range.length();
00235 #endif
00236 
00237   return MB_SUCCESS;
00238 }
00239 
00240 
00241 ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result,
00242                                                   Interface* instance,
00243                                                   const Range& vertices )
00244 {
00245   const Range::iterator begin = vertices.lower_bound( MBVERTEX );
00246   const Range::iterator end = vertices.upper_bound( MBVERTEX );
00247   size_t count = 0;
00248   
00249     // compute mean
00250   CartVect v;
00251   result.center = CartVect( 0, 0, 0 );
00252   for (Range::iterator i = begin; i != end; ++i)
00253   {
00254     ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
00255     if (MB_SUCCESS != rval)
00256       return rval;
00257     result.center += v;
00258     ++count;
00259   }
00260   result.center /= count;
00261   
00262     // compute covariance matrix
00263   Matrix3 a( 0.0 );
00264   for (Range::iterator i = begin; i != end; ++i)
00265   {
00266     ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
00267     if (MB_SUCCESS != rval)
00268       return rval;
00269   
00270     v -= result.center;
00271     a += outer_product( v, v );
00272   }
00273   a /= count;
00274 
00275     // Get axes (Eigenvectors) from covariance matrix
00276   double lambda[3];
00277   moab::Matrix::EigenDecomp( a, lambda, result.axis );
00278   
00279     // Calculate center and extents of box given orientation defined by axes
00280   return box_from_axes( result, instance, vertices );
00281 }
00282 
00283 ErrorCode OrientedBox::covariance_data_from_tris( CovarienceData& result,
00284                                                  Interface* instance,
00285                                                  const Range& elements )
00286 {
00287   ErrorCode rval;
00288   const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
00289   const Range::iterator end = elements.lower_bound( CN::TypeDimensionMap[3].first );
00290   
00291     // compute mean and moments
00292   result.matrix = Matrix3(0.0);
00293   result.center = CartVect(0.0);
00294   result.area = 0.0;
00295   for (Range::iterator i = begin; i != end; ++i)
00296   {
00297     const EntityHandle* conn;
00298     int conn_len;
00299     rval = instance->get_connectivity( *i, conn, conn_len );
00300     if (MB_SUCCESS != rval)
00301       return rval;
00302     
00303       // for each triangle in the 2-D cell
00304     for (int j = 2; j < conn_len; ++j)
00305     {
00306       EntityHandle vertices[3] = { conn[0], conn[j-1], conn[j] };
00307       CartVect coords[3];
00308       rval = instance->get_coords( vertices, 3, coords[0].array() );
00309       if (MB_SUCCESS != rval)
00310         return rval;
00311       
00312         // edge vectors
00313       const CartVect edge0 = coords[1] - coords[0];
00314       const CartVect edge1 = coords[2] - coords[0];
00315       const CartVect centroid = (coords[0] + coords[1] + coords[2]) / 3;
00316       const double tri_area2 = (edge0 * edge1).length();
00317       result.area += tri_area2;
00318       result.center += tri_area2 * centroid;
00319       
00320       result.matrix += tri_area2 * (9 * outer_product( centroid,  centroid  ) +
00321                                     outer_product( coords[0], coords[0] ) +
00322                                     outer_product( coords[1], coords[1] ) +
00323                                     outer_product( coords[2], coords[2] ));
00324     } // for each triangle
00325   } // for each element
00326 
00327   return MB_SUCCESS;
00328 }
00329 
00330 
00331 ErrorCode OrientedBox::compute_from_2d_cells( OrientedBox& result,
00332                                                   Interface* instance,
00333                                                   const Range& elements )
00334 {
00335     // Get orientation data from elements
00336   CovarienceData data;
00337   ErrorCode rval = covariance_data_from_tris( data, instance, elements );
00338   if (MB_SUCCESS != rval)
00339     return rval;
00340   
00341     // get vertices from elements
00342   Range points;
00343   rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
00344   if (MB_SUCCESS != rval)
00345     return rval;
00346     
00347     // Calculate box given points and orientation data
00348   return compute_from_covariance_data( result, instance, data, points );
00349 }
00350 
00351 ErrorCode OrientedBox::compute_from_covariance_data(
00352                                                 OrientedBox& result,
00353                                                 Interface* instance,
00354                                                 CovarienceData& data,
00355                                                 const Range& vertices )
00356 {
00357   if (data.area <= 0.0) {
00358     CartVect axis[3] = { CartVect(0.), CartVect(0.), CartVect(0.) };
00359     result = OrientedBox( axis, CartVect(0.) );
00360     return MB_SUCCESS;
00361   }
00362 
00363     // get center from sum
00364   result.center = data.center / data.area;
00365 
00366     // get covariance matrix from moments
00367   data.matrix /= 12 * data.area;
00368   data.matrix -= outer_product( result.center, result.center );
00369 
00370     // get axes (Eigenvectors) from covariance matrix
00371   double lamda[3];
00372   moab::Matrix::EigenDecomp( data.matrix, lamda, result.axis );
00373 
00374     // We now have only the axes.  Calculate proper center
00375     // and extents for enclosed points.
00376   return box_from_axes( result, instance, vertices );
00377 }      
00378 
00379 bool OrientedBox::contained( const CartVect& point, double tol ) const
00380 {
00381   CartVect from_center = point - center;
00382 #if MB_ORIENTED_BOX_UNIT_VECTORS
00383   return fabs(from_center % axis[0]) - length[0] <= tol &&
00384          fabs(from_center % axis[1]) - length[1] <= tol &&
00385          fabs(from_center % axis[2]) - length[2] <= tol ;
00386 #else
00387   for (int i = 0; i < 3; ++i) {
00388     double length = axis[i].length();
00389     if (fabs(from_center % axis[i]) - length*length > length*tol)
00390       return false;
00391   }
00392   return true;
00393 #endif
00394 }
00395 
00396 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
00397                                                 Interface* moab_instance,
00398                                                 const CovarienceData* data,
00399                                                 unsigned data_length,
00400                                                 const Range& vertices )
00401 {
00402     // Sum input CovarienceData structures
00403   CovarienceData data_sum( Matrix3(0.0), CartVect(0.0), 0.0 );
00404   for (const CovarienceData* const end = data+data_length; data != end; ++data) {
00405     data_sum.matrix += data->matrix;
00406     data_sum.center += data->center;
00407     data_sum.area += data->area;
00408   }
00409     // Compute box from sum of structs
00410   return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
00411 }
00412 
00413 
00414 
00415 //bool OrientedBox::contained( const OrientedBox& box, double tol ) const
00416 //{
00417 //  for (int i = -1; i < 2; i += 2) 
00418 //  {
00419 //    for (int j = -1; j < 2; j += 2) 
00420 //    {
00421 //      for (int k = -1; k < 2; k += 2) 
00422 //      {
00423 //        CartVect corner( center );
00424 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
00425 //        corner += i * box.length[0] * box.axis[0];
00426 //        corner += j * box.length[1] * box.axis[1];
00427 //        corner += k * box.length[2] * box.axis[2];
00428 //#else
00429 //        corner += i * box.axis[0];
00430 //        corner += j * box.axis[1];
00431 //        corner += k * box.axis[2];
00432 //#endif
00433 //        if (!contained( corner, tol ))
00434 //          return false;
00435 //      }
00436 //    }
00437 //  }
00438 //  return true;
00439 //}
00440 
00441 
00442 /* This is a helper function to check limits on ray length, turning the box-ray 
00443  * intersection test into a box-segment intersection test. Use this to test the
00444  * limits against one side (plane) of the box. The side of the box (plane) is
00445  * normal to an axis.
00446  *
00447  *   normal_par_pos  Coordinate of particle's position along axis normal to side of box
00448  *   normal_par_dir  Coordinate of particle's direction along axis normal to side of box
00449  *   half_extent     Distance between center of box and side of box
00450  *   nonneg_ray_len  Maximum ray length in positive direction (in front of origin)
00451  *   neg_ray_len     Maximum ray length in negative direction (behind origin)
00452  *   return          true if intersection with plane occurs within distance limits
00453  *
00454  * ray equation:   intersection = origin + dist*direction
00455  * plane equation: intersection.plane_normal = half_extent
00456  * 
00457  * Assume plane_normal and direction are unit vectors. Combine equations.
00458  *
00459  *     (origin + dist*direction).plane_normal = half_extent
00460  *     origin.plane_normal + dist*direction.plane_normal = half_extent
00461  *     dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
00462  *
00463  * Although this solves for distance, avoid floating point division.
00464  *
00465  *     dist*direction.plane_normal = half_extent - origin.plane_normal
00466  *
00467  * Use inequalities to test dist against ray length limits. Be aware that 
00468  * inequalities change due to sign of direction.plane_normal.
00469  */ 
00470 inline bool check_ray_limits(const double  normal_par_pos,
00471                              const double  normal_par_dir,
00472                              const double  half_extent,
00473                              const double* nonneg_ray_len,
00474                              const double* neg_ray_len ) {
00475 
00476   const double extent_pos_diff = half_extent - normal_par_pos;
00477 
00478   // limit in positive direction
00479   if(nonneg_ray_len) { // should be 0 <= t <= nonneg_ray_len
00480     assert(0 <= *nonneg_ray_len);
00481     if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
00482       if(*nonneg_ray_len*normal_par_dir>=extent_pos_diff && extent_pos_diff>=0) return true;
00483     } else if(normal_par_dir<0) {
00484       if(*nonneg_ray_len*normal_par_dir<=extent_pos_diff && extent_pos_diff<=0) return true;
00485     }
00486   } else {            // should be 0 <= t
00487     if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
00488       if(extent_pos_diff>=0) return true;
00489     } else if(normal_par_dir<0) {
00490       if(extent_pos_diff<=0) return true;
00491     }
00492   }
00493 
00494   // limit in negative direction
00495   if(neg_ray_len) {   // should be neg_ray_len <= t < 0
00496     assert(0 >= *neg_ray_len);
00497     if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
00498       if(*neg_ray_len*normal_par_dir<=extent_pos_diff && extent_pos_diff<0) return true;
00499     } else if(normal_par_dir<0) {
00500       if(*neg_ray_len*normal_par_dir>=extent_pos_diff && extent_pos_diff>0) return true;
00501     }
00502   }
00503 
00504   return false;
00505 }
00506 
00507 /* This implementation copied from cgmMC (overlap.C).
00508  * Original author:  Tim Tautges?
00509  */
00510 bool OrientedBox::intersect_ray( const CartVect& ray_origin,
00511                                  const CartVect& ray_direction,
00512                  const double    reps,
00513                  const double*   nonneg_ray_len,
00514                                  const double*   neg_ray_len ) const
00515 {
00516   // test distance from box center to line
00517   const CartVect cx       = center - ray_origin;
00518   const double dist_s     = cx % ray_direction;
00519   const double dist_sq    = cx % cx - (dist_s*dist_s);
00520   const double max_diagsq = outer_radius_squared(reps);
00521   
00522   // For the largest sphere, no intersections exist if discriminant is negative.
00523   // Geometrically, if distance from box center to line is greater than the 
00524   // longest diagonal, there is no intersection.
00525   // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
00526   if(dist_sq > max_diagsq) return false;
00527 
00528   // If the closest possible intersection must be closer than nonneg_ray_len. Be 
00529   // careful with absolute value, squaring distances, and subtracting squared
00530   // distances.
00531   if (nonneg_ray_len) {
00532     assert(0<=*nonneg_ray_len);
00533     double max_len;
00534     if(neg_ray_len) {
00535       assert(0>=*neg_ray_len);
00536       max_len = std::max(*nonneg_ray_len,-(*neg_ray_len));
00537     } else {
00538       max_len = *nonneg_ray_len;
00539     }
00540     const double temp = fabs(dist_s) - max_len;
00541     if(0.0<temp && temp*temp>max_diagsq) return false;
00542   } 
00543 
00544   // if smaller than shortest diagonal, we do hit
00545   if (dist_sq < inner_radius_squared(reps)) {
00546     // nonnegative direction
00547     if(dist_s>=0.0 ) {
00548       if(nonneg_ray_len) {
00549         if(*nonneg_ray_len>dist_s) return true;
00550       } else {
00551         return true;
00552       }
00553     // negative direction 
00554     } else {
00555       if(neg_ray_len && *neg_ray_len<dist_s) return true;
00556     }
00557   }
00558   
00559     // get transpose of axes
00560     // Note: if axes were stored as a matrix, could skip
00561     // transpose and just switch order of operands in
00562     // matrix-vector multiplies below. - J.K.
00563   //Matrix3 B( axis[0][0], axis[1][0], axis[2][0],
00564   //             axis[0][1], axis[1][1], axis[2][1],
00565   //             axis[0][2], axis[1][2], axis[2][2] );
00566   Matrix3 B( axis[0][0], axis[0][1], axis[0][2],
00567                axis[1][0], axis[1][1], axis[1][2],
00568                axis[2][0], axis[2][1], axis[2][2] );
00569   //CartVect T = B * -center;
00570   
00571     // transform ray to box coordintae system
00572   //CartVect par_pos = T + B * b;
00573   CartVect par_pos = B * (ray_origin - center);
00574   CartVect par_dir = B * ray_direction;
00575 
00576   // Fast Rejection Test: Ray will not intersect if it is going away from the box.
00577   // This will not work for rays with neg_ray_len. length[0] is half of box width 
00578   // along axis[0].
00579   const double half_x = length[0] + reps;
00580   const double half_y = length[1] + reps;
00581   const double half_z = length[2] + reps;
00582   if(!neg_ray_len) {
00583     if ((par_pos[0] >  half_x && par_dir[0] >= 0) ||
00584     (par_pos[0] < -half_x && par_dir[0] <= 0))
00585       return false;
00586   
00587     if ((par_pos[1] >  half_y && par_dir[1] >= 0) ||
00588     (par_pos[1] < -half_y && par_dir[1] <= 0))
00589       return false;
00590     
00591     if ((par_pos[2] >  half_z && par_dir[2] >= 0) ||
00592     (par_pos[2] < -half_z && par_dir[2] <= 0))
00593       return false;
00594   }
00595 
00596   // test if ray_origin is inside box
00597   if (par_pos[0] <= half_x && par_pos[0] >= -half_x &&
00598       par_pos[1] <= half_y && par_pos[1] >= -half_y &&
00599       par_pos[2] <= half_z && par_pos[2] >= -half_z)
00600     return true;
00601 
00602     //test two xy plane
00603   if (fabs(par_dir[0] * (half_z - par_pos[2]) + par_dir[2] * par_pos[0]) 
00604         <= fabs(par_dir[2] * half_x) &&  // test against x extents using z
00605       fabs(par_dir[1] * (half_z - par_pos[2]) + par_dir[2] * par_pos[1]) 
00606         <= fabs(par_dir[2] * half_y) &&  // test against y extents using z
00607       check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
00608     return true;
00609   if (fabs(par_dir[0] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[0]) 
00610         <= fabs(par_dir[2] * half_x) && 
00611       fabs(par_dir[1] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[1]) 
00612         <= fabs(par_dir[2] * half_y) &&
00613       check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
00614     return true;
00615 
00616     //test two xz plane
00617   if (fabs(par_dir[0] * (half_y - par_pos[1]) + par_dir[1] * par_pos[0]) 
00618         <= fabs(par_dir[1] * half_x) && 
00619       fabs(par_dir[2] * (half_y - par_pos[1]) + par_dir[1] * par_pos[2]) 
00620         <= fabs(par_dir[1] * half_z) &&
00621       check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
00622     return true;
00623   if (fabs(par_dir[0] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[0]) 
00624         <= fabs(par_dir[1] * half_x) && 
00625       fabs(par_dir[2] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[2])
00626         <= fabs(par_dir[1] * half_z) &&
00627       check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
00628     return true;
00629 
00630     //test two yz plane
00631   if (fabs(par_dir[1] * (half_x - par_pos[0]) + par_dir[0] * par_pos[1]) 
00632         <= fabs(par_dir[0] * half_y) &&
00633       fabs(par_dir[2] * (half_x - par_pos[0]) + par_dir[0] * par_pos[2]) 
00634         <= fabs(par_dir[0] * half_z) &&
00635       check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
00636     return true;
00637   if (fabs(par_dir[1] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[1])
00638         <= fabs(par_dir[0] * half_y) &&
00639       fabs(par_dir[2] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[2]) 
00640         <= fabs(par_dir[0] * half_z) &&
00641       check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
00642     return true;
00643 
00644   return false;
00645 }
00646 
00647 ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
00648 {
00649   ErrorCode rval;
00650   int signs[8][3] = { { -1, -1, -1 },
00651                       {  1, -1, -1 },
00652                       {  1,  1, -1 },
00653                       { -1,  1, -1 },
00654                       { -1, -1,  1 },
00655                       {  1, -1,  1 },
00656                       {  1,  1,  1 },
00657                       { -1,  1,  1 } };
00658                       
00659   std::vector<EntityHandle> vertices;
00660   for (int i = 0; i < 8; ++i)
00661   {
00662     CartVect coords(center);
00663     for (int j = 0; j < 3; ++j){
00664 #if MB_ORIENTED_BOX_UNIT_VECTORS
00665       coords += signs[i][j] * (axis[j]*length[j]);
00666 #else
00667       coords += signs[i][j] * axis[j];
00668 #endif
00669     }
00670     EntityHandle handle;
00671     rval = instance->create_vertex( coords.array(), handle );
00672     if (MB_SUCCESS != rval) {
00673       instance->delete_entities( &vertices[0], vertices.size() );
00674       return rval;
00675     }
00676     vertices.push_back( handle );
00677   }
00678   
00679   rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
00680   if (MB_SUCCESS != rval) {
00681     instance->delete_entities( &vertices[0], vertices.size() );
00682     return rval;
00683   }
00684   
00685   return MB_SUCCESS;
00686 }
00687   
00688 void OrientedBox::closest_location_in_box( 
00689                                     const CartVect& input_position,
00690                                     CartVect& output_position ) const
00691 {
00692     // get coordinates on box axes
00693   const CartVect from_center = input_position - center;
00694 
00695 #if MB_ORIENTED_BOX_UNIT_VECTORS
00696   CartVect local( from_center % axis[0],
00697                     from_center % axis[1],
00698                     from_center % axis[2] );
00699 
00700   for (int i = 0; i < 3; ++i) {
00701     if (local[i] < -length[i])
00702       local[i] = -length[i];
00703     else if (local[i] > length[i])
00704       local[i] =  length[i];
00705   }
00706 #else
00707   CartVect local( (from_center % axis[0]) / (axis[0] % axis[0]),
00708                     (from_center % axis[1]) / (axis[1] % axis[1]),
00709                     (from_center % axis[2]) / (axis[2] % axis[2]) );
00710 
00711   for (int i = 0; i < 3; ++i) {
00712     if (local[i] < -1.0)
00713       local[i] = -1.0;
00714     else if (local[i] > 1.0)
00715       local[i] = 1.0;
00716   }
00717 #endif
00718 
00719   output_position = center
00720                   + local[0] * axis[0] 
00721                   + local[1] * axis[1]
00722                   + local[2] * axis[2];
00723 }
00724   
00725 } // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines