File Coverage

blib/lib/Geo/Index.pm
Criterion Covered Total %
statement 810 1012 80.0
branch 342 534 64.0
condition 74 146 50.6
subroutine 41 54 75.9
pod 23 37 62.1
total 1290 1783 72.3


line stmt bran cond sub pod time code
1             package Geo::Index;
2              
3              
4             # Geo::Index
5             # Copyright 2019 Alexander Hajnal, All rights reserved
6             #
7             # Alex Kent Hajnal
8             # --------------------------------------
9             # akh@cpan.org
10             # https://alephnull.net/software
11             # https://github.com/Alex-Kent/Geo-Index
12             #
13             # This module is free software; you can redistribute it and/or modify it under
14             # the same terms as Perl itself. See the LICENSE file or perlartistic(1).
15              
16              
17             require 5.00405;
18              
19 11     11   81318 use warnings;
  11         23  
  11         361  
20 11     11   55 use strict;
  11         54  
  11         278  
21              
22 11     11   64 use Carp;
  11         21  
  11         666  
23 11     11   5613 use Math::Trig;
  11         149073  
  11         1668  
24 11     11   6245 use POSIX qw( ceil );
  11         68412  
  11         67  
25 11     11   16079 use Config;
  11         24  
  11         703  
26              
27             #. Note: Comments starting with #. #: and #> are used
28             #. by the author to trigger syntax highlighting rules.
29              
30              
31             #. Numeric keys are smaller and faster but require a Perl with 64-bit integer support
32             #. When debugging this module it can be useful to manually set this to 0
33 11 50   11   78 use constant USE_NUMERIC_KEYS => ( $Config{use64bitint} ) ? 1 : 0;
  11         20  
  11         1826  
34              
35             #use constant USE_NUMERIC_KEYS => 0; # Uncomment to force text keys
36             #use constant USE_NUMERIC_KEYS => 1; # Uncomment to force numeric keys
37 11     11   81 use constant USE_PACKED_KEYS => 0; # Change to 1 to pack numeric keys
  11         27  
  11         763  
38              
39              
40             #. Text keys
41             #. ------------------------------------------------------------------------------
42             #. Text keys have the format "level:lat_idx,lon_idx".
43             #.
44             #. Level is the level number, lat_idx and lon_idx are integer latitude and
45             #. longitude values scaled to the key's level. The ALL level, latitude, or
46             #. longitude is represented by the string "ALL".
47             #.
48             #. Indices using text keys require roughly twice as much memory as indices using
49             #. numeric keys. In addition, text keys are about 30% slower than numeric ones.
50              
51             #. Numeric keys
52             #. ------------------------------------------------------------------------------
53             #. 64-bit numeric keys are broken down into three parts:
54             #. Level -> bits 63..58
55             #. Latitude -> bits 29..57
56             #. Longitude -> bits 0..28
57             #.
58             #. The values of latitude and longitude in numeric keys are integer values
59             #. scaled to the key's level. The ALL level, latitude, or longitude is
60             #. represented by all 1 bits in the respective bitfield.
61              
62             #. Packed keys
63             #. ------------------------------------------------------------------------------
64             #. Packed numeric keys are numeric keys run through pack("Q", $key)
65             #.
66             #. There appears to be no performance benefit from using them.
67              
68              
69              
70             #. Value for ALL for Level
71             #. Can also be used for masking values
72 11     11   67 use constant MASK_LEVEL => ( 1 << 6 ) - 1;
  11         21  
  11         613  
73              
74             #. Value for ALL for Latitude or Longitude
75             #. Can also be used for masking values
76 11     11   86 use constant MASK_LATLON => ( 1 << 29 ) - 1;
  11         21  
  11         574  
77              
78             #. Used in key specifications for global [ ALL, ALL, ALL ] and polar areas [ ..., ..., ALL ]
79 11     11   82 use constant ALL => -1;
  11         23  
  11         578  
80              
81              
82              
83              
84              
85             =encoding utf8
86              
87             =head1 NAME
88              
89             Geo::Index - Geographic indexer
90              
91             =cut
92              
93 11     11   75 use vars qw ($VERSION);
  11         21  
  11         2647  
94             $VERSION = 'v0.0.8';
95              
96             =head1 VERSION
97              
98             This document describes Geo::Index version 0.0.8
99              
100             =cut
101              
102              
103              
104              
105             #. Attempt to load C low-level code library
106              
107              
108             #. The C code is located in Index.xs is the code's root directory
109              
110              
111             #. Boilerplate for compiled code
112             require Exporter;
113             *import = \&Exporter::import;
114             require DynaLoader;
115 11     11 0 2453 sub dl_load_flags {0} # Prevent DynaLoader from complaining and croaking
116              
117             #. Attempt to load the C low-level code library
118             eval { DynaLoader::bootstrap Geo::Index $VERSION; };
119              
120             #. Note whether C low-level code library is available
121             my $C_CODE_COMPILED;
122             if ($@) {
123             $C_CODE_COMPILED = 0;
124             } else {
125             $C_CODE_COMPILED = 1;
126             }
127            
128             #. Choose which C function to export
129             @Geo::Index::EXPORT = qw(); # Symbols to export by default
130             #@Geo::Index::EXPORT_OK = qw(
131             # GetCCodeVersion fast_log2_double fast_log2_float ComputeAreaExtrema_float
132             # ComputeAreaExtrema_double ComputeAreaExtrema_double SetUpDistance_float HaversineDistance_float
133             # SetUpDistance_double SetUpDistance_double HaversineDistance_double
134             # );
135             @Geo::Index::EXPORT_OK = qw(); # Symbols to export on request
136              
137              
138              
139              
140             =head1 SYNOPSIS
141              
142             # Create and populate a geographic index
143            
144             use Geo::Index;
145            
146             @points = (
147             { lat => 1.0, lon => 2.0 },
148             { lat => -90.0, lon => 0.0, name => 'South Pole' },
149             { lat => 30.0, lon => -20.0, ele => 123.4 }
150             );
151             $point = { lat=>10.0, lon=>20.0 };
152            
153             $index = Geo::Index->new();
154             $index->IndexPoints( \@points );
155             $index->Index( $point );
156             $index->Index( [ 30, 40 ] );
157            
158            
159             # Search index
160            
161             %search_options = ( sort_results => 1, radius=>5_000_000 );
162             $results = $index->Search( [ -80, 20 ], \%search_options );
163             print "$$results[0]{name}\n"; # Prints 'South Pole'
164            
165             # Get all points in the southern hemisphere
166             $results = $index->SearchByBounds( [ -180, -90, 180, 0 ] );
167             print "$$results[0]{name}\n"; # Also prints 'South Pole'
168            
169             ($closest) = $index->Closest( [ -80, 20 ] );
170             print "$$closest{name}\n"; # Also prints 'South Pole'
171            
172             ($closest) = $index->Closest( $points[1], { post_condition=>'NONE' } );
173             print "$$closest{name}\n"; # Also prints 'South Pole'
174            
175             ($farthest) = $index->Farthest( [ 90, 0 ] );
176             print "$$farthest{name}\n"; # Also prints 'South Pole'
177            
178             # Compute distance in meters between two points (using haversine formula)
179            
180             $m = $index->Distance( { lat=>51.507222, lon=>-0.1275 }, [ -6.2, 106.816667 ] );
181             printf("London to Jakarta: %i km\n", $m / 1000);
182            
183             $index->DistanceFrom( [ 90, 0 ] );
184             $m = $index->DistanceTo( $points[1] );
185             printf("Pole to pole: %i km\n", $m / 1000);
186            
187             =head1 DESCRIPTION
188              
189             Geo::Index is a Perl module for creating in-memory geographic points indices.
190             Once points have been indexed, fast searches can be run.
191              
192             Efficient searches methods include B>> to get all points
193             within a distance from a given point, B>> to get all
194             points in an given area, B>> to get the closest points to a
195             given point, and B>> to get the farthest points from a given
196             point.
197              
198             Additional methods are provided to compute distances between arbitrary points
199             (EB>>, B>>, and B>>E)
200             and to get the size in meters of one degree or the size in degrees of one meter
201             at a given point (B>> and B>>,
202             respectively).
203              
204             While by default computations are done for the Earth, other bodies can be used
205             by supplying appropriates radii and circumferences to Bnew( ... )>>>.
206              
207             =head1 POINTS
208              
209             Geo::Index works with points on a spherical body. Points are hash references
210             containing, at a minimum, C and C entries which give the point's
211             position in degrees. Additional hash entries can be present and will be both
212             ignored and preserved. The C>, C>,
213             C>, C>, C>,
214             C>, C>, and C>
215             methods add additional entries in point hashes.
216              
217             The hash entries used by Geo::Gpx are shown below. Apart from C and C
218             these values are created by Geo::Gpx. Unless noted, these values may be read but
219             should not be set, altered, or deleted.
220              
221             =over
222              
223             =item *
224             B> - Point's latitude in degrees [ -90 .. 90 ]
225              
226             =item *
227             B> - Point longitude in degrees [ -180 .. 180 )
228              
229             These two values may be changed but the altered point should then be re-indexed
230             using C> before further searches are run.
231              
232             =item *
233             B> - The optional user data supplied when a point was created
234             using the array shorthand. This contents of this field may be freely modified
235             by the user. See C> and C>, below.
236              
237             =item *
238             B> - The point's latitude in radians [ -pi/2 .. pi/2 ]
239              
240             =item *
241             B> - The point's longitude in radians [ -pi .. pi )
242              
243             =item *
244             B> - Circumference (in meters) of the circle of latitude
245             that the point falls on. This is computed from the body's equatorial
246             circumference assuming a spherical (not an oblate) body.
247              
248             =item *
249             B> - Distance (in meters) of point from search
250             point of previous search. The distance computation assumes a spherical body
251             and is computed using a ruggedized version of the haversine formula. This
252             value is only generated when C> is called with the C
253             or C option. See also C>, C>,
254             and C>.
255              
256             =item *
257             B> - Distance (in meters) of point from search
258             point's antipode as determined by a previous call to C>.
259             This distance is computed using a ruggedized version of the haversine formula.
260              
261             =back
262              
263             As a convenience, most methods allow points to be specified using a shorthand
264             notation S, I ]>> or S, I, I ]>>. Points
265             given in this notation will be converted to hash-based points. If a point
266             created using this notation is returned as a search result it will be as a
267             reference to the hash constructed by Geo::Index and not as a reference to the
268             original array. To access the data field of a point created using the shorthand
269             notation use C<$$point{'data'}> where C<$point> is a search result point.
270              
271             Any fields added to the indexed points by Geo::Index can be removed using
272             C> and C>.
273              
274             =head1 METHODS
275              
276             =cut
277              
278       11     BEGIN {
279              
280             } # END BEGIN
281              
282              
283 11     11   5208 use fields qw(index indices positions planetary_radius planetary_diameter polar_circumference equatorial_circumference levels max_level max_size quiet);
  11         17400  
  11         44  
