Merging in from trunk

(bzr r8254.1.37)

7 files changed, 145 insertions, 32 deletions

diff --git a/src/2geom/bezier.h b/src/2geom/bezier.h
index 4ab965f42..9e68d93ae 100644
--- a/src/2geom/bezier.h
+++ b/src/2geom/bezier.h
@@ -46,35 +46,32 @@ namespace Geom {
 inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, unsigned order) {
 /*
  *  Bernstein : 
- *	Evaluate a Bernstein function at a particular parameter value
+ *      Evaluate a Bernstein function at a particular parameter value
  *      Fill in control points for resulting sub-curves.
  * 
  */
 
     unsigned N = order+1;
-    std::valarray<Coord> vtemp(2*N);
+    std::valarray<Coord> row(N);
     for (unsigned i = 0; i < N; i++)
-        vtemp[i] = v[i];
+        row[i] = v[i];
 
     // Triangle computation
     const double omt = (1-t);
     if(left)
-        left[0] = vtemp[0];
+        left[0] = row[0];
     if(right)
-        right[order] = vtemp[order];
-    double *prev_row = &vtemp[0];
-    double *row = &vtemp[N];
+        right[order] = row[order];
     for (unsigned i = 1; i < N; i++) {
         for (unsigned j = 0; j < N - i; j++) {
-            row[j] = omt*prev_row[j] + t*prev_row[j+1];
+            row[j] = omt*row[j] + t*row[j+1];
         }
         if(left)
             left[i] = row[0];
         if(right)
             right[order-i] = row[order-i];
-        std::swap(prev_row, row);
     }
-    return (prev_row[0]);
+    return (row[0]);
 /*
     Coord vtemp[order+1][order+1];
 
@@ -97,6 +94,20 @@ inline Coord subdivideArr(Coord t, Coord const *v, Coord *left, Coord *right, un
             return (vtemp[order][0]);*/
 }
 
+template <typename T>
+inline T bernsteinValueAt(double t, T const *c_, unsigned n) {
+    double u = 1.0 - t;
+    double bc = 1;
+    double tn = 1;
+    T tmp = c_[0]*u;
+    for(unsigned i=1; i<n; i++){
+        tn = tn*t;
+        bc = bc*(n-i+1)/i;
+        tmp = (tmp + tn*bc*c_[i])*u;
+    }
+    return (tmp + tn*t*c_[n]);
+}
+
 
 class Bezier {
 private:
@@ -225,24 +236,34 @@ public:
     //inline Coord const &operator[](unsigned ix) const { return c_[ix]; }
     inline void setPoint(unsigned ix, double val) { c_[ix] = val; }
 
-    /* This is inelegant, as it uses several extra stores.  I think there might be a way to
-     * evaluate roughly in situ. */
-
+    /**
+    *  The size of the returned vector equals n_derivs+1.
+    */
     std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const {
-        std::vector<Coord> val_n_der;
+        /* This is inelegant, as it uses several extra stores.  I think there might be a way to
+         * evaluate roughly in situ. */
+
+         // initialize return vector with zeroes, such that we only need to replace the non-zero derivs
+        std::vector<Coord> val_n_der(n_derivs + 1, Coord(0.0));
+
+        // initialize temp storage variables
         std::valarray<Coord> d_(order()+1);
-        unsigned nn = n_derivs + 1; 	// the size of the result vector equals n_derivs+1 ...
-        if(nn > order())
-            nn = order()+1;		// .. but with a maximum of order() + 1!
-        for(unsigned i = 0; i < size(); i++)
+        for (unsigned i = 0; i < size(); i++) {
             d_[i] = c_[i];
-        for(unsigned di = 0; di < nn; di++) {
-            val_n_der.push_back(subdivideArr(t, &d_[0], NULL, NULL, order() - di));
-            for(unsigned i = 0; i < order() - di; i++) {
+        }
+
+        unsigned nn = n_derivs + 1;
+        if(n_derivs > order()) {
+            nn = order()+1; // only calculate the non zero derivs
+        }
+        for (unsigned di = 0; di < nn; di++) {
+            //val_n_der[di] = (subdivideArr(t, &d_[0], NULL, NULL, order() - di));
+            val_n_der[di] = bernsteinValueAt(t, &d_[0], order() - di);
+            for (unsigned i = 0; i < order() - di; i++) {
                 d_[i] = (order()-di)*(d_[i+1] - d_[i]);
             }
         }
-        val_n_der.resize(n_derivs);
+
         return val_n_der;
     }
 
@@ -257,6 +278,11 @@ public:
         find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, 0.0, 1.0);
         return solutions;
     }
+    std::vector<double> roots(Interval const ivl) const {
+        std::vector<double> solutions;
+        find_bernstein_roots(&const_cast<std::valarray<Coord>&>(c_)[0], order(), solutions, 0, ivl[0], ivl[1]);
+        return solutions;
+    }
 };
 
 
diff --git a/src/2geom/d2-sbasis.cpp b/src/2geom/d2-sbasis.cpp
index 2c52e4782..aef989fc7 100644
--- a/src/2geom/d2-sbasis.cpp
+++ b/src/2geom/d2-sbasis.cpp
@@ -58,8 +58,19 @@ Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){
   return result;
 }
 
-Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, 
-		      Piecewise<D2<SBasis> > const &b){
+/** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
+ * @return  \f[
+ *      f(t) \rightarrow \left\{ 
+ *      \begin{array}{c}
+ *      a_1 \bullet b_1 \\
+ *      a_2 \bullet b_2 \\
+ *      \ldots \\
+ *      a_n \bullet b_n \\
+ *      \end{array}\right.
+ * \f]
+ * @relates Piecewise */
+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b)
+{
   Piecewise<SBasis > result;
   if (a.empty() || b.empty()) return result;
   Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
@@ -72,6 +83,30 @@ Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a,
   return result;
 }
 
+/** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
+ * @return  \f[
+ *      f(t) \rightarrow \left\{ 
+ *      \begin{array}{c}
+ *      a_1 \bullet b \\
+ *      a_2 \bullet b \\
+ *      \ldots \\
+ *      a_n \bullet b \\
+ *      \end{array}\right.
+ * \f]
+ * @relates Piecewise */
+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Point const &b)
+{
+  Piecewise<SBasis > result;
+  if (a.empty()) return result;
+
+  result.push_cut(a.cuts.front());
+  for (unsigned i = 0; i < a.size(); ++i){
+    result.push(dot(a.segs[i],b), a.cuts[i+1]);
+  }
+  return result;
+}
+
+
 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, 
 			Piecewise<D2<SBasis> > const &b){
   Piecewise<SBasis > result;
diff --git a/src/2geom/d2-sbasis.h b/src/2geom/d2-sbasis.h
index c61052da7..d404e0618 100644
--- a/src/2geom/d2-sbasis.h
+++ b/src/2geom/d2-sbasis.h
@@ -73,6 +73,7 @@ Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a);
 D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a);
 Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &a);
 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
+Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Point const &b);
 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b);
 
 Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Matrix const &m);
diff --git a/src/2geom/d2.h b/src/2geom/d2.h
index afa00b40d..b2a0f8866 100644
--- a/src/2geom/d2.h
+++ b/src/2geom/d2.h
@@ -97,10 +97,10 @@ class D2{
     }
     std::vector<Point > valueAndDerivatives(double t, unsigned n) const {
         std::vector<Coord> x = f[X].valueAndDerivatives(t, n),
-                           y = f[Y].valueAndDerivatives(t, n);
-        std::vector<Point> res;
+                           y = f[Y].valueAndDerivatives(t, n); // always returns a vector of size n+1
+        std::vector<Point> res(n+1);
         for(unsigned i = 0; i <= n; i++) {
-            res.push_back(Point(x[i], y[i]));
+            res[i] = Point(x[i], y[i]);
         }
         return res;
     }
@@ -321,6 +321,25 @@ dot(D2<T> const & a, D2<T> const & b) {
     return r;
 }
 
