Adding upstream version 1.3.upstream/1.3upstream

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

1 files changed, 1135 insertions, 0 deletions

diff --git a/src/toys/sweeper.cpp b/src/toys/sweeper.cpp
new file mode 100644
index 0000000..7dae586
--- /dev/null
+++ b/src/toys/sweeper.cpp
@@ -0,0 +1,1135 @@
+#include <iostream>
+#include <2geom/path.h>
+#include <2geom/basic-intersection.h>
+#include <2geom/pathvector.h>
+#include <2geom/exception.h>
+
+#include <cstdlib>
+#include <cstdio>
+#include <set>
+#include <vector>
+#include <algorithm>
+
+#include <limits.h>
+#define NULL_IDX UINT_MAX
+
+#include <2geom/orphan-code/intersection-by-smashing.h>
+#include "../2geom/orphan-code/intersection-by-smashing.cpp"
+
+using namespace Geom;
+using namespace std;
+
+
+/*
+The sweeper class takes a PathVector as input and generates "events" to let clients construct the relevant graph.
+
+The basic strategy is the following:
+The path is split into "tiles": a tile consists in 2 boxes related by a (monotonic) curve.
+
+The tiles are created at the very beginning, using a sweep, but *no care* is taken to topology
+information at this step! All the boxes of all the tiles are then enlarged so that they are
+either equal or disjoint.
+[TODO: we should look for curves traversing boxes, split them and repeat the process...]
+
+The sweeper maintains a virtual sweepline, that is the limit of the "known area". The tiles can have 2 states:
+open if they have one end in the known area, and one in the unknown, closed otherwise.
+[TODO: open/close should belong to tiles pointers, not tiles...]
+
+The sorted list of open tiles intersecting the sweep line is called the "context".
+*!*WARNING*!*: because the sweep line is not straight, closed tiles can still be in the context!!
+they can only be removed once the end of the last box is reached.
+
+The events are changes in the context when the sweep line crosses boxes.
+They are obtained by sorting the tiles according to one or the other of theire end boxes depending
+on the open/close state.
+
+A "big" event happens when the sweep line reaches a new 'box'. After such a "big" event, the sweep
+line goes round the new box along it's 3 other sides.
+N.B.: in an ideal world, all tiles ending at one box would be on one side, all the tiles starting
+there on the other. Unfortunately, because we have boxes as vertices, things are not that nice:
+open/closed tiles can appear in any order around a vertex, even in the monotonic case(!). Morover,
+our fat vertices have a non zero "duration", during which many things can happen: this is why we
+have to keep closed edges in the context until both ends of theire boxes are reached...
+
+
+To keep things uniform, such "big" events are split into elementary ones: opening/closing of a single
+edge. One such event is generated for each tile around the current 'box', in CCW order (geometrically,
+the sweepline is deformed in a neighborhood of the box to go round it for a certain amount, enter the
+box and come back inside the box; the piece inside the box is a "virtual edge" that is not added for
+good but that we keep track of). The event knows if it's the last one in such a sequence, so that the
+client knows when to do the additional work required to "close" the vertex construction. Hmmm. It's
+hard to explain the moves without a drawing here...(see sweep.svg in the doc dir). There are
+
+*Closings: insert a new the relevant tile in the context with a "exit" flag.
+
+*Openings: insert a new the relevant tile in the context with a "entry" flag.
+
+At the end of a box, the relevant exit/entries are purged from the context.
+
+
+N.B. I doubt we can do boolops without building the full graph, i.e. having different clients to obtain
+different outputs. So splitting sweeper/grpah builder is maybe not so relevant w/r to functionality
+(only code organization).
+*/
+
+
+//TODO: decline intersections algorithms for each kind of curves...
+//TODO: write an intersector that can work on sub domains.
+//TODO: factor computation of derivative and the like out.
+std::vector<SmashIntersection> monotonic_smash_intersect( Curve const &a, Interval a_dom,
+			                                         Curve const &b, Interval b_dom, double tol){
+	std::vector<SmashIntersection> result;
+    D2<SBasis> asb = a.toSBasis();
+    asb = portion( asb, a_dom );
+    D2<SBasis> bsb = b.toSBasis();
+    bsb = portion( bsb, b_dom );
+    result = monotonic_smash_intersect(asb, bsb, tol );
+    for (auto & i : result){
+    	i.times[X] *= a_dom.extent();
+    	i.times[X] += a_dom.min();
+    	i.times[Y] *= b_dom.extent();
+    	i.times[Y] += b_dom.min();
+    }
+    return result;
+}
+
+
+
+class Sweeper{
+public:
+
+    //---------------------------
+    // utils...
+    //---------------------------
+
+    //near predicate utilized in process_splits
+    template<typename T>
+    struct NearPredicate {
+        double tol;
+        NearPredicate(double eps):tol(eps){}
+        NearPredicate(){tol = EPSILON;}//???
+        bool operator()(T x, T y) { return are_near(x, y, tol); } };
+
+    // ensures that f and t are elements of a vector, sorts and uniqueifies
+    // also asserts that no values fall outside of f and t
+    // if f is greater than t, the sort is in reverse
+    void process_splits(std::vector<double> &splits, double f, double t, double tol=EPSILON) {
+        //splits.push_back(f);
+        //splits.push_back(t);
+        std::sort(splits.begin(), splits.end());
+        std::vector<double>::iterator end = std::unique(splits.begin(), splits.end(), NearPredicate<double>(tol));
+        splits.resize(end - splits.begin());
+
+        //remove any splits which fall outside t / f
+        while(!splits.empty() && splits.front() < f+tol) splits.erase(splits.begin());
+        splits.insert(splits.begin(), f);
+        //splits[0] = f;
+        while(!splits.empty() && splits.back()  > t-tol) splits.erase(splits.end() - 1);
+        splits.push_back(t);
+        //splits.back() = t;
+    }
+
+    struct IntersectionMinTimeOrder {
+        unsigned which;
+        IntersectionMinTimeOrder (unsigned idx) : which(idx) {}
+        bool operator()(SmashIntersection const &a, SmashIntersection const &b) const {
+        	return a.times[which].min() < b.times[which].min();
+        }
+    };
+
+    // ensures that f and t are elements of a vector, sorts and uniqueifies
+    // also asserts that no values fall outside of f and t
+    // if f is greater than t, the sort is in reverse
+    std::vector<std::pair<Interval, Rect> >
+    process_intersections(std::vector<SmashIntersection> &inters, unsigned which, unsigned tileidx) {
+    	std::vector<std::pair<Interval, Rect> > result;
+    	std::pair<Interval, Rect> apair;
+    	Interval dom ( tiles_data[tileidx].f, tiles_data[tileidx].t );
+    	apair.first  = Interval( dom.min() );
+    	apair.second = tiles_data[tileidx].fbox;
+    	result.push_back( apair );
+
+    	std::sort(inters.begin(), inters.end(), IntersectionMinTimeOrder(which) );
+    	for (auto & inter : inters){
+    		if ( !inter.times[which].intersects( dom ) )//this should never happen.
+    			continue;
+    		if ( result.back().first.intersects( inter.times[which] ) ){
+    			result.back().first.unionWith( inter.times[which] );
+    			result.back().second.unionWith( inter.bbox );
+    		}else{
+    	    	apair.first  = inter.times[which];
+    	    	apair.second = inter.bbox;
+    	    	result.push_back( apair );
+    		}
+    	}
+    	apair.first  = Interval( dom.max() );
+    	apair.second = tiles_data[tileidx].tbox;
+		if ( result.size() > 1 && result.back().first.intersects( apair.first ) ){
+			result.back().first.unionWith( apair.first );
+			result.back().second.unionWith( apair.second );
+		}else{
+	    	result.push_back( apair );
+		}
+		return result;
+    }
+
+
+    //---------------------------
+    // Tiles.
+    //---------------------------
+
+    //A tile is a "light edge": just two boxes, joint by a curve.
+    //it is open iff intersected by the sweepline.
+    class Tile{
+    public:
+        unsigned path;
+        unsigned curve;
+        double f;
+        double t;
+        Rect fbox, tbox;
+        bool reversed;//with respect to sweep direction. Flip f/t instead?
+        bool open;//means sweepline currently cuts it (i.e. one end in the known area, the other in the unknown).
+        int state;//-1: both ends in unknown area, 0:one end in each, 1: both in known area.
+        //Warning: we can not delete a tile immediately when it's past(=closed again), only when the end of it's tbox is!.
+        Rect bbox(){Rect b = fbox; b.unionWith(tbox); return b;}
+        Point min(){return ( bbox().min() ); }
+        Point max(){return ( bbox().max() ); }
+//        Rect cur_box() const {return ((open)^(reversed) ) ? tbox : fbox; }
+        Rect cur_box() const { return ((state>=0)^(reversed) ) ? tbox : fbox; }
+        Rect first_box() const {return ( reversed ) ? tbox : fbox; }
+        Rect last_box() const {return ( reversed ) ? fbox : tbox; }
+    };
+
+    D2<SBasis> tileToSB(Tile const &tile){
+        //TODO: don't convert each time!!!!!!
+        assert( tile.path < paths.size() );
+        assert( tile.curve < paths[tile.path].size() );
+        D2<SBasis> c = paths[tile.path][tile.curve].toSBasis();
+        c = portion( c, Interval( tile.f, tile.t ) );
+        return c;
+    }
+
+    //SweepOrder for Rects or Tiles.
+    class SweepOrder{
+    public:
+        Dim2 dim;
+        SweepOrder(Dim2 d) : dim(d) {}
+        bool operator()(const Rect &a, const Rect &b) const {
+            return Point::LexLessRt(dim)(a.min(), b.min());
+        }
+        bool operator()(const Tile &a, const Tile &b) const {
+            return Point::LexLessRt(dim)(a.cur_box().min(), b.cur_box().min());
+        }
+    };
+
+    class PtrSweepOrder{
+    public:
+        Dim2 dim;
+        std::vector<Tile>::iterator const begin;
+        PtrSweepOrder(std::vector<Tile>::iterator const beg, Dim2 d) : dim(d), begin(beg){}
+        bool operator()(const unsigned a, const unsigned b) const {
+            return Point::LexLessRt(dim)((begin+a)->cur_box().min(), (begin+b)->cur_box().min());
+        }
+    };
+
+
+    //---------------------------
+    // Vertices.
+    //---------------------------
+
+    //A ray is nothing but an edge ending or starting at a given vertex, + some info about when/where it exited a "separating" box;
+    struct Ray{
+    public:
+        unsigned tile;
+        bool centrifuge;//true if the intrinsic orientation of curve points away from the vertex.
+        //exit info:
+        unsigned exit_side;//0:y=min; 1:x=max; 2:y=max; 3:x=min.
+        double exit_place; //x or y value on the exit line.
+        double exit_time;  //exit time on curve.
+        Ray(){tile = NULL_IDX; exit_side = 4;}
+        Ray(unsigned tile_idx,  unsigned s, double p, double t){
+            tile = tile_idx;
+            exit_side =s;
+            exit_place = p;
+            exit_time = t;
+        }
+        Ray(unsigned tile_idx, bool outward){
+            tile = tile_idx;
+            exit_side = 4;
+            centrifuge = outward;
+            exit_time = (centrifuge) ? 2 : -1 ;
+        }
+        void setExitInfo( unsigned side, double place, double time){
+            exit_side = side;
+            exit_place = place;
+            exit_time = time;
+        }
+    };
+
+    class FatVertex : public Rect{
+    public:
+        std::vector<Ray> rays;
+        FatVertex(const Rect &r, unsigned const tile, bool centrifuge) : Rect(r){
+            rays.push_back( Ray(tile, centrifuge) );
+        }
+        FatVertex(Rect r) : Rect(r){}
+        FatVertex() : Rect(){}
+        void erase(unsigned from, unsigned to){
+            unsigned size = to-from;
+            from = from % rays.size();
+            to = from + size;
+
+            if (to >= rays.size() ){
+                to  = to % rays.size();
+                rays.erase( rays.begin()+from, rays.end() );
+                rays.erase( rays.begin(), rays.begin()+to );
+            }else{
+                rays.erase( rays.begin()+from, rays.begin()+to );
+            }
+
+        }
+    };
+
+    //---------------------------
+    // Context related stuff.
+    //---------------------------
+
+    class Event{
+    public:
+        bool opening;//true means an edge is added, otherwise an edge is removed from context.
+        unsigned tile;//which tile to open/close.
+        unsigned insert_at;//where to insert the next tile in the context.
+        //unsigned erase_at;//idx of the tile to close in the context.  = context.find(tile).
+        bool to_be_continued;
+        bool empty(){
+            return  tile==NULL_IDX;
+        }
+        Event(){
+            opening = false;
+            insert_at = 0;
+            //erase_at = 0;
+            tile = NULL_IDX;
+            to_be_continued = false;
+        }
+    };
+
+    void printEvent(Event const &e){
+        std::printf("Event: ");
+        std::printf("%s, ", e.opening?"opening":"closing");
+        std::printf("insert_at:%u, ", e.insert_at);
+        //std::printf("erase_at:%u, ", e.erase_at);
+        std::printf("tile:%u.\n", e.tile);
+    }
+
+    class Context : public std::vector<std::pair<unsigned,bool> >{//first = tile, second = true if it's a birth (+).
+    public:
+        Point last_pos;
+        FatVertex pending_vertex;
+        Event pending_event;
+
+        unsigned find(unsigned const tile, bool positive_only=false){
+            for (unsigned i=0; i<size(); i++){
+                if ( (*this)[i].first == tile ){
+                    if ( (*this)[i].second || !positive_only ) return i;
+                }
+            }
+            return (*this).size();
+        }
+    };
+    void printContext(){
+        std::printf("context:[");
+        for (unsigned i=0; i<context.size(); i++){
+            unsigned tile = context[i].first;
+            assert( tile<tiles_data.size() );
+            std::printf(" %s%u%s", (tiles_data[ tile ].reversed)?"-":"+", tile, (context[i].second)?"o":"c");
+//            assert( context[i].second || !tiles_data[ tile ].open);
+            assert( context[i].second || tiles_data[ tile ].state==1);
+        }
+        std::printf("]\n");
+    }
+
+
+
+    //----
+    //This is the heart of it all!! Take particular care to non linear sweep line...
+    //----
+    //Given a point on the sweep line, (supposed to be the min() of a vertex not yet connected to the already known part),
+    //find the first edge "after" it in the context. Pretty tricky atm :-(
+    //TODO: implement this as a lower_bound (?).
+
+    unsigned contextRanking(Point const &pt){
+
+//        std::printf("contextRanking:------------------------------------\n");
+
+        unsigned rank = context.size();
+        std::vector<unsigned> unmatched_closed_tiles = std::vector<unsigned>();
+
+//        std::printf("Scan context.\n");
+
+        for (unsigned i=0; i<context.size(); i++){
+
+            unsigned tile_idx = context[i].first;
+            assert(  tile_idx < tiles_data.size() );
+//            std::printf("testing %u (e=%u),", i, tile_idx);
+
+            Tile tile = tiles_data[tile_idx];
+            assert( tile.state >= 0 );
+
+            //if the tile is open (i.e. not both ends in the known area) and point is below/above the tile's bbox:
+            if ( tile.state == 0 ){
+//                std::printf("opened tile, ");
+                if (pt[1-dim] < tile.min()[1-dim] ) {
+//                    printContext();
+//                    std::printf("below bbox %u!\n", i);
+                    rank = i;
+                    break;
+                }
+                if (pt[1-dim] > tile.max()[1-dim] ){
+//                    std::printf("above bbox %u!\n", i);
+                    continue;
+                }
+                
+                //TODO: don't convert each time!!!!!!
+                D2<SBasis> c = tileToSB( tile );
+            
+                std::vector<double> times = roots(c[dim]-pt[dim]);
+                if (times.size()==0){
+                    assert( tile.first_box()[dim].contains(pt[dim]) );
+                    if ( pt[1-dim] < tile.first_box()[1-dim].min() ){
+//                        std::printf("open+hit box %u!\n", i);
+                        rank = i;
+                        break;
+                    }else{
+                        continue;
+                    }
+                }
+                if ( pt[1-dim] < c[1-dim](times.front()) ){
+//                    std::printf("open+hit curve %u!\n", i);
+                    rank = i;
+                    break;
+                }            
+            }
+
+//            std::printf("closed tile, ");
+
+
+            //At this point, the tile is closed (i.e. both ends are in the known area)
+            //Such tiles do 'nested parens' like travels in the unknown area.
+            //We are interested in the second occurrence only (to give a chance to open tiles to exist in between).
+            if ( unmatched_closed_tiles.size()==0 || tile_idx != unmatched_closed_tiles.back() ){
+                unmatched_closed_tiles.push_back( tile_idx );
+//                std::printf("open paren %u\n",tile_idx);
+                continue;
+            }
+            unmatched_closed_tiles.pop_back();
+
+//            std::printf("close paren, ");
+            
+            if ( !tile.bbox().contains( pt ) ){
+                continue;
+            }
+            
+//            std::printf("in bbox, ");
+
+            //At least one of fbox[dim], tbox[dim] has to contain the pt[dim]: assert it?
+
+            //Find intersection with the hline(vline if dim=Y) through the point
+            double hit_place;
+            //TODO: don't convert each time!!!!!!
+           D2<SBasis> c = tileToSB( tile );
+           std::vector<double> times = roots(c[1-dim]-pt[1-dim]);
+           if ( times.size()>0 ){
+//               std::printf("hit curve,");
+               hit_place = c[dim](times.front());
+           }else{
+//               std::printf("hit box, ");
+               //if there was no intersection, the line went through the first_box
+               assert( tile.first_box()[1-dim].contains(pt[1-dim]) );
+               continue;
+           }
+            
+            if ( pt[dim] > hit_place ){
+//                std::printf("wrong side, ");
+                continue;
+            }
+//            std::printf("good side, ");
+            rank = i;
+            break;
+        }
+
+//        std::printf("rank %u.\n", rank);
+//        printContext();
+        assert( rank<=tiles_data.size() );
+        return rank;
+    }
+
+    //TODO: optimize this.
+    //it's done the slow way for debugging purpose...
+    void purgeDeadTiles(){
+        //std::printf("purge ");
+        //printContext();
+        for (unsigned i=0; i<context.size(); i++){
+            assert( context[i].first<tiles_data.size() );
+            Tile tile = tiles_data[context[i].first];
+            if (tile.state==1 && Point::LexLessRt(dim)( tile.fbox.max(), context.last_pos ) && Point::LexLessRt(dim)( tile.tbox.max(), context.last_pos ) ){
+//            if (!tile.open && Point::LexLessRt(dim)( tile.fbox.max(), context.last_pos ) && Point::LexLessRt(dim)( tile.tbox.max(), context.last_pos ) ){
+                unsigned j;
+                for (j=i+1; j<context.size() && context[j].first != context[i].first; j++){}
+                assert ( j < context.size() );
+                if ( context[j].first == context[i].first){
+                    assert ( context[j].second == !context[i].second );
+                    context.erase(context.begin()+j);
+                    context.erase(context.begin()+i);
+//                    printContext();
+                    i--;
+                }
+            }
+        }
+        return;
+    }
+
+    void applyEvent(Event event){
+//        std::printf("Apply event : ");
+        if(event.empty()){
+//            std::printf("empty event!\n");
+            return;
+        }
+
+//        printEvent(event);
+//        std::printf("    old ");
+//        printContext();
+
+        assert ( context.begin() + event.insert_at <= context.end() );
+
+        if (!event.opening){
+//             unsigned idx = event.erase_at;
+//             assert( idx == context.find(event.tile) );
+//             assert( context[idx].first == event.tile);
+            tiles_data[event.tile].open = false;
+            tiles_data[event.tile].state = 1;
+            //context.erase(context.begin()+idx);
+            unsigned idx = event.insert_at;
+            context.insert(context.begin()+idx, std::pair<unsigned, bool>(event.tile, false) );
+        }else{
+            unsigned idx = event.insert_at;
+            tiles_data[event.tile].open = true;
+            tiles_data[event.tile].state = 0;
+            context.insert(context.begin()+idx, std::pair<unsigned, bool>(event.tile, true) );
+            sortTiles();
+        }
+        context.last_pos = context.pending_vertex.min();
+        context.last_pos[1-dim] = context.pending_vertex.max()[1-dim];
+
+//        std::printf("    new ");
+//        printContext();
+//        std::printf("\n");
+        //context.pending_event = Event();is this a good idea?
+    }
+
+
+
+    //---------------------------
+    // Sweeper.
+    //---------------------------
+
+    PathVector paths;
+    std::vector<Tile> tiles_data;
+    std::vector<unsigned> tiles;
+    std::vector<Rect> vtxboxes;
+    Context context;
+    double tol;
+    Dim2 dim;
+
+
+    //-------------------------------
+    //-- Tiles preparation.
+    //-------------------------------
+
+    //split input paths into monotonic pieces...
+    void createMonotonicTiles(){
+        for ( unsigned i=0; i<paths.size(); i++){
+            for ( unsigned j=0; j<paths[i].size(); j++){
+                //find the points where slope is 0°, 45°, 90°, 135°...
+                D2<SBasis> deriv = derivative( paths[i][j].toSBasis() );
+                std::vector<double> splits0 = roots( deriv[X] );
+                std::vector<double> splits90 = roots( deriv[Y] );
+                std::vector<double> splits45 = roots( deriv[X]- deriv[Y] );
+                std::vector<double> splits135 = roots( deriv[X] + deriv[Y] );
+                std::vector<double> splits;
+                splits.insert(splits.begin(), splits0.begin(), splits0.end() );
+                splits.insert(splits.begin(), splits90.begin(), splits90.end() );
+                splits.insert(splits.begin(), splits45.begin(), splits45.end() );
+                splits.insert(splits.begin(), splits135.begin(), splits135.end() );
+                process_splits(splits,0,1);
+
+                for(unsigned k = 1; k < splits.size(); k++){
+                	Tile tile;
+                    tile.path = i;
+                    tile.curve = j;
+                    tile.f = splits[k-1];
+                    tile.t = splits[k];
+                    //TODO: use meaningful tolerance here!!
+                    Point fp = paths[i][j].pointAt(tile.f);
+                    Point tp = paths[i][j].pointAt(tile.t);
+                    tile.fbox = Rect(fp, fp );
+                    tile.tbox = Rect(tp, tp );
+                    tile.open = false;
+                    tile.state = -1;
+                    tile.reversed = Point::LexLessRt(dim)(tp, fp);
+
+                    tiles_data.push_back(tile);
+                }
+            }
+        }
+        std::sort(tiles_data.begin(), tiles_data.end(), SweepOrder(dim) );
+    }
+
+    void splitTile(unsigned i, double t, double tolerance=0, bool sort = true){
+        assert( i<tiles_data.size() );
+        Tile newtile = tiles_data[i];
+        assert( newtile.f < t && t < newtile.t );
+        newtile.f = t;
+        //newtile.fbox = fatPoint(paths[newtile.path][newtile.curve].pointAt(t), tolerance );
+        Point p = paths[newtile.path][newtile.curve].pointAt(t);
+        newtile.fbox = Rect(p, p);
+        newtile.fbox.expandBy( tolerance );
+        tiles_data[i].tbox = newtile.fbox;
+        tiles_data[i].t = t;
+        tiles_data.insert(tiles_data.begin()+i+1, newtile);
+        if (sort)
+            std::sort(tiles_data.begin()+i+1, tiles_data.end(), SweepOrder(dim) );
+    }
+    void splitTile(unsigned i, SmashIntersection inter, unsigned which,bool sort = true){
+    	double t = inter.times[which].middle();
+    	assert( i<tiles_data.size() );
+        Tile newtile = tiles_data[i];
+        assert( newtile.f < t && t < newtile.t );
+        newtile.f = t;
+        newtile.fbox = inter.bbox;
+        tiles_data[i].tbox = newtile.fbox;
+        tiles_data[i].t = t;
+        tiles_data.insert(tiles_data.begin()+i+1, newtile);
+        if (sort)
+            std::sort(tiles_data.begin()+i+1, tiles_data.end(), SweepOrder(dim) );
+    }
+#if 0
+    void splitTile(unsigned i, std::vector<double> const &times, double tolerance=0, bool sort = true){
+        if ( times.size()<3 ) return;
+        assert( i<tiles_data.size() );
+        std::vector<Tile> pieces ( times.size()-2, tiles_data[i] );
+        Rect prevbox = tiles_data[i].fbox;
+        for (unsigned k=0; k < times.size()-2; k++){
+                 pieces[k].f = times[k];
+                 pieces[k].t = times[k+1];
+                 pieces[k].fbox = prevbox;
+                 //TODO: use relevant precision here.
+                 prevbox = fatPoint(paths[tiles_data[i].path][tiles_data[i].curve].pointAt(times[k+1]), tolerance );
+                 pieces[k].tbox = prevbox;
+        }
+        tiles_data.insert(tiles_data.begin()+i, pieces.begin(), pieces.end() );
+        unsigned newi = i + times.size()-2;
+        assert( newi<tiles_data.size() );
+        assert( newi>=1 );
+        tiles_data[newi].f    = tiles_data[newi-1].t;
+        tiles_data[newi].fbox = tiles_data[newi-1].tbox;
+
+        if (sort)
+            std::sort(tiles_data.begin()+i, tiles_data.end(), SweepOrder(dim) );
+    }
+#else
+    void splitTile(unsigned i, std::vector<std::pair<Interval,Rect> > const &cuts, bool sort = true){
+        assert ( cuts.size() >= 2 );
+        assert( i<tiles_data.size() );
+        std::vector<Tile> pieces ( cuts.size()-1, tiles_data[i] );
+        for (unsigned k=1; k+1 < cuts.size(); k++){
+        	pieces[k-1].t = cuts[k].first.middle();
+        	pieces[k  ].f = cuts[k].first.middle();
+        	pieces[k-1].tbox = cuts[k].second;
+        	pieces[k  ].fbox = cuts[k].second;
+        }
+    	pieces.front().fbox.unionWith( cuts[0].second );
+    	pieces.back().tbox.unionWith( cuts.back().second );
+
+        tiles_data.insert(tiles_data.begin()+i, pieces.begin(), pieces.end()-1 );
+        unsigned newi = i + cuts.size()-2;
+        assert( newi < tiles_data.size() );
+        tiles_data[newi] = pieces.back();
+
+        if (sort)
+            std::sort(tiles_data.begin()+i, tiles_data.end(), SweepOrder(dim) );
+    }
+#endif
+
+    //TODO: maybe not optimal. For a fully optimized sweep, it would be nice to have
+    //an efficient way to way find *only the first* intersection (in sweep direction)...
+    void splitIntersectingTiles(){
+        //make sure it is sorted, but should be ok. (remove sorting at the end of monotonic tiles creation?
+        std::sort(tiles_data.begin(), tiles_data.end(), SweepOrder(dim) );
+
+//        std::printf("\nFind intersections: tiles_data.size():%u\n", tiles_data.size() );
+
+        for (unsigned i=0; i+1<tiles_data.size(); i++){
+            //std::printf("\ni=%u (%u([%f,%f]))\n", i, tiles_data[i].curve, tiles_data[i].f, tiles_data[i].t );
+            std::vector<SmashIntersection> inters_on_i;
+            for (unsigned j=i+1; j<tiles_data.size(); j++){
+                //std::printf("  j=%u (%u)\n", j,tiles_data[j].curve );
+                if ( Point::LexLessRt(dim)(tiles_data[i].max(), tiles_data[j].min()) ) break;
+
+                unsigned pi = tiles_data[i].path;
+                unsigned ci = tiles_data[i].curve;
+                unsigned pj = tiles_data[j].path;
+                unsigned cj = tiles_data[j].curve;
+                std::vector<SmashIntersection> intersections;
+
+                intersections = monotonic_smash_intersect(paths[pi][ci], Interval(tiles_data[i].f, tiles_data[i].t),
+                                                          paths[pj][cj], Interval(tiles_data[j].f, tiles_data[j].t), tol );
+                inters_on_i.insert( inters_on_i.end(), intersections.begin(), intersections.end() );
+                std::vector<std::pair<Interval, Rect> > cuts = process_intersections(intersections, 1, j);
+
+//                std::printf("  >|%u/%u|=%u. times_j:%u", i, j, crossings.size(),  times_j.size() );
+
+                splitTile(j, cuts, false);
+                j += cuts.size()-2;
+            }
+
+            //process_splits(times_i, tiles_data[i].f, tiles_data[i].t);
+            //assert(times_i.size()>=2);
+            //splitTile(i, times_i, tol, false);
+            //i+=times_i.size()-2;
+            std::vector<std::pair<Interval, Rect> > cuts_on_i = process_intersections(inters_on_i, 0, i);
+            splitTile(i, cuts_on_i, false);
+            i += cuts_on_i.size()-2;
+            //std::printf("new i:%u, tiles_data: %u\n",i ,tiles_data.size());
+            //std::sort(tiles_data.begin()+i+1, tiles_data.end(), SweepOrder(dim) );
+            std::sort(tiles_data.begin()+i+1, tiles_data.end(), SweepOrder(dim) );
+        }
+        //this last sorting should be useless!!
+        std::sort(tiles_data.begin(), tiles_data.end(), SweepOrder(dim) );
+    }
+
+    void sortTiles(){
+        std::sort(tiles.begin(), tiles.end(), PtrSweepOrder(tiles_data.begin(), dim) );
+    }
+
+
+    //-------------------------------
+    //-- Vertices boxes cookup.
+    //-------------------------------
+
+    void fuseInsert(const Rect &b,  std::vector<Rect> &boxes, Dim2 dim){
+        //TODO: this can be optimized...
+        for (unsigned i=0; i<boxes.size(); i++){
+            if ( Point::LexLessRt(dim)( b.max(), boxes[i].min() ) ) break;
+            if ( b.intersects( boxes[i] ) ){
+                Rect bigb = b;
+                bigb.unionWith( boxes[i] );
+                boxes.erase( boxes.begin()+i );
+                fuseInsert( bigb, boxes, dim);
+                return;
+            }
+        }
+        std::vector<Rect>::iterator pos = std::lower_bound(boxes.begin(), boxes.end(), b, SweepOrder(dim) );
+        boxes.insert( pos, b );
+    }
+
+    //debug only!!
+    bool isContained(Rect b,   std::vector<Rect> const &boxes ){
+        for (const auto & boxe : boxes){
+            if ( boxe.contains(b) ) return true;
+        }
+        return false;
+    }
+
+    //Collect vertex boxes. Fuse overlapping ones.
+    //NB: enlarging a vertex may create intersection with already scanned ones...
+    std::vector<Rect> collectBoxes(){
+        std::vector<Rect> ret;
+        for (auto & i : tiles_data){
+            fuseInsert(i.fbox, ret, dim);
+            fuseInsert(i.tbox, ret, dim);
+        }
+        return ret;
+    }
+
+    //enlarge tiles ends to match the vertices bounding boxes.
+    //remove edges fully contained in one vertex bbox.
+    void enlargeTilesEnds(const std::vector<Rect> &boxes ){
+        for (unsigned i=0; i<tiles_data.size(); i++){
+            std::vector<Rect>::const_iterator f_it;
+            f_it = std::lower_bound(boxes.begin(), boxes.end(), tiles_data[i].fbox, SweepOrder(dim) );
+            if ( f_it==boxes.end() ) f_it--;
+            while (!(*f_it).contains(tiles_data[i].fbox) && f_it != boxes.begin()){
+                f_it--;
+            }
+            assert( (*f_it).contains(tiles_data[i].fbox) );
+            tiles_data[i].fbox = *f_it;
+
+            std::vector<Rect>::const_iterator t_it;
+            t_it = std::lower_bound(boxes.begin(), boxes.end(), tiles_data[i].tbox, SweepOrder(dim) );
+            if ( t_it==boxes.end() ) t_it--;
+            while (!(*t_it).contains(tiles_data[i].tbox) && t_it != boxes.begin()){
+                t_it--;
+            }
+            assert( (*t_it).contains(tiles_data[i].tbox) );
+            tiles_data[i].tbox = *t_it;
+
+            //NB: enlarging the ends may swapp their sweep order!!!
+            tiles_data[i].reversed = Point::LexLessRt(dim)( tiles_data[i].tbox.min(), tiles_data[i].fbox.min());
+
+            if ( f_it==t_it ){
+                tiles_data.erase(tiles_data.begin()+i);
+                i-=1;
+            }
+        }
+    }
+
+    //Make sure tiles stop at vertices. Split them if needed.
+    //Returns true if at least one tile was split.
+    bool splitTilesThroughFatPoints(std::vector<Rect> &boxes ){
+        
+        std::sort(tiles.begin(), tiles.end(), PtrSweepOrder(tiles_data.begin(), dim) );
+        std::sort(boxes.begin(), boxes.end(), SweepOrder(dim) );
+
+        bool result = false;
+        for (unsigned i=0; i<tiles_data.size(); i++){
+            for (auto & boxe : boxes){
+                if ( Point::LexLessRt(dim)( tiles_data[i].max(), boxe.min()) ) break;
+                if ( Point::LexLessRt(dim)( boxe.max(), tiles_data[i].min()) ) continue;
+                if ( !boxe.intersects( tiles_data[i].bbox() ) ) continue;
+                if ( tiles_data[i].fbox.intersects( boxe ) ) continue;
+                if ( tiles_data[i].tbox.intersects( boxe ) ) continue;
+
+                //at this point box[k] intersects the curve bbox away from the fbox and tbox.
+
+                D2<SBasis> c = tileToSB( tiles_data[i] );
+//----------> use level-set!!
+                for (unsigned corner=0; corner<4; corner++){
+                    unsigned D = corner % 2;
+                    double val = boxe.corner(corner)[D];
+                    std::vector<double> times = roots( c[D] - val );
+                    if ( times.size()>0 ){
+                        double t = lerp( times.front(), tiles_data[i].f, tiles_data[i].t );
+                        double hit_place = c[1-D](times.front());
+                        if ( boxe[1-D].contains(hit_place) ){
+                            result = true;
+                            //taking a point on the boundary is dangerous!! 
+                            //Either use >0 tolerance here, or find 2 intersection points and split in between.
+                            splitTile( i, t, tol, false);
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+        return result;
+    }
+
+
+
+
+
+    //TODO: rewrite all this!...
+    //-------------------------------------------------------------------------------------------
+    //-------------------------------------------------------------------------------------------
+    //-------------------------------------------------------------------------------------------
+    //-------------------------------
+    //-- ccw Sorting of rays around a vertex.
+    //-------------------------------
+    //-------------------------------------------------------------------------------------------
+    //-------------------------------------------------------------------------------------------
+    //-------------------------------------------------------------------------------------------
+    //returns an (infinite) rect around "a" separating it from "b". Nota: 3 sides are infinite!
+    //TODO: place the cut where there is most space...
+    OptRect separate(Rect const &a, Rect const &b){
+        Rect ret ( Interval( -infinity(), infinity() ) , Interval(-infinity(), infinity() ) );
+        double gap = 0;
+        unsigned dir = 4;
+        if (b[X].min() - a[X].max() > gap){
+            gap = b[X].min() - a[X].max();
+            dir = 0;
+        }
+        if (a[X].min() - b[X].max() > gap){
+            gap = a[X].min() - b[X].max();
+            dir = 1;
+        }
+        if (b[Y].min() - a[Y].max() > gap){
+            gap = b[Y].min() - a[Y].max();
+            dir = 2;
+        }
+        if (a[Y].min() - b[Y].max() > gap){
+            gap = a[Y].min() - b[Y].max();
+            dir = 3;
+        }
+        switch (dir) {
+            case 0: ret[X].setMax(( a.max()[X] + b.min()[X] )/ 2); break;
+            case 1: ret[X].setMin(( b.max()[X] + a.min()[X] )/ 2); break;
+            case 2: ret[Y].setMax(( a.max()[Y] + b.min()[Y] )/ 2); break;
+            case 3: ret[Y].setMin(( b.max()[Y] + a.min()[Y] )/ 2); break;
+            case 4: return OptRect();
+        }
+        return OptRect(ret);
+    }
+
+    //Find 4 lines (returned as a Rect sides) that cut all the rays (=edges). *!* some side might be infinite.
+    OptRect isolateVertex(Rect const &box){
+        OptRect sep ( Interval( -infinity(), infinity() ) , Interval(-infinity(), infinity() ) );
+        //separate this vertex from the others. Find a better way.
+        for (auto & vtxboxe : vtxboxes){
+            if ( Point::LexLessRt(dim)( sep->max(), vtxboxe.min() ) ){
+                break;
+            }
+            if ( vtxboxe!=box ){//&& !vtxboxes[i].intersects(box) ){
+                OptRect sepi = separate(box, vtxboxe);
+                if ( sep && sepi ){
+                    sep = intersect(*sep, *sepi);
+                }else{
+                    std::cout<<"box="<<box<<"\n";
+                    std::cout<<"vtxboxes[i]="<<vtxboxe<<"\n";
+                    assert(sepi);
+                }
+            }
+        }
+        if (!sep) THROW_EXCEPTION("Invalid intersection data.");
+        return sep;
+    }
+
+    //TODO: argh... rewrite to have "dim"=min first place.
+    struct ExitPoint{
+    public:
+        unsigned side; //0:y=min; 1:x=max; 2:y=max; 3:x=min.
+        double place; //x or y value on the exit line.
+        unsigned ray_idx;
+        double time; //exit time on curve.
+        ExitPoint(){}
+        ExitPoint(unsigned s, double p, unsigned r, double t){
+            side =s;
+            place = p;
+            ray_idx = r;
+            time = t;
+        }
+    };
+
+
+    class ExitOrder{
+    public:
+        bool operator()(Ray a, Ray b) const {
+            if ( a.exit_side < b.exit_side ) return true;
+            if ( a.exit_side > b.exit_side ) return false;
+            if ( a.exit_side <= 1) {
+                return ( a.exit_place < b.exit_place );
+            }
+            return ( a.exit_place > b.exit_place );
+        }
+    };
+
+    void printRay(Ray const &r){
+        std::printf("Ray: tile=%u, centrifuge=%u, side=%u, place=%f\n",
+                    r.tile, r.centrifuge, r.exit_side, r.exit_place);
+    }
+    void printVertex(FatVertex const &v){
+        std::printf("Vertex: [%f,%f]x[%f,%f]\n", v[X].min(),v[X].max(),v[Y].min(),v[Y].max() );
+        for (const auto & ray : v.rays){
+            printRay(ray);
+        }
+    }
+
+    //TODO: use a partial order on input coming from the context + Try quadrant order just in case it's enough.
+    //TODO: use monotonic assumption.
+    void sortRays( FatVertex &v ){
+        OptRect sep = isolateVertex(v);
+
+        for (unsigned i=0; i < v.rays.size(); i++){
+            v.rays[i].centrifuge = (tiles_data[ v.rays[i].tile ].fbox == v);
+            v.rays[i].exit_time = v.rays[i].centrifuge ? 1 : 0 ;
+        }
+
+        for (unsigned i=0; i < v.rays.size(); i++){
+            //TODO: don't convert each time!!!
+            assert( v.rays[i].tile < tiles_data.size() );
+            D2<SBasis> c = tileToSB( tiles_data[ v.rays[i].tile ] );
+
+            for (unsigned side=0; side<4; side++){//scan X or Y direction, on level min or max...
+                double level = sep->corner( side )[1-side%2];
+                if (level != infinity() && level != -infinity() ){
+                    std::vector<double> times = roots(c[1-side%2]-level);
+                    if ( times.size() > 0 ) {
+                        double t;
+                        assert( v.rays[i].tile < tiles_data.size() );
+                        if (tiles_data[ v.rays[i].tile ].fbox == v){
+                            t = times.front();
+                            if ( v.rays[i].exit_side > 3 || v.rays[i].exit_time > t ){
+                                v.rays[i].setExitInfo( side, c[side%2](t), t);
+                            }
+                        }else{
+                            t = times.back();
+                            if ( v.rays[i].exit_side > 3 || v.rays[i].exit_time < t ){
+                                v.rays[i].setExitInfo( side, c[side%2](t), t);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        //Rk: at this point, side == 4 means the edge is contained in the intersection box (?)...;
+        std::sort( v.rays.begin(), v.rays.end(), ExitOrder() );
+    }
+
+
+    //-------------------------------
+    //-- initialize all data.
+    //-------------------------------
+
+    Sweeper(){}
+    Sweeper(PathVector const &input_paths, Dim2 sweep_dir, double tolerance=1e-5){
+        paths = input_paths;//use a ptr...
+        dim = sweep_dir;
+        tol = tolerance;
+
+        //split paths into monotonic tiles
+        createMonotonicTiles();
+
+        //split at tiles intersections
+        splitIntersectingTiles();
+
+        //handle overlapping end boxes/and tiles traversing boxes.
+        do{
+            vtxboxes = collectBoxes();
+        }while ( splitTilesThroughFatPoints(vtxboxes) );
+
+        enlargeTilesEnds(vtxboxes);
+
+        //now create the pointers to the tiles.
+        tiles = std::vector<unsigned>(tiles_data.size(), 0);
+        for (unsigned i=0; i<tiles_data.size(); i++){
+            tiles[i] = i;
+        }
+        sortTiles();
+
+        //initialize the context.
+        if (tiles_data.size()>0){
+            context.clear();
+            context.last_pos = tiles_data[tiles.front()].min();
+            context.last_pos[dim] -= 1;
+            context.pending_vertex = FatVertex();
+            context.pending_event = Event();
+        }
+
+//        std::printf("Sweeper initialized (%u tiles)\n", tiles_data.size());
+    }
+
+
+    //-------------------------------
+    //-- Event walk.
+    //-------------------------------
+
+
+    Event getNextEvent(){
+//        std::printf("getNextEvent():\n");
+
+//        std::printf("initial contex:\n");
+//        printContext();
+        Event old_event = context.pending_event, event;
+//        std::printf("apply old event\n");
+        applyEvent(context.pending_event);
+//        printContext();
+
+        if (context.pending_vertex.rays.size()== 0){
+//            std::printf("cook up a new vertex\n");
+
+            //find the edges at the next vertex.
+            //TODO: implement this as a lower bound!!
+            std::vector<unsigned>::iterator low, high;
+            //Warning: bad looking test, but make sure we advance even in case of 0 width boxes...
+            for ( low = tiles.begin(); low != tiles.end() && 
+                      ( tiles_data[*low].state==1 || Point::LexLessRt(dim)(tiles_data[*low].cur_box().min(), context.last_pos) ); low++){}
+
+            if ( low == tiles.end() ){
+//                std::printf("no more event found\n");
+                return(Event());
+            }
+            Rect pos = tiles_data[ *low ].cur_box();
+            context.last_pos = pos.min();
+            context.last_pos[1-dim] = pos.max()[1-dim];
+
+//            printContext();
+//            std::printf("purgeDeadTiles\n");
+            purgeDeadTiles();
+//            printContext();
+
+            FatVertex v(pos);
+            high = low;
+            do{
+//                v.rays.push_back( Ray(*high, !tiles_data[*high].open) );
+                v.rays.push_back( Ray(*high, tiles_data[*high].state!=0) );
+                high++;
+            }while( high != tiles.end() && tiles_data[ *high ].cur_box()==pos );
+
+//            std::printf("sortRays\n");
+            sortRays(v);
+
+            //Look for an opened tile
+            unsigned i=0;
+
+            for( i=0; i<v.rays.size(); ++i){assert( v.rays[i].tile<tiles_data.size() );}
+//            for( i=0; i<v.rays.size() && !tiles_data[ v.rays[i].tile ].open; ++i){}
+            for( i=0; i<v.rays.size() && tiles_data[ v.rays[i].tile ].state!=0; ++i){}
+
+            //if there are only openings:
+            if (i == v.rays.size() ){
+//                std::printf("only openings!\n");
+                event.insert_at = contextRanking(pos.min());
+            }
+            //if there is at least one closing: make it first, and catch 'insert_at'.
+            else{
+//                std::printf("not only openings\n");
+                if( i > 0 ){
+                    std::vector<Ray> head;
+                    head.assign  ( v.rays.begin(), v.rays.begin()+i );
+                    v.rays.erase ( v.rays.begin(), v.rays.begin()+i );
+                    v.rays.insert( v.rays.end(), head.begin(), head.end());
+                }
+                //assert( tiles_data[ v.rays[0].tile ].open );
+                assert( tiles_data[ v.rays[0].tile ].state==0 );
+                event.insert_at = context.find( v.rays.front().tile )+1;
+//                event.erase_at = context.find( v.rays.front().tile );
+//                std::printf("at least one closing!\n");
+            }
+            context.pending_vertex = v;
+        }else{
+//            std::printf("continue biting exiting vertex\n");
+            event.tile = context.pending_vertex.rays.front().tile;
+            event.insert_at = old_event.insert_at;
+//             if (old_event.opening){
+//                 event.insert_at++;
+//             }else{
+//                 if (old_event.erase_at < old_event.insert_at){
+//                     event.insert_at--;
+//                 }
+//             }
+            event.insert_at++;
+        }
+        event.tile = context.pending_vertex.rays.front().tile;
+//        event.opening = !tiles_data[ event.tile ].open;
+        event.opening = tiles_data[ event.tile ].state!=0;
+        event.to_be_continued = context.pending_vertex.rays.size()>1;
+//        if ( !event.opening ) event.erase_at = context.find(event.tile);
+
+        context.pending_vertex.rays.erase(context.pending_vertex.rays.begin());
+        context.pending_event = event;
+//        printEvent(event);
+        return event;
+    }
+};
+
+
+/*
+  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=4:softtabstop=4:fileencoding=utf-8:textwidth=99 :