Repository: josdejong/mathjs (js)

Issue:
math.intersect() is numerically unstable.
Example:
math.intersect([300, 90], [400, 97], [300, 130], [400, 125]) returns null.

The current code has a hard equivalence check which should probably be relaxed to a soft check within some epsilon.

Pull request:
Fix #540: math.intersect() is numerically unstable
This fixes the numerically unstable math.intersect() for the 2
dimensional case by computing whether the determinant is within
config.epsilon instead of a hard equality check.

Patch:
diff --git a/lib/function/geometry/intersect.js b/lib/function/geometry/intersect.js
index b248facd6..1ce3fa18a 100644
--- a/lib/function/geometry/intersect.js
+++ b/lib/function/geometry/intersect.js
@@ -1,7 +1,12 @@
 'use strict';
 
 function factory (type, config, load, typed) {
+
+  var abs = load(require('../arithmetic/abs'));
+  var add = load(require('../arithmetic/add'));
   var matrix = load(require('../../type/matrix/function/matrix'));
+  var multiply = load(require('../arithmetic/multiply'));
+  var subtract = load(require('../arithmetic/subtract'));
 
   /**
    * Calculates the point of intersection of two lines in two or three dimensions
@@ -46,7 +51,7 @@ function factory (type, config, load, typed) {
         if (!_2d(y)) { throw new TypeError('Array with 2 numbers expected for third argument'); }
         if (!_2d(z)) { throw new TypeError('Array with 2 numbers expected for fourth argument'); }
 
-        return _intersect2d(w[0], w[1], x[0], x[1], y[0], y[1], z[0], z[1]);
+        return _intersect2d(w, x, y, z);
       }
       else if (w.length === 3) {
         if (!_3d(w)) { throw new TypeError('Array with 3 numbers expected for first argument'); }
@@ -71,74 +76,66 @@ function factory (type, config, load, typed) {
     }
   });
 
-  return intersect;
-}
-
-function _2d(x) {
-  return x.length === 2 && typeof x[0] === 'number' && typeof x[1] === 'number';
-}
-
-function _3d(x) {
-  return x.length === 3 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number';
-}
+  function _2d(x) {
+    return x.length === 2 && typeof x[0] === 'number' && typeof x[1] === 'number';
+  }
 
-function _4d(x) {
-  return x.length === 4 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number' && typeof x[3] === 'number';
-}
+  function _3d(x) {
+    return x.length === 3 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number';
+  }
 
-function _intersect2d(x1, y1, x2, y2, x3, y3, x4, y4){
-  var d1343 = (x1 - x3)*(x4 - x3) + (y1 - y3)*(y4 - y3);
-  var d4321 = (x4 - x3)*(x2 - x1) + (y4 - y3)*(y2 - y1);
-  var d1321 = (x1 - x3)*(x2 - x1) + (y1 - y3)*(y2 - y1);
-  var d4343 = (x4 - x3)*(x4 - x3) + (y4 - y3)*(y4 - y3);
-  var d2121 = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
-  var ta = ( d1343*d4321 - d1321*d4343 ) / ( d2121*d4343 - d4321*d4321 );
-  var tb = ( d1343 + ta * d4321 ) / (d4343);
-
-  var pax = x1 + ta * (x2 - x1);
-  var pay = y1 + ta * (y2 - y1);
-  var pbx = x3 + tb * (x4 - x3);
-  var pby = y3 + tb * (y4 - y3);
-  if (pax === pbx && pay === pby){
-    return [pax, pay];
+  function _4d(x) {
+    return x.length === 4 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number' && typeof x[3] === 'number';
   }
-  else{
-    return null;
+
+  function _intersect2d(p1a, p1b, p2a, p2b){
+    var o1 = p1a;
+    var o2 = p2a;
+    var d1 = subtract(o1, p1b);
+    var d2 = subtract(o2, p2b);
+    var det = d1[0]*d2[1] - d2[0]*d1[1];
+    if (abs(det) < config.epsilon) {
+      return null;
+    }
+    var t = (d2[0]*o1[1] - d2[1]*o1[0] - d2[0]*o2[1] + d2[1]*o2[0]) / det;
+    return add(multiply(d1, t), o1); 
   }
-}
 
-function _intersect3d(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4){
-  var d1343 = (x1 - x3)*(x4 - x3) + (y1 - y3)*(y4 - y3) + (z1 - z3)*(z4 - z3);
-  var d4321 = (x4 - x3)*(x2 - x1) + (y4 - y3)*(y2 - y1) + (z4 - z3)*(z2 - z1);
-  var d1321 = (x1 - x3)*(x2 - x1) + (y1 - y3)*(y2 - y1) + (z1 - z3)*(z2 - z1);
-  var d4343 = (x4 - x3)*(x4 - x3) + (y4 - y3)*(y4 - y3) + (z4 - z3)*(z4 - z3);
-  var d2121 = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1);
-  var ta = ( d1343*d4321 - d1321*d4343 ) / ( d2121*d4343 - d4321*d4321 );
-  var tb = ( d1343 + ta * d4321 ) / (d4343);
-
-  var pax = x1 + ta * (x2 - x1);
-  var pay = y1 + ta * (y2 - y1);
-  var paz = z1 + ta * (z2 - z1);
-  var pbx = x3 + tb * (x4 - x3);
-  var pby = y3 + tb * (y4 - y3);
-  var pbz = z3 + tb * (z4 - z3);
-  if (pax === pbx && pay === pby && paz === pbz){
-    return [pax, pay, paz];
+  function _intersect3d(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4){
+    var d1343 = (x1 - x3)*(x4 - x3) + (y1 - y3)*(y4 - y3) + (z1 - z3)*(z4 - z3);
+    var d4321 = (x4 - x3)*(x2 - x1) + (y4 - y3)*(y2 - y1) + (z4 - z3)*(z2 - z1);
+    var d1321 = (x1 - x3)*(x2 - x1) + (y1 - y3)*(y2 - y1) + (z1 - z3)*(z2 - z1);
+    var d4343 = (x4 - x3)*(x4 - x3) + (y4 - y3)*(y4 - y3) + (z4 - z3)*(z4 - z3);
+    var d2121 = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1);
+    var ta = ( d1343*d4321 - d1321*d4343 ) / ( d2121*d4343 - d4321*d4321 );
+    var tb = ( d1343 + ta * d4321 ) / (d4343);
+
+    var pax = x1 + ta * (x2 - x1);
+    var pay = y1 + ta * (y2 - y1);
+    var paz = z1 + ta * (z2 - z1);
+    var pbx = x3 + tb * (x4 - x3);
+    var pby = y3 + tb * (y4 - y3);
+    var pbz = z3 + tb * (z4 - z3);
+    if (pax === pbx && pay === pby && paz === pbz){
+      return [pax, pay, paz];
+    }
+    else{
+      return null;
+    }
   }
-  else{
-    return null;
+
+  function _intersectLinePlane(x1, y1, z1, x2, y2, z2, x, y, z, c){
+    var t = (c - x1*x - y1*y - z1*z)/(x2*x + y2*y + z2*z - x1 - y1 - z1);
+    var px = x1 + t * (x2 - x1);
+    var py = y1 + t * (y2 - y1);
+    var pz = z1 + t * (z2 - z1);
+    return [px, py, pz];
+    // TODO: Add cases when line is parallel to the plane:
+    //       (a) no intersection,
+    //       (b) line contained in plane
   }
-}
 
-function _intersectLinePlane(x1, y1, z1, x2, y2, z2, x, y, z, c){
-  var t = (c - x1*x - y1*y - z1*z)/(x2*x + y2*y + z2*z - x1 - y1 - z1);
-  var px = x1 + t * (x2 - x1);
-  var py = y1 + t * (y2 - y1);
-  var pz = z1 + t * (z2 - z1);
-  return [px, py, pz];
-  // TODO: Add cases when line is parallel to the plane:
-  //       (a) no intersection,
-  //       (b) line contained in plane
+  return intersect;
 }
 
 exports.name = 'intersect';


Test patch:
diff --git a/test/function/geometry/intersect.test.js b/test/function/geometry/intersect.test.js
index eb46c88ea..c1a6b4790 100644
--- a/test/function/geometry/intersect.test.js
+++ b/test/function/geometry/intersect.test.js
@@ -7,6 +7,7 @@ describe('intersect', function() {
     assert.deepEqual(math.intersect([0, 0], [10, 10], [10, 0], [0, 10]), [5, 5]);
     assert.deepEqual(math.intersect(math.matrix([0, 0]), [10, 10], math.matrix([10, 0]), math.matrix([0, 10])), math.matrix([5, 5]));
     assert.deepEqual(math.intersect(math.matrix([0, 0]), math.matrix([10, 10]), math.matrix([10, 0]), math.matrix([0, 10])), math.matrix([5, 5]));
+    assert.deepEqual(math.intersect([300, 90], [400, 97], [300, 130], [400, 125]), [633.3333333333334, 113.33333333333334]);
   });
 
   it('should calculate the intersection point of two 3D lines', function() {