vmath2d/ybdecomposer.d

   1 /*
   2  * coded by Ketmar // Invisible Vector <ketmar@ketmar.no-ip.org>
   3  * Understanding is not required. Only obedience.
   4  *
   5  * This program is free software: you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation, version 3 of the License ONLY.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  16  */
  17 // convex decomposition algorithm originally created by Mark Bayazit (http://mnbayazit.com/)
  18 // for more information about this algorithm, see http://mnbayazit.com/406/bayazit
  19 // modified by Yogesh (http://yogeshkulkarni.com)
  20 // D port and further modifications by Ketmar // Invisible Vector
  21 //
  22 // decompose simple polygon (i.e. polygon without holes and self-intersections) into set of convex polygons
  23 module iv.vmath2d.ybdecomposer;
  24
  25 import iv.alice;
  26 import iv.vmath2d.math2d;
  27 import iv.vmath2d.vxstore;
  28
  29
  30 // ////////////////////////////////////////////////////////////////////////// //
  31 public struct YBDecomposer(VT) if (IsVectorDim!(VT, 2)) {
  32 private:
  33   // precondition: ccw
  34   static bool reflex() (in auto ref VT prev, in auto ref VT on, in auto ref VT next) {
  35     pragma(inline, true);
  36     // YOGESH: Added following condition of collinearity
  37     if (Math2D.isCollinear(prev, on, next)) return false;
  38     return right(prev, on, next);
  39   }
  40
  41   // if i is on the left of i-1 and i+1 line in the polygon vertices; checks area -ve
  42   static bool left() (in auto ref VT a, in auto ref VT b, in auto ref VT c) { pragma(inline, true); return (Math2D.area(a, b, c) > 0); }
  43
  44   // if i is on the left or ON of i-1 and i+1 line in the polygon vertices; checks area -ve and 0
  45   static bool leftOn() (in auto ref VT a, in auto ref VT b, in auto ref VT c) { pragma(inline, true); return (Math2D.area(a, b, c) >= 0); }
  46
  47   // if i is on the right of i-1 and i+1 line in the polygon vertices; checks area +ve
  48   static bool right() (in auto ref VT a, in auto ref VT b, in auto ref VT c) { pragma(inline, true); return (Math2D.area(a, b, c) < 0); }
  49
  50   // if i is on the right or ON of i-1 and i+1 line in the polygon vertices; checks area +ve and 0
  51   static bool rightOn() (in auto ref VT a, in auto ref VT b, in auto ref VT c) { pragma(inline, true); return (Math2D.area(a, b, c) <= 0); }
  52
  53 public:
  54   // decompose the polygon into several smaller non-concave polygon
  55   // if the polygon is already convex, it will return the original polygon
  56   // returned polygons are CCW
  57   static void convexPartition(VS) (ref VS[] accum, ref VS poly) if (IsGoodVertexStorage!(VS, VT)) {
  58     // can you see j from i in the polygon vertices without any obstructions?
  59     bool canSee (int i, int j) {
  60       VT prev = poly[i-1];
  61       VT on = poly[i];
  62       VT next = poly[i+1];
  63
  64       if (reflex(prev, on, next)) {
  65         if (leftOn(poly[i], poly[i-1], poly[j]) && rightOn(poly[i], poly[i+1], poly[j])) return false;
  66       } else {
  67         if (rightOn(poly[i], poly[i+1], poly[j]) || leftOn(poly[i], poly[i-1], poly[j])) return false;
  68       }
  69
  70       VT prevj = poly[j-1];
  71       VT onj = poly[j];
  72       VT nextj = poly[j+1];
  73
  74       if (reflex(prevj, onj, nextj)) {
  75         if (leftOn(poly[j], poly[j-1], poly[i]) && rightOn(poly[j], poly[j+1], poly[i])) return false;
  76       } else {
  77         if (rightOn(poly[j], poly[j+1], poly[i]) || leftOn(poly[j], poly[j-1], poly[i])) return false;
  78       }
  79
  80       foreach (immutable k; 0..poly.length) {
  81         //if ((k + 1) % vertices.Count == i || k == i || (k + 1) % vertices.Count == j || k == j) continue; // ignore incident edges
  82
  83         // YOGESH: changed from Line-Line intersection to Segment-Segment Intersection
  84         VT p1 = poly[i];
  85         VT p2 = poly[j];
  86         VT q1 = poly[k];
  87         VT q2 = poly[k+1];
  88
  89         // ignore incident edges
  90         if (p1 == q1 || p1 == q2 || p2 == q1 || p2 == q2) continue;
  91
  92         auto intPoint = Math2D.segIntersect(p1, p2, q1, q2);
  93         //Debug.Print("Line from point {0} to {1} is tested against [{2},{3} to see if they intersect]", i, j, k, k+1);
  94         if (intPoint.valid) {
  95           // intPoint is not any of the j line then false, else continue. Intersection has to be interior to qualify s 'false' from here
  96           //if (intPoint != poly.verts[k] || intPoint != poly.verts[k+1]) return false;
  97           if (intPoint.equals(poly[k]) || intPoint.equals(poly[k+1])) return false;
  98         }
  99       }
 100
 101       return true;
 102     }
 103
 104     void copyPolyPartTo (ref VS dest, int i, int j) {
 105       while (j < i) j += poly.length;
 106       for (; i <= j; ++i) dest ~= poly[i];
 107     }
 108
 109     // main
 110     poly.collinearSimplify();
 111     if (poly.length < 3) return;
 112
 113     // we force it to CCW as it is a precondition in this algorithm
 114     poly.forceCCW();
 115
 116     //VT.Float lowerDist = 0.0, upperDist = 0.0;
 117     //int lowerIndex = 0, upperIndex = 0;
 118     VT lowerInt, upperInt; // intersection points
 119
 120     // go thru all Verices until we find a reflex vertex i
 121     // extend the edges incident at i until they hit an edge
 122     // find BEST vertex within the range, that the partitioning chord
 123
 124     // a polygon can be broken into convex regions by eliminating all reflex vertices
 125     // eliminating two reflex vertices with one diagonal is better than eliminating just one
 126     // a reflex vertex can only be removed if the diagonal connecting to it is within the range given by extending its neighbouring edges;
 127     // otherwise, its angle is only reduced
 128     foreach (immutable i; 0..poly.length) {
 129       if (reflex(poly[i-1], poly[i], poly[i+1])) {
 130         VT.Float lowerDist = VT.Float.infinity;
 131         VT.Float upperDist = VT.Float.infinity;
 132         int lowerIndex = 0, upperIndex = 0;
 133         for (int j = 0; j < poly.length; ++j) {
 134           // YOGESH: if any of j line's endpoints matches with reflex i, skip
 135           if (poly.sameIndex(i, j) || poly.sameIndex(i, j-1) || poly.sameIndex(i, j+1)) continue; // no self and prev and next, for testing
 136
 137           // testing incoming edge:
 138           // if line coming into i vertex (i-1 to i) has j vertex of the test-line on left
 139           // AND have j-i on right, then they will be intersecting
 140
 141           VT iPrev = poly[i-1];
 142           VT iSelf = poly[i];
 143           VT jSelf = poly[j];
 144           VT jPrev = poly[j-1];
 145
 146           bool leftOK = left(iPrev, iSelf, jSelf);
 147           bool rightOK = right(iPrev, iSelf, jPrev);
 148
 149           bool leftOnOK = Math2D.isCollinear(iPrev, iSelf, jSelf); // YOGESH: cached into variables for better debugging
 150           bool rightOnOK = Math2D.isCollinear(iPrev, iSelf, jPrev); // YOGESH: cached into variables for better debugging
 151
 152           if (leftOnOK || rightOnOK) {
 153             // YOGESH: Checked "ON" condition as well, collinearity
 154             // lines are colinear, they can not be overlapping as polygon is simple
 155             // find closest point which is not internal to incoming line i , i -1
 156             VT.Float d = iSelf.distanceSquared(jSelf);
 157
 158             // this lower* is the point got from incoming edge into the i vertex,
 159             // lowerInt incoming edge intersection point
 160             // lowerIndex incoming edge intersection edge
 161             if (d < lowerDist) {
 162               // keep only the closest intersection
 163               lowerDist = d;
 164               lowerInt = jSelf;
 165               lowerIndex = j-1;
 166             }
 167
 168             d = iSelf.distanceSquared(jPrev);
 169
 170             // this lower* is the point got from incoming edge into the i vertex,
 171             // lowerInt incoming edge intersection point
 172             // lowerIndex incoming edge intersection edge
 173             if (d < lowerDist) {
 174               // keep only the closest intersection
 175               lowerDist = d;
 176               lowerInt = jPrev;
 177               lowerIndex = j;
 178             }
 179           } else if (leftOK && rightOK) {
 180             // YOGESH: Intersection in-between. Bayazit had ON condition in built here, which I have taken care above.
 181             // find the point of intersection
 182             VT p = Math2D.lineIntersect(poly[i-1], poly[i], poly[j], poly[j-1]);
 183             // make sure it's inside the poly,
 184             if (right(poly[i+1], poly[i], p)) {
 185               VT.Float d = poly[i].distanceSquared(p);
 186               // this lower* is the point got from incoming edge into the i vertex,
 187               // lowerInt incoming edge intersection point
 188               // lowerIndex incoming edge intersection edge
 189               if (d < lowerDist) {
 190                 // keep only the closest intersection
 191                 lowerDist = d;
 192                 lowerInt = p;
 193                 lowerIndex = j;
 194               }
 195             }
 196           }
 197
 198           // testing outgoing edge:
 199           // if line outgoing from i vertex (i to i+1) has j vertex of the test-line on right
 200           // AND has j+1 on left, they they will be intersecting
 201
 202           VT iNext = poly[i+1];
 203           VT jNext = poly[j+1];
 204
 205           bool leftOKn = left(iNext, iSelf, jNext);
 206           bool rightOKn = right(iNext, iSelf, jSelf);
 207
 208           bool leftOnOKn = Math2D.isCollinear(iNext, iSelf, jNext); // YOGESH: cached into variables for better debugging
 209           bool rightOnOKn = Math2D.isCollinear(iNext, iSelf, jSelf);
 210
 211           if (leftOnOKn || rightOnOKn) {
 212             // YOGESH: Checked "ON" condition as well, collinearity
 213             // lines are colinear, they can not be overlapping as polygon is simple
 214             // find closest point which is not internal to incoming line i , i -1
 215             VT.Float d = iSelf.distanceSquared(jNext);
 216
 217             // this upper* is the point got from outgoing edge into the i vertex,
 218             // upperInt outgoing edge intersection point
 219             // upperIndex outgoing edge intersection edge
 220             if (d < upperDist) {
 221               // keep only the closest intersection
 222               upperDist = d;
 223               upperInt = jNext;
 224               upperIndex = j+1;
 225             }
 226
 227             d = poly[i].distanceSquared(poly[j]);
 228
 229             // this upper* is the point got from outgoing edge into the i vertex,
 230             // upperInt outgoing edge intersection point
 231             // upperIndex outgoing edge intersection edge
 232             if (d < upperDist) {
 233               // keep only the closest intersection
 234               upperDist = d;
 235               upperInt = jSelf;
 236               upperIndex = j;
 237             }
 238           } else if (leftOKn && rightOKn) {
 239             // YOGESH: Intersection in-between. Bayazit had ON condition in built here, which I have taken care above.
 240             VT p = Math2D.lineIntersect(poly[i+1], poly[i], poly[j], poly[j+1]);
 241             if (left(poly[i-1], poly[i], p)) {
 242               VT.Float d = poly[i].distanceSquared(p);
 243               // this upper* is the point got from outgoing edge from the i vertex,
 244               // upperInt outgoing edge intersection point
 245               // upperIndex outgoing edge intersection edge
 246               if (d < upperDist) {
 247                 upperDist = d;
 248                 upperIndex = j;
 249                 upperInt = p;
 250               }
 251             }
 252           }
 253         }
 254
 255         VS lowerPoly, upperPoly;
 256         static if (is(VS == class)) {
 257           lowerPoly = new VS();
 258           upperPoly = new VS();
 259         }
 260
 261         // YOGESH: If no vertices in the range, lets not choose midpoint but closet point of that segment
 262         // if there are no vertices to connect to, choose a point in the middle
 263         if (lowerIndex == (upperIndex+1)%poly.length) {
 264           VT sp = ((lowerInt+upperInt)/cast(VT.Float)2);
 265           copyPolyPartTo(lowerPoly, i, upperIndex);
 266           lowerPoly ~= sp;
 267           copyPolyPartTo(upperPoly, lowerIndex, i);
 268           upperPoly ~= sp;
 269         } else {
 270           // find vertex to connect to
 271           VT.Float highestScore = 0, bestIndex = lowerIndex;
 272           while (upperIndex < lowerIndex) upperIndex += poly.length;
 273
 274           // go through all the vertices between the range of lower and upper
 275           for (int j = lowerIndex; j <= upperIndex; ++j) {
 276             if (canSee(i, j)) {
 277               VT.Float score = cast(VT.Float)1/(poly[i].distanceSquared(poly[j])+1);
 278
 279               // if another vertex is reflex, choosing it has highest score
 280               VT prevj = poly[j-1];
 281               VT onj = poly[j];
 282               VT nextj = poly[j+1];
 283
 284               if (reflex(prevj, onj, nextj)) {
 285                 if (rightOn(poly[j-1], poly[j], poly[i]) && leftOn(poly[j+1], poly[j], poly[i])) {
 286                   score += 3;
 287                 } else {
 288                   score += 2;
 289                 }
 290               } else {
 291                 score += 1;
 292               }
 293               if (score > highestScore) {
 294                 bestIndex = j;
 295                 highestScore = score;
 296               }
 297             }
 298           }
 299
 300           // YOGESH : Pending: if there are 2 vertices as 'bestIndex', its better to disregard both and put midpoint (M case)
 301           copyPolyPartTo(lowerPoly, i, cast(int)bestIndex);
 302           copyPolyPartTo(upperPoly, cast(int)bestIndex, i);
 303         }
 304
 305         // solve smallest poly first (SAW in Bayazit's C++ code)
 306         if (lowerPoly.length < upperPoly.length) {
 307           convexPartition(accum, lowerPoly);
 308           convexPartition(accum, upperPoly);
 309         } else {
 310           convexPartition(accum, upperPoly);
 311           convexPartition(accum, lowerPoly);
 312         }
 313         return;
 314       }
 315     }
 316
 317     // polygon is already convex
 318     accum ~= poly;
 319   }
 320 }