284              
285              
286              
287             =head2 Geo::Index-Enew( ... )
288              
289             =over
290              
291             C<$index = Geo::Index-Enew()>;
292              
293             =over
294              
295             Create a new empty index using default options: radius and circumferences are those of Earth, C is set to 20 (~40Em index resolution).
296              
297             =back
298              
299             C<$index = Geo::Index-Enew( \@points );>
300              
301             =over
302              
303             Create a new index using default options and populate it with the given points.
304              
305             The points in the array can be in either hash or array notation.
306              
307             =back
308              
309             C<$index = Geo::Index-Enew( \%options );>
310              
311             =over
312              
313             Create a new empty index using specified options.
314              
315             =back
316              
317             C<$index = Geo::Index-Enew( \@points, \%options );>
318              
319             =over
320              
321             Create a new index using specified options and populate it with the given points.
322              
323             The points in the array can be in either hash or array notation.
324              
325             =back
326              
327             B
328              
329             When a Geo::Index object is created, one can specify various options to fine-tune
330             its behavior. The default values are suitable for a high-resolution index of Earth.
331              
332             =over
333              
334             B>
335              
336             =over
337              
338             Average planetary radius (in meters). S<(default: 6371230)>
339              
340             If a C is specified but C or C
341             are not given then they will be calculated from the radius ( 2 * pi * radius )
342              
343             =back
344              
345             B>
346              
347             =over
348              
349             Polar (meridional) circumference of the object the points lie on (in meters). S<(default: 40007863)>
350              
351             =back
352              
353             B>
354              
355             =over
356              
357             Circumference at the equator of the object the points lie on (in meters). S<(default: 40075017)>
358              
359             =back
360              
361             B>
362              
363             =over
364              
365             Depth of index. S<(valid: E0, E31; default: 20)>
366              
367             Note that the C parameter specifies the number of non-full-globe index
368             levels to generate and NOT the deepest index level. (Level -1, covering the
369             entire globe, is always generated) For example, setting C to 20
370             generates indices at levels 0 through 19 (plus level -1).
371              
372             A summary of typical tile levels is shown below. To choose a value for the
373             C option using the table add 1 to the 'Level' shown for the desired
374             maximum level of detail. The 'Grid' column shows the north-south size of each
375             tile in meters at a each level. The 'Size' column shows the initial amount of
376             RAM needed for an indexed set of S<1Emillion> random points using numeric
377             keys on a 64-bit system when that level is the most detailed one (sizes may grow
378             moderately once searches are run).
379              
380             Level Grid Size Level Grid Size
381             ----- --------- --------------- ----- ------- -------
382             -1 ~40000 km (entire planet) 12 ~5 km ~2.4 GB
383             0 ~20000 km 13 ~2.5 km ~2.7 GB
384             1 ~10000 km ~1.0 GB 14 ~1.2 km ~3.1 GB
385             2 ~5000 km ~1.0 GB 15 ~600 m ~3.3 GB
386             3 ~2500 km ~1.1 GB 16 ~300 m ~3.6 GB
387             4 ~1250 km ~1.2 GB 17 ~150 m ~3.8 GB
388             5 ~625 km ~1.3 GB 18 ~75 m ~4.1 GB
389             6 ~315 km ~1.4 GB 19 ~40 m ~4.4 GB
390             7 ~155 km ~1.5 GB 20 ~20 m ~4.6 GB
391             8 ~80 km ~1.6 GB 21 ~10 m ~4.9 GB
392             9 ~40 km ~1.7 GB 22 ~5 m ~5.1 GB
393             10 ~20 km ~1.9 GB 23 ~2 m ~5.4 GB
394             11 ~10 km ~2.1 GB 24 ~1 m ~5.6 GB
395              
396             For reference, the memory usage of the array of S<1 million> random, unindexed
397             points is about S<440 MB>, growing to about S<540 MB> with index use (about 100
398             bytes per point); the former amount is included in the index memory usage shown
399             above.
400              
401             =back
402              
403             B>
404              
405             =over
406              
407             Choose the type of low-level functions to use.
408             S<(default: 'C' if available, 'C' otherwise)>
409              
410             Geo::Index will attempt to use compiled C code to speed up certain calculations.
411             If the compilation fails (or was blocked by the user) then equivalent (but
412             slower) Perl code will be used.
413              
414             This option can be used to explicitly request the type of code to use. When set
415             to 'C' then compiled C code using single-precision floating point will
416             be requested. When set to 'C' then compiled C code using double-precision
417             floating point will be requested. When set to 'C' then Perl code will be
418             used. If compiled code is unavailable then Perl code will be used regardless of
419             what was requested.
420              
421             Perl natively uses double-precision floating point. On modern hardware
422             double-precision is slightly faster than single-precision. On certain platforms,
423             however, it may be preferable to use single-precision instead of double-precision
424             floating point. When needed, using single-precision should not be an issue since
425             the minor errors introduced from loss of precision are drowned out by the errors
426             inherent in the haversine function that is used for distance calculations.
427              
428             =back
429              
430             =back
431              
432             =back
433              
434             =cut
435              
436              
437             #. Geo::Index uses C code to speed up distance computations.
438             #. Along with $C_CODE_COMPILED (set near top of this module), the following
439             #. variables hold the current state of the compiled code:
440              
441             my $ACTIVE_CODE = undef; #. Set to the type of low-level code currently being used:
442             #. 'perl, 'double', or 'float' (the latter two being C)
443             my @SUPPORTED_CODE = ( ); #. List of available low-level code types
444             my $C_CODE_ACTIVE = 0; #. Set true when compiled C code is being used
445              
446             sub new($;$$) {
447 19     19 1 12993 my ( $class, $_points, $_options ) = @_;
448            
449             #. Allow calling as Geo::Index->new( \%options )
450 19 100       98 if (ref $_points eq 'HASH') {
451 14         29 $_options = $_points;
452 14         31 $_points = undef;
453             }
454            
455             #. Initialize instance variables
456            
457 19   33     148 my Geo::Index $self = fields::new(ref $class || $class);
458            
459 19         35520 $self->{index} = { }; #. The points index
460 19         56 $self->{indices} = { }; #. Indices used for each point
461 19         41 $self->{positions} = { }; #. Each point's position when indexed
462            
463             #. Planetary parameters
464             #.
465             #. (Defaults are for the Earth)
466            
467 19         44 $self->{planetary_radius} = 6371230; #. Average radius of the object the points lie on (in m)
468 19         35 $self->{polar_circumference} = 40007863; #. Polar circumference of the object the points lie on (in m)
469 19         36 $self->{equatorial_circumference} = 40075017; #. Circumference at the equator of the object the points lie on (in m)
470            
471 19 100       78 if (ref $_options eq 'HASH') {
472            
473 14 50       50 if ($_options->{radius}) {
474             #. A custom planetary size is being used
475 0         0 $self->{planetary_radius} = $_options->{radius};
476            
477             #. If not specified, circumferences are calculated from radius
478 0 0       0 $self->{polar_circumference} = ( $_options->{polar_circumference} ) ? $_options->{polar_circumference} : 2 * pi * $self->{planetary_radius};
479 0 0       0 $self->{equatorial_circumference} = ( $_options->{equatorial_circumference} ) ? $_options->{equatorial_circumference} : 2 * pi * $self->{planetary_radius};
480             }
481            
482             }
483 19         54 $self->{planetary_diameter} = 2 * $self->{planetary_radius};
484            
485             #. Index parameters
486            
487             #. Level Grid Size Level Grid Size
488             #. ----- --------- --------------- ----- ------- -------
489             #. -1 ~40000 km (entire planet) 12 ~5 km ~2.4 GB
490             #. 0 ~20000 km 13 ~2.5 km ~2.7 GB
491             #. 1 ~10000 km ~1.0 GB 14 ~1.2 km ~3.1 GB
492             #. 2 ~5000 km ~1.0 GB 15 ~600 m ~3.3 GB
493             #. 3 ~2500 km ~1.1 GB 16 ~300 m ~3.6 GB
494             #. 4 ~1250 km ~1.2 GB 17 ~150 m ~3.8 GB
495             #. 5 ~625 km ~1.3 GB 18 ~75 m ~4.1 GB
496             #. 6 ~315 km ~1.4 GB 19 ~40 m ~4.4 GB
497             #. 7 ~155 km ~1.5 GB 20 ~20 m ~4.6 GB
498             #. 8 ~80 km ~1.6 GB 21 ~10 m ~4.9 GB
499             #. 9 ~40 km ~1.7 GB 22 ~5 m ~5.1 GB
500             #. 10 ~20 km ~1.9 GB 23 ~2 m ~5.4 GB
501             #. 11 ~10 km ~2.1 GB 24 ~1 m ~5.6 GB
502             #.
503             #. Memory usage shown above is for 1 million random points using numeric keys
504             #. at various settings of the 'levels' option ('Level' is one less than the
505             #. value of the 'levels' option with 'levels' values of -1 or 0 being invalid.)
506             #. For reference, the memory usage of the array of 1 million random, unindexed
507             #. points is about 440 MB; this amount is included in the index memory usage
508             #. shown above. The size of the points array will grow moderately (about 100
509             #. bytes per point) as the index is used.
510            
511             #. Depth of index
512 19         37 $self->{levels} = 20; #. ~40 m grid size at most detailed level of index
513 19 100       55 if (ref $_options eq 'HASH') {
514 14 100       56 $self->{levels} = int $_options->{levels} if ($_options->{levels});
515             }
516             #.Clip value
517 19 50       107 if ($self->{levels} > 30) {
    50          
518 0         0 $self->{levels} = 30;
519             } elsif ($self->{levels} < 1) {
520 0         0 $self->{levels} = 1;
521             }
522            
523             #. Number of grid tiles in each direction at most detailed level of index
524 19         58 $self->{max_size} = 2**$self->{levels};
525            
526             #. Index of the highest-resolution level
527 19         47 $self->{max_level} = $self->{levels} - 1;
528            
529             #. If possible, compiled C code will be used for the distance functions
530             #. If the code compiled and this is the first time new() is being called
531             #. we'll start using the compiled code by default.
532 19 100       56 unless ( defined $ACTIVE_CODE ) {
533 10         37 push @SUPPORTED_CODE, 'perl';
534            
535 10 50       33 if ( $C_CODE_COMPILED) {
536 10         24 push @SUPPORTED_CODE, 'float';
537 10         20 push @SUPPORTED_CODE, 'double';
538 10         32 SetDistanceFunctionType('double');
539             } else {
540 0         0 SetDistanceFunctionType('perl');
541             }
542             }
543            
544             #. Optionally switch function type to one of 'perl', 'double' or 'float'
545 19 100       63 if (ref $_options eq 'HASH') {
546 14 100       52 if ( defined $_options->{function_type} ) {
547 6         24 SetDistanceFunctionType( $_options->{function_type} );
548             }
549             }
550            
551             #. Optionally initialize the index with a set of points
552 19 0 50     59 if ( (defined $_points) && (scalar @$_points) ) {
553 0         0 $self->IndexPoints($_points);
554             }
555            
556 19 100       120 $self->{quiet} = ( defined $_options->{quiet} ) ? $_options->{quiet} : 0;
557            
558 19         303 return $self;
559             }
560              
561              
562              
563             =head2 IndexPoints( ... )
564              
565             =over
566              
567             C<$index-EIndexPoints( \@points );>
568              
569             Add points in list to the index
570              
571             If a point is added that already exists in the index and its position has
572             changed then the existing index entry will be deleted and the point will be
573             indexed again. If its position has not changed then no action will be taken.
574              
575             B>
576              
577             =over
578              
579             The points to add to the index
580              
581             Each point in the list is either a reference to a hash containing at a minimum
582             a C and a C value (both in degrees) or a reference to an array
583             giving the point. See the B> section above for details.
584              
585             =back
586              
587             =back
588              
589             =cut
590              
591              
592              
593             sub IndexPoints($$) {
594 15     15 1 5772 my ($self, $_points) = @_;
595            
596 15         59 foreach my $_point (@$_points) {
597 685         1286 my $type = ref $_point;
598            
599 685 50 33     2371 if ( ($type eq 'ARRAY') || ($type eq 'HASH') ) {
600 685         1486 $self->Index( $_point );
601            
602             } else {
603 0         0 croak "Unknown argument in Index: '$_point'";
604             }
605             }
606             }
607              
608              
609              
610             sub BuildPoints($) {
611 0     0 0 0 my (undef, $_in) = @_;
612 0         0 my @out = ( );
613            
614 0         0 foreach my $_point ( @$_in ) {
615 0         0 my $type = ref $_point;
616            
617 0 0       0 if ($type eq 'ARRAY') {
    0          
618             #. Got array; expand arguments into a full point
619 0         0 $_point = { 'lat'=>$$_point[0], 'lon'=>$$_point[1], 'data'=>$$_point[2] };
620            
621             } elsif ($type eq 'HASH') {
622             #. Got hash; no changes needed
623            
624             } else {
625 0         0 croak "Geo::Index::BuildPoints(...): Unknown point format '$_point'; maybe you passed a list of references instead of a reference to a list of references?\n";
626             }
627            
628 0         0 push @out, $_point;
629             }
630            
631 0 0       0 return ( wantarray ) ? @out : \@out;
632             }
633              
634              
635              
636             =head2 Index( ... )
637              
638             =over
639              
640             C<$index-EIndex( \%point );>
641              
642             C<$index-EIndex( \@point );>
643              
644             Add a single point to the index
645              
646             If the point being added already exists in the index and its position has
647             changed then the existing index entry will be deleted and the point will be
648             indexed again. If its position has not changed then no action will be taken.
649              
650             B> or B>
651              
652             =over
653              
654             The point to add to the index
655              
656             This can be either a hash containing at a minimum a C and a C value
657             (both in degrees) or an array giving the point. See the B> section above
658             for details.
659              
660             =back
661              
662             =back
663              
664             =cut
665              
666              
667             sub Index($$) {
668 687     687 1 1199 my ($self, $_point) = @_;
669            
670 687         983 my $type = ref $_point;
671            
672 687 100       1445 if ($type eq 'ARRAY') {
    50          
673             #. Got array; expand arguments into a full point
674 1         4 $_point = { 'lat'=>$$_point[0], 'lon'=>$$_point[1], 'data'=>$$_point[2] };
675            
676             } elsif ($type eq 'HASH') {
677             #. Got hash; no changes needed
678            
679             } else {
680 0         0 croak "Unknown argument in Index: '$_point'";
681             }
682            
683 687         1198 my $lat = $$_point{lat};
684 687         887 my $lon = $$_point{lon};
685            
686             #: Don't reindex points that are already in the index if they haven't moved
687 687 50       1574 if ( defined $$self{positions}{$_point} ) {
688            
689 0         0 my $indexed_position = $$self{positions}{$_point};
690            
691 0 0 0     0 if (
692             ( $lat == $$indexed_position[0] ) &&
693             ( $lon == $$indexed_position[1] )
694             ) {
695             #. Point is already indexed but hasn't moved; don't reindex
696 0         0 return;
697             }
698            
699             #. Point has moved; delete it so that it can be reindexed
700 0         0 Unindex($self, $_point);
701             }
702            
703 687         2111 $$self{positions}{$_point} = [ $lat, $lon ];
704            
705 687         1917 $$_point{lat_rad} = Math::Trig::deg2rad($lat);
706 687         6326 $$_point{lon_rad} = Math::Trig::deg2rad($lon);
707 687         5462 $$_point{circumference} = cos( $$_point{lat_rad} ) * $self->{equatorial_circumference};
708            
709 687         1443 my ($lat_int, $lon_int) = $self->GetIntLatLon($lat, $lon);
710            
711 687         1120 my $size = $self->{max_size};
712            
713 687         1052 my @to_index = ();
714 687         804 my $key;
715            
716 687         1380 my ($lat_idx_0, $lon_idx_0) = $self->GetIndices($self->{max_level}, $lat_int, $lon_int);
717            
718 687         1564 for (my $grid_level=$self->{max_level}; $grid_level>=0; $grid_level--) {
719            
720             #. Choose indices
721            
722 13740         16697 my $size_minus_one = $size - 1;
723 13740 100       23220 if ( $lat_idx_0 == 0 ) {
    100          
724             #. Near south pole
725            
726 464         537 if (USE_NUMERIC_KEYS) {
727 464         525 if (USE_PACKED_KEYS) {
728             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | MASK_LATLON );
729             } else {
730 464         682 $key = ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | MASK_LATLON;
731             }
732             } else {
733             $key = [ $grid_level, $lat_idx_0, ALL ];
734             }
735 464         697 push @to_index, $key;
736 464         533 if (USE_NUMERIC_KEYS) {
737 464         545 push @{ $self->{index}{$key} }, $_point;
  464         1199  
738             } else {
739             $self->AddValue( $key, $_point );
740             }
741            
742             } elsif ( $lat_idx_0 >= $size_minus_one ) {
743             #. Near north pole
744            
745             # Clip value
746 803         937 $lat_idx_0 = $size_minus_one;
747            
748 803         989 if (USE_NUMERIC_KEYS) {
749 803         902 if (USE_PACKED_KEYS) {
750             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | MASK_LATLON );
751             } else {
752 803         1180 $key = ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | MASK_LATLON;
753             }
754             } else {
755             $key = [ $grid_level, $size_minus_one, ALL ];
756             }
757 803         1258 push @to_index, $key;
758 803         947 if (USE_NUMERIC_KEYS) {
759 803         955 push @{ $self->{index}{$key} }, $_point;
  803         1820  
760             } else {
761             $self->AddValue( $key, $_point );
762             }
763            
764             } else {
765             #. Non-polar
766            
767 12473         13934 if (USE_NUMERIC_KEYS) {
768 12473         13741 if (USE_PACKED_KEYS) {
769             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | $lon_idx_0 );
770             } else {
771 12473         16758 $key = ( $grid_level << 58 ) | ( $lat_idx_0 << 29 ) | $lon_idx_0;
772             }
773             } else {
774             $key = [ $grid_level, $lat_idx_0, $lon_idx_0 ];
775             }
776            
777 12473         19694 push @to_index, $key;
778 12473         14223 if (USE_NUMERIC_KEYS) {
779 12473         13941 push @{ $self->{index}{$key} }, $_point;
  12473         38413  
780             } else {
781             $self->AddValue($key, $_point);
782             }
783             }
784            
785 13740         18008 $size >>= 1;
786 13740         16255 $lat_idx_0 >>= 1;
787 13740         23579 $lon_idx_0 >>= 1;
788             }
789            
790             #. All points in the world get this index
791 687         823 if (USE_NUMERIC_KEYS) {
792 687         798 if (USE_PACKED_KEYS) {
793             $key = pack("Q", ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON );
794             } else {
795 687         880 $key = ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON;
796             }
797             } else {
798             $key = [ ALL, ALL, ALL ];
799             }
800 687         1044 push @to_index, $key;
801 687         855 if (USE_NUMERIC_KEYS) {
802 687         808 push @{ $self->{index}{$key} }, $_point;
  687         1358  
803             } else {
804             $self->AddValue($key, $_point);
805             }
806            
807 687         3038 $$self{indices}{$_point} = \@to_index;
808             }
809              
810              
811              
812             =head2 Unindex( ... )
813              
814             =over
815              
816             C<$index-EUnindex( \%point );>
817              
818             Remove specified point from index
819              
820             This method will remove the point from the index but will not destroy the
821             actual point.
822              
823             B>
824              
825             =over
826              
827             The point to remove from the index.
828              
829             Note that this must be a reference to the actual point to remove and not to a
830             copy of it. Simply specifying a point's location as a new point hash will not
831             work.
832              
833             =back
834              
835             Added in 0.0.4 to replace the functionally identical (and now deprecated)
836             C.
837              
838             =back
839              
840             =cut
841              
842              
843             #. Trampoline to handle deprecated method name
844             sub DeletePointIndex($$) {
845 0     0 1 0 my ($self, $_point) = @_;
846            
847             print STDERR "DeletePointIndex(...) is deprecated and will be removed in future. " .
848             "Please update your code to use Unindex(...) instead.\n"
849 0 0       0 unless $self->{quiet};
850            
851             #. Update method pointer to point to new code
852 0         0 *Geo::Index::DeletePointIndex = *Geo::Index::Unindex;
853            
854             #. Fall through to new name
855 0 0       0 if (wantarray) {
856 0         0 return $self->Unindex($_point);
857             } else {
858 0         0 return scalar $self->Unindex($_point);
859             }
860             }
861              
862              
863             #. Remove specified point from index
864             # Used by Index
865             sub Unindex($$) {
866 5     5 1 14 my ($self, $_point) = @_;
867            
868             #. Remove the point from the index
869            
870             #. Get the full index
871 5         14 my $_index = $self->{index};
872            
873             #. Get the list of point's index keys
874 5         12 my $_indices = $$self{indices}{$_point};
875            
876             #. Loop through point's index keys...
877 5         10 foreach my $key ( @$_indices ) {
878            
879             #. Look up the key's index entry
880 105         168 my $_index_entry = $self->GetValue($key);
881            
882             #. Remove point from the index entry
883            
884 105         128 my $i = 0;
885             #. Loop through points lying in the index entry...
886 105         141 foreach my $_indexed ( @$_index_entry ) {
887 142 100       256 if ($_point eq $_indexed) {
888             #. Found point in index, delete it from the index
889 105         139 splice @$_index_entry, $i, 1;
890 105         151 last;
891             }
892 37         52 $i++;
893             }
894             #. Point has now been removed from index entry
895            
896             #. Delete the index entry if it's now empty
897            
898 105         120 my $formatted_key;
899 105         116 if (USE_NUMERIC_KEYS) {
900 105         125 $formatted_key = $key;
901             } else {
902             my $level = $$key[0];
903             my $lat = int $$key[1];
904             my $lon = int $$key[2];
905             if ( $lon == ALL ) {
906             if ( $level == ALL ) {
907             $formatted_key = 'ALL:ALL,ALL';
908             } else {
909             $formatted_key = "$level:$lat,ALL";
910             }
911             } else {
912             $formatted_key = "$level:$lat,$lon";
913             }
914            
915             }
916 105 100       115 delete $$_index{$formatted_key} unless (scalar @{$$_index{$formatted_key}});
  105         269  
917             }
918            
919             #. Delete the point from the index metadata
920 5         11 delete $$self{indices}{$_point};
921 5         31 delete $$self{positions}{$_point};
922             }
923              
924              
925              
926             #. Add a value to the index using a text key
927             # Used by Index
928             sub AddValue($$$) {
929 0     0 0 0 my ($self, $key, $value) = @_;
930            
931 0         0 my $_index = $self->{index};
932 0         0 if (USE_NUMERIC_KEYS) {
933 0         0 push @{ $$_index{$key} }, $value;
  0         0  
934             } else {
935             my $level = $$key[0];
936             my $lat = int $$key[1];
937             my $lon = int $$key[2];
938             if ( $lon == ALL ) {
939             if ( $level == ALL ) {
940             push @{ $$_index{'ALL:ALL,ALL'} }, $value;
941             } else {
942             push @{ $$_index{"$level:$lat,ALL"} }, $value;
943             }
944             } else {
945             push @{ $$_index{"$level:$lat,$lon"} }, $value;
946             }
947             }
948             }
949              
950              
951              
952             #. Return the index entry for a given key
953             #. Keys are either 64-bit integers or array
954             #. references: [ level, lat_idx, lon_idx ]
955             # Used by Search, Closest, Unindex, Sweep, Vacuum
956             sub GetValue($$) {
957 105     105 0 158 my ($self, $key) = @_;
958            
959 105         127 my $_index = $self->{index};
960 105         121 if (USE_NUMERIC_KEYS) {
961 105         188 return $$_index{$key};
962             } else {
963             my $level = $$key[0];
964             my $lat = int $$key[1];
965             my $lon = int $$key[2];
966             if ( $lon == ALL ) {
967             if ( $level == ALL ) {
968             return $$_index{"ALL:ALL,ALL"};
969             } else {
970             return $$_index{"$$key[0]:$$key[1],ALL"};
971             }
972             } else {
973             return $$_index{"$$key[0]:$$key[1],$$key[2]"};
974             }
975             }
976             }
977              
978              
979              
980              
981              
982              
983              
984             # ==============================================================================
985              
986              
987             =head2 Search( ... )
988              
989             =over
990              
991             C<@results = $index-ESearch( \%point, \%options );>
992              
993             C<$results = $index-ESearch( \%point, \%options );>
994              
995             Search index for points near a given point
996              
997             B>
998              
999             =over
1000              
1001             The point to search near
1002              
1003             This is either a reference to a hash containing at a minimum a C and a
1004             C value (both in degrees) or a reference to an array giving the point.
1005             See the B> section above for details.
1006              
1007             =back
1008              
1009             B>
1010              
1011             =over
1012              
1013             The parameters for the search (all are optional):
1014              
1015             Note that except for C, none of the below options have any effect when
1016             C is specified.
1017              
1018             B>
1019              
1020             =over
1021              
1022             Only return results within this distance (in meters) from search point.
1023              
1024             If no C is specified or the C is set to C then
1025             all points in the index will be returned.
1026              
1027             When C is specified then all points within the specified
1028             radius will be returned (additional points outside the radius may also
1029             be returned).
1030              
1031             =back
1032              
1033             B>
1034              
1035             =over
1036              
1037             Sort results by distance from search point before filtering and returning them.
1038              
1039             =back
1040              
1041             B>
1042              
1043             =over
1044              
1045             Return at most this many results.
1046              
1047             Unless sorting is also requested these are not guaranteed to be the closest
1048             results to the search point.
1049              
1050             =back
1051              
1052             B>
1053              
1054             =over
1055              
1056             Reference to additional user-supplied code to determine whether each point
1057             should be included in the results.
1058              
1059             This code is run before the distance from the search point to the result point
1060             has been calculated.
1061              
1062             See the B> section below for syntax.
1063              
1064             =back
1065              
1066             B>
1067              
1068             =over
1069              
1070             Reference to additional user-supplied code to determine whether each point
1071             should be included in the results.
1072              
1073             This code is run after the distance from the search point to the result point
1074             has been calculated.
1075              
1076             See the B> section below for syntax.
1077              
1078             =back
1079              
1080             B>
1081              
1082             =over
1083              
1084             Arbitrary user-supplied data that is passed to the condition functions.
1085              
1086             This can be used to allow the function access to additional data structures.
1087              
1088             See the B> section below for syntax.
1089              
1090             =back
1091              
1092             B>
1093              
1094             =over
1095              
1096             Return the raw preliminary results. (no result filtering is done)
1097              
1098             Normally when a search is performed the raw results have their distances
1099             computed and are filtered against the search radius and the condition functions.
1100             This process can be very slow when there are large numbers of points in a result
1101             set. Specifying this option skips those steps and can make searches several
1102             orders of magnitude faster.
1103              
1104             When this option is active then instead of returning a list of points the
1105             C method will return a list of lists of points (some of which
1106             may be C). If a C was specified then the results returned will
1107             contain all the points within that radius. Additional points outside of the
1108             search radius will likely be returned as well.
1109              
1110             When iterating over the arrays be sure to check whether a list element is C
1111             before trying to deference it.
1112              
1113             An example of returned quick results (in scalar context); B> are references
1114             to different points:
1115              
1116             C<[ [ POINT, POINT, POINT ], [ POINT, POINT ], undef, [ POINT, POINT ] ]>
1117              
1118             To be clear, when this option is active rough radius limiting is done but there
1119             is no filtering done, no distances are computed, and no sorting is performed.
1120              
1121             See the B> section below for a discussion of this option and when
1122             to S.
1123              
1124             =back
1125              
1126             B>
1127              
1128             =over
1129              
1130             Manual adjustment for tile size (signed integer, default is C<0>)
1131            
1132             Values C 0> use smaller tiles, values C 0> use larger tiles. Each
1133             increment of C<1> doubles or halves the tile size used. For example, set to
1134             C<-1> to use tiles half the normal size in each direction.
1135              
1136             This option can be a bit counter-intuitive. Although using smaller tiles will
1137             result in fewer points that need to be checked or returned it will also lead to
1138             a larger number of tiles that need to be processed. This can slow things down
1139             under some circumstances. Similarly using larger tiles results in more points
1140             spread over fewer tiles. What adjustment (if any) will result in the highest
1141             performance is highly dependant on both the search radius and on the number and
1142             distribution of the indexed points. If you adjust this value be sure to
1143             benchmark your application using a real dataset and the parameters (both typical
1144             and worst-case) that you expect to use.
1145              
1146             =back
1147              
1148             =back
1149              
1150             B
1151              
1152             =over
1153              
1154             In list context the return value is a list of references to the points found or
1155             an empty array if none were found.
1156              
1157             In scalar context the return value is a reference to the aforementioned list or
1158             C if no results were found.
1159              
1160             If either the C or C options were specified in the
1161             search options then for each point in the results the distance in meters from
1162             it to the search point will be stored in the C entry
1163             in the result point's hash. It can be retrieved using e.g.
1164             S>
1165              
1166             See above section for the results returned when the C option is
1167             active.
1168              
1169             =back
1170              
1171             =back
1172              
1173             =cut
1174              
1175              
1176             sub Search($$;$) {
1177 210     210 1 164994 my ($self, $_search_point, $_options) = @_;
1178            
1179 210         430 my $_points = $$self{index};
1180            
1181 210 100       606 if (ref $_search_point eq 'ARRAY') {
1182             #. Got array; expand arguments into a full point
1183 2         5 my $lat = $$_search_point[0];
1184 2         3 my $lon = $$_search_point[1];
1185            
1186 2         7 $_search_point = { 'lat'=>$lat, 'lon'=>$lon };
1187             }
1188            
1189 210         371 my $p_lat = $$_search_point{lat};
1190 210         290 my $p_lon = $$_search_point{lon};
1191            
1192             #. Search options; user should omit (or set to undef) inactive options:
1193            
1194 210         334 my $pre_condition = $$_options{pre_condition}; #. Reference to subroutine returning true if current point should be considered as
1195             #. a possible result, false otherwise. This subroutine should not modify any data.
1196             #. This subroutine is called before the distance from the search point to the
1197             #. result point has been calculated.
1198             #.
1199 210         287 my $post_condition = $$_options{post_condition}; #. Reference to subroutine returning true if current point should be considered as
1200             #. a possible result, false otherwise. This subroutine should not modify any data.
1201             #. This subroutine is called after the distance from the search point to the
1202             #. result point has been calculated.
1203             #.
1204             # $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
1205             #.
1206 210         636 my $search_radius = $$_options{radius}; #. Only points within radius (in meters) will be considered.
1207             #.
1208             # $$_options{sort_results}; #. Sort results by distance from point.
1209             #.
1210 210         273 my $max_results = $$_options{max_results}; #. Return at most this many results.
1211             #.
1212 210         275 my $quick_results = $$_options{quick_results}; #. Return preliminary results only. Do not compute distances or call the
1213             #. condition subroutines. Format returned is either a list of lists of points or
1214             #. a reference to a list of list of points (depending on how Search was called).
1215             #.
1216 210         271 my $tile_adjust = $$_options{tile_adjust}; #. Manual adjustment for tile size (signed integer, default is 0)
1217             #. Values <0 use smaller tiles, values >0 use larger tiles.
1218             #. Each increment of 1 doubles or halves the tile size used.
1219             #. For example, set to -1 to use tiles half the normal size in each direction.
1220             #.
1221             #. This option can be a bit counter-intuitive. Although using smaller tiles will
1222             #. result in fewer points that need to be checked or returned it will also lead to
1223             #. a larger number of tiles that need to be processed. This can slow things down
1224             #. under some circumstances. Similarly using larger tiles results in more points
1225             #. spread over fewer tiles. What adjustment (if any) will result in the highest
1226             #. performance is highly dependant on both the search radius and on the number
1227             #. and distribution of the indexed points. If you adjust this value be sure to
1228             #. benchmark your application using a real dataset and the parameters (both
1229             #. typical and worst-case) that you expect to use.
1230            
1231 210 100       450 $tile_adjust = ( defined $tile_adjust ) ? int $tile_adjust : 0;
1232            
1233 210 100       367 $quick_results = (defined $quick_results) ? 1 : 0;
1234            
1235 210 100       395 $search_radius = ALL unless defined ($search_radius);
1236            
1237 210         321 my $_results = [ ];
1238 210         354 my @result_set = ();
1239            
1240 210         266 my $p_lat_rad;
1241 210 100       395 if (defined $$_search_point{lat_rad}) {
1242 208         321 $p_lat_rad = $$_search_point{lat_rad};
1243             } else {
1244 2         7 $p_lat_rad = Math::Trig::deg2rad($p_lat);
1245 2         22 $$_search_point{lat_rad} = $p_lat_rad;
1246             }
1247            
1248 210         259 my $p_lon_rad;
1249 210 100       358 if (defined $$_search_point{lon_rad}) {
1250 208         325 $p_lon_rad = $$_search_point{lon_rad};
1251             } else {
1252 2         5 $p_lon_rad = Math::Trig::deg2rad($p_lon);
1253 2         15 $$_search_point{lon_rad} = $p_lon_rad;
1254             }
1255            
1256             #. Variables used while computing area extrema
1257 210         328 my $max_size = $self->{max_size};
1258 210         330 my $max_level = $self->{max_level};
1259             #. Set earlier: $p_lat, $p_lon, $self->{polar_circumference}, $search_radius
1260            
1261             #. Variables set/used while computing area extrema, used to select sets
1262 210         953 my $grid_level; #. Grid level to pull results from
1263             my $grid_size; #. Width or height of grid at chosen grid level
1264 210         0 my $max_grid_idx; #. Highest index in grid at chosen grid level
1265 210         0 my $lat_0_idx; #. Western extreme of search area (as grid index)
1266 210         0 my $lat_1_idx; #. Eastern extreme of search area (as grid index)
1267 210         0 my $lon_0_idx; #. Southern extreme of search area (as grid index)
1268 210         0 my $lon_1_idx; #. Northern extreme of search area (as grid index)
1269            
1270 210 100 100     1004 if (
1271             ( $search_radius == ALL ) || #. Explictly asked to search all points
1272             ( $search_radius > $self->{equatorial_circumference} / 4.0 ) #. A search radius over half the globe will search all points
1273             ) {
1274             #. Over half the globe covered by search so search all points
1275             #. Distances will be calculated but not checked
1276            
1277             #. KEY: [ ALL, ALL, ALL ]
1278            
1279 90         133 if (USE_NUMERIC_KEYS) {
1280 90         126 my $key;
1281 90         117 if (USE_PACKED_KEYS) {
1282             $key = pack("Q", ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON );
1283             } else {
1284 90         125 $key = ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON;
1285             }
1286 90         222 push @result_set, $$_points{$key};
1287            
1288             } else {
1289             #. Use split keys
1290            
1291             my $key = [ ALL, ALL, ALL ]; #. All points in the world get this index
1292             push @result_set, $self->GetValue($key);
1293             }
1294            
1295             } else {
1296             #. Less than half the globe being searched
1297            
1298 120 100       248 if ( $C_CODE_ACTIVE ) { # *** C ***
1299             #. The C code performs the same steps as the Perl version
1300            
1301             #. Dummy function that gets replace with the appropriate
1302             #. C version using either floats or doubles.
1303       0 0   sub ComputeAreaExtrema( $$$$$$$$ ) { }
1304            
1305 91         689 ( $grid_level, $grid_size, $max_grid_idx, $lat_0_idx, $lat_1_idx, $lon_0_idx, $lon_1_idx ) = ComputeAreaExtrema( $tile_adjust, $max_size, $max_level, $p_lat, $p_lat_rad, $p_lon, $self->{polar_circumference}, $search_radius );
1306            
1307             } else { # *** Perl ***
1308             # .--------------------------------------------.
1309             # | If you change code in this section be sure |
1310             # | to also update the C version of the code. |
1311             # '--------------------------------------------'
1312            
1313             #. Get search position in integer form
1314 29         77 my ($p_lat_int, $p_lon_int) = $self->GetIntLatLon($p_lat, $p_lon);
1315            
1316             #. Only return points within radius
1317            
1318             #. Determine grid level to search at
1319            
1320             #. Circumferences at this point (in meters)
1321 29         49 my $NS_circumference_in_meters = $self->{polar_circumference};
1322            
1323             #. Size of most detailed grid tile at this point (in meters)
1324             #. Distance is divided by two to get pole-to-pole distance
1325 29         50 my $ns_indexed_meters = ( $NS_circumference_in_meters / 2.0 ) / $max_size;
1326            
1327             #. The north-south size (in meters) of a grid tile is always
1328             #. larger than the east-west size (equal at the equator).
1329             #.
1330             #. Determine which levels tiles are larger than the search radius
1331             #.
1332             #. $search_radius / $ns_indexed_meters
1333             #. -> number of most detailed grid tiles at search radius
1334             #.
1335             #. ceil( log2( $search_radius / $ns_indexed_meters ) )
1336             #. -> relative grid level of tile no larger than search radius
1337             #.
1338             #. fast_log2 efficiently performs ceil( log2( n ) ); in C it calls those
1339             #. functions directly but a different method is used in the Perl code.
1340 29         64 my $shift = fast_log2( $search_radius / $ns_indexed_meters );
1341            
1342             #. Use a slightly higher-resolution grid
1343             #. We're betting that the time saved by searching a smaller physical area
1344             #. will outweigh the time needed for the additional index lookups
1345             #. In practice, this benchmarked over 40% slower than with no adjustment
1346             #$shift--;
1347            
1348 29         47 $shift += $tile_adjust;
1349            
1350             #. Use a slightly lower-resolution grid
1351             #. We're betting that the time saved by having to make fewer index lookups
1352             #. will outweigh the cost of having to search a larger physical area
1353             #.
1354             #. Using a lower resolution speeds things up when searching a small area
1355             #. (e.g. < 1 km radius at 1M random points) but makes things much slower
1356             #. when searching a large radius (e.g. > 10 km radius at 1M random points).
1357             #$shift += 3;
1358            
1359             #. Make sure the shift we computed lies within the index levels
1360 29 50       72 if ($shift < 0) {
    100          
1361 0         0 $shift = 0;
1362             } elsif ($shift >= $max_level) {
1363 1         3 $shift = $max_level - 1;
1364             }
1365            
1366             #. Shift is relative to the highest-resolution zoom level
1367            
1368             #. Determine grid level to use
1369 29         38 $grid_level = $max_level - $shift;
1370            
1371 29         49 $grid_size = 2**( $grid_level + 1 );
1372 29         43 $max_grid_idx = $grid_size - 1;
1373            
1374            
1375             #. Determine which grid tiles need to be checked
1376            
1377            
1378             #. Get search point's grid indices
1379            
1380             #. Determine number of degrees in one north-south meter
1381 29         41 my $lat_meter_in_degrees = 360.0 / $NS_circumference_in_meters;
1382            
1383             #. Searches are performed using angles (degrees or radians) for the radii, not meters
1384            
1385             #. Get north-south search radius
1386 29         41 my $lat_radius = $search_radius * $lat_meter_in_degrees;
1387 29         87 my $lat_radius_rad = deg2rad( $lat_radius );
1388            
1389             #. Get east-west search radius
1390             #. This is done as follows:
1391             #. o A point is placed on the equator at lat_radius east longitude.
1392             #. o Keeping the distance in meters from the the point to the prime
1393             #. meridian constant, the point is rotated north to the search point's
1394             #. latitude. The point's longitude will move east as this is done.
1395             #. o The point's new longitude is measured and this value is used for the
1396             #. lon_radius (east-west search radius).
1397 29         366 my $lon_radius = rad2deg( atan2( sin($lat_radius_rad), cos($lat_radius_rad) * cos($p_lat_rad) ) );
1398            
1399             #. The search radii have now been determined
1400            
1401             #. Determine the extreme latitudes and longitudes of the search circle
1402            
1403 29         220 my $lat_0 = $p_lat - $lat_radius;
1404 29         45 my $lat_1 = $p_lat + $lat_radius;
1405            
1406 29         45 my $lon_0 = $p_lon - $lon_radius;
1407 29         41 my $lon_1 = $p_lon + $lon_radius;
1408            
1409 29 100       79 if ( $lon_0 < -180.0 ) { $lon_0 += 360.0; }
  2 50       4  
1410 0         0 elsif ( $lon_0 > 180.0 ) { $lon_0 -= 360.0; }
1411            
1412 29 50       70 if ( $lon_1 < -180.0 ) { $lon_1 += 360.0; }
  0 100       0  
1413 6         11 elsif ( $lon_1 > 180.0 ) { $lon_1 -= 360.0; }
1414            
1415 29 100       76 if ( $lat_0 < -90.0 ) { $lat_0 = -90.0; }
  6 50       21  
1416 0         0 elsif ( $lat_0 > 90.0 ) { $lat_0 = 90.0; }
1417            
1418 29 50       66 if ( $lat_1 < -90.0 ) { $lat_1 = -90.0; }
  0 100       0  
1419 4         8 elsif ( $lat_1 > 90.0 ) { $lat_1 = 90.0; }
1420            
1421             #. Determine the grid indices for the search circle's extremes
1422            
1423             #. Inlined for speed:
1424             # $lat_0_idx = $self->GetIntLat($lat_0) >> $shift;
1425             # $lat_1_idx = $self->GetIntLat($lat_1) >> $shift;
1426             #
1427             # $lon_0_idx = $self->GetIntLon($lon_0) >> $shift;
1428             # $lon_1_idx = $self->GetIntLon($lon_1) >> $shift;
1429            
1430 29         54 $lat_0_idx = int( ( $lat_0 + 90.0 ) * $max_size / 180.0 );
1431 29 100       51 $lat_0_idx = $max_size - 1 if ($lat_0_idx >= $max_size);
1432 29         46 $lat_0_idx >>= $shift;
1433            
1434 29         47 $lat_1_idx = int( ( $lat_1 + 90.0 ) * $max_size / 180.0 );
1435 29 100       49 $lat_1_idx = $max_size - 1 if ($lat_1_idx >= $max_size);
1436 29         34 $lat_1_idx >>= $shift;
1437            
1438 29         55 $lon_0_idx = ( int( ( $lon_0 + 180.0 ) * $max_size / 360.0 ) % $max_size ) >> $shift;
1439            
1440 29         46 $lon_1_idx = ( int( ( $lon_1 + 180.0 ) * $max_size / 360.0 ) % $max_size ) >> $shift;
1441             #. END inline
1442            
1443             #. Make sure latitudes are ordered south then north
1444             #. (this is not always the case in polar regions)
1445             #. Longitudes can be in either order to allow straddling of antimeridian
1446 29 50       71 ($lat_0_idx, $lat_1_idx) = ($lat_1_idx, $lat_0_idx) if ($lat_0_idx > $lat_1_idx);
1447             } # END Perl
1448            
1449             #. Grid extrema have now been determined
1450            
1451             #. Gather preliminary search results
1452            
1453 120         202 if (USE_NUMERIC_KEYS) {
1454 120         175 my $seen_n_polar = 0;
1455 120         153 my $seen_s_polar = 0;
1456 120 100       215 if ($lon_0_idx <= $lon_1_idx) {
1457             #. Does not straddle antimeridian
1458 88         215 for (my $lat_idx = $lat_0_idx; $lat_idx <= $lat_1_idx; $lat_idx++) {
1459            
1460 123 100       259 if ( $lat_idx == 0 ) {
    100          
1461             #. Near south pole
1462 20 50       61 unless ( $seen_s_polar ) {
1463 20         28 $seen_s_polar = 1;
1464 20         26 my $key;
1465 20         56 if (USE_PACKED_KEYS) {
1466             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
1467             } else {
1468 20         52 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
1469             }
1470 20         85 push @result_set, $$_points{$key};
1471             }
1472            
1473             } elsif ( $lat_idx >= $max_grid_idx ) {
1474             #. Near north pole
1475 30 50       76 unless ( $seen_n_polar ) {
1476 30         46 $seen_n_polar = 1;
1477 30         36 my $key;
1478 30         43 if (USE_PACKED_KEYS) {
1479             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
1480             } else {
1481 30         62 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
1482             }
1483 30         99 push @result_set, $$_points{$key};
1484             }
1485            
1486             } else {
1487             #. Normal case
1488 73         153 for (my $lon_idx = $lon_0_idx; $lon_idx <= $lon_1_idx; $lon_idx++) {
1489 1391         1722 my $clipped_lon_idx = $lon_idx % $grid_size;
1490 1391         1569 my $key;
1491 1391         1538 if (USE_PACKED_KEYS) {
1492             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx );
1493             } else {
1494 1391         1912 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
1495             }
1496 1391         2706 push @result_set, $$_points{$key};
1497             }
1498             }
1499            
1500             }
1501             # END does not straddle antimeridian
1502            
1503             } else { # ($lon_0_idx > $lon_1_idx)
1504             #. Straddles antimeridian
1505 32         94 for (my $lat_idx = $lat_0_idx; $lat_idx <= $lat_1_idx; $lat_idx++) {
1506            
1507 72 100       162 if ( $lat_idx == 0 ) {
    50          
1508             #. Near south pole
1509 24 50       53 unless ( $seen_s_polar ) {
1510 24         34 $seen_s_polar = 1;
1511 24         41 my $key;
1512 24         33 if (USE_PACKED_KEYS) {
1513             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
1514             } else {
1515 24         69 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
1516             }
1517 24         77 push @result_set, $$_points{$key};
1518             }
1519            
1520             } elsif ( $lat_idx >= $max_grid_idx ) {
1521             #. Near north pole
1522 0 0       0 unless ( $seen_n_polar ) {
1523 0         0 $seen_n_polar = 1;
1524 0         0 my $key;
1525 0         0 if (USE_PACKED_KEYS) {
1526             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
1527             } else {
1528 0         0 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
1529             }
1530 0         0 push @result_set, $$_points{$key};
1531             }
1532            
1533             } else {
1534             #. Non-polar
1535            
1536             #. East side
1537 48         101 for (my $lon_idx = $lon_0_idx; $lon_idx <= $max_grid_idx; $lon_idx++) {
1538 21075         25778 my $clipped_lon_idx = $lon_idx % $grid_size;
1539 21075         23058 my $key;
1540 21075         24073 if (USE_PACKED_KEYS) {
1541             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx );
1542             } else {
1543 21075         27075 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
1544             }
1545 21075         40165 push @result_set, $$_points{$key};
1546             }
1547            
1548             #. West side
1549 48         91 for (my $lon_idx = 0; $lon_idx <= $lon_1_idx; $lon_idx++) {
1550 6247         7624 my $clipped_lon_idx = $lon_idx % $grid_size;
1551 6247         7095 my $key;
1552 6247         7170 if (USE_PACKED_KEYS) {
1553             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx );
1554             } else {
1555 6247         7847 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
1556             }
1557 6247         11858 push @result_set, $$_points{$key};
1558             }
1559             }
1560             }
1561             }
1562            
1563             } else {
1564             #. Use split keys
1565            
1566             my $seen_n_polar = 0;
1567             my $seen_s_polar = 0;
1568             if ($lon_0_idx <= $lon_1_idx) {
1569             #. Does not straddle antimeridian
1570             for (my $lat_idx = $lat_0_idx; $lat_idx <= $lat_1_idx; $lat_idx++) {
1571            
1572             if ( $lat_idx == 0 ) {
1573             #. Near south pole
1574             unless ( $seen_s_polar ) {
1575             $seen_s_polar = 1;
1576             my $key = [ $grid_level, $lat_idx, ALL ];
1577             push @result_set, $self->GetValue($key);
1578             }
1579            
1580             } elsif ( $lat_idx >= $max_grid_idx ) {
1581             #. Near north pole
1582             unless ( $seen_n_polar ) {
1583             $seen_n_polar = 1;
1584             my $key = [ $grid_level, $lat_idx, ALL ];
1585             push @result_set, $self->GetValue($key);
1586             }
1587            
1588             } else {
1589             #. Normal case
1590             for (my $lon_idx = $lon_0_idx; $lon_idx <= $lon_1_idx; $lon_idx++) {
1591             my $clipped_lon_idx = $lon_idx % $grid_size;
1592             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
1593             push @result_set, $self->GetValue($key);
1594             }
1595             }
1596            
1597             }
1598             # END does not straddle antimeridian
1599            
1600             } else { # ($lon_0_idx > $lon_1_idx)
1601             #. Straddles antimeridian
1602             for (my $lat_idx = $lat_0_idx; $lat_idx <= $lat_1_idx; $lat_idx++) {
1603            
1604             if ( $lat_idx == 0 ) {
1605             #. Near south pole
1606             unless ( $seen_s_polar ) {
1607             $seen_s_polar = 1;
1608             my $key = [ $grid_level, $lat_idx, ALL ];
1609             push @result_set, $self->GetValue($key);
1610             }
1611            
1612             } elsif ( $lat_idx >= $max_grid_idx ) {
1613             #. Near north pole
1614             unless ( $seen_n_polar ) {
1615             $seen_n_polar = 1;
1616             my $key = [ $grid_level, $lat_idx, ALL ];
1617             push @result_set, $self->GetValue($key);
1618             }
1619            
1620             } else {
1621             #. Non-polar
1622            
1623             #. East side
1624             for (my $lon_idx = $lon_0_idx; $lon_idx <= $max_grid_idx; $lon_idx++) {
1625             my $clipped_lon_idx = $lon_idx % $grid_size;
1626             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
1627             push @result_set, $self->GetValue($key);
1628             }
1629            
1630             #. West side
1631             for (my $lon_idx = 0; $lon_idx <= $lon_1_idx; $lon_idx++) {
1632             my $clipped_lon_idx = $lon_idx % $grid_size;
1633             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
1634             push @result_set, $self->GetValue($key);
1635             }
1636             }
1637             }
1638             }
1639            
1640             }
1641            
1642             }
1643            
1644 210 100       427 if ( $quick_results ) {
1645             #. Return preliminary results
1646             #. Format is a list of lists (some of which may be undef
1647             #. All points within the search radius will be returned
1648             #. possibly along with additional points outside the
1649             #. search radius.
1650             return ( wantarray )
1651             #. Return array:
1652 2 50       19 ? @result_set
1653             #. Return array reference:
1654             : \@result_set;
1655             }
1656            
1657            
1658             #. Gather results
1659            
1660             #. Prepare to compute result distances
1661            
1662 208         702 SetUpDistance($self->{planetary_diameter}, $p_lat_rad, $p_lon_rad);
1663            
1664 208 100       361 if ( $search_radius > 0 ) {
1665             #. Search for points within radius
1666            
1667 156 50 33     518 if (defined $pre_condition || defined $post_condition) {
1668             #. Filter specified
1669 0         0 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
1670            
1671 0 0       0 if (defined $pre_condition) {
1672 0 0       0 if (defined $post_condition) {
1673             #. Both pre- and post-distance calculation condition
1674 0         0 foreach my $_set ( @result_set ) {
1675 0 0       0 next unless (defined $_set);
1676 0         0 foreach my $_point (@$_set) {
1677 0 0       0 if ( &$pre_condition($_point, $_search_point, $user_data) ) {
1678 0         0 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1679 0 0       0 if ( $distance <= $search_radius ) {
1680 0         0 $$_point{search_result_distance} = $distance;
1681 0 0       0 if ( &$post_condition($_point, $_search_point, $user_data) ) {
1682 0         0 push @$_results, $_point;
1683             }
1684             }
1685             }
1686             }
1687             }
1688             } else {
1689             #. Only pre-distance calculation condition
1690 0         0 foreach my $_set ( @result_set ) {
1691 0 0       0 next unless (defined $_set);
1692 0         0 foreach my $_point (@$_set) {
1693 0 0       0 if ( &$pre_condition($_point, $_search_point, $user_data) ) {
1694 0         0 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1695 0 0       0 if ( $distance <= $search_radius ) {
1696 0         0 $$_point{search_result_distance} = $distance;
1697 0         0 push @$_results, $_point;
1698             }
1699             }
1700             }
1701             }
1702             }
1703             } else {
1704             #. Only post-distance calculation condition
1705 0         0 foreach my $_set ( @result_set ) {
1706 0 0       0 next unless (defined $_set);
1707 0         0 foreach my $_point (@$_set) {
1708 0         0 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1709 0 0       0 if ( $distance <= $search_radius ) {
1710 0         0 $$_point{search_result_distance} = $distance;
1711 0 0       0 if ( &$post_condition($_point, $_search_point, $user_data) ) {
1712 0         0 push @$_results, $_point;
1713             }
1714             }
1715             }
1716             }
1717             }
1718            
1719             # END conditions present
1720            
1721             } else {
1722             #. No filter
1723 156         286 foreach my $_set ( @result_set ) {
1724 28823 100       45783 next unless (defined $_set);
1725 214         331 foreach my $_point (@$_set) {
1726 642         1479 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1727 642 100       1823 if ( $distance <= $search_radius ) {
1728 494         629 $$_point{search_result_distance} = $distance;
1729 494         960 push @$_results, $_point;
1730             }
1731             }
1732             }
1733             }
1734            
1735             } else {
1736             #. Explictly asked to search all points
1737             #. Distances will be calculated but not checked
1738            
1739 52 100 100     216 if (defined $pre_condition || defined $post_condition) {
1740             #. Filter specified
1741 2         4 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
1742            
1743 2 100       6 if (defined $pre_condition) {
1744 1 50       3 if (defined $post_condition) {
1745             #. Both pre- and post-distance calculation condition
1746 0         0 foreach my $_set ( @result_set ) {
1747 0 0       0 next unless (defined $_set);
1748 0         0 foreach my $_point (@$_set) {
1749 0 0       0 if ( &$pre_condition($_point, $_search_point, $user_data) ) {
1750 0         0 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1751 0         0 $$_point{search_result_distance} = $distance;
1752 0 0       0 if ( &$post_condition($_point, $_search_point, $user_data) ) {
1753 0         0 push @$_results, $_point;
1754             }
1755             }
1756             }
1757             }
1758             } else {
1759             #. Only pre-distance calculation condition
1760 1         4 foreach my $_set ( @result_set ) {
1761 1 50       4 next unless (defined $_set);
1762 1         3 foreach my $_point (@$_set) {
1763 206 100       714 if ( &$pre_condition($_point, $_search_point, $user_data) ) {
1764 9         50 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1765 9         12 $$_point{search_result_distance} = $distance;
1766 9         21 push @$_results, $_point;
1767             }
1768             }
1769             }
1770             }
1771             } else {
1772             #. Only post-distance calculation condition
1773 1         3 foreach my $_set ( @result_set ) {
1774 1 50       4 next unless (defined $_set);
1775 1         3 foreach my $_point (@$_set) {
1776 206         911 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1777 206         270 $$_point{search_result_distance} = $distance;
1778 206 100       293 if ( &$post_condition($_point, $_search_point, $user_data) ) {
1779 9         47 push @$_results, $_point;
1780             }
1781             }
1782             }
1783             }
1784            
1785             # END conditions present
1786            
1787             } else {
1788             #. No filter
1789 50         95 foreach my $_set ( @result_set ) {
1790 50 50       97 next unless (defined $_set);
1791 50         91 foreach my $_point (@$_set) {
1792 556         1216 my $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
1793 556         937 $$_point{search_result_distance} = $distance;
1794 556         907 push @$_results, $_point;
1795             }
1796             }
1797             }
1798             }
1799            
1800 208 100       428 if (defined $$_options{sort_results}) {
1801             #. Return points sorted by distance
1802 205         350 my $_tmp = $_results; $_results = undef;
  205         282  
1803 205         715 @$_results = sort { $$a{search_result_distance} <=> $$b{search_result_distance} } @$_tmp;
  2157         3344  
1804             }
1805            
1806 208 100       460 if ( $max_results ) {
1807 1         4 my $result_count = scalar @$_results;
1808 1 50       6 $max_results = $result_count if ($result_count < $max_results);
1809             }
1810            
1811             return ( wantarray )
1812             #. Return array:
1813 208 0       1511 ? ( scalar @$_results )
    0          
    0          
    50          
    100          
    50          
    50          
