Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
OrientedBox.cpp
Go to the documentation of this file.
1 /*
2  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
3  * storing and accessing finite element mesh data.
4  *
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  *
9  * This library is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 2.1 of the License, or (at your option) any later version.
13  *
14  */
15 
16 /*
17  * The algorithms for the calculation of the oriented box from a
18  * set of points or a set of cells was copied from the implementation
19  " in the "Visualization Toolkit". J.K. - 2006-07-19
20  *
21  * Program: Visualization Toolkit
22  * Module: $RCSfile$
23  *
24  * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
25  * All rights reserved.
26  * See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
27  */
28 
29 /**\file OrientedBox.cpp
30  *\author Jason Kraftcheck ([email protected])
31  *\date 2006-07-18
32  */
33 
34 #include "moab/Interface.hpp"
35 #include "moab/CN.hpp"
36 #include "moab/OrientedBox.hpp"
37 #include "moab/Range.hpp"
38 #include <ostream>
39 #include <cassert>
40 #include <limits>
41 
42 namespace moab
43 {
44 
45 std::ostream& operator<<( std::ostream& s, const OrientedBox& b )
46 {
47  return s << b.center << " + " << b.axes.col( 0 )
48 #if MB_ORIENTED_BOX_UNIT_VECTORS
49  << ":" << b.length[0]
50 #endif
51  << " x " << b.axes.col( 1 )
52 #if MB_ORIENTED_BOX_UNIT_VECTORS
53  << ":" << b.length[1]
54 #endif
55  << " x " << b.axes.col( 2 )
56 #if MB_ORIENTED_BOX_UNIT_VECTORS
57  << ":" << b.length[2]
58 #endif
59  ;
60 }
61 
62 /**\brief Find closest point on line
63  *
64  * Find the point on the line for which a line trough the
65  * input point \a p and the result position is orthogonal to
66  * the input line.
67  * \param p The point for which to find the perpendicular
68  * \param b A point on the line
69  * \param m The direction of the line
70  * \return The location on the line specified as 't' in the
71  * formula t * m + b
72  */
73 static double point_perp( const CartVect& p, // closest to this point
74  const CartVect& b, // point on line
75  const CartVect& m ) // line direction
76 {
77 #if MB_ORIENTED_BOX_UNIT_VECTORS
78  double t = ( m % ( p - b ) );
79 #else
80  double t = ( m % ( p - b ) ) / ( m % m );
81 #endif
82  return Util::is_finite( t ) ? t : 0.0;
83 }
84 
85 void OrientedBox::order_axes_by_length( double ax1_len, double ax2_len, double ax3_len )
86 {
87 
88  CartVect len( ax1_len, ax2_len, ax3_len );
89 
90  if( len[2] < len[1] )
91  {
92  if( len[2] < len[0] )
93  {
94  std::swap( len[0], len[2] );
95  axes.swapcol( 0, 2 );
96  }
97  }
98  else if( len[1] < len[0] )
99  {
100  std::swap( len[0], len[1] );
101  axes.swapcol( 0, 1 );
102  }
103  if( len[1] > len[2] )
104  {
105  std::swap( len[1], len[2] );
106  axes.swapcol( 1, 2 );
107  }
108 
109 #if MB_ORIENTED_BOX_UNIT_VECTORS
110  length = len;
111  if( len[0] > 0.0 ) axes.colscale( 0, 1.0 / len[0] );
112  if( len[1] > 0.0 ) axes.colscale( 1, 1.0 / len[1] );
113  if( len[2] > 0.0 ) axes.colscale( 2, 1.0 / len[2] );
114 #endif
115 
116 #if MB_ORIENTED_BOX_OUTER_RADIUS
117  radius = len.length();
118 #endif
119 }
120 
121 OrientedBox::OrientedBox( const CartVect axes_in[3], const CartVect& mid ) : center( mid )
122 {
123 
124  axes = Matrix3( axes_in[0], axes_in[1], axes_in[2], false );
125 
126  order_axes_by_length( axes_in[0].length(), axes_in[1].length(), axes_in[2].length() );
127 }
128 
129 OrientedBox::OrientedBox( const Matrix3& axes_mat, const CartVect& mid ) : center( mid ), axes( axes_mat )
130 {
131  order_axes_by_length( axes.col( 0 ).length(), axes.col( 1 ).length(), axes.col( 2 ).length() );
132 }
133 
134 ErrorCode OrientedBox::tag_handle( Tag& handle_out, Interface* instance, const char* name )
135 {
136  // We're going to assume this when mapping the OrientedBox
137  // to tag data, so assert it.
138 #if MB_ORIENTED_BOX_OUTER_RADIUS
139  const int rad_size = 1;
140 #else
141  const int rad_size = 0;
142 #endif
143 #if MB_ORIENTED_BOX_UNIT_VECTORS
144  const int SIZE = rad_size + 15;
145 #else
146  const int SIZE = rad_size + 12;
147 #endif
148  assert( sizeof( OrientedBox ) == SIZE * sizeof( double ) );
149 
150  return instance->tag_get_handle( name, SIZE, MB_TYPE_DOUBLE, handle_out, MB_TAG_DENSE | MB_TAG_CREAT );
151 }
152 
153 /**\brief Common code for box calculation
154  *
155  * Given the orientation of the box and an approximate center,
156  * calculate the exact center and extents of the box.
157  *
158  *\param result.center As input, the approximate center of the box.
159  * As output, the exact center of the box.
160  *\param result.axes As input, directions of principal axes corresponding
161  * to the orientation of the box. Axes are assumed to
162  * be unit-length on input. Output will include extents
163  * of box.
164  *\param points The set of points the box should contain.
165  */
166 static ErrorCode box_from_axes( OrientedBox& result, Interface* instance, const Range& points )
167 {
168  ErrorCode rval;
169 
170  // project points onto axes to get box extents
171  CartVect min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() );
172  for( Range::iterator i = points.begin(); i != points.end(); ++i )
173  {
174  CartVect coords;
175  rval = instance->get_coords( &*i, 1, coords.array() );MB_CHK_ERR( rval );
176 
177  for( int d = 0; d < 3; ++d )
178  {
179  const double t = point_perp( coords, result.center, result.axes.col( d ) );
180  if( t < min[d] ) min[d] = t;
181  if( t > max[d] ) max[d] = t;
182  }
183  }
184 
185  // We now have a box defined by three orthogonal line segments
186  // that intersect at the center of the box. Each line segment
187  // is defined as result.center + t * result.axes[i], where the
188  // range of t is [min[i], max[i]].
189 
190  // Calculate new center
191  const CartVect mid = 0.5 * ( min + max );
192  result.center += mid[0] * result.axes.col( 0 ) + mid[1] * result.axes.col( 1 ) + mid[2] * result.axes.col( 2 );
193 
194  // reorder axes by length
195  CartVect range = 0.5 * ( max - min );
196  if( range[2] < range[1] )
197  {
198  if( range[2] < range[0] )
199  {
200  std::swap( range[0], range[2] );
201  result.axes.swapcol( 0, 2 );
202  }
203  }
204  else if( range[1] < range[0] )
205  {
206  std::swap( range[0], range[1] );
207  result.axes.swapcol( 0, 1 );
208  }
209  if( range[1] > range[2] )
210  {
211  std::swap( range[1], range[2] );
212  result.axes.swapcol( 1, 2 );
213  }
214 
215  // scale axis to encompass all points, divide in half
216 #if MB_ORIENTED_BOX_UNIT_VECTORS
217  result.length = range;
218 #else
219  result.axes.colscale( 0, range[0] );
220  result.axes.colscale( 1, range[1] );
221  result.axes.colscale( 2, range[2] );
222 #endif
223 
224 #if MB_ORIENTED_BOX_OUTER_RADIUS
225  result.radius = range.length();
226 #endif
227 
228  return MB_SUCCESS;
229 }
230 
231 ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result, Interface* instance, const Range& vertices )
232 {
233  const Range::iterator begin = vertices.lower_bound( MBVERTEX );
234  const Range::iterator end = vertices.upper_bound( MBVERTEX );
235  size_t count = 0;
236 
237  // compute mean
238  CartVect v;
239  result.center = CartVect( 0, 0, 0 );
240  for( Range::iterator i = begin; i != end; ++i )
241  {
242  ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
243  if( MB_SUCCESS != rval ) return rval;
244  result.center += v;
245  ++count;
246  }
247  result.center /= count;
248 
249  // compute covariance matrix
250  Matrix3 a( 0.0 );
251  for( Range::iterator i = begin; i != end; ++i )
252  {
253  ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
254  if( MB_SUCCESS != rval ) return rval;
255 
256  v -= result.center;
257  a += outer_product( v, v );
258  }
259  a /= count;
260 
261  // Get axes (Eigenvectors) from covariance matrix
262  CartVect lambda;
263  a.eigen_decomposition( lambda, result.axes );
264 
265  // Calculate center and extents of box given orientation defined by axes
266  return box_from_axes( result, instance, vertices );
267 }
268 
269 ErrorCode OrientedBox::covariance_data_from_tris( CovarienceData& result, Interface* instance, const Range& elements )
270 {
271  ErrorCode rval;
272  const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
273  const Range::iterator end = elements.lower_bound( CN::TypeDimensionMap[3].first );
274 
275  // compute mean and moments
276  result.matrix = Matrix3( 0.0 );
277  result.center = CartVect( 0.0 );
278  result.area = 0.0;
279  for( Range::iterator i = begin; i != end; ++i )
280  {
281  const EntityHandle* conn = NULL;
282  int conn_len = 0;
283  rval = instance->get_connectivity( *i, conn, conn_len );
284  if( MB_SUCCESS != rval ) return rval;
285 
286  // for each triangle in the 2-D cell
287  for( int j = 2; j < conn_len; ++j )
288  {
289  EntityHandle vertices[3] = { conn[0], conn[j - 1], conn[j] };
290  CartVect coords[3];
291  rval = instance->get_coords( vertices, 3, coords[0].array() );
292  if( MB_SUCCESS != rval ) return rval;
293 
294  // edge vectors
295  const CartVect edge0 = coords[1] - coords[0];
296  const CartVect edge1 = coords[2] - coords[0];
297  const CartVect centroid = ( coords[0] + coords[1] + coords[2] ) / 3;
298  const double tri_area2 = ( edge0 * edge1 ).length();
299  result.area += tri_area2;
300  result.center += tri_area2 * centroid;
301 
302  result.matrix +=
303  tri_area2 * ( 9 * outer_product( centroid, centroid ) + outer_product( coords[0], coords[0] ) +
304  outer_product( coords[1], coords[1] ) + outer_product( coords[2], coords[2] ) );
305  } // for each triangle
306  } // for each element
307 
308  return MB_SUCCESS;
309 }
310 
311 ErrorCode OrientedBox::compute_from_2d_cells( OrientedBox& result, Interface* instance, const Range& elements )
312 {
313  // Get orientation data from elements
314  CovarienceData data;
315  ErrorCode rval = covariance_data_from_tris( data, instance, elements );
316  if( MB_SUCCESS != rval ) return rval;
317 
318  // get vertices from elements
319  Range points;
320  rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
321  if( MB_SUCCESS != rval ) return rval;
322 
323  // Calculate box given points and orientation data
324  return compute_from_covariance_data( result, instance, data, points );
325 }
326 
327 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
328  Interface* instance,
329  CovarienceData& data,
330  const Range& vertices )
331 {
332  if( data.area <= 0.0 )
333  {
334  Matrix3 empty_axes( 0.0 );
335  result = OrientedBox( empty_axes, CartVect( 0. ) );
336  return MB_SUCCESS;
337  }
338 
339  // get center from sum
340  result.center = data.center / data.area;
341 
342  // get covariance matrix from moments
343  data.matrix /= 12 * data.area;
344  data.matrix -= outer_product( result.center, result.center );
345 
346  // get axes (Eigenvectors) from covariance matrix
347  CartVect lambda;
348  data.matrix.eigen_decomposition( lambda, result.axes );
349 
350  // We now have only the axes. Calculate proper center
351  // and extents for enclosed points.
352  return box_from_axes( result, instance, vertices );
353 }
354 
355 bool OrientedBox::contained( const CartVect& point, double tol ) const
356 {
357  CartVect from_center = point - center;
358 #if MB_ORIENTED_BOX_UNIT_VECTORS
359  return fabs( from_center % axes.col( 0 ) ) - length[0] <= tol &&
360  fabs( from_center % axes.col( 1 ) ) - length[1] <= tol &&
361  fabs( from_center % axes.col( 2 ) ) - length[2] <= tol;
362 #else
363  for( int i = 0; i < 3; ++i )
364  {
365  double length = axes.col( i ).length();
366  if( fabs( from_center % axes.col( i ) ) - length * length > length * tol ) return false;
367  }
368  return true;
369 #endif
370 }
371 
372 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
373  Interface* moab_instance,
374  const CovarienceData* data,
375  unsigned data_length,
376  const Range& vertices )
377 {
378  // Sum input CovarienceData structures
379  CovarienceData data_sum( Matrix3( 0.0 ), CartVect( 0.0 ), 0.0 );
380  for( const CovarienceData* const end = data + data_length; data != end; ++data )
381  {
382  data_sum.matrix += data->matrix;
383  data_sum.center += data->center;
384  data_sum.area += data->area;
385  }
386  // Compute box from sum of structs
387  return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
388 }
389 
390 // bool OrientedBox::contained( const OrientedBox& box, double tol ) const
391 //{
392 // for (int i = -1; i < 2; i += 2)
393 // {
394 // for (int j = -1; j < 2; j += 2)
395 // {
396 // for (int k = -1; k < 2; k += 2)
397 // {
398 // CartVect corner( center );
399 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
400 // corner += i * box.length[0] * box.axis.col(0);
401 // corner += j * box.length[1] * box.axis.col(1);
402 // corner += k * box.length[2] * box.axis.col(2);
403 //#else
404 // corner += i * box.axis.col(0);
405 // corner += j * box.axis.col(1);
406 // corner += k * box.axis.col(2);
407 //#endif
408 // if (!contained( corner, tol ))
409 // return false;
410 // }
411 // }
412 // }
413 // return true;
414 //}
415 
416 /* This is a helper function to check limits on ray length, turning the box-ray
417  * intersection test into a box-segment intersection test. Use this to test the
418  * limits against one side (plane) of the box. The side of the box (plane) is
419  * normal to an axis.
420  *
421  * normal_par_pos Coordinate of particle's position along axis normal to side of box
422  * normal_par_dir Coordinate of particle's direction along axis normal to side of box
423  * half_extent Distance between center of box and side of box
424  * nonneg_ray_len Maximum ray length in positive direction (in front of origin)
425  * neg_ray_len Maximum ray length in negative direction (behind origin)
426  * return true if intersection with plane occurs within distance limits
427  *
428  * ray equation: intersection = origin + dist*direction
429  * plane equation: intersection.plane_normal = half_extent
430  *
431  * Assume plane_normal and direction are unit vectors. Combine equations.
432  *
433  * (origin + dist*direction).plane_normal = half_extent
434  * origin.plane_normal + dist*direction.plane_normal = half_extent
435  * dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
436  *
437  * Although this solves for distance, avoid floating point division.
438  *
439  * dist*direction.plane_normal = half_extent - origin.plane_normal
440  *
441  * Use inequalities to test dist against ray length limits. Be aware that
442  * inequalities change due to sign of direction.plane_normal.
443  */
444 inline bool check_ray_limits( const double normal_par_pos,
445  const double normal_par_dir,
446  const double half_extent,
447  const double* nonneg_ray_len,
448  const double* neg_ray_len )
449 {
450 
451  const double extent_pos_diff = half_extent - normal_par_pos;
452 
453  // limit in positive direction
454  if( nonneg_ray_len )
455  { // should be 0 <= t <= nonneg_ray_len
456  assert( 0 <= *nonneg_ray_len );
457  if( normal_par_dir > 0 )
458  { // if/else if needed for pos/neg divisor
459  if( *nonneg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff >= 0 ) return true;
460  }
461  else if( normal_par_dir < 0 )
462  {
463  if( *nonneg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff <= 0 ) return true;
464  }
465  }
466  else
467  { // should be 0 <= t
468  if( normal_par_dir > 0 )
469  { // if/else if needed for pos/neg divisor
470  if( extent_pos_diff >= 0 ) return true;
471  }
472  else if( normal_par_dir < 0 )
473  {
474  if( extent_pos_diff <= 0 ) return true;
475  }
476  }
477 
478  // limit in negative direction
479  if( neg_ray_len )
480  { // should be neg_ray_len <= t < 0
481  assert( 0 >= *neg_ray_len );
482  if( normal_par_dir > 0 )
483  { // if/else if needed for pos/neg divisor
484  if( *neg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff < 0 ) return true;
485  }
486  else if( normal_par_dir < 0 )
487  {
488  if( *neg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff > 0 ) return true;
489  }
490  }
491 
492  return false;
493 }
494 
495 /* This implementation copied from cgmMC (overlap.C).
496  * Original author: Tim Tautges?
497  */
498 bool OrientedBox::intersect_ray( const CartVect& ray_origin,
499  const CartVect& ray_direction,
500  const double reps,
501  const double* nonneg_ray_len,
502  const double* neg_ray_len ) const
503 {
504  // test distance from box center to line
505  const CartVect cx = center - ray_origin;
506  const double dist_s = cx % ray_direction;
507  const double dist_sq = cx % cx - ( dist_s * dist_s );
508  const double max_diagsq = outer_radius_squared( reps );
509 
510  // For the largest sphere, no intersections exist if discriminant is negative.
511  // Geometrically, if distance from box center to line is greater than the
512  // longest diagonal, there is no intersection.
513  // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
514  if( dist_sq > max_diagsq ) return false;
515 
516  // If the closest possible intersection must be closer than nonneg_ray_len. Be
517  // careful with absolute value, squaring distances, and subtracting squared
518  // distances.
519  if( nonneg_ray_len )
520  {
521  assert( 0 <= *nonneg_ray_len );
522  double max_len;
523  if( neg_ray_len )
524  {
525  assert( 0 >= *neg_ray_len );
526  max_len = std::max( *nonneg_ray_len, -( *neg_ray_len ) );
527  }
528  else
529  {
530  max_len = *nonneg_ray_len;
531  }
532  const double temp = fabs( dist_s ) - max_len;
533  if( 0.0 < temp && temp * temp > max_diagsq ) return false;
534  }
535 
536  // if smaller than shortest diagonal, we do hit
537  if( dist_sq < inner_radius_squared( reps ) )
538  {
539  // nonnegative direction
540  if( dist_s >= 0.0 )
541  {
542  if( nonneg_ray_len )
543  {
544  if( *nonneg_ray_len > dist_s ) return true;
545  }
546  else
547  {
548  return true;
549  }
550  // negative direction
551  }
552  else
553  {
554  if( neg_ray_len && *neg_ray_len < dist_s ) return true;
555  }
556  }
557 
558  // get transpose of axes
559  Matrix3 B = axes.transpose();
560 
561  // transform ray to box coordintae system
562  CartVect par_pos = B * ( ray_origin - center );
563  CartVect par_dir = B * ray_direction;
564 
565  // Fast Rejection Test: Ray will not intersect if it is going away from the box.
566  // This will not work for rays with neg_ray_len. length[0] is half of box width
567  // along axes.col(0).
568  const double half_x = length[0] + reps;
569  const double half_y = length[1] + reps;
570  const double half_z = length[2] + reps;
571  if( !neg_ray_len )
572  {
573  if( ( par_pos[0] > half_x && par_dir[0] >= 0 ) || ( par_pos[0] < -half_x && par_dir[0] <= 0 ) ) return false;
574 
575  if( ( par_pos[1] > half_y && par_dir[1] >= 0 ) || ( par_pos[1] < -half_y && par_dir[1] <= 0 ) ) return false;
576 
577  if( ( par_pos[2] > half_z && par_dir[2] >= 0 ) || ( par_pos[2] < -half_z && par_dir[2] <= 0 ) ) return false;
578  }
579 
580  // test if ray_origin is inside box
581  if( par_pos[0] <= half_x && par_pos[0] >= -half_x && par_pos[1] <= half_y && par_pos[1] >= -half_y &&
582  par_pos[2] <= half_z && par_pos[2] >= -half_z )
583  return true;
584 
585  // test two xy plane
586  if( fabs( par_dir[0] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <=
587  fabs( par_dir[2] * half_x ) && // test against x extents using z
588  fabs( par_dir[1] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <=
589  fabs( par_dir[2] * half_y ) && // test against y extents using z
590  check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
591  return true;
592  if( fabs( par_dir[0] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <= fabs( par_dir[2] * half_x ) &&
593  fabs( par_dir[1] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <= fabs( par_dir[2] * half_y ) &&
594  check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
595  return true;
596 
597  // test two xz plane
598  if( fabs( par_dir[0] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
599  fabs( par_dir[2] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
600  check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
601  return true;
602  if( fabs( par_dir[0] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
603  fabs( par_dir[2] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
604  check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
605  return true;
606 
607  // test two yz plane
608  if( fabs( par_dir[1] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
609  fabs( par_dir[2] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
610  check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
611  return true;
612  if( fabs( par_dir[1] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
613  fabs( par_dir[2] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
614  check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
615  return true;
616 
617  return false;
618 }
619 
620 ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
621 {
622  ErrorCode rval;
623  int signs[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
624  { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
625 
626  std::vector< EntityHandle > vertices;
627  for( int i = 0; i < 8; ++i )
628  {
629  CartVect coords( center );
630  for( int j = 0; j < 3; ++j )
631  {
632 #if MB_ORIENTED_BOX_UNIT_VECTORS
633  coords += signs[i][j] * ( axes.col( j ) * length[j] );
634 #else
635  coords += signs[i][j] * axes.col( j );
636 #endif
637  }
638  EntityHandle handle;
639  rval = instance->create_vertex( coords.array(), handle );
640  if( MB_SUCCESS != rval )
641  {
642  instance->delete_entities( &vertices[0], vertices.size() );
643  return rval;
644  }
645  vertices.push_back( handle );
646  }
647 
648  rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
649  if( MB_SUCCESS != rval )
650  {
651  instance->delete_entities( &vertices[0], vertices.size() );
652  return rval;
653  }
654 
655  return MB_SUCCESS;
656 }
657 
658 void OrientedBox::closest_location_in_box( const CartVect& input_position, CartVect& output_position ) const
659 {
660  // get coordinates on box axes
661  const CartVect from_center = input_position - center;
662 
663 #if MB_ORIENTED_BOX_UNIT_VECTORS
664  CartVect local( from_center % axes.col( 0 ), from_center % axes.col( 1 ), from_center % axes.col( 2 ) );
665 
666  for( int i = 0; i < 3; ++i )
667  {
668  if( local[i] < -length[i] )
669  local[i] = -length[i];
670  else if( local[i] > length[i] )
671  local[i] = length[i];
672  }
673 #else
674  CartVect local( ( from_center % axes.col( 0 ) ) / ( axes.col( 0 ) % axes.col( 0 ) ),
675  ( from_center % axes.col( 1 ) ) / ( axes.col( 1 ) % axes.col( 1 ) ),
676  ( from_center % axes.col( 2 ) ) / ( axes.col( 2 ) % axes.col( 2 ) ) );
677 
678  for( int i = 0; i < 3; ++i )
679  {
680  if( local[i] < -1.0 )
681  local[i] = -1.0;
682  else if( local[i] > 1.0 )
683  local[i] = 1.0;
684  }
685 #endif
686 
687  output_position = center + local[0] * axes.col( 0 ) + local[1] * axes.col( 1 ) + local[2] * axes.col( 2 );
688 }
689 
690 } // namespace moab