OrientedBox.cpp

Go to the documentation of this file.

/*
 * 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 ([email protected])
 *\date 2006-07-18
 */

#include "moab/Interface.hpp"
#include "moab/CN.hpp"
#include "moab/OrientedBox.hpp"
#include "moab/Range.hpp"
#include <ostream>
#include <cassert>
#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