moab
|
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