davideberly / GeometricTools

A collection of source code for computing in the fields of mathematics, geometry, graphics, image analysis and physics.
Boost Software License 1.0
1.08k stars 202 forks source link

Plane3Plane3 intersection algebraic mistake #86

Closed Greezzee closed 5 months ago

Greezzee commented 5 months ago

Plane specification formula is $Ax + By + Cz + D = 0$ Written using normals is $Dot(\vec{r}, \vec{N}) + D = 0$ So, in this function $Dot(\vec{N_0}, \vec{L}(t)) = -D_0$ and $Dot(\vec{N_1}, \vec{L}(t)) = -D_1$, so minus is missing

davideberly commented 5 months ago

Please read the comments immediately after the pragma-once statement in Hyperplane.h. My implementation of a plane with unit-length normal U and origin P is Dot(U,X-P) = 0. In terms of a plane constant c rather than specifying a plane origin, the equation is Dot(U,X) = c, where c = Dot(U,P). If U = (A,B,C) and X=(x,y,z), my plane is Ax+By+C*z-c = 0.

I believe my plane-plane intersection code is correct. For example,

    Plane3<double> plane0({ 0.0, 0.0, 1.0 }, 1.0);  // z = 1
    Plane3<double> plane1({ 0.0, 1.0, 0.0 }, -1.0);  // y = -1
    FIQuery<double, Plane3<double>, Plane3<double>> query{};
    auto result = query(plane0, plane1);
    // result.intersect = true
    // result.isLine = true
    // result.line.origin = (0, -1, 1) 
    // result.line.direction = (-1, 0, 0)
    // result.plane = <invalid, all members are zero>
    // NOTE: Points on the line of intersection have y = -1 and z = +1.

    // If I negate the right-hand sides of c0 and c1 as you suggest,
    // the line of intersection has points with y = +1 and z = -1,
    // but these are not my planes.
    double dot = Dot(plane0.normal, plane1.normal);
    double invDet = (double)1 / ((double)1 - dot * dot);
    double c0 = -(plane0.constant - dot * plane1.constant) * invDet;
    double c1 = -(plane1.constant - dot * plane0.constant) * invDet;
    result.intersect = true;
    result.isLine = true;
    result.line.origin = c0 * plane0.normal + c1 * plane1.normal;
    result.line.direction = UnitCross(plane0.normal, plane1.normal);
    // result.intersect = true
    // result.isLine = true
    // result.line.origin = (0, 1, -1)
    // result.line.direction = (-1, 0, 0)
    // result.plane = <invalid, all members are zero>