From 61f3ab8f23f4c924d455757bf3e65f8487521b5a Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 13 Apr 2024 13:57:42 +0200 Subject: Adding upstream version 1.3. Signed-off-by: Daniel Baumann --- src/toys/smash-intersector.cpp | 583 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 583 insertions(+) create mode 100644 src/toys/smash-intersector.cpp (limited to 'src/toys/smash-intersector.cpp') diff --git a/src/toys/smash-intersector.cpp b/src/toys/smash-intersector.cpp new file mode 100644 index 0000000..b01acfb --- /dev/null +++ b/src/toys/smash-intersector.cpp @@ -0,0 +1,583 @@ +/* + * Diffeomorphism-based intersector: given two curves + * M(t)=(x(t),y(t)) and N(u)=(X(u),Y(u)) + * and supposing M is a graph over the x-axis, we compute y(x) and solve + * Y(u) - y(X(u)) = 0 + * to get the intersections of the two curves... + * + * Notice the result can be far from intuitive because of the choice we have + * to make to consider a curve as a graph over x or y. For instance the two + * branches of xy=eps are never close from this point of view (!)... + * + * Authors: + * J.-F. Barraud + * Copyright 2010 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 <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/path.h> +#include <2geom/bezier-to-sbasis.h> +#include <2geom/sbasis-geometric.h> +#include +#include + +#include +#include +#include +#include +#include + +#include <2geom/orphan-code/intersection-by-smashing.h> +#include "../2geom/orphan-code/intersection-by-smashing.cpp" + +using namespace Geom; + +#define VERBOSE 0 + +static double exp_rescale(double x){ return ::pow(10, x);} +std::string exp_formatter(double x){ return default_formatter(exp_rescale(x));} + + + +#if 0 +//useless here; +Piecewise > linearizeCusps( D2 f, double tol){ + D2 df = derivative( f ); + std::vector xdoms = level_set( df[X], 0., tol); + std::vector ydoms = level_set( df[Y], 0., tol); + std::vector doms; + //TODO: use order!! + for ( unsigned i=0; i > result; + if (doms.size() == 0 ) return Piecewise >(f); + if (doms[0].min() > 0 ){ + result.cuts.push_back( 0 ); + result.cuts.push_back( doms[0].min() ); + result.segs.push_back( portion( f, Interval( 0, doms[0].min() ) ) ); + } + for ( unsigned i=0; i( Linear( a[X], b[X] ), Linear( a[Y], b[Y] ) ) ); + result.cuts.push_back( doms[i].max() ); + result.segs.push_back( D2( Linear( b[X], c[X] ), Linear( b[Y], c[Y] ) ) ); + double t = ( i+1 == doms.size() )? 1 : doms[i+1].min(); + result.cuts.push_back( t ); + result.segs.push_back( portion( f, Interval( doms[i].max(), t ) ) ); + } + return result; +} +#endif + +#if 0 +/* Computes the intersection of two sets given as (ordered) union intervals. + */ +std::vector intersect( std::vector const &a, std::vector const &b){ + std::vector result; + //TODO: use order! + for (unsigned i=0; i < a.size(); i++){ + for (unsigned j=0; j < b.size(); j++){ + OptInterval c( a[i] ); + c &= b[j]; + if ( c ) { + result.push_back( *c ); + } + } + } + return result; +} + +/* Computes the top and bottom boundaries of the L_\infty neighborhood + * of a curve. The curve is supposed to be a graph over the x-axis. + */ +void computeLinfinityNeighborhood( D2 const &f, double tol, D2 > &topside, D2 > &botside ){ + double signx = ( f[X].at0() > f[X].at1() )? -1 : 1; + double signy = ( f[Y].at0() > f[Y].at1() )? -1 : 1; + + Piecewise > top, bot; + top = Piecewise > (f); + top.cuts.insert( top.cuts.end(), 2); + top.segs.insert( top.segs.end(), D2(Linear( f[X].at1(), f[X].at1()+2*tol*signx), + Linear( f[Y].at1() )) ); + bot = Piecewise >(f); + bot.cuts.insert( bot.cuts.begin(), - 1 ); + bot.segs.insert( bot.segs.begin(), D2(Linear( f[X].at0()-2*tol*signx, f[X].at0()), + Linear( f[Y].at0() )) ); + top += Point(-tol*signx, tol); + bot += Point( tol*signx, -tol); + + if ( signy < 0 ){ + swap( top, bot ); + top += Point( 0, 2*tol); + bot += Point( 0, -2*tol); + } + topside = make_cuts_independent(top); + botside = make_cuts_independent(bot); +} + + +/*Compute top and bottom boundaries of the L^infty nbhd of the graph of a *monotonic* function f. + * if f is increasing, it is given by [f(t-tol)-tol, f(t+tol)+tol]. + * if not, it is [f(t+tol)-tol, f(t-tol)+tol]. + */ +void computeLinfinityNeighborhood( Piecewise const &f, double tol, Piecewise &top, Piecewise &bot){ + top = f + tol; + top.offsetDomain( - tol ); + top.cuts.insert( top.cuts.end(), f.domain().max() + tol); + top.segs.insert( top.segs.end(), SBasis(Linear( f.lastValue() + tol )) ); + + bot = f - tol; + bot.offsetDomain( tol ); + bot.cuts.insert( bot.cuts.begin(), f.domain().min() - tol); + bot.segs.insert( bot.segs.begin(), SBasis(Linear( f.firstValue() - tol )) ); + + if ( f.firstValue() > f.lastValue() ){ + swap( top, bot ); + top += 2*tol; + bot -= 2*tol; + } +} + +std::vector level_set( D2 const &f, Rect region){ + std::vector x_in_reg = level_set( f[X], region[X] ); + std::vector y_in_reg = level_set( f[Y], region[Y] ); + std::vector result = intersect ( x_in_reg, y_in_reg ); + return result; +} + +void prolongateByConstants( Piecewise &f, double paddle_width ){ + if ( f.size() == 0 ) return; //do we have a covention about the domain of empty pwsb? + f.cuts.insert( f.cuts.begin(), f.cuts.front() - paddle_width ); + f.segs.insert( f.segs.begin(), SBasis( f.segs.front().at0() ) ); + f.cuts.insert( f.cuts.end(), f.cuts.back() + paddle_width ); + f.segs.insert( f.segs.end(), SBasis( f.segs.back().at1() ) ); +} + + + +/* Returns the intervals over which the curve keeps its slope + * in one of the 8 sectors delimited by x=0, y=0, y=x, y=-x. + * WARNING: both curves are supposed to be a graphs over x or y axis, + * and the smaller the slopes the better (typically <=45°). + */ +std::vector > smash_intersect( D2 const &a, D2 const &b, + double tol, cairo_t *cr , bool draw_more_stuff=false ){ + + std::vector > res; + + // a and b or X and Y may have to be exchanged, so make local copies. + D2 aa = a; + D2 bb = b; + bool swapresult = false; + bool swapcoord = false;//debug only! + + if ( draw_more_stuff ){ + cairo_set_line_width (cr, 3); + cairo_set_source_rgba(cr, .5, .9, .7, 1 ); + cairo_d2_sb(cr, aa); + cairo_d2_sb(cr, bb); + cairo_stroke(cr); + } + +#if 1 + //if the (enlarged) bounding boxes don't intersect, stop. + if ( !draw_more_stuff ){ + OptRect abounds = bounds_fast( a ); + OptRect bbounds = bounds_fast( b ); + if ( !abounds || !bbounds ) return res; + abounds->expandBy(tol); + if ( !(abounds->intersects(*bbounds))){ + return res; + } + } +#endif + + //Choose the best curve to be re-parametrized by x or y values. + OptRect dabounds = bounds_exact(derivative(a)); + OptRect dbbounds = bounds_exact(derivative(b)); + if ( dbbounds->min().length() > dabounds->min().length() ){ + aa=b; + bb=a; + swap( dabounds, dbbounds ); + swapresult = true; + } + + //Choose the best coordinate to use as new parameter + double dxmin = std::min( abs((*dabounds)[X].max()), abs((*dabounds)[X].min()) ); + double dymin = std::min( abs((*dabounds)[Y].max()), abs((*dabounds)[Y].min()) ); + if ( (*dabounds)[X].max()*(*dabounds)[X].min() < 0 ) dxmin=0; + if ( (*dabounds)[Y].max()*(*dabounds)[Y].min() < 0 ) dymin=0; + assert (dxmin>=0 && dymin>=0); + + if (dxmin < dymin) { + aa = D2( aa[Y], aa[X] ); + bb = D2( bb[Y], bb[X] ); + swapcoord = true; + } + + //re-parametrize aa by the value of x. + Interval x_range_strict( aa[X].at0(), aa[X].at1() ); + Piecewise y_of_x = pw_compose_inverse(aa[Y],aa[X], 2, 1e-5); + + //Compute top and bottom boundaries of the L^infty nbhd of aa. + Piecewise top_ay, bot_ay; + computeLinfinityNeighborhood( y_of_x, tol, top_ay, bot_ay); + + Interval ax_range = top_ay.domain();//i.e. aa[X] domain ewpanded by tol. + + if ( draw_more_stuff ){ + Piecewise dbg_x( SBasis( Linear( top_ay.domain().min(), top_ay.domain().max() ) ) ); + dbg_x.setDomain( top_ay.domain() ); + D2 > dbg_side ( Piecewise( SBasis( Linear( 0 ) ) ), + Piecewise( SBasis( Linear( 0, 2*tol) ) ) ); + + D2 > dbg_rgn; + unsigned h = ( swapcoord ) ? Y : X; + dbg_rgn[h].concat ( dbg_x ); + dbg_rgn[h].concat ( dbg_side[X] + dbg_x.lastValue() ); + dbg_rgn[h].concat ( reverse(dbg_x) ); + dbg_rgn[h].concat ( dbg_side[X] + dbg_x.firstValue() ); + + dbg_rgn[1-h].concat ( bot_ay ); + dbg_rgn[1-h].concat ( dbg_side[Y] + bot_ay.lastValue() ); + dbg_rgn[1-h].concat ( reverse(top_ay) ); + dbg_rgn[1-h].concat ( reverse( dbg_side[Y] ) + bot_ay.firstValue() ); + + cairo_set_line_width (cr, 1.); + cairo_set_source_rgba(cr, 0., 1., 0., .75 ); + cairo_d2_pw_sb(cr, dbg_rgn ); + cairo_stroke(cr); + + D2 bbb = bb; + if ( swapcoord ) swap( bbb[X], bbb[Y] ); + //Piecewise > dbg_rgnB = neighborhood( bbb, tol ); + D2 > dbg_topB, dbg_botB; + computeLinfinityNeighborhood( bbb, tol, dbg_topB, dbg_botB ); + cairo_set_line_width (cr, 1.); + cairo_set_source_rgba(cr, .2, 8., .2, .4 ); +// cairo_pw_d2_sb(cr, dbg_rgnB ); + cairo_d2_pw_sb(cr, dbg_topB ); + cairo_d2_pw_sb(cr, dbg_botB ); + cairo_stroke(cr); + } + + std::vector bx_in_ax_range = level_set(bb[X], ax_range ); + + // find times when bb is in the neighborhood of aa. + std::vector tbs; + for (unsigned i=0; i > bb_in; + bb_in[X] = Piecewise ( portion( bb[X], bx_in_ax_range[i] ) ); + bb_in[Y] = Piecewise ( portion( bb[Y], bx_in_ax_range[i]) ); + bb_in[X].setDomain( bx_in_ax_range[i] ); + bb_in[Y].setDomain( bx_in_ax_range[i] ); + + Piecewise h; + Interval level; + h = bb_in[Y] - compose( top_ay, bb_in[X] ); + level = Interval( -infinity(), 0 ); + std::vector rts_lo = level_set( h, level); + h = bb_in[Y] - compose( bot_ay, bb_in[X] ); + level = Interval( 0, infinity()); + std::vector rts_hi = level_set( h, level); + + std::vector rts = intersect( rts_lo, rts_hi ); + tbs.insert(tbs.end(), rts.begin(), rts.end() ); + } + + std::vector > result(tbs.size(),std::pair()); + + /* for each solution I, find times when aa is in the neighborhood of bb(I). + * (Note: the preimage of bb[X](I) by aa[X], enlarged by tol, is a good approximation of this: + * it would give points in the 2*tol neighborhood of bb (if the slope of aa is never more than 1). + * + faster computation. + * - implies little jumps depending on the subdivision of the input curve into monotonic pieces + * and on the choice of preferred axis. If noticeable, these jumps would feel random to the user :-( + */ + for (unsigned j=0; j tas; + Piecewise fat_y_of_x = y_of_x; + prolongateByConstants( fat_y_of_x, 100*(1+tol) ); + + D2 > top_b, bot_b; + D2 bbj = portion( bb, tbs[j] ); + computeLinfinityNeighborhood( bbj, tol, top_b, bot_b ); + + Piecewise h; + Interval level; + h = top_b[Y] - compose( fat_y_of_x, top_b[X] ); + level = Interval( +infinity(), 0 ); + std::vector rts_top = level_set( h, level); + for (unsigned idx=0; idx < rts_top.size(); idx++){ + rts_top[idx] = Interval( top_b[X].valueAt( rts_top[idx].min() ), + top_b[X].valueAt( rts_top[idx].max() ) ); + } + assert( rts_top.size() == 1 ); + + h = bot_b[Y] - compose( fat_y_of_x, bot_b[X] ); + level = Interval( 0, -infinity()); + std::vector rts_bot = level_set( h, level); + for (unsigned idx=0; idx < rts_bot.size(); idx++){ + rts_bot[idx] = Interval( bot_b[X].valueAt( rts_bot[idx].min() ), + bot_b[X].valueAt( rts_bot[idx].max() ) ); + } + assert( rts_bot.size() == 1 ); + +#if VERBOSE + printf("range(aa[X]) = [%f, %f];\n", y_of_x.domain().min(), y_of_x.domain().max()); + printf("range(bbj[X]) = [%f, %f]; tol= %f\n", bbj[X].at0(), bbj[X].at1(), tol); + + printf("rts_top = "); + for (unsigned dbgi=0; dbgi= x_range_strict.max() ){ + tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 1 : 0 ) ); + }else{ + tas = level_set(aa[X], x_dom ); + } + +#if VERBOSE + if ( tas.size() != 1 ){ + printf("Error: preimage of [%f, %f] by x:[0,1]->[%f, %f] is ", + x_dom.min(), x_dom.max(), x_range_strict.min(), x_range_strict.max()); + if ( tas.size() == 0 ){ + printf( "empty.\n"); + }else{ + printf("\n [%f,%f]", tas[0].min(), tas[0].max() ); + for (unsigned toto=1; toto A = handles_to_sbasis(psh.pts.begin(), A_bez_ord-1); + D2 B = handles_to_sbasis(psh.pts.begin()+A_bez_ord, B_bez_ord-1); + cairo_set_line_width (cr, .8); + cairo_set_source_rgba(cr,0.,0.,0.,.6); + cairo_d2_sb(cr, A); + cairo_d2_sb(cr, B); + cairo_stroke(cr); + + Rect tolbytol( anchor.pos, anchor.pos ); + tolbytol.expandBy( tol ); + cairo_rectangle(cr, tolbytol); + cairo_stroke(cr); +/* + Piecewise > smthA = linearizeCusps(A+Point(0,10), tol); + cairo_set_line_width (cr, 1.); + cairo_set_source_rgba(cr, 1., 0., 1., 1. ); + cairo_pw_d2_sb(cr, smthA); + cairo_stroke(cr); +*/ + + std::vector Acuts = monotonicSplit(A); + std::vector Bcuts = monotonicSplit(B); + +#if 0 + for (unsigned i=0; i Ai = portion( A, Acuts[i]); + cairo_set_line_width (cr, .2); + cairo_set_source_rgba(cr, 0., 0., 0., 1. ); + draw_cross(cr, Ai.at0()); + cairo_stroke(cr); + for (unsigned j=0; j > my_intersections; + D2 Bj = portion( B, Bcuts[j]); + cairo_set_line_width (cr, .2); + cairo_set_source_rgba(cr, 0., 0., 0., 1. ); + draw_cross(cr, Bj.at0()); + cairo_stroke(cr); + } + } +#endif + + std::vector my_intersections; + my_intersections = smash_intersect( A, B, tol ); + + for (auto & my_intersection : my_intersections){ + cairo_set_line_width (cr, 2.5); + cairo_set_source_rgba(cr, 1., 0., 0., .8 ); + cairo_d2_sb(cr, portion( A, my_intersection.times[X])); + cairo_stroke(cr); + cairo_set_line_width (cr, 2.5); + cairo_set_source_rgba(cr, 0., 0., 1., .8 ); + cairo_d2_sb(cr, portion( B, my_intersection.times[Y])); + cairo_stroke(cr); + } +#if 0 + + unsigned apiece( slidera.value()/100. * Acuts.size() ); + unsigned bpiece( sliderb.value()/100. * Bcuts.size() ); + + + for (unsigned i=0; i Ai = portion( A, Acuts[i]); + for (unsigned j=0; j my_intersections; + D2 Bj = portion( B, Bcuts[j]); + bool draw_more = toggle.on && i == apiece && j == bpiece; +// my_intersections = smash_intersect( Ai, Bj, tol, cr, draw_more ); + my_intersections = monotonic_smash_intersect( Ai, Bj, tol ); + + for (unsigned k=0; k 2) + sscanf(argv[2], "%d", &B_bez_ord); + if(argc > 1) + sscanf(argv[1], "%d", &A_bez_ord); + + init( argc, argv, new Intersector(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 : -- cgit v1.2.3