I strongly refer to that book or any other good book on computer algorithms for more details.

The method should be O(N*log N), both for the intitial point sorting, and for the building/traversing of the tree. Note that you could also just compare each and every line, but that's O(N^2). I don't think you can cut this down any further while some additional huristics that are based on human observation.

For graphics algorithms, check the Graphic Gems series of books. Many are raster-rendering, but vector-polygon processing is covered too.

The old standby was the CALGO Collected Algorithms of the ACM. I used to have that on Microfiche subscription, but I think it's online now.)

-- Bill

I like that way. Upon reading the question, I realized that the normal edge-sorting technique used in rendering does tell when lines cross, but it only samples the shape at the rastors. It makes sence, though, that you can compare intervals instead.

I think you might be able to simplify this problem by preprocessing your data set. You shouldn't be trying to do some 200 million intersection tests. A (human generated) map is usually made of lots of relatively small line segments, mainly because it has to convey a good amount of information while remaining readable. I would start by giving every line segment an id number and put them in a database (or flat file, hash is fine). You should be able to quickly find what length line segment is most popular and what range of lengths around that covers say 80% of all the lines. Also, if you search starting with the longer ones, it should save you time, that is if long segments intersect more lines than do short segments.

This characterization would give you a guide to creating a grid of cells across the map, so that you can sort your segments by spatial distribution, in effect reducing the number of pairs. It is possible that you could define some bounding boxes of your own (say around cities?) which would do this efficiently for many of the longer line segments as well.

You can probably also do away with the bounding box calls and that epsilon stuff. You can also mark each segment in the database which has already been found to intersecting some other line, in effect taking that line segment out of circulation.

One thing to note is your definition of intersection. I think the algorithm below will find two lines which touch at their endpoints to be intersecting lines, but your example of the lake seems to indicate that such a pair of line segments should not be said to intersect. In that case you can eliminate a lot of computation by doing such a test first.

There is another approach which might be interesting. You could select a minimum resolution which would be needed to represent your map reasonably well, and provide a two dimensional array to represent this canvas. By drawing all the segments on this map with a routine which draws differently upon intersection, you can "download" the final answer by simply reading off the canvas at the end. Possibly you could do something with C-based drawing functions and various paint modes (and GIMP might be cool for that indeed, using an additive drawing mode) or collision detection functions (say in OpenGL), but I'm thinking here of just writing a line segment rendering function in Perl which "draws" to a data structure.

The neat part is that each "pixel", being an array element, could contain a list of the id numbers of all lines which crossed that pixel, which would provide a very rich output. You only have to execute the drawing code 14000 times, which sounds pretty good compared to 200 million of anything, but I'm not sure which approach is faster. Perhaps the most interesting thing is that you can run it at arbitrary resolutions, so you could for example run it once at a low resolution, then divide the canvas into an arbitrary number of cells and then perhaps use mathematical approach afterwards, or perhaps skip making cells and just use math to check the questionable pairs.

The pixel rendering approach is more interesting if you can grasp the entire image visually in a gestalt, for example by using colored filters. You might like to check out Gimp even just for fun since it is quite Perl-compatible.

Lastly I think there is one more approach to consider, if this is a one-shot job. If you just draw all the lines and open up the image in Gimp or Photoshop, you can use a combination of using the magic wand, paint bucket, and zoom tool to try solving this by hand.

You can also assign each line a unique RGB color combination and read its id with an eye-dropper (or Perl call if using Gimp).

Finding the exact point of intersection is too much work if
all we care about is whether two lines intersect at all. 
The intersection, if any, can be found by examining the 
signs of the two cross products 
(p2 - p0) × (p1 - p0) and (p3 - p0)× (p1 - p0).
The line_intersect() subroutine returns a simple true or
false value indicating whether two lines intersect:

# line_intersect( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 )
# Returns true if the two lines defined by these points intersect.
# In borderline cases, it relies on epsilon to decide.
 sub line_intersect {
  my ( $x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3 ) = @_;
  my @box_a = bounding_box( 2, $x0, $y0, $x1, $y1 );
  my @box_b = bounding_box( 2, $x2, $y2, $x3, $y3 );
  # If even the bounding boxes do not intersect, give up right now.
  return 0 unless bounding_box_intersect( 2, @box_a, @box_b );
  # If the signs of the two determinants (absolute values or lengths
  # of the cross products, actually) are different, the lines
  # intersect.
  my $dx10 = $x1 - $x0;
  my $dy10 = $y1 - $y0;
  my $det_a = determinant( $x2 - $x0, $y2 - $y0, $dx10, $dy10 );
  my $det_b = determinant( $x3 - $x0, $y3 - $y0, $dx10, $dy10 );
  return 1 if $det_a < 0 and $det_b > 0 or
  $det_a > 0 and $det_b < 0;
  if ( abs( $det_a ) < epsilon ) {
   if ( abs( $det_b ) < epsilon ) {
    # Both cross products are "zero".
    return 1;
   } elsif (abs($x3-$x2)<epsilon and abs($y3-$y2)<epsilon){
     # The other cross product is "zero" and
     # the other vector (from (x2,y2) to (x3,y3))
     # is also "zero".
    return 1;
   }
  } elsif ( abs( $det_b < epsilon ) ) {
  # The other cross product is "zero" and
  # the other vector is also "zero".
  return 1 if abs( $dx10 ) < epsilon and abs( $dy10 ) < epsilon;
 }
 return 0; # Default is no intersection.
}

# determinant( $x0, $y0, $x1, $y1 )    [from p. 432]
# Computes the determinant given the four elements of a matrix
# as arguments.
#
sub determinant { $_[0] * $_[3] - $_[1] * $_[2] }
[download]

Fascinating question, but I don't have any hard answers.

I assume you're looking at the simpler line_intersect function on page 439 (1st ed). Have you run Devel::DProf on a semi-reasonable sample to see where the time is being burnt?

Glancing at the code, I would assume that bounding_box and determinant would be good candidats for memoization, which will get you the on the RAM side of the space/time trade-off. You do have a lot of RAM, right?

g r i n d e r

PDL (``Perl Data Language'') gives standard Perl the ability to compactly store and speedily manipulate the large N-dimensional data arrays which are the bread and butter of scientific computing.

Well, the first part of this problem to consider, that makes it much more difficult, is that you are actually talking about the intersection of line segments, not lines themselves. This defeats the elementary proposition that two lines in a plane WILL intersect somewhere if their slopes are inequal.

However, there is a "brute-force" approach that will eliminate most of the possible pairs to consider in a pairwise approach: instead of looking at the segments pairwise to determine whether they intersect, look at the bounding rectangles of the segments pairwise to see if they overlap. Then, run the pairwise segment comparison only on this smaller subset.

For the purpose of this proposition, suppose that in some orthogonal, two-dimensional coordinate system (can you tell I really was a geometer at one point in my life?), the endpoints of a line segment are described by (a,b) and (p,q). Then the bounding rectangle of the line segment {(a,b),(p,q)} is the rectangle described by the point-set {(a,b),(a,q),(p,q),(p,b)}. It is a trivial exercise to demonstrate that if two line segments intersect, their bounding rectangles overlap; what we are actually interested in then is the contrapositive of that result: if their bounding rectangles do not overlap, two line segments do not intersect.

Spud Zeppelin * spud@spudzeppelin.com

If your shoreline data already knows which point is next down the beach, this becomes a problem in graph topology - grouping connected pieces of the shoreline "network". Pick a point and walk the list till you see the same point again. Pick an unvisited point, either an island or a neighboring lake, and repeat. This algorithm is linear in the number of points.

Chapter 8 of the wolf book covers this, but is geared to more complex graphs.

After Compline,
Zaxo

Chapter 33.2 of the Cormen–Leiserson–Rivest–Stein: Introduction to Algorithms book gives a good algorithm for solving this.