fathomteam-moab-7bde9dfb84a8/src/OrientedBox.cpp

/*
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

/*
 * The algorithms for the calculation of the oriented box from a
 * set of points or a set of cells was copied from the implementation
 " in the "Visualization Toolkit".  J.K. - 2006-07-19
 *
 * Program:   Visualization Toolkit
 * Module:    $RCSfile$
 *
 * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 * All rights reserved.
 * See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
 */

/**\file OrientedBox.cpp
 *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
 *\date 2006-07-18
 */

#include "moab/Interface.hpp"
#include "moab/CN.hpp"
#include "moab/OrientedBox.hpp"
#include "moab/Range.hpp"
#include <ostream>
#include <assert.h>
#include <limits>

namespace moab {

std::ostream& operator<<( std::ostream& s, const OrientedBox& b )
{
  return s << b.center
           << " + "
           << b.axes.col(0)
#if MB_ORIENTED_BOX_UNIT_VECTORS
           << ":" << b.length[0]
#endif
           << " x "
           << b.axes.col(1)
#if MB_ORIENTED_BOX_UNIT_VECTORS
           << ":" << b.length[1]
#endif
           << " x "
           << b.axes.col(2)
#if MB_ORIENTED_BOX_UNIT_VECTORS
           << ":" << b.length[2]
#endif
            ;
}

/**\brief Find closest point on line
 *
 * Find the point on the line for which a line trough the
 * input point \a p and the result position is orthogonal to
 * the input line.
 * \param p  The point for which to find the perpendicular
 * \param b  A point on the line
 * \param m  The direction of the line
 * \return   The location on the line specified as 't' in the
 *           formula t * m + b
 */
static double point_perp( const CartVect& p,   // closest to this point
                          const CartVect& b,   // point on line
                          const CartVect& m )  // line direction
{
#if MB_ORIENTED_BOX_UNIT_VECTORS
  double t = (m % (p - b));
#else
  double t = (m % (p - b)) / (m % m);
#endif
  return Util::is_finite(t) ? t : 0.0;
}

void OrientedBox::order_axes_by_length(double ax1_len, double ax2_len, double ax3_len)
{

  CartVect len( ax1_len, ax2_len, ax3_len );

  if (len[2] < len[1])
    {
      if (len[2] < len[0]) {
	std::swap( len[0], len[2] );
	axes.swapcol( 0, 2 );
      }
    }
  else if (len[1] < len[0]) {
    std::swap( len[0], len[1] );
    axes.swapcol( 0, 1 );
  }
  if (len[1] > len[2]) {
    std::swap( len[1], len[2] );
    axes.swapcol( 1, 2 );
  }

#if MB_ORIENTED_BOX_UNIT_VECTORS
    length = len;
    if (len[0] > 0.0)
      axes.colscale(0, 1.0/len[0]);
    if (len[1] > 0.0)
      axes.colscale(1, 1.0/len[1]);
    if (len[2] > 0.0)
      axes.colscale(2 ,1.0/len[2]);
#endif

#if MB_ORIENTED_BOX_OUTER_RADIUS
    radius = len.length();
#endif
}

OrientedBox::OrientedBox( const CartVect axes_in[3], const CartVect& mid )
  : center(mid)
{

  axes = Matrix3( axes_in[0], axes_in[1], axes_in[2], false );

  order_axes_by_length( axes_in[0].length(),axes_in[1].length(),axes_in[2].length() );

}

OrientedBox::OrientedBox( const Matrix3& axes_mat, const CartVect& mid )
  : center(mid), axes(axes_mat)
{
  order_axes_by_length( axes.col(0).length(), axes.col(1).length(), axes.col(2).length() );
}

ErrorCode OrientedBox::tag_handle( Tag& handle_out,
                                       Interface* instance,
                                       const char* name)
{
    // We're going to assume this when mapping the OrientedBox
    // to tag data, so assert it.
#if MB_ORIENTED_BOX_OUTER_RADIUS
  const int rad_size = 1;
#else
  const int rad_size = 0;
#endif
#if MB_ORIENTED_BOX_UNIT_VECTORS
  const int SIZE = rad_size + 15;
#else
  const int SIZE = rad_size + 12;
#endif
  assert( sizeof(OrientedBox) == SIZE*sizeof(double) );

  return instance->tag_get_handle( name, SIZE, MB_TYPE_DOUBLE,
                                   handle_out, MB_TAG_DENSE|MB_TAG_CREAT );
}

/**\brief Common code for box calculation
 *
 * Given the orientation of the box and an approximate center,
 * calculate the exact center and extents of the box.
 *
 *\param result.center  As input, the approximate center of the box.
 *                      As output, the exact center of the box.
 *\param result.axes    As input, directions of principal axes corresponding
 *                      to the orientation of the box.  Axes are assumed to
 *                      be unit-length on input.  Output will include extents
 *                      of box.
 *\param points  The set of points the box should contain.
 */
static ErrorCode box_from_axes( OrientedBox& result,
                                  Interface* instance,
                                  const Range& points )
{
  ErrorCode rval;

    // project points onto axes to get box extents
  CartVect min(std::numeric_limits<double>::max()),
             max(-std::numeric_limits<double>::max());
  for (Range::iterator i = points.begin(); i != points.end(); ++i)
  {
    CartVect coords;
    rval = instance->get_coords( &*i, 1, coords.array() );MB_CHK_ERR(rval);

    for (int d = 0; d < 3; ++d)
    {
      const double t = point_perp( coords, result.center, result.axes.col(d) );
      if (t < min[d])
        min[d] = t;
      if (t > max[d])
        max[d] = t;
    }
  }

    // We now have a box defined by three orthogonal line segments
    // that intersect at the center of the box.  Each line segment
    // is defined as result.center + t * result.axes[i], where the
    // range of t is [min[i], max[i]].

    // Calculate new center
  const CartVect mid = 0.5 * (min + max);
  result.center += mid[0] * result.axes.col(0) +
                   mid[1] * result.axes.col(1) +
                   mid[2] * result.axes.col(2);

    // reorder axes by length
  CartVect range = 0.5 * (max - min);
  if (range[2] < range[1])
  {
    if (range[2] < range[0]) {
      std::swap( range[0], range[2] );
      result.axes.swapcol( 0, 2 );
    }
  }
  else if (range[1] < range[0]) {
    std::swap( range[0], range[1] );
    result.axes.swapcol( 0, 1 );
  }
  if (range[1] > range[2]) {
    std::swap( range[1], range[2] );
    result.axes.swapcol( 1, 2 );
  }

    // scale axis to encompass all points, divide in half
#if MB_ORIENTED_BOX_UNIT_VECTORS
  result.length = range;
#else
  result.axes.colscale(0, range[0]);
  result.axes.colscale(1, range[1]);
  result.axes.colscale(2, range[2]);
#endif

#if MB_ORIENTED_BOX_OUTER_RADIUS
  result.radius = range.length();
#endif

  return MB_SUCCESS;
}


ErrorCode OrientedBox::compute_from_vertices( OrientedBox& result,
                                                  Interface* instance,
                                                  const Range& vertices )
{
  const Range::iterator begin = vertices.lower_bound( MBVERTEX );
  const Range::iterator end = vertices.upper_bound( MBVERTEX );
  size_t count = 0;

    // compute mean
  CartVect v;
  result.center = CartVect( 0, 0, 0 );
  for (Range::iterator i = begin; i != end; ++i)
  {
    ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
    if (MB_SUCCESS != rval)
      return rval;
    result.center += v;
    ++count;
  }
  result.center /= count;

    // compute covariance matrix
  Matrix3 a( 0.0 );
  for (Range::iterator i = begin; i != end; ++i)
  {
    ErrorCode rval = instance->get_coords( &*i, 1, v.array() );
    if (MB_SUCCESS != rval)
      return rval;

    v -= result.center;
    a += outer_product( v, v );
  }
  a /= count;

    // Get axes (Eigenvectors) from covariance matrix
  CartVect lambda;
  a.eigen_decomposition(lambda, result.axes);

    // Calculate center and extents of box given orientation defined by axes
  return box_from_axes( result, instance, vertices );
}

ErrorCode OrientedBox::covariance_data_from_tris( CovarienceData& result,
                                                 Interface* instance,
                                                 const Range& elements )
{
  ErrorCode rval;
  const Range::iterator begin = elements.lower_bound( CN::TypeDimensionMap[2].first );
  const Range::iterator end = elements.lower_bound( CN::TypeDimensionMap[3].first );

    // compute mean and moments
  result.matrix = Matrix3(0.0);
  result.center = CartVect(0.0);
  result.area = 0.0;
  for (Range::iterator i = begin; i != end; ++i)
  {
    const EntityHandle* conn = NULL;
    int conn_len = 0;
    rval = instance->get_connectivity( *i, conn, conn_len );
    if (MB_SUCCESS != rval)
      return rval;

      // for each triangle in the 2-D cell
    for (int j = 2; j < conn_len; ++j)
    {
      EntityHandle vertices[3] = { conn[0], conn[j-1], conn[j] };
      CartVect coords[3];
      rval = instance->get_coords( vertices, 3, coords[0].array() );
      if (MB_SUCCESS != rval)
        return rval;

        // edge vectors
      const CartVect edge0 = coords[1] - coords[0];
      const CartVect edge1 = coords[2] - coords[0];
      const CartVect centroid = (coords[0] + coords[1] + coords[2]) / 3;
      const double tri_area2 = (edge0 * edge1).length();
      result.area += tri_area2;
      result.center += tri_area2 * centroid;

      result.matrix += tri_area2 * (9 * outer_product( centroid,  centroid  ) +
                                    outer_product( coords[0], coords[0] ) +
                                    outer_product( coords[1], coords[1] ) +
                                    outer_product( coords[2], coords[2] ));
    } // for each triangle
  } // for each element

  return MB_SUCCESS;
}


ErrorCode OrientedBox::compute_from_2d_cells( OrientedBox& result,
                                                  Interface* instance,
                                                  const Range& elements )
{
    // Get orientation data from elements
  CovarienceData data;
  ErrorCode rval = covariance_data_from_tris( data, instance, elements );
  if (MB_SUCCESS != rval)
    return rval;

    // get vertices from elements
  Range points;
  rval = instance->get_adjacencies( elements, 0, false, points, Interface::UNION );
  if (MB_SUCCESS != rval)
    return rval;

    // Calculate box given points and orientation data
  return compute_from_covariance_data( result, instance, data, points );
}

ErrorCode OrientedBox::compute_from_covariance_data(
                                                OrientedBox& result,
                                                Interface* instance,
                                                CovarienceData& data,
                                                const Range& vertices )
{
  if (data.area <= 0.0) {
    Matrix3 empty_axes (0.0);
    result = OrientedBox( empty_axes, CartVect(0.) );
    return MB_SUCCESS;
  }

    // get center from sum
  result.center = data.center / data.area;

    // get covariance matrix from moments
  data.matrix /= 12 * data.area;
  data.matrix -= outer_product( result.center, result.center );

    // get axes (Eigenvectors) from covariance matrix
  CartVect lambda;
  data.matrix.eigen_decomposition(lambda, result.axes);

    // We now have only the axes.  Calculate proper center
    // and extents for enclosed points.
  return box_from_axes( result, instance, vertices );
}

bool OrientedBox::contained( const CartVect& point, double tol ) const
{
  CartVect from_center = point - center;
#if MB_ORIENTED_BOX_UNIT_VECTORS
  return fabs(from_center % axes.col(0)) - length[0] <= tol &&
         fabs(from_center % axes.col(1)) - length[1] <= tol &&
         fabs(from_center % axes.col(2)) - length[2] <= tol ;
#else
  for (int i = 0; i < 3; ++i) {
    double length = axes.col(i).length();
    if (fabs(from_center % axes.col(i)) - length*length > length*tol)
      return false;
  }
  return true;
#endif
}

ErrorCode OrientedBox::compute_from_covariance_data( OrientedBox& result,
                                                Interface* moab_instance,
                                                const CovarienceData* data,
                                                unsigned data_length,
                                                const Range& vertices )
{
    // Sum input CovarienceData structures
  CovarienceData data_sum( Matrix3(0.0), CartVect(0.0), 0.0 );
  for (const CovarienceData* const end = data+data_length; data != end; ++data) {
    data_sum.matrix += data->matrix;
    data_sum.center += data->center;
    data_sum.area += data->area;
  }
    // Compute box from sum of structs
  return compute_from_covariance_data( result, moab_instance, data_sum, vertices );
}


//bool OrientedBox::contained( const OrientedBox& box, double tol ) const
//{
//  for (int i = -1; i < 2; i += 2)
//  {
//    for (int j = -1; j < 2; j += 2)
//    {
//      for (int k = -1; k < 2; k += 2)
//      {
//        CartVect corner( center );
//#ifdef MB_ORIENTED_BOX_UNIT_VECTORS
//        corner += i * box.length[0] * box.axis.col(0);
//        corner += j * box.length[1] * box.axis.col(1);
//        corner += k * box.length[2] * box.axis.col(2);
//#else
//        corner += i * box.axis.col(0);
//        corner += j * box.axis.col(1);
//        corner += k * box.axis.col(2);
//#endif
//        if (!contained( corner, tol ))
//          return false;
//      }
//    }
//  }
//  return true;
//}


/* This is a helper function to check limits on ray length, turning the box-ray
 * intersection test into a box-segment intersection test. Use this to test the
 * limits against one side (plane) of the box. The side of the box (plane) is
 * normal to an axis.
 *
 *   normal_par_pos  Coordinate of particle's position along axis normal to side of box
 *   normal_par_dir  Coordinate of particle's direction along axis normal to side of box
 *   half_extent     Distance between center of box and side of box
 *   nonneg_ray_len  Maximum ray length in positive direction (in front of origin)
 *   neg_ray_len     Maximum ray length in negative direction (behind origin)
 *   return          true if intersection with plane occurs within distance limits
 *
 * ray equation:   intersection = origin + dist*direction
 * plane equation: intersection.plane_normal = half_extent
 *
 * Assume plane_normal and direction are unit vectors. Combine equations.
 *
 *     (origin + dist*direction).plane_normal = half_extent
 *     origin.plane_normal + dist*direction.plane_normal = half_extent
 *     dist = (half_extent - origin.plane_normal)/(direction.plane_normal)
 *
 * Although this solves for distance, avoid floating point division.
 *
 *     dist*direction.plane_normal = half_extent - origin.plane_normal
 *
 * Use inequalities to test dist against ray length limits. Be aware that
 * inequalities change due to sign of direction.plane_normal.
 */
inline bool check_ray_limits(const double  normal_par_pos,
                             const double  normal_par_dir,
                             const double  half_extent,
                             const double* nonneg_ray_len,
                             const double* neg_ray_len ) {

  const double extent_pos_diff = half_extent - normal_par_pos;

  // limit in positive direction
  if(nonneg_ray_len) { // should be 0 <= t <= nonneg_ray_len
    assert(0 <= *nonneg_ray_len);
    if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
      if(*nonneg_ray_len*normal_par_dir>=extent_pos_diff && extent_pos_diff>=0) return true;
    } else if(normal_par_dir<0) {
      if(*nonneg_ray_len*normal_par_dir<=extent_pos_diff && extent_pos_diff<=0) return true;
    }
  } else {            // should be 0 <= t
    if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
      if(extent_pos_diff>=0) return true;
    } else if(normal_par_dir<0) {
      if(extent_pos_diff<=0) return true;
    }
  }

  // limit in negative direction
  if(neg_ray_len) {   // should be neg_ray_len <= t < 0
    assert(0 >= *neg_ray_len);
    if       (normal_par_dir>0) { // if/else if needed for pos/neg divisor
      if(*neg_ray_len*normal_par_dir<=extent_pos_diff && extent_pos_diff<0) return true;
    } else if(normal_par_dir<0) {
      if(*neg_ray_len*normal_par_dir>=extent_pos_diff && extent_pos_diff>0) return true;
    }
  }

  return false;
}

