A script to compute intersections between unit disks

Given a set of unit disks with known coordinates (centers) and ranges, the objective of the algorithm is to find all the partition areas (unit disk intersections) in the Euclidean plane. The following Perl script applies Helly's theorem on disks and prints a collection of unique segment areas as well as their coverage status.
The complexity of the approach is O(n3), where n the number of disks.

Input:
a hash table %ncoords containing the coordinates (values) of the centers of the disks (keys)
a hash table %range with keys the disk members and values the disk range

...
my $pt = "A"; # a letter is assigned to each segment area
my %point = (); 
foreach my $n (keys %ncoords){ 
    my ($x, $y) = @{$ncoords{$n}};
    push (@{$point{$pt}}, $n);
    foreach my $n_ (keys %ncoords){
        next if ($n_ == $n);
        my ($x_, $y_) = @{$ncoords{$n_}};
        if (distance($x, $x_, $y, $y_) < ($range{$n} + $range{$n_})){ 
            my ($x1, $y1, $x2, $y2) = circles_inters_points($n, $n_); 
            my $p1 = increment($pt); # subroutine to increase letters
            my $p2 = increment($p1);
            push (@{$point{$p1}}, ($n, $n_));
            push (@{$point{$p2}}, ($n, $n_));
            foreach my $nd (keys %ncoords){
                next if (($nd == $n) || ($nd == $n_));
                    my ($ndx, $ndy) = @{$ncoords{$nd}};
                if (distance($ndx, $x1, $ndy, $y1) < $range{$nd}){
                    push (@{$point{$p1}}, $nd);
                }
                if (distance($ndx, $x2, $ndy, $y2) < $range{$nd}){
                    push (@{$point{$p2}}, $nd);
                }
            }
            $pt = $p2++;
        }
        $pt++;
    }
    $pt++;
}

# remove redundant points
foreach my $pt (keys %point){
    next if (!exists $point{$pt});
    foreach my $pt_ (keys %point){
        next if ($pt_ eq $pt);
        next if ((scalar @{$point{$pt_}}) != (scalar @{$point{$pt}}));
        my $check = 1;
        foreach my $n (@{$point{$pt}}){
            if (!grep {$_ == $n} @{$point{$pt_}}){
                $check = 0;
                last;
            }
        }
        if ($check == 1){
            delete $point{$pt_};
        }
    }
}

foreach my $pt (keys %point){
    print "# $pt:";
    foreach my $n (@{$point{$pt}}){
        print "$n ";
    }
    print "\n";
}
sub increment {
    my $self = shift;
    $self++;
    return $self;
}

sub circles_inters_points{
    my ($n, $n_) = @_;
    my ($x0, $x1, $y0, $y1, $r0, $r1) = ($ncoords{$n}[0],$ncoords{$n_}[0],$ncoords{$n}[1],$ncoords{$n_}[1], $range{$n}, $range{$n_});

    my $d = distance($x0, $x1, $y0, $y1);
    my $a = ($r0**2 - $r1**2 + $d**2)/(2*$d);
    my ($x2, $y2) = ($x0 + $a*($x1-$x0)/$d, $y0 + $a*($y1-$y0)/$d);
    my $h = sqrt(abs($r0**2 - $a**2));
    my ($x3, $y3) = ($x2 + $h*($y1-$y0)/$d, $y2 - $h*($x1-$x0)/$d);
    my ($x4, $y4) = ($x2 - $h*($y1-$y0)/$d, $y2 + $h*($x1-$x0)/$d);

    return ($x3, $y3, $x4, $y4);
}

sub distance {
    my ($x1, $x2, $y1, $y2) = @_;
    return sqrt( ($x1-$x2)*($x1-$x2) + ($y1-$y2)*($y1-$y2) );
}
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License