Code

update 2geom
authorjohanengelen <johanengelen@users.sourceforge.net>
Sat, 14 Jun 2008 15:01:19 +0000 (15:01 +0000)
committerjohanengelen <johanengelen@users.sourceforge.net>
Sat, 14 Jun 2008 15:01:19 +0000 (15:01 +0000)
20 files changed:
src/2geom/convex-cover.cpp
src/2geom/curve.h
src/2geom/elliptical-arc.h
src/2geom/hvlinesegment.h
src/2geom/nearest-point.cpp
src/2geom/nearest-point.h
src/2geom/path.cpp
src/2geom/path.h
src/2geom/pathvector.cpp
src/2geom/pathvector.h
src/2geom/piecewise.cpp
src/2geom/piecewise.h
src/2geom/poly-laguerre-solve.cpp
src/2geom/poly.cpp
src/2geom/quadtree.cpp
src/2geom/quadtree.h
src/2geom/sbasis-geometric.cpp
src/2geom/sbasis-geometric.h
src/2geom/sbasis.h
src/2geom/solve-bezier-one-d.cpp

index 2d11b74190226f6d6361cc0fae4e450267eb0418..50c43e6d144de46e67904a21b7cf27156533dca5 100644 (file)
@@ -114,7 +114,7 @@ ConvexHull::angle_sort() {
 void
 ConvexHull::graham_scan() {
     unsigned stac = 2;
-    for(int i = 2; i < boundary.size(); i++) {
+    for(unsigned int i = 2; i < boundary.size(); i++) {
         double o = SignedTriangleArea(boundary[stac-2], 
                                       boundary[stac-1], 
                                       boundary[i]);
@@ -350,6 +350,7 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) {
     // al and bl now point to the top of the first pair of edges that overlap in y value
     double sweep_y = std::min(a.boundary[al][Y],
                               b.boundary[bl][Y]);
+    return ret;
 }
 
 /*** ConvexHull intersection(ConvexHull a, ConvexHull b);
@@ -359,10 +360,11 @@ ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) {
  */
 ConvexHull intersection(ConvexHull a, ConvexHull b) {
     ConvexHull ret;
+    /*
     int ai = 0, bi = 0;
     int aj = a.boundary.size() - 1;
     int bj = b.boundary.size() - 1;
-    
+    */
     /*while (true) {
         if(a[ai]
     }*/
index b76c0d328cdc4fb6defe90142c23a47b816b8c57..d201d5d5e0e1d1aecd8419652117e37e96355ee3 100644 (file)
@@ -105,10 +105,16 @@ public:
 
   virtual Point pointAt(Coord t) const { return pointAndDerivatives(t, 1).front(); }
   virtual Coord valueAt(Coord t, Dim2 d) const { return pointAt(t)[d]; }
+  virtual Point operator() (double t)  const { return pointAt(t); }
   virtual std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const = 0;
   virtual D2<SBasis> toSBasis() const = 0;
 };
 
+inline
+Coord nearest_point(Point const& p, Curve const& c)
+{
+       return c.nearestPoint(p);
+}
 
 } // end namespace Geom
 
index f3d4f49312f519b1312cf5d7db2715c022478b4e..3a3b90670199cdcda2db4bab342799d532cf280f 100644 (file)
@@ -166,7 +166,7 @@ class EllipticalArc : public Curve
         return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );
     }
     