1814             ? ( $max_results )
1815             ? ( scalar @$_results > $max_results )
1816             ? @$_results[0..($max_results-1)]
1817             : @$_results
1818             : @$_results
1819             : ()
1820             #. Return array reference:
1821             : ( scalar @$_results )
1822             ? ( $max_results )
1823             ? ( scalar @$_results > $max_results )
1824             ? [ @$_results[0..($max_results-1)] ]
1825             : $_results
1826             : $_results
1827             : undef; #. undef == No result found
1828             }
1829              
1830              
1831              
1832              
1833              
1834              
1835              
1836              
1837              
1838             # ==============================================================================
1839              
1840              
1841             =head2 SearchByBounds( ... )
1842              
1843             =over
1844              
1845             C<@results = $index-ESearchByBounds( \@bounds, \%options );>
1846              
1847             C<@results = $index-ESearchByBounds( \%bounds, \%options );>
1848              
1849             C<$results = $index-ESearchByBounds( \@bounds, \%options );>
1850              
1851             C<$results = $index-ESearchByBounds( \%bounds, \%options );>
1852              
1853             Search index for points within a given bounding box
1854              
1855             The points returned are those that lie between the specified latitudes and
1856             longitudes. The four corners form a rectangle only when using certain map
1857             projections such as L
1858             or L
1859             (including L
1860             a.k.a. web mercator as used by slippy maps). If you are using a projection
1861             that does not have horizontal lines of latitude and vertical lines of longitude
1862             and you want the results to lie within and/or fill a rectangle on your map then
1863             your code will need to perform the appropriate bounds adjustment and point
1864             filtering itself.
1865              
1866             B> or B>
1867              
1868             =over
1869              
1870             The search boundaries
1871              
1872             This can be specified either as a list or a hash. Any of the following are
1873             acceptable:
1874              
1875             =over
1876              
1877             C<( I, I, I, I )>
1878              
1879             C<( 'south' =E I, 'north' =E I, 'west' =E I, 'east' =E I )>
1880              
1881             C<( 's' =E I, 'n' =E I, 'w' =E I, 'e' =E I )>
1882              
1883             C<( 'S' =E I, 'N' =E I, 'W' =E I, 'E' =E I )>
1884              
1885             =back
1886              
1887             C> and C> are the south and north latitudes of the bounding
1888             box and C> and C> are its west and east longitudes (all values
1889             are in degrees). For the list form the order matches that used by PostGIS,
1890             shapefiles, and GeoJSON but be aware that the order of the fields is not
1891             standardized.
1892              
1893             =back
1894              
1895             B>
1896              
1897             =over
1898              
1899             The parameters for the search (all are optional):
1900              
1901             Note that none of the below options have any effect when C is
1902             specified.
1903              
1904             =over
1905              
1906             B>
1907              
1908             =over
1909              
1910             Return at most this many results.
1911              
1912             =back
1913              
1914             B>
1915              
1916             =over
1917              
1918             Reference to additional user-supplied code to determine whether each point
1919             should be included in the results.
1920              
1921             Note that unlike with the other search methods there is only a single condition
1922             function. Instead of the C<$_search_point>, the second parameter to the
1923             condition function is a reference to the bounding box in list form (as described
1924             above for C<@bounds>). Lastly, since "distance from search point" makes no
1925             sense in the context of a bounding box, none is provided to the function.
1926              
1927             See the B> section below for syntax.
1928              
1929             =back
1930              
1931             B>
1932              
1933             =over
1934              
1935             Arbitrary user-supplied data that is passed to the condition function.
1936              
1937             This can be used to allow the function access to additional data structures.
1938              
1939             See the B> section below for syntax.
1940              
1941             =back
1942              
1943             B>
1944              
1945             =over
1946              
1947             Return the raw preliminary results. (no result filtering is done)
1948              
1949             Normally when a search is performed the raw results are filtered to ensure that
1950             they lie within the bounds and the condition function is applied to each result
1951             point. This process can be slow when there are large numbers of points in a
1952             result set. Specifying this option skips those steps and can make searches
1953             faster.
1954              
1955             When this option is active then instead of returning a list of points the
1956             C method will return a list of lists of points (some
1957             of which may be C). The results are guaranteed to contain all points
1958             within the requested bounds but additional points outside of the requested
1959             bounds will likely be returned as well.
1960              
1961             When iterating over the arrays be sure to check whether a list element is C
1962             before trying to deference it.
1963              
1964             An example of returned quick results (in scalar context); B> are references
1965             to different points:
1966              
1967             C<[ [ POINT, POINT, POINT ], [ POINT, POINT ], undef, [ POINT, POINT ] ]>
1968              
1969             To be clear, when this option is active rough bounds limiting is done but there
1970             is no filtering done and no bound checks are actually performed.
1971              
1972             See the B> section below for a discussion of this
1973             option and when to use it.
1974              
1975             =back
1976              
1977              
1978             B>
1979              
1980             =over
1981              
1982             Manual adjustment for tile size (signed integer, default is C<0>)
1983              
1984             Values C 0> use smaller tiles, values C 0> use larger tiles. Each
1985             increment of C<1> doubles or halves the tile size used. For example, set to
1986             C<-1> to use tiles half the normal size in each direction.
1987              
1988             This option can be a bit counter-intuitive. Although using smaller tiles will
1989             result in fewer points that need to be checked or returned it will also lead to
1990             a larger number of tiles that need to be processed. This can slow things down
1991             under some circumstances. Similarly using larger tiles results in more points
1992             spread over fewer tiles. What adjustment (if any) will result in the highest
1993             performance is highly dependant on both the search radius and on the number and
1994             distribution of the indexed points. If you adjust this value be sure to
1995             benchmark your application using a real dataset and the parameters (both typical
1996             and worst-case) that you expect to use.
1997              
1998             =back
1999              
2000             =back
2001              
2002             B
2003              
2004             =over
2005              
2006             In list context the return value is a list of references to the points found or
2007             an empty array if none were found.
2008              
2009             In scalar context the return value is a reference to the aforementioned list or
2010             C if no results were found.
2011              
2012             See above section for the results returned when the C option is
2013             active.
2014              
2015             =back
2016              
2017             =back
2018              
2019             =back
2020              
2021             =cut
2022              
2023              
2024 11     11   51422 use constant TRACE_BOUNDS => 0; # Uncomment for no tracing of this method
  11         48  
  11         98075  
2025             #use constant TRACE_BOUNDS => 1; # Uncomment to enable tracing of this method
2026              
2027             sub SearchByBounds($$;$) {
2028 21     21 1 21658 my ($self, $_bounding_box, $_options) = @_;
2029            
2030 21         47 my $_points = $$self{index};
2031            
2032             #. Determine bounding box
2033            
2034 21         39 my ( $south, $north, $west, $east );
2035            
2036 21 100       55 if (ref $_bounding_box eq 'HASH') {
2037             #. Got hash; convert bounds to array
2038            
2039             $west = $$_bounding_box{w} ? $$_bounding_box{w}
2040             : $$_bounding_box{W} ? $$_bounding_box{W}
2041             : $$_bounding_box{west} ? $$_bounding_box{west}
2042 14 50       53 : croak "Could not find west value for bounds";
    50          
    50          
2043            
2044             $south = $$_bounding_box{s} ? $$_bounding_box{s}
2045             : $$_bounding_box{S} ? $$_bounding_box{S}
2046             : $$_bounding_box{south} ? $$_bounding_box{south}
2047 14 50       41 : croak "Could not find south value for bounds";
    50          
    50          
2048            
2049             $east = $$_bounding_box{e} ? $$_bounding_box{e}
2050             : $$_bounding_box{E} ? $$_bounding_box{E}
2051             : $$_bounding_box{east} ? $$_bounding_box{east}
2052 14 50       46 : croak "Could not find east value for bounds";
    50          
    50          
2053            
2054             $north = $$_bounding_box{n} ? $$_bounding_box{n}
2055             : $$_bounding_box{N} ? $$_bounding_box{N}
2056             : $$_bounding_box{north} ? $$_bounding_box{north}
2057 14 50       49 : croak "Could not find north value for bounds";
    50          
    50          
2058            
2059 14         33 $_bounding_box = [ $west, $south, $east, $north ];
2060            
2061             } else {
2062             #. Bounds in array form
2063            
2064 7         13 $south = $$_bounding_box[1];
2065 7         11 $north = $$_bounding_box[3];
2066            
2067 7         12 $west = $$_bounding_box[0];
2068 7         10 $east = $$_bounding_box[2];
2069             }
2070            
2071             #. Make sure bounds are valid
2072            
2073 21         28 my $has_errors = 0;
2074            
2075 21 100       63 if ($west < -180.0) {
2076 1 50       4 carp "In SearchByBounds: west bound $west is out of range ( < -180 )" unless $self->{quiet};
2077 1         2 $has_errors = 1;
2078             }
2079            
2080 21 100       48 if ($east > 180.0) {
2081 1 50       4 carp "In SearchByBounds: east bound $east is out of range ( > 180 )" unless $self->{quiet};
2082 1         2 $has_errors = 1;
2083             }
2084            
2085 21 100       59 if ($south < -90.0) {
2086 1 50       2 carp "In SearchByBounds: south bound $south is out of range ( < -90 )" unless $self->{quiet};
2087 1         2 $has_errors = 1;
2088             }
2089            
2090 21 100       49 if ($north > 90.0) {
2091 1 50       5 carp "In SearchByBounds: north bound $north is out of range ( > 90 )" unless $self->{quiet};
2092 1         2 $has_errors = 1;
2093             }
2094            
2095 21 100       51 if ($south > $north) {
    50          
2096 1 50       4 carp "In SearchByBounds: south bound greater than north bound ( $south > $north )" unless $self->{quiet};
2097 1         4 $has_errors = 1;
2098             } elsif ($south == $north) {
2099 0 0       0 carp "In SearchByBounds: bounds cover no area ( south > north: $south > $north )" unless $self->{quiet};
2100 0         0 $has_errors = 1;
2101             }
2102            
2103 21 100       43 if ($has_errors) {
2104 2 50       12 return (wantarray) ? ( ) : undef;
2105             }
2106            
2107             #. Search options; user should omit (or set to undef) inactive options:
2108            
2109 19         35 my $condition = $$_options{condition}; #. Reference to subroutine returning true if current point should be considered as
2110             #. a possible result, false otherwise. This subroutine should not modify any data.
2111             #. This subroutine is called before the distance from the search point to the
2112             #. result point has been calculated.
2113             #.
2114             # $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
2115             #.
2116 19         32 my $max_results = $$_options{max_results}; #. Return at most this many results.
2117             #.
2118 19         30 my $quick_results = $$_options{quick_results}; #. Return preliminary results only. Do not compute distances or call the
2119             #. condition subroutines. Format returned is either a list of lists of points or
2120             #. a reference to a list of list of points (depending on how Search was called).
2121             #.
2122 19         33 my $tile_adjust = $$_options{tile_adjust}; #. Manual adjustment for tile size (signed integer, default is 0)
2123             #. Values <0 use smaller tiles, values >0 use larger tiles.
2124             #. Each increment of 1 doubles or halves the tile size used.
2125             #. For example, set to -1 to use tiles half the normal size in each direction.
2126             #.
2127             #. This option can be a bit counter-intuitive. Although using smaller tiles will
2128             #. result in fewer points that need to be checked or returned it will also lead to
2129             #. a larger number of tiles that need to be processed. This can slow things down
2130             #. under some circumstances. Similarly using larger tiles results in more points
2131             #. spread over fewer tiles. What adjustment (if any) will result in the highest
2132             #. performance is highly dependent on both the search radius and on the number
2133             #. and distribution of the indexed points. If you adjust this value be sure to
2134             #. benchmark your application using a real dataset and the parameters (both
2135             #. typical and worst-case) that you expect to use.
2136            
2137 19 100       42 $tile_adjust = ( defined $tile_adjust ) ? int $tile_adjust : 0;
2138            
2139 19 100       33 $quick_results = (defined $quick_results) ? 1 : 0;
2140            
2141 19         37 if ( TRACE_BOUNDS ) {
2142             print "Bounds: W: $west, S: $south, E: $east, N: $north\n";
2143             }
2144            
2145 19         29 my @keys; #. The index keys to return points from
2146 19         22 if ( TRACE_BOUNDS) {
2147             @keys = ();
2148             }
2149            
2150 19         51 my $_results = [ ];
2151 19         31 my @result_set = ();
2152            
2153 19         36 my $max_level = $self->{max_level};
2154            
2155             #. Variables set/used while computing area extrema, used to select sets
2156 19         52 my $grid_level; #. Grid level to pull results from
2157             my $grid_size; #. Width or height of grid at chosen grid level
2158 19         0 my $max_grid_idx; #. Highest index in grid at chosen grid level
2159 19         0 my $shift; #. Number of bits to shift integer latitudes and longitudes to yield valid indices at the active grid level
2160            
2161             #. Determine how many degrees the bounds cover in each direction
2162            
2163 19         32 my $lat_degrees = $north - $south;
2164 19 50       40 if ($lat_degrees <= 0) {
2165 0 0       0 carp "In SearchByBounds north ($north) is less than or equal to south ($south)!\n" unless $self->{quiet};
2166 0 0       0 return (wantarray) ? ( ) : undef;
2167             }
2168            
2169 19 100       49 my $lon_degrees = ($west <= $east) ? $east - $west #. Normal case
2170             : ( $east + 360.0 ) - $west; #. Straddles antimeridian
2171            
2172 19         52 if ( TRACE_BOUNDS) {
2173             print "my \$lon_degrees = $lon_degrees = ($west <= $east) ? $east - $west : ( $east + 360.0 ) - $west;\n";
2174             }
2175            
2176             #. Determine grid level to use
2177            
2178 19         98 my $lat_best_level = fast_log2( 180.0 / $lat_degrees );
2179 19         48 my $lon_best_level = fast_log2( 360.0 / $lon_degrees );
2180            
2181 19 100       42 $grid_level = ( $lat_best_level > $lon_best_level ) ? $lat_best_level : $lon_best_level;
2182 19         27 $grid_level -= $tile_adjust;
2183            
2184             #. Determine shift and grid size for the chosen level
2185            
2186 19         24 $shift = $max_level - $grid_level;
2187            
2188             #$shift += $tile_adjust;
2189            
2190 19         32 $grid_size = 2**( $grid_level + 1 );
2191 19         29 $max_grid_idx = $grid_size - 1;
2192            
2193             #. Compute grid indices for bounds extrema
2194            
2195 19         76 my ($south_int, $west_int) = $self->GetIntLatLon($south, $west);
2196 19         41 my ($north_int, $east_int) = $self->GetIntLatLon($north, $east);
2197            
2198 19         37 my $north_idx = $north_int >> $shift; #. Northern extreme of search area (as grid index)
2199 19         29 my $south_idx = $south_int >> $shift; #. Southern extreme of search area (as grid index)
2200 19         30 my $east_idx = $east_int >> $shift; #. Eastern extreme of search area (as grid index)
2201 19         28 my $west_idx = $west_int >> $shift; #. Western extreme of search area (as grid index)
2202            
2203 19 100       40 $east_idx = $max_grid_idx if ($east == 180.0); #. Special case for antimeridian as east bound
2204 19 100       37 $north_idx = $max_grid_idx if ($north == 90.0); #. Special case for north pole as north bound
2205            
2206 19 100       49 my $include_south_pole = ( $south == -90.0 ) ? 1 : 0;
2207 19 100       42 my $include_north_pole = ( $north == 90.0 ) ? 1 : 0;
2208            
2209 19         28 if ( TRACE_BOUNDS) {
2210             print "Poles:\tN: $include_north_pole\tS: $include_south_pole\n";
2211            
2212             print "Integer: W: $west_int, S: $south_int, E: $east_int, N: $north_int\n";
2213             print "Index ($grid_level): W: $west_idx, S: $south_idx, E: $east_idx, N: $north_idx\n";
2214             print "Shift: $shift\tgrid_size: $grid_size\tmax_grid_idx: $max_grid_idx\n";
2215             print "---\n\n";
2216             }
2217            
2218             #. Gather preliminary search results
2219            
2220 19         32 if (USE_NUMERIC_KEYS) {
2221 19         28 my $seen_n_polar = 0;
2222 19         24 my $seen_s_polar = 0;
2223 19 100       34 if ( $west <= $east ) {
2224             #. Does not straddle antimeridian
2225 17         19 if ( TRACE_BOUNDS) {
2226             print "NORMAL\n";
2227             }
2228 17         54 for (my $lat_idx = $south_idx; $lat_idx <= $north_idx; $lat_idx++) {
2229 58 100       115 if ( $lat_idx == 0 ) {
    100          
2230             #. Near south pole
2231 9 50       20 unless ( $seen_s_polar ) {
2232 9         11 $seen_s_polar = 1;
2233 9         12 my $key;
2234 9         12 if (USE_PACKED_KEYS) {
2235             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
2236             } else {
2237 9         21 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
2238             }
2239 9         13 push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2240 9         33 push @result_set, $$_points{$key};
2241             }
2242            
2243             } elsif ( $lat_idx >= $max_grid_idx ) {
2244             #. Near north pole
2245 8 50       24 unless ( $seen_n_polar ) {
2246 8         13 $seen_n_polar = 1;
2247 8         10 my $key;
2248 8         10 if (USE_PACKED_KEYS) {
2249             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
2250             } else {
2251 8         15 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
2252             }
2253 8         12 push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2254 8         24 push @result_set, $$_points{$key};
2255             }
2256            
2257             } else {
2258             #. Normal case
2259 41         89 for (my $lon_idx = $west_idx; $lon_idx <= $east_idx; $lon_idx++) {
2260 214         319 my $clipped_lon_idx = $lon_idx % $grid_size;
2261 214         230 my $key;
2262 214         232 if (USE_PACKED_KEYS) {
2263             $key = pack("Q", ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx );
2264             } else {
2265 214         284 $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
2266             }
2267 214         461 push @result_set, $$_points{$key};
2268 214         389 push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ]" if ( TRACE_BOUNDS);
2269             }
2270             }
2271            
2272             }
2273             # END does not straddle antimeridian
2274            
2275             } else { # ($west_idx > $east_idx)
2276             #. Straddles antimeridian
2277 2         5 if ( TRACE_BOUNDS) {
2278             print "STRADDLES ANTIMERIDIAN\n";
2279             }
2280 2         13 for (my $lat_idx = $south_idx; $lat_idx <= $north_idx; $lat_idx++) {
2281            
2282 9 100       24 if ( $lat_idx == 0 ) {
    100          
2283             #. Near south pole
2284 1 50       4 unless ( $seen_s_polar ) {
2285 1         1 $seen_s_polar = 1;
2286 1         4 my $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
2287 1         3 push @result_set, $$_points{$key};
2288 1         3 push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2289             }
2290            
2291             } elsif ( $lat_idx >= $max_grid_idx ) {
2292             #. Near north pole
2293 1 50       3 unless ( $seen_n_polar ) {
2294 1         2 $seen_n_polar = 1;
2295 1         3 my $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
2296 1         11 push @result_set, $$_points{$key};
2297 1         4 push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2298             }
2299            
2300             } else {
2301             #. Non-polar
2302            
2303             #. East side
2304 7         31 for (my $lon_idx = $west_idx; $lon_idx <= $max_grid_idx; $lon_idx++) {
2305 22         32 my $clipped_lon_idx = $lon_idx % $grid_size;
2306 22         31 my $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
2307 22         75 push @result_set, $$_points{$key};
2308 22         43 push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ]" if ( TRACE_BOUNDS);
2309             }
2310            
2311             #. West side
2312 7         22 for (my $lon_idx = 0; $lon_idx < $east_idx; $lon_idx++) {
2313 2         3 my $clipped_lon_idx = $lon_idx % $grid_size;
2314 2         4 my $key = ( $grid_level << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
2315 2         5 push @result_set, $$_points{$key};
2316 2         6 push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ]" if ( TRACE_BOUNDS);
2317             }
2318             }
2319             }
2320             } # END straddles antimeridian
2321            
2322             } else {
2323             #. Use split keys
2324            
2325             my $seen_n_polar = 0;
2326             my $seen_s_polar = 0;
2327             if ($west <= $east) {
2328             #. Does not straddle antimeridian
2329             if ( TRACE_BOUNDS) {
2330             print "NORMAL\n";
2331             }
2332             for (my $lat_idx = $south_idx; $lat_idx <= $north_idx; $lat_idx++) {
2333            
2334             if ( $lat_idx == 0 ) {
2335             #. Near south pole
2336             unless ( $seen_s_polar ) {
2337             $seen_s_polar = 1;
2338             my $key = [ $grid_level, $lat_idx, ALL ];
2339             push @result_set, $self->GetValue($key);
2340             push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2341             }
2342            
2343             } elsif ( $lat_idx >= $max_grid_idx ) {
2344             #. Near north pole
2345             unless ( $seen_n_polar ) {
2346             $seen_n_polar = 1;
2347             my $key = [ $grid_level, $lat_idx, ALL ];
2348             push @result_set, $self->GetValue($key);
2349             push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2350             }
2351            
2352             } else {
2353             #. Normal case
2354             for (my $lon_idx = $west_idx; $lon_idx <= $east_idx; $lon_idx++) {
2355             my $clipped_lon_idx = $lon_idx % $grid_size;
2356             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
2357             push @result_set, $self->GetValue($key);
2358             push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ($lon_idx) ]" if ( TRACE_BOUNDS);
2359             }
2360             }
2361            
2362             }
2363             # END does not straddle antimeridian
2364            
2365             } else { # ($west_idx > $east_idx)
2366             #. Straddles antimeridian
2367             if ( TRACE_BOUNDS) {
2368             print "STRADDLES ANTIMERIDIAN\n";
2369             }
2370             for (my $lat_idx = $south_idx; $lat_idx <= $north_idx; $lat_idx++) {
2371            
2372             if ( $lat_idx == 0 ) {
2373             #. Near south pole
2374             unless ( $seen_s_polar ) {
2375             $seen_s_polar = 1;
2376             my $key = [ $grid_level, $lat_idx, ALL ];
2377             push @result_set, $self->GetValue($key);
2378             push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2379             }
2380            
2381             } elsif ( $lat_idx >= $max_grid_idx ) {
2382             #. Near north pole
2383             unless ( $seen_n_polar ) {
2384             $seen_n_polar = 1;
2385             my $key = [ $grid_level, $lat_idx, ALL ];
2386             push @result_set, $self->GetValue($key);
2387             push @keys, "[ $grid_level, $lat_idx, ALL ]" if ( TRACE_BOUNDS);
2388             }
2389            
2390             } else {
2391             #. Non-polar
2392            
2393             #. East side
2394             for (my $lon_idx = $west_idx; $lon_idx <= $max_grid_idx; $lon_idx++) {
2395             my $clipped_lon_idx = $lon_idx % $grid_size;
2396             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
2397             push @result_set, $self->GetValue($key);
2398             push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ($lon_idx, E) ]" if ( TRACE_BOUNDS);
2399             }
2400            
2401             #. West side
2402             for (my $lon_idx = 0; $lon_idx < $east_idx; $lon_idx++) {
2403             my $clipped_lon_idx = $lon_idx % $grid_size;
2404             my $key = [ $grid_level, $lat_idx, $clipped_lon_idx ];
2405             push @result_set, $self->GetValue($key);
2406             push @keys, "[ $grid_level, $lat_idx, $clipped_lon_idx ($lon_idx, W) ]" if ( TRACE_BOUNDS);
2407             }
2408             }
2409             }
2410             } # END straddles antimeridian
2411            
2412             } # END using split keys
2413            
2414 19         59 print join("\n", @keys); print "\n" if ( TRACE_BOUNDS);
  19         28  
2415            
2416 19 100       37 if ( $quick_results ) {
2417             #. Return preliminary results
2418             #. Format is a list of lists (some of which may be undef
2419             #. All points within the search radius will be returned
2420             #. possibly along with additional points outside the
2421             #. search radius.
2422             return ( wantarray )
2423             #. Return array:
2424 1 50       12 ? @result_set
2425             #. Return array reference:
2426             : \@result_set;
2427             }
2428            
2429            
2430             #. Gather results
2431            
2432 18 100       43 if ($west <= $east) {
2433             #. Normal case
2434 16         21 if ( TRACE_BOUNDS) {
2435             print "NORMAL\n";
2436             }
2437            
2438 16 100       34 if (defined $condition) {
2439             #. Filter specified
2440 1         3 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
2441 1         4 foreach my $_set ( @result_set ) {
2442 12 100       22 next unless (defined $_set);
2443 9         15 foreach my $_point (@$_set) {
2444 49         67 my $p_lat = $$_point{lat};
2445 49         55 my $p_lon = $$_point{lon};
2446            
2447 49         56 if ( TRACE_BOUNDS) {
2448             print "--------------------------------------------------------------------------------\n";
2449             print " ( p_lat:$p_lat >= south:$south ) &&\n";
2450             print " ( p_lat:$p_lat <= north:$north ) &&\n";
2451             print " (\n";
2452             print " ( p_lon:$p_lon >= west:$west ) ||\n";
2453             print " ( p_lon:$p_lon <= east:$east ) ||\n";
2454             print " ( include_south_pole:$include_south_pole && p_lat:$p_lat == -90.0 ) ||\n";
2455             print " ( include_north_pole:$include_north_pole && p_lat:$p_lat == 90.0 )\n";
2456             print " )\n";
2457             print "--------------------------------------------------------------------------------\n";
2458             }
2459            
2460 49 50 66     244 if (
      66        
      66        
      33        
      66        
      33        
      33        
2461             (
2462             ( $p_lat >= $south ) &&
2463             ( $p_lat <= $north ) &&
2464             ( $p_lon >= $west ) &&
2465             ( $p_lon <= $east )
2466             ) ||
2467             ( $include_south_pole && $p_lat == -90.0 ) ||
2468             ( $include_north_pole && $p_lat == 90.0 )
2469             ) {
2470 44 100       75 if ( &$condition($_point, $_bounding_box, $user_data) ) {
2471 5         25 push @$_results, $_point;
2472 5         10 print "POINT: $$_point{name}\n" if ( TRACE_BOUNDS);
2473             } else {
2474 39         169 print "-----: $$_point{name}\n" if ( TRACE_BOUNDS);
2475             }
2476             }
2477             }
2478             }
2479             } else {
2480             #. No filter
2481 15         31 foreach my $_set ( @result_set ) {
2482 217 100       346 next unless (defined $_set);
2483 108         186 foreach my $_point (@$_set) {
2484 734         944 my $p_lat = $$_point{lat};
2485 734         863 my $p_lon = $$_point{lon};
2486            
2487 734         807 if ( TRACE_BOUNDS) {
2488             print "--------------------------------------------------------------------------------\n";
2489             print " ( p_lat:$p_lat >= south:$south ) &&\n";
2490             print " ( p_lat:$p_lat <= north:$north ) &&\n";
2491             print " (\n";
2492             print " ( p_lon:$p_lon >= west:$west ) ||\n";
2493             print " ( p_lon:$p_lon <= east:$east ) ||\n";
2494             print " ( include_south_pole:$include_south_pole && p_lat:$p_lat == -90.0 ) ||\n";
2495             print " ( include_north_pole:$include_north_pole && p_lat:$p_lat == 90.0 )\n";
2496             print " )\n";
2497             print "--------------------------------------------------------------------------------\n";
2498             }
2499            
2500 734 100 100     3343 if (
      100        
      100        
      66        
      66        
      66        
      33        
2501             (
2502             ( $p_lat >= $south ) &&
2503             ( $p_lat <= $north ) &&
2504             ( $p_lon >= $west ) &&
2505             ( $p_lon <= $east )
2506             ) ||
2507             ( $include_south_pole && $p_lat == -90.0 ) ||
2508             ( $include_north_pole && $p_lat == 90.0 )
2509             ) {
2510 670         911 push @$_results, $_point;
2511 670         1004 print "POINT: $$_point{name}\n" if ( TRACE_BOUNDS);
2512             } else {
2513 64         101 print "-----: $$_point{name}\n" if ( TRACE_BOUNDS);
2514             }
2515             }
2516             }
2517             }
2518            
2519             } else {
2520             #. Straddles antimeridian (west > east)
2521 2         4 if ( TRACE_BOUNDS) {
2522             print "STRADDLES ANTIMERIDIAN\n";
2523             }
2524            
2525 2 50       6 if (defined $condition) {
2526             #. Filter specified
2527 0         0 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
2528 0         0 foreach my $_set ( @result_set ) {
2529 0 0       0 next unless (defined $_set);
2530 0         0 foreach my $_point (@$_set) {
2531 0         0 my $p_lat = $$_point{lat};
2532 0         0 my $p_lon = $$_point{lon};
2533            
2534 0         0 if ( TRACE_BOUNDS) {
2535             print "--------------------------------------------------------------------------------\n";
2536             print " ( p_lat:$p_lat >= south:$south ) &&\n";
2537             print " ( p_lat:$p_lat <= north:$north ) &&\n";
2538             print " (\n";
2539             print " ( p_lon:$p_lon >= west:$west ) ||\n";
2540             print " ( p_lon:$p_lon <= east:$east ) ||\n";
2541             print " ( include_south_pole:$include_south_pole && p_lat:$p_lat == -90.0 ) ||\n";
2542             print " ( include_north_pole:$include_north_pole && p_lat:$p_lat == 90.0 )\n";
2543             print " )\n";
2544             print "--------------------------------------------------------------------------------\n";
2545             }
2546            
2547 0 0 0     0 if (
      0        
      0        
2548             ( $p_lat >= $south ) &&
2549             ( $p_lat <= $north ) &&
2550             (
2551             ( $p_lon >= $west ) ||
2552             ( $p_lon <= $east ) ||
2553             ( $include_south_pole && $p_lat == -90.0 ) ||
2554             ( $include_north_pole && $p_lat == 90.0 )
2555             )
2556             ) {
2557 0 0       0 if ( &$condition($_point, $_bounding_box, $user_data) ) {
2558 0         0 push @$_results, $_point;
2559 0         0 print "POINT: $$_point{name}\n" if ( TRACE_BOUNDS);
2560             } else {
2561 0         0 print "-----: $$_point{name}\n" if ( TRACE_BOUNDS);
2562             }
2563             }
2564             }
2565             }
2566             } else {
2567             #. No filter
2568 2         6 foreach my $_set ( @result_set ) {
2569 26 100       47 next unless (defined $_set);
2570 12         22 foreach my $_point (@$_set) {
2571 20         29 my $p_lat = $$_point{lat};
2572 20         27 my $p_lon = $$_point{lon};
2573            
2574 20         30 if ( TRACE_BOUNDS) {
2575             print "--------------------------------------------------------------------------------\n";
2576             print " ( p_lat:$p_lat >= south:$south ) &&\n";
2577             print " ( p_lat:$p_lat <= north:$north ) &&\n";
2578             print " (\n";
2579             print " ( p_lon:$p_lon >= west:$west ) ||\n";
2580             print " ( p_lon:$p_lon <= east:$east ) ||\n";
2581             print " ( include_south_pole:$include_south_pole && p_lat:$p_lat == -90.0 ) ||\n";
2582             print " ( include_north_pole:$include_north_pole && p_lat:$p_lat == 90.0 )\n";
2583             print " )\n";
2584             print "--------------------------------------------------------------------------------\n";
2585             }
2586            
2587 20 100 66     137 if (
      100        
      100        
2588             ( $p_lat >= $south ) &&
2589             ( $p_lat <= $north ) &&
2590             (
2591             ( $p_lon >= $west ) ||
2592             ( $p_lon <= $east ) ||
2593             ( $include_south_pole && $p_lat == -90.0 ) ||
2594             ( $include_north_pole && $p_lat == 90.0 )
2595             )
2596             ) {
2597 11         20 push @$_results, $_point;
2598 11         24 print "POINT: $$_point{name}\n" if ( TRACE_BOUNDS);
2599             } else {
2600 9         18 print "-----: $$_point{name}\n" if ( TRACE_BOUNDS);
2601             }
2602             }
2603             }
2604             }
2605            
2606             }
2607            
2608 18 100       39 if ( $max_results ) {
2609 1         3 my $results_count = scalar @$_results;
2610 1 50       5 $max_results = $results_count if ($results_count < $max_results);
2611             }
2612            
2613             return ( wantarray )
2614             #. Return array:
2615 18 0       181 ? ( scalar @$_results )
    50          
    50          
    50          
    100          
    50          
    100          