/* This implementation copied from cgmMC (overlap.C).
 * Original author:  Tim Tautges?
 */
bool OrientedBox::intersect_ray( const CartVect& ray_origin,
                                 const CartVect& ray_direction,
				 const double    reps,
				 const double*   nonneg_ray_len,
                                 const double*   neg_ray_len ) const
{
  // test distance from box center to line
  const CartVect cx       = center - ray_origin;
  const double dist_s     = cx % ray_direction;
  const double dist_sq    = cx % cx - (dist_s*dist_s);
  const double max_diagsq = outer_radius_squared(reps);

  // For the largest sphere, no intersections exist if discriminant is negative.
  // Geometrically, if distance from box center to line is greater than the
  // longest diagonal, there is no intersection.
  // manipulate the discriminant: 0 > dist_s*dist_s - cx%cx + max_diagsq
  if(dist_sq > max_diagsq) return false;

  // If the closest possible intersection must be closer than nonneg_ray_len. Be
  // careful with absolute value, squaring distances, and subtracting squared
  // distances.
  if (nonneg_ray_len) {
    assert(0<=*nonneg_ray_len);
    double max_len;
    if(neg_ray_len) {
      assert(0>=*neg_ray_len);
      max_len = std::max(*nonneg_ray_len,-(*neg_ray_len));
    } else {
      max_len = *nonneg_ray_len;
    }
    const double temp = fabs(dist_s) - max_len;
    if(0.0<temp && temp*temp>max_diagsq) return false;
  }

  // if smaller than shortest diagonal, we do hit
  if (dist_sq < inner_radius_squared(reps)) {
    // nonnegative direction
    if(dist_s>=0.0 ) {
      if(nonneg_ray_len) {
        if(*nonneg_ray_len>dist_s) return true;
      } else {
        return true;
      }
    // negative direction
    } else {
      if(neg_ray_len && *neg_ray_len<dist_s) return true;
    }
  }

  // get transpose of axes
  Matrix3 B = axes.transpose();

  // transform ray to box coordintae system
  CartVect par_pos = B * (ray_origin - center);
  CartVect par_dir = B * ray_direction;

  // Fast Rejection Test: Ray will not intersect if it is going away from the box.
  // This will not work for rays with neg_ray_len. length[0] is half of box width
  // along axes.col(0).
  const double half_x = length[0] + reps;
  const double half_y = length[1] + reps;
  const double half_z = length[2] + reps;
  if(!neg_ray_len) {
    if ((par_pos[0] >  half_x && par_dir[0] >= 0) ||
	(par_pos[0] < -half_x && par_dir[0] <= 0))
      return false;

    if ((par_pos[1] >  half_y && par_dir[1] >= 0) ||
	(par_pos[1] < -half_y && par_dir[1] <= 0))
      return false;

    if ((par_pos[2] >  half_z && par_dir[2] >= 0) ||
	(par_pos[2] < -half_z && par_dir[2] <= 0))
      return false;
  }

  // test if ray_origin is inside box
  if (par_pos[0] <= half_x && par_pos[0] >= -half_x &&
      par_pos[1] <= half_y && par_pos[1] >= -half_y &&
      par_pos[2] <= half_z && par_pos[2] >= -half_z)
    return true;

    //test two xy plane
  if (fabs(par_dir[0] * (half_z - par_pos[2]) + par_dir[2] * par_pos[0])
        <= fabs(par_dir[2] * half_x) &&  // test against x extents using z
      fabs(par_dir[1] * (half_z - par_pos[2]) + par_dir[2] * par_pos[1])
        <= fabs(par_dir[2] * half_y) &&  // test against y extents using z
      check_ray_limits( par_pos[2], par_dir[2], half_z, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[0] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[0])
        <= fabs(par_dir[2] * half_x) &&
      fabs(par_dir[1] * (-half_z - par_pos[2]) + par_dir[2] * par_pos[1])
        <= fabs(par_dir[2] * half_y) &&
      check_ray_limits( par_pos[2], par_dir[2], -half_z, nonneg_ray_len, neg_ray_len ) )
    return true;

    //test two xz plane
  if (fabs(par_dir[0] * (half_y - par_pos[1]) + par_dir[1] * par_pos[0])
        <= fabs(par_dir[1] * half_x) &&
      fabs(par_dir[2] * (half_y - par_pos[1]) + par_dir[1] * par_pos[2])
        <= fabs(par_dir[1] * half_z) &&
      check_ray_limits( par_pos[1], par_dir[1], half_y, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[0] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[0])
        <= fabs(par_dir[1] * half_x) &&
      fabs(par_dir[2] * (-half_y - par_pos[1]) + par_dir[1] * par_pos[2])
        <= fabs(par_dir[1] * half_z) &&
      check_ray_limits( par_pos[1], par_dir[1], -half_y, nonneg_ray_len, neg_ray_len ) )
    return true;

    //test two yz plane
  if (fabs(par_dir[1] * (half_x - par_pos[0]) + par_dir[0] * par_pos[1])
        <= fabs(par_dir[0] * half_y) &&
      fabs(par_dir[2] * (half_x - par_pos[0]) + par_dir[0] * par_pos[2])
        <= fabs(par_dir[0] * half_z) &&
      check_ray_limits( par_pos[0], par_dir[0], half_x, nonneg_ray_len, neg_ray_len ) )
    return true;
  if (fabs(par_dir[1] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[1])
        <= fabs(par_dir[0] * half_y) &&
      fabs(par_dir[2] * (-half_x - par_pos[0]) + par_dir[0] * par_pos[2])
        <= fabs(par_dir[0] * half_z) &&
      check_ray_limits( par_pos[0], par_dir[0], -half_x, nonneg_ray_len, neg_ray_len ) )
    return true;

  return false;
}

