Function to extract intersections of crossing lines on axes

Question

Code in Perl is 5.18.2.

sub extract_crossing {
    my @x = @{ $_[0] }; my @y = @{ $_[1] };
    my @xcross =(); my @ycross =();
    for (my $i=0; $i<$#x; $i++) {
        my $k = ($y[$i] - $y[$i+1]) / ($x[$i] - $x[$i+1]);
        if($y[$i+1] * $y[$i] < 0) {
            my $xc = $x[$i+1] - $y[$i+1] / $k;
            push(@xcross, $xc);
        }
        if($x[$i+1] * $x[$i] < 0) {
            my $yc = $y[$i+1] - $x[$i+1] * $k;
            push(@ycross, $yc);
        }
    }
    return (\@xcross, \@ycross);
}

which extracts successfully crossing points with x- and y-axes. It looks first point where the product of two subsequent points is negative. If so, intersection point with corresponding axis.

However, I feel this function extraneous, since it is so basic operation.

How can you do this extraction better with default tools in Perl?

Answer 1

If List::MoreUtils qualifies as one of Perl's "default tools" as you say in a comment, Math::Geometry::Planar should qualify as well. Math::Geometry::Planar provides a number of handy functions for calculating the intersection of segments, rays, and lines, as well as functions for manipulating polygons, calculating distances, and other goodies.

When evaluating any solution, you should make sure it generates correct results for a number of inputs, including edge cases. Your original code has at least one bug (a divide-by-zero error for vertical line segments)...let's make sure SegmentLineIntersection from Math::Geometry::Planar works as expected:

use strict;
use warnings;

use Math::Geometry::Planar qw(SegmentLineIntersection);
use Test::More tests => 8;

my @x_axis = ( [0, 0], [1, 0] );
my @y_axis = ( [0, 0], [0, 1] );

is_deeply(
    SegmentLineIntersection([ [-1, 2], [2, -1], @x_axis ]),
    [1, 0],
    'Segment (-1, 2), (2, -1) intersects x-axis once at (1, 0)'
);

is_deeply(
    SegmentLineIntersection([ [-1, 2], [2, -1], @y_axis ]),
    [0, 1],
    'Segment (-1, 2), (2, -1) intersects y-axis once at (0, 1)'
);

is(
    SegmentLineIntersection([ [0, 1], [1, 1], @x_axis ]),
    0,
    'Horizontal segment above x-axis never intersects x-axis'
);

is(
    SegmentLineIntersection([ [1, 0], [1, 1], @y_axis ]),
    0,
    'Vertical segment to the right of y-axis never intersects y-axis'
);

is(
    SegmentLineIntersection([ [0, 0], [1, 0], @x_axis ]),
    0,
    'Horizontal segment on x-axis returns false (intersects infinite times)'
);

is(
    SegmentLineIntersection([ [0, 0], [0, 1], @y_axis ]),
    0,
    'Vertical segment on y-axis returns false (intersects infinite times)'
);

is_deeply(
    SegmentLineIntersection([ [0, 0], [1, 1], @x_axis ]),
    [0, 0],
    'Segment beginning at origin intersects x-axis at (0, 0)'
);

is_deeply(
    SegmentLineIntersection([ [0, 0], [1, 1], @y_axis ]),
    [0, 0],
    'Segment beginning at origin intersects y-axis at (0, 0)'
);

Output:

1..8
ok 1 - Segment (-1, 2), (2, -1) intersects x-axis once at (1, 0)
ok 2 - Segment (-1, 2), (2, -1) intersects y-axis once at (0, 1)
ok 3 - Horizontal segment above x-axis never intersects x-axis
ok 4 - Vertical segment to the right of y-axis never intersects y-axis
ok 5 - Horizontal segment on x-axis returns false (intersects infinite times)
ok 6 - Vertical segment on y-axis returns false (intersects infinite times)
not ok 7 - Segment beginning at origin intersects x-axis at (0, 0)
#   Failed test 'Segment beginning at origin intersects x-axis at (0, 0)'
#   at geometry line 49.
#     Structures begin differing at:
#          $got = '0'
#     $expected = ARRAY(0x1b1f088)
not ok 8 - Segment beginning at origin intersects y-axis at (0, 0)
#   Failed test 'Segment beginning at origin intersects y-axis at (0, 0)'
#   at geometry line 55.
#     Structures begin differing at:
#          $got = '0'
#     $expected = ARRAY(0x1b1f010)
# Looks like you failed 2 tests of 8.