-    // TODO: native implementation of the following methods
+    
     Rect boundsFast() const
     {
        return boundsExact();
@@ -174,6 +174,7 @@ class EllipticalArc : public Curve
   
     Rect boundsExact() const;
     
+    // TODO: native implementation of the following methods
     Rect boundsLocal(Interval i, unsigned int deg) const
     {
        return SBasisCurve(toSBasis()).boundsLocal(i, deg);
@@ -193,6 +194,7 @@ class EllipticalArc : public Curve
        return allNearestPoints(p, from, to).front();
     }
     
+    // TODO: native implementation of the following methods
     int winding(Point p) const
     {
        return SBasisCurve(toSBasis()).winding(p);
@@ -200,6 +202,7 @@ class EllipticalArc : public Curve
     
     Curve *derivative() const;
     
+    // TODO: native implementation of the following methods
     Curve *transformed(Matrix const &m) const
     {
        return SBasisCurve(toSBasis()).transformed(m);
index 8c479858bd48ab39ee82aec4455892a188264322..b3c20037897aece708fcd63014d61d7d20251ef2 100644 (file)
@@ -170,6 +170,11 @@ class HLineSegment : public Curve
                if ( from > to ) std::swap(from, to);
                double xfrom = pointAt(from)[X];
                double xto = pointAt(to)[X];
+               if ( xfrom > xto )
+               {
+                   std::swap(xfrom, xto);
+                   std::swap(from, to);
+               }
                if ( p[X] > xfrom && p[X] < xto )
                {
                        return (p[X] - initialPoint()[X]) / (finalPoint()[X] - initialPoint()[X]);
@@ -395,6 +400,11 @@ class VLineSegment : public Curve
                if ( from > to ) std::swap(from, to);
                double yfrom = pointAt(from)[Y];
                double yto = pointAt(to)[Y];
+               if (yfrom > yto)
+               {
+                   std::swap(yfrom, yto);
+                   std::swap(from, to);
+               }
                if ( p[Y] > yfrom && p[Y] < yto )
                {
                        return (p[Y] - initialPoint()[Y]) / (finalPoint()[Y] - initialPoint()[Y]);
index de880713f651bcc8feb710c1bc3d9f364c5f4c83..074de1cb3364c54780731ffb0588c23e33eb9be3 100644 (file)
-/*
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
- *
- * Authors:
- *             
- *             Marco Cecchetti <mrcekets at gmail.com>
- * 
- * Copyright 2007-2008  authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-#include "nearest-point.h"
-
-namespace Geom
-{
-
-////////////////////////////////////////////////////////////////////////////////
-// D2<SBasis> versions
-
-/*
- * Return the parameter t of a nearest point on the portion of the curve "c", 
- * related to the interval [from, to], to the point "p".
- * The needed curve derivative "dc" is passed as parameter.
- * The function return the first nearest point to "p" that is found.
- */
-
-double nearest_point( Point const& p, 
-                                         D2<SBasis> const& c, 
-                                         D2<SBasis> const& dc, 
-                             double from, double to )
-{
-       if ( from > to ) std::swap(from, to);
-       if ( from < 0 || to > 1 )
-       {
-               THROW_RANGEERROR("[from,to] interval out of bounds");
-       }
-
-       SBasis dd = dot(c - p, dc);     
-       std::vector<double> zeros = Geom::roots(dd);
-       
-       double closest = from;
-       double min_dist_sq = L2sq(c(from) - p);
-       double distsq;
-       for ( unsigned int i = 0; i < zeros.size(); ++i )
-       {
-               distsq = L2sq(c(zeros[i]) - p);
-               if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
-               {
-                       closest = zeros[i];
-                       min_dist_sq = distsq;
-               }
-       }
-       if ( min_dist_sq > L2sq( c(to) - p ) )
-               closest = to;
-       return closest;
-
-}
-
-/*
- * Return the parameters t of all the nearest points on the portion of 
- * the curve "c", related to the interval [from, to], to the point "p".
- * The needed curve derivative "dc" is passed as parameter.
- */
-
-std::vector<double> 
-all_nearest_points( Point const& p, 
-                           D2<SBasis> const& c, 
-                           D2<SBasis> const& dc, 
-                           double from, double to )
-{
-       std::swap(from, to);
-       if ( from > to ) std::swap(from, to);
-       if ( from < 0 || to > 1 )
-       {
-               THROW_RANGEERROR("[from,to] interval out of bounds");
-       }
-
-       std::vector<double> result;
-       SBasis dd = dot(c - p, Geom::derivative(c));
-       
-       std::vector<double> zeros = Geom::roots(dd);
-       std::vector<double> candidates;
-       candidates.push_back(from);
-       candidates.insert(candidates.end(), zeros.begin(), zeros.end());
-       candidates.push_back(to);
-       std::vector<double> distsq;
-       distsq.reserve(candidates.size());
-       for ( unsigned int i = 0; i < candidates.size(); ++i )
-       {
-               distsq.push_back( L2sq(c(candidates[i]) - p) );
-       }
-       unsigned int closest = 0;
-       double dsq = distsq[0];
-       for ( unsigned int i = 1; i < candidates.size(); ++i )
-       {
-               if ( dsq > distsq[i] )
-               {
-                       closest = i;
-                       dsq = distsq[i];
-               }
-       }
-       for ( unsigned int i = 0; i < candidates.size(); ++i )
-       {
-               if( distsq[closest] == distsq[i] )
-               {
-                       result.push_back(candidates[i]);
-               }
-       }
-       return result;
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// Piecewise< D2<SBasis> > versions
-
-
-double nearest_point( Point const& p,  
-                                         Piecewise< D2<SBasis> > const& c,
-                             double from, double to )
-{
-       if ( from > to ) std::swap(from, to);
-       if ( from < c.cuts[0] || to > c.cuts[c.size()] )
-       {
-               THROW_RANGEERROR("[from,to] interval out of bounds");
-       }
-       
-       unsigned int si = c.segN(from);
-       unsigned int ei = c.segN(to);
-       if ( si == ei )
-       {
-               double nearest=
-                       nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
-               return c.mapToDomain(nearest, si);
-       }
-       double t;
-       double nearest = nearest_point(p, c[si], c.segT(from, si));
-       unsigned int ni = si;
-       double dsq;
-       double mindistsq = distanceSq(p, c[si](nearest));
-       Rect bb;
-       for ( unsigned int i = si + 1; i < ei; ++i )
-       {
-               bb = bounds_fast(c[i]);
-               dsq = distanceSq(p, bb);
-               if ( mindistsq <= dsq ) continue;
-               t = nearest_point(p, c[i]);
-               dsq = distanceSq(p, c[i](t));
-               if ( mindistsq > dsq )
-               {
-                       nearest = t;
-                       ni = i;
-                       mindistsq = dsq;
-               }
-       }
-       bb = bounds_fast(c[ei]);
-       dsq = distanceSq(p, bb);
-       if ( mindistsq > dsq )
-       {
-               t = nearest_point(p, c[ei], 0, c.segT(to, ei));
-               dsq = distanceSq(p, c[ei](t));
-               if ( mindistsq > dsq )
-               {
-                       nearest = t;
-                       ni = ei;
-               }
-       }
-       return c.mapToDomain(nearest, ni);
-}
-
-std::vector<double> 
-all_nearest_points( Point const& p, 
-                    Piecewise< D2<SBasis> > const& c, 
-                           double from, double to )
-{
-       if ( from > to ) std::swap(from, to);
-       if ( from < c.cuts[0] || to > c.cuts[c.size()] )
-       {
-               THROW_RANGEERROR("[from,to] interval out of bounds");
-       }
-       
-       unsigned int si = c.segN(from);
-       unsigned int ei = c.segN(to);
-       if ( si == ei )
-       {
-               std::vector<double>     all_nearest = 
-                       all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
-               for ( unsigned int i = 0; i < all_nearest.size(); ++i )
-               {
-                       all_nearest[i] = c.mapToDomain(all_nearest[i], si);
-               }
-               return all_nearest;
-       }
-       std::vector<double> all_t;
-       std::vector< std::vector<double> > all_np;
-       all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
-       std::vector<unsigned int> ni;
-       ni.push_back(si);
-       double dsq;
-       double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
-       Rect bb;
-       for ( unsigned int i = si + 1; i < ei; ++i )
-       {
-               bb = bounds_fast(c[i]);
-               dsq = distanceSq(p, bb);
-               if ( mindistsq < dsq ) continue;
-               all_t = all_nearest_points(p, c[i]);
-               dsq = distanceSq( p, c[i](all_t.front()) );
-               if ( mindistsq > dsq )
-               {
-                       all_np.clear();
-                       all_np.push_back(all_t);
-                       ni.clear();
-                       ni.push_back(i);
-                       mindistsq = dsq;
-               }
-               else if ( mindistsq == dsq )
-               {
-                       all_np.push_back(all_t);
-                       ni.push_back(i);
-               }
-       }
-       bb = bounds_fast(c[ei]);
-       dsq = distanceSq(p, bb);
-       if ( mindistsq >= dsq )
-       {
-               all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
-               dsq = distanceSq( p, c[ei](all_t.front()) );
-               if ( mindistsq > dsq )
-               {
-                       for ( unsigned int i = 0; i < all_t.size(); ++i )
-                       {
-                               all_t[i] = c.mapToDomain(all_t[i], ei);
-                       }
-                       return all_t;
-               }
-               else if ( mindistsq == dsq )
-               {
-                       all_np.push_back(all_t);
-                       ni.push_back(ei);
-               }
-       }
-       std::vector<double> all_nearest;
-       for ( unsigned int i = 0; i < all_np.size(); ++i )
-       {
-               for ( unsigned int j = 0; j < all_np[i].size(); ++j )
-               {
-                       all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );
-               }
-       }
-       return all_nearest;
-}
-
-} // end namespace Geom
-
-
+/*\r
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
+ *\r
+ * Authors:\r
+ *             \r
+ *             Marco Cecchetti <mrcekets at gmail.com>\r
+ * \r
+ * Copyright 2007-2008  authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#include "nearest-point.h"\r
+\r
+namespace Geom\r
+{\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// D2<SBasis> versions\r
+\r
+/*\r
+ * Return the parameter t of a nearest point on the portion of the curve "c", \r
+ * related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ * The function return the first nearest point to "p" that is found.\r
+ */\r
+\r
+double nearest_point( Point const& p, \r
+                                         D2<SBasis> const& c, \r
+                                         D2<SBasis> const& dc, \r
+                             double from, double to )\r
+{\r
+       if ( from > to ) std::swap(from, to);\r
+       if ( from < 0 || to > 1 )\r
+       {\r
+               THROW_RANGEERROR("[from,to] interval out of bounds");\r
+       }\r
+\r
+       SBasis dd = dot(c - p, dc);     \r
+       std::vector<double> zeros = Geom::roots(dd);\r
+       \r
+       double closest = from;\r
+       double min_dist_sq = L2sq(c(from) - p);\r
+       double distsq;\r
+       for ( unsigned int i = 0; i < zeros.size(); ++i )\r
+       {\r
+               distsq = L2sq(c(zeros[i]) - p);\r
+               if ( min_dist_sq > L2sq(c(zeros[i]) - p) )\r
+               {\r
+                       closest = zeros[i];\r
+                       min_dist_sq = distsq;\r
+               }\r
+       }\r
+       if ( min_dist_sq > L2sq( c(to) - p ) )\r
+               closest = to;\r
+       return closest;\r
+\r
+}\r
+\r
+/*\r
+ * Return the parameters t of all the nearest points on the portion of \r
+ * the curve "c", related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ */\r
+\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+                           D2<SBasis> const& c, \r
+                           D2<SBasis> const& dc, \r
+                           double from, double to )\r
+{\r
+       std::swap(from, to);\r
+       if ( from > to ) std::swap(from, to);\r
+       if ( from < 0 || to > 1 )\r
+       {\r
+               THROW_RANGEERROR("[from,to] interval out of bounds");\r
+       }\r
+\r
+       std::vector<double> result;\r
+       SBasis dd = dot(c - p, Geom::derivative(c));\r
+       \r
+       std::vector<double> zeros = Geom::roots(dd);\r
+       std::vector<double> candidates;\r
+       candidates.push_back(from);\r
+       candidates.insert(candidates.end(), zeros.begin(), zeros.end());\r
+       candidates.push_back(to);\r
+       std::vector<double> distsq;\r
+       distsq.reserve(candidates.size());\r
+       for ( unsigned int i = 0; i < candidates.size(); ++i )\r
+       {\r
+               distsq.push_back( L2sq(c(candidates[i]) - p) );\r
+       }\r
+       unsigned int closest = 0;\r
+       double dsq = distsq[0];\r
+       for ( unsigned int i = 1; i < candidates.size(); ++i )\r
+       {\r
+               if ( dsq > distsq[i] )\r
+               {\r
+                       closest = i;\r
+                       dsq = distsq[i];\r
+               }\r
+       }\r
+       for ( unsigned int i = 0; i < candidates.size(); ++i )\r
+       {\r
+               if( distsq[closest] == distsq[i] )\r
+               {\r
+                       result.push_back(candidates[i]);\r
+               }\r
+       }\r
+       return result;\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// Piecewise< D2<SBasis> > versions\r
+\r
+\r
+double nearest_point( Point const& p,  \r
+                                         Piecewise< D2<SBasis> > const& c,\r
+                             double from, double to )\r
+{\r
+       if ( from > to ) std::swap(from, to);\r
+       if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
+       {\r
+               THROW_RANGEERROR("[from,to] interval out of bounds");\r
+       }\r
+       \r
+       unsigned int si = c.segN(from);\r
+       unsigned int ei = c.segN(to);\r
+       if ( si == ei )\r
+       {\r
+               double nearest=\r
+                       nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));\r
+               return c.mapToDomain(nearest, si);\r
+       }\r
+       double t;\r
+       double nearest = nearest_point(p, c[si], c.segT(from, si));\r
+       unsigned int ni = si;\r
+       double dsq;\r
+       double mindistsq = distanceSq(p, c[si](nearest));\r
+       Rect bb;\r
+       for ( unsigned int i = si + 1; i < ei; ++i )\r
+       {\r
+               bb = bounds_fast(c[i]);\r
+               dsq = distanceSq(p, bb);\r
+               if ( mindistsq <= dsq ) continue;\r
+               t = nearest_point(p, c[i]);\r
+               dsq = distanceSq(p, c[i](t));\r
+               if ( mindistsq > dsq )\r
+               {\r
+                       nearest = t;\r
+                       ni = i;\r
+                       mindistsq = dsq;\r
+               }\r
+       }\r
+       bb = bounds_fast(c[ei]);\r
+       dsq = distanceSq(p, bb);\r
+       if ( mindistsq > dsq )\r
+       {\r
+               t = nearest_point(p, c[ei], 0, c.segT(to, ei));\r
+               dsq = distanceSq(p, c[ei](t));\r
+               if ( mindistsq > dsq )\r
+               {\r
+                       nearest = t;\r
+                       ni = ei;\r
+               }\r
+       }\r
+       return c.mapToDomain(nearest, ni);\r
+}\r
+\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+                    Piecewise< D2<SBasis> > const& c, \r
+                           double from, double to )\r
+{\r
+       if ( from > to ) std::swap(from, to);\r
+       if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
+       {\r
+               THROW_RANGEERROR("[from,to] interval out of bounds");\r
+       }\r
+       \r
+       unsigned int si = c.segN(from);\r
+       unsigned int ei = c.segN(to);\r
+       if ( si == ei )\r
+       {\r
+               std::vector<double>     all_nearest = \r
+                       all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));\r
+               for ( unsigned int i = 0; i < all_nearest.size(); ++i )\r
+               {\r
+                       all_nearest[i] = c.mapToDomain(all_nearest[i], si);\r
+               }\r
+               return all_nearest;\r
+       }\r
+       std::vector<double> all_t;\r
+       std::vector< std::vector<double> > all_np;\r
+       all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );\r
+       std::vector<unsigned int> ni;\r
+       ni.push_back(si);\r
+       double dsq;\r
+       double mindistsq = distanceSq( p, c[si](all_np.front().front()) );\r
+       Rect bb;\r
+       for ( unsigned int i = si + 1; i < ei; ++i )\r
+       {\r
+               bb = bounds_fast(c[i]);\r
+               dsq = distanceSq(p, bb);\r
+               if ( mindistsq < dsq ) continue;\r
+               all_t = all_nearest_points(p, c[i]);\r
+               dsq = distanceSq( p, c[i](all_t.front()) );\r
+               if ( mindistsq > dsq )\r
+               {\r
+                       all_np.clear();\r
+                       all_np.push_back(all_t);\r
+                       ni.clear();\r
+                       ni.push_back(i);\r
+                       mindistsq = dsq;\r
+               }\r
+               else if ( mindistsq == dsq )\r
+               {\r
+                       all_np.push_back(all_t);\r
+                       ni.push_back(i);\r
+               }\r
+       }\r
+       bb = bounds_fast(c[ei]);\r
+       dsq = distanceSq(p, bb);\r
+       if ( mindistsq >= dsq )\r
+       {\r
+               all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));\r
+               dsq = distanceSq( p, c[ei](all_t.front()) );\r
+               if ( mindistsq > dsq )\r
+               {\r
+                       for ( unsigned int i = 0; i < all_t.size(); ++i )\r
+                       {\r
+                               all_t[i] = c.mapToDomain(all_t[i], ei);\r
+                       }\r
+                       return all_t;\r
+               }\r
+               else if ( mindistsq == dsq )\r
+               {\r
+                       all_np.push_back(all_t);\r
+                       ni.push_back(ei);\r
+               }\r
+       }\r
+       std::vector<double> all_nearest;\r
+       for ( unsigned int i = 0; i < all_np.size(); ++i )\r
+       {\r
+               for ( unsigned int j = 0; j < all_np[i].size(); ++j )\r
+               {\r
+                       all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );\r
+               }\r
+       }\r
+       return all_nearest;\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
index 8c1a6f7d7be90dde4374ea5e63b2cb89f0ca2efd..73ac0c3ce73f9785a68d19d3afb9f83bdf310f79 100644 (file)
-/*
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
- *
- *
- * Authors:
- *             
- *             Marco Cecchetti <mrcekets at gmail.com>
- * 
- * Copyright 2007-2008  authors
- *
- * This library is free software; you can redistribute it and/or
- * modify it either under the terms of the GNU Lesser General Public
- * License version 2.1 as published by the Free Software Foundation
- * (the "LGPL") or, at your option, under the terms of the Mozilla
- * Public License Version 1.1 (the "MPL"). If you do not alter this
- * notice, a recipient may use your version of this file under either
- * the MPL or the LGPL.
- *
- * You should have received a copy of the LGPL along with this library
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- * You should have received a copy of the MPL along with this library
- * in the file COPYING-MPL-1.1
- *
- * The contents of this file are subject to the Mozilla Public License
- * Version 1.1 (the "License"); you may not use this file except in
- * compliance with the License. You may obtain a copy of the License at
- * http://www.mozilla.org/MPL/
- *
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for
- * the specific language governing rights and limitations.
- */
-
-
-#ifndef _NEAREST_POINT_H_
-#define _NEAREST_POINT_H_
-
-
-#include <vector>
-
-#include "d2.h"
-#include "piecewise.h"
-#include "exception.h"
-
-
-
-namespace Geom
-{
-
-/*
- * Given a line L specified by a point A and direction vector v,
- * return the point on L nearest to p. Note that the returned value
- * is with respect to the _normalized_ direction of v!
- */
-inline double nearest_point(Point const &p, Point const &A, Point const &v)
-{
-    Point d(p - A);
-    return d[0] * v[0] + d[1] * v[1];
-}
-
-////////////////////////////////////////////////////////////////////////////////
-// D2<SBasis> versions
-
-/*
- * Return the parameter t of a nearest point on the portion of the curve "c", 
- * related to the interval [from, to], to the point "p".
- * The needed curve derivative "dc" is passed as parameter.
- * The function return the first nearest point to "p" that is found.
- */
-double nearest_point( Point const& p,
-                             D2<SBasis> const& c, D2<SBasis> const& dc, 
-                             double from = 0, double to = 1 );
-
-inline
-double nearest_point( Point const& p, 
-                             D2<SBasis> const& c, 
-                             double from = 0, double to = 1 )
-{
-       return nearest_point(p, c, Geom::derivative(c), from, to);
-}
-
-/*
- * Return the parameters t of all the nearest points on the portion of 
- * the curve "c", related to the interval [from, to], to the point "p".
- * The needed curve derivative "dc" is passed as parameter.
- */
-std::vector<double> 
-all_nearest_points( Point const& p, 
-                           D2<SBasis> const& c, D2<SBasis> const& dc, 
-                           double from = 0, double to = 1 );
-
-inline
-std::vector<double> 
-all_nearest_points( Point const& p, 
-                           D2<SBasis> const& c, 
-                           double from = 0, double to = 1 )
-{
-       return all_nearest_points(p, c,  Geom::derivative(c), from, to);
-}
-
-
-////////////////////////////////////////////////////////////////////////////////
-// Piecewise< D2<SBasis> > versions
-
-double nearest_point( Point const& p, 
-                             Piecewise< D2<SBasis> > const& c, 
-                             double from, double to );
-
-inline
-double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) 
-{
-       return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);
-}
-
-
-std::vector<double> 
-all_nearest_points( Point const& p, 
-                                       Piecewise< D2<SBasis> > const& c, 
-                           double from, double to );
-
-inline
-std::vector<double> 
-all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) 
-{
-       return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);
-}
-
-} // end namespace Geom
-
-
-
-#endif /*_NEAREST_POINT_H_*/
+/*\r
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
+ *\r
+ *\r
+ * Authors:\r
+ *             \r
+ *             Marco Cecchetti <mrcekets at gmail.com>\r
+ * \r
+ * Copyright 2007-2008  authors\r
+ *\r
+ * This library is free software; you can redistribute it and/or\r
+ * modify it either under the terms of the GNU Lesser General Public\r
+ * License version 2.1 as published by the Free Software Foundation\r
+ * (the "LGPL") or, at your option, under the terms of the Mozilla\r
+ * Public License Version 1.1 (the "MPL"). If you do not alter this\r
+ * notice, a recipient may use your version of this file under either\r
+ * the MPL or the LGPL.\r
+ *\r
+ * You should have received a copy of the LGPL along with this library\r
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
+ * You should have received a copy of the MPL along with this library\r
+ * in the file COPYING-MPL-1.1\r
+ *\r
+ * The contents of this file are subject to the Mozilla Public License\r
+ * Version 1.1 (the "License"); you may not use this file except in\r
+ * compliance with the License. You may obtain a copy of the License at\r
+ * http://www.mozilla.org/MPL/\r
+ *\r
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
+ * the specific language governing rights and limitations.\r
+ */\r
+\r
+\r
+#ifndef _NEAREST_POINT_H_\r
+#define _NEAREST_POINT_H_\r
+\r
+\r
+#include <vector>\r
+\r
+#include "d2.h"\r
+#include "piecewise.h"\r
+#include "exception.h"\r
+\r
+\r
+\r
+namespace Geom\r
+{\r
+\r
+/*\r
+ * Given a line L specified by a point A and direction vector v,\r
+ * return the point on L nearest to p. Note that the returned value\r
+ * is with respect to the _normalized_ direction of v!\r
+ */\r
+inline double nearest_point(Point const &p, Point const &A, Point const &v)\r
+{\r
+    Point d(p - A);\r
+    return d[0] * v[0] + d[1] * v[1];\r
+}\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// D2<SBasis> versions\r
+\r
+/*\r
+ * Return the parameter t of a nearest point on the portion of the curve "c", \r
+ * related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ * The function return the first nearest point to "p" that is found.\r
+ */\r
+double nearest_point( Point const& p,\r
+                             D2<SBasis> const& c, D2<SBasis> const& dc, \r
+                             double from = 0, double to = 1 );\r
+\r
+inline\r
+double nearest_point( Point const& p, \r
+                             D2<SBasis> const& c, \r
+                             double from = 0, double to = 1 )\r
+{\r
+       return nearest_point(p, c, Geom::derivative(c), from, to);\r
+}\r
+\r
+/*\r
+ * Return the parameters t of all the nearest points on the portion of \r
+ * the curve "c", related to the interval [from, to], to the point "p".\r
+ * The needed curve derivative "dc" is passed as parameter.\r
+ */\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+                           D2<SBasis> const& c, D2<SBasis> const& dc, \r
+                           double from = 0, double to = 1 );\r
+\r
+inline\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+                           D2<SBasis> const& c, \r
+                           double from = 0, double to = 1 )\r
+{\r
+       return all_nearest_points(p, c,  Geom::derivative(c), from, to);\r
+}\r
+\r
+\r
+////////////////////////////////////////////////////////////////////////////////\r
+// Piecewise< D2<SBasis> > versions\r
+\r
+double nearest_point( Point const& p, \r
+                             Piecewise< D2<SBasis> > const& c, \r
+                             double from, double to );\r
+\r
+inline\r
+double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
+{\r
+       return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);\r
+}\r
+\r
+\r
+std::vector<double> \r
+all_nearest_points( Point const& p, \r
+                                       Piecewise< D2<SBasis> > const& c, \r
+                           double from, double to );\r
+\r
+inline\r
+std::vector<double> \r
+all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
+{\r
+       return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);\r
+}\r
+\r
+} // end namespace Geom\r
+\r
+\r
+\r
+#endif /*_NEAREST_POINT_H_*/\r
index 35c7b76bce80201f6ebd03f2e6e66bf94fc58c40..0ff1f3b0c65a624fec53fc4f644096be5b245aff 100644 (file)
@@ -50,8 +50,11 @@ void Path::swap(Path &other) {
 
 Rect Path::boundsFast() const {
   Rect bounds=front().boundsFast();
-  for ( const_iterator iter=++begin(); iter != end() ; ++iter ) {
-    bounds.unionWith(iter->boundsFast());
+  const_iterator iter = begin();
+  if ( iter != end() ) {
+         for ( ++iter; iter != end() ; ++iter ) {
+           bounds.unionWith(iter->boundsFast());
+         }
   }
   return bounds;
 }
@@ -236,6 +239,56 @@ double Path::nearestPoint(Point const& _point, double from, double to) const
        return ni + nearest;
 }
 
+Rect Path::boundsFast()
+{
+    Rect bound;
+    if (empty()) return bound;
+    
+    bound = begin()->boundsFast();
+    double top = bound.top();
+    double bottom = bound.bottom();
+    double left = bound.left();
+    double right = bound.right();
+    for (iterator it = ++begin(); it != end(); ++it)
+    {
+        bound = it->boundsFast();
+        if ( top > bound.top() )           top = bound.top();
+        if ( bottom < bound.bottom() )     bottom = bound.bottom();
+        if ( left > bound.left() )         left = bound.left();
+        if ( right < bound.right() )       right = bound.right();
+    }
+    bound[0][0] = left;
+    bound[0][1] = right;
+    bound[1][0] = top;
+    bound[1][1] = bottom;
+    return bound;
+}
+
+Rect Path::boundsExact()
+{
+    Rect bound;
+    if (empty()) return bound;
+    
+    bound = begin()->boundsExact();
+    double top = bound.top();
+    double bottom = bound.bottom();
+    double left = bound.left();
+    double right = bound.right();
+    for (iterator it = ++begin(); it != end(); ++it)
+    {
+        bound = it->boundsExact();
+        if ( top > bound.top() )           top = bound.top();
+        if ( bottom < bound.bottom() )     bottom = bound.bottom();
+        if ( left > bound.left() )         left = bound.left();
+        if ( right < bound.right() )       right = bound.right();
+    }
+    bound[0][0] = left;
+    bound[0][1] = right;
+    bound[1][0] = top;
+    bound[1][1] = bottom;
+    return bound;
+}
+
 //This assumes that you can't be perfect in your t-vals, and as such, tweaks the start
 void Path::appendPortionTo(Path &ret, double from, double to) const {
   assert(from >= 0 && to >= 0);
index 670d4c1251eeb0c23fe97195a49055712290ee60..8f51c46c20b7ac26e58655706171326a1fd81e06 100644 (file)
@@ -266,6 +266,12 @@ public:
          return (*this)[i].valueAt(lt, d);
   }
 
+  
+  Point operator() (double t) const
+  {
+         return pointAt(t);
+  }
+  
   std::vector<double> roots(double v, Dim2 d) const {
     std::vector<double> res;
     for(unsigned i = 0; i <= size(); i++) {
@@ -297,6 +303,9 @@ public:
          return nearestPoint(_point, 0, sz);
   }
    
+  Rect boundsFast();
+  Rect boundsExact();
+  
   void appendPortionTo(Path &p, double f, double t) const;
 
   Path portion(double f, double t) const {
@@ -543,6 +552,13 @@ inline static Piecewise<D2<SBasis> > paths_to_pw(std::vector<Path> paths) {
     return ret;
 }
 
+inline
+Coord nearest_point(Point const& p, Path const& c)
+{
+       return c.nearestPoint(p);
+}
+
+
 /*
 class PathPortion : public Curve {
   Path *source;
index 696ccad7752349554e250bbae1f50f6f376cabcf..9c8111767df1c37d6c000b34b6ecad183cdee583 100644 (file)
@@ -65,6 +65,59 @@ PathVector reverse_paths_and_order (PathVector const & path_in)
     return path_out;
 }
 
+Rect bounds_fast( PathVector const& pv )
+{
+    typedef PathVector::const_iterator const_iterator;
+    
+    Rect bound;
+    if (pv.empty()) return bound;
+    
+    bound = (pv.begin())->boundsFast();
+    double top = bound.top();
+    double bottom = bound.bottom();
+    double left = bound.left();
+    double right = bound.right();
+    for (const_iterator it = ++(pv.begin()); it != pv.end(); ++it)
+    {
+        bound = it->boundsFast();
+        if ( top > bound.top() )           top = bound.top();
+        if ( bottom < bound.bottom() )     bottom = bound.bottom();
+        if ( left > bound.left() )         left = bound.left();
+        if ( right < bound.right() )       right = bound.right();
+    }
+    bound[0][0] = left;
+    bound[0][1] = right;
+    bound[1][0] = top;
+    bound[1][1] = bottom;
+    return bound;
+}
+
+Rect bounds_exact( PathVector const& pv )
+{
+    typedef PathVector::const_iterator const_iterator;
+    
+    Rect bound;
+    if (pv.empty()) return bound;
+    
+    bound = (pv.begin())->boundsExact();
+    double top = bound.top();
+    double bottom = bound.bottom();
+    double left = bound.left();
+    double right = bound.right();
+    for (const_iterator it = ++(pv.begin()); it != pv.end(); ++it)
+    {
+        bound = it->boundsExact();
+        if ( top > bound.top() )           top = bound.top();
+        if ( bottom < bound.bottom() )     bottom = bound.bottom();
+        if ( left > bound.left() )         left = bound.left();
+        if ( right < bound.right() )       right = bound.right();
+    }
+    bound[0][0] = left;
+    bound[0][1] = right;
+    bound[1][0] = top;
+    bound[1][1] = bottom;
+    return bound;
+}
 
 } // namespace Geom
 
index 97a0d19f19ff33fc15a948d3e5eaf2adef25fca3..a9f4d88ed328b86bea3043e57829768182032498 100644 (file)
@@ -46,6 +46,8 @@ PathVector operator* (PathVector const & path_in, Matrix const &m);
 
 PathVector reverse_paths_and_order (PathVector const & path_in);
 
+Rect bounds_fast( PathVector const & pv );
+Rect bounds_exact( PathVector const & pv );
 }
 
 #endif // SEEN_GEOM_PATHVECTOR_H
index dc91ab4a93fcd047417dabb6197689f68c3a5763..222152aa36fd8209c39f7f9fb2bbc56ca22a9ca1 100644 (file)
@@ -167,6 +167,22 @@ std::vector<double> roots(Piecewise<SBasis> const &f){
     return result;
 }
 
+Piecewise<SBasis> reverse(Piecewise<SBasis> const &f) {
+    Piecewise<SBasis> ret = Piecewise<SBasis>();
+    ret.cuts.resize(f.cuts.size());
+    ret.segs.resize(f.segs.size());
+    double start = f.cuts[0];
+    double end = f.cuts.back();
+    for (unsigned i = 0; i < f.cuts.size(); i++) {
+        double x = f.cuts[f.cuts.size() - 1 - i];
+        ret.cuts[i] = end - (x - start);
+    }
+    for (unsigned i = 0; i < f.segs.size(); i++)
+        ret.segs[i] = reverse(f[f.segs.size() - i - 1]);
+    return ret;
+}
+
+
 }
 /*
   Local Variables:
index 574d6de62f29011eab0e63f5e53e5d8140669c7f..3f13d15b01ef9215ba7b16262d6bbce8c70bf863 100644 (file)
@@ -704,6 +704,7 @@ Piecewise<T> derivative(Piecewise<T> const &a) {
 }
 
 std::vector<double> roots(Piecewise<SBasis> const &f);
+Piecewise<SBasis> reverse(Piecewise<SBasis> const &f);
 
 }
 
index 766f16edafa1e2621e9d7e24bfa0c3da9cb2f3a1..266e24c2b93c690babeb0f12891fb3f4fbd76f81 100644 (file)
@@ -12,7 +12,7 @@ cdouble laguerre_internal_complex(Poly const & p,
     double n = p.degree();
     quad_root = false;
     const unsigned shuffle_rate = 10;
-    static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
+    //static double shuffle[] = {0, 0.5, 0.25, 0.75, 0.125, 0.375, 0.625, 0.875, 1.0};
     unsigned shuffle_counter = 0;
     while(std::norm(a) > (tol*tol)) {
         //std::cout << "xk = " << xk << std::endl;
@@ -132,7 +132,12 @@ laguerre(Poly p, const double tol) {
 std::vector<double>
 laguerre_real_interval(Poly const & ply,
                        const double lo, const double hi,
-                       const double tol) {
+                       const double tol) 
+{
+    /* not implemented*/
+    assert(false);
+    std::vector<double> result;
+    return result;
 }
 
 /*
index 4f5ba6c557ea25227c6672c4dd9327ad2a200447..27f98596b531eb7820ab8091e284cd851ea84837 100644 (file)
@@ -45,7 +45,7 @@ std::vector<std::complex<double> > solve(Poly const & pp) {
        
     gsl_complex_packed_ptr z = new double[p.degree()*2];
     double* a = new double[p.size()];
-    for(int i = 0; i < p.size(); i++)
+    for(unsigned int i = 0; i < p.size(); i++)
         a[i] = p[i];
     std::vector<std::complex<double> > roots;
     //roots.resize(p.degree());
@@ -55,7 +55,7 @@ std::vector<std::complex<double> > solve(Poly const & pp) {
      
     gsl_poly_complex_workspace_free (w);
      
-    for (int i = 0; i < p.degree(); i++) {
+    for (unsigned int i = 0; i < p.degree(); i++) {
         roots.push_back(std::complex<double> (z[2*i] ,z[2*i+1]));
         //printf ("z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]);
     }    
@@ -67,7 +67,7 @@ std::vector<double > solve_reals(Poly const & p) {
     std::vector<std::complex<double> > roots = solve(p);
     std::vector<double> real_roots;
     
-    for(int i = 0; i < roots.size(); i++) {
+    for(unsigned int i = 0; i < roots.size(); i++) {
         if(roots[i].imag() == 0) // should be more lenient perhaps
             real_roots.push_back(roots[i].real());
     }
index a84a5a7d45f5bc2ce1db4994227d2e6f1bc8b77d..bb041edbe1fd7d77359da7c5857a06bf7b8cf4ab 100644 (file)
@@ -1,5 +1,17 @@
 #include "quadtree.h"
 
+namespace Geom{
+Quad* QuadTree::search(Rect const &r) {
+    return search(r[0].min(), r[1].min(),
+                  r[0].max(), r[1].max());
+}
+
+void QuadTree::insert(Rect const &r, int shape) {
+    insert(r[0].min(), r[1].min(),
+           r[0].max(), r[1].max(), shape);
+}
+
+
 Quad* QuadTree::search(double x0, double y0, double x1, double y1) {
     Quad *q = root;
         
@@ -117,6 +129,8 @@ void QuadTree::erase(Quad *q, int shape) {
     return;
 }
 
+};
+
 /*
   Local Variables:
   mode:c++
index 9b5b75e62344335652758e41f8b50a6656e59b1e..716c56673661cf080bb290e4ff7b5d372ef2b8ca 100644 (file)
@@ -1,6 +1,10 @@
 #include <vector>
 #include <cassert>
 
+#include "d2.h"
+
+namespace Geom{
+
 class Quad{
 public:
     Quad* children[4];
@@ -10,6 +14,23 @@ public:
             children[i] = 0;
     }
     typedef std::vector<int>::iterator iterator;
+    Rect bounds(unsigned i, double x, double y, double d) {
+        double dd = d/2;
+        switch(i % 4) {
+            case 0:
+                return Rect(Interval(x, x+dd), Interval(y, y+dd));
+            case 1:
+                return Rect(Interval(x+dd, x+d), Interval(y, y+dd));
+            case 2:
+                return Rect(Interval(x, x+dd), Interval(y+dd, y+d));
+            case 3:
+                return Rect(Interval(x+dd, x+d), Interval(y+dd, y+d));
+            default: 
+                /* just to suppress warning message
+                 * this case should be never reached */
+                assert(false);
+        }        
+    }
 };
 
 class QuadTree{
@@ -23,9 +44,12 @@ public:
 
     Quad* search(double x0, double y0, double x1, double y1);
     void insert(double x0, double y0, double x1, double y1, int shape);
+    Quad* search(Geom::Rect const &r);
+    void insert(Geom::Rect const &r, int shape);
     void erase(Quad *q, int shape);
 };
 
