diff options
Diffstat (limited to 'include/2geom/symbolic/implicit.h')
| -rw-r--r-- | include/2geom/symbolic/implicit.h | 353 |
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 : |