Looks like our last two tests failed: apparently segments with one end on a line don't count as intersecting (this is the case in your original algorithm as well). I'm not a geometry expert, so I can't evaluate whether this is a bug or mathematically correct.

Calculating intercepts for multiple segments

The following function returns the x-intercepts for multiple connected line segments. The implementation for calculating y-intercepts would be nearly identical. Note that if a pair of segments meet exactly on the axis, it doesn't count as an intercept, as in your original function. That may or may not be desirable.

use strict;
use warnings;

use Math::Geometry::Planar qw(SegmentLineIntersection);
use Test::Exception;
use Test::More tests => 3;

sub x_intercepts {
    my ($points) = @_;

    die 'Must pass at least 2 points' unless @$points >= 2;

    my @intercepts;
    my @x_axis = ( [0, 0], [1, 0] );

    foreach my $i (0 .. $#$points - 1) {
        my $intersect = SegmentLineIntersection([ @$points[$i, $i + 1], @x_axis ]);
        push @intercepts, $intersect if $intersect;
    }

    return \@intercepts;
}

dies_ok { x_intercepts([ [0, 0] ]) } 'Dies with < 2 points';

is_deeply(
    x_intercepts([ [-1, -1], [1, 1], [1, -1] ]),
    [ [0, 0], [1, 0] ],
    'Intersects x-axis at (0, 0) and (1, 0)'
);

is_deeply(
    x_intercepts([ [-1, -1], [0, 0], [1, 1] ]),
    [],
    "No intercept when segments start or end on x-axis but don't cross it"
);

Output:

1..3
ok 1 - Dies with < 2 points
ok 2 - Intersects x-axis at (0, 0) and (1, 0)
ok 3 - No intercept when segments start or end on x-axis but don't cross it

Note that this implementation accepts a single array reference of points, where a point is a reference to a two-element array, instead of separate array references of x- and y-coordinates. I think this is a little bit more intuitive.

Answer 2

This is one alternative solution. The code in the question calculates the value of $k every iteration, that's way too often, since you'll only need it if you want to store the value in one of the return arrays. This code loops over the indices from 0 to $#x-1 and makes use of the fact that the index is always stored in $_. This way, you some sort of lazy evaluation. Also I'm not sure if it is a bug, but in the first iteration of the for loop in the OP the comparison is $y[-1] * $y[0] > 0. This code doesn't do that.

use 5.010;              # for the // iterator
sub extract_crossing2 {
  my @x = @{ $_[0] }; my @y = @{ $_[1] };
  my (@xcross, @ycross);
  # "lazily" calculate $k,
  # there is a possibility of division by zero here! maybe catch that 
  # with if (defined $@){…}
  my $get_k = sub {
    eval {($y[$_] - $y[$_+1]) / ($x[$_] - $x[$_+1])}
  };
  foreach (0..$#x-1){
    my $k;        # only gets set if needed
    push @xcross, ($x[$_] - $y[$_]) / ($k  = $get_k->()) if $y[$_] * $y[$_+1] < 0;
    push @ycross, ($y[$_] - $x[$_]) * ($k // $get_k->()) if $x[$_] * $x[$_+1] < 0;
  }
  return \(@xcross, @ycross);
}

There might be a nice solution with the pair* routines provided by List::MoreUtils and List::Util.

Edit: As ThisSuitIsBlackNot has pointed out, there is a possibility of a division by zero error in the question's code. I didn't work around that error.