Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 353 insertions, 0 deletions

diff --git a/include/2geom/symbolic/implicit.h b/include/2geom/symbolic/implicit.h
new file mode 100644
index 0000000..82d77cd
--- /dev/null
+++ b/include/2geom/symbolic/implicit.h
@@ -0,0 +1,353 @@
+/*
+ * Routines to compute the implicit equation of a parametric polynomial curve
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 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 _GEOM_SL_IMPLICIT_H_
+#define _GEOM_SL_IMPLICIT_H_
+
+
+
+#include <2geom/symbolic/multipoly.h>
+#include <2geom/symbolic/matrix.h>
+
+
+#include <2geom/exception.h>
+
+#include <array>
+
+
+namespace Geom { namespace SL {
+
+typedef MultiPoly<1, double> MVPoly1;
+typedef MultiPoly<2, double> MVPoly2;
+typedef MultiPoly<3, double> MVPoly3;
+typedef std::array<MVPoly1, 3> poly_vector_type;
+typedef std::array<poly_vector_type, 2> basis_type;
+typedef std::array<double, 3> coeff_vector_type;
+
+namespace detail {
+
+/*
+ *  transform a univariate polynomial f(t) in a 3-variate polynomial
+ *  p(t, x, y) = f(t) * x^i * y^j
+ */
+inline
+void poly1_to_poly3(MVPoly3 & p3, MVPoly1 const& p1, size_t i, size_t j)
+{
+    multi_index_type I = make_multi_index(0, i, j);
+    for (; I[0] < p1.get_poly().size(); ++I[0])
+    {
+        p3.coefficient(I, p1[I[0]]);
+    }
+}
+
+/*
+ *  evaluates the degree of a poly_vector_type, such a degree is defined as:
+ *  deg({p[0](t), p[1](t), p[2](t)}) := {max(deg(p[i](t)), i = 0, 1, 2), k}
+ *  here k is the index where the max is achieved,
+ *  if deg(p[i](t)) == deg(p[j](t)) and i < j then k = i
+ */
+inline
+std::pair<size_t, size_t> deg(poly_vector_type const& p)
+{
+    std::pair<size_t, size_t> d;
+    d.first = p[0].get_poly().real_degree();
+    d.second = 0;
+    size_t k = p[1].get_poly().real_degree();
+    if (d.first < k)
+    {
+        d.first = k;
+        d.second = 1;
+    }
+    k = p[2].get_poly().real_degree();
+    if (d.first < k)
+    {
+        d.first = k;
+        d.second = 2;
+    }
+    return d;
+}
+
+} // end namespace detail
+
+
+/*
+ * A polynomial parametrization could be seen as 1-variety V in R^3,
+ * intersection of two surfaces x = f(t), y = g(t), this variety V has
+ * attached an ideal I in the ring of polynomials in t, x, y with coefficients
+ * on reals; a basis of generators for I is given by p(t, x, y) = x - f(t),
+ * q(t, x, y) = y - g(t); such a basis has the nice property that could be
+ * written as a couple of vectors of dim 3 with entries in R[t]; the original
+ * polinomials p and q can be obtained by doing a dot product between each
+ * vector and the vector {x, y, 1}
+ * As reference you can read the text book:
+ * Ideals, Varieties and Algorithms by Cox, Little, O'Shea
+ */
+inline
+void make_initial_basis(basis_type& b, MVPoly1 const& p, MVPoly1 const& q)
+{
+    // first basis vector
+    b[0][0] = 1;
+    b[0][1] = 0;
+    b[0][2] = -p;
+
+    // second basis vector
+    b[1][0] = 0;
+    b[1][1] = 1;
+    b[1][2] = -q;
+}
+
+/*
+ * Starting from the initial basis for the ideal I is possible to make up
+ * a new basis, still showing off the nice property that each generator is
+ * a moving line that is a linear combination of x, y, 1 where the coefficients
+ * are polynomials in R[t], and moreover each generator is of minimal degree.
+ * Can be proved that given a polynomial parametrization f(t), g(t)
+ * we are able to make up a "micro" basis of generators p(t, x, y), q(t, x, y)
+ * for the ideal I such that the deg(p, t) = m <= n/2 and deg(q, t) = n - m,
+ * where n = max(deg(f(t)), deg(g(t))); this let us halve the order of
+ * the Bezout matrix.
+ * Reference:
+ * Zheng, Sederberg - A Direct Approach to Computing the micro-basis
+ *                    of a Planar Rational Curves
+ * Deng, Chen, Shen - Computing micro-Basis of Rational Curves and Surfaces
+ *                    Using Polynomial Matrix Factorization
+ */
+inline
+void microbasis(basis_type& b, MVPoly1 const& p, MVPoly1 const& q)
+{
+    typedef std::pair<size_t, size_t> degree_pair_t;
+
+    size_t n = std::max(p.get_poly().real_degree(), q.get_poly().real_degree());
+    make_initial_basis(b, p, q);
+    degree_pair_t n0 = detail::deg(b[0]);
+    degree_pair_t n1 = detail::deg(b[1]);
+    size_t d;
+    double r0, r1;
+    //size_t iter = 0;
+    while ((n0.first + n1.first) > n)// && iter < 30)
+    {
+//        ++iter;
+//        std::cout << "iter = " << iter << std::endl;
+//        for (size_t i= 0; i < 2; ++i)
+//            for (size_t j= 0; j < 3; ++j)
+//                std::cout << b[i][j] << std::endl;
+//        std::cout << n0.first << ", " << n0.second << std::endl;
+//        std::cout << n1.first << ", " << n1.second << std::endl;
+//        std::cout << "-----" << std::endl;
+//        if (n0.first < n1.first)
+//        {
+//            d = n1.first - n0.first;
+//            r = b[1][n1.second][n1.first] / b[0][n1.second][n0.first];
+//            for (size_t i = 0; i < b[0].size(); ++i)
+//                b[1][i] -= ((r * b[0][i]).get_poly() << d);
+//            b[1][n1.second][n1.first] = 0;
+//            n1 = detail::deg(b[1]);
+//        }
+//        else
+//        {
+//            d = n0.first - n1.first;
+//            r = b[0][n0.second][n0.first] / b[1][n0.second][n1.first];
+//            for (size_t i = 0; i < b[0].size(); ++i)
+//                b[0][i] -= ((r * b[1][i]).get_poly() << d);
+//            b[0][n0.second][n0.first] = 0;
+//            n0 = detail::deg(b[0]);
+//        }
+
+        // this version shouldn't suffer of ill-conditioning due to
+        // cancellation issue
+        if (n0.first < n1.first)
+        {
+            d = n1.first - n0.first;
+            r0 = b[0][n1.second][n0.first];
+            r1 = b[1][n1.second][n1.first];
+            for (size_t i = 0; i < b[0].size(); ++i)
+            {
+                b[1][i] *= r0;
+                b[1][i] -= ((r1 * b[0][i]).get_poly() << d);
+                // without the following division the modulus grows
+                // beyond the limit of the double type
+                b[1][i] /= r0;
+            }
+            n1 = detail::deg(b[1]);
+        }
+        else
+        {
+            d = n0.first - n1.first;
+            r0 = b[0][n1.second][n0.first];
+            r1 = b[1][n1.second][n1.first];
+
+            for (size_t i = 0; i < b[0].size(); ++i)
+            {
+                b[0][i] *= r1;
+                b[0][i] -= ((r0 * b[1][i]).get_poly() << d);
+                b[0][i] /= r1;
+            }
+            n0 = detail::deg(b[0]);
+        }
+
+    }
+}
+
+/*
+ *  computes the dot product:
+ *  p(t, x, y) = {p0(t), p1(t), p2(t)} . {x, y, 1}
+ */
+inline
+void basis_to_poly(MVPoly3 & p0, poly_vector_type const& v)
+{
+    MVPoly3 p1, p2;
+    detail::poly1_to_poly3(p0, v[0], 1,0);
+    detail::poly1_to_poly3(p1, v[1], 0,1);
+    detail::poly1_to_poly3(p2, v[2], 0,0);
+    p0 += p1;
+    p0 += p2;
+}
+
+
+/*
+ * Make up a Bezout matrix with two basis genarators as input.
+ *
+ * A Bezout matrix is the matrix related to the symmetric bilinear form
+ * (f,g) -> B[f,g] where B[f,g](s,t) = (f(t)*g(s) - f(s)*g(t))/(s-t)
+ * where f, g are polynomials, this function is called a bezoutian.
+ * Given a basis of generators {p(t, x, y), q(t, x, y)} for the ideal I
+ * related to our parametrization x = f(t), y = g(t), we are able to prove
+ * that the implicit equation of such polynomial parametrization can be
+ * evaluated computing the determinant of the Bezout matrix made up using
+ * the polinomial p and q as univariate polynomials in t with coefficients
+ * in R[x,y], so the resulting Bezout matrix will be a matrix with bivariate
+ * polynomials as entries. A Bezout matrix is always symmetric.
+ * Reference:
+ * Sederberg, Zheng - Algebraic Methods for Computer Aided Geometric Design
+ */
+Matrix<MVPoly2>
+make_bezout_matrix (MVPoly3 const& p, MVPoly3 const& q)
+{
+    size_t pdeg = p.get_poly().real_degree();
+    size_t qdeg = q.get_poly().real_degree();
+    size_t n = std::max(pdeg, qdeg);
+
+    Matrix<MVPoly2> BM(n, n);
+    //std::cerr << "rows, columns " << BM.rows() << " , " << BM.columns() << std::endl;
+    for (size_t i = n; i >= 1; --i)
+    {
+        for (size_t j = n; j >= i; --j)
+        {
+            size_t m = std::min(i, n + 1 - j);
+            //std::cerr << "m = " << m << std::endl;
+            for (size_t k = 1; k <= m; ++k)
+            {
+                //BM(i-1,j-1) += (p[j-1+k] * q[i-k] - p[i-k] * q[j-1+k]);
+                BM(n-i,n-j) += (p.coefficient(j-1+k) * q.coefficient(i-k)
+                                - p.coefficient(i-k) * q.coefficient(j-1+k));
+            }
+        }
+    }
+
+    for (size_t i = 0; i < n; ++i)
+    {
+        for (size_t j = 0; j < i; ++j)
+            BM(j,i) = BM(i,j);
+    }
+    return BM;
+}
+
+/*
+ *  Make a matrix that represents a main minor (i.e. with the diagonal
+ *  on the diagonal of the matrix to which it owns) of the Bezout matrix
+ *  with order n-1 where n is the order of the Bezout matrix.
+ *  The minor is obtained by removing the "h"-th row and the "h"-th column,
+ *  and as the Bezout matrix is symmetric.
+ */
+Matrix<MVPoly2>
+make_bezout_main_minor (MVPoly3 const& p, MVPoly3 const& q, size_t h)
+{
+    size_t pdeg = p.get_poly().real_degree();
+    size_t qdeg = q.get_poly().real_degree();
+    size_t n = std::max(pdeg, qdeg);
+
+    Matrix<MVPoly2> BM(n-1, n-1);
+    size_t u = 0, v;
+    for (size_t i = 1; i <= n; ++i)
+    {
+        v = 0;
+        if (i == h)
+        {
+            u = 1;
+            continue;
+        }
+        for (size_t j = 1; j <= i; ++j)
+        {
+            if (j == h)
+            {
+                v = 1;
+                continue;
+            }
+            size_t m = std::min(i, n + 1 - j);
+            for (size_t k = 1; k <= m; ++k)
+            {
+                //BM(i-u-1,j-v-1) += (p[j-1+k] * q[i-k] - p[i-k] * q[j-1+k]);
+                BM(i-u-1,j-v-1) += (p.coefficient(j-1+k) * q.coefficient(i-k)
+                                  - p.coefficient(i-k) * q.coefficient(j-1+k));
+            }
+        }
+    }
+
+    --n;
+    for (size_t i = 0; i < n; ++i)
+    {
+        for (size_t j = 0; j < i; ++j)
+            BM(j,i) = BM(i,j);
+    }
+    return BM;
+}
+
+
+} /*end namespace Geom*/  } /*end namespace SL*/
+
+
+
+
+#endif // _GEOM_SL_IMPLICIT_H_
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :