#!/usr/bin/perl -w

use constant epsilon => 1e-14;

# bounding_box($d, @p [,@b])
#   Return the bounding box of the points @p in $d dimensions.
#   The @b is an optional initial bounding box: we can use this
#   to create a cumulative bounding box that includes boxes found
#   by earlier runs of the subroutine (this feature is used by
#   bounding_box_of_points()).
#
#   The bounding box is returned as a list.  The first $d elements
#   are the minimum coordinates, the last $d elements are the
#   maximum coordinates.

sub bounding_box {
    my ( $d, @bb ) = @_; # $d is the number of dimensions.
    # Extract the points, leave the bounding box.
    my @p = splice( @bb, 0, @bb - 2 * $d );

    @bb = ( @p, @p ) unless @bb;

    # Scan each coordinate and remember the extrema.
    for ( my $i = 0; $i < $d; $i++ ) {
        for ( my $j = 0; $j < @p; $j += $d ) {
            my $ij = $i + $j;
            # The minima.
            $bb[ $i      ] = $p[ $ij ] if $p[ $ij ] < $bb[ $i      ];
            # The maxima.
            $bb[ $i + $d ] = $p[ $ij ] if $p[ $ij ] > $bb[ $i + $d ];
        }
    }

    return @bb;
}


# bounding_box_intersect($d, @a, @b)
#   Return true if the given bounding boxes @a and @b intersect
#   in $d dimensions.  Used by line_intersection().

sub bounding_box_intersect {
    my ( $d, @bb ) = @_; # Number of dimensions and box coordinates.
    my @aa = splice( @bb, 0, 2 * $d ); # The first box.
    # (@bb is the second one.)

    # Must intersect in all dimensions.
    for ( my $i_min = 0; $i_min < $d; $i_min++ ) {
        my $i_max = $i_min + $d; # The index for the maximum.
        return 0 if ( $aa[ $i_max ] + epsilon ) < $bb[ $i_min ];
        return 0 if ( $bb[ $i_max ] + epsilon ) < $aa[ $i_min ];
    }

    return 1;
}

