http://qs1969.pair.com?node_id=480510


in reply to Millions of line segment intersection calcs: Looking for speed tips

How does the set of "unintersected edges" computed by your algorithm differ from a Delaunay triangulation? If the difference is insignificant for your application (or can be quickly computed), you might be better off starting with a fast Delaunay triangulation (e.g., using Qhull) and adapting its output to your needs.

Update: Take a look at Aurenhammer and Xu's "Optimal Triangulations". It points out that the greedy triangulation – what you are computing – is not necessarily the minimum-weight triangulation (MWT). The paper also notes, however, that for uniformly distributed point sets both the greedy and Delaunay triangulations yield satisfactory approximations to the MWT.

Second update: The following code uses Qhull and computes an (approximate) DT-derived solution almost instantly:

#!/usr/bin/perl

# solve.pl

use warnings;
use strict;

use YAML;
use File::Temp 'tempfile';
use List::Util 'sum';

my $points        = YAML::Load( do { local $/; <> } );

my $triangles     = compute_delaunay_triangulation( $points );
my $segs          = unique_segments( $triangles );
my @seg_endpoints = map [ @{$points}[@$_] ], @$segs;
my @seg_lengths   = map seg_length($_), @seg_endpoints;

my $count = @$segs;
my $mean  = sum(@seg_lengths) / $count;

print "edge count          = $count\n";
print "average edge length = $mean\n\n";
print "edge endpoints follow in YAML format\n\n";

print YAML::Dump( @seg_endpoints );

sub compute_delaunay_triangulation {

    # use Qhull's qdelaunay for this

    my ($points) = @_;

    my $fh = tempfile();
    die "cannot open tempfile: $!" unless $fh;

    my $pid = open my $write_to_child, "|-";
    die "cannot fork: $!" unless defined $pid;

    if ($pid) {  # parent
        print $write_to_child "2\n", scalar @$points, "\n",
            map "@$_\n", @$points;
        close $write_to_child || warn "kid exited $?";
    } else {     # child
        open STDOUT, '>&', $fh  or die "could not redirect stdout: $!"
+;
        close $fh               or die "error on close of tempfile: $!
+";
        exec qw(qdelaunay Qt i) or die "could not exec qdelaunay: $!";
    }

    wait;

    seek $fh, 0, 0 or die "could not seek to start of tempfile: $!";
    my $results = do { local $/; <$fh> };
    close $fh;

    (undef, my @triangles) = split /\n/, $results;

    return [ map [split], @triangles ];
}

sub unique_segments {
    my ($triangles) = @_;
    my %segs;
    for my $tri (@$triangles) {
        for ([0,1],[1,2],[2,0]) {
            my @points = sort @{$tri}[@$_];
            $segs{"@points"} = \@points;
        }
    }
    [ values %segs ];
}

sub seg_length {
    my ($seg) = @_;
    my ($x,$y,$x1,$y1) = ( $seg->[0][0], $seg->[0][1]
                         , $seg->[1][0], $seg->[1][1] );
    sqrt( ($x - $x1) ** 2 + ($y - $y1) ** 2);
}

[download]

Example: solve dump_301.yml. Output: 

edge count          = 889
average edge length = 123.835184709118

edge endpoints follow in YAML format

---
-
  - 2590
  - 3183
-
  - 2707
  - 3257
---
-
  - 2560
  - 2850
-
  - 2586
  - 2869

... other edges omitted ...

[download]

Cheers,
Tom

Tom Moertel : Blog / Talks / CPAN / LectroTest / PXSL / Coffee / Movie Rating Decoder

Replies are listed 'Best First'.
Re: Delaunay triangulation
by drewhead (Beadle) on Aug 04, 2005 at 02:32 UTC
    How does the set of "unintersected edges" computed by your algorithm differ from a Delaunay triangulation?
    That's a good question. I'd never heard of Delaunay triangulation until your post, so I'm going to need to do some reading. If the answer is that it doesn't or maybe if it's even a better way to show what I'm really after (the idea of remoteness of a set of coordinates... vague enough? ;-) I may just go this route.

    I admit to being somewhat of a perl biggot, however; speed in this case is really more important than my love affair with perl. This code will be triggered by a web app submittal of a file with the coords in it. Infact, speed is probally more important than finite accuracy of the stats I'm trying to produce. The end result will be used to relativally diffrenciate sets of coordinates. Or I should say it is my hope that the people who will use these number see them that way. :)

    Thanks for the links and the example. I've grabbed qhull and installed it so I can play around in the next few days. I'll update when I learn some more so the next guy that does a search.

[reply]
      If you end up implementing Delaunay Triangulation, you may as well throw in Voronoi diagrams :-)
[reply]
Re: Delaunay triangulation
by anonymized user 468275 (Curate) on Aug 03, 2005 at 16:16 UTC
    I for one would agree you are right on both counts, but to paraphrase Shakespeare, 'Does a problem by any other name solve more sweetly?'

    Update: To clarify what I mean: the OP knows the problem and its solution but is looking for the fastest way to achieve it with perl.

    One world, one people

[reply]
[reply]
      To clarify what I mean: the OP knows the problem and its solution but is looking for the fastest way to achieve it with perl.

      How do you know what the OP knows and is looking for? Do you have some knowledge beyond what is in the original post? The OP says he is "willing to entertain ideas of differing algorithms." What makes you think otherwise?

      Cheers,
      Tom

      Tom Moertel : Blog / Talks / CPAN / LectroTest / PXSL / Coffee / Movie Rating Decoder

[reply]
      I tend to side with tmoertel. I once in ages past implemented my own naive algorithm for this very thing ( O(n^3) or O(n^4), probably the same way that the OP is doing, but in AutoLISP using AutoCAD), and although I knew there must be a better way, I didn't know at the time that just knowing the name Delaunay Triangulation would've helped in finding better algorithms (and it was long before Google and probably before Al Gore invented the internet). I've since been unable to find the time to implement this in perl :-(

      There was a start made on this in perl utilizing a C library, but all I see is a readme file.

[reply]

In Section Seekers of Perl Wisdom