+};
 
 /*
   Local Variables:
index b30c3a655c99b71f64920bb4be5c6289bce8d6c3..6558303891278c6ef0ea4e71e420c01cf4d3dfab 100644 (file)
@@ -157,6 +157,17 @@ Geom::atan2(D2<SBasis> const &vect, double tol, unsigned order){
     return atan2(Piecewise<D2<SBasis> >(vect),tol,order);
 }
 
+//tan2 is the pseudo-inverse of atan2.  It takes an angle and returns a unit_vector that points in the direction of angle.
+D2<Piecewise<SBasis> >
+Geom::tan2(SBasis const &angle, double tol, unsigned order){
+    return tan2(Piecewise<SBasis>(angle), tol, order);
+}
+
+D2<Piecewise<SBasis> >
+Geom::tan2(Piecewise<SBasis> const &angle, double tol, unsigned order){
+    return D2<Piecewise<SBasis> >(cos(angle, tol, order), sin(angle, tol, order));
+}
+
 //unitVector(x,y) is computed as (b,-a) where a and b are solutions of:
 //     ax+by=0 (eqn1)   and   a^2+b^2=1 (eqn2)
 Piecewise<D2<SBasis> >
index c4f139aa6a742507cd84d8f4aed1a445d8478fda..6e21a6395c51aab5b52212608821feaa72d798fd 100644 (file)
@@ -33,6 +33,14 @@ Piecewise<SBasis>
 atan2(Piecewise<D2<SBasis> >const &vect, 
            double tol=.01, unsigned order=3);
 
+D2<Piecewise<SBasis> >
+tan2(SBasis const &angle, 
+           double tol=.01, unsigned order=3);
+
+D2<Piecewise<SBasis> >
+tan2(Piecewise<SBasis> const &angle, 
+           double tol=.01, unsigned order=3);
+
 Piecewise<D2<SBasis> >
 unitVector(D2<SBasis> const &vect, 
            double tol=.01, unsigned order=3);
index 0a48ff1cf78039b0c4fe22d686a6ea8ca51b1dc5..d6f55c8d8c0a8dde929e6bd43f85578f43b957a7 100644 (file)
@@ -296,6 +296,7 @@ inline std::ostream &operator<< (std::ostream &out_file, const SBasis & p) {
     return out_file;
 }
 
+// These are deprecated, use sbasis-math versions if possible
 SBasis sin(Linear bo, int k);
 SBasis cos(Linear bo, int k);
 
index cd7c36cc33e7a71bcff15e9e34f504ccfda2ba6b..3ec738df1d1be1d65a7d58d4ac33dd0911ab6150 100644 (file)
@@ -86,6 +86,11 @@ Bernsteins::find_bernstein_roots(double const *w, /* The control points  */
         /* if deep enough, return 1 solution at midpoint  */
         if (depth >= MAXDEPTH) {
             //printf("bottom out %d\n", depth);
+            const double Ax = right_t - left_t;
+            const double Ay = w[degree] - w[0];
+            
+            solutions.push_back(left_t - Ax*w[0] / Ay);
+            return;
             solutions.push_back((left_t + right_t) / 2.0);
             return;
         }
@@ -101,6 +106,7 @@ Bernsteins::find_bernstein_roots(double const *w, /* The control points  */
                 solutions.push_back(left_t - Ax*w[0] / Ay);
                 return;
             }
+
     }
 
     /* Otherwise, solve recursively after subdividing control polygon  */