Intersect of two planes - divide by zero

Question

I have following alghoritm to find line intersection of two planes:

public static function getIntersectOf2Planes ( self $P1 , self $P2 )
{
    /*  Line equation in parametric form:
            x = x0 + t*a
            y = y0 + t*b
            z = z0 + t*c
        */
    $x0 = ( $P1->B * $P2->D - $P2->B * $P1->D ) / ( $P1->A * $P2->B - $P2->A * $P1->B ) ;
    $a = ( $P1->B * $P2->C - $P2->B * $P1->C );
    $y0 = ( $P2->A * $P1->D - $P1->A * $P2->D ) / ( $P1->A * $P2->B - $P2->A * $P1->B ) ;
    $b = ( $P2->A * $P1->C - $P1->A * $P2->C );
    $z0 = 0;
    $c = (  $P1->A * $P2->B - $P2->A * $P1->B );

    $IntersectionLine = new Line3D( $x0, $a, $y0, $b, $z0, $c );

    return $IntersectionLine;
}

and it works fine, but to compute $x0 and $y0 i have to divide by:

( $P1->A * $P2->B - $P2->A * $P1->B )

and in some cases, the value of this expression is equal to zero, so I get an "dividing by zero" error :(

What should I do in this case?

I know, that the case when this expression is equal to zero, doesn't mean that there is no intersection, because it's happen when I have planes perpendicular to one of the axies.

For example for:

Plane1:
A = 0
B = 0
C = 100
D = 0

Plane2:
A = 50
B = 0
C = 0
D = -250

so the equation of line should exists.

PS I wrote this code with a view to: https://math.stackexchange.com/questions/2766615/line-by-two-planes-intersection?noredirect=1#comment5706281_2766615

Answer 1

In short, you have to implement the intersection algorithm for the case when (a1*b2 - a2*b1) = 0 (ie. when the planes are not independent when you set z=0).

To expand on that, first we need to understand how you got this far. First let us write down the equation of two planes:

a1x + b1y + c1z + d1 = 0

and

a2x + b2y + c2z + d2 = 0

When two planes intersect, the intersection is a line. So, the most usual way to solve that is by finding a point first on such a line and then figuring out its orientation (a, b, c) in your case. The orientation is a straight forward cross product. The intersection point is typically calculated by setting one of the co-ordinates to 0 and then solving the 2 linear equations that you get. In your code, this is done by setting:

z = 0

But this only works when the equations

a1x + b1y + d1 = 0 and a2x + b2y + d2 = 0

are capable of giving a solution for x and y, which is not the case when a1b2-a2b1=0. So In such cases, you can solve the same by setting either x or y to 0 which again gives two linear equations that you can solve to obtain a point on the line. Then you can compute the parametric form much like how you did. For example (setting y to 0):

x0 = (c1d2 - c2d1)/(a1c2 - a2c1)
y0 = 0
z0 = (a2d1 - a1d2)/(a1c2 - a2c1)

But for this to be a defined value you need to have (a1c2 - a2c1) to be non-zero. Does this help?