mirror of
https://github.com/josdejong/mathjs.git
synced 2026-01-25 15:07:57 +00:00
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.
This commit is contained in:
Keunhong Park
parent
e36c46d3ff
commit
f6d021e2f9
@ -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,75 +76,67 @@ function factory (type, config, load, typed) {
|
||||
}
|
||||
});
|
||||
|
||||
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 _4d(x) {
|
||||
return x.length === 4 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number' && typeof x[3] === 'number';
|
||||
}
|
||||
|
||||
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];
|
||||
}
|
||||
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
|
||||
}
|
||||
|
||||
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 _4d(x) {
|
||||
return x.length === 4 && typeof x[0] === 'number' && typeof x[1] === 'number' && typeof x[2] === 'number' && typeof x[3] === '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];
|
||||
}
|
||||
else{
|
||||
return null;
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
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
|
||||
}
|
||||
|
||||
exports.name = 'intersect';
|
||||
exports.factory = factory;
|
||||
|
||||
@ -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() {
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user