| 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 | |
|---|
| 15 | use lib './lib'; |
|---|
| 16 | use Math::Trig; |
|---|
| 17 | use Math::Complex; |
|---|
| 18 | use tahproject; |
|---|
| 19 | use strict; |
|---|
| 20 | |
|---|
| 21 | my $nodecount = 0; |
|---|
| 22 | my $waycount = 0; |
|---|
| 23 | my $relcount = 0; |
|---|
| 24 | my $halfpi = pi()/2; |
|---|
| 25 | my $twopi = pi()*2; |
|---|
| 26 | |
|---|
| 27 | my $LimitY = pi(); |
|---|
| 28 | my $LimitY2 = -pi(); |
|---|
| 29 | my $RangeY = $twopi; |
|---|
| 30 | my $debug=0; |
|---|
| 31 | my $make_open_ways=0; |
|---|
| 32 | my $segments; |
|---|
| 33 | my $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. |
|---|
| 41 | my $direction="cw"; |
|---|
| 42 | |
|---|
| 43 | my @coastline_segments; |
|---|
| 44 | |
|---|
| 45 | my $tilex; |
|---|
| 46 | my $tiley; |
|---|
| 47 | my $zoom; |
|---|
| 48 | my $minlat; |
|---|
| 49 | my $minlon; |
|---|
| 50 | my $maxlat; |
|---|
| 51 | my $maxlon; |
|---|
| 52 | my $border_crossed; |
|---|
| 53 | |
|---|
| 54 | if (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 | } |
|---|
| 61 | elsif (scalar(@ARGV) == 4) |
|---|
| 62 | { |
|---|
| 63 | ($minlat, $minlon, $maxlat, $maxlon) = @ARGV; |
|---|
| 64 | } |
|---|
| 65 | else |
|---|
| 66 | { |
|---|
| 67 | die "need either 'tilex tiley [zoom]' (defaults to zoom=12) or 'minlat minlon maxlat maxlon' on cmdline"; |
|---|
| 68 | } |
|---|
| 69 | |
|---|
| 70 | my $sigma = 0.01; |
|---|
| 71 | |
|---|
| 72 | my $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 | |
|---|
| 78 | my $copy = 1; |
|---|
| 79 | my $waybuf; |
|---|
| 80 | my $is_coastline; |
|---|
| 81 | my @seglist; |
|---|
| 82 | my $last_node_ref; |
|---|
| 83 | my $segcount = 0; |
|---|
| 84 | |
|---|
| 85 | while(<STDIN>) |
|---|
| 86 | { |
|---|
| 87 | while(/(<[^'"<>]*((["'])[^\3]*?\3[^<>"']*)*>)/og) |
|---|
| 88 | { |
|---|
| 89 | my $xmltag = $1; |
|---|
| 90 | if ($xmltag =~ /^\s*<node.*\sid=["'](\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 | |
|---|
| 164 | foreach 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 | |
|---|
| 231 | if (!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 |
|---|
| 259 | my $available_segs; |
|---|
| 260 | grep { $available_segs->{$_} = 1 } keys(%$segments); |
|---|
| 261 | |
|---|
| 262 | while (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 | |
|---|
| 274 | STRING: |
|---|
| 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 | last if (check_segments_intersect($helpernodes->{$helper}->[1], $helpernodes->{$helper}->[0], $nodes->{$currentnode}->{"lon"}, $nodes->{$currentnode}->{"lat"}, $segments) ); # new |
|---|
| 473 | my $newnode = make_node($helpernodes->{$helper}->[0], $helpernodes->{$helper}->[1]); |
|---|
| 474 | my $newseg = make_seg($currentnode, $newnode); |
|---|
| 475 | $currentnode = $newnode; |
|---|
| 476 | push(@seglist, $newseg); |
|---|
| 477 | } |
|---|
| 478 | } |
|---|
| 479 | |
|---|
| 480 | my $newseg = make_seg($currentnode, $segments->{$found}->{"from"}); |
|---|
| 481 | push(@seglist, $newseg); |
|---|
| 482 | |
|---|
| 483 | # if the segment we have found is our start segment (which we added |
|---|
| 484 | # as a special case earlier!), we have a closed way. |
|---|
| 485 | if ($found == $seglist[0]) |
|---|
| 486 | { |
|---|
| 487 | printf("CLOSED\n") if ($debug); |
|---|
| 488 | push(@coastline_segments, \@seglist); |
|---|
| 489 | next; |
|---|
| 490 | } |
|---|
| 491 | |
|---|
| 492 | # else we just plod on. |
|---|
| 493 | push (@seglist, $found); |
|---|
| 494 | $currentseg = $found; |
|---|
| 495 | $currentnode = $segments->{$found}->{"to"}; |
|---|
| 496 | delete $available_segs->{$found}; |
|---|
| 497 | #nprintf "REMAINING: %d = %s\n", scalar(keys(%$available_segs)), join(",", keys(%$available_segs)) if ($debug); |
|---|
| 498 | goto STRING; |
|---|
| 499 | } |
|---|
| 500 | |
|---|
| 501 | ENDOSM: |
|---|
| 502 | |
|---|
| 503 | # if we had no coastline intersection with the bounding box, but coastline |
|---|
| 504 | # was present, we have an island situation and need to add a blue background. |
|---|
| 505 | # FIXME what if we have a "little lake" situation? |
|---|
| 506 | |
|---|
| 507 | unless( $border_crossed ) |
|---|
| 508 | { |
|---|
| 509 | my $state = lookup_handler($helpernodes, $tilex, $tiley, $zoom); |
|---|
| 510 | if( $state eq "10" ) |
|---|
| 511 | { |
|---|
| 512 | # sea |
|---|
| 513 | addBlueRectangle($helpernodes); |
|---|
| 514 | } |
|---|
| 515 | elsif ( $state eq "01" ) |
|---|
| 516 | { |
|---|
| 517 | #land |
|---|
| 518 | } |
|---|
| 519 | else # state 00 or 11 (unknown or mixed) |
|---|
| 520 | { |
|---|
| 521 | my %temp = ("00"=>0, "10"=>0, "01"=>0, "11"=>0);; |
|---|
| 522 | $temp{lookup_handler($helpernodes, $tilex-1, $tiley, $zoom)}++; |
|---|
| 523 | $temp{lookup_handler($helpernodes, $tilex+1, $tiley, $zoom)}++; |
|---|
| 524 | $temp{lookup_handler($helpernodes, $tilex, $tiley-1, $zoom)}++; |
|---|
| 525 | $temp{lookup_handler($helpernodes, $tilex, $tiley+1, $zoom)}++; |
|---|
| 526 | |
|---|
| 527 | if( $temp{"10"} + $temp{"11"} > $temp{"01"} ) # if more sea/coast than land around. |
|---|
| 528 | { |
|---|
| 529 | addBlueRectangle($helpernodes); |
|---|
| 530 | } |
|---|
| 531 | elsif ( ($state eq "11") and ( $temp{"01"} == 0 ) ) |
|---|
| 532 | # if the tile is marked coast but no land near, assume it's a group of islands instead of lakes. |
|---|
| 533 | { |
|---|
| 534 | # coast |
|---|
| 535 | addBlueRectangle($helpernodes); |
|---|
| 536 | } |
|---|
| 537 | else |
|---|
| 538 | { |
|---|
| 539 | #land |
|---|
| 540 | } |
|---|
| 541 | } |
|---|
| 542 | } |
|---|
| 543 | |
|---|
| 544 | if (scalar(@coastline_segments) == 1) |
|---|
| 545 | { |
|---|
| 546 | make_way($coastline_segments[0]); |
|---|
| 547 | } |
|---|
| 548 | else |
|---|
| 549 | { |
|---|
| 550 | make_multipolygon(\@coastline_segments); |
|---|
| 551 | } |
|---|
| 552 | print "</osm>\n"; |
|---|
| 553 | |
|---|
| 554 | |
|---|
| 555 | |
|---|
| 556 | |
|---|
| 557 | # checks segment from ($x0,$y0) to ($x1,$y1) for intersection with all segments |
|---|
| 558 | # in $s |
|---|
| 559 | sub check_segments_intersect |
|---|
| 560 | { |
|---|
| 561 | my ($x0, $y0, $x1, $y1, $s)= @_; |
|---|
| 562 | |
|---|
| 563 | printf("check_segments_intersect: (%f,%f) (%f,%f) %s\n", $x0, $y0, $x1, $y1, join(", ",keys(%$s)) ) if ($debug); |
|---|
| 564 | foreach my $seg(keys(%$s)) |
|---|
| 565 | { |
|---|
| 566 | my $from= $s->{$seg}->{"from"}; |
|---|
| 567 | my $to = $s->{$seg}->{"to"}; |
|---|
| 568 | my $x2= $nodes->{$from}->{"lon"}; |
|---|
| 569 | my $y2= $nodes->{$from}->{"lat"}; |
|---|
| 570 | my $x3= $nodes->{$to} ->{"lon"}; |
|---|
| 571 | my $y3= $nodes->{$to} ->{"lat"}; |
|---|
| 572 | my ($ret)= intersect($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3); |
|---|
| 573 | if ($ret>0) # found intersection |
|---|
| 574 | { |
|---|
| 575 | printf(" SEG %3d: %d (%f,%f) --> %d (%f,%f) \n", $seg, $from, $x2, $y2, $to, $x3, $y3 ) if ($debug); |
|---|
| 576 | return 1; |
|---|
| 577 | } |
|---|
| 578 | } |
|---|
| 579 | return 0; # no intersection |
|---|
| 580 | } |
|---|
| 581 | |
|---|
| 582 | |
|---|
| 583 | # Perl implementation of intersect algorithm for segments |
|---|
| 584 | # first segments is from ($x0,$y0) to ($x1,$y1), second segment is from ($x2,$y2) to ($x3,$y3) |
|---|
| 585 | # source: http://softsurfer.com/Archive/algorithm_0104/algorithm_0104B.htm#Line%20Intersections |
|---|
| 586 | # Modification: This implementation accepts only _true_ intersections, i.e. if the two segments |
|---|
| 587 | # only "touch" at the ends it's _not_ detected as an intersection _intentionally_ |
|---|
| 588 | sub intersect |
|---|
| 589 | { |
|---|
| 590 | my ($x0, $y0, $x1, $y1, $x2, $y2, $x3, $y3)= @_; |
|---|
| 591 | |
|---|
| 592 | my $cpproduct= ($x1-$x0)*($y3-$y2) - ($y1-$y0)*($x3-$x2); |
|---|
| 593 | if (abs($cpproduct)>0.000000000001) |
|---|
| 594 | { |
|---|
| 595 | my $t1= (($x1-$x0)*($y0-$y2) - ($y1-$y0)*($x0-$x2)) / $cpproduct; |
|---|
| 596 | return 0 if ( ($t1<=0) || ($t1>=1) ); |
|---|
| 597 | my $t2= (($x3-$x2)*($y0-$y2) - ($y3-$y2)*($x0-$x2)) / $cpproduct; |
|---|
| 598 | return 0 if ( ($t2<=0) || ($t2>=1) ); |
|---|
| 599 | # a true intersection |
|---|
| 600 | my $xs= $x0 + $t2*($x1-$x0); |
|---|
| 601 | my $ys= $y0 + $t2*($y1-$y0); |
|---|
| 602 | return (1, $xs, $ys); |
|---|
| 603 | } |
|---|
| 604 | else |
|---|
| 605 | { |
|---|
| 606 | # are the two segments collinear ? |
|---|
| 607 | my $cpp1= ($x1-$x0)*($y0-$y2) -($y1-$y0)*($x0-$x2); |
|---|
| 608 | my $cpp2= ($x3-$x2)*($y0-$y2) -($y3-$y2)*($x0-$x2); |
|---|
| 609 | return 0 if ( ($cpp1!=0) || ($cpp2!=0) ); # they are _not_ collinear |
|---|
| 610 | |
|---|
| 611 | return -1 if ( ($x0==$x1) && ($y0==$y1) ); # error: first segment is only a single point |
|---|
| 612 | return -1 if ( ($x2==$x3) && ($y2==$y3) ); # error: second segment is only a single point |
|---|
| 613 | |
|---|
| 614 | my ($t0, $t1); |
|---|
| 615 | if (($x3-$x2)!=0) |
|---|
| 616 | { |
|---|
| 617 | $t0= ($x0-$x2) / ($x3-$x2); |
|---|
| 618 | $t1= ($x1-$x2) / ($x3-$x2); |
|---|
| 619 | } |
|---|
| 620 | else |
|---|
| 621 | { |
|---|
| 622 | $t0= ($y0-$y2) / ($y3-$y2); |
|---|
| 623 | $t1= ($y1-$y2) / ($y3-$y2); |
|---|
| 624 | } |
|---|
| 625 | if ($t0>$t1) |
|---|
| 626 | { |
|---|
| 627 | my $tmp= $t0; $t0= $t1; $t1= $tmp; |
|---|
| 628 | } |
|---|
| 629 | return 0 if ( ($t0>1) || ($t1<0) ); |
|---|
| 630 | $t0 = ($t0<0)? 0 : $t0; # clip to min 0 |
|---|
| 631 | $t1 = ($t1>1)? 1 : $t1; # clip to max 1 |
|---|
| 632 | if ($t0==$t1) |
|---|
| 633 | { |
|---|
| 634 | return (1, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2) ); |
|---|
| 635 | } |
|---|
| 636 | return (2, $x2 + $t0*($x3-$x2), $y2 + $t0*($y3-$y2), $x2 + $t1*($x3-$x2), $y2 + $t1*($y3-$y2)); |
|---|
| 637 | } |
|---|
| 638 | } |
|---|
| 639 | |
|---|
| 640 | |
|---|
| 641 | |
|---|
| 642 | sub make_node |
|---|
| 643 | { |
|---|
| 644 | my ($lat, $lon) = @_; |
|---|
| 645 | my $id = --$nodecount; |
|---|
| 646 | print "<node id='$id' lat='$lat' lon='$lon' />\n"; |
|---|
| 647 | # save node in array for later processing |
|---|
| 648 | $nodes->{$id}={"lat"=>$lat, "lon"=>$lon}; |
|---|
| 649 | return $id; |
|---|
| 650 | } |
|---|
| 651 | |
|---|
| 652 | sub make_seg |
|---|
| 653 | { |
|---|
| 654 | my ($from, $to) = @_; |
|---|
| 655 | my $id = ++$segcount; |
|---|
| 656 | $segments->{$id} = { "from" => $from, "to" => $to }; |
|---|
| 657 | return $id; |
|---|
| 658 | } |
|---|
| 659 | |
|---|
| 660 | sub make_way |
|---|
| 661 | { |
|---|
| 662 | my ($seglist, $open) = @_; |
|---|
| 663 | my $id = --$waycount; |
|---|
| 664 | print "<way id='$id'>\n"; |
|---|
| 665 | print " <tag k='natural' v='coastline' />\n"; |
|---|
| 666 | print " <tag k='created-with' v='close-areas.pl' />\n"; |
|---|
| 667 | print " <tag k='close-areas.pl:debug' v='open way' />\n" if ($open); |
|---|
| 668 | |
|---|
| 669 | my $first = 1; |
|---|
| 670 | foreach my $seg(@$seglist) |
|---|
| 671 | { |
|---|
| 672 | printf " <nd ref='%s' />\n", $segments->{$seg}->{"from"} if ($first); |
|---|
| 673 | printf " <nd ref='%s' />\n", $segments->{$seg}->{"to"}; |
|---|
| 674 | $first = 0; |
|---|
| 675 | } |
|---|
| 676 | |
|---|
| 677 | print "</way>\n"; |
|---|
| 678 | return $id; |
|---|
| 679 | } |
|---|
| 680 | |
|---|
| 681 | sub make_multipolygon |
|---|
| 682 | { |
|---|
| 683 | my ($seglistlist) = @_; |
|---|
| 684 | return unless scalar(@$seglistlist); |
|---|
| 685 | |
|---|
| 686 | # we will create a list of multipolygons. each multipolygon is a list of |
|---|
| 687 | # ways, with the first being the "outer" polygon. |
|---|
| 688 | my $multipolygons = []; |
|---|
| 689 | |
|---|
| 690 | # first create all the ways. |
|---|
| 691 | my $way_ids = []; |
|---|
| 692 | for (my $i=0; $i<scalar(@$seglistlist); $i++) |
|---|
| 693 | { |
|---|
| 694 | $way_ids->[$i] = make_way($seglistlist->[$i]); |
|---|
| 695 | } |
|---|
| 696 | |
|---|
| 697 | # now find out which way is contained in which. since our polygons do not |
|---|
| 698 | # intersect, we simply take one node of a way and check whether this node |
|---|
| 699 | # is contained in any of the other ways; if yes, the whole way is contained. |
|---|
| 700 | my $contained_in = []; |
|---|
| 701 | my $children = []; |
|---|
| 702 | for (my $i=0; $i<scalar(@$seglistlist); $i++) |
|---|
| 703 | { |
|---|
| 704 | my $pt = $nodes->{$segments->{$seglistlist->[$i]->[0]}->{"from"}}; |
|---|
| 705 | for (my $j = 0; $j < scalar(@$seglistlist); $j++) # fixme really need to check all? |
|---|
| 706 | { |
|---|
| 707 | next if ($j==$i); |
|---|
| 708 | if (polygon_contains_point($seglistlist->[$j], $pt)) |
|---|
| 709 | { |
|---|
| 710 | $contained_in->[$i] = $j; |
|---|
| 711 | $children->[$j]++; |
|---|
| 712 | last; |
|---|
| 713 | } |
|---|
| 714 | } |
|---|
| 715 | } |
|---|
| 716 | |
|---|
| 717 | # now write multipolygons |
|---|
| 718 | for (my $i=0; $i<scalar(@$seglistlist); $i++) |
|---|
| 719 | { |
|---|
| 720 | # ways that don't have children do not trigger a relation |
|---|
| 721 | next unless $children->[$i]; |
|---|
| 722 | my $id = --$relcount; |
|---|
| 723 | print "<relation id='$id'>\n"; |
|---|
| 724 | print " <tag k='type' v='multipolygon' />\n"; |
|---|
| 725 | printf " <member type='way' ref='%s' role='outer' />\n", $way_ids->[$i]; |
|---|
| 726 | for (my $j = 0; $j < scalar(@$seglistlist); $j++) |
|---|
| 727 | { |
|---|
| 728 | if (defined($contained_in->[$j]) && $contained_in->[$j] == $i) |
|---|
| 729 | { |
|---|
| 730 | printf " <member type='way' ref='%s' role='inner' />\n", $way_ids->[$j]; |
|---|
| 731 | } |
|---|
| 732 | } |
|---|
| 733 | print "</relation>\n"; |
|---|
| 734 | } |
|---|
| 735 | } |
|---|
| 736 | |
|---|
| 737 | # index lookup by Martijn van Oosterhout |
|---|
| 738 | sub lookup_handler |
|---|
| 739 | { |
|---|
| 740 | my ($helpernodes, $x, $y, $zoom) = @_; |
|---|
| 741 | # make it use z12 x,y coordinates |
|---|
| 742 | # this looks up the most upper left z12 tile in zoom<12. This probably |
|---|
| 743 | # needs to be made smarter for islands etc. |
|---|
| 744 | my $tilex = $x*(2**(12-$zoom)); |
|---|
| 745 | my $tiley = $y*(2**(12-$zoom)); |
|---|
| 746 | my $tileoffset = ($tiley * (2**12)) + $tilex; |
|---|
| 747 | my $fh; |
|---|
| 748 | open($fh, "<", "png2tileinfo/oceantiles_12.dat") or die; |
|---|
| 749 | seek $fh, int($tileoffset / 4), 0; |
|---|
| 750 | my $buffer; |
|---|
| 751 | read $fh, $buffer, 1; |
|---|
| 752 | $buffer = substr( $buffer."\0", 0, 1 ); |
|---|
| 753 | $buffer = unpack "B*", $buffer; |
|---|
| 754 | my $str = substr( $buffer, 2*($tileoffset % 4), 2 ); |
|---|
| 755 | close($fh); |
|---|
| 756 | |
|---|
| 757 | print("lookup handler finds: $str\n") if $debug; |
|---|
| 758 | |
|---|
| 759 | # $str eq "00" => unknown (not yet checked) |
|---|
| 760 | # $str eq "01" => known land |
|---|
| 761 | # $str eq "10" => known sea |
|---|
| 762 | # $str eq "11" => known edge tile |
|---|
| 763 | |
|---|
| 764 | return $str; |
|---|
| 765 | } |
|---|
| 766 | |
|---|
| 767 | sub addBlueRectangle |
|---|
| 768 | { |
|---|
| 769 | my $helpernodes = shift; |
|---|
| 770 | my @n; |
|---|
| 771 | my @s; |
|---|
| 772 | for (my $i=0; $i<4; $i++) |
|---|
| 773 | { |
|---|
| 774 | $n[$i] = make_node($helpernodes->{$i}->[0], |
|---|
| 775 | $helpernodes->{$i}->[1]); |
|---|
| 776 | } |
|---|
| 777 | for (my $i=0; $i<4; $i++) |
|---|
| 778 | { |
|---|
| 779 | if ($direction eq "ccw") |
|---|
| 780 | { |
|---|
| 781 | $s[$i] = make_seg($n[$i], $n[($i+1)%4]); |
|---|
| 782 | } |
|---|
| 783 | else |
|---|
| 784 | { |
|---|
| 785 | $s[3-$i] = make_seg($n[($i+3)%4], $n[($i+2)%4]); |
|---|
| 786 | } |
|---|
| 787 | } |
|---|
| 788 | push(@coastline_segments, \@s); |
|---|
| 789 | } |
|---|
| 790 | |
|---|
| 791 | # this takes two points (hashes with lat/lon keys) as arguments and returns |
|---|
| 792 | # an array reference to an array containing up to four points (hashes with |
|---|
| 793 | # lat/lon keys and an added "side" key, where 0=right 1=top 2=left 3=bottom) |
|---|
| 794 | # denoting the intersections of the line formed by the input points with the |
|---|
| 795 | # bounding box. |
|---|
| 796 | # |
|---|
| 797 | # 0 results - segment is fully inside or fully outside bounding box |
|---|
| 798 | # 1 result - segment begins inside, ends outside or vice versa |
|---|
| 799 | # 2 results - segment begins outside and ends outside, but cuts through |
|---|
| 800 | # bounding box |
|---|
| 801 | # 3 results - can't think how this can happen |
|---|
| 802 | # 4 results - as case 2 but segment is exactly at the bounding box diagonal, |
|---|
| 803 | # thus intersecting with all four edges (freak case) |
|---|
| 804 | # |
|---|
| 805 | # If there is more than one result, results are sorted by distance from the |
|---|
| 806 | # first point. |
|---|
| 807 | sub compute_bbox_intersections |
|---|
| 808 | { |
|---|
| 809 | my ($from, $to) = @_; |
|---|
| 810 | my $result = []; |
|---|
| 811 | my $point; |
|---|
| 812 | |
|---|
| 813 | # div by zero FIXME |
|---|
| 814 | my $latd = ($to->{"lat"} - $from->{"lat"}); |
|---|
| 815 | my $lond = ($to->{"lon"} - $from->{"lon"}); |
|---|
| 816 | |
|---|
| 817 | # segment's bbox |
|---|
| 818 | my $s_minlon = $from->{"lon"}; |
|---|
| 819 | my $s_minlat = $from->{"lat"}; |
|---|
| 820 | $s_minlon = $to->{"lon"} if ($to->{"lon"} < $s_minlon); |
|---|
| 821 | $s_minlat = $to->{"lat"} if ($to->{"lat"} < $s_minlat); |
|---|
| 822 | my $s_maxlon = $from->{"lon"}; |
|---|
| 823 | my $s_maxlat = $from->{"lat"}; |
|---|
| 824 | $s_maxlon = $to->{"lon"} if ($to->{"lon"} > $s_maxlon); |
|---|
| 825 | $s_maxlat = $to->{"lat"} if ($to->{"lat"} > $s_maxlat); |
|---|
| 826 | |
|---|
| 827 | printf "BBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n", |
|---|
| 828 | $minlat, $minlon, $maxlat, $maxlon if ($debug); |
|---|
| 829 | |
|---|
| 830 | printf "SBBOX:\n minlat %f\n minlon %f\n maxlat %f\n maxlon %f\n", |
|---|
| 831 | $s_minlat, $s_minlon, $s_maxlat, $s_maxlon if ($debug); |
|---|
| 832 | |
|---|
| 833 | # only if the segment is not horizontal |
|---|
| 834 | if ($latd != 0) |
|---|
| 835 | { |
|---|
| 836 | # intersection with top of bounding box |
|---|
| 837 | $point = { |
|---|
| 838 | "side" => "1", |
|---|
| 839 | "lat" => $maxlat, |
|---|
| 840 | "lon" => $from->{"lon"} + ($maxlat - $from->{"lat"}) * $lond / $latd |
|---|
| 841 | }; |
|---|
| 842 | push (@$result, $point) |
|---|
| 843 | if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) && |
|---|
| 844 | ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) && |
|---|
| 845 | ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon); |
|---|
| 846 | |
|---|
| 847 | |
|---|
| 848 | # intersection with bottom of bounding box |
|---|
| 849 | $point = { |
|---|
| 850 | "side" => "3", |
|---|
| 851 | "lat" => $minlat, |
|---|
| 852 | "lon" => $from->{"lon"} + ($minlat - $from->{"lat"}) * $lond / $latd |
|---|
| 853 | }; |
|---|
| 854 | push (@$result, $point) |
|---|
| 855 | if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) && |
|---|
| 856 | ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) && |
|---|
| 857 | ($point->{"lon"} >= $minlon && $point->{"lon"} <= $maxlon); |
|---|
| 858 | } |
|---|
| 859 | |
|---|
| 860 | # only if the segment is not vertical |
|---|
| 861 | if ($lond != 0) |
|---|
| 862 | { |
|---|
| 863 | # intersection with left of bounding box |
|---|
| 864 | $point = { |
|---|
| 865 | "side" => "2", |
|---|
| 866 | "lat" => $from->{"lat"} + $latd / $lond * ($minlon - $from->{"lon"}), |
|---|
| 867 | "lon" => $minlon |
|---|
| 868 | }; |
|---|
| 869 | push (@$result, $point) |
|---|
| 870 | if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) && |
|---|
| 871 | ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) && |
|---|
| 872 | ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat); |
|---|
| 873 | |
|---|
| 874 | # intersection with right of bounding box |
|---|
| 875 | $point = { |
|---|
| 876 | "side" => "0", |
|---|
| 877 | "lat" => $from->{"lat"} + $latd / $lond * ($maxlon - $from->{"lon"}), |
|---|
| 878 | "lon" => $maxlon |
|---|
| 879 | }; |
|---|
| 880 | push (@$result, $point) |
|---|
| 881 | if ($point->{"lat"} >= $s_minlat && $point->{"lat"} <= $s_maxlat) && |
|---|
| 882 | ($point->{"lon"} >= $s_minlon && $point->{"lon"} <= $s_maxlon) && |
|---|
| 883 | ($point->{"lat"} >= $minlat && $point->{"lat"} <= $maxlat); |
|---|
| 884 | } |
|---|
| 885 | |
|---|
| 886 | # if more than 1 result, sort by distance from origin of segment |
|---|
| 887 | # (strictly speaking this sorts by distance squared but why waste |
|---|
| 888 | # a sqrt call) |
|---|
| 889 | if (scalar(@$result) > 1) |
|---|
| 890 | { |
|---|
| 891 | my @tmp = sort |
|---|
| 892 | { (($a->{"lat"} - $from->{"lat"})**2 + ($a->{"lon"} - $from->{"lon"})**2) <=> |
|---|
| 893 | (($b->{"lat"} - $from->{"lat"})**2 + ($b->{"lon"} - $from->{"lon"})**2) } @$result; |
|---|
| 894 | $result = \@tmp; |
|---|
| 895 | } |
|---|
| 896 | |
|---|
| 897 | #printf "intersections for segment %f,%f - %f,%f:\n", $from->{"lat"}, $from->{"lon"}, $to->{"lat"}, $to->{"lon"}; |
|---|
| 898 | #print Dumper($result); |
|---|
| 899 | |
|---|
| 900 | return $result; |
|---|
| 901 | } |
|---|
| 902 | |
|---|
| 903 | # expects point as lat/lon hash, and polygon as a list of segment ids |
|---|
| 904 | # which will be looked up in global $segments to get from/to node ids |
|---|
| 905 | # which will in turn be looked up in $nodes to get lat/lon |
|---|
| 906 | sub polygon_contains_point |
|---|
| 907 | { |
|---|
| 908 | my ($seglist, $point) = @_; |
|---|
| 909 | my $p1 = $nodes->{$segments->{$seglist->[0]}->{"from"}}; |
|---|
| 910 | my $counter = 0; |
|---|
| 911 | foreach my $seg(@$seglist) |
|---|
| 912 | { |
|---|
| 913 | my $p2 = $nodes->{$segments->{$seg}->{"to"}}; |
|---|
| 914 | if ($point->{"lat"} >= $p1->{"lat"} || $point->{"lat"} >= $p2->{"lat"}) |
|---|
| 915 | { |
|---|
| 916 | if ($point->{"lat"} < $p1->{"lat"} || $point->{"lat"} < $p2->{"lat"}) |
|---|
| 917 | { |
|---|
| 918 | if ($point->{"lon"} < $p1->{"lon"} || $point->{"lon"} < $p2->{"lon"}) |
|---|
| 919 | { |
|---|
| 920 | if ($p1->{"lat"} != $p2->{"lat"}) |
|---|
| 921 | { |
|---|
| 922 | my $xint = ($point->{"lat"}-$p1->{"lat"})*($p2->{"lon"}-$p1->{"lon"})/($p2->{"lat"}-$p1->{"lat"})+$p1->{"lon"}; |
|---|
| 923 | if ($p1->{"lon"} == $p2->{"lon"} || $point->{"lon"} <= $xint) |
|---|
| 924 | { |
|---|
| 925 | $counter++; |
|---|
| 926 | } |
|---|
| 927 | } |
|---|
| 928 | } |
|---|
| 929 | } |
|---|
| 930 | } |
|---|
| 931 | $p1 = $p2; |
|---|
| 932 | } |
|---|
| 933 | return ($counter%2); |
|---|
| 934 | } |
|---|
| 935 | |
|---|
| 936 | sub node_is_inside |
|---|
| 937 | { |
|---|
| 938 | my $point = shift; |
|---|
| 939 | return 0 if ($point->{"lat"} > $maxlat); |
|---|
| 940 | return 0 if ($point->{"lat"} < $minlat); |
|---|
| 941 | return 0 if ($point->{"lon"} > $maxlon); |
|---|
| 942 | return 0 if ($point->{"lon"} < $minlon); |
|---|
| 943 | return 1; |
|---|
| 944 | } |
|---|
| 945 | |
|---|
| 946 | |
|---|
| 947 | sub compute_angle_from_bbox_center |
|---|
| 948 | { |
|---|
| 949 | my ($node) = @_; |
|---|
| 950 | my $opposite_leg = ($node->{"lat"}-($maxlat-$minlat)/2-$minlat); |
|---|
| 951 | my $adjacent_leg = ($node->{"lon"}-($maxlon-$minlon)/2-$minlon); |
|---|
| 952 | my $z = cplx($adjacent_leg, $opposite_leg); |
|---|
| 953 | return (arg($z) < 0) ? arg($z) + 2*pi : arg($z); |
|---|
| 954 | } |
|---|
| 955 | |
|---|
