@@ -1265,6 +1265,58 @@ namespace sgl {
1265
1265
1266
1266
namespace {
1267
1267
1268
+
1269
+ static bool feq (double x, double y) {
1270
+ return !((x < y) | (x > y));
1271
+ }
1272
+
1273
+ static bool eq_zero (double x) {
1274
+ return feq (x, 0 );
1275
+ }
1276
+
1277
+ static bool collinear (
1278
+ double x1, double y1, // point 1
1279
+ double x2, double y2, // point 2
1280
+ double x3, double y3 // point 3
1281
+ ) {
1282
+ bool x1x2 = feq (x1, x2);
1283
+ bool x1x3 = feq (x1, x3);
1284
+ bool x2x3 = feq (x2, x3);
1285
+ bool y1y2 = feq (y1, y2);
1286
+ bool y1y3 = feq (y1, y3);
1287
+ bool y2y3 = feq (y2, y3);
1288
+ if (x1x2) {
1289
+ return x1x3;
1290
+ }
1291
+ if (y1y2) {
1292
+ return y1y3;
1293
+ }
1294
+ if ((x1x2 & y1y2) | (x1x3 & y1y3) | (x2x3 & y2y3)) {
1295
+ return true ;
1296
+ }
1297
+ double cx1 = x3 - x1;
1298
+ double cy1 = y3 - y1;
1299
+ double cx2 = x2 - x1;
1300
+ double cy2 = y2 - y1;
1301
+ double s1 = cx1 * cy2;
1302
+ double s2 = cy1 * cx2;
1303
+ // Check if precision was lost.
1304
+ double s3 = (s1 / cy2) - cx1;
1305
+ double s4 = (s2 / cx2) - cy1;
1306
+ if (s3 < 0 ) {
1307
+ s1 = nexttoward (s1, -INFINITY);
1308
+ } else if (s3 > 0 ) {
1309
+ s1 = nexttoward (s1, +INFINITY);
1310
+ }
1311
+ if (s4 < 0 ) {
1312
+ s2 = nexttoward (s2, -INFINITY);
1313
+ } else if (s4 > 0 ) {
1314
+ s2 = nexttoward (s2, +INFINITY);
1315
+ }
1316
+ return eq_zero (s1 - s2);
1317
+ }
1318
+
1319
+
1268
1320
// TODO: Make robust
1269
1321
// Returns the orientation of the triplet (p, q, r)
1270
1322
// 0 if collinear, >0 if clockwise, <0 if counter-clockwise
@@ -1280,6 +1332,93 @@ enum class raycast_result {
1280
1332
CROSS,
1281
1333
BOUNDARY
1282
1334
};
1335
+ static bool pteq (const vertex_xy &a, const vertex_xy &b) {
1336
+ return feq (a.x , b.x ) && feq (a.y , b.y );
1337
+ }
1338
+
1339
+ static raycast_result raycast_robust (const vertex_xy &prev, const vertex_xy &next, const vertex_xy &vert) {
1340
+ extent_xy r;
1341
+ r.min .x = std::min (prev.x , next.x );
1342
+ r.min .y = std::min (prev.y , next.y );
1343
+ r.max .x = std::max (prev.x , next.x );
1344
+ r.max .y = std::max (prev.y , next.y );
1345
+
1346
+ if (vert.y < r.min .y || vert.y > r.max .y ) {
1347
+ return raycast_result::NONE;
1348
+ }
1349
+ if (vert.x < r.min .x ) {
1350
+ if (vert.y != r.min .y && vert.y != r.max .y ) {
1351
+ return raycast_result::CROSS;
1352
+ }
1353
+ } else if (vert.x > r.max .x ) {
1354
+ if (r.min .y != r.max .y && r.min .x != r.max .x ) {
1355
+ return raycast_result::NONE;
1356
+ }
1357
+ }
1358
+ vertex_xy p = vert;
1359
+ vertex_xy a = prev;
1360
+ vertex_xy b = next;
1361
+
1362
+ if (a.y < b.y ) {
1363
+ vertex_xy t = a;
1364
+ a = b;
1365
+ b = t;
1366
+ }
1367
+ if (pteq (vert, a) || pteq (vert, b)) {
1368
+ return raycast_result::BOUNDARY;
1369
+ }
1370
+ if (a.y == b.y ) {
1371
+ if (a.x == b.x ) {
1372
+ return raycast_result::NONE;
1373
+ }
1374
+ if (p.y == b.y ) {
1375
+ if (!(p.x < r.min .x || p.x > r.max .x )) {
1376
+ return raycast_result::BOUNDARY;
1377
+ }
1378
+ }
1379
+ }
1380
+ if (a.x == b.x && vert.x == b.x ) {
1381
+ if (vert.y >= a.y && vert.y <= b.y ) {
1382
+ return raycast_result::BOUNDARY;
1383
+ }
1384
+ }
1385
+ if (collinear (a.x , a.y , b.x , b.y , p.x , p.y )) {
1386
+ if (p.x < r.min .x ) {
1387
+ if (r.min .y == r.max .y ) {
1388
+ return raycast_result::NONE;
1389
+ }
1390
+ } else if (p.x > r.max .x ) {
1391
+ return raycast_result::NONE;
1392
+ }
1393
+ return raycast_result::BOUNDARY;
1394
+ }
1395
+ if (p.y == a.y || p.y == b.y ) {
1396
+ p.y = nexttoward (p.y , INFINITY);
1397
+ }
1398
+ if (p.y < a.y || p.y > b.y ) {
1399
+ return raycast_result::NONE;
1400
+ }
1401
+ if (a.x > b.x ) {
1402
+ if (p.x >= a.x ) {
1403
+ return raycast_result::NONE;
1404
+ }
1405
+ if (p.x <= b.x ) {
1406
+ return raycast_result::CROSS;
1407
+ }
1408
+ } else {
1409
+ if (p.x >= b.x ) {
1410
+ return raycast_result::NONE;
1411
+ }
1412
+ if (p.x <= a.x ) {
1413
+ return raycast_result::CROSS;
1414
+ }
1415
+ }
1416
+ if ((p.y -a.y )/(p.x -a.x ) >= (b.y -a.y )/(b.x -a.x )) {
1417
+ return raycast_result::CROSS;
1418
+ }
1419
+ return raycast_result::NONE;
1420
+ }
1421
+
1283
1422
1284
1423
// TODO: Make robust
1285
1424
raycast_result raycast_fast (const vertex_xy &prev, const vertex_xy &next, const vertex_xy &vert) {
@@ -2410,7 +2549,7 @@ point_in_polygon_result prepared_geometry::contains(const vertex_xy &vert) const
2410
2549
}
2411
2550
2412
2551
// The end of this node is either the end of the current node, or the end of the level
2413
- const auto node_end = (stack[depth - 1 ] + 1 ) * NODE_SIZE;
2552
+ const auto node_end = (( stack[depth - 1 ] + 1 ) * NODE_SIZE) - 1 ;
2414
2553
const auto levl_end = index.level_array [depth].entry_count - 1 ;
2415
2554
2416
2555
const auto end = math::min (node_end, levl_end);
0 commit comments