src/2geom/shape.cpp

   1 /**
   2  * \brief  Shapes are special paths on which boolops can be performed
   3  *
   4  * Authors:
   5  *      Michael G. Sloan <mgsloan@gmail.com>
   6  *      Nathan Hurst <njh@mail.csse.monash.edu.au>
   7  *      MenTaLguY <mental@rydia.net>
   8  *
   9  * Copyright 2007-2009  Authors
  10  *
  11  * This library is free software; you can redistribute it and/or
  12  * modify it either under the terms of the GNU Lesser General Public
  13  * License version 2.1 as published by the Free Software Foundation
  14  * (the "LGPL") or, at your option, under the terms of the Mozilla
  15  * Public License Version 1.1 (the "MPL"). If you do not alter this
  16  * notice, a recipient may use your version of this file under either
  17  * the MPL or the LGPL.
  18  *
  19  * You should have received a copy of the LGPL along with this library
  20  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  21  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  22  * You should have received a copy of the MPL along with this library
  23  * in the file COPYING-MPL-1.1
  24  *
  25  * The contents of this file are subject to the Mozilla Public License
  26  * Version 1.1 (the "License"); you may not use this file except in
  27  * compliance with the License. You may obtain a copy of the License at
  28  * http://www.mozilla.org/MPL/
  29  *
  30  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  31  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  32  * the specific language governing rights and limitations.
  33  *
  34  */
  35
  36 #include <2geom/shape.h>
  37
  38 #include <2geom/utils.h>
  39 #include <2geom/sweep.h>
  40 #include <2geom/ord.h>
  41
  42 #include <iostream>
  43 #include <algorithm>
  44 #include <cstdlib>
  45
  46 //#define SHAPE_DEBUG  // turns on debug outputting to cout.
  47
  48 namespace Geom {
  49
  50 // A little sugar for appending a list to another
  51 template<typename T>
  52 void append(T &a, T const &b) {
  53     a.insert(a.end(), b.begin(), b.end());
  54 }
  55
  56 //Orders a list of indices according to their containment within eachother.
  57 struct ContainmentOrder {
  58     std::vector<Region> const *rs;
  59     explicit ContainmentOrder(std::vector<Region> const *r) : rs(r) {}
  60     bool operator()(unsigned a, unsigned b) const { return (*rs)[b].contains((*rs)[a]); }
  61 };
  62
  63 //Returns the list of regions containing a particular point.  Useful in tandem with ContainmentOrder
  64 std::vector<unsigned> Shape::containment_list(Point p) const {
  65     std::vector<Rect> pnt;
  66     pnt.push_back(Rect(p, p));
  67     std::vector<std::vector<unsigned> > cull = sweep_bounds(pnt, bounds(*this));
  68     std::vector<unsigned> containers;
  69     if(cull[0].size() == 0) return containers;
  70     for(unsigned i = 0; i < cull[0].size(); i++)
  71         if(content[cull[0][i]].contains(p)) containers.push_back(cull[0][i]);
  72     return containers;
  73 }
  74
  75 /* Used within shape_boolean and related functions, as the name describes, finds the
  76  * first false within the list of lists of booleans.
  77  */
  78 void first_false(std::vector<std::vector<bool> > visited, unsigned &i, unsigned &j) {
  79     for(i = 0, j = 0; i < visited.size(); i++) {
  80         std::vector<bool>::iterator unvisited = std::find(visited[i].begin(), visited[i].end(), false);
  81         if(unvisited != visited[i].end()) {
  82             j = unvisited - visited[i].begin();
  83             break;
  84         }
  85     }
  86 }
  87
  88 // Finds a crossing in a list of them, given the sorting index.
  89 unsigned find_crossing(Crossings const &cr, Crossing x, unsigned i) {
  90     return std::lower_bound(cr.begin(), cr.end(), x, CrossingOrder(i)) - cr.begin();
  91 }
  92
  93 /* This function handles boolean ops on shapes.  The first parameter is a bool
  94  * which determines its behavior in each combination of cases.  For proper
  95  * fill information and noncrossing behavior, the fill data of the regions
  96  * must be correct.  The boolean parameter determines whether the operation
  97  * is a union or a subtraction.  Reversed paths represent inverse regions,
  98  * where everything is included in the fill except for the insides.
  99  *
 100  * Here is a chart of the behavior under various circumstances:
 101  *
 102  * rev = false (union)
 103  *            A
 104  *        F         H
 105  * F  A+B -> F  A-B -> H
 106  *B
 107  * H  B-A -> H  AxB -> H
 108  *
 109  * rev = true (intersect)
 110  *            A
 111  *        F         H
 112  * F  AxB -> F  B-A -> F
 113  *B
 114  * H  A-B -> F  A+B -> H
 115  *
 116  * F/H = Fill outer / Hole outer
 117  * A/B specify operands
 118  * + = union, - = subtraction, x = intersection
 119  * -> read as "produces"
 120  *
 121  * This is the main function of boolops, yet its operation isn't very complicated.
 122  * It traverses the crossings, and uses the crossing direction to decide whether
 123  * the next segment should be taken from A or from B.  The second half of the
 124  * function deals with figuring out what to do with bits that have no intersection.
 125  */
 126 Shape shape_boolean(bool rev, Shape const & a, Shape const & b, CrossingSet const & crs) {
 127     const Regions ac = a.content, bc = b.content;
 128
 129     //Keep track of which crossings we've hit.
 130     std::vector<std::vector<bool> > visited;
 131     for(unsigned i = 0; i < crs.size(); i++)
 132         visited.push_back(std::vector<bool>(crs[i].size(), false));
 133
 134     //Traverse the crossings, creating chunks
 135     Regions chunks;
 136     while(true) {
 137         unsigned i, j;
 138         first_false(visited, i, j);
 139         if(i == visited.size()) break;
 140
 141         Path res;
 142         do {
 143             Crossing cur = crs[i][j];
 144             visited[i][j] = true;
 145
 146             //get indices of the dual:
 147             unsigned io = cur.getOther(i), jo = find_crossing(crs[io], cur, io);
 148             if(jo < visited[io].size()) visited[io][jo] = true;
 149
 150             //main driving logic
 151             if(logical_xor(cur.dir, rev)) {
 152                 if(i >= ac.size()) { i = io; j = jo; }
 153                 j++;
 154                 if(j >= crs[i].size()) j = 0;
 155                 Crossing next = crs[i][j];
 156                 ac[next.a].boundary.appendPortionTo(res, cur.ta, next.ta);
 157             } else {
 158                 if(i < ac.size()) { i = io; j = jo; }
 159                 j++;
 160                 if(j >= crs[i].size()) j = 0;
 161                 Crossing next = crs[i][j];
 162                 bc[next.b - ac.size()].boundary.appendPortionTo(res, cur.tb, next.tb);
 163             }
 164         } while (!visited[i][j]);
 165         if(res.size() > 0) chunks.push_back(Region(res));
 166     }
 167
 168     //If true, then we are on the 'subtraction diagonal'
 169     bool const on_sub = logical_xor(a.fill, b.fill);
 170     //If true, outer paths are filled
 171     bool const res_fill = rev ? (on_sub || (a.fill && b.fill)) : (a.fill && b.fill);
 172
 173     //Handle unintersecting portions
 174     for(unsigned i = 0; i < crs.size(); i++) {
 175         if(crs[i].size() == 0) {
 176             bool env;
 177             bool on_a = i < ac.size();
 178             Region const & r(on_a ? ac[i] : bc[i - ac.size()]);
 179             Shape const & other(on_a ? b : a);
 180
 181             std::vector<unsigned> containers = other.containment_list(r.boundary.initialPoint());
 182             if(containers.empty()) {
 183                 //not included in any container, the environment fill is the opposite of the outer fill
 184                 env = !res_fill;
 185                 if(on_sub && logical_xor(other.fill, res_fill)) env = !env;  //If on the subtractor, invert the environment fill
 186             } else {
 187                 //environment fill is the same as the inner-most container
 188                 std::vector<unsigned>::iterator cit = std::min_element(containers.begin(), containers.end(), ContainmentOrder(&other.content));
 189                 env = other[*cit].isFill();
 190             }
 191             if(!logical_xor(rev, env)) chunks.push_back(r); //When unioning, environment must be hole for inclusion, when intersecting, it must be filled
 192         }
 193     }
 194
 195     return Shape(chunks, res_fill);
 196 }
 197
 198 // Just a convenience wrapper for shape_boolean, which handles the crossings
 199 Shape shape_boolean(bool rev, Shape const & a, Shape const & b) {
 200     CrossingSet crs = crossings_between(a, b);
 201
 202     return shape_boolean(rev, a, b, crs);
 203 }
 204
 205
 206 // Some utility functions for boolop:
 207
 208 std::vector<double> region_sizes(Shape const &a) {
 209     std::vector<double> ret;
 210     for(unsigned i = 0; i < a.size(); i++) {
 211         ret.push_back(double(a[i].size()));
 212     }
 213     return ret;
 214 }
 215
 216 Shape shape_boolean_ra(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) {
 217     return shape_boolean(rev, a.inverse(), b, reverse_ta(crs, a.size(), region_sizes(a)));
 218 }
 219
 220 Shape shape_boolean_rb(bool rev, Shape const &a, Shape const &b, CrossingSet const &crs) {
 221     return shape_boolean(rev, a, b.inverse(), reverse_tb(crs, a.size(), region_sizes(b)));
 222 }
 223
 224 /* This is a function based on shape_boolean which allows boolean operations
 225  * to be specified as a logic table.  This logic table is 4 bit-flags, which
 226  * correspond to the elements of the 'truth table' for a particular operation.
 227  * These flags are defined with the enums starting with BOOLOP_ .
 228  *
 229  * NOTE: currently doesn't work, as the CrossingSet reversal functions crash
 230  */
 231 Shape boolop(Shape const &a, Shape const &b, unsigned flags, CrossingSet const &crs) {
 232     THROW_NOTIMPLEMENTED();
 233     flags &= 15;
 234     if(flags <= BOOLOP_UNION) {
 235         switch(flags) {
 236             case BOOLOP_INTERSECT:    return shape_boolean(true, a, b, crs);
 237             case BOOLOP_SUBTRACT_A_B: return shape_boolean_rb(true, a, b, crs);
 238             case BOOLOP_IDENTITY_A:   return a;
 239             case BOOLOP_SUBTRACT_B_A: return shape_boolean_ra(true, a, b, crs);
 240             case BOOLOP_IDENTITY_B:   return b;
 241             case BOOLOP_EXCLUSION: {
 242                 Shape res = shape_boolean_rb(true, a, b, crs);
 243                 append(res.content, shape_boolean_ra(true, a, b, crs).content);
 244                 return res;
 245             }
 246             case BOOLOP_UNION:        return shape_boolean(false, a, b);
 247         }
 248     } else {
 249         flags = ~flags & 15;
 250         switch(flags - BOOLOP_NEITHER) {
 251             case BOOLOP_SUBTRACT_A_B: return shape_boolean_ra(false, a, b, crs);
 252             case BOOLOP_SUBTRACT_B_A: return shape_boolean_rb(false, a, b, crs);
 253             case BOOLOP_EXCLUSION: {
 254                 Shape res = shape_boolean_ra(false, a, b, CrossingSet(crs));
 255                 append(res.content, shape_boolean_rb(false, a, b, CrossingSet(crs)).content);
 256                 return res;
 257             }
 258         }
 259         return boolop(a, b, flags, crs).inverse();
 260     }
 261     return Shape();
 262 }
 263
 264 /* This version of the boolop function doesn't require a set of crossings, as
 265  * it computes them for you.  This is more efficient in some cases, as the
 266  * shape can be inverted before finding crossings.  In the special case of
 267  * exclusion it uses the other version of boolop.
 268  */
 269 Shape boolop(Shape const &a, Shape const &b, unsigned flags) {
 270     flags &= 15;
 271     if(flags <= BOOLOP_UNION) {
 272         switch(flags) {
 273             case BOOLOP_INTERSECT:    return shape_boolean(true, a, b);
 274             case BOOLOP_SUBTRACT_A_B: return shape_boolean(true, a, b.inverse());
 275             case BOOLOP_IDENTITY_A:   return a;
 276             case BOOLOP_SUBTRACT_B_A: return shape_boolean(true, b, a.inverse());
 277             case BOOLOP_IDENTITY_B:   return b;
 278             case BOOLOP_EXCLUSION: {
 279                 Shape res = shape_boolean(true, a, b.inverse());
 280                 append(res.content, shape_boolean(true, b, a.inverse()).content);
 281                 return res;
 282             } //return boolop(a, b, flags,  crossings_between(a, b));
 283             case BOOLOP_UNION:        return shape_boolean(false, a, b);
 284         }
 285     } else {
 286         flags = ~flags & 15;
 287         switch(flags) {
 288             case BOOLOP_SUBTRACT_A_B: return shape_boolean(false, b, a.inverse());
 289             case BOOLOP_SUBTRACT_B_A: return shape_boolean(false, a, b.inverse());
 290             case BOOLOP_EXCLUSION: {
 291                 Shape res = shape_boolean(false, a, b.inverse());
 292                 append(res.content, shape_boolean(false, b, a.inverse()).content);
 293                 return res;
 294             } //return boolop(a, b, flags, crossings_between(a, b));
 295         }
 296         return boolop(a, b, flags).inverse();
 297     }
 298     return Shape();
 299 }
 300
 301 int paths_winding(std::vector<Path> const &ps, Point p) {
 302     int ret = 0;
 303     for(unsigned i = 0; i < ps.size(); i++)
 304         ret += winding(ps[i], p);
 305     return ret;
 306 }
 307
 308 void add_to_shape(Shape &s, Path const &p, bool fill) {
 309     if(fill)
 310         s.content.push_back(Region(p).asFill());
 311     else
 312         s.content.push_back(Region(p).asHole());
 313 }
 314
 315 int inner_winding(Path const & p, std::vector<Path> const &ps) {
 316     Point pnt = p.initialPoint();
 317     return paths_winding(ps, pnt) - winding(p, pnt) + 1;
 318 }
 319
 320 double fudgerize(double d, bool rev) {
 321     double ret = rev ? d - 0.01 : d + 0.01;
 322     if(ret < 0) ret = 0;
 323     return ret;
 324 }
 325
 326 unsigned pick_coincident(unsigned ix, unsigned jx, bool &rev, std::vector<Path> const &ps, CrossingSet const &crs) {
 327     unsigned ex_jx = jx;
 328     unsigned oix = crs[ix][jx].getOther(ix);
 329     double otime = crs[ix][jx].getTime(oix);
 330     Point cross_point = ps[oix].pointAt(otime),
 331           along = ps[oix].pointAt(fudgerize(otime, rev)) - cross_point,
 332           prev = -along;
 333     bool ex_dir = rev;
 334     for(unsigned k = jx; k < crs[ix].size(); k++) {
 335         unsigned koix = crs[ix][k].getOther(ix);
 336         if(koix == oix) {
 337             if(!are_near(otime, crs[ix][k].getTime(oix))) break;
 338             for(unsigned dir = 0; dir < 2; dir++) {
 339                 Point val = ps[ix].pointAt(fudgerize(crs[ix][k].getTime(ix), dir)) - cross_point;
 340                 Cmp to_prev = cmp(cross(val, prev), 0);
 341                 Cmp from_along = cmp(cross(along, val), 0);
 342                 Cmp c = cmp(from_along, to_prev);
 343                 if(c == EQUAL_TO && from_along == LESS_THAN) {
 344                     ex_jx = k;
 345                     prev = val;
 346                     ex_dir = dir;
 347                 }
 348             }
 349         }
 350     }
 351     rev = ex_dir;
 352     return ex_jx;
 353 }
 354
 355 unsigned crossing_along(double t, unsigned ix, unsigned jx, bool dir, Crossings const & crs) {
 356     Crossing cur = Crossing(t, t, ix, ix, false);
 357     if(jx < crs.size()) {
 358         double ct = crs[jx].getTime(ix);
 359         if(t == ct) {
 360             cur = crs[jx];
 361             if(cur.a == cur.b) {
 362                 if(jx+1 <= crs.size() && crs[jx+1].getOther(ix) == ix) return jx+1;
 363                 if(jx > 0 && crs[jx-1].getOther(ix) == ix) return jx-1;
 364             }
 365         }
 366     }
 367     if(!dir) {
 368         jx = std::upper_bound(crs.begin(), crs.end(), cur, CrossingOrder(ix)) - crs.begin();
 369     } else {
 370         jx = std::lower_bound(crs.begin(), crs.end(), cur, CrossingOrder(ix)) - crs.begin();
 371         if(jx == 0) jx = crs.size() - 1; else jx--;
 372         jx = std::lower_bound(crs.begin(), crs.end(), crs[jx], CrossingOrder(ix)) - crs.begin();
 373     }
 374     if(jx >= crs.size()) jx = 0;
 375     return jx;
 376 }
 377
 378 void crossing_dual(unsigned &i, unsigned &j, CrossingSet const & crs) {
 379     Crossing cur = crs[i][j];
 380     i = cur.getOther(i);
 381 #ifdef SHAPE_DEBUG
 382     std::cout << i << "\n";
 383 #endif
 384     if(crs[i].empty())
 385         j = 0;
 386     else
 387         j = std::lower_bound(crs[i].begin(), crs[i].end(), cur, CrossingOrder(i)) - crs[i].begin();
 388 }
 389
 390 //locate a crossing on the outside, by casting a ray through the middle of the bbox
 391 void outer_crossing(unsigned &ix, unsigned &jx, bool & dir, std::vector<Path> const & ps, CrossingSet const & crs) {
 392     Rect bounds = *(ps[ix].boundsFast());
 393     double ry = bounds[Y].middle();
 394     double max_val = bounds.left(), max_t = 0;
 395     ix = ps.size();
 396     for(unsigned i = 0; i < ps.size(); i++) {
 397         if(!crs[i].empty()) {
 398             std::vector<double> rts = ps[i].roots(ry, Y);
 399             for(unsigned j = 0; j < rts.size(); j++) {
 400                 double val = ps[i].valueAt(rts[j], X);
 401                 if(val > max_val) {
 402                     ix = i;
 403                     max_val = val;
 404                     max_t = rts[j];
 405                 }
 406             }
 407         }
 408     }
 409     if(ix != ps.size()) {
 410         dir = ps[ix].valueAt(max_t + 0.01, Y) >
 411               ps[ix].valueAt(max_t - 0.01, Y);
 412         jx = crossing_along(max_t, ix, jx, dir, crs[ix]);
 413     }
 414 }
 415
 416 std::vector<Path> inner_sanitize(std::vector<Path> const & ps) {
 417     CrossingSet crs(crossings_among(ps));
 418
 419     Regions chunks;
 420
 421     std::vector<bool> used_path(ps.size(), false);
 422     std::vector<std::vector<bool> > visited;
 423     for(unsigned i = 0; i < crs.size(); i++)
 424         visited.push_back(std::vector<bool>(crs[i].size(), false));
 425
 426     std::vector<Path> result_paths;
 427
 428     while(true) {
 429         unsigned ix = 0, jx = 0;
 430         bool dir = false;
 431
 432         //find an outer crossing by trying various paths and checking if the crossings are used
 433         for(; ix < crs.size(); ix++) {
 434             //TODO: optimize so it doesn't unecessarily check stuff
 435             bool cont = true;
 436             for(unsigned j = 0; j < crs[ix].size(); j++) {
 437                 if(!visited[ix][j]) { cont = false; break; }
 438             }
 439             if(cont) continue;
 440             unsigned rix = ix, rjx = jx;
 441             outer_crossing(rix, rjx, dir, ps, crs);
 442             if(rix >= crs.size() || visited[rix][rjx]) continue;
 443             ix = rix; jx = rjx;
 444             break;
 445         }
 446         if(ix == crs.size()) break;
 447         crossing_dual(ix, jx, crs);
 448
 449         dir = !dir;
 450
 451         Path res;
 452         do {
 453             visited[ix][jx] = true;
 454             //unsigned nix = ix, njx = jx;
 455             //crossing_dual(nix, njx, crs);
 456             //visited[nix][njx] = true;
 457             unsigned fix = ix, fjx = jx;
 458
 459             bool new_dir = dir;
 460
 461             jx = crossing_along(crs[ix][jx].getTime(ix), ix, jx, dir, crs[ix]);
 462             if(crs[ix][jx].a != crs[ix][jx].b) crossing_dual(ix, jx, crs); else new_dir = !new_dir;
 463             jx = pick_coincident(ix, jx, new_dir, ps, crs);
 464
 465             //unsigned nix = ix, njx = jx;
 466             //crossing_dual(nix, njx, crs);
 467
 468             Crossing from = crs[fix][fjx],
 469                      to = crs[ix][jx];
 470             if(dir) {
 471                 // backwards
 472 #ifdef SHAPE_DEBUG
 473                 std::cout << "r" << ix << "[" << from.getTime(ix)  << ", " << to.getTime(ix) << "]\n";
 474 #endif
 475                 Path p = ps[ix].portion(from.getTime(ix), to.getTime(ix)).reverse();
 476                 for(unsigned i = 0; i < p.size(); i++)
 477                     res.append(p[i], Path::STITCH_DISCONTINUOUS);
 478             } else {
 479                 // forwards
 480 #ifdef SHAPE_DEBUG
 481                 std::cout << "f" << ix << "[" << from.getTime(ix) << ", " << to.getTime(ix) << "]\n";
 482 #endif
 483                 ps[ix].appendPortionTo(res, from.getTime(ix), to.getTime(ix));
 484             }
 485             dir = new_dir;
 486         } while(!visited[ix][jx]);
 487 #ifdef SHAPE_DEBUG
 488         std::cout << "added " << res.size() << "\n";
 489 #endif
 490         result_paths.push_back(res);
 491     }
 492     for(unsigned i = 0; i < crs.size(); i++) {
 493         if(crs[i].empty() && !used_path[i])
 494             result_paths.push_back(ps[i]);
 495     }
 496     return result_paths;
 497 }
 498
 499 Shape sanitize(std::vector<Path> const & ps) {
 500     std::vector<Path> res;
 501     for(unsigned i = 0; i < ps.size(); i++) {
 502         append(res, inner_sanitize(std::vector<Path>(1, ps[i])));
 503     }
 504     return stopgap_cleaner(res);
 505 }
 506
 507 /*  WIP sanitizer:
 508 unsigned pick_coincident(unsigned ix, unsigned jx, bool pref, bool &rev, std::vector<Path> const &ps, CrossingSet const &crs) {
 509     unsigned ex_jx = jx;
 510     unsigned oix = crs[ix][jx].getOther(ix);
 511     double otime = crs[ix][jx].getTime(oix);
 512     Point cross_point = ps[oix].pointAt(otime),
 513           along = ps[oix].pointAt(otime + (rev ? -0.01 : 0.01)) - cross_point,
 514           prev = -along;
 515     bool ex_dir = rev;
 516     for(unsigned k = jx; k < crs[ix].size(); k++) {
 517         unsigned koix = crs[ix][k].getOther(ix);
 518         if(koix == oix) {
 519             if(!are_near(otime, crs[ix][k].getTime(oix))) break;
 520             for(unsigned dir = 0; dir < 2; dir++) {
 521                 Point val = ps[ix].pointAt(crs[ix][k].getTime(ix) + (dir ? -0.01 : 0.01)) - cross_point;
 522                 Cmp to_prev = cmp(cross(val, prev), 0);
 523                 Cmp from_along = cmp(cross(along, val), 0);
 524                 Cmp c = cmp(from_along, to_prev);
 525                 if(c == EQUAL_TO && (from_along == LESS_THAN) == pref) {
 526                     ex_jx = k;
 527                     prev = val;
 528                     ex_dir = dir;
 529                 }
 530             }
 531         }
 532     }
 533     rev = ex_dir;
 534     return ex_jx;
 535 }
 536
 537 unsigned corner_index(unsigned &i) {
 538     div_t div_res = div(i, 4);
 539     i = div_res.quot;
 540     return div_res.rem;
 541 }
 542
 543 bool corner_direction(unsigned ix, unsigned jc, unsigned corner, CrossingSet const &crs) {
 544     if(crs[ix][jc].a == ix) return corner > 1; else return corner %2 == 1;
 545 }
 546
 547 Shape sanitize(std::vector<Path> const & ps) {
 548     CrossingSet crs = crossings_among(ps);
 549
 550     //Keep track of which CORNERS we've hit.
 551     // FF FR RF RR, first is A dir, second B dir
 552     std::vector<std::vector<bool> > visited;
 553     for(unsigned i = 0; i < crs.size(); i++)
 554         visited.push_back(std::vector<bool>(crs[i].size()*4, false));
 555
 556     Regions chunks;
 557     while(true) {
 558         unsigned i, j;
 559         first_false(visited, i, j);
 560         unsigned corner = corner_index(j);
 561
 562         if(i == visited.size()) break;
 563
 564         bool dir = corner_direction(i, j, corner, crs);
 565
 566         //Figure out whether we hug the path cw or ccw, based on the orientation of the initial corner:
 567         unsigned oix = crs[i][j].getOther(i);
 568         double otime = crs[i][j].getTime(oix);
 569         bool odir = (oix == crs[i][j].a) ? corner > 1 : corner % 2 == 1;
 570         Point cross_point = ps[oix].pointAt(otime),
 571               along = ps[oix].pointAt(otime + (odir ? -0.01 : 0.01)) - cross_point,
 572                 val = ps[i].pointAt(crs[i][j].getTime(i) + (dir ? -0.01 : 0.01)) - cross_point;
 573
 574         Cmp from_along = cmp(cross(along, val), 0);
 575         bool cw = from_along == LESS_THAN;
 576         std::cout << "cw = " << cw << "\n";
 577         Path res;
 578         do {
 579             Crossing cur = crs[i][j];
 580             visited[i][j*4+corner] = true;
 581
 582             unsigned fix = i, fjx = j;
 583             crossing_dual(i, j, crs);
 584             visited[i][j*4+corner] = true;
 585             i = fix; j = fjx;
 586
 587             j = crossing_along(crs[i][j].getTime(i), i, j, dir, crs[i]);
 588
 589             crossing_dual(i, j, crs);
 590
 591             bool new_dir = dir;
 592             pick_coincident(i, j, cw, new_dir, ps, crs);
 593
 594             Crossing from = crs[fix][fjx],
 595                      to = crs[i][j];
 596             if(dir) {
 597                 // backwards
 598                 std::cout << "r" << i << "[" << to.getTime(i)  << ", " << from.getTime(i) << "]\n";
 599                 Path p = ps[i].portion(to.getTime(i) + 0.001, from.getTime(i)).reverse();
 600                 for(unsigned k = 0; k < p.size(); k++)
 601                     res.append(p[k]);
 602             } else {
 603                 // forwards
 604                 std::cout << "f" << i << "[" << from.getTime(i) << ", " << to.getTime(i) << "]\n";
 605                 ps[i].appendPortionTo(res, from.getTime(i) + 0.001, to.getTime(i));
 606             }
 607             if(i == to.a)
 608                 corner = (new_dir ? 2 : 0) + (dir ? 1 : 0);
 609             else
 610                 corner = (new_dir ? 1 : 0) + (dir ? 2 : 0);
 611             dir = new_dir;
 612         } while(!visited[i][j*4+corner]);
 613         chunks.push_back(Region(res));
 614 //        if(use) {
 615 //            chunks.push_back(Region(res, true));
 616 //        }
 617     }
 618     return Shape(chunks);
 619 //    return ret;
 620 } */
 621
 622 /* This transforms a shape by a matrix.  In the case that the matrix flips
 623  * the shape, it reverses the paths in order to preserve the fill.
 624  */
 625 Shape Shape::operator*(Matrix const &m) const {
 626     Shape ret;
 627     for(unsigned i = 0; i < size(); i++)
 628         ret.content.push_back(content[i] * m);
 629     ret.fill = fill;
 630     return ret;
 631 }
 632
 633 // Inverse is a boolean not, and simply reverses all the paths & fill flags
 634 Shape Shape::inverse() const {
 635     Shape ret;
 636     for(unsigned i = 0; i < size(); i++)
 637         ret.content.push_back(content[i].inverse());
 638     ret.fill = !fill;
 639     return ret;
 640 }
 641
 642 bool Shape::contains(Point const &p) const {
 643     std::vector<unsigned> containers = containment_list(p);
 644     if(containers.empty()) return !isFill();
 645     unsigned ix = *min_element(containers.begin(), containers.end(), ContainmentOrder(&content));
 646     return content[ix].isFill();
 647 }
 648
 649 Shape stopgap_cleaner(std::vector<Path> const &ps) {
 650     if(ps.empty()) return Shape(false);
 651     Shape ret;
 652     for(unsigned i = 0; i < ps.size(); i++)
 653         add_to_shape(ret, ps[i], inner_winding(ps[i], ps) % 2 != 0);
 654     return ret;
 655 }
 656
 657 bool Shape::inside_invariants() const {  //semi-slow & easy to violate
 658     for(unsigned i = 0; i < size(); i++)
 659         if( logical_xor(content[i].isFill(), contains(content[i].boundary.initialPoint())) ) return false;
 660     return true;
 661 }
 662 bool Shape::region_invariants() const { //semi-slow
 663     for(unsigned i = 0; i < size(); i++)
 664         if(!content[i].invariants()) return false;
 665     return true;
 666 }
 667 bool Shape::cross_invariants() const { //slow
 668     CrossingSet crs; // = crossings_among(paths_from_regions(content));
 669     for(unsigned i = 0; i < crs.size(); i++)
 670         if(!crs[i].empty()) return false;
 671     return true;
 672 }
 673
 674 bool Shape::invariants() const {
 675     return inside_invariants() && region_invariants() && cross_invariants();
 676 }
 677
 678 }
 679
 680 /*
 681   Local Variables:
 682   mode:c++
 683   c-file-style:"stroustrup"
 684   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 685   indent-tabs-mode:nil
 686   fill-column:99
 687   End:
 688 */
 689 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :