diff options
Diffstat (limited to 'src/toys/curve-intersection-by-implicitization.cpp')
| -rw-r--r-- | src/toys/curve-intersection-by-implicitization.cpp | 300 |
1 files changed, 300 insertions, 0 deletions
diff --git a/src/toys/curve-intersection-by-implicitization.cpp b/src/toys/curve-intersection-by-implicitization.cpp new file mode 100644 index 0000000..e071911 --- /dev/null +++ b/src/toys/curve-intersection-by-implicitization.cpp @@ -0,0 +1,300 @@ +/* + * Show off crossings between two D2<SBasis> curves. + * The intersection points are found by using implicitization tecnique. + * + * 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. + */ + + +#include <toys/path-cairo.h> +#include <toys/toy-framework-2.h> +#include <2geom/d2.h> +#include <2geom/sbasis-poly.h> +#include <2geom/numeric/linear_system.h> +#include <2geom/symbolic/implicit.h> + + +using namespace Geom; + +/* + * helper routines + */ +void poly_to_mvpoly1(SL::MVPoly1 & p, Geom::Poly const& q) +{ + for (size_t i = 0; i < q.size(); ++i) + { + p.coefficient(i, q[i]); + } + p.normalize(); +} + +void mvpoly1_to_poly(Geom::Poly & p, SL::MVPoly1 const& q) +{ + p.resize(q.get_poly().size()); + for (size_t i = 0; i < q.get_poly().size(); ++i) + { + p[i] = q[i]; + } +} + + +/* + * intersection_info + * structure utilized to store intersection info + * + * p - the intersection point + * t0 - the parameter t value at which the first curve pass through p + * t1 - the parameter t value at which the first curve pass through p + */ +struct intersection_info +{ + intersection_info() + {} + + intersection_info(Point const& _p, Coord _t0, Coord _t1) + : p(_p), t0(_t0), t1(_t1) + {} + + Point p; + Coord t0, t1; +}; + +typedef std::vector<intersection_info> intersections_info; + + + +/* + * intersection algorithm + */ +void intersect(intersections_info& xs, D2<SBasis> const& A, D2<SBasis> const& B) +{ + using std::swap; + + // supposing implicitization the most expensive step + // we perform a call to intersect with curve arguments swapped + if (A[0].size() > B[0].size()) + { + intersect(xs, B, A); + for (auto & x : xs) + swap(x.t0, x.t1); + + return; + } + + // convert A from symmetric power basis to power basis + Geom::Poly A0 = sbasis_to_poly(A[0]); + Geom::Poly A1 = sbasis_to_poly(A[1]); + + // convert to MultiPoly type + SL::MVPoly1 Af, Ag; + poly_to_mvpoly1(Af, A0); + poly_to_mvpoly1(Ag, A1); + + // compute a basis of the ideal related to the curve A + // in vector form + Geom::SL::basis_type b; + // if we compute the micro-basis the bezout matrix is made up + // by one only entry so we can't do the inversion step. + if (A0.size() == 3) + { + make_initial_basis(b, Af, Ag); + } + else + { + microbasis(b, Af, Ag); + } + + // we put the basis in of the form of two independent moving line + Geom::SL::MVPoly3 p, q; + basis_to_poly(p, b[0]); + basis_to_poly(q, b[1]); + + // compute the Bezout matrix and the implicit equation of the curve A + Geom::SL::Matrix<Geom::SL::MVPoly2> BZ = make_bezout_matrix(p, q); + SL::MVPoly2 ic = determinant_minor(BZ); + ic.normalize(); + + + // convert B from symmetric power basis to power basis + Geom::Poly B0 = sbasis_to_poly(B[0]); + Geom::Poly B1 = sbasis_to_poly(B[1]); + + // convert to MultiPoly type + SL::MVPoly1 Bf, Bg; + poly_to_mvpoly1(Bf, B0); + poly_to_mvpoly1(Bg, B1); + + // evaluate the implicit equation of A on B + // so we get an s(t) polynomial that give us + // the t values for B at which intersection happens + SL::MVPoly1 s = ic(Bf, Bg); + + // convert s(t) to Poly type, in order to use the real_solve function + Geom::Poly z; + mvpoly1_to_poly(z, s); + + // compute t values for the curve B at which intersection happens + std::vector<double> sol = solve_reals(z); + + // filter the found solutions wrt the domain interval [0,1] of B + // and compute the related point coordinates + std::vector<double> pt; + pt.reserve(sol.size()); + std::vector<Point> points; + points.reserve(sol.size()); + for (double & i : sol) + { + if (i >= 0 && i <= 1) + { + pt.push_back(i); + points.push_back(B(pt.back())); + } + } + + // case: A is parametrized by polynomial of degree 1 + // we compute the t values of A at the intersection points + // and filter the results wrt the domain interval [0,1] + double t; + xs.clear(); + xs.reserve(pt.size()); + if (A0.size() == 2) + { + for (size_t i = 0; i < points.size(); ++i) + { + t = (points[i][X] - A0[0]) / A0[1]; + if (t >= 0 && t <= 1) + { + xs.push_back(intersection_info(points[i], t, pt[i])); + } + } + return; + } + + // general case + // we compute the value of the parameter t of A at each intersection point + // and we filter the final result wrt the domain interval [0,1] + // the computation is performed by using the inversion formula for each point + // As reference see: + // Sederberger - Computer Aided Geometric Design + // par 16.5 - Implicitization and Inversion + size_t n = BZ.rows(); + Geom::NL::Matrix BZN(n, n); + Geom::NL::MatrixView BZV(BZN, 0, 0, n-1, n-1); + Geom::NL::VectorView cv = BZN.column_view(n-1); + Geom::NL::VectorView bv(cv, n-1); + Geom::NL::LinearSystem ls(BZV, bv); + for (size_t i = 0; i < points.size(); ++i) + { + // evaluate the first main minor of order n-1 at each intersection point + polynomial_matrix_evaluate(BZN, BZ, points[i]); + // solve the linear system with the powers of t as unknowns + ls.SV_solve(); + // the last element contains the t value + t = -ls.solution()[n-2]; + // filter with respect to the domain of A + if (t >= 0 && t <= 1) + { + xs.push_back(intersection_info(points[i], t, pt[i])); + } + } +} + + + +class IntersectImplicit : public Toy +{ + + void draw( cairo_t *cr, std::ostringstream *notify, + int width, int height, bool save, std::ostringstream *timer_stream) override + { + cairo_set_line_width (cr, 0.3); + cairo_set_source_rgba (cr, 0.8, 0., 0, 1); + D2<SBasis> A = pshA.asBezier(); + cairo_d2_sb(cr, A); + cairo_stroke(cr); + cairo_set_source_rgba (cr, 0.0, 0., 0, 1); + D2<SBasis> B = pshB.asBezier(); + cairo_d2_sb(cr, B); + cairo_stroke(cr); + + intersect(xs, A, B); + for (auto & x : xs) + { + draw_handle(cr, x.p); + } + + Toy::draw(cr, notify, width, height, save,timer_stream); + } + + +public: + IntersectImplicit(unsigned int _A_bez_ord, unsigned int _B_bez_ord) + : A_bez_ord(_A_bez_ord), B_bez_ord(_B_bez_ord) + { + handles.push_back(&pshA); + for (unsigned int i = 0; i <= A_bez_ord; ++i) + pshA.push_back(Geom::Point(uniform()*400, uniform()*400)); + handles.push_back(&pshB); + for (unsigned int i = 0; i <= B_bez_ord; ++i) + pshB.push_back(Geom::Point(uniform()*400, uniform()*400)); + + } + +private: + unsigned int A_bez_ord, B_bez_ord; + PointSetHandle pshA, pshB; + intersections_info xs; +}; + + +int main(int argc, char **argv) +{ + unsigned int A_bez_ord = 4; + unsigned int B_bez_ord = 6; + if(argc > 1) + sscanf(argv[1], "%d", &A_bez_ord); + if(argc > 2) + sscanf(argv[2], "%d", &B_bez_ord); + + + init( argc, argv, new IntersectImplicit(A_bez_ord, B_bez_ord)); + return 0; +} + + +/* + 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 : |