2616             ? ( $max_results )
2617             ? ( scalar @$_results > $max_results )
2618             ? @$_results[0..($max_results-1)]
2619             : @$_results
2620             : @$_results
2621             : ()
2622             #. Return array reference:
2623             : ( scalar @$_results )
2624             ? ( $max_results )
2625             ? ( scalar @$_results > $max_results )
2626             ? [ @$_results[0..($max_results-1)] ]
2627             : $_results
2628             : $_results
2629             : undef; #. undef == No result found
2630             }
2631              
2632              
2633             # ==============================================================================
2634              
2635              
2636              
2637              
2638              
2639              
2640              
2641              
2642              
2643              
2644              
2645             =head2 Closest( ... )
2646              
2647             =over
2648              
2649             C<@results = $index-EClosest( \%point, $number_of_points_desired, \%options );>
2650              
2651             C<$results = $index-EClosest( \%point, $number_of_points_desired, \%options );>
2652              
2653             Find the point or points closest to a given point
2654              
2655             Note that if you want to find the closest points within a given radius it may be
2656             faster to use C> instead. See the B>
2657             section below for more details.
2658              
2659             B>
2660              
2661             =over
2662              
2663             The point to search near
2664              
2665             This is either a reference to a hash containing at a minimum a C and a
2666             C value (both in degrees) or a reference to an array giving the point.
2667             See the B> section above for details.
2668              
2669             =back
2670              
2671             B>
2672              
2673             =over
2674              
2675             The number of points that should be returned.
2676              
2677             Set to C<0> to not restrict the number of points returned or set it S C<0>>
2678             to set the maximum number of points to return.
2679              
2680             If omitted then this will default to C<1>.
2681              
2682             =back
2683              
2684             B>
2685              
2686             =over
2687              
2688             The parameters for the search (all are optional):
2689              
2690             B>
2691              
2692             =over
2693              
2694             Only return results within this distance (in meters) from search point.
2695              
2696             If no C is specified or the C is set to C then
2697             all points in the index may potentially be returned.
2698              
2699             =back
2700              
2701             B>
2702              
2703             =over
2704              
2705             Sort results by distance from point
2706              
2707             By default points returned are sorted by distance. Set this to C<0> to not
2708             sort the returned points.
2709              
2710             Although sorting is not mandatory, performing it is strongly recommended since
2711             otherwise the set of points returned are not guaranteed to be the closest.
2712              
2713             =back
2714              
2715             B>
2716              
2717             =over
2718              
2719             Reference to additional user-supplied code to determine whether each point
2720             should be included in the results.
2721              
2722             This code is run before the distance from the search point to the result point
2723             has been calculated.
2724              
2725             See the B> section below for syntax.
2726              
2727             =back
2728              
2729             B>
2730              
2731             =over
2732              
2733             Reference to additional user-supplied code to determine whether each point
2734             should be included in the results.
2735              
2736             This code is run after the distance from the search point to the result point
2737             has been calculated.
2738              
2739             By default, a C function that filters out the search point
2740             is used. To remove this default function either specify a new one or set
2741             C to "C".
2742              
2743             See the B> section below for syntax.
2744              
2745             =back
2746              
2747             B>
2748              
2749             =over
2750              
2751             Arbitrary user-supplied data that is passed to the condition functions.
2752              
2753             This can be used to allow the function access to additional data structures.
2754              
2755             See the B> section below for syntax.
2756              
2757             =back
2758              
2759             =back
2760              
2761             B
2762              
2763             =over
2764              
2765             In list context the return value is a list of references to the points found or
2766             an empty array if none were found.
2767              
2768             In scalar context the return value is a reference to the aforementioned list or
2769             C if no results were found.
2770              
2771             For each point in the results the distance in meters from it to the search point
2772             will be stored in the C entry in the result point's
2773             hash. It can be retrieved using e.g. S>
2774              
2775             =back
2776              
2777             =back
2778              
2779             =cut
2780              
2781              
2782              
2783             #. Return closest point or points to search point
2784              
2785             sub Closest($$;$$) {
2786 16     16 1 8199 my ($self, $_search_point, $number_of_points_desired, $_options) = @_;
2787            
2788             # Allow calling as Closest( POINT, OPTIONS ) when only a single point is desired.
2789 16 100       42 if (ref $number_of_points_desired) {
2790 2         3 $_options = $number_of_points_desired;
2791 2         2 $number_of_points_desired = 1;
2792             }
2793            
2794             #. Get the point index
2795 16         35 my $_points = $$self{index};
2796            
2797            
2798 16 100       44 if (ref $_search_point eq 'ARRAY') {
2799             #. Got array; expand arguments into a full point
2800 13         25 my $lat = $$_search_point[0];
2801 13         15 my $lon = $$_search_point[1];
2802            
2803 13         41 $_search_point = { 'lat'=>$lat, 'lon'=>$lon };
2804             }
2805            
2806             #. Get search point's position
2807 16         26 my $p_lat = $$_search_point{lat};
2808 16         26 my $p_lon = $$_search_point{lon};
2809            
2810             #. Search options; user should omit (or set to undef) inactive options:
2811            
2812 16         24 my $pre_condition = $$_options{pre_condition}; #. Reference to subroutine returning true if current point should be considered as
2813             #. a possible result, false otherwise. This subroutine should not modify any data.
2814             #. This subroutine is called before the distance from the search point to the
2815             #. result point has been calculated.
2816             #.
2817 16         23 my $post_condition = $$_options{post_condition}; #. Reference to subroutine returning true if current point should be considered as
2818             #. a possible result, false otherwise. This subroutine should not modify any data.
2819             #. This subroutine is called after the distance from the search point to the
2820             #. result point has been calculated.
2821             #.
2822             #. If no post_condition is specified then the default function (given below, omits
2823             #. the search point from the results) will be used. To override this behavior
2824             #. either define your own post_condition function or set post_condition to "NONE".
2825             #.
2826             #. Default post_condition:
2827             #.
2828             #. sub {
2829             #. my ( $result_point, $search_point, $user_data ) = @_;
2830             #. return ( $result_point != $search_point );
2831             #. }
2832             #.
2833 16         22 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
2834             #.
2835 16         20 my $search_radius = $$_options{radius}; #. Only points within radius (in meters) will be considered.
2836             #. Default: No distance restriction
2837             #.
2838 16         24 my $sort_results = $$_options{sort_results}; #. Sort results by distance from point.
2839             #. Set to 0 to not sort results
2840             #. Default: Points are sorted by distance
2841            
2842             #. Maximum number of results to return.
2843             #. Set to 0 to return all matching results (use with care; specifying a radius is strongly suggested)
2844             #. Default: Only one point is returned
2845 16 100       40 $number_of_points_desired = 1 unless (defined $number_of_points_desired);
2846            
2847 16 100       42 if ( ! defined $post_condition ) {
    100          
2848 10     67   45 $post_condition = sub { my ( $result_point, $search_point, $user_data ) = @_; return ( $result_point != $_search_point ); };
  67         116  
  67         218  
2849             } elsif ($post_condition eq 'NONE') {
2850 2         10 $post_condition = undef;
2851             }
2852            
2853             #. Used to speed up inner loops:
2854            
2855 16 100       41 my $no_pre_condition = ( defined $pre_condition ) ? 0 : 1;
2856 16 100       39 my $no_post_condition = ( defined $post_condition ) ? 0 : 1;
2857            
2858 16 50       34 my $no_search_radius = ( defined $search_radius ) ? 0 : 1;
2859            
2860             #. Always sort unless explicitly told not to
2861 16 50       46 $sort_results = 1 unless ( defined $sort_results );
2862            
2863             #. Get parameters for level one past the most detailed level in the index
2864            
2865 16         28 my $max_level = $self->{max_level};
2866            
2867 16         24 my $cur_level = $max_level + 1;
2868 16         26 my $cur_size = 2**$cur_level;
2869 16         29 my $cur_max_idx = $cur_size - 1;
2870            
2871             #. Get the integer forms of the search point's latitude and longitude
2872            
2873 16         35 my $p_lat_idx = int( ( $p_lat + 90.0 ) * $cur_size / 180.0 );
2874 16 50       34 $p_lat_idx = $cur_size if ($p_lat_idx > $cur_size); # This includes the north pole
2875            
2876 16         44 my $p_lon_idx = ( int( ( $p_lon + 180.0 ) * $cur_size / 360.0 ) % $cur_size );
2877            
2878             #. Determine the low bit of the integer latitude and longitude
2879             #. This is used to determine which tile edge the search point is closest to.
2880            
2881 16         22 my $lat_low_bit = $p_lat_idx & 1;
2882 16         24 my $lon_low_bit = $p_lon_idx & 1;
2883            
2884             #. Initialize the bit shift to one past the end of the index
2885            
2886 16         20 my $shift = -1;
2887            
2888             #. Note the search point's position in radians
2889            
2890 16         19 my $p_lat_rad;
2891 16 100       35 if (defined $$_search_point{lat_rad}) {
2892 3         4 $p_lat_rad = $$_search_point{lat_rad};
2893             } else {
2894 13         34 $p_lat_rad = Math::Trig::deg2rad($p_lat);
2895 13         127 $$_search_point{lat_rad} = $p_lat_rad;
2896             }
2897            
2898 16         21 my $p_lon_rad;
2899 16 100       31 if (defined $$_search_point{lon_rad}) {
2900 3         4 $p_lon_rad = $$_search_point{lon_rad};
2901             } else {
2902 13         24 $p_lon_rad = Math::Trig::deg2rad($p_lon);
2903 13         105 $$_search_point{lon_rad} = $p_lon_rad;
2904             }
2905            
2906             #. Determine grid sizes, etc. in meters
2907            
2908 16         25 my $NS_circumference_in_meters = $self->{polar_circumference};
2909            
2910 16         43 my $lat_meter_in_degrees = 360.0 / $NS_circumference_in_meters;
2911            
2912 16         47 my $EW_circumference_in_meters = $self->LongitudeCircumference($p_lat_rad);
2913            
2914 16         34 my $lon_meter_in_degrees = 360.0 / $EW_circumference_in_meters;
2915            
2916 16         24 my $lat_grid_in_meters = ( $NS_circumference_in_meters / 2.0 ) / $cur_size;
2917            
2918 16         28 my $lon_grid_in_meters = $EW_circumference_in_meters / $cur_size;
2919            
2920 16         27 my $p_lat_m = $NS_circumference_in_meters / 360.0 * ( $p_lat + 90.0 );
2921 16         30 my $p_lon_m = $EW_circumference_in_meters / 360.0 * ( $p_lon + 180.0 );
2922            
2923             #. Initialize the scratch pads and results
2924            
2925 16         21 my %distances = (); #. Distance to each point in current set of results
2926 16         21 my %considered = (); #. Points that have been considered
2927             #. (used to skip results that are already in @valid)
2928 16         28 my @valid = (); #. Points that meet the search criteria
2929            
2930             #. Set the distance origin to the search point
2931 16         74 SetUpDistance($self->{planetary_diameter}, $p_lat_rad, $p_lon_rad);
2932            
2933 16         23 my $_results; #. Holds the points found in the most recent result tile
2934            
2935             #. Used to exit loop early when a search radius is specified
2936 16         17 my $largest_distance_seen = 0;
2937            
2938             #. Loop through zoom levels from most zoomed in to least...
2939 16         40 while ( $cur_level >= 0 ) {
2940            
2941 267         359 my $adj_lat_idx; #. Adjacent grid index (latitude)
2942             my $adj_lon_idx;
2943            
2944 267         344 my $valid_radius_m = 0; #. Distance from search point to closest edge of grid tiles
2945            
2946             #. Set up for current grid level
2947 267         335 $cur_level--; #. Zoom out one level
2948 267         318 $cur_size >>= 1; #. Grid is now half the size
2949 267         336 my $max_grid_idx = $cur_size - 1;
2950 267         333 $lat_grid_in_meters *= 2.0; #. Grid tiles now have twice the height
2951 267         306 $lon_grid_in_meters *= 2.0; #. Grid tiles now have twice the width
2952            
2953 267 100       395 if ($cur_level < 0) {
2954             #. World-wide
2955            
2956 9         13 if (USE_NUMERIC_KEYS) {
2957            
2958 9         12 my $key;
2959 9         10 if (USE_PACKED_KEYS) {
2960             $key = pack("Q", ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON );
2961             } else {
2962 9         12 $key = ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON;
2963             }
2964 9         17 $_results = $$_points{$key};
2965 9 50       19 next unless defined $_results;
2966            
2967 9         14 foreach my $_point ( @$_results ) {
2968 99 100       193 next if ($considered{$_point});
2969            
2970 81 100 100     180 if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
2971            
2972 73         111 my $distance;
2973 73 100       127 if (defined $distances{$_point}) {
2974 31         45 $distance = $distances{$_point};
2975             } else {
2976 42         124 $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
2977 42         94 $distances{$_point} = $distance;
2978 42         62 $$_point{search_result_distance} = $distance;
2979             }
2980            
2981 73         118 $considered{$_point} = 1;
2982 73 50 33     136 if ( $no_search_radius || ($distance <= $search_radius) ) {
2983 73 50 33     137 if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
2984 73         148 push @valid, $_point;
2985             }
2986             }
2987             } else { #. Pre-condition failed
2988 8         45 $considered{$_point} = 1;
2989             }
2990             }
2991            
2992             } else {
2993             #. Use split keys
2994            
2995             my $key = [ ALL, ALL, ALL ];
2996             $_results = $self->GetValue($key);
2997             next unless defined $_results;
2998            
2999             foreach my $_point ( @$_results ) {
3000             next if ($considered{$_point});
3001            
3002             if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3003            
3004             my $distance;
3005             if (defined $distances{$_point}) {
3006             $distance = $distances{$_point};
3007             } else {
3008             $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3009             $distances{$_point} = $distance;
3010             $$_point{search_result_distance} = $distance;
3011             }
3012            
3013             $considered{$_point} = 1;
3014             if ( $no_search_radius || ($distance <= $search_radius) ) {
3015             if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3016             push @valid, $_point;
3017             }
3018             }
3019             } else { #. Pre-condition failed
3020             $considered{$_point} = 1;
3021             }
3022             }
3023            
3024             } # END if split keys
3025            
3026             } else {
3027             #. Normal case
3028            
3029             #. Determine which grid tile edge the point is closest to
3030             #. This is done by looking at LSB of the indices of the previous (more zoomed in) level
3031 258         328 my $lat_low_bit = $p_lat_idx & 1; #. Latitude previous LSB
3032 258         316 my $lon_low_bit = $p_lon_idx & 1; #. Longitude previous LSB
3033            
3034 258         296 $p_lat_idx >>= 1; #. Latitude index is now half the previous value
3035 258         288 $p_lon_idx >>= 1; #. Longitude index is now half the previous value
3036            
3037 258         368 my $valid_radius_lat_m;
3038             my $valid_radius_lon_m;
3039            
3040 258 100       381 if ($lat_low_bit) {
3041             #. Closer to top
3042            
3043             #. Get adjacent tile's latitude grid index
3044 70         89 $adj_lat_idx = $p_lat_idx + 1;
3045            
3046 70 100       120 if ( $adj_lat_idx == $cur_size ) {
3047             #. We're abutting the the north pole
3048 10         48 $adj_lat_idx = undef;
3049             }
3050            
3051             #. Current search radius is the distance to bottom edge of point's grid tile
3052 70         98 my $lower_edge_m = $p_lat_idx * $lat_grid_in_meters;
3053 70         98 $valid_radius_lat_m = $p_lat_m - $lower_edge_m;
3054            
3055             } else {
3056             #. Closer to bottom
3057            
3058             #. Get adjacent tile's latitude grid index
3059 188         223 $adj_lat_idx = $p_lat_idx - 1;
3060            
3061 188 100       293 if ($adj_lat_idx < 0) {
3062             #. South polar
3063 76         109 $adj_lat_idx = undef;
3064             }
3065             #. Current search radius is the distance to upper edge of point's grid tile
3066 188         272 my $upper_edge_m = ($p_lat_idx + 1) * $lat_grid_in_meters;
3067 188         243 $valid_radius_lat_m = $upper_edge_m - $p_lat_m;
3068             }
3069            
3070 258 100       375 if ($lon_low_bit) {
3071             #. Closer to right
3072 71         90 $adj_lon_idx = ( $p_lon_idx + 1 ) % $cur_size;
3073            
3074             #. Current search radius is the distance to left edge of point's grid tile
3075 71         93 my $left_edge_m = $p_lon_idx * $lon_grid_in_meters;
3076 71         93 $valid_radius_lon_m = $p_lon_m - $left_edge_m;
3077            
3078             } else {
3079             #. Closer to left
3080 187         241 $adj_lon_idx = ( $p_lon_idx - 1 ) % $cur_size;
3081            
3082             #. Current search radius is the distance to right edge of point's grid tile
3083 187         249 my $right_edge_m = ($p_lon_idx + 1) * $lon_grid_in_meters;
3084 187         224 $valid_radius_lon_m = $right_edge_m - $p_lon_m;
3085            
3086             }
3087            
3088 258 50       481 $valid_radius_m = ( $valid_radius_lat_m < $valid_radius_lon_m ) ? $valid_radius_lat_m : $valid_radius_lon_m;
3089             }
3090            
3091             #. Oddly it's actually slightly faster to NOT split this code into four versions
3092             #. (pre- and post-condition, pre- only, post- only and no conditions) and instead
3093             #. do the checks inline as coded below.
3094            
3095 267         308 if (USE_NUMERIC_KEYS) {
3096            
3097 267         380 foreach my $lat_idx ( $p_lat_idx, $adj_lat_idx) {
3098 534 100       912 next unless defined $lat_idx;
3099            
3100 439 100       719 if ( $lat_idx == 0 ) {
    100          
3101             #. Near south pole
3102 91         109 my $key;
3103 91         101 if (USE_PACKED_KEYS) {
3104             $key = pack("Q", ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
3105             } else {
3106 91         140 $key = ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
3107             }
3108 91         158 $_results = $$_points{$key};
3109 91 100       166 next unless defined $_results;
3110            
3111 75         107 foreach my $_point ( @$_results ) {
3112 100 100       203 next if ($considered{$_point});
3113            
3114 75 100 100     146 if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3115            
3116 69         133 my $distance;
3117 69 100       127 if (defined $distances{$_point}) {
3118 45         77 $distance = $distances{$_point};
3119             } else {
3120 24         97 $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3121 24         46 $distances{$_point} = $distance;
3122 24         44 $$_point{search_result_distance} = $distance;
3123             }
3124            
3125 69 100       147 if ($distance <= $valid_radius_m) {
3126 9         18 $considered{$_point} = 1;
3127 9 50 33     27 $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3128 9 50 33     20 if ( $no_search_radius || ($distance <= $search_radius) ) {
3129 9 50 66     25 if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3130 9         27 push @valid, $_point;
3131             }
3132             }
3133             }
3134            
3135             } else { #. Pre-condition failed
3136 6         75 $considered{$_point} = 1;
3137             }
3138             }
3139            
3140             } elsif ( $lat_idx >= $max_grid_idx ) {
3141             #. Near north pole
3142 130         148 my $key;
3143 130         141 if (USE_PACKED_KEYS) {
3144             $key = pack("Q", ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON );
3145             } else {
3146 130         221 $key = ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | MASK_LATLON ;
3147             }
3148 130         231 $_results = $$_points{$key};
3149 130 100       228 next unless defined $_results;
3150            
3151 68         110 foreach my $_point ( @$_results ) {
3152 115 100       243 next if ($considered{$_point});
3153            
3154 86 50 33     178 if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3155            
3156 86         100 my $distance;
3157 86 100       153 if (defined $distances{$_point}) {
3158 58         97 $distance = $distances{$_point};
3159             } else {
3160 28         86 $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3161 28         58 $distances{$_point} = $distance;
3162 28         42 $$_point{search_result_distance} = $distance;
3163             }
3164            
3165 86 100       185 if ($distance <= $valid_radius_m) {
3166 8         16 $considered{$_point} = 1;
3167 8 50 33     20 $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3168 8 50 33     23 if ( $no_search_radius || ($distance <= $search_radius) ) {
3169            
3170 8 50 33     24 if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3171 8         24 push @valid, $_point;
3172             }
3173             }
3174             }
3175            
3176             } else { #. Pre-condition failed
3177 0         0 $considered{$_point} = 1;
3178             }
3179             }
3180            
3181             } else {
3182             #. Normal case
3183 218         291 foreach my $lon_idx ( $p_lon_idx, $adj_lon_idx ) {
3184 436         571 my $clipped_lon_idx = $lon_idx % $cur_size;
3185 436         511 my $key;
3186 436         478 if (USE_PACKED_KEYS) {
3187             $key = pack("Q", ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx );
3188             } else {
3189 436         609 $key = ( ($cur_level-1) << 58 ) | ( $lat_idx << 29 ) | $clipped_lon_idx ;
3190             }
3191 436         644 $_results = $$_points{$key};
3192 436 100       798 next unless defined $_results;
3193            
3194 21         29 foreach my $_point ( @$_results ) {
3195 21 100       57 next if ($considered{$_point});
3196            
3197 4 50 33     14 if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3198            
3199 4         5 my $distance;
3200 4 50       16 if (defined $distances{$_point}) {
3201 0         0 $distance = $distances{$_point};
3202             } else {
3203 4         24 $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3204 4         9 $distances{$_point} = $distance;
3205 4         10 $$_point{search_result_distance} = $distance;
3206             }
3207            
3208 4 100       12 if ($distance <= $valid_radius_m) {
3209 3         16 $considered{$_point} = 1;
3210 3 50 33     25 $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3211 3 50 33     10 if ( $no_search_radius || ($distance <= $search_radius) ) {
3212 3 100 66     49 if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3213 2         7 push @valid, $_point;
3214             }
3215             }
3216             }
3217             } else { #. Pre-condition failed
3218 0         0 $considered{$_point} = 1;
3219             }
3220             }
3221             }
3222             }
3223            
3224             }
3225            
3226             } else {
3227             #. Use split keys
3228            
3229             foreach my $lat_idx ( $p_lat_idx, $adj_lat_idx) {
3230             next unless defined $lat_idx;
3231            
3232             if ( $lat_idx == 0 ) {
3233             #. Near south pole
3234             my $key = [ $cur_level-1, $lat_idx, ALL ];
3235             $_results = $self->GetValue($key);
3236             next unless defined $_results;
3237            
3238             foreach my $_point ( @$_results ) {
3239             next if ($considered{$_point});
3240            
3241             if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3242            
3243             my $distance;
3244             if (defined $distances{$_point}) {
3245             $distance = $distances{$_point};
3246             } else {
3247             $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3248             $distances{$_point} = $distance;
3249             $$_point{search_result_distance} = $distance;
3250             }
3251            
3252             if ($distance <= $valid_radius_m) {
3253             $considered{$_point} = 1;
3254             $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3255             if ( $no_search_radius || ($distance <= $search_radius) ) {
3256             if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3257             push @valid, $_point;
3258             }
3259             }
3260             }
3261            
3262             } else { #. Pre-condition failed
3263             $considered{$_point} = 1;
3264             }
3265             }
3266            
3267             } elsif ( $lat_idx >= $max_grid_idx ) {
3268             #. Near north pole
3269             my $key = [ $cur_level-1, $lat_idx, ALL ];
3270             $_results = $self->GetValue($key);
3271             next unless defined $_results;
3272            
3273             foreach my $_point ( @$_results ) {
3274             next if ($considered{$_point});
3275            
3276             if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3277            
3278             my $distance;
3279             if (defined $distances{$_point}) {
3280             $distance = $distances{$_point};
3281             } else {
3282             $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3283             $distances{$_point} = $distance;
3284             $$_point{search_result_distance} = $distance;
3285             }
3286            
3287             if ($distance <= $valid_radius_m) {
3288             $considered{$_point} = 1;
3289             $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3290             if ( $no_search_radius || ($distance <= $search_radius) ) {
3291             if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3292             push @valid, $_point;
3293             }
3294             }
3295             }
3296            
3297             } else { #. Pre-condition failed
3298             $considered{$_point} = 1;
3299             }
3300             }
3301            
3302             } else {
3303             #. Normal case
3304             foreach my $lon_idx ( $p_lon_idx, $adj_lon_idx ) {
3305             my $clipped_lon_idx = $lon_idx % $cur_size;
3306             my $key = [ $cur_level-1, $lat_idx, $clipped_lon_idx ];
3307             $_results = $self->GetValue($key);
3308             next unless defined $_results;
3309            
3310             foreach my $_point ( @$_results ) {
3311             next if ($considered{$_point});
3312            
3313             if ( $no_pre_condition || &$pre_condition($_point, $_search_point, $user_data) ) {
3314            
3315             my $distance;
3316             if (defined $distances{$_point}) {
3317             $distance = $distances{$_point};
3318             } else {
3319             $distance = HaversineDistance($$_point{lat_rad}, $$_point{lon_rad});
3320             $distances{$_point} = $distance;
3321             $$_point{search_result_distance} = $distance;
3322             }
3323            
3324             if ($distance <= $valid_radius_m) {
3325             $considered{$_point} = 1;
3326             $largest_distance_seen = $distance unless ( $no_search_radius || ($distance < $largest_distance_seen) );
3327             if ( $no_search_radius || ($distance <= $search_radius) ) {
3328             if ( $no_post_condition || &$post_condition($_point, $_search_point, $user_data) ) {
3329             push @valid, $_point;
3330             }
3331             }
3332             }
3333             } else { #. Pre-condition failed
3334             $considered{$_point} = 1;
3335             }
3336             }
3337             }
3338             }
3339            
3340             }
3341            
3342             } # END if split keys
3343            
3344             #. Stop searching if we have found sufficient points
3345 267 100 100     576 last if ( ($number_of_points_desired) && (scalar @valid >= $number_of_points_desired) );
3346            
3347             #. Stop searching if search radius exceeded
3348 259 50 33     601 last unless ( $no_search_radius || ( $largest_distance_seen < $search_radius) );
3349             } # END loop through levels
3350            
3351             #. Sort results by distance
3352 16 50       65 @valid = sort { $$a{search_result_distance} <=> $$b{search_result_distance} } @valid if ( $sort_results );
  195         275  
