Common/Geom/CurveAlg.h

Go to the documentation of this file.

//-------------------------------------------------------------------
// Copyright (c) 2007 Celeritive Technologies, Inc.
//
// 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.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
//-------------------------------------------------------------------
#pragma once
#ifndef VOLUMILL_CURVEALG_H
#define VOLUMILL_CURVEALG_H

#include <Geom/Curve.h>
#include <Util/TypeStripper.h>
#include <boost/variant.hpp>

namespace geom
{

template<class P> class BoundingBox;
typedef boost::variant<geom::Point2d, CurvePtr2d> CurveIntersectionResult;
typedef std::vector<CurveIntersectionResult> CurveIntersectionResults;

namespace detail
{
   void getBoundingBoxImpl (const geom::ParametricCurve<2>& c, BoundingBox<geom::Point2d>* pBox);
   void intersectCurvesImpl (
      const geom::ParametricCurve<2>& c1,
      const geom::ParametricCurve<2>& c2,
      double tol,
      CurveIntersectionResults* pIntersections);
   bool projectPointToCurveImpl (
      const geom::ParametricCurve<2>& curve, const Point2d& pt, double tol, Point2d* pProjPt=0, double* pParam=0);
   void tessellateCurveImpl (
      const geom::ParametricCurve<2>& curve,
      double tol,
      std::vector<Point2d>* pPoints);
   double getLengthImpl (
      const geom::ParametricCurve<2>& curve);
}

template<class Curve>
inline void getBoundingBox (const Curve& c, BoundingBox<geom::Point2d>* pBox)
{
   detail::getBoundingBoxImpl (util::TypeStripper<Curve>::cstrip (c), pBox);
}

template<class Curve>
inline void intersectCurves (const Curve& c1, const Curve& c2, double tol, CurveIntersectionResults* pIntersections)
{
   detail::intersectCurvesImpl (
      util::TypeStripper<Curve>::cstrip (c1),
      util::TypeStripper<Curve>::cstrip (c2),
      tol,
      pIntersections);
}

template<class Curve>
inline double getLength (const Curve& curve)
{
   return detail::getLengthImpl (util::TypeStripper<Curve>::cstrip (curve));
}

CurvePtr2d createSubCurve (const CurvePtr2d& curve, double startParam, double endParam);

template<class Curve, class Point>
inline bool projectPointToCurve (
   const Curve& curve, const Point& pt, double tol, Point* pProjPt=0, double* pParam=0)
{
   return detail::projectPointToCurveImpl (
      util::TypeStripper<Curve>::cstrip (curve), pt, tol, pProjPt, pParam);
}

template<class Curve, class Point>
inline void tessellateCurve (const Curve& curve, double tol, std::vector<Point>* pPoints)
{
   detail::tessellateCurveImpl (
      util::TypeStripper<Curve>::cstrip (curve), tol, pPoints);
}

template<class CurveIter, class Point>
inline void tessellateLoop (CurveIter begin, CurveIter end, double tol, std::vector<Point>* pPoints)
{
   for (; begin != end; ++begin)
   {
      if (pPoints->size())
         pPoints->pop_back();
      tessellateCurve (*begin, tol, pPoints);
   }

   if (pPoints->size() >= 2)
   {
      // Ensure the loop really is closed if it's nearly closed
      if (geom::distance (pPoints->front(), pPoints->back()) < util::pointResolution)
         pPoints->back() = pPoints->front();
   }
   else
      WARN (pPoints->size() >= 2);
}

template<class Curves, class Point>
inline void tessellateLoop (const Curves& curves, double tol, std::vector<Point>* pPoints)
{
   tessellateLoop (curves.begin(), curves.end(), tol, pPoints);
}

template<class Curves>
inline void reverseLoop (Curves& curves)
{
   std::reverse (curves.begin(), curves.end());
   for (typename Curves::iterator it = curves.begin(); it != curves.end(); ++it)
   {
      util::TypeStripper<typename Curves::value_type>::strip (*it).reverse();
   }
}

template<class CurveIter, class KnotIter>
inline void createKnots (CurveIter begin, CurveIter end, KnotIter outIt)
{
   double lastKnot = 0.0;
   *outIt++ = lastKnot;
   for (; begin != end; ++begin)
   {
      double knot = lastKnot + getLength (*begin);
      *outIt++ = knot;
      lastKnot = knot;
   }
}

template<class KnotIter>
int getSpan (KnotIter begin, KnotIter end, double param)
{
   KnotIter it = std::upper_bound (begin, end, param);
   int idx = (int) (it - begin);
   if (it == end)
      idx -= 2;
   else if (idx != 0)
      --idx;
   return idx;
}

template<class Arc, class P>
bool projectPointToArc (const Arc& arc, const P& pt, double tol, P* pProjPt = 0, double* pParam = 0)
{
   const P& center = arc.getCenter();
   VecNd<P::Dim> dir (pt - center);
   double dirLen = length (dir);
   if (length (dir) < util::doubleEps)
   {
      WARN (length (dir) >= util::doubleEps)("Cannot project the center point of an arc to the arc");
      util::assignPtr (pParam, 1e30);
      return false;
   }
   double r = arc.getRadius();
   // Scale the tolerance so that it is valid in the unitless world of angles
   tol /= r;
   const std::pair<double, double>& angleRange = arc.getAngleRange();
   VecNd<P::Dim> startDir;
   bool cw = angleRange.first > angleRange.second;
   double totalAngle;
   if (cw)
   {
      startDir[0] = cos (angleRange.second);
      startDir[1] = sin (angleRange.second);
      totalAngle = angleRange.first - angleRange.second;
   }
   else
   {
      startDir[0] = cos (angleRange.first);
      startDir[1] = sin (angleRange.first);
      totalAngle = angleRange.second - angleRange.first;
   }
   double cosAngle = startDir * dir;
   double sinAngle = dir[1] * startDir[0] - dir[0] * startDir[1];
   double angle = atan2 (sinAngle, cosAngle);
   if (angle < -tol)
      angle += 2 * util::pi;
   if (angle < tol)
   {
      util::assignPtr (pProjPt, cw ? getEndPoint (&arc) : getStartPoint (&arc));
      util::assignPtr (pParam, cw ? 1.0 : 0.0);
      return true;
   }
   if (fabs (angle - totalAngle) < tol)
   {
      util::assignPtr (pProjPt, cw ? getStartPoint (&arc) : getEndPoint (&arc));
      util::assignPtr (pParam, cw ? 0.0 : 1.0);
      return true;
   }
   if (angle > totalAngle)
   {
      P startPt (getStartPoint (&arc));
      P endPt (getEndPoint (&arc));
      double startDist = geom::distance (pt, startPt);
      double endDist = geom::distance (pt, endPt);
      if (startDist < endDist)
      {
         util::assignPtr (pProjPt, startPt);
         util::assignPtr (pParam, 0.0);
      }
      else
      {
         util::assignPtr (pProjPt, endPt);
         util::assignPtr (pParam, 1.0);
      }
      return false;
   }
   util::assignPtr (pProjPt, center + (r / dirLen) * dir);
   if (cw)
      util::assignPtr (pParam, 1 - angle / totalAngle);
   else
      util::assignPtr (pParam, angle / totalAngle);
   VALIDATE (angle > 0 && angle < totalAngle);
   return true;
}
}

#endif

---