summaryrefslogtreecommitdiffstats
path: root/src/3rdparty/2geom/src/toys/sweeper.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/3rdparty/2geom/src/toys/sweeper.cpp')
-rw-r--r--src/3rdparty/2geom/src/toys/sweeper.cpp1135
1 files changed, 1135 insertions, 0 deletions
diff --git a/src/3rdparty/2geom/src/toys/sweeper.cpp b/src/3rdparty/2geom/src/toys/sweeper.cpp
new file mode 100644
index 0000000..7dae586
--- /dev/null
+++ b/src/3rdparty/2geom/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 :
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-15 13:04:31 +0000