File Coverage

blib/lib/Math/Geometry/Voronoi.pm

Criterion	Covered	Total	%
statement	120	127	94.4
branch	45	56	80.3
condition	21	36	58.3
subroutine	14	15	93.3
pod	3	8	37.5
total	203	242	83.8

package Math::Geometry::Voronoi;

use 5.008;
use strict;
use warnings;

our $VERSION = '1.3';

require XSLoader;
XSLoader::load('Math::Geometry::Voronoi', $VERSION);

use Params::Validate qw(validate ARRAYREF CODEREF);
use List::Util qw(min max sum);

use base 'Class::Accessor::Fast';
__PACKAGE__->mk_accessors(
                         qw(points lines edges vertices xmin ymin xmax ymax));

sub new {
    my $pkg  = shift;
    my %args = validate(@_, {points => {type => ARRAYREF}});
    my $self = bless({points => $args{points}}, $pkg);

    $self->sort_points();

    return $self;
}

# C code needs points sorted by y then by x and needs min and max for
# both - should provide a way for the client to provide this
sub sort_points {
    my $self   = shift;
    my $points = $self->points();

    @$points =
      sort { $a->[1] <=> $b->[1] || $a->[0] <=> $b->[0] } @$points;

    $self->ymin($points->[0][1]);
    $self->ymax($points->[-1][1]);

    my @x = map { $_->[0] } @$points;
    $self->xmin(min(@x));
    $self->xmax(max(@x));

    return;
}

sub compute {
    my $self = shift;

    my $result = compute_voronoi_xs($self->points,
                                    $self->xmin,
                                    $self->xmax,
                                    $self->ymin,
                                    $self->ymax);

    $self->lines($result->{lines});
    $self->vertices($result->{vertices});
    $self->edges($result->{edges});

    return;
}

sub cmp_verts_ab { return cmp_verts($a,$b); }

# a low x value is placed before a high x value. if both x values
# are the same, a high y value is placed before a low y value.
sub cmp_verts {
    return ($_[0]->[0] <=> $_[1]->[0] || $_[1]->[1] <=> $_[0]->[1] );
}


sub vert_inside_bounds {
    my ($self, $x,$y) = @_;
    return (
        $x >= $self->xmin and $x <= $self->xmax and
        $y >= $self->ymin and $y <= $self->ymax
    );
}


sub boundry_interesection_verts {
    my($self, $a,$b,$c) = @_;
    my $verts = [];

    if($b){
        my $v1 = [$self->xmin,($c-$a*$self->xmin)/$b];
        my $v2 = [$self->xmax,($c-$a*$self->xmax)/$b];
        push ( @$verts, $v1 ) if ( $self->vert_inside_bounds( @$v1 ) );
        push ( @$verts, $v2 ) if ( $self->vert_inside_bounds( @$v2 ) );
    }

    if($a){
        my $v1 = [($c-$b*$self->ymax)/$a,$self->ymax];
        my $v2 = [($c-$b*$self->ymin)/$a,$self->ymin];
        push ( @$verts, $v1 ) if ( $self->vert_inside_bounds( @$v1 ) );
        push ( @$verts, $v2 ) if ( $self->vert_inside_bounds( @$v2 ) );
    }

    $verts;
}


