Mesh Oriented datABase  (version 5.6.0)
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  // project points onto axes to get box extents
169  CartVect min( std::numeric_limits< double >::max() ), max( -std::numeric_limits< double >::max() );
170  for( Range::iterator i = points.begin(); i != points.end(); ++i )
171  {
172  CartVect coords;
173  MB_CHK_ERR( instance->get_coords( &*i, 1, coords.array() ) );
174 
175  for( int d = 0; d < 3; ++d )
176  {
177  const double t = point_perp( coords, result.center, result.axes.col( d ) );
178  if( t < min[d] ) min[d] = t;
179  if( t > max[d] ) max[d] = t;
180  }
181  }
182 
183  // We now have a box defined by three orthogonal line segments
184  // that intersect at the center of the box. Each line segment
185  // is defined as result.center + t * result.axes[i], where the
186  // range of t is [min[i], max[i]].
187 
188  // Calculate new center
189  const CartVect mid = 0.5 * ( min + max );
190  result.center += mid[0] * result.axes.col( 0 ) + mid[1] * result.axes.col( 1 ) + mid[2] * result.axes.col( 2 );
191 
192  // reorder axes by length
193  CartVect range = 0.5 * ( max - min );
194  if( range[2] < range[1] )
195  {
196  if( range[2] < range[0] )
197  {
198  std::swap( range[0], range[2] );
199  result.axes.swapcol( 0, 2 );
200  }
201  }
202  else if( range[1] < range[0] )
203  {
204  std::swap( range[0], range[1] );
205  result.axes.swapcol( 0, 1 );
206  }
207  if( range[1] > range[2] )
208  {
209  std::swap( range[1], range[2] );
210  result.axes.swapcol( 1, 2 );
211  }
212 
213  // scale axis to encompass all points, divide in half
214 #if MB_ORIENTED_BOX_UNIT_VECTORS
215  result.length = range;
216 #else
217  result.axes.colscale( 0, range[0] );
218  result.axes.colscale( 1, range[1] );
219  result.axes.colscale( 2, range[2] );
220 #endif
221 
222 #if MB_ORIENTED_BOX_OUTER_RADIUS
223  result.radius = range.length();
224 #endif
225 
226  return MB_SUCCESS;
227 }
228 
229 ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result, Interface* instance, const Range& vertices )
230 {
231  const Range::iterator begin = vertices.lower_bound( MBVERTEX );
232  const Range::iterator end = vertices.upper_bound( MBVERTEX );
233  size_t count = 0;
234 
235  // compute mean
236  CartVect v;
237  result.center = CartVect( 0, 0, 0 );
238  for( Range::iterator i = begin; i != end; ++i )
239  {
240  ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
241  if( MB_SUCCESS != rval ) return rval;
242  result.center += v;
243  ++count;
244  }
245  result.center /= count;
246 
247  // compute covariance matrix
248  Matrix3 a( 0.0 );
249  for( Range::iterator i = begin; i != end; ++i )
250  {
251  ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
252  if( MB_SUCCESS != rval ) return rval;
253 
254  v -= result.center;
255  a += outer_product( v, v );
256  }
257  a /= count;
258 
259  // Get axes (Eigenvectors) from covariance matrix
260  CartVect lambda;
261  a.eigen_decomposition( lambda, result.axes );
262 
263  // Calculate center and extents of box given orientation defined by axes
264  return box_from_axes( result, instance, vertices );
265 }
266 
267 ErrorCode OrientedBox::covariance_data_from_tris( CovarienceData& result, Interface* instance, const Range& elements )
268 {
269  ErrorCode rval;
270  const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
271  const Range::iterator end = elements.lower_bound( CN::TypeDimensionMap[3].first );
272 
273  // compute mean and moments
274  result.matrix = Matrix3( 0.0 );
275  result.center = CartVect( 0.0 );
276  result.area = 0.0;
277  for( Range::iterator i = begin; i != end; ++i )
278  {
279  const EntityHandle* conn = NULL;
280  int conn_len = 0;
281  rval = instance->get_connectivity( *i, conn, conn_len );
282  if( MB_SUCCESS != rval ) return rval;
283 
284  // for each triangle in the 2-D cell
285  for( int j = 2; j < conn_len; ++j )
286  {
287  EntityHandle vertices[3] = { conn[0], conn[j - 1], conn[j] };
288  CartVect coords[3];
289  rval = instance->get_coords( vertices, 3, coords[0].array() );
290  if( MB_SUCCESS != rval ) return rval;
291 
292  // edge vectors
293  const CartVect edge0 = coords[1] - coords[0];
294  const CartVect edge1 = coords[2] - coords[0];
295  const CartVect centroid = ( coords[0] + coords[1] + coords[2] ) / 3;
296  const double tri_area2 = ( edge0 * edge1 ).length();
297  result.area += tri_area2;
298  result.center += tri_area2 * centroid;
299 
300  result.matrix +=
301  tri_area2 * ( 9 * outer_product( centroid, centroid ) + outer_product( coords[0], coords[0] ) +
302  outer_product( coords[1], coords[1] ) + outer_product( coords[2], coords[2] ) );
303  } // for each triangle
304  } // for each element
305 
306  return MB_SUCCESS;
307 }
308 
309 ErrorCode OrientedBox::compute_from_2d_cells( OrientedBox& result, Interface* instance, const Range& elements )
310 {
311  // Get orientation data from elements
312  CovarienceData data;
313  ErrorCode rval = covariance_data_from_tris( data, instance, elements );
314  if( MB_SUCCESS != rval ) return rval;
315 
316  // get vertices from elements
317  Range points;
318  rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
319  if( MB_SUCCESS != rval ) return rval;
320 
321  // Calculate box given points and orientation data
322  return compute_from_covariance_data( result, instance, data, points );
323 }
324 
325 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
326  Interface* instance,
327  CovarienceData& data,
328  const Range& vertices )
329 {
330  if( data.area <= 0.0 )
331  {
332  Matrix3 empty_axes( 0.0 );
333  result = OrientedBox( empty_axes, CartVect( 0. ) );
334  return MB_SUCCESS;
335  }
336 
337  // get center from sum
338  result.center = data.center / data.area;
339 
340  // get covariance matrix from moments
341  data.matrix /= 12 * data.area;
342  data.matrix -= outer_product( result.center, result.center );
343 
344  // get axes (Eigenvectors) from covariance matrix
345  CartVect lambda;
346  data.matrix.eigen_decomposition( lambda, result.axes );
347 
348  // We now have only the axes. Calculate proper center
349  // and extents for enclosed points.
350  return box_from_axes( result, instance, vertices );
351 }
352 
353 bool OrientedBox::contained( const CartVect& point, double tol ) const
354 {
355  CartVect from_center = point - center;
356 #if MB_ORIENTED_BOX_UNIT_VECTORS
357  return fabs( from_center % axes.col( 0 ) ) - length[0] <= tol &&
358  fabs( from_center % axes.col( 1 ) ) - length[1] <= tol &&
359  fabs( from_center % axes.col( 2 ) ) - length[2] <= tol;
360 #else
361  for( int i = 0; i < 3; ++i )
362  {
363  double length = axes.col( i ).length();
364  if( fabs( from_center % axes.col( i ) ) - length * length > length * tol ) return false;
365  }
366  return true;
367 #endif
368 }
369 
370 ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
371  Interface* moab_instance,
372  const CovarienceData* data,
373  unsigned data_length,
374  const Range& vertices )
375 {
376  // Sum input CovarienceData structures
377  CovarienceData data_sum( Matrix3( 0.0 ), CartVect( 0.0 ), 0.0 );
378  for( const CovarienceData* const end = data + data_length; data != end; ++data )
379  {
380  data_sum.matrix += data->matrix;
381  data_sum.center += data->center;
382  data_sum.area += data->area;
383  }
384  // Compute box from sum of structs
385  return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
386 }
387 
388 // bool OrientedBox::contained( const OrientedBox& box, double tol ) const
389 //{
390 // for (int i = -1; i < 2; i += 2)
391 // {
392 // for (int j = -1; j < 2; j += 2)
393 // {
394 // for (int k = -1; k < 2; k += 2)
395 // {
396 // CartVect corner( center );
397 //#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
398 // corner += i * box.length[0] * box.axis.col(0);
399 // corner += j * box.length[1] * box.axis.col(1);
400 // corner += k * box.length[2] * box.axis.col(2);
401 //#else
402 // corner += i * box.axis.col(0);
403 // corner += j * box.axis.col(1);
404 // corner += k * box.axis.col(2);
405 //#endif
406 // if (!contained( corner, tol ))
407 // return false;
408 // }
409 // }
410 // }
411 // return true;
412 //}
413 
414 /* This is a helper function to check limits on ray length, turning the box-ray
415  * intersection test into a box-segment intersection test. Use this to test the
416  * limits against one side (plane) of the box. The side of the box (plane) is
417  * normal to an axis.
418  *
419  * normal_par_pos Coordinate of particle's position along axis normal to side of box
420  * normal_par_dir Coordinate of particle's direction along axis normal to side of box
421  * half_extent Distance between center of box and side of box
422  * nonneg_ray_len Maximum ray length in positive direction (in front of origin)
423  * neg_ray_len Maximum ray length in negative direction (behind origin)
424  * return true if intersection with plane occurs within distance limits
425  *
426  * ray equation: intersection = origin + dist*direction
427  * plane equation: intersection.plane_normal = half_extent
428  *
429  * Assume plane_normal and direction are unit vectors. Combine equations.
430  *
431  * (origin + dist*direction).plane_normal = half_extent
432  * origin.plane_normal + dist*direction.plane_normal = half_extent
433  * dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
434  *
435  * Although this solves for distance, avoid floating point division.
436  *
437  * dist*direction.plane_normal = half_extent - origin.plane_normal
438  *
439  * Use inequalities to test dist against ray length limits. Be aware that
440  * inequalities change due to sign of direction.plane_normal.
441  */
442 inline bool check_ray_limits( const double normal_par_pos,
443  const double normal_par_dir,
444  const double half_extent,
445  const double* nonneg_ray_len,
446  const double* neg_ray_len )
447 {
448 
449  const double extent_pos_diff = half_extent - normal_par_pos;
450 
451  // limit in positive direction
452  if( nonneg_ray_len )
453  { // should be 0 <= t <= nonneg_ray_len
454  assert( 0 <= *nonneg_ray_len );
455  if( normal_par_dir > 0 )
456  { // if/else if needed for pos/neg divisor
457  if( *nonneg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff >= 0 ) return true;
458  }
459  else if( normal_par_dir < 0 )
460  {
461  if( *nonneg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff <= 0 ) return true;
462  }
463  }
464  else
465  { // should be 0 <= t
466  if( normal_par_dir > 0 )
467  { // if/else if needed for pos/neg divisor
468  if( extent_pos_diff >= 0 ) return true;
469  }
470  else if( normal_par_dir < 0 )
471  {
472  if( extent_pos_diff <= 0 ) return true;
473  }
474  }
475 
476  // limit in negative direction
477  if( neg_ray_len )
478  { // should be neg_ray_len <= t < 0
479  assert( 0 >= *neg_ray_len );
480  if( normal_par_dir > 0 )
481  { // if/else if needed for pos/neg divisor
482  if( *neg_ray_len * normal_par_dir <= extent_pos_diff && extent_pos_diff < 0 ) return true;
483  }
484  else if( normal_par_dir < 0 )
485  {
486  if( *neg_ray_len * normal_par_dir >= extent_pos_diff && extent_pos_diff > 0 ) return true;
487  }
488  }
489 
490  return false;
491 }
492 
493 /* This implementation copied from cgmMC (overlap.C).
494  * Original author: Tim Tautges?
495  */
496 bool OrientedBox::intersect_ray( const CartVect& ray_origin,
497  const CartVect& ray_direction,
498  const double reps,
499  const double* nonneg_ray_len,
500  const double* neg_ray_len ) const
501 {
502  // test distance from box center to line
503  const CartVect cx = center - ray_origin;
504  const double dist_s = cx % ray_direction;
505  const double dist_sq = cx % cx - ( dist_s * dist_s );
506  const double max_diagsq = outer_radius_squared( reps );
507 
508  // For the largest sphere, no intersections exist if discriminant is negative.
509  // Geometrically, if distance from box center to line is greater than the
510  // longest diagonal, there is no intersection.
511  // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
512  if( dist_sq > max_diagsq ) return false;
513 
514  // If the closest possible intersection must be closer than nonneg_ray_len. Be
515  // careful with absolute value, squaring distances, and subtracting squared
516  // distances.
517  if( nonneg_ray_len )
518  {
519  assert( 0 <= *nonneg_ray_len );
520  double max_len;
521  if( neg_ray_len )
522  {
523  assert( 0 >= *neg_ray_len );
524  max_len = std::max( *nonneg_ray_len, -( *neg_ray_len ) );
525  }
526  else
527  {
528  max_len = *nonneg_ray_len;
529  }
530  const double temp = fabs( dist_s ) - max_len;
531  if( 0.0 < temp && temp * temp > max_diagsq ) return false;
532  }
533 
534  // if smaller than shortest diagonal, we do hit
535  if( dist_sq < inner_radius_squared( reps ) )
536  {
537  // nonnegative direction
538  if( dist_s >= 0.0 )
539  {
540  if( nonneg_ray_len )
541  {
542  if( *nonneg_ray_len > dist_s ) return true;
543  }
544  else
545  {
546  return true;
547  }
548  // negative direction
549  }
550  else
551  {
552  if( neg_ray_len && *neg_ray_len < dist_s ) return true;
553  }
554  }
555 
556  // get transpose of axes
557  Matrix3 B = axes.transpose();
558 
559  // transform ray to box coordintae system
560  CartVect par_pos = B * ( ray_origin - center );
561  CartVect par_dir = B * ray_direction;
562 
563  // Fast Rejection Test: Ray will not intersect if it is going away from the box.
564  // This will not work for rays with neg_ray_len. length[0] is half of box width
565  // along axes.col(0).
566  const double half_x = length[0] + reps;
567  const double half_y = length[1] + reps;
568  const double half_z = length[2] + reps;
569  if( !neg_ray_len )
570  {
571  if( ( par_pos[0] > half_x && par_dir[0] >= 0 ) || ( par_pos[0] < -half_x && par_dir[0] <= 0 ) ) return false;
572 
573  if( ( par_pos[1] > half_y && par_dir[1] >= 0 ) || ( par_pos[1] < -half_y && par_dir[1] <= 0 ) ) return false;
574 
575  if( ( par_pos[2] > half_z && par_dir[2] >= 0 ) || ( par_pos[2] < -half_z && par_dir[2] <= 0 ) ) return false;
576  }
577 
578  // test if ray_origin is inside box
579  if( par_pos[0] <= half_x && par_pos[0] >= -half_x && par_pos[1] <= half_y && par_pos[1] >= -half_y &&
580  par_pos[2] <= half_z && par_pos[2] >= -half_z )
581  return true;
582 
583  // test two xy plane
584  if( fabs( par_dir[0] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <=
585  fabs( par_dir[2] * half_x ) && // test against x extents using z
586  fabs( par_dir[1] * ( half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <=
587  fabs( par_dir[2] * half_y ) && // test against y extents using z
588  check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
589  return true;
590  if( fabs( par_dir[0] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[0] ) <= fabs( par_dir[2] * half_x ) &&
591  fabs( par_dir[1] * ( -half_z - par_pos[2] ) + par_dir[2] * par_pos[1] ) <= fabs( par_dir[2] * half_y ) &&
592  check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
593  return true;
594 
595  // test two xz plane
596  if( fabs( par_dir[0] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
597  fabs( par_dir[2] * ( half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
598  check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
599  return true;
600  if( fabs( par_dir[0] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[0] ) <= fabs( par_dir[1] * half_x ) &&
601  fabs( par_dir[2] * ( -half_y - par_pos[1] ) + par_dir[1] * par_pos[2] ) <= fabs( par_dir[1] * half_z ) &&
602  check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
603  return true;
604 
605  // test two yz plane
606  if( fabs( par_dir[1] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
607  fabs( par_dir[2] * ( half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
608  check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
609  return true;
610  if( fabs( par_dir[1] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[1] ) <= fabs( par_dir[0] * half_y ) &&
611  fabs( par_dir[2] * ( -half_x - par_pos[0] ) + par_dir[0] * par_pos[2] ) <= fabs( par_dir[0] * half_z ) &&
612  check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
613  return true;
614 
615  return false;
616 }
617 
618 ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
619 {
620  ErrorCode rval;
621  int signs[8][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, { -1, 1, -1 },
622  { -1, -1, 1 }, { 1, -1, 1 }, { 1, 1, 1 }, { -1, 1, 1 } };
623 
624  std::vector< EntityHandle > vertices;
625  for( int i = 0; i < 8; ++i )
626  {
627  CartVect coords( center );
628  for( int j = 0; j < 3; ++j )
629  {
630 #if MB_ORIENTED_BOX_UNIT_VECTORS
631  coords += signs[i][j] * ( axes.col( j ) * length[j] );
632 #else
633  coords += signs[i][j] * axes.col( j );
634 #endif
635  }
636  EntityHandle handle;
637  rval = instance->create_vertex( coords.array(), handle );
638  if( MB_SUCCESS != rval )
639  {
640  instance->delete_entities( &vertices[0], vertices.size() );
641  return rval;
642  }
643  vertices.push_back( handle );
644  }
645 
646  rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
647  if( MB_SUCCESS != rval )
648  {
649  instance->delete_entities( &vertices[0], vertices.size() );
650  return rval;
651  }
652 
653  return MB_SUCCESS;
654 }
655 
656 void OrientedBox::closest_location_in_box( const CartVect& input_position, CartVect& output_position ) const
657 {
658  // get coordinates on box axes
659  const CartVect from_center = input_position - center;
660 
661 #if MB_ORIENTED_BOX_UNIT_VECTORS
662  CartVect local( from_center % axes.col( 0 ), from_center % axes.col( 1 ), from_center % axes.col( 2 ) );
663 
664  for( int i = 0; i < 3; ++i )
665  {
666  if( local[i] < -length[i] )
667  local[i] = -length[i];
668  else if( local[i] > length[i] )
669  local[i] = length[i];
670  }
671 #else
672  CartVect local( ( from_center % axes.col( 0 ) ) / ( axes.col( 0 ) % axes.col( 0 ) ),
673  ( from_center % axes.col( 1 ) ) / ( axes.col( 1 ) % axes.col( 1 ) ),
674  ( from_center % axes.col( 2 ) ) / ( axes.col( 2 ) % axes.col( 2 ) ) );
675 
676  for( int i = 0; i < 3; ++i )
677  {
678  if( local[i] < -1.0 )
679  local[i] = -1.0;
680  else if( local[i] > 1.0 )
681  local[i] = 1.0;
682  }
683 #endif
684 
685  output_position = center + local[0] * axes.col( 0 ) + local[1] * axes.col( 1 ) + local[2] * axes.col( 2 );
686 }
687 
688 } // namespace moab