# line_intersection( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 )
#
#    Compute the intersection point of the line segments
#    (x0,y0)-(x1,y1) and (x2,y2)-(x3,y3).
#
#    Or, if given four arguments, they should be the slopes of the
#    two lines and their crossing points at the y-axis.  That is,
#    if you express both lines as y = ax+b, you should provide the
#    two 'a's and then the two 'b's.
#
#    line_intersection() returns either a triplet ($x, $y, $s) for the
#    intersection point, where $x and $y are the coordinates, and $s
#    is true when the line segments cross and false when the line
#    segments don't cross (but their extrapolated lines would).
#
#    Otherwise, it's a string describing a non-intersecting situation:
#      "out of bounding box"
#      "parallel"
#      "parallel collinear"
#      "parallel horizontal"
#      "parallel vertical"
#    Because of the bounding box checks, the cases "parallel horizontal"
#    and "parallel vertical" never actually happen.  (Bounding boxes
#    are discussed later in the chapter.)
#
sub line_intersection {
 my ( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 );

 if ( @_ == 8 ) {
     ( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 ) = @_;

     # The bounding boxes chop the lines into line segments.
     # bounding_box() is defined later in this chapter.
     my @box_a = bounding_box( 2, $x0, $y0, $x1, $y1 );
     my @box_b = bounding_box( 2, $x2, $y2, $x3, $y3 );

     # Take this test away and the line segments are
     # turned into lines going from infinite to another.
     # bounding_box_intersect() defined later in this chapter.
     return "out of bounding box"
         unless bounding_box_intersect( 2, @box_a, @box_b );
 } elsif ( @_ == 4 ) { # The parametric form.
     $x0 = $x2 = 0;
     ( $y0, $y2 ) = @_[ 1, 3 ];
     # Need to multiply by 'enough' to get 'far enough'.
     my $abs_y0 = abs $y0;
     my $abs_y2 = abs $y2;
     my $enough = 10 * ( $abs_y0 > $abs_y2 ? $abs_y0 : $abs_y2 );
     $x1 = $x3 = $enough;
     $y1 = $_[0] * $x1 + $y0;
     $y3 = $_[2] * $x2 + $y2;
 }

 my ($x, $y);  # The as-yet-undetermined intersection point.

 my $dy10 = $y1 - $y0; # dyPQ, dxPQ are the coordinate differences
 my $dx10 = $x1 - $x0; # between the points P and Q.
 my $dy32 = $y3 - $y2;
 my $dx32 = $x3 - $x2;

 my $dy10z = abs( $dy10 ) < epsilon; # Is the difference $dy10 "zero"?
 my $dx10z = abs( $dx10 ) < epsilon;
 my $dy32z = abs( $dy32 ) < epsilon;
 my $dx32z = abs( $dx32 ) < epsilon;

 my $dyx10;                            # The slopes.
 my $dyx32;

 $dyx10 = $dy10 / $dx10 unless $dx10z;
 $dyx32 = $dy32 / $dx32 unless $dx32z;

 # Now we know all differences and the slopes;
 # we can detect horizontal/vertical special cases.
 # E.g., slope = 0 means a horizontal line.

 unless ( defined $dyx10 or defined $dyx32 ) {
     return "parallel vertical";
 } elsif ( $dy10z and not $dy32z ) { # First line horizontal.
     $y = $y0;
     $x = $x2 + ( $y - $y2 ) * $dx32 / $dy32;
 } elsif ( not $dy10z and $dy32z ) { # Second line horizontal.
     $y = $y2;
     $x = $x0 + ( $y - $y0 ) * $dx10 / $dy10;
 } elsif ( $dx10z and not $dx32z ) { # First line vertical.
     $x = $x0;
     $y = $y2 + $dyx32 * ( $x - $x2 );
 } elsif ( not $dx10z and $dx32z ) { # Second line vertical.
     $x = $x2;
     $y = $y0 + $dyx10 * ( $x - $x0 );
 } elsif ( abs( $dyx10 - $dyx32 ) < epsilon ) {
     # The slopes are suspiciously close to each other.
     # Either we have parallel collinear or just parallel lines.

     # The bounding box checks have already weeded the cases
     # "parallel horizontal" and "parallel vertical" away.

     my $ya = $y0 - $dyx10 * $x0;
     my $yb = $y2 - $dyx32 * $x2;

     return "parallel collinear" if abs( $ya - $yb ) < epsilon;
     return "parallel";
 } else {
     # None of the special cases matched.
     # We have a "honest" line intersection.

     $x = ($y2 - $y0 + $dyx10*$x0 - $dyx32*$x2)/($dyx10 - $dyx32);
     $y = $y0 + $dyx10 * ($x - $x0);
 }

 my $h10 = $dx10 ? ($x - $x0) / $dx10 : ($dy10 ? ($y - $y0) / $dy10 : 1);
 my $h32 = $dx32 ? ($x - $x2) / $dx32 : ($dy32 ? ($y - $y2) / $dy32 : 1);

 return ($x, $y, $h10 >= 0 && $h10 <= 1 && $h32 >= 0 && $h32 <= 1);
}

print "@{[line_intersection( 1, 1,  5, 5,  1, 4,  4, 1 )]}\n";
print "@{[line_intersection( 1, 1,  5, 5,  2, 4,  7, 4 )]}\n";
print "@{[line_intersection( 1, 1,  5, 5,  3, 0,  3, 6 )]}\n";
print "@{[line_intersection( 1, 1,  5, 5,  5, 2,  7, 2 )]}\n";
print     line_intersection( 1, 1,  5, 5,  4, 2,  7, 5 ), "\n";
print     line_intersection( 1, 1,  5, 5,  3, 3,  6, 6 ), "\n";