sub polygons {
    my $self = shift;
    my %args = validate(@_,
                        {normalize_vertices => {type     => CODEREF,
                                                optional => 1
                                               },
                        });
    my $points   = $self->points;
    my $lines    = $self->lines;
    my $edges    = $self->edges;
    my $vertices = $self->vertices;

    if (my $norm = $args{normalize_vertices}) {
        $vertices = [map { [$norm->($_->[0]), $norm->($_->[1])] } @$vertices];
    }

    my @edges_by_point;
    EDGE: foreach my $edge (@$edges) {
        my ($l, $v1, $v2) = @$edge;
        next EDGE if( $v1 == -1 and $v2 == -1 );
        my ($lon1, $lat1, $lon2, $lat2);

        my $ivs = $self->boundry_interesection_verts(@{$lines->[$l]});
        $ivs = [sort cmp_verts_ab @$ivs];

        if( my $norm = $args{normalize_vertices}) {
            $ivs = [map { [$norm->($_->[0]), $norm->($_->[1])] } @$ivs];
        }

        ($lat1,$lon1) = @{$vertices->[$v1]} if( $v1 != -1 );
        ($lat2,$lon2) = @{$vertices->[$v2]} if( $v2 != -1 );

        if( $v1 == -1 ) {
            next EDGE unless( @$ivs and $lat2 +0 == $lat2 and $lon2 +0 == $lon2 );
            if( cmp_verts( [$lat2,$lon2], $ivs->[0] ) > 0 ) {
                ($lat1,$lon1) = @{$ivs->[0]};
            } elsif( cmp_verts( [$lat2,$lon2], $ivs->[1] ) > 0 ) {
                ($lat1,$lon1) = @{$ivs->[1]};
            } else {
                next EDGE;
            }
        }
        if( $v2 == -1 ) {
            next EDGE unless( @$ivs and $lat1 +0 == $lat1 and $lon1 +0 == $lon1 );
            if( cmp_verts( [$lat1,$lon1], $ivs->[1] ) < 0 ) {
                ($lat2,$lon2) = @{$ivs->[1]};
            } elsif( cmp_verts( [$lat1,$lon1], $ivs->[0] ) < 0 ) {
                ($lat2,$lon2) = @{$ivs->[0]};
            } else {
                next EDGE;
            }
        }

        # if any of the coords are NaN things break.
        next EDGE if( grep {$_ +0 != $_ } ($lat1,$lon1,$lat2,$lon2));

        my ($p1, $p2) = ($lines->[$l][3], $lines->[$l][4]);

        if ($p1 != -1 and $p2 != -1) {
            foreach my $p ($p1, $p2) {
                push @{$edges_by_point[$p]}, [$lat1, $lon1, $lat2, $lon2];
            }
        }
    }

    my @polygons;
    foreach my $p (0 .. $#$points) {
        my $stack = $edges_by_point[$p];
        next unless $stack;
        # can't make a polygon with less than 2 edges
        next unless @$stack >= 2;

        my @poly = ();
        foreach my $this ( @$stack ) {
            if(
                !grep { $_->[0] == $this->[0] && $_->[1] == $this->[1] } @poly
                and $this->[0] +0 == $this->[0]
                and $this->[1] +0 == $this->[1]
            ) {
                push @poly, [$this->[0],$this->[1]];
            }
            if(
                !grep { $_->[0] == $this->[2] && $_->[1] == $this->[3] } @poly
                and $this->[2] +0 == $this->[2]
                and $this->[3] +0 == $this->[3]
            ) {
                push @poly, [$this->[2],$this->[3]];
            }
        }

        #TODO: if this point is the closest point to a corner...
        # add that corner as a vert on this poly

        # sort poly's verts (anti?) clockwise around the point $points->[$p];
        @poly = sort {
            my($lat1,$lon1) = (
                $a->[0] - $points->[$p]->[0],
                $a->[1] - $points->[$p]->[1]
            );

            my($lat2,$lon2) = (
                $b->[0] - $points->[$p]->[0],
                $b->[1] - $points->[$p]->[1]
            );
            return atan2($lon1,$lat1) <=> atan2($lon2,$lat2);
        } @poly;

        # make a list of the first points
        push @polygons, [$p, map { [$_->[0], $_->[1]] } @poly];
    }

    return @polygons;
}

sub _dump_poly {
    my $poly = shift;
    return
      "[ \n\t" . join(", \n\t", map { "[$_->[0],$_->[1]]" } @$poly) . " ]\n";
}

1;
__END__