ClippingAndIntersection.hpp

Go to the documentation of this file.

// Copyright (c) 2012-2019 University of Lyon and CNRS (France).
// All rights reserved.
//
// This file is part of MEPP2; 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 3 of
// the License, or (at your option) any later version.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
#pragma once
 
#include "FEVV/Tools/Math/MatrixOperations.hpp"
 
namespace FEVV {
namespace Operators {
 
namespace Geometry {
 
template< typename GeometryTraits >
static bool
sphere_clip_vector(
    const typename GeometryTraits::Point &center, 
    double                               r,                                     
    const typename GeometryTraits::Point &p,
    typename GeometryTraits::Vector      &v,
    const GeometryTraits &gt)
{
  typedef typename GeometryTraits::Vector Vector;
  if(r < 0)
    r = -r; // the radius cannot be negative
 
  Vector w = gt.sub_p(p, center); // direction towards P from the sphere center
  double a = gt.dot_product(v, v);
 
  if(fabs(a) < std::numeric_limits< double >::epsilon())
    return false;
 
  double b = 2.0 * gt.dot_product(v, w);
  double c = gt.dot_product(w, w) - r * r;
  double delta = b * b - 4. * a * c;
  if(delta < 0.)
  {
    // Should not happen when p is inside the sphere, but happens sometimes 
  // (numerical precision)
    return true;
  }
 
  double t1 = (-b + ::sqrt(delta)) /
              (2.0 * a) /*, t2 = (-b - ::sqrt(delta)) / (2.0 * a)*/;
 
  if(t1 >= 1.)
  {
    // Inside the sphere
    return false;
  }
 
  // if (t1 < 0.) {
  // Should not happen, but happens sometimes (numerical precision or P not
  // inside the sphere)
  //  return true;
  //}
 
  v = gt.scalar_mult(v, t1);
 
  return true;
}
} // namespace Geometry
 
} // namespace Operators
} // namespace FEVV