diff options
Diffstat (limited to '')
| -rw-r--r-- | src/2geom/orphan-code/intersection-by-smashing.cpp | 349 |
1 files changed, 349 insertions, 0 deletions
diff --git a/src/2geom/orphan-code/intersection-by-smashing.cpp b/src/2geom/orphan-code/intersection-by-smashing.cpp new file mode 100644 index 0000000..02e44b1 --- /dev/null +++ b/src/2geom/orphan-code/intersection-by-smashing.cpp @@ -0,0 +1,349 @@ +#include <2geom/d2.h> +#include <2geom/sbasis.h> +#include <2geom/sbasis-geometric.h> +#include <2geom/orphan-code/intersection-by-smashing.h> + +#include <cstdlib> +#include <cstdio> +#include <vector> +#include <algorithm> + +namespace Geom { +using namespace Geom; + +/* + * 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. + */ +static +void computeLinfinityNeighborhood( D2<SBasis > const &f, double tol, D2<Piecewise<SBasis> > &topside, D2<Piecewise<SBasis> > &botside ){ + double signx = ( f[X].at0() > f[X].at1() )? -1 : 1; + double signy = ( f[Y].at0() > f[Y].at1() )? -1 : 1; + + Piecewise<D2<SBasis> > top, bot; + top = Piecewise<D2<SBasis> > (f); + top.cuts.insert( top.cuts.end(), 2); + top.segs.insert( top.segs.end(), D2<SBasis>(SBasis(Linear( f[X].at1(), f[X].at1()+2*tol*signx)), + SBasis(Linear( f[Y].at1() )) )); + bot = Piecewise<D2<SBasis> >(f); + bot.cuts.insert( bot.cuts.begin(), - 1 ); + bot.segs.insert( bot.segs.begin(), D2<SBasis>(SBasis(Linear( f[X].at0()-2*tol*signx, f[X].at0())), + SBasis(Linear( f[Y].at0() )) )); + top += Point(-tol*signx, tol); + bot += Point( tol*signx, -tol); + + if ( signy < 0 ){ + std::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]. + */ +static +void computeLinfinityNeighborhood( Piecewise<SBasis> const &f, double tol, Piecewise<SBasis> &top, Piecewise<SBasis> &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()) { + std::swap(top, bot); + top += 2 * tol; + bot -= 2 * tol; + } +} + +/* + * 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. + */ +std::vector<Interval> monotonicSplit(D2<SBasis> const &p){ + std::vector<Interval> result; + + D2<SBasis> v = derivative(p); + + std::vector<double> someroots; + std::vector<double> cuts (2,0.); + cuts[1] = 1.; + + someroots = roots(v[X]); + cuts.insert( cuts.end(), someroots.begin(), someroots.end() ); + + someroots = roots(v[Y]); + cuts.insert( cuts.end(), someroots.begin(), someroots.end() ); + + //we could split in the middle to avoid computing roots again... + someroots = roots(v[X]-v[Y]); + cuts.insert( cuts.end(), someroots.begin(), someroots.end() ); + + someroots = roots(v[X]+v[Y]); + cuts.insert( cuts.end(), someroots.begin(), someroots.end() ); + + sort(cuts.begin(),cuts.end()); + unique(cuts.begin(), cuts.end() ); + + for (unsigned i=1; i<cuts.size(); i++){ + result.push_back( Interval( cuts[i-1], cuts[i] ) ); + } + return result; +} + +//std::vector<Interval> level_set( D2<SBasis> const &f, Rect region){ +// std::vector<Interval> x_in_reg = level_set( f[X], region[X] ); +// std::vector<Interval> y_in_reg = level_set( f[Y], region[Y] ); +// std::vector<Interval> result = intersect ( x_in_reg, y_in_reg ); +// return result; +//} + +/*TODO: remove this!!! + * the minimum would be to move it to piecewise.h but this would be stupid. + * The best would be to let 'compose' be aware of extension modes (constant, linear, polynomial..) + * (I think the extension modes (at start and end) should be properties of the pwsb). + */ +static +void prolongateByConstants( Piecewise<SBasis> &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() ) ); +} + +static +bool compareIntersectionsTimesX( SmashIntersection const &inter1, SmashIntersection const &inter2 ){ + return inter1.times[X].min() < inter2.times[Y].min(); +} +/*Fuse contiguous intersection domains + * + */ +static +void cleanup_and_fuse( std::vector<SmashIntersection> &inters ){ + std::sort( inters.begin(), inters.end(), compareIntersectionsTimesX); + for (unsigned i=0; i < inters.size(); i++ ){ + for (unsigned j=i+1; j < inters.size() && inters[i].times[X].intersects( inters[j].times[X]) ; j++ ){ + if (inters[i].times[Y].intersects( inters[j].times[Y] ) ){ + inters[i].times.unionWith(inters[j].times); + inters[i].bbox.unionWith(inters[j].bbox); + inters.erase( inters.begin() + j ); + } + } + } +} + +/* Computes the intersection of two sets given as (ordered) union intervals. + */ +static +std::vector<Interval> intersect( std::vector<Interval> const &a, std::vector<Interval> const &b){ + std::vector<Interval> result; + //TODO: use order to optimize this! + for (auto i : a){ + for (auto j : b){ + OptInterval c( i ); + c &= j; + if ( c ) { + result.push_back( *c ); + } + } + } + return result; +} + +/* Returns the intervals over which the curves are in the + * tol-neighborhood one of the other for the L_\infty metric. + * WARNING: each curve is supposed to be a graph over x or y axis + * (but not necessarily the same axis for both) and the smaller + * the slope the better (typically <=45°). + */ +std::vector<SmashIntersection> monotonic_smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){ + using std::swap; + + // a and b or X and Y may have to be exchanged, so make local copies. + D2<SBasis> aa = a; + D2<SBasis> bb = b; + bool swapresult = false; + bool swapcoord = false;//debug only! + + //if the (enlarged) bounding boxes don't intersect, stop. + OptRect abounds = bounds_fast( a ); + OptRect bbounds = bounds_fast( b ); + if ( !abounds || !bbounds ) return std::vector<SmashIntersection>(); + abounds->expandBy(tol); + if ( !(abounds->intersects(*bbounds))){ + return std::vector<SmashIntersection>(); + } + + //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( std::abs((*dabounds)[X].max()), std::abs((*dabounds)[X].min()) ); + double dymin = std::min( std::abs((*dabounds)[Y].max()), std::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<SBasis>( aa[Y], aa[X] ); + bb = D2<SBasis>( 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<SBasis> 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<SBasis> 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. + std::vector<Interval> bx_in_ax_range = level_set(bb[X], ax_range ); + + // find times when bb is in the neighborhood of aa. + std::vector<Interval> tbs; + for (auto & i : bx_in_ax_range){ + D2<Piecewise<SBasis> > bb_in; + bb_in[X] = Piecewise<SBasis> ( portion( bb[X], i ) ); + bb_in[Y] = Piecewise<SBasis> ( portion( bb[Y], i) ); + bb_in[X].setDomain( i ); + bb_in[Y].setDomain( i ); + + Piecewise<SBasis> h; + Interval level; + h = bb_in[Y] - compose( top_ay, bb_in[X] ); + level = Interval( -infinity(), 0 ); + std::vector<Interval> rts_lo = level_set( h, level); + h = bb_in[Y] - compose( bot_ay, bb_in[X] ); + level = Interval( 0, infinity()); + std::vector<Interval> rts_hi = level_set( h, level); + + std::vector<Interval> rts = intersect( rts_lo, rts_hi ); + tbs.insert(tbs.end(), rts.begin(), rts.end() ); + } + + std::vector<SmashIntersection> result(tbs.size(), SmashIntersection()); + + /* 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<tbs.size(); j++){ + result[j].times[Y] = tbs[j]; + std::vector<Interval> tas; + //TODO: replace this by some option in the "compose(pw,pw)" method! + Piecewise<SBasis> fat_y_of_x = y_of_x; + prolongateByConstants( fat_y_of_x, 100*(1+tol) ); + + D2<Piecewise<SBasis> > top_b, bot_b; + D2<SBasis> bbj = portion( bb, tbs[j] ); + computeLinfinityNeighborhood( bbj, tol, top_b, bot_b ); + + Piecewise<SBasis> h; + Interval level; + h = top_b[Y] - compose( fat_y_of_x, top_b[X] ); + level = Interval( +infinity(), 0 ); + std::vector<Interval> rts_top = level_set( h, level); + for (auto & idx : rts_top){ + idx = Interval( top_b[X].valueAt( idx.min() ), + top_b[X].valueAt( 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<Interval> rts_bot = level_set( h, level); + for (auto & idx : rts_bot){ + idx = Interval( bot_b[X].valueAt( idx.min() ), + bot_b[X].valueAt( idx.max() ) ); + } + assert( rts_bot.size() == 1 ); + rts_top = intersect( rts_top, rts_bot ); + assert (rts_top.size() == 1); + Interval x_dom = rts_top[0]; + + if ( x_dom.max() <= x_range_strict.min() ){ + tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 0 : 1 ) ); + }else if ( x_dom.min() >= x_range_strict.max() ){ + tas.push_back( Interval ( ( aa[X].at0() < aa[X].at1() ) ? 1 : 0 ) ); + }else{ + tas = level_set(aa[X], x_dom ); + } + assert( tas.size()==1 ); + result[j].times[X] = tas.front(); + + result[j].bbox = Rect( bbj.at0(), bbj.at1() ); + Interval y_dom( aa[Y](result[j].times[X].min()), aa[Y](result[j].times[X].max()) ); + result[j].bbox.unionWith( Rect( x_dom, y_dom ) ); + } + + if (swapresult) { + for (auto & i : result){ + swap( i.times[X], i.times[Y]); + } + } + if (swapcoord) { + for (auto & i : result){ + swap( i.bbox[X], i.bbox[Y] ); + } + } + + //TODO: cleanup result? fuse contiguous intersections? + return result; +} + +std::vector<SmashIntersection> smash_intersect( D2<SBasis> const &a, D2<SBasis> const &b, double tol){ + std::vector<SmashIntersection> result; + + std::vector<Interval> acuts = monotonicSplit(a); + std::vector<Interval> bcuts = monotonicSplit(b); + for (auto & acut : acuts){ + D2<SBasis> ai = portion( a, acut); + for (auto & bcut : bcuts){ + D2<SBasis> bj = portion( b, bcut); + std::vector<SmashIntersection> ai_cap_bj = monotonic_smash_intersect( ai, bj, tol ); + for (auto & k : ai_cap_bj){ + k.times[X] = k.times[X] * acut.extent() + acut.min(); + k.times[Y] = k.times[Y] * bcut.extent() + bcut.min(); + } + result.insert( result.end(), ai_cap_bj.begin(), ai_cap_bj.end() ); + } + } + cleanup_and_fuse( result ); + return result; +} + +} + +/* + 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 : |