source: subversion/applications/rendering/tahNG/development/close-areas.pl @ 26211

Last change on this file since 26211 was 7475, checked in by deelkar, 12 years ago

be more verbose in close-areas, improve cutting process

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