#!/usr/bin/perl sub closest_points { my ( @p ) = @_; return () unless @p and @p % 2 == 0; my $unsorted_x = [ map { $p[ 2 * $_ ] } 0..$#p/2 ]; my $unsorted_y = [ map { $p[ 2 * $_ + 1 ] } 0..$#p/2 ]; # Compute the permutation and ordinal indices. # X Permutation Index. # # If @$unsorted_x is (18, 7, 25, 11), @xpi will be (1, 3, 0, 2), # e.g., $xpi[0] == 1 meaning that the $sorted_x[0] is in # $unsorted_x->[1]. # # We do this because we may now sort @$unsorted_x to @sorted_x # and can still restore the original ordering as @sorted_x[@xpi]. # This is needed because we will want to sort the points by x and y # but might also want to identify the result by the original point # indices: "the 12th point and the 45th point are the closest pair". my @xpi = sort { $unsorted_x->[ $a ] <=> $unsorted_x->[ $b ] } 0..$#$unsorted_x; # Y Permutation Index. # my @ypi = sort { $unsorted_y->[ $a ] <=> $unsorted_y->[ $b ] } 0..$#$unsorted_y; # Y Ordinal Index. # # The ordinal index is the inverse of the permutation index: If # @$unsorted_y is (16, 3, 42, 10) and @ypi is (1, 3, 0, 2), @yoi # will be (2, 0, 3, 1), e.g. $yoi[0] == 1 meaning that # $unsorted_y->[0] is the $sorted_y[1]. my @yoi; @yoi[ @ypi ] = 0..$#ypi; # Recurse to find the closest points. my ( $p, $q, $d ) = _closest_points_recurse( [ @$unsorted_x[@xpi] ], [ @$unsorted_y[@xpi] ], \@xpi, \@yoi, 0, $#xpi ); my $pi = $xpi[ $p ]; # Permute back. my $qi = $xpi[ $q ]; ( $pi, $qi ) = ( $qi, $pi ) if $pi > $qi; # Smaller id first. return ( $pi, $qi, $d ); } sub _closest_points_recurse { my ( $x, $y, $xpi, $yoi, $x_l, $x_r ) = @_; # $x, $y: array references to the x- and y-coordinates of the points # $xpi: x permutation indices: computed by closest_points_recurse() # $yoi: y ordering indices: computed by closest_points_recurse() # $x_l: the left bound of the currently interesting point set # $x_r: the right bound of the currently interesting point set # That is, only points $x->[$x_l..$x_r] and $y->[$x_l..$x_r] # will be inspected. my $d; # The minimum distance found. my $p; # The index of the other end of the minimum distance. my $q; # Ditto. my $N = $x_r - $x_l + 1; # Number of interesting points. if ( $N > 3 ) { # We have lots of points. Recurse! my $x_lr = int( ( $x_l + $x_r ) / 2 ); # Right bound of left half. my $x_rl = $x_lr + 1; # Left bound of right half. # First recurse to find out how the halves do. my ( $p1, $q1, $d1 ) = _closest_points_recurse( $x, $y, $xpi, $yoi, $x_l, $x_lr ); my ( $p2, $q2, $d2 ) = _closest_points_recurse( $x, $y, $xpi, $yoi, $x_rl, $x_r ); # Then merge the halves' results. # Update the $d, $p, $q to be the closest distance # and the indices of the closest pair of points so far. if ( $d1 < $d2 ) { $d = $d1; $p = $p1; $q = $q1 } else { $d = $d2; $p = $p2; $q = $q2 } # Then check the straddling area. # The x-coordinate halfway between the left and right halves. my $x_d = ( $x->[ $x_lr ] + $x->[ $x_rl ] ) / 2; # The indices of the "potential" points: those point pairs # that straddle the area and have the potential to be closer # to each other than the closest pair so far. # my @xi; # Find the potential points from the left half. # The left bound of the left segment with potential points. my $x_ll; if ( $x_lr == $x_l ) { $x_ll = $x_l } else { # Binary search. my $x_ll_lo = $x_l; my $x_ll_hi = $x_lr; do { $x_ll = int( ( $x_ll_lo + $x_ll_hi ) / 2 ); if ( $x_d - $x->[ $x_ll ] > $d ) { $x_ll_lo = $x_ll + 1; } elsif ( $x_d - $x->[ $x_ll ] < $d ) { $x_ll_hi = $x_ll - 1; } } until $x_ll_lo > $x_ll_hi or ( $x_d - $x->[ $x_ll ] < $d and ( $x_ll == 0 or $x_d - $x->[ $x_ll - 1 ] > $d ) ); } push @xi, $x_ll..$x_lr; # Find the potential points from the right half. # The right bound of the right segment with potential points. my $x_rr; if ( $x_rl == $x_r ) { $x_rr = $x_r } else { # Binary search. my $x_rr_lo = $x_rl; my $x_rr_hi = $x_r; do { $x_rr = int( ( $x_rr_lo + $x_rr_hi ) / 2 ); if ( $x->[ $x_rr ] - $x_d > $d ) { $x_rr_hi = $x_rr - 1; } elsif ( $x->[ $x_rr ] - $x_d < $d ) { $x_rr_lo = $x_rr + 1; } } until $x_rr_hi < $x_rr_lo or ( $x->[ $x_rr ] - $x_d < $d and ( $x_rr == $x_r or $x->[ $x_rr + 1 ] - $x_d > $d ) ); } push @xi, $x_rl..$x_rr; # Now we know the potential points. Are they any good? # This gets kind of intense. # First sort the points by their original indices. my @x_by_y = @$yoi[ @$xpi[ @xi ] ]; my @i_x_by_y = sort { $x_by_y[ $a ] <=> $x_by_y[ $b ] } 0..$#x_by_y; my @xi_by_yi; @xi_by_yi[ 0..$#xi ] = @xi[ @i_x_by_y ]; my @xi_by_y = @$yoi[ @$xpi[ @xi_by_yi ] ]; my @x_by_yi = @$x[ @xi_by_yi ]; my @y_by_yi = @$y[ @xi_by_yi ]; # Inspect each potential pair of points (the first point # from the left half, the second point from the right). for ( my $i = 0; $i <= $#xi_by_yi; $i++ ) { my $i_i = $xi_by_y[ $i ]; my $x_i = $x_by_yi[ $i ]; my $y_i = $y_by_yi[ $i ]; for ( my $j = $i + 1; $j <= $#xi_by_yi; $j++ ) { # Skip over points that can't be closer # to each other than the current best pair. last if $xi_by_y[ $j ] - $i_i > 7; # Too far? my $y_j = $y_by_yi[ $j ]; my $dy = $y_j - $y_i; last if $dy > $d; # Too tall? my $x_j = $x_by_yi[ $j ]; my $dx = $x_j - $x_i; next if abs( $dx ) > $d; # Too wide? # Still here? We may have a winner. # Check the distance and update if so. my $d3 = sqrt( $dx**2 + $dy**2 ); if ( $d3 < $d ) { $d = $d3; $p = $xi_by_yi[ $i ]; $q = $xi_by_yi[ $j ]; } } } } elsif ( $N == 3 ) { # Just three points? No need to recurse. my $x_m = $x_l + 1; # Compare the square sums and leave the sqrt for later. my $s1 = ($x->[ $x_l ]-$x->[ $x_m ])**2 + ($y->[ $x_l ]-$y->[ $x_m ])**2; my $s2 = ($x->[ $x_m ]-$x->[ $x_r ])**2 + ($y->[ $x_m ]-$y->[ $x_r ])**2; my $s3 = ($x->[ $x_l ]-$x->[ $x_r ])**2 + ($y->[ $x_l ]-$y->[ $x_r ])**2; if ( $s1 < $s2 ) { if ( $s1 < $s3 ) { $d = $s1; $p = $x_l; $q = $x_m } else { $d = $s3; $p = $x_l; $q = $x_r } } elsif ( $s2 < $s3 ) { $d = $s2; $p = $x_m; $q = $x_r } else { $d = $s3; $p = $x_l; $q = $x_r } $d = sqrt $d; } elsif ( $N == 2 ) { # Just two points? No need to recurse. $d = sqrt(($x->[ $x_l ]-$x->[ $x_r ])**2 + ($y->[ $x_l ]-$y->[ $x_r ])**2); $p = $x_l; $q = $x_r; } else { # Less than two points? Strange. return ( ); } return ( $p, $q, $d ); } @clopo = closest_points( 1, 2, 2, 5, 3, 1, 3, 3, 4, 5, 5, 1, 5, 6, 6, 4, 7, 4, 8, 1 ), "\n"; print "@clopo\n";