ErrorCode OrientedBox::make_hex( EntityHandle& hex, Interface* instance )
{
  ErrorCode rval;
  int signs[8][3] = { { -1, -1, -1 },
                      {  1, -1, -1 },
                      {  1,  1, -1 },
                      { -1,  1, -1 },
                      { -1, -1,  1 },
                      {  1, -1,  1 },
                      {  1,  1,  1 },
                      { -1,  1,  1 } };

  std::vector<EntityHandle> vertices;
  for (int i = 0; i < 8; ++i)
  {
    CartVect coords(center);
    for (int j = 0; j < 3; ++j){
#if MB_ORIENTED_BOX_UNIT_VECTORS
      coords += signs[i][j] * (axes.col(j)*length[j]);
#else
      coords += signs[i][j] * axes.col(j);
#endif
    }
    EntityHandle handle;
    rval = instance->create_vertex( coords.array(), handle );
    if (MB_SUCCESS != rval) {
      instance->delete_entities( &vertices[0], vertices.size() );
      return rval;
    }
    vertices.push_back( handle );
  }

  rval = instance->create_element( MBHEX, &vertices[0], vertices.size(), hex );
  if (MB_SUCCESS != rval) {
    instance->delete_entities( &vertices[0], vertices.size() );
    return rval;
  }

  return MB_SUCCESS;
}

void OrientedBox::closest_location_in_box(
                                    const CartVect& input_position,
                                    CartVect& output_position ) const
{
    // get coordinates on box axes
  const CartVect from_center = input_position - center;

#if MB_ORIENTED_BOX_UNIT_VECTORS
  CartVect local( from_center % axes.col(0),
                    from_center % axes.col(1),
                    from_center % axes.col(2) );

  for (int i = 0; i < 3; ++i) {
    if (local[i] < -length[i])
      local[i] = -length[i];
    else if (local[i] > length[i])
      local[i] =  length[i];
  }
#else
  CartVect local( (from_center % axes.col(0)) / (axes.col(0) % axes.col(0)),
                    (from_center % axes.col(1)) / (axes.col(1) % axes.col(1)),
                    (from_center % axes.col(2)) / (axes.col(2) % axes.col(2)) );

  for (int i = 0; i < 3; ++i) {
    if (local[i] < -1.0)
      local[i] = -1.0;
    else if (local[i] > 1.0)
      local[i] = 1.0;
  }
#endif

  output_position = center
                  + local[0] * axes.col(0)
                  + local[1] * axes.col(1)
                  + local[2] * axes.col(2);
}

} // namespace moab