From 09b657b95d9d62079dca21908127e157c84703fc Mon Sep 17 00:00:00 2001
From: Manuel Baclet <mbaclet@gmail.com>
Date: Mon, 24 Jul 2017 15:30:34 +0200
Subject: [PATCH] Improve epsilon logic

---
 lib/epsilon.js | 74 ++++++++++++++++++++++++++++----------------------
 1 file changed, 41 insertions(+), 33 deletions(-)

diff --git a/lib/epsilon.js b/lib/epsilon.js
index 2d7d68d..5815d16 100644
--- a/lib/epsilon.js
+++ b/lib/epsilon.js
@@ -24,29 +24,28 @@ function Epsilon(eps){
 			var By = right[1];
 			var Cx = pt[0];
 			var Cy = pt[1];
-			return (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax) >= -eps;
+			var ABx = Bx - Ax;
+			var ABy = By - Ay;
+			var AB = Math.sqrt(ABx * ABx + ABy * ABy);
+			// algebraic distance of 'pt' to ('left', 'right') line is:
+			// [ABx * (Cy - Ay) - ABy * (Cx - Ax)] / AB
+			return ABx * (Cy - Ay) - ABy * (Cx - Ax) >= -eps * AB;
 		},
 		pointBetween: function(p, left, right){
 			// p must be collinear with left->right
 			// returns false if p == left, p == right, or left == right
+			if (my.pointsSame(p, left) || my.pointsSame(p, right)) return false;
 			var d_py_ly = p[1] - left[1];
 			var d_rx_lx = right[0] - left[0];
 			var d_px_lx = p[0] - left[0];
 			var d_ry_ly = right[1] - left[1];
 
 			var dot = d_px_lx * d_rx_lx + d_py_ly * d_ry_ly;
-			// if `dot` is 0, then `p` == `left` or `left` == `right` (reject)
-			// if `dot` is less than 0, then `p` is to the left of `left` (reject)
-			if (dot < eps)
-				return false;
-
+			// dot < 0 is p is to the left of 'left'
+			if (dot < 0) return false;
 			var sqlen = d_rx_lx * d_rx_lx + d_ry_ly * d_ry_ly;
-			// if `dot` > `sqlen`, then `p` is to the right of `right` (reject)
-			// therefore, if `dot - sqlen` is greater than 0, then `p` is to the right of `right` (reject)
-			if (dot - sqlen > -eps)
-				return false;
-
-			return true;
+			// dot <= sqlen is p is to the left of 'right'
+			return dot <= sqlen;
 		},
 		pointsSameX: function(p1, p2){
 			return Math.abs(p1[0] - p2[0]) < eps;
@@ -71,7 +70,18 @@ function Epsilon(eps){
 			var dy1 = pt1[1] - pt2[1];
 			var dx2 = pt2[0] - pt3[0];
 			var dy2 = pt2[1] - pt3[1];
-			return Math.abs(dx1 * dy2 - dx2 * dy1) < eps;
+			var n1 = Math.sqrt(dx1 * dx1 + dy1 * dy1);
+			var n2 = Math.sqrt(dx2 * dx2 + dy2 * dy2);
+			// Assuming det(u, v) = 0, we have:
+			// |det(u + u_err, v + v_err)| = |det(u + u_err, v + v_err) - det(u,v)|
+			// =|det(u, v_err) + det(u_err. v) + det(u_err, v_err)|
+			// <= |det(u, v_err)| + |det(u_err, v)| + |det(u_err, v_err)|
+			// <= N(u)N(v_err) + N(u_err)N(v) + N(u_err)N(v_err)
+			// <= eps * (N(u) + N(v) + eps)
+			// We have N(u) ~ N(u + u_err) and N(v) ~ N(v + v_err).
+			// Assuming eps << N(u) and eps << N(v), we end with:
+			// |det(u + u_err, v + v_err)| <= eps * (N(u + u_err) + N(v + v_err))
+			return Math.abs(dx1 * dy2 - dx2 * dy1) <= eps * (n1 + n2);
 		},
 		linesIntersect: function(a0, a1, b0, b1){
 			// returns false if the lines are coincident (e.g., parallel or on top of each other)
@@ -98,7 +108,9 @@ function Epsilon(eps){
 			var bdy = b1[1] - b0[1];
 
 			var axb = adx * bdy - ady * bdx;
-			if (Math.abs(axb) < eps)
+			var n1 = Math.sqrt(adx * adx + ady * ady);
+			var n2 = Math.sqrt(bdx * bdx + bdy * bdy);
+			if (Math.abs(axb) <= eps * (n1 + n2))
 				return false; // lines are coincident
 
 			var dx = a0[0] - b0[0];
@@ -106,38 +118,34 @@ function Epsilon(eps){
 
 			var A = (bdx * dy - bdy * dx) / axb;
 			var B = (adx * dy - ady * dx) / axb;
+			var pt = [
+				a0[0] + A * adx,
+				a0[1] + A * ady
+			];
 
 			var ret = {
 				alongA: 0,
 				alongB: 0,
-				pt: [
-					a0[0] + A * adx,
-					a0[1] + A * ady
-				]
+				pt: pt
 			};
 
 			// categorize where intersection point is along A and B
-
-			if (A <= -eps)
-				ret.alongA = -2;
-			else if (A < eps)
+			if (my.pointsSame(pt, a0))
 				ret.alongA = -1;
-			else if (A - 1 <= -eps)
-				ret.alongA = 0;
-			else if (A - 1 < eps)
+			else if (my.pointsSame(pt, a1))
 				ret.alongA = 1;
-			else
+			else if (A < 0)
+				ret.alongA = -2;
+			else if (A > 1)
 				ret.alongA = 2;
 
-			if (B <= -eps)
-				ret.alongB = -2;
-			else if (B < eps)
+			if (my.pointsSame(pt, b0))
 				ret.alongB = -1;
-			else if (B - 1 <= -eps)
-				ret.alongB = 0;
-			else if (B - 1 < eps)
+			else if (my.pointsSame(pt, b1))
 				ret.alongB = 1;
-			else
+			else if (B < 0)
+				ret.alongB = -2;
+			else if (B > 1)
 				ret.alongB = 2;
 
 			return ret;