source: subversion/applications/rendering/tilesAtHome/close-areas.pl @ 11553

Last change on this file since 11553 was 11110, checked in by matthiasj, 11 years ago

move all *.pm files into lib and add lib to the library search paths

File size: 27.8 KB
Line 
1#!/usr/bin/perl
2
3# A program to form closed areas from incomplete sets of bordering
4# segments, by making an assumption about which side of the segment
5# the enclosed area lies on.
6#
7# Currently only supports "coastline" ways but may be extended in
8# the future.
9#
10# A detailed discussion should be in the Wiki:
11# http://wiki.openstreetmap.org/index.php/Tiles%40home/Dev/Interim_Coastline_Support
12#
13# Written by Frederik Ramm <frederik@remote.org> - public domain
14
15use lib './lib';
16use Math::Trig;
17use Math::Complex;
18use tahproject;
19use strict;
20
21my $nodecount = 0;
22my $waycount = 0;
23my $relcount = 0;
24my $halfpi = pi()/2;
25my $twopi = pi()*2;
26
27my $LimitY = pi();
28my $LimitY2 = -pi();
29my $RangeY = $twopi;
30my $debug=0;
31my $make_open_ways=0;
32my $segments;
33my $nodes;
34
35# The direction in which areas are closed - cw (clockwise), or ccw
36# (counter-clockwise). For OSM we normally use cw, which means that
37# area boundaries are drawn such that the area is "to the right" of
38# the segment (i.e. coastline segments have water on their right).
39#
40# "ccw" is unused but supported in case somebody needs it.
41my $direction="cw";
42
43my @coastline_segments;
44
45my $tilex;
46my $tiley; 
47my $zoom;
48my $minlat;
49my $minlon;
50my $maxlat;
51my $maxlon;
52my $border_crossed;
53
54if (scalar(@ARGV) == 2 or scalar(@ARGV) == 3)
55{
56    ($tilex, $tiley, $zoom) = @ARGV;
57    if (!$zoom) {$zoom = 12;}
58    ($maxlat, $minlat) = Project($tiley, $zoom);
59    ($minlon, $maxlon) = ProjectL($tilex, $zoom);
60}
61elsif (scalar(@ARGV) == 4)
62{
63    ($minlat, $minlon, $maxlat, $maxlon) = @ARGV;
64}
65else
66{
67    die "need either 'tilex tiley [zoom]' (defaults to zoom=12) or 'minlat minlon maxlat maxlon' on cmdline";
68}
69
70my $sigma = 0.01;
71
72my $helpernodes =
73    { 0 => [ $maxlat + $sigma, $maxlon + $sigma], 
74      1 => [ $maxlat + $sigma, $minlon - $sigma ], 
75      2 => [ $minlat - $sigma, $minlon - $sigma ], 
76      3 => [ $minlat - $sigma, $maxlon + $sigma ] };
77
78my $copy = 1;
79my $waybuf;
80my $is_coastline;
81my @seglist;
82my $last_node_ref; 
83my $segcount = 0; 
84
85while(<STDIN>)
86{
87    while(/(<[^'"<>]*((["'])[^\3]*?\3[^<>"']*)*>)/og)
88    {
89        my $xmltag = $1;
90        if ($xmltag =~ /^\s*<node.*id=["'](\d+)['"].*lat=["']([0-9.Ee-]+)["'].*lon=["']([0-9.Ee-]+)["']/)
91        {
92            $nodes->{$1} = { "lat" => $2, "lon" => $3 };
93        }
94        elsif ($xmltag =~ /^\s*<way /)
95        {
96            $copy = 0;
97            undef $waybuf;
98            undef @seglist;
99            undef $is_coastline;
100            undef $last_node_ref;
101        }
102        elsif($xmltag =~ /^\s*<\/osm/)
103        {
104            last;
105        }
106        elsif($xmltag =~ /^\s*<(osm.*)>/)
107        {
108            # If frollo encounters an empty file, it outputs <osm foo />. Detect this and exit
109            if (substr($1, -1) eq "/" )
110            { print $xmltag; exit }
111           
112            if (/version=['"](.*?)['"]/)
113            {
114                die ("close-areas.pl does not support version $1") unless ($1 eq "0.5");
115            }
116        }
117
118        if ($copy)
119        {
120            print $xmltag."\n";
121        }
122        elsif ($xmltag =~ /^\s*<tag k=['"]natural["'].*v=["']coastline['"]/)
123        {
124            $is_coastline = 1; 
125        }
126        elsif ($xmltag =~ /^\s*<nd ref=['"](\d+)["']/)
127        {
128            $waybuf .= $xmltag . "\n";
129            if (defined($last_node_ref))
130            {
131                push(@seglist, { "from" => $last_node_ref, "to" => $1 });
132            }
133            $last_node_ref = $1;
134        }
135        elsif ($xmltag =~ /^\s*<\/way/)
136        {
137            $copy = 1;
138            # non-coastal ways are written and forgotten
139            if (!$is_coastline)
140            {
141                print $waybuf;
142                print $xmltag . "\n";
143            }
144            # coastal ways are forgotten too, but their segments are kept.
145            else
146            {
147                foreach my $seg(@seglist)
148                {
149                    $segments->{++$segcount} = $seg;
150                }
151            }
152        }
153        else
154        {
155            $waybuf .= $xmltag . "\n";
156        }
157    }
158}
159
160# file fully read. delete non-coastal segments (we have printed them
161# all, so no need for us to keep them), and check coastal segments for
162# intersection with the bounding box.
163
164foreach my $seg(keys(%$segments))
165{
166    my $fromnode = $nodes->{$segments->{$seg}->{"from"}};
167    my $tonode = $nodes->{$segments->{$seg}->{"to"}};
168
169    if (!(defined($fromnode) && defined($tonode)))
170    {
171        delete $segments->{$seg};
172        print "delete segment - incomplete $seg\n" if ($debug);
173        next;
174    }
175
176    # returns 0, 1, or 2 points.
177    # (may probably return 3 or 4 points in freak cases)
178    my $intersect = compute_bbox_intersections($fromnode, $tonode);
179    printf "intersections for $seg: %d\n", scalar(@$intersect) if ($debug);
180
181    if (!scalar(@$intersect))
182    {
183        # this segment has no intersections with the bounding box.
184        if (!node_is_inside($fromnode))
185        {
186            # this segment is fully outside the bounding box, and of
187            # no concern to us.
188            delete $segments->{$seg};
189            print "delete $seg fully out\n" if ($debug);
190        }
191    }
192    elsif (scalar(@$intersect) == 1)
193    {
194        # this segments enters OR exits the bounding box. find out which,
195        # and tag accordingly.
196        if (node_is_inside($fromnode))
197        {
198            $segments->{$seg}->{"exit_intersect"} = $intersect->[0];
199            $segments->{$seg}->{"exit_intersect_angle"} =
200                compute_angle_from_bbox_center($intersect->[0]);
201        }
202        else
203        {
204            $segments->{$seg}->{"entry_intersect"} = $intersect->[0];
205            $segments->{$seg}->{"entry_intersect_angle"} =
206                compute_angle_from_bbox_center($intersect->[0]);
207        }
208    }
209    else
210    {
211        # this segment enters AND exits the bounding box. as intersection
212        # points are ordered by distance from the segment's origin, we
213        # assume that the very first intersection is the entry point and
214        # the very last is the exit point.
215        #   
216        # FIXME: segments like this - probably a very long one cutting right
217        # through the box or a short one cutting diagonally at one edge -
218        # are very little tested.
219        $segments->{$seg}->{"entry_intersect_angle"} =
220            compute_angle_from_bbox_center($intersect->[0]);
221        $segments->{$seg}->{"exit_intersect_angle"} =
222            compute_angle_from_bbox_center($intersect->[scalar(@$intersect)-1]);
223        $segments->{$seg}->{"entry_intersect"} = $intersect->[0];
224        $segments->{$seg}->{"exit_intersect"} = $intersect->[scalar(@$intersect)-1];
225    }
226}
227
228# if no coastline segments are present, switch over to
229# special handler to decide whether to draw a blue tile.
230
231if (!scalar(%$segments))
232{
233    ## instead of filling the tile here, we do it at the end...
234    $border_crossed = 0;
235    goto ENDOSM;
236}
237
238# now start building artificial ways for coastline segments.
239#
240# strategy:
241# 1. grab any one segment. if none available, we're done.
242# 2. find segment that starts where previous segment ends,
243#    repeat until none found OR back where we started.
244# 3. if back where we started, write the circular list of segments,
245#    delete them, and repeat from 1.
246# 4. if none found AND last position is inside drawing area,
247#    remove processed segments without writing, and repeat from 1.
248# 5. if none found AND last position is outside drawing area,
249#    search clockwise outside of drawing area for another
250#    segment that begins here, create articifial segment joining
251#    latest position and segment's start node, and continue with
252#    step 2.
253#
254# Special rule for creating artificial segments: artificial segments
255# must never intersect with drawing area. use artificial nodes outside
256# the four corners of the drawing area to route segments if necessary.
257
258# copy list of segment ids
259my $available_segs;
260grep { $available_segs->{$_} = 1 } keys(%$segments);
261
262while (scalar(%$available_segs))
263{
264    my @seglist;
265    # grab any segment as the first one
266    my @tmp = keys(%$available_segs);
267    my $currentseg = shift(@tmp);
268    printf "GRAB %d (from %d to %d)\n", $currentseg, $segments->{$currentseg}->{"from"}, $segments->{$currentseg}->{"to"} if ($debug);
269    delete $available_segs->{$currentseg};
270    #printf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
271    my $currentnode = $segments->{$currentseg}->{"to"};
272    push (@seglist, $currentseg);
273
274STRING:
275    # find a segment that begins where the previous one ended
276    while ($currentseg)
277    {
278        printf "SEARCH begin at %d\n", $currentnode if ($debug);
279        undef $currentseg;
280        foreach my $seg(keys(%$available_segs))
281        {
282            #printf "  TEST seg %d begin at %d\n", $seg, $segments->{$seg}->{"from"} if ($debug);
283            if ($segments->{$seg}->{"from"} == $currentnode)
284            {
285                printf "  FOUND seg %d begin at %d\n", $seg, $segments->{$seg}->{"from"} if ($debug);
286                push (@seglist, $seg);
287                $currentseg = $seg;
288                $currentnode = $segments->{$currentseg}->{"to"};
289                delete $available_segs->{$seg};
290                #printf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
291                last;
292            }
293        }
294    }
295
296    # no more segments found. do we have a loop?
297    if ($currentnode == $segments->{$seglist[0]}->{"from"})
298    {
299        printf("LOOP\n") if ($debug);
300        # loop detected. store the segment list for later output,
301        # and move on trying to connect the rest.
302        push(@coastline_segments, \@seglist);
303        next;
304    }
305
306    my $lastseg = @seglist[scalar(@seglist)-1];
307    my $exit_angle = $segments->{$lastseg}->{"exit_intersect_angle"};
308
309    # are we inside the drawable area? (we are if the last segment did not
310    # exit the drawable area.)
311    # if yes, give up as we obviously have an imcomplete piece of coastline.
312    if (!defined($exit_angle))
313    {
314        printf("NOT OUTSIDE\n") if ($debug);
315        # this is a debug option that allows one to detect the incomplete
316        # way by looking at the output file with JOSM etc.; it is not intended
317        # for production use!
318        if ($make_open_ways)
319        {
320            make_way(\@seglist, 1);
321        }
322        next;
323    }
324
325    # "else" case: we are outside the drawable area and want
326    # to find another segment that begins outside the drawable
327    # area and enters it.
328   
329    my $segs_to_check;
330    $segs_to_check = [];
331    foreach my $seg(keys(%$available_segs))
332    {
333        push(@$segs_to_check, $seg) if defined($segments->{$seg}->{"entry_intersect"});
334    }
335
336    # we will also accept the first segment of the current way
337    # if it is outside. this is a special case as, being used
338    # already, it isn't in the $segments list any
339    # more
340   
341    push(@$segs_to_check, $seglist[0]) 
342        if defined($segments->{$seglist[0]}->{"entry_intersect"});
343
344    # sort all segments entering the drawable area by angle from area
345    # centrepoint to the point where they enter the area.
346    my @sorted_segs_to_check;
347    @sorted_segs_to_check = sort
348        { $segments->{$a}->{"entry_intersect_angle"} <=> 
349          $segments->{$b}->{"entry_intersect_angle"} } @$segs_to_check;
350
351    @sorted_segs_to_check = reverse @sorted_segs_to_check if ($direction eq "cw");
352
353    # find the nearest entering segment.
354    my $found;
355    $found = 0;
356    foreach my $seg(@sorted_segs_to_check)
357    {
358        if ($direction eq "cw")
359        {
360            next if ($segments->{$seg}->{"entry_intersect_angle"}) > $exit_angle;
361        }
362        else
363        {
364            next if ($segments->{$seg}->{"entry_intersect_angle"}) < $exit_angle;
365        }
366        printf("use seg %d angle %f\n", $seg, $segments->{$seg}->{"entry_intersect_angle"}) if ($debug);
367        $found = $seg;
368        last;
369    }
370    if (!$found)
371    {
372        foreach my $seg(@sorted_segs_to_check)
373        {
374            printf("use (2) seg %d angle %f\n", $seg, $segments->{$seg}->{"entry_intersect_angle"}) if ($debug);
375            $found = $seg;
376            last;
377        }
378    }
379
380    # if no segment remains outside, give up
381    if (!$found)
382    {
383        printf("NO SEG OUTSIDE\n") if ($debug);
384        # this is a debug option that allows one to detect the incomplete
385        # way by looking at the output file with JOSM etc.; it is not intended
386        # for production use!
387        if ($make_open_ways)
388        {
389            make_way(\@seglist, 1);
390        }
391        next;
392    }
393   
394    # at this point, we have a list of segments that ends with a segment
395    # leaving the drawable area, and we also have the next segment where
396    # the coastline comes back into the drawable area. we need to connect
397    # them with an artifical segment - or more than one.
398
399    # If a coastline leaves the visible area at the top (side 1) and comes
400    # back in from the right (side 0), then we must not make a direct connection
401    # between the points for fear of clipping our viewport; instead, we must
402    # "hop" over the top right corner. Same for other sides. The "helpernodes"
403    # hash contains nodes to be used for hopping from one side to the next.
404    #
405    # (an extreme case of this is coastline leaving and entering at the same
406    # side but leaving south of where it enters - need to go around all four
407    # corners then, for a total of 5 artificial segments.)
408
409    $border_crossed = 1;
410    my $height = $maxlat-$minlat;
411    my $width = $maxlon-$minlon;
412
413    # exit_angle is already defined
414    my $entry_angle = $segments->{$found}->{"entry_intersect_angle"};
415
416    my $exit_side = $segments->{$lastseg}->{"exit_intersect"}->{"side"};
417    my $entry_side = $segments->{$found}->{"entry_intersect"}->{"side"};
418
419    printf("exit angle %s entry angle %s\n", $exit_angle, $entry_angle) if $debug;
420    printf("exit side %d entry side %d\n", $exit_side, $entry_side) if $debug;
421
422    # the following two blocks (similar but not identical for clockwise/
423    # counter-clockwise) will find out whether we need to "go around corners"
424    # and add 0-4 segments if needed.
425    if ($direction eq "ccw")
426    {
427        # $min_once is for the special case where entry and exit side are
428        # identical but we still need to go all the way around the box.
429        my $min_once;
430        if ($exit_side == 0 and $entry_side == 0) 
431        {
432            # Take into account that the angle flips from 2*pi() to zero at this side
433            $min_once = (($exit_angle > $entry_angle and ($exit_angle - $entry_angle) < pi()) ||
434                         ($exit_angle < pi() and $entry_angle > pi()));
435        }
436        else 
437        {
438            $min_once = ($exit_angle > $entry_angle); 
439        }
440        printf("min_once=%d\n", $min_once) if $debug;
441        for (my $i = $exit_side;; $i++)
442        {
443            $i=0 if ($i==4);
444            last if ($i==$entry_side && !$min_once);
445            $min_once = 0;
446            my $newnode = make_node($helpernodes->{$i}->[0], $helpernodes->{$i}->[1]);
447            my $newseg = make_seg($currentnode, $newnode);
448            $currentnode = $newnode;
449            push(@seglist, $newseg);
450        }
451    }
452    else
453    {
454        my $min_once;
455        if ($exit_side == 0 and $entry_side == 0)
456        {
457            $min_once = (($exit_angle < $entry_angle and ($entry_angle - $exit_angle) < pi()) ||
458                         ($exit_angle > pi() and $entry_angle < pi()));
459        }
460        else
461        {
462            $min_once = ($exit_angle < $entry_angle);
463        }
464        printf("min_once=%d\n", $min_once) if $debug;
465        for (my $i = $exit_side;; $i--)
466        {
467            $i=3 if ($i==-1);
468            last if ($i==$entry_side && !$min_once);
469            $min_once = 0;
470            my $helper = $i-1;
471            $helper = 3 if ($helper == -1);
472            my $newnode = make_node($helpernodes->{$helper}->[0], $helpernodes->{$helper}->[1]);
473            my $newseg = make_seg($currentnode, $newnode);
474            $currentnode = $newnode;
475            push(@seglist, $newseg);
476        }
477    }
478
479    my $newseg = make_seg($currentnode, $segments->{$found}->{"from"});
480    push(@seglist, $newseg);
481
482    # if the segment we have found is our start segment (which we added
483    # as a special case earlier!), we have a closed way.
484    if ($found == $seglist[0])
485    {
486        printf("CLOSED\n") if ($debug);
487        push(@coastline_segments, \@seglist);
488        next;
489    }
490
491    # else we just plod on.
492    push (@seglist, $found);
493    $currentseg = $found;
494    $currentnode = $segments->{$found}->{"to"};
495    delete $available_segs->{$found};
496    #nprintf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug);
497    goto STRING;
498}
499
500ENDOSM:
501
502# if we had no coastline intersection with the bounding box, but coastline
503# was present, we have an island situation and need to add a blue background.
504# FIXME what if we have a "little lake" situation?
505
506unless( $border_crossed )
507{
508  my $state = lookup_handler($helpernodes, $tilex, $tiley, $zoom);
509  if( $state eq "10" )
510  {
511    # sea
512    addBlueRectangle($helpernodes);
513  }
514  elsif ( $state eq "01" )
515  {
516    #land
517  }
518  else
519  { 
520    my %temp = ("00"=>0, "10"=>0, "01"=>0, "11"=>0);;
521    $temp{lookup_handler($helpernodes, $tilex-1, $tiley, $zoom)}++;
522    $temp{lookup_handler($helpernodes, $tilex+1, $tiley, $zoom)}++;
523    $temp{lookup_handler($helpernodes, $tilex, $tiley-1, $zoom)}++;
524    $temp{lookup_handler($helpernodes, $tilex, $tiley+1, $zoom)}++;
525
526    if( $temp{"10"} > $temp{"01"} )
527    {
528      addBlueRectangle($helpernodes);
529    }
530    elsif ( ($state eq "11") and ( $temp{"01"} == 0 ) ) 
531    # if the tile is marked coast but no land near, assume it's a group of islands instead of lakes.
532    {
533      # coast
534      addBlueRectangle($helpernodes);
535    }
536    else
537    {
538      #land
539    }
540  }
541}
542
543if (scalar(@coastline_segments) == 1)
544{
545    make_way($coastline_segments[0]);
546}
547else
548{
549    make_multipolygon(\@coastline_segments);
550}
551print "</osm>\n";
552
553sub make_node
554{
555    my ($lat, $lon) = @_;
556    my $id = --$nodecount;
557    print "<node id='$id' lat='$lat' lon='$lon' />\n";
558    # save node in array for later processing
559    $nodes->{$id}={"lat"=>$lat, "lon"=>$lon};
560    return $id;
561}
562
563sub make_seg
564{
565    my ($from, $to) = @_;
566    my $id = ++$segcount;
567    $segments->{$id} = { "from" => $from, "to" => $to };
568    return $id;
569}
570
571sub make_way
572{
573    my ($seglist, $open) = @_;
574    my $id = --$waycount;
575    print "<way id='$id'>\n";
576    print "  <tag k='natural' v='coastline' />\n";
577    print "  <tag k='created-with' v='close-areas.pl' />\n";
578    print "  <tag k='close-areas.pl:debug' v='open way' />\n" if ($open);
579
580    my $first = 1;
581    foreach my $seg(@$seglist) 
582    { 
583        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"from"} if ($first);
584        printf "  <nd ref='%s' />\n", $segments->{$seg}->{"to"};
585        $first = 0;
586    }
587
588    print "</way>\n";
589    return $id;
590}
591
592sub make_multipolygon
593{
594    my ($seglistlist) = @_;
595    return unless scalar(@$seglistlist);
596
597    # we will create a list of multipolygons. each multipolygon is a list of
598    # ways, with the first being the "outer" polygon.
599    my $multipolygons = [];
600       
601    # first create all the ways.
602    my $way_ids = [];
603    for (my $i=0; $i<scalar(@$seglistlist); $i++)
604    {
605        $way_ids->[$i] = make_way($seglistlist->[$i]);
606    }
607
608    # now find out which way is contained in which. since our polygons do not
609    # intersect, we simply take one node of a way and check whether this node
610    # is contained in any of the other ways; if yes, the whole way is contained.
611    my $contained_in = [];
612    my $children = [];
613    for (my $i=0; $i<scalar(@$seglistlist); $i++)
614    {
615        my $pt = $nodes->{$segments->{$seglistlist->[$i]->[0]}->{"from"}};
616        for (my $j = 0; $j < scalar(@$seglistlist); $j++) # fixme really need to check all?
617        {
618            next if ($j==$i);
619            if (polygon_contains_point($seglistlist->[$j], $pt))
620            {
621                $contained_in->[$i] = $j; 
622                $children->[$j]++;
623                last;
624            }
625        }
626    }
627
628    # now write multipolygons
629    for (my $i=0; $i<scalar(@$seglistlist); $i++)
630    {
631        # ways that don't have children do not trigger a relation
632        next unless $children->[$i];
633        my $id = --$relcount;
634        print "<relation id='$id'>\n";
635        print "  <tag k='type' v='multipolygon' />\n";
636        printf "  <member type='way' ref='%s' role='outer' />\n", $way_ids->[$i];
637        for (my $j = 0; $j < scalar(@$seglistlist); $j++) 
638        {
639            if (defined($contained_in->[$j]) && $contained_in->[$j] == $i)
640            {
641                printf "  <member type='way' ref='%s' role='inner' />\n", $way_ids->[$j];
642            }
643        }
644        print "</relation>\n";
645    }
646}
647
648# index lookup by Martijn van Oosterhout
649sub lookup_handler
650{
651    my ($helpernodes, $x, $y, $zoom) = @_;
652    # make it use z12 x,y coordinates
653    # this looks up the most upper left z12 tile in zoom<12. This probably
654    # needs to be made smarter for islands etc.
655    my $tilex = $x*(2**(12-$zoom));
656    my $tiley = $y*(2**(12-$zoom));
657    my $tileoffset = ($tiley * (2**12)) + $tilex;
658    my $fh;
659    open($fh, "<", "png2tileinfo/oceantiles_12.dat") or die;
660    seek $fh, int($tileoffset / 4), 0; 
661    my $buffer;
662    read $fh, $buffer, 1;
663    $buffer = substr( $buffer."\0", 0, 1 );
664    $buffer = unpack "B*", $buffer;
665    my $str = substr( $buffer, 2*($tileoffset % 4), 2 );
666    close($fh);
667
668    print("lookup handler finds: $str\n") if $debug;
669
670    # $str eq "00" => unknown (not yet checked)
671    # $str eq "01" => known land
672    # $str eq "10" => known sea
673    # $str eq "11" => known edge tile
674
675    return $str; 
676}
677
678sub addBlueRectangle
679{
680    my $helpernodes = shift;
681    my @n;
682    my @s;
683    for (my $i=0; $i<4; $i++)
684    {
685        $n[$i] = make_node($helpernodes->{$i}->[0], 
686                         $helpernodes->{$i}->[1]);
687    }
688    for (my $i=0; $i<4; $i++)
689    {
690        if ($direction eq "ccw")
691        {
692            $s[$i] = make_seg($n[$i], $n[($i+1)%4]);
693        }
694        else
695        {
696            $s[3-$i] = make_seg($n[($i+3)%4], $n[($i+2)%4]);
697        }
698    }
699    push(@coastline_segments, \@s);
700}
701
702# this takes two points (hashes with lat/lon keys) as arguments and returns
703# an array reference to an array containing up to four points (hashes with
704# lat/lon keys and an added "side" key, where 0=right 1=top 2=left 3=bottom)
705# denoting the intersections of the line formed by the input points with the
706# bounding box.
707#
708# 0 results - segment is fully inside or fully outside bounding box
709# 1 result  - segment begins inside, ends outside or vice versa
710# 2 results - segment begins outside and ends outside, but cuts through
711#             bounding box
712# 3 results - can't think how this can happen
713# 4 results - as case 2 but segment is exactly at the bounding box diagonal,
714#             thus intersecting with all four edges (freak case)
715#
716# If there is more than one result, results are sorted by distance from the
717# first point.
718sub compute_bbox_intersections
719{
720    my ($from, $to) = @_;
721    my $result = [];
722    my $point;
723
724    # div by zero FIXME
725    my $latd = ($to->{"lat"} - $from->{"lat"});
726    my $lond = ($to->{"lon"} - $from->{"lon"});
727
728    # segment's bbox
729    my $s_minlon = $from->{"lon"};
730    my $s_minlat = $from->{"lat"};
731    $s_minlon = $to->{"lon"} if ($to->{"lon"} < $s_minlon);
732    $s_minlat = $to->{"lat"} if ($to->{"lat"} < $s_minlat);
733    my $s_maxlon = $from->{"lon"};
734    my $s_maxlat = $from->{"lat"};
735    $s_maxlon = $to->{"lon"} if ($to->{"lon"} > $s_maxlon);
736    $s_maxlat = $to->{"lat"} if ($to->{"lat"} > $s_maxlat);
737
738    printf "BBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
739        $minlat, $minlon, $maxlat, $maxlon if ($debug);
740
741    printf "SBBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n",
742        $s_minlat, $s_minlon, $s_maxlat, $s_maxlon if ($debug);
743
744    # only if the segment is not horizontal
745    if ($latd != 0)
746    {
747        # intersection with top of bounding box
748        $point = { 
749            "side" => "1",
750            "lat" => $maxlat, 
751            "lon" => $from->{"lon"} + ($maxlat - $from->{"lat"}) * $lond / $latd 
752        };
753        push (@$result, $point) 
754            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
755               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
756               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
757
758
759        # intersection with bottom of bounding box
760        $point = { 
761            "side" => "3",
762            "lat" => $minlat, 
763            "lon" => $from->{"lon"} + ($minlat - $from->{"lat"}) * $lond / $latd 
764        };
765        push (@$result, $point) 
766            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
767               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
768               ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon);
769    }
770
771    # only if the segment is not vertical
772    if ($lond != 0)
773    {
774        # intersection with left of bounding box
775        $point = { 
776            "side" => "2",
777            "lat" => $from->{"lat"} + $latd / $lond * ($minlon - $from->{"lon"}),
778            "lon" => $minlon
779        };
780        push (@$result, $point) 
781            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
782               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
783               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
784
785        # intersection with right of bounding box
786        $point = { 
787            "side" => "0",
788            "lat" => $from->{"lat"} + $latd / $lond * ($maxlon - $from->{"lon"}),
789            "lon" => $maxlon
790        };
791        push (@$result, $point) 
792            if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) &&
793               ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) &&
794               ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat);
795    }
796
797    # if more than 1 result, sort by distance from origin of segment
798    # (strictly speaking this sorts by distance squared but why waste
799    # a sqrt call)
800    if (scalar(@$result) > 1)
801    {
802        my @tmp = sort 
803            { (($a->{"lat"} - $from->{"lat"})**2 + ($a->{"lon"} - $from->{"lon"})**2) <=> 
804              (($b->{"lat"} - $from->{"lat"})**2 + ($b->{"lon"} - $from->{"lon"})**2) } @$result;
805        $result = \@tmp;
806    }
807
808    #printf "intersections for segment %f,%f - %f,%f:\n", $from->{"lat"}, $from->{"lon"}, $to->{"lat"}, $to->{"lon"};
809    #print Dumper($result);
810
811    return $result;
812}
813
814# expects point as lat/lon hash, and polygon as a list of segment ids
815# which will be looked up in global $segments to get from/to node ids
816# which will in turn be looked up in $nodes to get lat/lon
817sub polygon_contains_point
818{
819    my ($seglist, $point) = @_;
820    my $p1 =  $nodes->{$segments->{$seglist->[0]}->{"from"}};
821    my $counter = 0;
822    foreach my $seg(@$seglist)
823    {
824        my $p2 = $nodes->{$segments->{$seg}->{"to"}};
825        if ($point->{"lat"} > $p1->{"lat"} || $point->{"lat"} > $p2->{"lat"})
826        {
827            if ($point->{"lat"} < $p1->{"lat"} || $point->{"lat"} < $p2->{"lat"})
828            {
829                if ($point->{"lon"} < $p1->{"lon"} || $point->{"lon"} < $p2->{"lon"})
830                {
831                    if ($p1->{"lat"} != $p2->{"lat"})
832                    {
833                        my $xint = ($point->{"lat"}-$p1->{"lat"})*($p2->{"lon"}-$p1->{"lon"})/($p2->{"lat"}-$p1->{"lat"})+$p1->{"lon"};
834                        if ($p1->{"lon"} == $p2->{"lon"} || $point->{"lon"} <= $xint)
835                        {
836                            $counter++;
837                        }
838                    }
839                }
840            }
841        }
842        $p1 = $p2;
843    }
844    return ($counter%2);
845}
846
847sub node_is_inside
848{
849    my $point = shift;
850    return 0 if ($point->{"lat"} > $maxlat);
851    return 0 if ($point->{"lat"} < $minlat);
852    return 0 if ($point->{"lon"} > $maxlon);
853    return 0 if ($point->{"lon"} < $minlon);
854    return 1;
855}
856
857
858sub compute_angle_from_bbox_center
859{
860    my ($node) = @_;
861    my $opposite_leg = ($node->{"lat"}-($maxlat-$minlat)/2-$minlat);
862    my $adjacent_leg = ($node->{"lon"}-($maxlon-$minlon)/2-$minlon);
863    my $z = cplx($adjacent_leg, $opposite_leg);
864    return (arg($z) < 0) ? arg($z) + 2*pi : arg($z);
865}
866
Note: See TracBrowser for help on using the repository browser.