#!/usr/bin/perl

######################################################################
# CPAN Modules
######################################################################

use YAML qw(Load);
use Data::Dumper;
use Carp;
use Benchmark;

use strict;

######################################################################
# A few globals
######################################################################

my $points; # hash ref holding the raw data {'1' => [X,Y], '2' => [X,Y], ...}
my $edges; # Array ref to keep refrences to all edges [edge, X1, Y1, X2, Y2] ...
my $neighborEdges; # hash ref to hold sorted neighbor edges {'1' => edge1, '2' => edge2, ...}

my $precision = 7;
my $delta = 10 ** (-$precision);

######################################################################
#
#       MAIN
#
######################################################################

# Read in the X,Y coordinates from supplied data
{
  local $/;
  open(my $fh, "dump_301.dat") or die "Unable to open file: $!\n";
  $points = Load( <$fh> );
  croak($@) if $@;
  close($fh);
}

my $start_pairs = Benchmark->new;
my $edges = find_unique_edges($points);
print "Unique edges ", scalar(@{$edges}), "\n";
print "Pairs parsed in :", timestr(timediff(Benchmark->new, $start_pairs)), "\n";

my $start_intersects = Benchmark->new;
my $neighborEdges = find_unintersected_edges($edges);
print "Unintersected edges ", scalar(keys %{$neighborEdges}), "\n";
print "Intersects parsed in :", timestr(timediff(Benchmark->new, $start_intersects)), "\n";

######################################################################
######################################################################

sub find_unique_edges
{

  my $MAP = shift;

  my @edges; # returnable data structure

  my $start_distance = Benchmark->new;
  my $found = {};

  # p1 is starting or anchor point of the line segment
  foreach my $p1 (@{$points}) {

    # p2 is end point of the line segment
    foreach my $p2 (@{$points}) {

      # We don't need to caculate if anchor and end are the same
      # or we have already seen these two pairs [reversed] before    
      unless ( ($p1 == $p2) || $found->{$p2}->{$p1} ) {

        # Compute the edge
        my $edge = sqrt(
          ($p1->[0] - $p2->[0])**2 +
          ($p1->[1] - $p2->[1])**2
        );

        # Push the whole thing on a stack
        push (@edges, [ $edge, $p1->[0], $p1->[1], $p2->[0], $p2->[1] ]);

        # Keep track of the edges we've already computed (no need to do so twice)
        $found->{$p2}->{$p1} = 1;

      }
    }
  }

  return(\@edges);
}

######################################################################
######################################################################

sub find_unintersected_edges
{

  my $PAIRS = shift;

  my $neighborEdges; # returnable data structure
  my $neighboorCNT = 1;
  foreach my $aref ( sort {$a->[0] <=> $b->[0] } @{$PAIRS}) {

    my $neighbor = 1;

    my $line_1_ref = [];
    push (@{$line_1_ref}, [$aref->[1],$aref->[2]]); # line 1 point A X,Y
    push (@{$line_1_ref}, [$aref->[3],$aref->[4]]); # line 1 point B X,Y

    for (my $n=1;$n<=$neighboorCNT;$n++) {
      if ($neighborEdges->{$n}) {

        my $points_ref;
        push (@{$points_ref}, @{$line_1_ref}); # line 1 points
        push (@{$points_ref}, [$neighborEdges->{$n}->[1],$neighborEdges->{$n}->[2]]); # line 2 point A X,Y
        push (@{$points_ref}, [$neighborEdges->{$n}->[3],$neighborEdges->{$n}->[4]]); # line 2 point B X,Y

        # If a intersect is found, set false and quit checking
        if (SegmentIntersection($points_ref)) { 
          $neighbor = 0;
          last;
        }

      }
    }


    if ($neighbor) { 
      $neighborEdges->{$neighboorCNT} = $aref;
      $neighboorCNT++;
    }
  }

  return($neighborEdges);
}

################################################################################
#
# Shamlessly stolen from Math::Geometry::Planar

sub SegmentIntersection {
  my $pointsref = $_[0];
  my @points = @$pointsref;
  my @p1 = @{$points[0]}; # p1,p2 = segment 1
  my @p2 = @{$points[1]};
  my @p3 = @{$points[2]}; # p3,p4 = segment 2
  my @p4 = @{$points[3]};

  my $n1 = Determinant(($p3[0]-$p1[0]),($p3[0]-$p4[0]),($p3[1]-$p1[1]),($p3[1]-$p4[1]));
  my $n2 = Determinant(($p2[0]-$p1[0]),($p3[0]-$p1[0]),($p2[1]-$p1[1]),($p3[1]-$p1[1]));
  my $d  = Determinant(($p2[0]-$p1[0]),($p3[0]-$p4[0]),($p2[1]-$p1[1]),($p3[1]-$p4[1]));
  if (abs($d) < $delta) {
    return 0; # parallel
  }
  if (!(($n1/$d < 1) && ($n2/$d < 1) &&
        ($n1/$d > 0) && ($n2/$d > 0))) {
    return 0;
  }
  return(1);
}

################################################################################
#  
#  The determinant for the matrix  | x1 y1 |
#                                  | x2 y2 |
#
# args : x1,y1,x2,y2
#

sub Determinant {
  my ($x1,$y1,$x2,$y2) = @_;
  return ($x1*$y2 - $x2*$y1);
}

##</code><code>##

%Time    Sec.     #calls   sec/call  F  name
58.21 1889.8621  28882029   0.000065     main::SegmentIntersection
30.60  993.5924        1  993.592408     main::find_unintersected_edges
11.12  361.0416  86646087   0.000004     main::Determinant
 0.05    1.7613        1   1.761285     main::find_unique_edges