3353            
3354             #. Only include requested number of points
3355 16 100       55 if ( $number_of_points_desired ) {
3356 8         16 my $count_to_return = scalar @valid;
3357 8 100       18 $count_to_return = $number_of_points_desired if ($number_of_points_desired < $count_to_return);
3358 8         40 @valid = @valid[0 .. $count_to_return-1];
3359             }
3360            
3361             #. Return points found
3362             return ( wantarray )
3363             #. Return array:
3364 16 100       161 ? @valid
3365             #. Return array reference:
3366             : \@valid;
3367             }
3368              
3369              
3370              
3371             # ==============================================================================
3372              
3373              
3374              
3375              
3376              
3377              
3378              
3379              
3380              
3381              
3382              
3383             =head2 Farthest( ... )
3384              
3385             =over
3386              
3387             C<@results = $index-EFarthest( \%point, $number_of_points_desired, \%options );>
3388              
3389             C<$results = $index-EFarthest( \%point, $number_of_points_desired, \%options );>
3390              
3391             Find the point or points farthest from a given point
3392              
3393             In other words, find the points closest to a given point's antipode.
3394              
3395             B>
3396              
3397             =over
3398              
3399             The point to search relative to
3400              
3401             This is either a reference to a hash containing at a minimum a C and a
3402             C value (both in degrees) or a reference to an array giving the point.
3403             See the B> section above for details.
3404              
3405             =back
3406              
3407             B>
3408              
3409             =over
3410              
3411             The number of points that should be returned.
3412              
3413             Set to C<0> to not restrict the number of points returned or set it SC<0>>
3414             to set the maximum number of points to return.
3415              
3416             If omitted then this will default to C<1>.
3417              
3418             =back
3419              
3420             B>
3421              
3422             =over
3423              
3424             The parameters for the search (all are optional):
3425              
3426             B>
3427              
3428             =over
3429              
3430             Only return results within this distance (in meters) from search point.
3431              
3432             If no C is specified or the C is set to C then
3433             all points in the index may potentially be returned.
3434              
3435             =back
3436              
3437             B>
3438              
3439             =over
3440              
3441             Sort results by distance from point
3442              
3443             By default points returned are sorted by distance. Set this to C<0> to not
3444             sort the returned points.
3445              
3446             Although sorting is not mandatory, performing it is strongly recommended since
3447             otherwise the set of points returned are not guaranteed to be the farthest.
3448              
3449             =back
3450              
3451             B>
3452              
3453             =over
3454              
3455             Reference to additional user-supplied code to determine whether each point
3456             should be included in the results.
3457              
3458             This code is run before the distance from the search point to the result point
3459             has been calculated.
3460              
3461             See the B> section below for syntax.
3462              
3463             =back
3464              
3465             B>
3466              
3467             =over
3468              
3469             Reference to additional user-supplied code to determine whether each point
3470             should be included in the results.
3471              
3472             This code is run after the distance from the search point to the result point
3473             has been calculated.
3474              
3475             By default, a C function that filters out the search point is
3476             used. To remove this default function either specify a new one, set a value for
3477             C, or set C to "C".
3478              
3479             See the B> section below for syntax.
3480              
3481             =back
3482              
3483             B>
3484              
3485             =over
3486              
3487             Arbitrary user-supplied data that is passed to the condition functions.
3488              
3489             This can be used to allow the function access to additional data structures.
3490              
3491             If the default C is active and no C value has been
3492             provided by the caller then this is set to the actual (non-antipodal) search
3493             point.
3494              
3495             See the B> section below for syntax.
3496              
3497             =back
3498              
3499             =back
3500              
3501             B
3502              
3503             =over
3504              
3505             In list context the return value is a list of references to the points found or
3506             an empty array if none were found.
3507              
3508             In scalar context the return value is a reference to the aforementioned list or
3509             C if no results were found.
3510              
3511             For each point in the results the distance in meters from it to the search point
3512             will be stored in the C entry in the result point's
3513             hash. In addition, the distance from a result point to the search point's
3514             antipode will be stored in the C entry. These can be
3515             retrieved using e.g.:
3516              
3517             $meters_from_search_point = $$point{search_result_distance};
3518             $meters_to_antipodal_point = $$point{antipode_distance};
3519              
3520             =back
3521              
3522             =back
3523              
3524             =cut
3525              
3526              
3527              
3528             #. Return farthest point or points from search point
3529              
3530             sub Farthest($$;$$) {
3531 4     4 1 4422 my ($self, $_search_point, $number_of_points_desired, $_options) = @_;
3532            
3533             # Allow calling as Farthest( POINT, OPTIONS ) when only a single point is desired.
3534 4 50       15 if (ref $number_of_points_desired) {
3535 0         0 $_options = $number_of_points_desired;
3536 0         0 $number_of_points_desired = 1;
3537             }
3538            
3539             #. Get the point index
3540 4         8 my $_points = $$self{index};
3541            
3542            
3543 4 50       11 if (ref $_search_point eq 'ARRAY') {
3544             #. Got array; expand arguments into a full point
3545 4         8 my $lat = $$_search_point[0];
3546 4         7 my $lon = $$_search_point[1];
3547            
3548 4         11 $_search_point = { 'lat'=>$lat, 'lon'=>$lon };
3549             }
3550            
3551             #. Get the point's position
3552 4         8 my $p_lat = $$_search_point{lat};
3553 4         7 my $p_lon = $$_search_point{lon};
3554            
3555             #. We'll be using the antipodal point as the center of the search
3556 4         6 my $antipode_lat = -1 * $p_lat;
3557 4         8 my $antipode_lon = $p_lon + 180.0;
3558 4 50       11 $antipode_lon -= 360.0 if ( $antipode_lon >= 180.0 );
3559            
3560             #. Search options; user should omit (or set to undef) inactive options:
3561            
3562 4         7 my $pre_condition = $$_options{pre_condition}; #. Reference to subroutine returning true if current point should be considered as
3563             #. a possible result, false otherwise. This subroutine should not modify any data.
3564             #. This subroutine is called before the distance from the search point to the
3565             #. result point has been calculated.
3566             #.
3567 4         15 my $post_condition = $$_options{post_condition}; #. Reference to subroutine returning true if current point should be considered as
3568             #. a possible result, false otherwise. This subroutine should not modify any data.
3569             #. This subroutine is called after the distance from the search point to the
3570             #. result point has been calculated.
3571             #.
3572             #. If no post_condition and no user_data is specified then the default function
3573             #. (given below, omits the search point from the results) will be used. To over-
3574             #. ride this behavior either define your own post_condition function, specify your
3575             #. own user_data, or set post_condition to "NONE".
3576             #.
3577             #. Default post_condition:
3578             #.
3579             #. sub {
3580             #. my ( $result_point, $search_point, $user_data ) = @_;
3581             #. return ( $result_point != $user_data );
3582             #. }
3583             #.
3584 4         8 my $user_data = $$_options{user_data}; #. User-defined data that is passed on to the condition subroutine.
3585             #.
3586             #. Default user_data: The actual (non-antipodal) search point
3587             #.
3588 4         8 my $search_radius = $$_options{radius}; #. Only points within radius (in meters) will be considered.
3589             #. Default: No distance restriction
3590             #.
3591 4         7 my $sort_results = $$_options{sort_results}; #. Sort results by distance from point.
3592             #. Set to 0 to not sort results
3593             #. Default: Points are sorted by distance
3594            
3595 4 50 33     19 if ( ( ! defined $post_condition ) && ( ! defined $user_data ) ) {
    0          
3596 4     24   18 $post_condition = sub { my ( $result_point, $search_point, $user_data ) = @_; return ( $result_point != $user_data ); };
  24         43  
  24         86  
3597 4         9 $user_data = $_search_point;
3598            
3599             } elsif ($post_condition eq 'NONE') {
3600 0         0 $post_condition = undef;
3601             }
3602            
3603             #. Maximum number of results to return.
3604             #. Set to 0 to return all matching results (use with care; specifying a radius is strongly suggested)
3605             #. Default: Only one point is returned
3606 4 100       16 $number_of_points_desired = 1 unless (defined $number_of_points_desired);
3607            
3608 4         16 my %options = (
3609             pre_condition => $pre_condition,
3610             post_condition => $post_condition,
3611             user_data => $user_data,
3612             radius => $search_radius,
3613             sort_results => $sort_results
3614             );
3615            
3616 4         14 my $_results = $self->Closest( [ $antipode_lat, $antipode_lon ], $number_of_points_desired, \%options);
3617            
3618 4         10 my $p_lat_rad;
3619 4 50       9 if (defined $$_search_point{lat_rad}) {
3620 0         0 $p_lat_rad = $$_search_point{lat_rad};
3621             } else {
3622 4         14 $p_lat_rad = Math::Trig::deg2rad($p_lat);
3623 4         40 $$_search_point{lat_rad} = $p_lat_rad;
3624             }
3625            
3626 4         7 my $p_lon_rad;
3627 4 50       12 if (defined $$_search_point{lon_rad}) {
3628 0         0 $p_lon_rad = $$_search_point{lon_rad};
3629             } else {
3630 4         21 $p_lon_rad = Math::Trig::deg2rad($p_lon);
3631 4         33 $$_search_point{lon_rad} = $p_lon_rad;
3632             }
3633            
3634 4         19 SetUpDistance($self->{planetary_diameter}, $p_lat_rad, $p_lon_rad);
3635            
3636 4         9 foreach my $_antipode (@$_results) {
3637             #. Distance found is actually the distance from the search point's antipodal point
3638 24         59 $$_antipode{antipode_distance} = $$_antipode{search_result_distance};
3639            
3640             #. Calculate and record the actual distance from the search point
3641 24         61 $$_antipode{search_result_distance} = HaversineDistance($$_antipode{lat_rad}, $$_antipode{lon_rad});
3642             }
3643            
3644             #. Return points found
3645             return ( wantarray )
3646             #. Return array:
3647 4 50       40 ? @$_results
3648             #. Return array reference:
3649             : $_results;
3650             }
3651              
3652              
3653              
3654             # ==============================================================================
3655              
3656              
3657              
3658              
3659              
3660             =head2 Distance( ... )
3661              
3662             =over
3663              
3664             C<$meters = $index-EDistance( \%point_1, \%point_2 );>
3665              
3666             C<$meters = $index-EDistance( \@point_1, \@point_2 );>
3667              
3668             C<$meters = $index-EDistance( \%point_1, \@point_2 );>
3669              
3670             C<$meters = $index-EDistance( \@point_1, \%point_2 );>
3671              
3672             Returns the distance in meters between two points
3673              
3674             The haversine function is used to compute the distance. As this assumes a
3675             spherical body the distances returned may show errors. Using the default
3676             options, these errors are up to 0.056% (north - south) or 1.12% (east - west).
3677             Such errors typically start becoming significant at distances over S<20 km>.
3678              
3679             B> or B>, B> or B>
3680              
3681             =over
3682              
3683             The points to measure the distance between
3684              
3685             These can be either hashes containing at a minimum a C and a C value
3686             (both in degrees) or arrays giving each point. See the B>
3687             section above for details.
3688              
3689             =back
3690              
3691             =back
3692              
3693             =cut
3694              
3695              
3696             # Not used internally
3697             sub Distance($$$) {
3698 14     14 1 3868 my ($self, $p0, $p1) = @_;
3699            
3700 14         42 $self->DistanceFrom($p0);
3701 14         36 return $self->DistanceTo($p1);
3702             }
3703              
3704              
3705              
3706              
3707             =head2 DistanceFrom( ... )
3708              
3709             =over
3710              
3711             C<$meters = $index-EDistanceFrom( \%point_1 );>
3712              
3713             C<$meters = $index-EDistanceFrom( \@point_1 );>
3714              
3715             Set an initial point to measure distances from
3716              
3717             Note that any call to C> and some calls to
3718             C> (those using the C or C
3719             options) will overwrite the initial point set with this method.
3720              
3721             B> or B>
3722              
3723             =over
3724              
3725             The point to measure distances from
3726              
3727             This can be either a hash containing at a minimum a C and a C value
3728             (both in degrees) or an array giving the point. See the B>
3729             section above for details.
3730              
3731              
3732             =back
3733              
3734             =back
3735              
3736             =cut
3737              
3738              
3739              
3740             # Not used internally
3741             sub DistanceFrom($$) {
3742 14     14 1 27 my ($self, $p0) = @_;
3743            
3744 14 100       39 if (ref $p0 eq 'ARRAY') {
3745             #. Got array; expand arguments into a full point
3746 5         15 $p0 = { 'lat'=>$$p0[0], 'lon'=>$$p0[1] };
3747             }
3748            
3749 14 100       47 $$p0{lat_rad} = Math::Trig::deg2rad($$p0{lat}) unless ($$p0{lat_rad});
3750 14 100       123 $$p0{lon_rad} = Math::Trig::deg2rad($$p0{lon}) unless ($$p0{lon_rad});
3751            
3752 14         121 SetUpDistance($self->{planetary_diameter}, $$p0{lat_rad}, $$p0{lon_rad});
3753             }
3754              
3755              
3756              
3757              
3758             =head2 DistanceTo( ... )
3759              
3760             =over
3761              
3762             C<$meters = $index-EDistanceTo( \%point_2 );>
3763              
3764             C<$meters = $index-EDistanceTo( \@point_2 );>
3765              
3766             Returns the distance in meters between the specified point and the one set
3767             earlier with C>.
3768              
3769             The haversine function is used to compute the distance. As this assumes a
3770             spherical body the distances returned may show errors. Using the default
3771             options, these errors are up to 0.056% (north - south) or 1.12% (east - west).
3772             Such errors typically start becoming significant at distances over S<20 km>.
3773              
3774             B> or B>
3775              
3776             =over
3777              
3778             The point to measure distances to
3779              
3780             This can be either a hash containing at a minimum a C and a C value
3781             (both in degrees) or an array giving the point. See the B>
3782             section above for details.
3783              
3784             =back
3785              
3786             =back
3787              
3788             =cut
3789              
3790              
3791              
3792             # Used by Distance(...)
3793             sub DistanceTo($$) {
3794 14     14 1 29 my ($self, $p1) = @_;
3795            
3796 14 100       37 if (ref $p1 eq 'ARRAY') {
3797             #. Got array; expand arguments into a full point
3798 5         13 $p1 = { 'lat'=>$$p1[0], 'lon'=>$$p1[1] };
3799             }
3800            
3801 14 100       44 $$p1{lat_rad} = Math::Trig::deg2rad($$p1{lat}) unless ($$p1{lat_rad});
3802 14 100       79 $$p1{lon_rad} = Math::Trig::deg2rad($$p1{lon}) unless ($$p1{lon_rad});
3803            
3804 14         121 return HaversineDistance($$p1{lat_rad}, $$p1{lon_rad});
3805             }
3806              
3807              
3808              
3809             #. Distance functions
3810             #.
3811             #. Geo::Index uses the haversine formula to compute great circle distances
3812             #. between points.
3813             #.
3814             #. Three versions are supported: a fallback version written in Perl (used if the
3815             #. C versions fail to compile) and two accelerated versions written in C, one
3816             #. using floats and the other using doubles. By default the C float version is
3817             #. used; if it fails to compile then the Perl version is used. Use of a specific
3818             #. version can also be requested with Geo::Index->SetDistanceFunctionType(...).
3819             #.
3820             #. The Perl version uses doubles. When using floats instead of doubles the loss
3821             #. of precision is typically under a meter (about 2 meters in the worst case).
3822             #. #. Compared to the errors inherent to the haversine function, this loss of
3823             #. precision is negligable.
3824              
3825             #. Here are the results of benchmarking the three versions on a fairly high-end
3826             #. workstation (higher numbers are better). The test dataset is 1 million random
3827             #. points and each search type was performed once for each point in random order.
3828             #. The same points were used for each test and they were in the same order. All
3829             #. searches returned results as lists except for the 'all points' search which
3830             #. returned a list reference. The default options (Earth, 20-level index) were
3831             #. used for each test. Each version's benchmark was run 32 times; some jitter
3832             #. was observed.
3833             #.
3834             #. Average number of operations per second using each version (rounded):
3835             #. Percentages (in parentheses) are relative to the pure-Perl version.
3836             #.
3837             #. Results for 32 iterations of each test:
3838             #.
3839             #. Operation Perl C double (%) C float (%)
3840             #. -------------------------------- ------ ------------ ------------
3841             #. Add point to index 35861 36067 (101) 36060 (101)
3842             #. Search: return all points 256127 256877 (101) 252116 (98)
3843             #. Search: sort, max 5 6733 8718 (129) 8860 (132)
3844             #. Search: sort, radius 1000, max 5 45364 49063 (108) 49831 (110)
3845             #. Search: sort, radius 1000 45902 49418 (108) 51673 (113)
3846             #. Search: max 5 198499 190404 (96) 204905 (103)
3847             #. Search: radius 1000, max 5 46942 51295 (109) 54554 (116)
3848             #. Search: radius 1000 47908 51941 (108) 55522 (116)
3849              
3850              
3851             #. These will be set to references to the Perl versions of the functions
3852             my $SetUpDistance_perl = undef;
3853             my $HaversineDistance_perl = undef;
3854             my $ComputeAreaExtrema_perl = undef;
3855             my $fast_log2_perl = undef;
3856              
3857             #. Choose whether to use Perl or C code for distance and log2 calculations
3858             #. Default is to use the Perl functions
3859              
3860             # Used by new(...)
3861             sub SetDistanceFunctionType($) {
3862 16     16 0 43 my ($type) = @_;
3863            
3864             #. Get function pointers for the Perl versions (if not already recorded)
3865 16 100       74 $SetUpDistance_perl = *SetUpDistance unless (defined $SetUpDistance_perl);
3866 16 100       56 $HaversineDistance_perl = *HaversineDistance unless (defined $HaversineDistance_perl);
3867 16 100       54 $ComputeAreaExtrema_perl = *ComputeAreaExtrema unless (defined $ComputeAreaExtrema_perl);
3868 16 100       48 $fast_log2_perl = *fast_log2 unless (defined $fast_log2_perl);
3869            
3870             #. Choose the type of functions to use:
3871            
3872 16 100 66     115 if ( $type eq 'perl' ) {
    100 33        
    50          
3873             #. Switch to using Perl code for distance and log2 calculations
3874            
3875 2         10 *Geo::Index::SetUpDistance = $SetUpDistance_perl;
3876 2         6 *Geo::Index::HaversineDistance = $HaversineDistance_perl;
3877 2         6 *Geo::Index::fast_log2 = $fast_log2_perl;
3878 2         6 *Geo::Index::ComputeAreaExtrema = $ComputeAreaExtrema_perl;
3879 2         5 $ACTIVE_CODE = 'perl';
3880            
3881 2         4 $C_CODE_ACTIVE = 0;
3882            
3883 2         4 return 1; # success
3884            
3885             } elsif ( $C_CODE_COMPILED && $type eq 'double' ) {
3886             #. Switch to using C double code for distance and log2 calculations
3887            
3888 12         42 *Geo::Index::SetUpDistance = *Geo::Index::SetUpDistance_double;
3889 12         30 *Geo::Index::HaversineDistance = *Geo::Index::HaversineDistance_double;
3890 12         25 *Geo::Index::fast_log2 = *Geo::Index::fast_log2_double;
3891 12         26 *Geo::Index::ComputeAreaExtrema = *Geo::Index::ComputeAreaExtrema_double;
3892 12         24 $ACTIVE_CODE = 'double';
3893            
3894 12         21 $C_CODE_ACTIVE = 1;
3895            
3896 12         32 return 1; # success
3897            
3898             } elsif ( $C_CODE_COMPILED && $type eq 'float' ) {
3899             #. Switch to using C float code for distance and log2 calculations
3900            
3901 2         11 *Geo::Index::SetUpDistance = *Geo::Index::SetUpDistance_float;
3902 2         7 *Geo::Index::HaversineDistance = *Geo::Index::HaversineDistance_float;
3903 2         6 *Geo::Index::ComputeAreaExtrema = *Geo::Index::ComputeAreaExtrema_float;
3904 2         6 *Geo::Index::fast_log2 = *Geo::Index::fast_log2_float;
3905            
3906 2         4 $ACTIVE_CODE = 'float';
3907            
3908 2         6 $C_CODE_ACTIVE = 1;
3909            
3910 2         5 return 1; # success
3911             }
3912            
3913 0         0 return undef; # Failed, no change
3914             }
3915              
3916              
3917             #. Returns the type of low-level functions that is active
3918             #. (one of 'perl', 'float', or 'double')
3919             # used by GetConfiguration, t/low-level.t
3920             sub GetLowLevelCodeType() {
3921 5     5 0 28 return $ACTIVE_CODE;
3922             }
3923              
3924             #. Returns reference to list of the supported low-level function types
3925             #. (list values as per GetLowLevelCodeType)
3926             # used by GetConfiguration, t/low-level.t
3927             sub GetSupportedLowLevelCodeTypes() {
3928 3     3 0 12 return [ @SUPPORTED_CODE ];
3929             }
3930              
3931             #. Perl version of the distance functions
3932              
3933             # Used by Search(...)
3934             #. For the values that the module is interested in the
3935             #. return value is the same as ceil(log2(n))
3936             sub fast_log2($) {
3937 0     0 0 0 my ($n) = @_;
3938 0         0 my $i = 0;
3939 0         0 my $c = 1;
3940 0         0 for ( $n = ceil( $n );
3941             $n > $c;
3942             $c<<=1, $i++ ) { }
3943 0         0 return $i;
3944             }
3945              
3946             #. Perl doesn't have a log2(n) function; if one wants
3947             #. to use it the following performs it:
3948              
3949             # Not used internally
3950             sub log2($) {
3951 0     0 0 0 my ($n) = @_;
3952 0         0 return log($n) / log(2);
3953             }
3954              
3955             #. Used internally by the Perl versions of SetUpDistance and HaversineDistance
3956             my ( $DistanceFrom_diameter, $DistanceFrom_lat_1, $DistanceFrom_lon_1 );
3957             my $DistanceFrom_cos_lat_1;
3958              
3959             #. Specify the point to get distances from
3960             #. Diameter is in meters, Lat and Lon are in radians
3961             #.
3962             #. If possible, this function will be replaced by an equivalent written in C.
3963             sub SetUpDistance($$$) {
3964 51     51 1 96 my ($new_diameter, $new_lat_1, $new_lon_1) = @_;
3965 51         75 $DistanceFrom_diameter = $new_diameter;
3966 51         62 $DistanceFrom_lat_1 = $new_lat_1;
3967 51         76 $DistanceFrom_lon_1 = $new_lon_1;
3968 51         109 $DistanceFrom_cos_lat_1 = cos( $DistanceFrom_lat_1 );
3969             }
3970              
3971             #. Returns the approximate distance from previously-set point to specified point
3972             #. Lat and Lon are in radians, return value is in meters
3973             #.
3974             #. If possible, this function will be replaced by an equivalent written in C.
3975             sub HaversineDistance($$) {
3976 161     161 0 237 my ($lat_0, $lon_0)= @_;
3977            
3978 161         271 my $sin_lat_diff_over_2 = sin( ( $lat_0 - $DistanceFrom_lat_1 ) / 2.0 );
3979 161         249 my $sin_lon_diff_over_2 = sin( ( $lon_0 - $DistanceFrom_lon_1 ) / 2.0 );
3980            
3981 161         347 my $n = ( $sin_lat_diff_over_2 * $sin_lat_diff_over_2 )
3982             + (
3983             ( $sin_lon_diff_over_2 * $sin_lon_diff_over_2 )
3984             * $DistanceFrom_cos_lat_1
3985             * cos( $lat_0 )
3986             );
3987            
3988             #. The haversine formula may get messy around antipodal points so clip to the largest sane value.
3989 161 50       295 if ( $n < 0.0 ) { $n = 0.0; }
  0         0  
3990            
3991 161         371 return $DistanceFrom_diameter * asin( sqrt($n) );
3992             }
3993              
3994              
3995              
3996              
3997             =head2 GetConfiguration( )
3998              
3999             =over
4000              
4001             C<%configuration = $index-EGetConfiguration( );>
4002              
4003             Returns the running configuration of the Geo::Index object.
4004              
4005             See also C> and C
4006              
4007             The return value is a hash with the following entries:
4008              
4009             =over
4010              
4011             B> - The key type in use:
4012              
4013             =over
4014              
4015             'C' for text keys (e.g. 'C<12:345,6789>')
4016              
4017             'C' for 64-bit numeric keys
4018              
4019             'C' for 64-bit numeric keys packed into an 8-byte string
4020              
4021             =back
4022              
4023             B> - The types of keys that can be used
4024              
4025             =over
4026              
4027             Value is a reference to a list of supported key types (as given above).
4028              
4029             =back
4030              
4031             B> - The type of low-level code in use:
4032              
4033             =over
4034              
4035             'C' for Perl functions
4036              
4037             'C' for C functions mostly using C values.
4038              
4039             'C' for C functions mostly using C values.
4040              
4041             =back
4042              
4043             B> - The types of low-level code that can be used
4044              
4045             =over
4046              
4047             Value is a reference to a list of supported code types (as given above).
4048              
4049             =back
4050              
4051             B> - Number of levels in index (excluding the global level)
4052            
4053             B> - Average planetary radius (in meters)
4054              
4055             B> - Polar circumference (in meters)
4056              
4057             B> - Equatorial circumference (in meters)
4058            
4059             B> - Number of points currently indexed
4060              
4061             B> - Size in meters (at the equator) of each tile the at
4062             most-detailed level of index
4063              
4064             =back
4065              
4066             =back
4067              
4068             =cut
4069              
4070              
4071             #. Returns the index's current configuration
4072             # not used internally
4073             sub GetConfiguration($) {
4074 2     2 1 899 my ($self) = @_;
4075 2         4 my %config = ();
4076            
4077             #. Low-level configuration
4078 2         6 $config{key_type} = ( USE_NUMERIC_KEYS ) ? ( USE_PACKED_KEYS ) ? 'packed' : 'numeric' : 'text';
4079 2         5 $config{supported_key_types} = [ 'text', 'numeric', 'packed' ];
4080 2         6 $config{code_type} = $self->GetLowLevelCodeType();
4081 2         7 $config{supported_code_types} = $self->GetSupportedLowLevelCodeTypes();
4082            
4083 2         6 $config{module_version} = "$VERSION";
4084 2         11 $config{module_version} =~ s/^v//;
4085            
4086 2 50       8 if ($C_CODE_COMPILED == 1) {
4087             #. C low-level function library is loaded
4088            
4089 2         10 my $c_code_version = Geo::Index::GetCCodeVersion();
4090 2         4 my $mask = ( 1 << 10 ) - 1;
4091 2         5 my $major_version = ( $c_code_version >> 20 ) & $mask;
4092 2         4 my $minor_version = ( $c_code_version >> 10 ) & $mask;
4093 2         11 my $sub_version = $c_code_version & $mask;
4094 2         8 $config{c_code_version} = "${major_version}.${minor_version}.${sub_version}";
4095            
4096             } else {
4097             #. No C low-level function library
4098 0         0 $config{c_code_version} = undef;
4099             }
4100            
4101             #. Index depth
4102 2         6 $config{levels} = $self->{levels};
4103            
4104             #. Planery size
4105 2         6 $config{planetary_radius} = $self->{planetary_radius};
4106 2         3 $config{polar_circumference} = $self->{polar_circumference};
4107 2         5 $config{equatorial_circumference} = $self->{equatorial_circumference};
4108            
4109             #. Number of points in index
4110 2         3 $config{size} = scalar keys %{$self->{indices}};
  2         7  
4111            
4112             #. Width in meters of each tile at most-detailed level of index
4113 2         6 $config{tile_meters} = $config{equatorial_circumference} / ( 2**$config{levels} );
4114            
4115 2         68 return %config;
4116             }
4117              
4118              
4119              
4120              
4121             =head2 GetStatistics( )
4122              
4123             =over
4124              
4125             C<@stats = $index-EGetStatistics( );>
4126              
4127             Returns statistics regarding the Geo::Index object.
4128              
4129             See also C> and C
4130              
4131             The return value is a list with one entry per level. Each list entry is a hash
4132             reference giving statistics for a single level of the index and contains the
4133             following entries:
4134              
4135             =over
4136              
4137             B> - The level number the statistics are for
4138              
4139             B> - Total number of points indexed in this level
4140              
4141             B> - Number of tiles containing at least one point
4142              
4143             B> - Minimum number of points in a non-empty tile
4144              
4145             B> - Maximum number of points in a non-empty tile
4146              
4147             B> - Average number of points in a non-empty tile
4148              
4149             =back
4150              
4151             =back
4152              
4153             =cut
4154              
4155              
4156              
4157             #. Returns statistics regarding the current index
4158             # not used internally
4159             sub GetStatistics($) {
4160 0     0 1 0 my ($self) = @_;
4161            
4162 0         0 my $_index = $self->{index};
4163 0         0 my $levels = $self->{levels};
4164            
4165 0         0 my @stats = ();
4166            
4167 0         0 foreach my $key (keys %$_index) {
4168 0         0 my ($level, $lat, $lon);
4169            
4170 0         0 if ( USE_NUMERIC_KEYS ) {
4171 0         0 if ( USE_PACKED_KEYS ) {
4172             #. packed numeric key
4173             $key = unpack("Q", $key);
4174             }
4175             #. numeric key
4176            
4177 0         0 $level = $key >> 58;
4178 0         0 $lat = ($key >> 29 ) & MASK_LATLON;
4179 0         0 $lon = $key & MASK_LATLON;
4180            
4181 0 0       0 $lat = 'ALL' if ($lat == MASK_LATLON);
4182 0 0       0 $lon = 'ALL' if ($lon == MASK_LATLON);
4183            
4184 0 0       0 next if ($level > $levels);
4185             } else {
4186             #. text key
4187             ($level,$lat,$lon) = split /[:,]/, $key;
4188            
4189             next if ($level eq 'ALL');
4190             }
4191            
4192 0         0 $stats[$level]->{level} = $level;
4193            
4194 0         0 my $count = scalar @{$$_index{$key}};
  0         0  
4195            
4196 0         0 $stats[$level]->{points} += $count;
4197 0         0 $stats[$level]->{tiles}++;
4198 0         0 $stats[$level]->{avg_tile_points} += $count;
4199            
4200 0 0 0     0 if (
4201             ( ! defined $stats[$level]->{min_tile_points} ) ||
4202             ( $count < $stats[$level]->{min_tile_points} )
4203             ) {
4204 0         0 $stats[$level]->{min_tile_points} = $count;
4205             }
4206 0 0 0     0 if (
4207             ( ! defined $stats[$level]->{max_tile_points} ) ||
4208             ( $count > $stats[$level]->{max_tile_points} )
4209             ) {
4210 0         0 $stats[$level]->{max_tile_points} = $count;
4211             }
4212            
4213             }
4214            
4215 0         0 for (my $level=0; $level<$levels; $level++) {
4216 0         0 my $tiles = $stats[$level]->{tiles};
4217 0 0       0 if ( $tiles ) {
4218 0         0 $stats[$level]->{avg_tile_points} /= $tiles;
4219             } else {
4220 0         0 $stats[$level]->{avg_tile_points} = undef;
4221             }
4222             }
4223            
4224 0         0 return @stats;
4225             }
4226              
4227              
4228              
4229              
4230             =head2 Sweep( ... )
4231              
4232             =over
4233              
4234             C<$index-ESweep( );>
4235              
4236             C<$index-ESweep( \%point );>
4237              
4238             C<$index-ESweep( \@points );>
4239              
4240             C<$index-ESweep( undef, \@extra_keys );>
4241              
4242             C<$index-ESweep( \%point, \@extra_keys );>
4243              
4244             C<$index-ESweep( \@points, \@extra_keys );>
4245              
4246             Remove data generated by searches from some or all points
4247              
4248             The fields that will be removed are C and C.
4249              
4250             Called on its own (with no point or points specified) this method will remove
4251             data generated by searches from all points.
4252              
4253             A list of additional keys to remove can optionally be supplied. To request
4254             vacuuming of all points with additional keys specified, use C instead
4255             of C<\%point> or C<\@points>.
4256              
4257             See also C>.
4258              
4259             B> or B>
4260              
4261             =over
4262              
4263             The point or a list of points to remove metadata from.
4264              
4265             =back
4266              
4267             B>
4268              
4269             =over
4270              
4271             List of additional keys to remove
4272              
4273             =back
4274              
4275             =back
4276              
4277             =cut
4278              
4279              
4280             #. Remove data generated search methods from some or all points
4281             # not used internally
4282             sub Sweep($;$$) {
4283 0     0 1 0 my ($self, $_points, $_extra_keys) = @_;
4284            
4285 0 0       0 if ( ! defined $_points ) {
    0          
4286             # Use all points in index if none were specified
4287 0         0 my $key;
4288 0         0 if (USE_NUMERIC_KEYS) {
4289 0         0 if (USE_PACKED_KEYS) {
4290             $key = pack("Q", ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON );
4291             } else {
4292 0         0 $key = ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON;
4293             }
4294             } else {
4295             $key = [ ALL, ALL, ALL ];
4296             }
4297 0         0 $_points = $self->GetValue($key);
4298            
4299             } elsif ( ref $_points eq 'HASH' ) {
4300             # Build list if passed a single point
4301 0         0 $_points = [ $_points ];
4302             }
4303            
4304 0 0       0 if ( ref $_extra_keys eq 'ARRAY' ) {
4305 0         0 foreach my $_point ( @$_points ) {
4306 0         0 delete $$_point{search_result_distance};
4307 0         0 delete $$_point{antipode_distance};
4308 0         0 foreach my $key (@$_extra_keys) {
4309 0         0 delete $$_point{$key};
4310             }
4311             }
4312             } else {
4313 0         0 foreach my $_point ( @$_points ) {
4314 0         0 delete $$_point{search_result_distance};
4315 0         0 delete $$_point{antipode_distance};
4316             }
4317             }
4318             }
4319              
4320              
4321              
4322              
4323             =head2 Vacuum( ... )
4324              
4325             =over
4326              
4327             C<$index-EVacuum( );>
4328              
4329             C<$index-EVacuum( \%point );>
4330              
4331             C<$index-EVacuum( \@points );>
4332              
4333             C<$index-EVacuum( undef, \@extra_keys );>
4334              
4335             C<$index-EVacuum( \%point, \@extra_keys );>
4336              
4337             C<$index-EVacuum( \@points, \@extra_keys );>
4338              
4339             Remove all data generated by Geo::Index from some or all points
4340              
4341             The fields that will be removed are: C, C, C,
4342             C, C.
4343              
4344             Called on its own (with no point or points specified) this method will remove
4345             all generated data from all points.
4346              
4347             A list of additional keys to remove can optionally be supplied. To request
4348             vacuuming of all points with additional keys specified, use C instead
4349             of C<\%point> or C<\@points>.
4350              
4351             See also C>.
4352              
4353             B> or B>
4354              
4355             =over
4356              
4357             The point or a list of points to remove metadata from.
4358              
4359             =back
4360              
4361             B>
4362              
4363             =over
4364              
4365             List of additional keys to remove
4366              
4367             =back
4368              
4369             =back
4370              
4371             =cut
4372              
4373              
4374             #. Remove all data generated by Geo::Index from some or all points
4375             # not used internally
4376             sub Vacuum($;$$) {
4377 0     0 1 0 my ($self, $_points, $_extra_keys) = @_;
4378            
4379 0 0       0 if ( ! defined $_points ) {
    0          
4380             # Use all points in index if none were specified
4381 0         0 my $key;
4382 0         0 if (USE_NUMERIC_KEYS) {
4383 0         0 if (USE_PACKED_KEYS) {
4384             $key = pack("Q", ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON );
4385             } else {
4386 0         0 $key = ( MASK_LEVEL << 58 ) | ( MASK_LATLON << 29 ) | MASK_LATLON;
4387             }
4388             } else {
4389             $key = [ ALL, ALL, ALL ];
4390             }
4391 0         0 $_points = $self->GetValue($key);
4392            
4393             } elsif ( ref $_points eq 'HASH' ) {
4394             # Build list if passed a single point
4395 0         0 $_points = [ $_points ];
4396             }
4397            
4398 0 0       0 if ( ref $_extra_keys eq 'ARRAY' ) {
4399 0         0 foreach my $_point ( @$_points ) {
4400 0         0 delete $$_point{lat_rad};
4401 0         0 delete $$_point{lon_rad};
4402 0         0 delete $$_point{circumference};
4403 0         0 delete $$_point{search_result_distance};
4404 0         0 delete $$_point{antipode_distance};
4405 0         0 foreach my $key (@$_extra_keys) {
4406 0         0 delete $$_point{$key};
4407             }
4408             }
4409             } else {
4410 0         0 foreach my $_point ( @$_points ) {
4411 0         0 delete $$_point{lat_rad};
4412 0         0 delete $$_point{lon_rad};
4413 0         0 delete $$_point{circumference};
4414 0         0 delete $$_point{search_result_distance};
4415 0         0 delete $$_point{antipode_distance};
4416             }
4417             }
4418             }
4419              
4420              
4421              
4422              
4423             =head2 PointCount( )
4424              
4425             =over
4426              
4427             C<$count = $index-EPointCount( );>
4428              
4429             Returns the number of points currently in index
4430              
4431             =back
4432              
4433             =cut
4434              
4435              
4436             #. Returns the number of points currently in index
4437             # Used by t/trampoline.t
4438             sub PointCount($) {
4439 7     7 1 33 my ($self) = @_;
4440 7         16 return scalar keys %{$self->{indices}}
  7         43  
4441             }
4442              
4443              
4444              
4445              
4446             =head2 AllPoints( )
4447              
4448             =over
4449              
4450             C<@all_points = $index-EAllPoints( );>
4451              
4452             C<$all_points = $index-EAllPoints( );>
4453              
4454             Returns all points currently in index
4455              
4456             B
4457              
4458             =over
4459              
4460             In list context the return value is a list of references to all points in the
4461             index or an empty array if there are no points in the index.
4462              
4463             In scalar context the return value is a reference to the aforementioned list or
4464             a reference to an empty list if there are no points in the index.
4465              
4466             =back
4467              
4468             =back
4469              
4470             =cut
4471              
4472              
4473             #. Returns a list of all points currently in index
4474             sub AllPoints($) {
4475 1     1 1 865 my ($self) = @_;
4476              
4477 1         9 my $_results = $self->Search( [0,0], { radius => ALL, quick_results => 1 } );
4478            
4479 1 50       5 if (wantarray) {
4480 0         0 return @{$$_results[0]};
  0         0  
4481             } else {
4482 1         4 return $$_results[0];
4483             }
4484            
4485             # Slower method:
4486            
4487             # my @return = ();
4488             # foreach my $_set ( @$_results ) {
4489             # push @return, @$_set if ( defined $_set );
4490             # }
4491             #
4492             # if (wantarray) {
4493             # return @return;
4494             # } else {
4495             # return \@return;
4496             # }
4497             }
4498              
4499              
4500              
4501              
4502             # ==============================================================================
4503              
4504             #. These methods are experimental and may be removed in future
4505              
4506              
4507             #. END Experimental methods
4508              
4509             # ==============================================================================
4510              
4511             #. The code below is only used internally by Geo::Index
4512              
4513              
4514             #: Proximity and indexing helper functions
4515              
4516             #. Note that for performance reasons the code from these functions may have been
4517             #. inlined into the various methods and C functions.
4518              
4519              
4520             #: Convert latitude/longitude pair from degrees to integers suitable for indexing
4521             #.
4522             #> IN: lat, lon -> The point's latitude and longitude in degrees
4523             #>
4524             #> OUT: lat_int, lon_int -> The point's latitude and longitude as scaled integers
4525             #.
4526             #. Input ranges are:
4527             #. latitude: [ -90 degrees .. 90 degrees ]
4528             #. longitude: [ -180 degrees .. 180 degrees )
4529             #.
4530             #. Output ranges are:
4531             #. latitude: [ 0 .. max_size ] = [ 0 .. 2**levels ]
4532             #. longitude: [ 0 .. max_size - 1 ] = [ 0 .. 2**levels - 1 ]
4533             #.
4534             #. Note that the latitude range includes both poles. The longitude range only
4535             #. includes the antimeridian once and thus its output integer range is one
4536             #. smaller. Code that converts the returned integers to index keys must take
4537             #. this into account.
4538              
4539             # Used by Index and Search
4540             sub GetIntLatLon($$) {
4541 754     754 1 1274 my ($self, $lat, $lon) = @_;
4542            
4543 754         1421 my $lat_int = int( ( $lat + 90.0 ) * $self->{max_size} / 180.0 );
4544 754 50       1528 $lat_int = $self->{max_size} if ($lat_int > $self->{max_size});
4545            
4546 754         1356 my $lon_int = int( ( $lon + 180.0 ) * $self->{max_size} / 360.0 ) % $self->{max_size};
4547            
4548 754         1405 return ($lat_int, $lon_int);
4549             }
4550              
4551              
4552             #: Convert latitude from degrees to integers suitable for indexing
4553             #.
4554             #> IN: lat -> The point's latitude in degrees
4555             #>
4556             #> OUT: lat_int -> The point's latitude as scaled integers
4557             #.
4558             #. Input ranges are:
4559             #. latitude: [ -90 degrees .. 90 degrees ]
4560             #.
4561             #. Output ranges are:
4562             #. latitude: [ 0 .. max_size ] = [ 0 .. 2**levels ]
4563             #.
4564             #. Note that the latitude range includes both poles. Code that converts the
4565             #. returned integer to part of an index keys must take this into account.
4566              
4567             # Not used
4568             sub GetIntLat($$) {
4569 0     0 1 0 my ($self, $lat) = @_;
4570            
4571 0         0 my $lat_int = int( ( $lat + 90.0 ) * $self->{max_size} / 180.0 );
4572 0 0       0 $lat_int = $self->{max_size} if ($lat_int > $self->{max_size});
4573            
4574 0         0 return ($lat_int);
4575             }
4576              
4577              
4578             #: Convert longitude from degrees to integers suitable for indexing
4579             #.
4580             #> IN: lon -> The point's longitude in degrees
4581             #>
4582             #> OUT: lon_int -> The point's longitude as scaled integers
4583             #.
4584             #. Input ranges are:
4585             #. longitude: [ -180 degrees .. 180 degrees )
4586             #.
4587             #. Output ranges are:
4588             #. longitude: [ 0 .. max_size - 1 ] = [ 0 .. 2**levels - 1 ]
4589              
4590             # Not used
4591             sub GetIntLon($$) {
4592 0     0 0 0 my ($self, $lon) = @_;
4593 0         0 return int( ( $lon + 180.0 ) * $self->{max_size} / 360.0 ) % $self->{max_size};
4594             }
4595              
4596              
4597             #: Return length of circle of latitude (in meters) at given latitude (in radians)
4598             #.
4599             #. Circles of latitude run east-west (e.g. the equator is a circle of latitude).
4600             #. Values are approximate. The diameters are those for an oblate spheroid but
4601             #. the math assumes a sphere.
4602              
4603             # Use by Closest, OneMeterInDegrees
4604             sub LongitudeCircumference($$) {
4605 16     16 0 33 my ($self, $radians) = @_;
4606 16         57 return abs(cos( $radians )) * $self->{equatorial_circumference};
4607             }
4608              
4609              
4610             =head2 OneMeterInDegrees( ... )
4611              
4612             =over
4613              
4614             C<$index-EOneMeterInDegrees( $latitude );>
4615              
4616             Returns length in degrees of one meter (N/S and E/W) at given latitude
4617              
4618             Values are approximate. The diameters are those for an oblate spheroid but
4619             the math assumes a sphere. As one approaches the poles these values get
4620             heavily distorted; code that uses them needs to take this into account.
4621              
4622             See also C>.
4623              
4624             B>
4625              
4626             =over
4627              
4628             The latitude in radians
4629              
4630             =back
4631              
4632             B
4633              
4634             =over
4635              
4636             An two-element list containing the width and height in meters of one degree
4637             at the given latitude:
4638              
4639             C<( I, I )>
4640              
4641             =back
4642              
4643             =back
4644              
4645             =cut
4646              
4647              
4648             #: Returns length in degrees of one meter (N/S and E/W) at given latitude
4649             #.
4650             #> IN: latitude -> The latitude in radians
4651             #>
4652             #> OUT: An array containing the width and height in meters of one degree at the
4653             #> given latitude: ( DEGREES_OF_LATITUDE, DEGREES_OF_LONGITUDE )
4654             #.
4655             #. Values are approximate. The diameters are those for an oblate spheroid but
4656             #. the math assumes a sphere. As one approaches the poles these values get
4657             #. heavily distorted; code that uses them needs to take this into account.
4658              
4659             # Not used
4660             sub OneMeterInDegrees($$) {
4661 0     0 1 0 my ($self, $latitude) = @_;
4662            
4663 0         0 my $NS_circumference_in_meters = $self->{polar_circumference};
4664 0         0 my $EW_circumference_in_meters = LongitudeCircumference( $self, $latitude );
4665            
4666 0 0       0 if ( $EW_circumference_in_meters ) {
4667 0         0 return ( ( 360.0 / $NS_circumference_in_meters ), ( 360.0 / $EW_circumference_in_meters ) )
4668             } else {
4669 0         0 return ( ( 360.0 / $NS_circumference_in_meters ), undef )
4670             }
4671             }
4672              
4673              
4674             =head2 OneDegreeInMeters( ... )
4675              
4676             =over
4677              
4678             C<$index-EOneMeterInDegrees( $latitude );>
4679              
4680             Returns length in meters of one degree (N/S and E/W) at given latitude
4681              
4682             Values are approximate. The diameters are those for an oblate spheroid but
4683             the math assumes a sphere. As one approaches the poles these values get
4684             heavily distorted; code that uses them needs to take this into account.
4685              
4686             See also C>.
4687              
4688             B>
4689              
4690             =over
4691              
4692             The latitude in radians
4693              
4694             =back
4695              
4696             B
4697              
4698             =over
4699              
4700             An two-element list containing the width and height in degrees of one meter
4701             at the given latitude:
4702              
4703             C<( I, I )>
4704              
4705             =back
4706              
4707             =back
4708              
4709             =cut
4710              
4711              
4712             #: Returns length in meters of one degree (N/S and E/W) at given latitude
4713             #.
4714             #> IN: latitude -> The latitude in radians
4715             #>
4716             #> OUT: An array containing the width and height in degrees of one meter at the
4717             #> given latitude: ( NORTH_SOUTH_METERS, EAST_WEST_METERS )
4718             #.
4719             #. Values are approximate. The diameters are those for an oblate spheroid but
4720             #. the math assumes a sphere. As one approaches the poles these values get
4721             #. heavily distorted; code that uses them needs to take this into account.
4722              
4723             # Not used
4724             sub OneDegreeInMeters($$) {
4725 0     0 1 0 my ($self, $latitude) = @_;
4726            
4727 0         0 my $NS_circumference_in_meters = $self->{polar_circumference};
4728 0         0 my $EW_circumference_in_meters = LongitudeCircumference( $self, $latitude );
4729            
4730 0         0 return ( ( $NS_circumference_in_meters / 360.0 ), ( $EW_circumference_in_meters / 360.0 ) )
4731             }
4732              
4733              
4734             #: Returns the latitude and longitude indices for a point at a given level
4735             #. This method is used by Index(...)
4736             #.
4737             #> IN: level -> The level number to get indices for
4738             #> lat_int -> Point's integer latitude (as returned by GetIntLatLon(...)
4739             #> lon_int -> Point's integer longitude (as returned by GetIntLatLon(...)
4740             #>
4741             #> OUT: lat_idx -> Point's integer latitude within the requested level
4742             #> lon_idx -> Point's integer longitude within the requested level
4743              
4744             # Used by Index
4745             sub GetIndices($$$) {
4746 687     687 0 1189 my ($self, $level, $lat_int, $lon_int) = @_;
4747            
4748 687         970 my $shift = $self->{max_level} - $level;
4749            
4750 687         954 my $lat_idx = $lat_int >> $shift;
4751 687         863 my $lon_idx = $lon_int >> $shift;
4752            
4753 687         1179 return ($lat_idx, $lon_idx);
4754             }
4755              
4756              
4757              
4758              
4759             # *** Method aliases ***
4760              
4761             sub index; *index = *Index;
4762             sub index_points; *index_points = *IndexPoints;
4763             sub unindex; *unindex = *Unindex;
4764              
4765             sub build_points; *build_points = *BuildPoints;
4766             sub add_value; *add_value = *AddValue;
4767             sub get_value; *get_value = *GetValue;
4768              
4769             sub search; *search = *Search;
4770             sub search_by_bounds; *search_by_bounds = *SearchByBounds;
4771             sub closest; *closest = *Closest;
4772             sub farthest; *farthest = *Farthest;
4773             sub all_points; *all_points = *AllPoints;
4774              
4775             sub distance; *distance = *Distance;
4776             sub distance_from; *distance_from = *DistanceFrom;
4777             sub distance_to; *distance_to = *DistanceTo;
4778              
4779             sub one_meter_in_degrees; *one_meter_in_degrees = *OneMeterInDegrees;
4780             sub one_degree_in_meters; *one_degree_in_meters = *OneDegreeInMeters;
4781              
4782             sub get_configuration; *get_configuration = *GetConfiguration;
4783             sub get_statistics; *get_statistics = *GetStatistics;
4784              
4785             sub sweep; *sweep = *Sweep;
4786             sub vacuum; *vacuum = *Vacuum;
4787             sub point_count; *point_count = *PointCount;
4788              
4789             # Internal methods:
4790              
4791             sub set_distance_function_type; *set_distance_function_type = *SetDistanceFunctionType;
4792             sub get_low_level_code_type; *get_low_level_code_type = *GetLowLevelCodeType;
4793             sub get_supported_low_level_code_types; *get_supported_low_level_code_types = *GetSupportedLowLevelCodeTypes;
4794             sub set_up_distance; *set_up_distance = *SetUpDistance;
4795             sub haversine_distance; *haversine_distance = *HaversineDistance;
4796             sub get_int_lat_lon; *get_int_lat_lon = *GetIntLatLon;
4797             sub get_int_lat; *get_int_lat = *GetIntLat;
4798             sub get_int_lon; *get_int_lon = *GetIntLon;
4799             sub longitude_circumference; *longitude_circumference = *LongitudeCircumference;
4800             sub get_indices; *get_indices = *GetIndices;
4801              
4802             =head2 Alternate method names
4803              
4804             Geo::Index uses CamelCase for its method names.
4805              
4806             For those who prefer using snake case, alternate method names are provided:
4807              
4808             Method Alternate name
4809             ----------------- --------------------
4810             Index index
4811             IndexPoints index_points
4812             Unindex unindex
4813            
4814             BuildPoints build_points
4815             AddValue add_value
4816             GetValue get_value
4817            
4818             Search search
4819             SearchByBounds search_by_bounds
4820             Closest closest
4821             Farthest farthest
4822             AllPoints all_points
4823            
4824             Distance distance
4825             DistanceFrom distance_from
4826             DistanceTo distance_to
4827            
4828             OneMeterInDegrees one_meter_in_degrees
4829             OneDegreeInMeters one_degree_in_meters
4830            
4831             GetConfiguration get_configuration
4832             GetStatistics get_statistics
4833            
4834             Sweep sweep
4835             Vacuum vacuum
4836             PointCount point_count
4837              
4838              
4839             =head1 CONDITION FUNCTIONS
4840              
4841             The C>, C>,
4842             C>, and C>
4843             methods allow a user-supplied condition function to filter potential results.
4844              
4845             If present, these condition functions are called for each potential search
4846             result. They should be idempotent* and could potentially be called multiple
4847             times for a given point. The code should return B (e.g. C<1>) if a potential
4848             point should be included in the results or B (e.g. C<0> or C) if the
4849             point should be excluded.
4850              
4851             For C>, C>, and
4852             C>, the C function runs before
4853             the distance to the result point has been calculated and the C
4854             function runs after it has been calculated. For C>
4855             no distances are calculated and the function is simply called once per point.
4856              
4857             * Functions can set outside values provided they do not affect any values
4858             used internally by C> and so long as those
4859             outside values have no effect on the condition's outcome. Such behavior is,
4860             of course, frowned upon.
4861              
4862             The parameters to the condition function are, in order:
4863              
4864             =over
4865              
4866             B>
4867              
4868             =over
4869              
4870             Reference to the potential search result being checked
4871              
4872             =back
4873              
4874             B>
4875              
4876             =over
4877              
4878             Reference to the point at the center of the search
4879              
4880             For C this is instead the bounding box:
4881             S, I, I, I ]>>
4882              
4883             =back
4884              
4885             B>
4886              
4887             =over
4888              
4889             Arbitrary user-supplied data
4890              
4891             =back
4892              
4893             =back
4894              
4895             For example, the options set in the following code allows all points in the
4896             results except for the one named 'S':
4897              
4898             $options{pre_condition} =
4899             sub {
4900             my ( $_result_point, $_search_point, $user_data ) = @_;
4901             if ( $$_result_point{name} eq $user_data ) {
4902             return 0; # Exclude result
4903             }
4904             return 1; # Point is a valid search result
4905             };
4906             $options{user_data} = "Point Nada";
4907              
4908             To exclude the search point from the search results use:
4909              
4910             $options{post_condition} =
4911             sub {
4912             my ( $_result_point, $_search_point, $user_data ) = @_;
4913             return ( $_result_point != $_search_point );
4914             };
4915              
4916             or more concisely
4917              
4918             $options{post_condition} = sub { return $_[0] != $_[1]; };
4919              
4920             In general, C functions should be preferred since the overhead
4921             of the Perl function call is typically larger than that of the distance
4922             calculation. By checking the distance first, running the C
4923             function might not be necessary.
4924              
4925              
4926              
4927              
4928             =head1 PERFORMANCE
4929              
4930             =head2 Overview
4931              
4932             Geo::Index is intended for stand-alone applications that need a way to quickly
4933             perform proximity searches on relatively small datasets (at most a few million
4934             points). Typical search speeds are three to five orders of magnitude faster
4935             than a linear search. For larger datasets and for applications running in a
4936             server environment using something like PostGIS is more appropriate.
4937              
4938             Indexing speed is about 50,000 points per second when C is 20. Search
4939             speeds are highly dependent on the data indexed and on search parameters but are
4940             typically in the neighborhood of a few thousand searches per second.
4941              
4942             Memory usage tends to be rather high; for 1,000,000 points the index is S<~3.2 GB>
4943             for tightly clustered points or S<~4.6 GB> when spread evenly world-wide. The
4944             size of an index grows linearly with each added point at a rate of about S<4 kB>
4945             per point. When a point is first encountered whilst searching its size will
4946             increase by about 100 bytes (this only happens once per point).
4947              
4948             Since performance is so dependant on data and usage, the the user is encouraged
4949             to test all available options while developing their application before choosing
4950             the one that works fastest. The C script included with
4951             this module may be helpful for measuring this module's performance.
4952              
4953             =head2 General tips
4954              
4955             Here are some guidelines for best results:
4956              
4957             =over
4958              
4959             =item * B
4960              
4961             That is, e.g., S> is faster than C<@results = Search(...);>
4962              
4963             =item * B
4964              
4965             Benchmarking has shown that the cost of the Perl function call is higher than
4966             that of the distance-related code. Thus there is probably no reason to use
4967             pre-conditions. Put concisely,
4968              
4969             =over
4970              
4971             $results = $index->Search( $point, { post_condition => $code_ref } );
4972              
4973             =back
4974              
4975             is faster than
4976              
4977             =over
4978              
4979             $results = $index->Search( $point, { pre_condition => $code_ref } );
4980              
4981             =back
4982              
4983             =item * B when creating the index>
4984              
4985             The C> method has best performance when the size of the most
4986             detailed level of the index has a smaller physical size than the radius of a
4987             typical search. For example, if your searches are typically for points within
4988             100 meters then an index with C should be set to at least 18 (~75 meters
4989             at the equator) to yield best results; if typical searches have 10 meter radius
4990             then C should be 22.
4991              
4992             The C> method works best when the most detailed level of the
4993             index contains a single point per tile and search points lie close to potential
4994             result points.
4995              
4996             To help tune the C value, the C> method can
4997             be used to find out the physical size of the most detailed level along with
4998             statistics on the number of points per index tile.
4999              
5000             =item * B option when possible.>
5001              
5002             Filtering points and combining them into a single, flat list can be very
5003             expensive. Many applications can tolerate getting additional points beyond
5004             those matching the search criteria. An example of this is drawing points on
5005             a map; if points are clipped to the visible area when they are
5006             drawn it may not matter if some of them lie outside of it.
5007              
5008             =item * B> instead of C> when you have a search radius.>
5009              
5010             The C> function is most efficient when no search radius is
5011             specified or when result points lie very close to the search point. Closeness
5012             is relative to the tile size of the most detailed index level; for the default
5013             index depth (C<20>), "very close" is roughly within about 100 meters.
5014              
5015             When clipping results to a maximal radius it is typically much faster to use
5016             C> with the C and C options*.
5017              
5018             For example, to find the closest C<$n> points within distance C<$d> of a point
5019             C<$p> it is usually much faster to use
5020              
5021             =over
5022              
5023             %options = (
5024             max_results => $n,
5025             radius => $d,
5026             sort_results => 1,
5027             post_condition => sub { return $_[0] != $_[1]; }
5028             );
5029             $results = $index->Search( $p, \%options );
5030              
5031             =back
5032              
5033             instead of
5034              
5035             =over
5036              
5037             $results = $index->Closest( $p, $n { radius => $d } );
5038              
5039             =back
5040              
5041             * The C shown in the example omits the search point from
5042             the results and is needed to fully emulate the behavior of C.
5043              
5044             =back
5045              
5046             =head2 Technical discussion
5047              
5048             Both C> and C> are very fast since
5049             they can find the relevant index tiles in linear time. Since the time needed to
5050             filter the results is directly proportional to the number of points retrieved
5051             from the index, best performance occurs when the size of the most detailed tiles
5052             is smaller than that of the typical search radius or search bounds.
5053              
5054             Searches run using C> are done starting from the most
5055             detailed level and work upwards. Best performance occurs when a result is found
5056             in the first few iterations. If the first iteration that finds points yields a
5057             large number of points then performance will suffer since the distance to each
5058             of these points will need to be measured to find the closest. For similar
5059             reasons, requesting a large number of closest points in a single call will also
5060             impact performance. The C> method is largely a wrapper for
5061             C> and thus exhibits similar behavior.
5062              
5063             Some functions within Geo::Index have optional implementations written in C. If
5064             these are active (by default they are whenever possible) searches typically run
5065             25% to 50% faster.
5066              
5067             Whenever possible Geo::Index uses numeric index keys. Compared to text index
5068             keys, numeric keys improve performance with about 30% faster speed and about 50%
5069             smaller index memory usage. The downside to numeric keys is that they are less
5070             legible to humans while debugging. (Whether numeric or text keys are used can
5071             be changed by setting the appropriate value at the top of C)
5072              
5073             =head2 Benchmark results
5074              
5075             Typical benchmark results run on a modern workstation using numeric keys and
5076             double-precision C low-level code with the index containing 1,000,000 points
5077             are as follows:
5078              
5079             =over
5080              
5081             =item
5082              
5083             B>>
5084              
5085             Points can be added to an index at the rate of about 50,000 per second.
5086              
5087             =item
5088              
5089             B>>
5090              
5091             Typical searches returning values run at about 25,000 to 50,000 searches per
5092             second. Worst-case performance is under 50 searches per second and searches
5093             returning no results run at over 100,000 searches per second. The overhead of
5094             traversing the results is fairly negligable.
5095              
5096             Quick searches run at 120,000 to 150,000 searches per second. Actually doing
5097             anything with the results slows things down a lot. Including traversal of the
5098             results, a typical quick search runs at 40,000 to 100,000 searches per second
5099             with the worst-case being about 80 searches per second.
5100              
5101             If distances to the result points are not needed, quick searches are typically
5102             about 75% faster than normal ones albeit with about 5 times as many results
5103             being returned.
5104              
5105             =item
5106              
5107             B>>
5108              
5109             For the C> method run time correlates with the
5110             size of the bounding box with smaller bounding boxes typically yielding faster
5111             run times.
5112              
5113             A fairly typical search yielding about 50 results runs at about 10,000 searches
5114             per second in normal mode and about 30,000 searches per second in quick mode.
5115             A nearly worst case example is a search returning 100,000 points; this will run
5116             at about 5 searches per second in normal mode or about 8,000 searches per second
5117             in quick mode.
5118              
5119              
5120             =item
5121              
5122             B>>
5123              
5124             For the Closest(...) method the highest performance is seen when there are
5125             result points close to the search point. Search speeds for the single closest
5126             point are typically in excess of 20,000 per second for close-by results or
5127             about 8,000 per second when results are far away. Worst case speeds of about
5128             1,000 searches per second occur when all indexed points are in the hemisphere
5129             opposite the search point.
5130              
5131             =item
5132              
5133             B>>
5134              
5135             For the Farthest(...) method the highest performance is seen when there are
5136             result points nearly antipodal to the search point. Search speeds for the
5137             single farthest point are typically in excess of 20,000 per second when
5138             nearly-antipodal points exist. Worst case speeds of about 1,000 searches per
5139             second occur when all indexed points are in the same hemisphere as the search
5140             point.
5141              
5142             =back
5143              
5144             Note that the numbers above are approximate and are highly dependant on the
5145             data being searched, the type of search being run, and on the number of results
5146             returned. Actual searches may be an order of magnitude or more slower.
5147              
5148             A sample benchmark run can be found in C To run the
5149             benchmarks yourself you can run C It needs the
5150             Devel::Size and Time::HiRes modules installed and a single run takes about 8
5151             minutes.
5152              
5153             Since Perl constants cannot be changed from the commandline you will need to
5154             edit the C to force the use of numeric keys, packed numeric keys,
5155             or text keys. This can be done by uncommenting the appropriate lines at the
5156             head of the file (look for C and C). When
5157             running C, the various other options can be found at the top of
5158             the script. When writing your own programs you can switch between the Perl and
5159             C single- or double-precision low-level code by using the C
5160             option when calling C.
5161              
5162             =head2 Potential optimizations
5163              
5164             The high cost of constructing and traversing the results seems inherent to the
5165             Perl language and there does not seem to be any way to avoid it. The is some
5166             potential for optimization though:
5167              
5168             =over
5169              
5170             =item *
5171              
5172             The C and C function calls might be sped up by
5173             assigning them to function handles (much as is done with C,
5174             etc.) instead of making the calls by dereferencing the variables.
5175              
5176             =item *
5177              
5178             Performance could potentially be improved by splitting the current combined
5179             index into individual indices for each level. Having smaller keys and indices
5180             should result in higher performance but the additional layer of indirection
5181             could slow things down in some circumstances.
5182              
5183             =item *
5184              
5185             Improvements might be possible to the performance of C, ...)>
5186             where S>E1> and per-point distances are not needed by the caller. At
5187             each iteration of the algorithm the previously-used search radius gives the
5188             maximal distance to all points already found, obviating the need to calculate
5189             every point's distance. Only points in the final level of the search would need
5190             to have their distances calculated. The downside to this method is that the
5191             point distances would not be available for all points in a result set (only for
5192             those found in the final search level).
5193              
5194             =item *
5195              
5196             A number of alternative datastructures were explored for the point index but
5197             benchmarking showed plain Perl hashes to be by far the most efficient. It is
5198             possible, though in my opinion unlikely, that a faster data structure choice
5199             exists that is suitable for use in this module.
5200              
5201             =back
5202              
5203             =head1 THEORY OF OPERATION
5204              
5205             =head2 Overview
5206              
5207             A given index comprises sets of tiles at various zoom levels with each tile
5208             containing a list of the points that lie within it. The lowest level of the
5209             index covers the entire globe. Each higher index level contains twice as many
5210             tiles in each direction. At each zoom level points are linearly mapped to
5211             grid tiles based on their latitudes and longitudes using an equirectangular
5212             projection. This is fairly analogous to how typical web slippy maps are
5213             organized (though they use a pseudo-mercator projection).
5214              
5215             As one approaches the poles the tiles become increasingly distorted with the
5216             area (in square meters) covered by each tile becoming progressively smaller.
5217             The distance in meters for one degree of longitude gets smaller as one moves
5218             away from the equator. The distance in meters for one degree of latitude,
5219             however, remains constant at all latitudes.
5220              
5221             Each tile has a name that gives its zoom level and position. These names are
5222             used as keys into a Perl hash allowing the quick retrieval of the points that
5223             lie within a given tile. The various search methods are designed to efficiently
5224             pull points from this index and filter them in various ways. The format used
5225             for the keys is described in the B section below.
5226              
5227             Additional datastructures (e.g. the list of all points in the index) are also
5228             present but knowing their details is not needed to understand how the index
5229             functions. In the descriptions below, some minor (but often critical) details
5230             have been omitted and some simplifications have been made; these details (mainly
5231             edge cases) are discussed in the code comments.
5232              
5233              
5234             =head2 Populating the index
5235              
5236             When a point is added to the index it is stored multiple times in the index
5237             hash, once for each level. This is done as follows:
5238              
5239             =over
5240              
5241             =item *
5242              
5243             The point's latitude and longitude are converted to integers. This is done
5244             using a simple linear mapping. In pseudo-code, the equations used are:
5245              
5246             max_size = 2**levels
5247            
5248             integer_latitude = floor( ( latitude + 90.0 ) * max_size / 180.0 )
5249             integer_latitude = max_size - 1 if (integer_latitude == max_size)
5250            
5251             integer_longitude = floor( ( longitude + 180.0 ) * max_size / 360.0 ) % max_size
5252            
5253             The values for C and C are in degrees and C is
5254             the number of levels in the index (not counting the S
5255              
5256             =item *
5257              
5258             Each index level is looped through from the index's maximum level to C<0>. At
5259             each level, the key (comprised of C, C, and
5260             C, see also below) is used to retrieve the corresponding
5261             value from the index hash. This value is a reference to the list of points that
5262             lie within the grid tile named by the key. The point being indexed is added to
5263             the retrieved list. If there is no list stored in the index for the current key
5264             then a new list is created and added. As a special case, all points adjacent to
5265             the poles (that is points with integer latitudes of C<0> or C) use
5266             the longitude C in their keys.
5267              
5268             Once the point has been added, the integer latitudes and longitudes as well as
5269             the C are shifted right by one bit in preparation for the the next
5270             level.
5271              
5272             =item *
5273              
5274             Once a the point has been added to the index at each level, the point is added to
5275             the global index entry using the key S, C, C>. (All indexed
5276             points can be found under this key.)
5277              
5278             =back
5279              
5280             =head2 Basic searching
5281              
5282             The C> method is typically used to find all points lying
5283             within a given radius of a search point. Two steps are performed by this
5284             method: retrieval of preliminary results and filtering of the results based on
5285             the search criteria.
5286              
5287             If no search radius was specified, if a global search was requested, or if the
5288             search radius covers more than half the globe then the preliminary results are
5289             all points in the index. Otherwise, the preliminary results are gathered as
5290             follows:
5291              
5292             The appropriate tile zoom level to use is determined using:
5293              
5294             shift = ceil( log2( search_radius / half_circumference ) )
5295             level = max_level - shift
5296              
5297             This results in the smallest level that could potentially contain all result
5298             points within a single tile.
5299              
5300             The search radius (in meters) is converted to two search angular radii, one for
5301             latitude and one for longitude. This is done since the number of meters per
5302             degree longitude decreases as one approaches the poles. Thus the north-south
5303             (latitude) search radius remains constant at all latitudes but the east-west
5304             (longitude) search radius increases as one nears the poles.
5305              
5306             Each extreme is converted to an integer and shifted right by the determined
5307             C, The preliminary results are retrieved from the index by iterating
5308             over the keys for the computed level, bounded by the integer extrema. This
5309             typically, but not always, results in a list of pointers to four tiles' points.
5310              
5311             If the C option is active then this preliminary list of lists of
5312             points is returned. If not then the points are filtered to only include those
5313             matching the search criteria. The filtered points are optionally sorted and
5314             then returned. Note that when large numbers of points have been found this
5315             filtering can be very slow; see B> above for details.
5316              
5317             =head2 Proximity searching
5318              
5319             The C> and C> methods find the points
5320             closest to (or farthest from) a search point. The C> method
5321             works as follows:
5322              
5323             The search starts at the most detailed level of the index and proceeds to the
5324             least detailed (C<0>). At each level, the grid tile that the search point lies
5325             in along with the three closest grid squares are identified. The method used
5326             for selecting the adjacent tiles is to look at the least-significant bits of the
5327             integer position at the previous (more detailed) level. A C<1> bit for the
5328             latitude selects tiles to the north, a C<0> bit the ones to the south. Likewise
5329             a C<1> for the longitude selects the ones east and a C<0> the ones west.
5330              
5331             Now that the four tiles have been identified, the largest radius from the search
5332             point to the tile edges is determined. The distance from the search point to
5333             each point within the four tiles is measured. If the point is within the radius
5334             computed and it passes any pre- or post-condition tests it is added to the
5335             results list. To speed up processing, points that have already been rejected
5336             along with the distances so far measured are cached. As a convenience, by
5337             default a filter is applied that omits the search point from the results.
5338              
5339             Once all points within a level's four chosen tiles have been gathered a check is
5340             done to see whether at least the requested number of points have been found. If
5341             they have then the loop ends, if not then the next (less-detailed) level is
5342             processed.
5343              
5344             By default, the results are then sorted by distance which ensures that the
5345             closest results are earliest in the list. This is necessary since although the
5346             nature of the algorithm tends to place closer points earlier in the results
5347             there is no inherent order to the points added from a particular index level.
5348             Lastly, the requested number of result points is returned.
5349              
5350             The C> method is largely implemented as a wrapper for
5351             C>. It functions by finding the closest points to the search
5352             point's antipode.
5353              
5354             =head2 Searching by bounding box
5355              
5356             The C> method works much the same as C>
5357             method. Instead of computing extrema of a search circle, those of the supplied
5358             bounding box are used. The tile level used is C, I )>
5359             where I> and I> are the most detailed levels that
5360             could potentially (given their extrema's angular distances) contain their
5361             respective extrema within a single tile in each direction. The remainder of the
5362             method is identical to that of C> albeit with all
5363             distance-related code removed.
5364              
5365             =head2 Tile naming (key generation)
5366              
5367             As mentioned earlier, keys consist of a zoom level, a latitude, and a longitude.
5368             Each key uniquely names a given tile.
5369              
5370             Zoom levels are either integers between C<0> and the maximum zoom level or the
5371             special zoom level C (with the value C<-1>) that covers the entire globe.
5372             Latitudes and longitudes are integers between C<0> and one less than the maximum
5373             grid size for the level. The tiles immediately adjacent to the poles are
5374             treated differently. In these areas the coverage of each tile is quite small
5375             and the algorithm around the poles would normally be complex. To accommodate
5376             these issues, the special value C (with the value C<-1>) is used for the
5377             longitude of the polar tiles (those areas with the lowest or highest latitude
5378             value for the key's level). All points lying in a polar region are assigned to
5379             that region's overarching C tile. At the global level all three components
5380             are set to C.
5381              
5382             If Perl has been compiled with 64-bit support then each key is packed into a
5383             64 bit integer. The level is stored in the upper 6 bits (bits 58 .. 63), the
5384             integer latitude in the next 29 bits (bits 29 .. 57), and the integer longitude
5385             in the low 29 bits (bits 0 .. 28). To represent the C value all bits in
5386             the relevant field are set to C<1>. Note that even on 32-bit systems Perl is
5387             often compiled with 64-bit support.
5388              
5389             If Perl does not have 64-bit support then a different format is used. In most
5390             places within Geo::Index, keys are stored as three-element array references.
5391             The first field contains the level, the second the integer latitude and the
5392             third the integer longitude. If present, C values are stored as-is as
5393             their integer value (C<-1>). For accessing the index, keys are converted to
5394             strings with the format "IC<:>IC<,>I" with the
5395             literal string "C" being used for C values.
5396              
5397             =head2 Object structure
5398              
5399             Each index object is a hash containing a number of entries These are:
5400              
5401             =over
5402              
5403             B{index}>> - The points index
5404              
5405             Entry is a hash reference. Keys are tile names (as discussed above), values are lists of
5406             point references.
5407              
5408             B{indices}>> - Indices used for each point
5409              
5410             Entry is a hash reference. Keys are point references, values are lists of tile names.
5411              
5412             B{positions}>> - Each point's position when indexed
5413              
5414             Entry is a hash reference. Keys are point references, values are two-element lists giving
5415             each point's latitude and longitude at the time it was indexed.
5416              
5417             B{levels}>> - Number of levels in the index (excluding the
5418             global level)
5419              
5420             B{max_level}>> - The highest-resolution level number (i.e. C - 1)
5421              
5422             B{max_size}>> - Number of grid tiles in each direction at most detailed level of index
5423              
5424             B{planetary_radius}>> - The planetary radius used by the index
5425             (in meters)
5426              
5427             B{polar_circumference}>> - The polar circumference used by the
5428             index (in meters)
5429              
5430             B{equatorial_circumference}>> - The equatorial circumference
5431             used by the index (in meters)
5432              
5433             =back
5434              
5435              
5436              
5437             =head1 BUGS AND DEFICIENCIES
5438              
5439             =head3 Known issues
5440              
5441             =over
5442              
5443             =item * This module is not believed to be thread-safe. In specific:
5444              
5445             =over
5446              
5447             =item *
5448             The C function stores the first point's position in
5449             global variables.
5450              
5451             To fix this, C and C would
5452             need to be removed plus C and
5453             C would need to be combined into a single
5454             4-argument C function. Calling code would need
5455             to be modified as appropriate. In terms of performance, the overall cost of
5456             doing this is likely quite low.
5457              
5458             =item *
5459             The search code stores distances computed for a specific search into the point
5460             datastructures. If multiple concurrent searches are run against a single index
5461             then distances computed by one search may overwrite those from another search.
5462             This can lead to inconsistent results.
5463              
5464             To fix this a per-search distance hash would need to be maintained. This could
5465             have serious performance implications and would preclude returning the point
5466             distances within the point hashes. The distances could, however, be returned in
5467             an additional datastructure.
5468              
5469             =item *
5470             Adding and deleting points to/from the index is not atomic. Running e.g. a
5471             search while points are being added or deleted can lead to unpredictable
5472             behavior (up to and including the program crashing).
5473              
5474             One could fix this by adding object-level locks:
5475              
5476             =over
5477              
5478             =item * Block concurrent calls to the C and Cmethods
5479              
5480             =item * Block calls to the C and C methods while searches are running
5481              
5482             =item * Block calls to C I when the C or C methods are active
5483              
5484             =back
5485              
5486             =back
5487              
5488             =back
5489              
5490             =over
5491              
5492             =item *
5493             Including the same point in multiple indices or searches at the same time could
5494             lead to interesting results.
5495              
5496             As mentioned above, this is due to the storage of search result distances within
5497             the points and not within the index object. Each search that involves a given
5498             point will likely overwrite its C value.
5499              
5500             This could be encountered in a number of ways. For example, a search using a
5501             condition function that itself runs a search against the second index could be
5502             problematic. This could be encountered even when using a single index. For
5503             example, if code relies on the distances values from a search it should save a
5504             copy of them before running subsearches against the same set of points. If this
5505             is not done then the distance values from the first search may be overwritten by
5506             those of the subsequent searches.
5507              
5508             =item *
5509              
5510             Geo::Index uses the spherical haversine formula to compute distances. While
5511             quite fast, its accuracy over long distances is poor, with a worst case error
5512             of about 0.1% (22 km). Since the module already has provision for changing the
5513             backends used for the distance methods, adding a new backend to, for example,
5514             compute accurate distances on e.g. a WGS-84 spheroid would be simple and
5515             straight-forward.
5516              
5517             =item *
5518              
5519             In places the code can be repetitious or awkward in style. This was done
5520             because, especially in the inner loops, speed has been favoured over clarity.
5521              
5522             =back
5523              
5524             =head3 Reporting bugs
5525              
5526             Please submit any bugs or feature requests either to C,
5527             through L,
5528             or through L. In any case I will
5529             receive notification when you do and you will be automatically notified of progress
5530             on your submission as it takes place. Any other comments can be sent to C.
5531              
5532             =head1 VERSION HISTORY
5533              
5534             B<0.0.8> (2019-04-10) - Corrected internal links, CPAN release
5535              
5536             =over
5537              
5538             =item * Corrected POD links to use spaces instead of underscores
5539              
5540             =item * Uploaded to CPAN
5541              
5542             =back
5543              
5544              
5545              
5546             B<0.0.7> (2019-04-08) - Methods can now be called using snake case
5547              
5548             =over
5549              
5550             =item * Added method aliases so that CamelCase methods can also be called using snake_case.
5551              
5552             =back
5553              
5554              
5555              
5556             B<0.0.6> (2019-04-05) - Bug fixes, additional tests
5557              
5558             =over
5559              
5560             =item * B>: Fixed off-by-one error at the north pole
5561              
5562             This affected B>> and B>>.
5563              
5564             =item * B>: Fixed off-by-one error at the north pole
5565              
5566             =item * More thorough tests
5567              
5568             =item * Improved documentation
5569              
5570             =back
5571              
5572              
5573              
5574             B<0.0.5> (2019-04-04) - Added methods, enhancements
5575              
5576             =over
5577              
5578             =item * B>>: New method
5579              
5580             =item * B>>: New method
5581              
5582             =item * B>>: Added C option
5583              
5584             =item * B>>: Added C option
5585              
5586             =back
5587              
5588              
5589              
5590             B<0.0.4> (2019-04-03) - Switched from Inline::C to XS for low-level C functions, minor restructuring
5591              
5592             =over
5593              
5594             =item * Low-level C code is now in C.
5595              
5596             All references to Inline::C have been removed.
5597              
5598             =item * Deprecated B> and replaced it with B>>
5599              
5600             =back
5601              
5602              
5603              
5604             B<0.0.3> (2019-04-01) - Added Vacuum(...), Sweep(...), and tests plus bug fixes and minor enhancements
5605              
5606             =over
5607              
5608             =item * B>>: New method
5609              
5610             =item * B>>: New method
5611              
5612             =item * Added tests
5613              
5614             =item * B>>: Bug fixes
5615              
5616             =item * Bnew(_..._)>>>: Added C option
5617              
5618             =back
5619              
5620              
5621             B<0.0.2> (2019-03-31) - Bug fixes and minor enhancements
5622              
5623             =over
5624              
5625             =item * B>>: Fixed bug for points added near (but not at) the north pole
5626              
5627             =item * B>>: Added C, C, and C values>
5628              
5629             =back
5630              
5631              
5632              
5633             B<0.0.1> (2018-03-30) - Initial release
5634              
5635             =head1 AUTHOR
5636              
5637             Alex Kent Hajnal S< > C S< > L
5638              
5639              
5640             =head1 COPYRIGHT
5641              
5642             Geo::Index
5643              
5644             Copyright 2019 Alexander Hajnal, All rights reserved
5645              
5646             This module is free software; you can redistribute it and/or modify it under
5647             the same terms as Perl itself. See L.
5648              
5649              
5650             =cut
5651              
5652              
5653             1;