+/** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
+ * @return \f$a \bullet b = a_X b_X + a_Y b_Y\f$.
+ * @relates D2 */
+template <typename T>
+inline T
+dot(D2<T> const & a, Point const & b) {
+    boost::function_requires<AddableConcept<T> >();
+    boost::function_requires<ScalableConcept<T> >();
+
+    T r;
+    for(unsigned i = 0; i < 2; i++) {
+        r += a[i] * b[i];
+    }
+    return r;
+}
+
+/** @brief Calculates the 'cross product' or 'outer product' of \c a and \c b
+ * @return \f$a \times b = a_Y b_X - a_X b_Y\f$.
+ * @relates D2 */
 template <typename T>
 inline T
 cross(D2<T> const & a, D2<T> const & b) {
diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp
index 04a21d624..cc91743b1 100644
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
@@ -179,7 +179,8 @@ Matrix Matrix::inverse() const {
     Matrix d;
 
     Geom::Coord const determ = det();
-    if (!are_near(determ, 0.0)) {
+    // the numerical precision of the determinant must be significant
+    if (fabs(determ) > 1e-18) {
         Geom::Coord const ideterm = 1.0 / determ;
 
         d._c[0] =  _c[3] * ideterm;
diff --git a/src/2geom/piecewise.h b/src/2geom/piecewise.h
index a5be42587..a0628daf1 100644
--- a/src/2geom/piecewise.h
+++ b/src/2geom/piecewise.h
@@ -58,7 +58,7 @@ namespace Geom {
  *      \begin{array}{cc}
  *      s_1,& t <= c_2 \\
  *      s_2,& c_2 <= t <= c_3\\
- *      \ldots
+ *      \ldots \\
  *      s_n,& c_n <= t
  *      \end{array}\right.
  * \f]
@@ -105,6 +105,10 @@ class Piecewise {
     inline output_type lastValue() const {
         return valueAt(cuts.back());
     }
+
+    /**
+    *  The size of the returned vector equals n_derivs+1.
+    */
     std::vector<output_type> valueAndDerivatives(double t, unsigned n_derivs) const {
         unsigned n = segN(t);
         std::vector<output_type> ret, val = segs[n].valueAndDerivatives(segT(t, n), n_derivs);
@@ -115,6 +119,7 @@ class Piecewise {
         }
         return ret;
     }
+
     //TODO: maybe it is not a good idea to have this?
     Piecewise<T> operator()(SBasis f);
     Piecewise<T> operator()(Piecewise<SBasis>f);
@@ -773,10 +778,10 @@ Piecewise<T> reverse(Piecewise<T> const &f) {
     return ret;
 }
 
-
 /**
  *  Interpolates between a and b.
  *  \return a if t = 0, b if t = 1, or an interpolation between a and b for t in [0,1]
+ *  \relates Piecewise
  */
 template<typename T>
 Piecewise<T> lerp(double t, Piecewise<T> const &a, Piecewise<T> b) {
diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h
index f2902750c..f1fd67867 100644
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
@@ -64,33 +64,58 @@ public:
     void moveTo(Point p) {
         finish();
         _path.start(p);
+        _start_p = p;
         _in_path = true;
     }
 //TODO: what if _in_path = false?
 
     void hlineTo(Coord v) {
-    	_path.template appendNew<HLineSegment>(Point(v, _path.finalPoint()[Y]));
+    // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
+        _path.template appendNew<HLineSegment>(Point(v, _path.finalPoint()[Y]));
     }
 
     void vlineTo(Coord v) {
-    	_path.template appendNew<VLineSegment>(Point(_path.finalPoint()[X], v));
+    // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
+        _path.template appendNew<VLineSegment>(Point(_path.finalPoint()[X], v));
     }
 
     void lineTo(Point p) {
+        // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
         _path.template appendNew<LineSegment>(p);
     }
 
     void quadTo(Point c, Point p) {
+        // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
         _path.template appendNew<QuadraticBezier>(c, p);
     }
 
     void curveTo(Point c0, Point c1, Point p) {
+        // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
         _path.template appendNew<CubicBezier>(c0, c1, p);
     }
 
     void arcTo(double rx, double ry, double angle,
                bool large_arc, bool sweep, Point p)
     {
+        // check for implicit moveto, like in: "M 1,1 L 2,2 z l 2,2 z"
+        if (!_in_path) {
+            moveTo(_start_p);
+        }
         _path.template appendNew<SVGEllipticalArc>(rx, ry, angle,
                                                  large_arc, sweep, p);
     }
@@ -113,6 +138,7 @@ protected:
     bool _in_path;
     OutputIterator _out;
     Path _path;
+    Point _start_p;
 };
 
 typedef std::back_insert_iterator<std::vector<Path> > iter;