Groonga 3.0.9 Source Code Document
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
geo.c
Go to the documentation of this file.
1 /* -*- c-basic-offset: 2 -*- */
2 /* Copyright(C) 2009-2012 Brazil
3 
4  This library is free software; you can redistribute it and/or
5  modify it under the terms of the GNU Lesser General Public
6  License version 2.1 as published by the Free Software Foundation.
7 
8  This library is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11  Lesser General Public License for more details.
12 
13  You should have received a copy of the GNU Lesser General Public
14  License along with this library; if not, write to the Free Software
15  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
16 */
17 
18 #include "geo.h"
19 #include "pat.h"
20 #include "util.h"
21 
22 #include <string.h>
23 #include <stdlib.h>
24 
25 typedef struct {
26  grn_id id;
27  double d;
28 } geo_entry;
29 
30 typedef struct
31 {
32  grn_geo_point key;
33  int key_size;
34 } mesh_entry;
35 
36 typedef struct {
37  grn_obj *pat;
38  grn_obj top_left_point_buffer;
39  grn_obj bottom_right_point_buffer;
40  grn_geo_point *top_left;
41  grn_geo_point *bottom_right;
42  grn_geo_point min;
43  grn_geo_point max;
44  int rectangle_common_bit;
45  uint8_t rectangle_common_key[sizeof(grn_geo_point)];
46 } in_rectangle_data;
47 
48 static int
49 compute_diff_bit(uint8_t *geo_key1, uint8_t *geo_key2)
50 {
51  int i, j, diff_bit = 0;
52 
53  for (i = 0; i < sizeof(grn_geo_point); i++) {
54  if (geo_key1[i] != geo_key2[i]) {
55  diff_bit = 8;
56  for (j = 0; j < 8; j++) {
57  if ((geo_key1[i] & (1 << (7 - j))) != (geo_key2[i] & (1 << (7 - j)))) {
58  diff_bit = j;
59  break;
60  }
61  }
62  break;
63  }
64  }
65 
66  return i * 8 + diff_bit;
67 }
68 
69 static void
70 compute_min_and_max_key(uint8_t *key_base, int diff_bit,
71  uint8_t *key_min, uint8_t *key_max)
72 {
73  int diff_byte, diff_bit_mask;
74 
75  diff_byte = diff_bit / 8;
76  diff_bit_mask = 0xff >> (diff_bit % 8);
77 
78  if (diff_byte == sizeof(grn_geo_point)) {
79  if (key_min) { memcpy(key_min, key_base, diff_byte); }
80  if (key_max) { memcpy(key_max, key_base, diff_byte); }
81  } else {
82  if (key_min) {
83  memcpy(key_min, key_base, diff_byte + 1);
84  key_min[diff_byte] &= ~diff_bit_mask;
85  memset(key_min + diff_byte + 1, 0,
86  sizeof(grn_geo_point) - diff_byte - 1);
87  }
88 
89  if (key_max) {
90  memcpy(key_max, key_base, diff_byte + 1);
91  key_max[diff_byte] |= diff_bit_mask;
92  memset(key_max + diff_byte + 1, 0xff,
93  sizeof(grn_geo_point) - diff_byte - 1);
94  }
95  }
96 }
97 
98 static void
99 compute_min_and_max(grn_geo_point *base_point, int diff_bit,
100  grn_geo_point *geo_min, grn_geo_point *geo_max)
101 {
102  uint8_t geo_key_base[sizeof(grn_geo_point)];
103  uint8_t geo_key_min[sizeof(grn_geo_point)];
104  uint8_t geo_key_max[sizeof(grn_geo_point)];
105 
106  grn_gton(geo_key_base, base_point, sizeof(grn_geo_point));
107  compute_min_and_max_key(geo_key_base, diff_bit,
108  geo_min ? geo_key_min : NULL,
109  geo_max ? geo_key_max : NULL);
110  if (geo_min) {
111  grn_ntog((uint8_t *)geo_min, geo_key_min, sizeof(grn_geo_point));
112  }
113  if (geo_max) {
114  grn_ntog((uint8_t *)geo_max, geo_key_max, sizeof(grn_geo_point));
115  }
116 }
117 
118 /* #define GEO_DEBUG */
119 
120 #ifdef GEO_DEBUG
121 #include <stdio.h>
122 
123 static void
124 inspect_mesh(grn_ctx *ctx, grn_geo_point *key, int key_size, int n)
125 {
126  grn_geo_point min, max;
127 
128  printf("mesh: %d:%d\n", n, key_size);
129 
130  printf("key: ");
131  grn_p_geo_point(ctx, key);
132 
133  compute_min_and_max(key, key_size, &min, &max);
134  printf("min: ");
135  grn_p_geo_point(ctx, &min);
136  printf("max: ");
137  grn_p_geo_point(ctx, &max);
138 }
139 
140 static void
141 inspect_mesh_entry(grn_ctx *ctx, mesh_entry *entries, int n)
142 {
143  mesh_entry *entry;
144 
145  entry = entries + n;
146  inspect_mesh(ctx, &(entry->key), entry->key_size, n);
147 }
148 
149 static void
150 inspect_tid(grn_ctx *ctx, grn_id tid, grn_geo_point *point, double d)
151 {
152  printf("tid: %d:%g", tid, d);
153  grn_p_geo_point(ctx, point);
154 }
155 
156 static void
157 inspect_key(grn_ctx *ctx, uint8_t *key)
158 {
159  int i;
160  for (i = 0; i < 8; i++) {
161  int j;
162  for (j = 0; j < 8; j++) {
163  printf("%d", (key[i] & (1 << (7 - j))) >> (7 - j));
164  }
165  printf(" ");
166  }
167  printf("\n");
168 }
169 
170 static void
171 print_key_mark(grn_ctx *ctx, int target_bit)
172 {
173  int i;
174 
175  for (i = 0; i < target_bit; i++) {
176  printf(" ");
177  if (i > 0 && i % 8 == 0) {
178  printf(" ");
179  }
180  }
181  if (i > 0 && i % 8 == 0) {
182  printf(" ");
183  }
184  printf("^\n");
185 }
186 
187 static void
188 inspect_cursor_entry(grn_ctx *ctx, grn_geo_cursor_entry *entry)
189 {
190  grn_geo_point point;
191 
192  printf("entry: ");
193  grn_ntog((uint8_t *)&point, entry->key, sizeof(grn_geo_point));
194  grn_p_geo_point(ctx, &point);
195  inspect_key(ctx, entry->key);
196  print_key_mark(ctx, entry->target_bit);
197 
198  printf(" target bit: %d\n", entry->target_bit);
199 
200 #define INSPECT_STATUS_FLAG(name) \
201  ((entry->status_flags & GRN_GEO_CURSOR_ENTRY_STATUS_ ## name) ? "true" : "false")
202 
203  printf(" top included: %s\n", INSPECT_STATUS_FLAG(TOP_INCLUDED));
204  printf("bottom included: %s\n", INSPECT_STATUS_FLAG(BOTTOM_INCLUDED));
205  printf(" left included: %s\n", INSPECT_STATUS_FLAG(LEFT_INCLUDED));
206  printf(" right included: %s\n", INSPECT_STATUS_FLAG(RIGHT_INCLUDED));
207  printf(" latitude inner: %s\n", INSPECT_STATUS_FLAG(LATITUDE_INNER));
208  printf("longitude inner: %s\n", INSPECT_STATUS_FLAG(LONGITUDE_INNER));
209 
210 #undef INSPECT_STATUS_FLAG
211 }
212 
213 static void
214 inspect_cursor_entry_targets(grn_ctx *ctx, grn_geo_cursor_entry *entry,
215  uint8_t *top_left_key, uint8_t *bottom_right_key,
216  grn_geo_cursor_entry *next_entry0,
217  grn_geo_cursor_entry *next_entry1)
218 {
219  printf("entry: ");
220  inspect_key(ctx, entry->key);
221  printf("top-left: ");
222  inspect_key(ctx, top_left_key);
223  printf("bottom-right: ");
224  inspect_key(ctx, bottom_right_key);
225  printf("next-entry-0: ");
226  inspect_key(ctx, next_entry0->key);
227  printf("next-entry-1: ");
228  inspect_key(ctx, next_entry1->key);
229  printf(" ");
230  print_key_mark(ctx, entry->target_bit + 1);
231 }
232 #else
233 # define inspect_mesh(...)
234 # define inspect_mesh_entry(...)
235 # define inspect_tid(...)
236 # define inspect_key(...)
237 # define print_key_mark(...)
238 # define inspect_cursor_entry(...)
239 # define inspect_cursor_entry_targets(...)
240 #endif
241 
242 static int
243 grn_geo_table_sort_detect_far_point(grn_ctx *ctx, grn_obj *table, grn_obj *index,
244  grn_pat *pat, geo_entry *entries,
245  grn_pat_cursor *pc, int n,
246  grn_bool accessorp,
247  grn_geo_point *base_point,
248  double *d_far, int *diff_bit)
249 {
250  int i = 0, diff_bit_prev, diff_bit_current;
251  grn_id tid;
252  geo_entry *ep, *p;
253  double d;
254  uint8_t geo_key_prev[sizeof(grn_geo_point)];
255  uint8_t geo_key_curr[sizeof(grn_geo_point)];
256  grn_geo_point point;
257 
258  *d_far = 0.0;
259  grn_gton(geo_key_curr, base_point, sizeof(grn_geo_point));
260  *diff_bit = sizeof(grn_geo_point) * 8;
261  diff_bit_current = sizeof(grn_geo_point) * 8;
262  memcpy(&point, base_point, sizeof(grn_geo_point));
263  ep = entries;
264  inspect_mesh(ctx, &point, *diff_bit, -1);
265  while ((tid = grn_pat_cursor_next(ctx, pc))) {
266  grn_ii_cursor *ic = grn_ii_cursor_open(ctx, (grn_ii *)index, tid, 0, 0, 1, 0);
267  if (ic) {
268  grn_ii_posting *posting;
269  grn_gton(geo_key_prev, &point, sizeof(grn_geo_point));
270  grn_pat_get_key(ctx, pat, tid, &point, sizeof(grn_geo_point));
271  grn_gton(geo_key_curr, &point, sizeof(grn_geo_point));
272  d = grn_geo_distance_rectangle_raw(ctx, base_point, &point);
273  inspect_tid(ctx, tid, &point, d);
274 
275  diff_bit_prev = diff_bit_current;
276  diff_bit_current = compute_diff_bit(geo_key_curr, geo_key_prev);
277 #ifdef GEO_DEBUG
278  printf("diff: %d:%d:%d\n", *diff_bit, diff_bit_prev, diff_bit_current);
279 #endif
280  if ((diff_bit_current % 2) == 1) {
281  diff_bit_current--;
282  }
283  if (diff_bit_current < diff_bit_prev && *diff_bit > diff_bit_current) {
284  if (i == n) {
285  grn_ii_cursor_close(ctx, ic);
286  break;
287  }
288  *diff_bit = diff_bit_current;
289  }
290 
291  if (d > *d_far) {
292  *d_far = d;
293  }
294  while ((posting = grn_ii_cursor_next(ctx, ic))) {
295  grn_id rid = accessorp
296  ? grn_table_get(ctx, table, &posting->rid, sizeof(grn_id))
297  : posting->rid;
298  if (rid) {
299  for (p = ep; entries < p && p[-1].d > d; p--) {
300  p->id = p[-1].id;
301  p->d = p[-1].d;
302  }
303  p->id = rid;
304  p->d = d;
305  if (i < n) {
306  ep++;
307  i++;
308  }
309  }
310  }
311  grn_ii_cursor_close(ctx, ic);
312  }
313  }
314 
315  return i;
316 }
317 
318 typedef enum {
319  MESH_LEFT_TOP,
320  MESH_RIGHT_TOP,
321  MESH_RIGHT_BOTTOM,
322  MESH_LEFT_BOTTOM
323 } mesh_position;
324 
325 /*
326  meshes should have
327  86 >= spaces when include_base_point_hash == GRN_FALSE,
328  87 >= spaces when include_base_point_hash == GRN_TRUE.
329 */
330 static int
331 grn_geo_get_meshes_for_circle(grn_ctx *ctx, grn_geo_point *base_point,
332  double d_far, int diff_bit,
333  int include_base_point_mesh,
334  mesh_entry *meshes)
335 {
336  double d;
337  int n_meshes;
338  int lat_diff, lng_diff;
339  mesh_position position;
340  grn_geo_point geo_base, geo_min, geo_max;
341 
342  compute_min_and_max(base_point, diff_bit - 2, &geo_min, &geo_max);
343 
344  lat_diff = (geo_max.latitude - geo_min.latitude + 1) / 2;
345  lng_diff = (geo_max.longitude - geo_min.longitude + 1) / 2;
346  geo_base.latitude = geo_min.latitude + lat_diff;
347  geo_base.longitude = geo_min.longitude + lng_diff;
348  if (base_point->latitude >= geo_base.latitude) {
349  if (base_point->longitude >= geo_base.longitude) {
350  position = MESH_RIGHT_TOP;
351  } else {
352  position = MESH_LEFT_TOP;
353  }
354  } else {
355  if (base_point->longitude >= geo_base.longitude) {
356  position = MESH_RIGHT_BOTTOM;
357  } else {
358  position = MESH_LEFT_BOTTOM;
359  }
360  }
361  /*
362  base_point: b
363  geo_min: i
364  geo_max: a
365  geo_base: x: must be at the left bottom in the top right mesh.
366 
367  e.g.: base_point is at the left bottom mesh case:
368  +------+------+
369  | | a|
370  | |x |
371  ^+------+------+
372  || | |
373  lng_diff || b | |
374  \/i------+------+
375  <------>
376  lat_diff
377 
378  grn_min + lat_diff -> the right mesh.
379  grn_min + lng_diff -> the top mesh.
380  */
381 #ifdef GEO_DEBUG
382  grn_p_geo_point(ctx, base_point);
383  printf("base: ");
384  grn_p_geo_point(ctx, &geo_base);
385  printf("min: ");
386  grn_p_geo_point(ctx, &geo_min);
387  printf("max: ");
388  grn_p_geo_point(ctx, &geo_max);
389  printf("diff: %d (%d, %d)\n", diff_bit, lat_diff, lng_diff);
390  switch (position) {
391  case MESH_LEFT_TOP :
392  printf("position: left-top\n");
393  break;
394  case MESH_RIGHT_TOP :
395  printf("position: right-top\n");
396  break;
397  case MESH_RIGHT_BOTTOM :
398  printf("position: right-bottom\n");
399  break;
400  case MESH_LEFT_BOTTOM :
401  printf("position: left-bottom\n");
402  break;
403  }
404 #endif
405 
406  n_meshes = 0;
407 
408 #define add_mesh(lat_diff_,lng_diff_,key_size_) do {\
409  meshes[n_meshes].key.latitude = geo_base.latitude + (lat_diff_);\
410  meshes[n_meshes].key.longitude = geo_base.longitude + (lng_diff_);\
411  meshes[n_meshes].key_size = key_size_;\
412  n_meshes++;\
413 } while (0)
414 
415  if (include_base_point_mesh || position != MESH_LEFT_TOP) {
416  add_mesh(0, -lng_diff, diff_bit);
417  }
418  if (include_base_point_mesh || position != MESH_RIGHT_TOP) {
419  add_mesh(0, 0, diff_bit);
420  }
421  if (include_base_point_mesh || position != MESH_RIGHT_BOTTOM) {
422  add_mesh(-lat_diff, 0, diff_bit);
423  }
424  if (include_base_point_mesh || position != MESH_LEFT_BOTTOM) {
425  add_mesh(-lat_diff, -lng_diff, diff_bit);
426  }
427 
428  /*
429  b: base_point
430  x: geo_base
431  0-83: sub meshes. 0-83 are added order.
432 
433  j: -5 -4 -3 -2 -1 0 1 2 3 4
434  +---+---+---+---+---+---+---+---+---+---+
435  |74 |75 |76 |77 |78 |79 |80 |81 |82 |83 | 4
436  +---+---+---+---+---+---+---+---+---+---+
437  |64 |65 |66 |67 |68 |69 |70 |71 |72 |73 | 3
438  +---+---+---+---+---+---+---+---+---+---+
439  |54 |55 |56 |57 |58 |59 |60 |61 |62 |63 | 2
440  +---+---+---+---+---+---+---+---+---+---+
441  |48 |49 |50 | b | |51 |52 |53 | 1
442  +---+---+---+ | +---+---+---+
443  |42 |43 |44 | |x |45 |46 |47 | 0
444  +---+---+---+-------+-------+---+---+---+
445  |36 |37 |38 | | |39 |40 |41 | -1
446  +---+---+---+ base meshes +---+---+---+
447  |30 |31 |32 | | |33 |34 |35 | -2
448  +---+---+---+---+---+---+---+---+---+---+
449  |20 |21 |22 |23 |24 |25 |26 |27 |28 |29 | -3
450  +---+---+---+---+---+---+---+---+---+---+
451  |10 |11 |12 |13 |14 |15 |16 |17 |18 |19 | -4
452  +---+---+---+---+---+---+---+---+---+---+
453  | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | -5
454  +---+---+---+---+---+---+---+---+---+---+
455  i
456  */
457  {
458  int i, j, n_sub_meshes, lat, lat_min, lat_max, lng, lng_min, lng_max;
459  n_sub_meshes = 0;
460  for (i = -5; i < 5; i++) {
461  lat_min = ((lat_diff + 1) / 2) * i;
462  lat_max = ((lat_diff + 1) / 2) * (i + 1) - 1;
463  for (j = -5; j < 5; j++) {
464  if (-3 < i && i < 2 && -3 < j && j < 2) {
465  continue;
466  }
467  lng_min = ((lng_diff + 1) / 2) * j;
468  lng_max = ((lng_diff + 1) / 2) * (j + 1) - 1;
469  if (base_point->latitude <= geo_base.latitude + lat_min) {
470  lat = geo_base.latitude + lat_min;
471  } else if (geo_base.latitude + lat_max < base_point->latitude) {
472  lat = geo_base.latitude + lat_max;
473  } else {
474  lat = base_point->latitude;
475  }
476  if (base_point->longitude <= geo_base.longitude + lng_min) {
477  lng = geo_base.longitude + lng_min;
478  } else if (geo_base.longitude + lng_max < base_point->longitude) {
479  lng = geo_base.longitude + lng_max;
480  } else {
481  lng = base_point->longitude;
482  }
483  meshes[n_meshes].key.latitude = lat;
484  meshes[n_meshes].key.longitude = lng;
485  d = grn_geo_distance_rectangle_raw(ctx, base_point,
486  &(meshes[n_meshes].key));
487  if (d < d_far) {
488 #ifdef GEO_DEBUG
489  printf("sub-mesh: %d: (%d,%d): (%d,%d;%d,%d)\n",
490  n_sub_meshes, base_point->latitude, base_point->longitude,
491  geo_base.latitude + lat_min,
492  geo_base.latitude + lat_max,
493  geo_base.longitude + lng_min,
494  geo_base.longitude + lng_max);
495  grn_p_geo_point(ctx, &(meshes[n_meshes].key));
496 #endif
497  meshes[n_meshes].key_size = diff_bit + 2;
498  n_meshes++;
499  }
500  n_sub_meshes++;
501  }
502  }
503  }
504 
505 #undef add_mesh
506 
507  return n_meshes;
508 }
509 
510 static int
511 grn_geo_table_sort_collect_points(grn_ctx *ctx, grn_obj *table, grn_obj *index,
512  grn_pat *pat,
513  geo_entry *entries, int n_entries,
514  int n, grn_bool accessorp,
515  grn_geo_point *base_point,
516  double d_far, int diff_bit)
517 {
518  int n_meshes;
519  mesh_entry meshes[86];
520  geo_entry *ep, *p;
521 
522  n_meshes = grn_geo_get_meshes_for_circle(ctx, base_point, d_far, diff_bit,
523  GRN_FALSE, meshes);
524 
525  ep = entries + n_entries;
526  while (n_meshes--) {
527  grn_id tid;
528  grn_pat_cursor *pc = grn_pat_cursor_open(ctx, pat,
529  &(meshes[n_meshes].key),
530  meshes[n_meshes].key_size,
531  NULL, 0,
532  0, -1,
533  GRN_CURSOR_PREFIX|GRN_CURSOR_SIZE_BY_BIT);
534  inspect_mesh_entry(ctx, meshes, n_meshes);
535  if (pc) {
536  while ((tid = grn_pat_cursor_next(ctx, pc))) {
537  grn_ii_cursor *ic = grn_ii_cursor_open(ctx, (grn_ii *)index, tid, 0, 0, 1, 0);
538  if (ic) {
539  double d;
540  grn_geo_point pos;
541  grn_ii_posting *posting;
542  grn_pat_get_key(ctx, pat, tid, &pos, sizeof(grn_geo_point));
543  d = grn_geo_distance_rectangle_raw(ctx, base_point, &pos);
544  inspect_tid(ctx, tid, &pos, d);
545  while ((posting = grn_ii_cursor_next(ctx, ic))) {
546  grn_id rid = accessorp
547  ? grn_table_get(ctx, table, &posting->rid, sizeof(grn_id))
548  : posting->rid;
549  if (rid) {
550  for (p = ep; entries < p && p[-1].d > d; p--) {
551  p->id = p[-1].id;
552  p->d = p[-1].d;
553  }
554  p->id = rid;
555  p->d = d;
556  if (n_entries < n) {
557  ep++;
558  n_entries++;
559  }
560  }
561  }
562  grn_ii_cursor_close(ctx, ic);
563  }
564  }
565  grn_pat_cursor_close(ctx, pc);
566  }
567  }
568  return n_entries;
569 }
570 
571 static inline grn_obj *
572 find_geo_sort_index(grn_ctx *ctx, grn_obj *key)
573 {
574  grn_obj *index = NULL;
575 
576  if (GRN_ACCESSORP(key)) {
577  grn_accessor *accessor = (grn_accessor *)key;
578  if (accessor->action != GRN_ACCESSOR_GET_KEY) {
579  return NULL;
580  }
581  if (!(DB_OBJ(accessor->obj)->id & GRN_OBJ_TMP_OBJECT)) {
582  return NULL;
583  }
584  if (accessor->obj->header.type != GRN_TABLE_HASH_KEY) {
585  return NULL;
586  }
587  if (!accessor->next) {
588  return NULL;
589  }
590  grn_column_index(ctx, accessor->next->obj, GRN_OP_LESS, &index, 1, NULL);
591  } else {
592  grn_column_index(ctx, key, GRN_OP_LESS, &index, 1, NULL);
593  }
594 
595  return index;
596 }
597 
598 static inline int
599 grn_geo_table_sort_by_distance(grn_ctx *ctx,
600  grn_obj *table,
601  grn_obj *index,
602  grn_pat *pat,
603  grn_pat_cursor *pc,
604  grn_bool accessorp,
605  grn_geo_point *base_point,
606  int offset,
607  int limit,
608  grn_obj *result)
609 {
610  int n_entries = 0, e = offset + limit;
611  geo_entry *entries;
612 
613  if ((entries = GRN_MALLOC(sizeof(geo_entry) * (e + 1)))) {
614  int n, diff_bit;
615  double d_far;
616  geo_entry *ep;
617  grn_bool need_not_indexed_records;
618  grn_hash *indexed_records = NULL;
619 
620  n = grn_geo_table_sort_detect_far_point(ctx, table, index, pat,
621  entries, pc, e, accessorp,
622  base_point,
623  &d_far, &diff_bit);
624  if (diff_bit > 0) {
625  n = grn_geo_table_sort_collect_points(ctx, table, index, pat,
626  entries, n, e, accessorp,
627  base_point, d_far, diff_bit);
628  }
629  need_not_indexed_records = offset + limit > n;
630  if (need_not_indexed_records) {
631  indexed_records = grn_hash_create(ctx, NULL, sizeof(grn_id), 0,
632  GRN_OBJ_TABLE_HASH_KEY|GRN_HASH_TINY);
633  }
634  for (ep = entries + offset;
635  n_entries < limit && ep < entries + n;
636  n_entries++, ep++) {
637  grn_id *sorted_id;
638  if (!grn_array_add(ctx, (grn_array *)result, (void **)&sorted_id)) {
639  if (indexed_records) {
640  grn_hash_close(ctx, indexed_records);
641  indexed_records = NULL;
642  }
643  break;
644  }
645  *sorted_id = ep->id;
646  if (indexed_records) {
647  grn_hash_add(ctx, indexed_records, &(ep->id), sizeof(grn_id),
648  NULL, NULL);
649  }
650  }
651  GRN_FREE(entries);
652  if (indexed_records) {
653  GRN_TABLE_EACH(ctx, table, GRN_ID_NIL, GRN_ID_MAX, id, NULL, NULL, NULL, {
654  if (!grn_hash_get(ctx, indexed_records, &id, sizeof(grn_id), NULL)) {
655  grn_id *sorted_id;
656  if (grn_array_add(ctx, (grn_array *)result, (void **)&sorted_id)) {
657  *sorted_id = id;
658  }
659  n_entries++;
660  if (n_entries == limit) {
661  break;
662  }
663  };
664  });
665  grn_hash_close(ctx, indexed_records);
666  }
667  }
668 
669  return n_entries;
670 }
671 
672 int
673 grn_geo_table_sort(grn_ctx *ctx, grn_obj *table, int offset, int limit,
674  grn_obj *result, grn_table_sort_key *keys, int n_keys)
675 {
676  grn_obj *index;
677  int i = 0, e = offset + limit;
678  grn_bool accessorp = GRN_ACCESSORP(keys->key);
679  if (n_keys == 2 && (index = find_geo_sort_index(ctx, keys->key))) {
680  grn_id tid;
681  grn_obj *arg = keys[1].key;
682  grn_pat *pat = (grn_pat *)grn_ctx_at(ctx, index->header.domain);
683  grn_id domain = pat->obj.header.domain;
684  grn_pat_cursor *pc = grn_pat_cursor_open(ctx, pat, NULL, 0,
685  GRN_BULK_HEAD(arg), GRN_BULK_VSIZE(arg),
686  0, -1, GRN_CURSOR_PREFIX);
687  if (pc) {
688  if (domain != GRN_DB_TOKYO_GEO_POINT && domain != GRN_DB_WGS84_GEO_POINT) {
689  while (i < e && (tid = grn_pat_cursor_next(ctx, pc))) {
690  grn_ii_cursor *ic = grn_ii_cursor_open(ctx, (grn_ii *)index, tid, 0, 0, 1, 0);
691  if (ic) {
692  grn_ii_posting *posting;
693  while (i < e && (posting = grn_ii_cursor_next(ctx, ic))) {
694  if (offset <= i) {
695  grn_id *v;
696  if (!grn_array_add(ctx, (grn_array *)result, (void **)&v)) { break; }
697  *v = posting->rid;
698  }
699  i++;
700  }
701  grn_ii_cursor_close(ctx, ic);
702  }
703  }
704  } else {
705  grn_geo_point *base_point = (grn_geo_point *)GRN_BULK_HEAD(arg);
706  i = grn_geo_table_sort_by_distance(ctx, table, index, pat,
707  pc, accessorp, base_point,
708  offset, limit, result);
709  }
710  grn_pat_cursor_close(ctx, pc);
711  }
712  }
713  return i;
714 }
715 
716 grn_rc
717 grn_geo_resolve_approximate_type(grn_ctx *ctx, grn_obj *type_name,
718  grn_geo_approximate_type *type)
719 {
720  grn_rc rc;
721  grn_obj approximate_type;
722 
723  GRN_TEXT_INIT(&approximate_type, 0);
724  rc = grn_obj_cast(ctx, type_name, &approximate_type, GRN_FALSE);
725  if (rc == GRN_SUCCESS) {
726  const char *name;
727  unsigned int size;
728  name = GRN_TEXT_VALUE(&approximate_type);
729  size = GRN_TEXT_LEN(&approximate_type);
730  if ((strncmp("rectangle", name, size) == 0) ||
731  (strncmp("rect", name, size) == 0)) {
732  *type = GRN_GEO_APPROXIMATE_RECTANGLE;
733  } else if ((strncmp("sphere", name, size) == 0) ||
734  (strncmp("sphr", name, size) == 0)) {
735  *type = GRN_GEO_APPROXIMATE_SPHERE;
736  } else if ((strncmp("ellipsoid", name, size) == 0) ||
737  (strncmp("ellip", name, size) == 0)) {
738  *type = GRN_GEO_APPROXIMATE_ELLIPSOID;
739  } else {
740  ERR(GRN_INVALID_ARGUMENT,
741  "geo distance approximate type must be one of "
742  "[rectangle, rect, sphere, sphr, ellipsoid, ellip]"
743  ": <%.*s>",
744  size, name);
745  }
746  }
747  GRN_OBJ_FIN(ctx, &approximate_type);
748 
749  return rc;
750 }
751 
752 typedef double (*grn_geo_distance_raw_func)(grn_ctx *ctx,
753  grn_geo_point *point1,
754  grn_geo_point *point2);
755 
756 grn_rc
757 grn_selector_geo_in_circle(grn_ctx *ctx, grn_obj *table, grn_obj *index,
758  int nargs, grn_obj **args,
759  grn_obj *res, grn_operator op)
760 {
761  grn_geo_approximate_type type = GRN_GEO_APPROXIMATE_RECTANGLE;
762 
763  if (!(nargs == 4 || nargs == 5)) {
764  ERR(GRN_INVALID_ARGUMENT,
765  "geo_in_circle(): requires 3 or 4 arguments but was <%d> arguments",
766  nargs - 1);
767  return ctx->rc;
768  }
769 
770  if (!index) {
771  grn_obj *point_column;
772  char column_name[GRN_TABLE_MAX_KEY_SIZE];
773  int column_name_size;
774  point_column = args[1];
775  column_name_size = grn_obj_name(ctx, point_column,
776  column_name, GRN_TABLE_MAX_KEY_SIZE);
777  ERR(GRN_INVALID_ARGUMENT,
778  "geo_in_circle(): index for <%.*s> is missing",
779  column_name_size, column_name);
780  return ctx->rc;
781  }
782 
783  if (nargs == 5) {
784  if (grn_geo_resolve_approximate_type(ctx, args[4], &type) != GRN_SUCCESS) {
785  return ctx->rc;
786  }
787  }
788 
789  {
790  grn_obj *center_point, *distance;
791  center_point = args[2];
792  distance = args[3];
793  grn_geo_select_in_circle(ctx, index, center_point, distance, type, res, op);
794  }
795 
796  return ctx->rc;
797 }
798 
799 static grn_geo_distance_raw_func
800 grn_geo_resolve_distance_raw_func (grn_ctx *ctx,
801  grn_geo_approximate_type approximate_type,
802  grn_id domain)
803 {
804  grn_geo_distance_raw_func distance_raw_func = NULL;
805 
806  switch (approximate_type) {
807  case GRN_GEO_APPROXIMATE_RECTANGLE :
808  distance_raw_func = grn_geo_distance_rectangle_raw;
809  break;
810  case GRN_GEO_APPROXIMATE_SPHERE :
811  distance_raw_func = grn_geo_distance_sphere_raw;
812  break;
813  case GRN_GEO_APPROXIMATE_ELLIPSOID :
814  if (domain == GRN_DB_WGS84_GEO_POINT) {
815  distance_raw_func = grn_geo_distance_ellipsoid_raw_wgs84;
816  } else {
817  distance_raw_func = grn_geo_distance_ellipsoid_raw_tokyo;
818  }
819  break;
820  default :
821  break;
822  }
823 
824  return distance_raw_func;
825 }
826 
827 grn_rc
828 grn_geo_select_in_circle(grn_ctx *ctx, grn_obj *index,
829  grn_obj *center_point, grn_obj *distance,
830  grn_geo_approximate_type approximate_type,
831  grn_obj *res, grn_operator op)
832 {
833  grn_id domain;
834  double center_longitude, center_latitude;
835  double d;
836  grn_obj *pat, *point_on_circle = NULL, center_point_, point_on_circle_;
837  grn_geo_point *center, on_circle;
838  grn_geo_distance_raw_func distance_raw_func;
839  pat = grn_ctx_at(ctx, index->header.domain);
840  domain = pat->header.domain;
841  if (domain != GRN_DB_TOKYO_GEO_POINT && domain != GRN_DB_WGS84_GEO_POINT) {
842  char name[GRN_TABLE_MAX_KEY_SIZE];
843  int name_size = 0;
844  grn_obj *domain_object;
845  domain_object = grn_ctx_at(ctx, domain);
846  if (domain_object) {
847  name_size = grn_obj_name(ctx, domain_object, name, GRN_TABLE_MAX_KEY_SIZE);
848  grn_obj_unlink(ctx, domain_object);
849  } else {
850  strcpy(name, "(null)");
851  name_size = strlen(name);
852  }
853  ERR(GRN_INVALID_ARGUMENT,
854  "geo_in_circle(): index table must be "
855  "TokyoGeoPoint or WGS84GeoPoint key type table: <%.*s>",
856  name_size, name);
857  goto exit;
858  }
859 
860  if (center_point->header.domain != domain) {
861  GRN_OBJ_INIT(&center_point_, GRN_BULK, 0, domain);
862  if (grn_obj_cast(ctx, center_point, &center_point_, GRN_FALSE)) { goto exit; }
863  center_point = &center_point_;
864  }
865  center = GRN_GEO_POINT_VALUE_RAW(center_point);
866  center_longitude = GRN_GEO_INT2RAD(center->longitude);
867  center_latitude = GRN_GEO_INT2RAD(center->latitude);
868 
869  distance_raw_func = grn_geo_resolve_distance_raw_func(ctx,
870  approximate_type,
871  domain);
872  if (!distance_raw_func) {
873  ERR(GRN_INVALID_ARGUMENT,
874  "unknown approximate type: <%d>", approximate_type);
875  goto exit;
876  }
877 
878  switch (distance->header.domain) {
879  case GRN_DB_INT32 :
880  d = GRN_INT32_VALUE(distance);
881  on_circle.latitude = center->latitude + GRN_GEO_RAD2INT(d / (double)GRN_GEO_RADIUS);
882  on_circle.longitude = center->longitude;
883  break;
884  case GRN_DB_UINT32 :
885  d = GRN_UINT32_VALUE(distance);
886  on_circle.latitude = center->latitude + GRN_GEO_RAD2INT(d / (double)GRN_GEO_RADIUS);
887  on_circle.longitude = center->longitude;
888  break;
889  case GRN_DB_INT64 :
890  d = GRN_INT64_VALUE(distance);
891  on_circle.latitude = center->latitude + GRN_GEO_RAD2INT(d / (double)GRN_GEO_RADIUS);
892  on_circle.longitude = center->longitude;
893  break;
894  case GRN_DB_UINT64 :
895  d = GRN_UINT64_VALUE(distance);
896  on_circle.latitude = center->latitude + GRN_GEO_RAD2INT(d / (double)GRN_GEO_RADIUS);
897  on_circle.longitude = center->longitude;
898  break;
899  case GRN_DB_FLOAT :
900  d = GRN_FLOAT_VALUE(distance);
901  on_circle.latitude = center->latitude + GRN_GEO_RAD2INT(d / (double)GRN_GEO_RADIUS);
902  on_circle.longitude = center->longitude;
903  break;
904  case GRN_DB_SHORT_TEXT :
905  case GRN_DB_TEXT :
906  case GRN_DB_LONG_TEXT :
907  GRN_OBJ_INIT(&point_on_circle_, GRN_BULK, 0, domain);
908  if (grn_obj_cast(ctx, point_on_circle, &point_on_circle_, GRN_FALSE)) { goto exit; }
909  point_on_circle = &point_on_circle_;
910  /* fallthru */
911  case GRN_DB_TOKYO_GEO_POINT :
912  case GRN_DB_WGS84_GEO_POINT :
913  if (domain != distance->header.domain) { /* todo */ goto exit; }
914  if (!point_on_circle) { point_on_circle = distance; }
915  GRN_GEO_POINT_VALUE(point_on_circle,
916  on_circle.latitude, on_circle.longitude);
917  d = distance_raw_func(ctx, center, &on_circle);
918  if (point_on_circle == &point_on_circle_) {
919  grn_obj_unlink(ctx, point_on_circle);
920  }
921  break;
922  default :
923  goto exit;
924  }
925  {
926  int n_meshes, diff_bit;
927  double d_far;
928  mesh_entry meshes[87];
929  uint8_t geo_key1[sizeof(grn_geo_point)];
930  uint8_t geo_key2[sizeof(grn_geo_point)];
931 
932  d_far = grn_geo_distance_rectangle_raw(ctx, center, &on_circle);
933  grn_gton(geo_key1, center, sizeof(grn_geo_point));
934  grn_gton(geo_key2, &on_circle, sizeof(grn_geo_point));
935  diff_bit = compute_diff_bit(geo_key1, geo_key2);
936 #ifdef GEO_DEBUG
937  printf("center point: ");
938  grn_p_geo_point(ctx, center);
939  printf("point on circle: ");
940  grn_p_geo_point(ctx, &on_circle);
941  printf("diff: %d\n", diff_bit);
942 #endif
943  if ((diff_bit % 2) == 1) {
944  diff_bit--;
945  }
946  n_meshes = grn_geo_get_meshes_for_circle(ctx, center,
947  d_far, diff_bit, GRN_TRUE,
948  meshes);
949  while (n_meshes--) {
950  grn_table_cursor *tc;
951  tc = grn_table_cursor_open(ctx, pat,
952  &(meshes[n_meshes].key),
953  meshes[n_meshes].key_size,
954  NULL, 0,
955  0, -1,
956  GRN_CURSOR_PREFIX|GRN_CURSOR_SIZE_BY_BIT);
957  inspect_mesh_entry(ctx, meshes, n_meshes);
958  if (tc) {
959  grn_id tid;
960  grn_geo_point point;
961  while ((tid = grn_table_cursor_next(ctx, tc))) {
962  double point_distance;
963  grn_table_get_key(ctx, pat, tid, &point, sizeof(grn_geo_point));
964  point_distance = distance_raw_func(ctx, &point, center);
965  if (point_distance <= d) {
966  inspect_tid(ctx, tid, &point, point_distance);
967  grn_ii_at(ctx, (grn_ii *)index, tid, (grn_hash *)res, op);
968  }
969  }
970  grn_table_cursor_close(ctx, tc);
971  }
972  }
973  }
974 exit :
975  grn_ii_resolve_sel_and(ctx, (grn_hash *)res, op);
976  return ctx->rc;
977 }
978 
979 grn_rc
980 grn_selector_geo_in_rectangle(grn_ctx *ctx, grn_obj *table, grn_obj *index,
981  int nargs, grn_obj **args,
982  grn_obj *res, grn_operator op)
983 {
984  if (nargs == 4) {
985  grn_obj *top_left_point, *bottom_right_point;
986  top_left_point = args[2];
987  bottom_right_point = args[3];
988  grn_geo_select_in_rectangle(ctx, index,
989  top_left_point, bottom_right_point,
990  res, op);
991  } else {
992  ERR(GRN_INVALID_ARGUMENT,
993  "geo_in_rectangle(): requires 3 arguments but was <%d> arguments",
994  nargs - 1);
995  }
996  return ctx->rc;
997 }
998 
999 static grn_rc
1000 in_rectangle_data_prepare(grn_ctx *ctx, grn_obj *index,
1001  grn_obj *top_left_point,
1002  grn_obj *bottom_right_point,
1003  const char *process_name,
1004  in_rectangle_data *data)
1005 {
1006  grn_id domain;
1007 
1008  if (!index) {
1009  ERR(GRN_INVALID_ARGUMENT, "%s: index column is missing", process_name);
1010  goto exit;
1011  }
1012 
1013  data->pat = grn_ctx_at(ctx, index->header.domain);
1014  domain = data->pat->header.domain;
1015  if (domain != GRN_DB_TOKYO_GEO_POINT && domain != GRN_DB_WGS84_GEO_POINT) {
1016  char name[GRN_TABLE_MAX_KEY_SIZE];
1017  int name_size = 0;
1018  grn_obj *domain_object;
1019  domain_object = grn_ctx_at(ctx, domain);
1020  if (domain_object) {
1021  name_size = grn_obj_name(ctx, domain_object, name, GRN_TABLE_MAX_KEY_SIZE);
1022  grn_obj_unlink(ctx, domain_object);
1023  } else {
1024  strcpy(name, "(null)");
1025  name_size = strlen(name);
1026  }
1027  ERR(GRN_INVALID_ARGUMENT,
1028  "%s: index table must be "
1029  "TokyoGeoPoint or WGS84GeoPoint key type table: <%.*s>",
1030  process_name, name_size, name);
1031  goto exit;
1032  }
1033 
1034  if (top_left_point->header.domain != domain) {
1035  grn_obj_reinit(ctx, &(data->top_left_point_buffer), domain, GRN_BULK);
1036  if (grn_obj_cast(ctx, top_left_point, &(data->top_left_point_buffer),
1037  GRN_FALSE)) {
1038  goto exit;
1039  }
1040  top_left_point = &(data->top_left_point_buffer);
1041  }
1042  data->top_left = GRN_GEO_POINT_VALUE_RAW(top_left_point);
1043 
1044  if (bottom_right_point->header.domain != domain) {
1045  grn_obj_reinit(ctx, &(data->bottom_right_point_buffer), domain, GRN_BULK);
1046  if (grn_obj_cast(ctx, bottom_right_point, &(data->bottom_right_point_buffer),
1047  GRN_FALSE)) {
1048  goto exit;
1049  }
1050  bottom_right_point = &(data->bottom_right_point_buffer);
1051  }
1052  data->bottom_right = GRN_GEO_POINT_VALUE_RAW(bottom_right_point);
1053 
1054  {
1055  grn_geo_point *top_left, *bottom_right;
1056 
1057  top_left = data->top_left;
1058  bottom_right = data->bottom_right;
1059 
1060  if (top_left->latitude < 0 || top_left->longitude < 0 ||
1061  bottom_right->latitude < 0 || bottom_right->longitude < 0) {
1062  ERR(GRN_FUNCTION_NOT_IMPLEMENTED,
1063  "%s: negative coordinate is not implemented.", process_name);
1064  goto exit;
1065  }
1066 
1067  if (top_left->latitude >= GRN_GEO_MAX_LATITUDE) {
1068  ERR(GRN_INVALID_ARGUMENT,
1069  "%s: top left point's latitude is too big: "
1070  "<%d>(max:%d): (%d,%d) (%d,%d)",
1071  process_name,
1072  GRN_GEO_MAX_LATITUDE, top_left->latitude,
1073  top_left->latitude, top_left->longitude,
1074  bottom_right->latitude, bottom_right->longitude);
1075  goto exit;
1076  }
1077 
1078  if (top_left->longitude >= GRN_GEO_MAX_LONGITUDE) {
1079  ERR(GRN_INVALID_ARGUMENT,
1080  "%s: top left point's longitude is too big: "
1081  "<%d>(max:%d): (%d,%d) (%d,%d)",
1082  process_name,
1083  GRN_GEO_MAX_LONGITUDE, top_left->longitude,
1084  top_left->latitude, top_left->longitude,
1085  bottom_right->latitude, bottom_right->longitude);
1086  goto exit;
1087  }
1088 
1089  if (bottom_right->latitude >= GRN_GEO_MAX_LATITUDE) {
1090  ERR(GRN_INVALID_ARGUMENT,
1091  "%s: bottom right point's latitude is too big: "
1092  "<%d>(max:%d): (%d,%d) (%d,%d)",
1093  process_name,
1094  GRN_GEO_MAX_LATITUDE, bottom_right->latitude,
1095  top_left->latitude, top_left->longitude,
1096  bottom_right->latitude, bottom_right->longitude);
1097  goto exit;
1098  }
1099 
1100  if (bottom_right->longitude >= GRN_GEO_MAX_LONGITUDE) {
1101  ERR(GRN_INVALID_ARGUMENT,
1102  "%s: bottom right point's longitude is too big: "
1103  "<%d>(max:%d): (%d,%d) (%d,%d)",
1104  process_name,
1105  GRN_GEO_MAX_LONGITUDE, bottom_right->longitude,
1106  top_left->latitude, top_left->longitude,
1107  bottom_right->latitude, bottom_right->longitude);
1108  goto exit;
1109  }
1110  }
1111 
1112  {
1113  int latitude_distance, longitude_distance;
1114  int diff_bit;
1115  grn_geo_point base;
1116  grn_geo_point *top_left, *bottom_right;
1117  grn_geo_point *geo_point_input;
1118  uint8_t geo_key_input[sizeof(grn_geo_point)];
1119  uint8_t geo_key_base[sizeof(grn_geo_point)];
1120  uint8_t geo_key_top_left[sizeof(grn_geo_point)];
1121  uint8_t geo_key_bottom_right[sizeof(grn_geo_point)];
1122 
1123  top_left = data->top_left;
1124  bottom_right = data->bottom_right;
1125 
1126  latitude_distance = top_left->latitude - bottom_right->latitude;
1127  longitude_distance = bottom_right->longitude - top_left->longitude;
1128  if (latitude_distance > longitude_distance) {
1129  geo_point_input = bottom_right;
1130  base.latitude = bottom_right->latitude;
1131  base.longitude = bottom_right->longitude - longitude_distance;
1132  } else {
1133  geo_point_input = top_left;
1134  base.latitude = top_left->latitude - latitude_distance;
1135  base.longitude = top_left->longitude;
1136  }
1137  grn_gton(geo_key_input, geo_point_input, sizeof(grn_geo_point));
1138  grn_gton(geo_key_base, &base, sizeof(grn_geo_point));
1139  diff_bit = compute_diff_bit(geo_key_input, geo_key_base);
1140  compute_min_and_max(&base, diff_bit, &(data->min), &(data->max));
1141 
1142  grn_gton(geo_key_top_left, top_left, sizeof(grn_geo_point));
1143  grn_gton(geo_key_bottom_right, bottom_right, sizeof(grn_geo_point));
1144  data->rectangle_common_bit =
1145  compute_diff_bit(geo_key_top_left, geo_key_bottom_right) - 1;
1146  compute_min_and_max_key(geo_key_top_left, data->rectangle_common_bit + 1,
1147  data->rectangle_common_key, NULL);
1148 
1149 #ifdef GEO_DEBUG
1150  printf("base: ");
1151  grn_p_geo_point(ctx, &base);
1152  printf("min: ");
1153  grn_p_geo_point(ctx, &(data->min));
1154  printf("max: ");
1155  grn_p_geo_point(ctx, &(data->max));
1156  printf("top-left: ");
1157  grn_p_geo_point(ctx, top_left);
1158  printf("bottom-right: ");
1159  grn_p_geo_point(ctx, bottom_right);
1160  printf("rectangle-common-bit:%10d\n", data->rectangle_common_bit);
1161  printf("distance(latitude): %10d\n", latitude_distance);
1162  printf("distance(longitude): %10d\n", longitude_distance);
1163 #endif
1164  }
1165 
1166 exit :
1167  return ctx->rc;
1168 }
1169 
1170 #define SAME_BIT_P(a, b, n_bit)\
1171  ((((uint8_t *)(a))[(n_bit) / 8] & (1 << (7 - ((n_bit) % 8)))) ==\
1172  (((uint8_t *)(b))[(n_bit) / 8] & (1 << (7 - ((n_bit) % 8)))))
1173 
1174 #define CURSOR_ENTRY_UPDATE_STATUS(entry, name, other_key) do {\
1175  if (SAME_BIT_P((entry)->key, (other_key), (entry)->target_bit)) {\
1176  (entry)->status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_ ## name;\
1177  } else {\
1178  (entry)->status_flags &= ~GRN_GEO_CURSOR_ENTRY_STATUS_ ## name;\
1179  }\
1180 } while (0)
1181 
1182 #define CURSOR_ENTRY_CHECK_STATUS(entry, name)\
1183  ((entry)->status_flags & GRN_GEO_CURSOR_ENTRY_STATUS_ ## name)
1184 #define CURSOR_ENTRY_IS_INNER(entry)\
1185  (((entry)->status_flags &\
1186  (GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER |\
1187  GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER)) ==\
1188  (GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER |\
1189  GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER))
1190 #define CURSOR_ENTRY_INCLUDED_IN_LATITUDE_DIRECTION(entry)\
1191  ((entry)->status_flags &\
1192  (GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER |\
1193  GRN_GEO_CURSOR_ENTRY_STATUS_TOP_INCLUDED |\
1194  GRN_GEO_CURSOR_ENTRY_STATUS_BOTTOM_INCLUDED))
1195 #define CURSOR_ENTRY_INCLUDED_IN_LONGITUDE_DIRECTION(entry)\
1196  ((entry)->status_flags &\
1197  (GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER |\
1198  GRN_GEO_CURSOR_ENTRY_STATUS_LEFT_INCLUDED |\
1199  GRN_GEO_CURSOR_ENTRY_STATUS_RIGHT_INCLUDED))
1200 
1201 #define SET_N_BIT(a, n_bit)\
1202  (((uint8_t *)(a))[((n_bit) / 8)] ^= (1 << (7 - ((n_bit) % 8))))
1203 
1204 #define N_BIT(a, n_bit)\
1205  ((((uint8_t *)(a))[((n_bit) / 8)] &\
1206  (1 << (7 - ((n_bit) % 8)))) >> (1 << (7 - ((n_bit) % 8))))
1207 
1208 grn_obj *
1209 grn_geo_cursor_open_in_rectangle(grn_ctx *ctx,
1210  grn_obj *index,
1211  grn_obj *top_left_point,
1212  grn_obj *bottom_right_point,
1213  int offset,
1214  int limit)
1215 {
1216  grn_geo_cursor_in_rectangle *cursor = NULL;
1217  in_rectangle_data data;
1218 
1219  GRN_API_ENTER;
1220  GRN_VOID_INIT(&(data.top_left_point_buffer));
1221  GRN_VOID_INIT(&(data.bottom_right_point_buffer));
1222  if (in_rectangle_data_prepare(ctx, index, top_left_point, bottom_right_point,
1223  "geo_in_rectangle()", &data)) {
1224  goto exit;
1225  }
1226 
1227  cursor = GRN_MALLOCN(grn_geo_cursor_in_rectangle, 1);
1228  if (!cursor) {
1229  ERR(GRN_NO_MEMORY_AVAILABLE,
1230  "[geo][cursor][in-rectangle] failed to allocate memory for geo cursor");
1231  goto exit;
1232  }
1233 
1234  cursor->pat = data.pat;
1235  cursor->index = index;
1236  memcpy(&(cursor->top_left), data.top_left, sizeof(grn_geo_point));
1237  memcpy(&(cursor->bottom_right), data.bottom_right, sizeof(grn_geo_point));
1238  grn_gton(cursor->top_left_key, data.top_left, sizeof(grn_geo_point));
1239  grn_gton(cursor->bottom_right_key, data.bottom_right, sizeof(grn_geo_point));
1240  cursor->pat_cursor = NULL;
1241  cursor->ii_cursor = NULL;
1242  cursor->offset = offset;
1243  cursor->rest = limit;
1244 
1245  cursor->current_entry = 0;
1246  {
1247  grn_geo_cursor_entry *entry;
1248 
1249  entry = &(cursor->entries[cursor->current_entry]);
1250  entry->target_bit = data.rectangle_common_bit;
1251  memcpy(entry->key, data.rectangle_common_key, sizeof(grn_geo_point));
1252  entry->status_flags =
1253  GRN_GEO_CURSOR_ENTRY_STATUS_TOP_INCLUDED |
1254  GRN_GEO_CURSOR_ENTRY_STATUS_BOTTOM_INCLUDED |
1255  GRN_GEO_CURSOR_ENTRY_STATUS_LEFT_INCLUDED |
1256  GRN_GEO_CURSOR_ENTRY_STATUS_RIGHT_INCLUDED;
1257  if (data.min.latitude == data.bottom_right->latitude &&
1258  data.max.latitude == data.top_left->latitude) {
1259  entry->status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER;
1260  }
1261  if (data.min.longitude == data.top_left->longitude &&
1262  data.max.longitude == data.bottom_right->longitude) {
1263  entry->status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER;
1264  }
1265  }
1266  {
1267  const char *minimum_reduce_bit_env;
1268  cursor->minimum_reduce_bit = 0;
1269  minimum_reduce_bit_env = getenv("GRN_GEO_IN_RECTANGLE_MINIMUM_REDUCE_BIT");
1270  if (minimum_reduce_bit_env) {
1271  cursor->minimum_reduce_bit = atoi(minimum_reduce_bit_env);
1272  }
1273  if (cursor->minimum_reduce_bit < 1) {
1274  cursor->minimum_reduce_bit = 1;
1275  }
1276  }
1277  GRN_DB_OBJ_SET_TYPE(cursor, GRN_CURSOR_COLUMN_GEO_INDEX);
1278  {
1279  grn_obj *db;
1280  grn_id id;
1281  db = grn_ctx_db(ctx);
1282  id = grn_obj_register(ctx, db, NULL, 0);
1283  DB_OBJ(cursor)->header.domain = GRN_ID_NIL;
1284  DB_OBJ(cursor)->range = GRN_ID_NIL;
1285  grn_db_obj_init(ctx, db, id, DB_OBJ(cursor));
1286  }
1287 
1288 exit :
1289  grn_obj_unlink(ctx, &(data.top_left_point_buffer));
1290  grn_obj_unlink(ctx, &(data.bottom_right_point_buffer));
1291  GRN_API_RETURN((grn_obj *)cursor);
1292 }
1293 
1294 static inline grn_bool
1295 grn_geo_cursor_entry_next_push(grn_ctx *ctx,
1296  grn_geo_cursor_in_rectangle *cursor,
1297  grn_geo_cursor_entry *entry)
1298 {
1299  grn_geo_cursor_entry *next_entry;
1300  grn_geo_point entry_base;
1301  grn_table_cursor *pat_cursor;
1302  grn_bool pushed = GRN_FALSE;
1303 
1304  grn_ntog((uint8_t*)(&entry_base), entry->key, sizeof(grn_geo_point));
1305  pat_cursor = grn_table_cursor_open(ctx,
1306  cursor->pat,
1307  &entry_base,
1308  entry->target_bit + 1,
1309  NULL, 0,
1310  0, -1,
1311  GRN_CURSOR_PREFIX|GRN_CURSOR_SIZE_BY_BIT);
1312  if (pat_cursor) {
1313  if (grn_table_cursor_next(ctx, pat_cursor)) {
1314  next_entry = &(cursor->entries[++cursor->current_entry]);
1315  memcpy(next_entry, entry, sizeof(grn_geo_cursor_entry));
1316  pushed = GRN_TRUE;
1317  }
1318  grn_table_cursor_close(ctx, pat_cursor);
1319  }
1320 
1321  return pushed;
1322 }
1323 
1324 static inline grn_bool
1325 grn_geo_cursor_entry_next(grn_ctx *ctx,
1326  grn_geo_cursor_in_rectangle *cursor,
1327  grn_geo_cursor_entry *entry)
1328 {
1329  uint8_t *top_left_key = cursor->top_left_key;
1330  uint8_t *bottom_right_key = cursor->bottom_right_key;
1331  int max_target_bit = GRN_GEO_KEY_MAX_BITS - cursor->minimum_reduce_bit;
1332 
1333  if (cursor->current_entry < 0) {
1334  return GRN_FALSE;
1335  }
1336 
1337  memcpy(entry,
1338  &(cursor->entries[cursor->current_entry--]),
1339  sizeof(grn_geo_cursor_entry));
1340  while (GRN_TRUE) {
1341  grn_geo_cursor_entry next_entry0, next_entry1;
1342  grn_bool pushed = GRN_FALSE;
1343 
1344  /*
1345  top_left_key: tl
1346  bottom_right_key: br
1347 
1348  e.g.: top_left_key is at the top left sub mesh and
1349  bottom_right_key is at the bottom right sub mesh.
1350  top_left_key is also at the top left - bottom right
1351  sub-sub mesh and
1352  bottom_right_key is at the bottom right - bottom left
1353  sub-sub mesh.
1354 
1355  ^latitude +----+----+----+----+
1356  | 1 |1010|1011|1110|1111|
1357  | | | | | |
1358  | 1 +----+----+----+----+
1359  \/ 0 |1000|1001|1100|1101|
1360  | | tl | | |
1361  +----+----+----+----+
1362  1 |0010|0011|0110|0111|
1363  | | | | |
1364  0 +----+----+----+----+
1365  0 |0000|0001|0100|0101|
1366  | | | br | |
1367  +----+----+----+----+
1368  0 1 0 1
1369  |-------| |-------|
1370  0 1
1371  <------>
1372  longitude
1373 
1374  entry.target_bit + 1 -> next_entry0
1375  entry.target_bit + 1 and entry.key ^ (entry.target_bit + 1) in bit
1376  -> next_entry1
1377 
1378  entry: represents the biggest mesh.
1379  (1010, 1011, 1110, 1111,
1380  1000, 1001, 1100, 1101,
1381  0010, 0011, 0110, 0111,
1382  0000, 0001, 0100, 0101)
1383  next_entry0: represents bottom sub-mesh.
1384  (0010, 0011, 0110, 0111,
1385  0000, 0001, 0100, 0101)
1386  next_entry1: represents top sub-mesh.
1387  (1010, 1011, 1110, 1111,
1388  1000, 1001, 1100, 1101)
1389 
1390  entry->status_flags = TOP_INCLUDED |
1391  BOTTOM_INCLUDED |
1392  LEFT_INCLUDED |
1393  RIGHT_INCLUDED
1394  next_entry0->status_flags = BOTTOM_INCLUDED |
1395  LEFT_INCLUDED |
1396  RIGHT_INCLUDED
1397  next_entry1->status_flags = TOP_INCLUDED |
1398  LEFT_INCLUDED |
1399  RIGHT_INCLUDED
1400 
1401  Both next_entry1 and next_entry0 are pushed to the stack in cursor.
1402  */
1403 
1404 #ifdef GEO_DEBUG
1405  inspect_cursor_entry(ctx, entry);
1406 #endif
1407 
1408  if (entry->target_bit >= max_target_bit) {
1409 #ifdef GEO_DEBUG
1410  printf("%d: force stopping to reduce a mesh\n", entry->target_bit);
1411 #endif
1412  break;
1413  }
1414 
1415  if (CURSOR_ENTRY_IS_INNER(entry)) {
1416 #ifdef GEO_DEBUG
1417  printf("%d: inner entries\n", entry->target_bit);
1418 #endif
1419  break;
1420  }
1421 
1422  memcpy(&next_entry0, entry, sizeof(grn_geo_cursor_entry));
1423  next_entry0.target_bit++;
1424  memcpy(&next_entry1, entry, sizeof(grn_geo_cursor_entry));
1425  next_entry1.target_bit++;
1426  SET_N_BIT(next_entry1.key, next_entry1.target_bit);
1427 
1428 #ifdef GEO_DEBUG
1429  inspect_cursor_entry_targets(ctx, entry, top_left_key, bottom_right_key,
1430  &next_entry0, &next_entry1);
1431 #endif
1432 
1433  if ((entry->target_bit + 1) % 2 == 0) {
1434  if (CURSOR_ENTRY_CHECK_STATUS(entry, TOP_INCLUDED)) {
1435  CURSOR_ENTRY_UPDATE_STATUS(&next_entry0, TOP_INCLUDED, top_left_key);
1436  CURSOR_ENTRY_UPDATE_STATUS(&next_entry1, TOP_INCLUDED, top_left_key);
1437  }
1438  if (CURSOR_ENTRY_CHECK_STATUS(entry, BOTTOM_INCLUDED)) {
1439  CURSOR_ENTRY_UPDATE_STATUS(&next_entry0, BOTTOM_INCLUDED,
1440  bottom_right_key);
1441  CURSOR_ENTRY_UPDATE_STATUS(&next_entry1, BOTTOM_INCLUDED,
1442  bottom_right_key);
1443  }
1444 
1445  if (CURSOR_ENTRY_CHECK_STATUS(entry, TOP_INCLUDED) &&
1446  !CURSOR_ENTRY_CHECK_STATUS(entry, BOTTOM_INCLUDED) &&
1447  CURSOR_ENTRY_CHECK_STATUS(&next_entry1, TOP_INCLUDED)) {
1448  next_entry0.status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER;
1449  } else if (!CURSOR_ENTRY_CHECK_STATUS(entry, TOP_INCLUDED) &&
1450  CURSOR_ENTRY_CHECK_STATUS(entry, BOTTOM_INCLUDED) &&
1451  CURSOR_ENTRY_CHECK_STATUS(&next_entry0, BOTTOM_INCLUDED)) {
1452  next_entry1.status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LATITUDE_INNER;
1453  }
1454 
1455  if (CURSOR_ENTRY_INCLUDED_IN_LATITUDE_DIRECTION(&next_entry1)) {
1456  if (grn_geo_cursor_entry_next_push(ctx, cursor, &next_entry1)) {
1457  pushed = GRN_TRUE;
1458 #ifdef GEO_DEBUG
1459  printf("%d: latitude: push 1\n", next_entry1.target_bit);
1460 #endif
1461  }
1462  }
1463  if (CURSOR_ENTRY_INCLUDED_IN_LATITUDE_DIRECTION(&next_entry0)) {
1464  if (grn_geo_cursor_entry_next_push(ctx, cursor, &next_entry0)) {
1465  pushed = GRN_TRUE;
1466 #ifdef GEO_DEBUG
1467  printf("%d: latitude: push 0\n", next_entry0.target_bit);
1468 #endif
1469  }
1470  }
1471  } else {
1472  if (CURSOR_ENTRY_CHECK_STATUS(entry, RIGHT_INCLUDED)) {
1473  CURSOR_ENTRY_UPDATE_STATUS(&next_entry0, RIGHT_INCLUDED,
1474  bottom_right_key);
1475  CURSOR_ENTRY_UPDATE_STATUS(&next_entry1, RIGHT_INCLUDED,
1476  bottom_right_key);
1477  }
1478  if (CURSOR_ENTRY_CHECK_STATUS(entry, LEFT_INCLUDED)) {
1479  CURSOR_ENTRY_UPDATE_STATUS(&next_entry0, LEFT_INCLUDED, top_left_key);
1480  CURSOR_ENTRY_UPDATE_STATUS(&next_entry1, LEFT_INCLUDED, top_left_key);
1481  }
1482 
1483  if (CURSOR_ENTRY_CHECK_STATUS(entry, LEFT_INCLUDED) &&
1484  !CURSOR_ENTRY_CHECK_STATUS(entry, RIGHT_INCLUDED) &&
1485  CURSOR_ENTRY_CHECK_STATUS(&next_entry0, LEFT_INCLUDED)) {
1486  next_entry1.status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER;
1487  } else if (!CURSOR_ENTRY_CHECK_STATUS(entry, LEFT_INCLUDED) &&
1488  CURSOR_ENTRY_CHECK_STATUS(entry, RIGHT_INCLUDED) &&
1489  CURSOR_ENTRY_CHECK_STATUS(&next_entry1, RIGHT_INCLUDED)) {
1490  next_entry0.status_flags |= GRN_GEO_CURSOR_ENTRY_STATUS_LONGITUDE_INNER;
1491  }
1492 
1493  if (CURSOR_ENTRY_INCLUDED_IN_LONGITUDE_DIRECTION(&next_entry1)) {
1494  if (grn_geo_cursor_entry_next_push(ctx, cursor, &next_entry1)) {
1495  pushed = GRN_TRUE;
1496 #ifdef GEO_DEBUG
1497  printf("%d: longitude: push 1\n", next_entry1.target_bit);
1498 #endif
1499  }
1500  }
1501  if (CURSOR_ENTRY_INCLUDED_IN_LONGITUDE_DIRECTION(&next_entry0)) {
1502  if (grn_geo_cursor_entry_next_push(ctx, cursor, &next_entry0)) {
1503  pushed = GRN_TRUE;
1504 #ifdef GEO_DEBUG
1505  printf("%d: longitude: push 0\n", next_entry0.target_bit);
1506 #endif
1507  }
1508  }
1509  }
1510 
1511  if (pushed) {
1512 #ifdef GEO_DEBUG
1513  int i;
1514 
1515  printf("%d: pushed\n", entry->target_bit);
1516  printf("stack:\n");
1517  for (i = cursor->current_entry; i >= 0; i--) {
1518  grn_geo_cursor_entry *stack_entry;
1519  stack_entry = &(cursor->entries[i]);
1520  printf("%2d: ", i);
1521  inspect_key(ctx, stack_entry->key);
1522  printf(" ");
1523  print_key_mark(ctx, stack_entry->target_bit);
1524  }
1525 #endif
1526  memcpy(entry,
1527  &(cursor->entries[cursor->current_entry--]),
1528  sizeof(grn_geo_cursor_entry));
1529 #ifdef GEO_DEBUG
1530  printf("%d: pop entry\n", entry->target_bit);
1531 #endif
1532  } else {
1533  break;
1534  }
1535  }
1536 
1537 #ifdef GEO_DEBUG
1538  printf("found:\n");
1539  inspect_cursor_entry(ctx, entry);
1540 #endif
1541 
1542  return GRN_TRUE;
1543 }
1544 
1545 typedef grn_bool (*grn_geo_cursor_callback)(grn_ctx *ctx, grn_ii_posting *posting, void *user_data);
1546 
1547 static void
1548 grn_geo_cursor_each(grn_ctx *ctx, grn_obj *geo_cursor,
1549  grn_geo_cursor_callback callback, void *user_data)
1550 {
1551  grn_geo_cursor_in_rectangle *cursor;
1552  grn_obj *pat;
1553  grn_table_cursor *pat_cursor;
1554  grn_ii *ii;
1555  grn_ii_cursor *ii_cursor;
1556  grn_ii_posting *posting = NULL;
1557  grn_geo_point *current, *top_left, *bottom_right;
1558  grn_id index_id;
1559 
1560  cursor = (grn_geo_cursor_in_rectangle *)geo_cursor;
1561  if (cursor->rest == 0) {
1562  return;
1563  }
1564 
1565  pat = cursor->pat;
1566  pat_cursor = cursor->pat_cursor;
1567  ii = (grn_ii *)(cursor->index);
1568  ii_cursor = cursor->ii_cursor;
1569  current = &(cursor->current);
1570  top_left = &(cursor->top_left);
1571  bottom_right = &(cursor->bottom_right);
1572 
1573  while (GRN_TRUE) {
1574  if (!pat_cursor) {
1575  grn_geo_cursor_entry entry;
1576  grn_geo_point entry_base;
1577  if (!grn_geo_cursor_entry_next(ctx, cursor, &entry)) {
1578  cursor->rest = 0;
1579  return;
1580  }
1581  grn_ntog((uint8_t*)(&entry_base), entry.key, sizeof(grn_geo_point));
1582  if (!(cursor->pat_cursor = pat_cursor =
1583  grn_table_cursor_open(ctx,
1584  pat,
1585  &entry_base,
1586  entry.target_bit + 1,
1587  NULL, 0,
1588  0, -1,
1589  GRN_CURSOR_PREFIX|GRN_CURSOR_SIZE_BY_BIT))) {
1590  cursor->rest = 0;
1591  return;
1592  }
1593 #ifdef GEO_DEBUG
1594  inspect_mesh(ctx, &entry_base, entry.target_bit, 0);
1595 #endif
1596  }
1597 
1598  while (ii_cursor || (index_id = grn_table_cursor_next(ctx, pat_cursor))) {
1599  if (!ii_cursor) {
1600  grn_table_get_key(ctx, pat, index_id, current, sizeof(grn_geo_point));
1601  if (grn_geo_in_rectangle_raw(ctx, current, top_left, bottom_right)) {
1602  inspect_tid(ctx, index_id, current, 0);
1603  if (!(cursor->ii_cursor = ii_cursor =
1604  grn_ii_cursor_open(ctx,
1605  ii,
1606  index_id,
1607  GRN_ID_NIL,
1608  GRN_ID_MAX,
1609  ii->n_elements,
1610  0))) {
1611  continue;
1612  }
1613  } else {
1614  continue;
1615  }
1616  }
1617 
1618  while ((posting = grn_ii_cursor_next(ctx, ii_cursor))) {
1619  if (cursor->offset == 0) {
1620  grn_bool keep_each;
1621  keep_each = callback(ctx, posting, user_data);
1622  if (cursor->rest > 0) {
1623  if (--(cursor->rest) == 0) {
1624  keep_each = GRN_FALSE;
1625  }
1626  }
1627  if (!keep_each) {
1628  return;
1629  }
1630  } else {
1631  cursor->offset--;
1632  }
1633  }
1634  grn_ii_cursor_close(ctx, ii_cursor);
1635  cursor->ii_cursor = ii_cursor = NULL;
1636  }
1637  grn_table_cursor_close(ctx, pat_cursor);
1638  cursor->pat_cursor = pat_cursor = NULL;
1639  }
1640 }
1641 
1642 static grn_bool
1643 grn_geo_cursor_next_callback(grn_ctx *ctx, grn_ii_posting *posting,
1644  void *user_data)
1645 {
1646  grn_ii_posting **return_posting = user_data;
1647  *return_posting = posting;
1648  return GRN_FALSE;
1649 }
1650 
1651 grn_posting *
1652 grn_geo_cursor_next(grn_ctx *ctx, grn_obj *geo_cursor)
1653 {
1654  grn_ii_posting *posting = NULL;
1655  grn_geo_cursor_each(ctx, geo_cursor, grn_geo_cursor_next_callback, &posting);
1656  return (grn_posting *)posting;
1657 }
1658 
1659 grn_rc
1660 grn_geo_cursor_close(grn_ctx *ctx, grn_obj *geo_cursor)
1661 {
1662  grn_geo_cursor_in_rectangle *cursor;
1663 
1664  if (!geo_cursor) { return GRN_INVALID_ARGUMENT; }
1665 
1666  cursor = (grn_geo_cursor_in_rectangle *)geo_cursor;
1667  if (cursor->pat) { grn_obj_unlink(ctx, cursor->pat); }
1668  if (cursor->index) { grn_obj_unlink(ctx, cursor->index); }
1669  if (cursor->pat_cursor) { grn_table_cursor_close(ctx, cursor->pat_cursor); }
1670  if (cursor->ii_cursor) { grn_ii_cursor_close(ctx, cursor->ii_cursor); }
1671  GRN_FREE(cursor);
1672 
1673  return GRN_SUCCESS;
1674 }
1675 
1676 typedef struct {
1677  grn_hash *res;
1678  grn_operator op;
1679 } grn_geo_select_in_rectangle_data;
1680 
1681 static grn_bool
1682 grn_geo_select_in_rectangle_callback(grn_ctx *ctx, grn_ii_posting *posting,
1683  void *user_data)
1684 {
1685  grn_geo_select_in_rectangle_data *data = user_data;
1686  grn_ii_posting_add(ctx, posting, data->res, data->op);
1687  return GRN_TRUE;
1688 }
1689 
1690 grn_rc
1691 grn_geo_select_in_rectangle(grn_ctx *ctx, grn_obj *index,
1692  grn_obj *top_left_point,
1693  grn_obj *bottom_right_point,
1694  grn_obj *res, grn_operator op)
1695 {
1696  grn_obj *cursor;
1697 
1698  cursor = grn_geo_cursor_open_in_rectangle(ctx, index,
1699  top_left_point, bottom_right_point,
1700  0, -1);
1701  if (cursor) {
1702  grn_geo_select_in_rectangle_data data;
1703  data.res = (grn_hash *)res;
1704  data.op = op;
1705  grn_geo_cursor_each(ctx, cursor, grn_geo_select_in_rectangle_callback,
1706  &data);
1707  grn_obj_unlink(ctx, cursor);
1708  grn_ii_resolve_sel_and(ctx, (grn_hash *)res, op);
1709  }
1710 
1711  return ctx->rc;
1712 }
1713 
1714 static grn_rc
1715 geo_point_get(grn_ctx *ctx, grn_obj *pat, int flags, grn_geo_point *geo_point)
1716 {
1717  grn_rc rc = GRN_SUCCESS;
1718  grn_id id;
1719  grn_table_cursor *cursor = NULL;
1720 
1721  cursor = grn_table_cursor_open(ctx, pat,
1722  NULL, 0,
1723  NULL, 0,
1724  0, 1,
1725  GRN_CURSOR_BY_KEY | flags);
1726  if (!cursor) {
1727  rc = ctx->rc;
1728  goto exit;
1729  }
1730 
1731  id = grn_table_cursor_next(ctx, cursor);
1732  if (id == GRN_ID_NIL) {
1733  rc = GRN_END_OF_DATA;
1734  } else {
1735  void *key;
1736  int key_size;
1737  key_size = grn_table_cursor_get_key(ctx, cursor, &key);
1738  memcpy(geo_point, key, key_size);
1739  }
1740 
1741 exit:
1742  if (cursor) {
1743  grn_table_cursor_close(ctx, cursor);
1744  }
1745  return rc;
1746 }
1747 
1748 int
1749 grn_geo_estimate_in_rectangle(grn_ctx *ctx,
1750  grn_obj *index,
1751  grn_obj *top_left_point,
1752  grn_obj *bottom_right_point)
1753 {
1754  int n = 0;
1755  int total_records;
1756  grn_rc rc;
1757  in_rectangle_data data;
1758 
1759  GRN_VOID_INIT(&(data.top_left_point_buffer));
1760  GRN_VOID_INIT(&(data.bottom_right_point_buffer));
1761  if (in_rectangle_data_prepare(ctx, index, top_left_point, bottom_right_point,
1762  "grn_geo_estimate_in_rectangle()", &data)) {
1763  n = -1;
1764  goto exit;
1765  }
1766 
1767  total_records = grn_table_size(ctx, data.pat);
1768  if (total_records > 0) {
1769  grn_geo_point min, max;
1770  int select_latitude_distance, select_longitude_distance;
1771  int total_latitude_distance, total_longitude_distance;
1772  double select_ratio;
1773  double estimated_n_records;
1774 
1775  rc = geo_point_get(ctx, data.pat, GRN_CURSOR_ASCENDING, &min);
1776  if (!rc) {
1777  rc = geo_point_get(ctx, data.pat, GRN_CURSOR_DESCENDING, &max);
1778  }
1779  if (rc) {
1780  if (rc == GRN_END_OF_DATA) {
1781  n = total_records;
1782  rc = GRN_SUCCESS;
1783  } else {
1784  n = -1;
1785  }
1786  goto exit;
1787  }
1788 
1789  select_latitude_distance = abs(data.max.latitude - data.min.latitude);
1790  select_longitude_distance = abs(data.max.longitude - data.min.longitude);
1791  total_latitude_distance = abs(max.latitude - min.latitude);
1792  total_longitude_distance = abs(max.longitude - min.longitude);
1793 
1794  select_ratio = 1.0;
1795  if (select_latitude_distance < total_latitude_distance) {
1796  select_ratio *= ((double)select_latitude_distance /
1797  (double)total_latitude_distance);
1798  }
1799  if (select_longitude_distance < total_longitude_distance) {
1800  select_ratio *= ((double)select_longitude_distance /
1801  (double)total_longitude_distance);
1802  }
1803  estimated_n_records = ceil(total_records * select_ratio);
1804  n = (int)estimated_n_records;
1805  }
1806 
1807 exit :
1808  grn_obj_unlink(ctx, &(data.top_left_point_buffer));
1809  grn_obj_unlink(ctx, &(data.bottom_right_point_buffer));
1810  return n;
1811 }
1812 
1813 grn_bool
1814 grn_geo_in_circle(grn_ctx *ctx, grn_obj *point, grn_obj *center,
1815  grn_obj *radius_or_point,
1816  grn_geo_approximate_type approximate_type)
1817 {
1818  grn_bool r = GRN_FALSE;
1819  grn_obj center_, radius_or_point_;
1820  grn_id domain = point->header.domain;
1821  if (domain == GRN_DB_TOKYO_GEO_POINT || domain == GRN_DB_WGS84_GEO_POINT) {
1822  grn_geo_distance_raw_func distance_raw_func;
1823  double d;
1824  if (center->header.domain != domain) {
1825  GRN_OBJ_INIT(&center_, GRN_BULK, 0, domain);
1826  if (grn_obj_cast(ctx, center, &center_, GRN_FALSE)) { goto exit; }
1827  center = &center_;
1828  }
1829 
1830  distance_raw_func = grn_geo_resolve_distance_raw_func(ctx,
1831  approximate_type,
1832  domain);
1833  if (!distance_raw_func) {
1834  ERR(GRN_INVALID_ARGUMENT,
1835  "unknown approximate type: <%d>", approximate_type);
1836  goto exit;
1837  }
1838  d = distance_raw_func(ctx,
1839  GRN_GEO_POINT_VALUE_RAW(point),
1840  GRN_GEO_POINT_VALUE_RAW(center));
1841  switch (radius_or_point->header.domain) {
1842  case GRN_DB_INT32 :
1843  r = d <= GRN_INT32_VALUE(radius_or_point);
1844  break;
1845  case GRN_DB_UINT32 :
1846  r = d <= GRN_UINT32_VALUE(radius_or_point);
1847  break;
1848  case GRN_DB_INT64 :
1849  r = d <= GRN_INT64_VALUE(radius_or_point);
1850  break;
1851  case GRN_DB_UINT64 :
1852  r = d <= GRN_UINT64_VALUE(radius_or_point);
1853  break;
1854  case GRN_DB_FLOAT :
1855  r = d <= GRN_FLOAT_VALUE(radius_or_point);
1856  break;
1857  case GRN_DB_SHORT_TEXT :
1858  case GRN_DB_TEXT :
1859  case GRN_DB_LONG_TEXT :
1860  GRN_OBJ_INIT(&radius_or_point_, GRN_BULK, 0, domain);
1861  if (grn_obj_cast(ctx, radius_or_point, &radius_or_point_, GRN_FALSE)) { goto exit; }
1862  radius_or_point = &radius_or_point_;
1863  /* fallthru */
1864  case GRN_DB_TOKYO_GEO_POINT :
1865  case GRN_DB_WGS84_GEO_POINT :
1866  if (domain != radius_or_point->header.domain) { /* todo */ goto exit; }
1867  r = d <= distance_raw_func(ctx,
1868  GRN_GEO_POINT_VALUE_RAW(radius_or_point),
1869  GRN_GEO_POINT_VALUE_RAW(center));
1870  break;
1871  default :
1872  goto exit;
1873  }
1874  } else {
1875  /* todo */
1876  }
1877 exit :
1878  return r;
1879 }
1880 
1881 grn_bool
1882 grn_geo_in_rectangle_raw(grn_ctx *ctx, grn_geo_point *point,
1883  grn_geo_point *top_left, grn_geo_point *bottom_right)
1884 {
1885  return ((top_left->longitude <= point->longitude) &&
1886  (point->longitude <= bottom_right->longitude) &&
1887  (bottom_right->latitude <= point->latitude) &&
1888  (point->latitude <= top_left->latitude));
1889 }
1890 
1891 grn_bool
1892 grn_geo_in_rectangle(grn_ctx *ctx, grn_obj *point,
1893  grn_obj *top_left, grn_obj *bottom_right)
1894 {
1895  grn_bool r = GRN_FALSE;
1896  grn_obj top_left_, bottom_right_;
1897  grn_id domain = point->header.domain;
1898  if (domain == GRN_DB_TOKYO_GEO_POINT || domain == GRN_DB_WGS84_GEO_POINT) {
1899  if (top_left->header.domain != domain) {
1900  GRN_OBJ_INIT(&top_left_, GRN_BULK, 0, domain);
1901  if (grn_obj_cast(ctx, top_left, &top_left_, GRN_FALSE)) { goto exit; }
1902  top_left = &top_left_;
1903  }
1904  if (bottom_right->header.domain != domain) {
1905  GRN_OBJ_INIT(&bottom_right_, GRN_BULK, 0, domain);
1906  if (grn_obj_cast(ctx, bottom_right, &bottom_right_, GRN_FALSE)) { goto exit; }
1907  bottom_right = &bottom_right_;
1908  }
1909  r = grn_geo_in_rectangle_raw(ctx,
1910  GRN_GEO_POINT_VALUE_RAW(point),
1911  GRN_GEO_POINT_VALUE_RAW(top_left),
1912  GRN_GEO_POINT_VALUE_RAW(bottom_right));
1913  } else {
1914  /* todo */
1915  }
1916 exit :
1917  return r;
1918 }
1919 
1920 typedef enum {
1921  LONGITUDE_SHORT,
1922  LONGITUDE_LONG,
1923 } distance_type;
1924 
1925 typedef enum {
1926  QUADRANT_1ST,
1927  QUADRANT_2ND,
1928  QUADRANT_3RD,
1929  QUADRANT_4TH,
1930  QUADRANT_1ST_TO_2ND,
1931  QUADRANT_1ST_TO_3RD,
1932  QUADRANT_1ST_TO_4TH,
1933  QUADRANT_2ND_TO_1ST,
1934  QUADRANT_2ND_TO_3RD,
1935  QUADRANT_2ND_TO_4TH,
1936  QUADRANT_3RD_TO_1ST,
1937  QUADRANT_3RD_TO_2ND,
1938  QUADRANT_3RD_TO_4TH,
1939  QUADRANT_4TH_TO_1ST,
1940  QUADRANT_4TH_TO_2ND,
1941  QUADRANT_4TH_TO_3RD,
1942 } quadrant_type;
1943 
1944 static distance_type
1945 geo_longitude_distance_type(int start_longitude, int end_longitude)
1946 {
1947  int diff_longitude;
1948  int east_to_west;
1949  int west_to_east;
1950  if (start_longitude >= 0) {
1951  diff_longitude = abs(start_longitude - end_longitude);
1952  } else {
1953  diff_longitude = abs(end_longitude - start_longitude);
1954  }
1955  east_to_west = start_longitude > 0 && end_longitude < 0;
1956  west_to_east = start_longitude < 0 && end_longitude > 0;
1957  if (start_longitude != end_longitude &&
1958  (east_to_west || west_to_east) &&
1959  diff_longitude > 180 * GRN_GEO_RESOLUTION) {
1960  return LONGITUDE_LONG;
1961  } else {
1962  return LONGITUDE_SHORT;
1963  }
1964 }
1965 
1966 static inline quadrant_type
1967 geo_quadrant_type(grn_geo_point *point1, grn_geo_point *point2)
1968 {
1969 #define QUADRANT_1ST_WITH_AXIS(point) \
1970  (point->longitude >= 0) && (point->latitude >= 0)
1971 #define QUADRANT_2ND_WITH_AXIS(point) \
1972  (point->longitude <= 0) && (point->latitude >= 0)
1973 #define QUADRANT_3RD_WITH_AXIS(point) \
1974  (point->longitude <= 0) && (point->latitude <= 0)
1975 #define QUADRANT_4TH_WITH_AXIS(point) \
1976  (point->longitude >= 0) && (point->latitude <= 0)
1977 
1978  if (QUADRANT_1ST_WITH_AXIS(point1) && QUADRANT_1ST_WITH_AXIS(point2)) {
1979  return QUADRANT_1ST;
1980  } else if (QUADRANT_2ND_WITH_AXIS(point1) && QUADRANT_2ND_WITH_AXIS(point2)) {
1981  return QUADRANT_2ND;
1982  } else if (QUADRANT_3RD_WITH_AXIS(point1) && QUADRANT_3RD_WITH_AXIS(point2)) {
1983  return QUADRANT_3RD;
1984  } else if (QUADRANT_4TH_WITH_AXIS(point1) && QUADRANT_4TH_WITH_AXIS(point2)) {
1985  return QUADRANT_4TH;
1986  } else {
1987  if (point1->longitude > 0 && point2->longitude < 0 &&
1988  point1->latitude >= 0 && point2->latitude >= 0) {
1989  return QUADRANT_1ST_TO_2ND;
1990  } else if (point1->longitude < 0 && point2->longitude > 0 &&
1991  point1->latitude >= 0 && point2->latitude >= 0) {
1992  return QUADRANT_2ND_TO_1ST;
1993  } else if (point1->longitude < 0 && point2->longitude > 0 &&
1994  point1->latitude <= 0 && point2->latitude <= 0) {
1995  return QUADRANT_3RD_TO_4TH;
1996  } else if (point1->longitude > 0 && point2->longitude < 0 &&
1997  point1->latitude <= 0 && point2->latitude <= 0) {
1998  return QUADRANT_4TH_TO_3RD;
1999  } else if (point1->longitude >= 0 && point2->longitude >= 0 &&
2000  point1->latitude > 0 && point2->latitude < 0) {
2001  return QUADRANT_1ST_TO_4TH;
2002  } else if (point1->longitude >= 0 && point2->longitude >= 0 &&
2003  point1->latitude < 0 && point2->latitude > 0) {
2004  return QUADRANT_4TH_TO_1ST;
2005  } else if (point1->longitude <= 0 && point2->longitude <= 0 &&
2006  point1->latitude > 0 && point2->latitude < 0) {
2007  return QUADRANT_2ND_TO_3RD;
2008  } else if (point1->longitude <= 0 && point2->longitude <= 0 &&
2009  point1->latitude < 0 && point2->latitude > 0) {
2010  return QUADRANT_3RD_TO_2ND;
2011  } else if (point1->longitude >= 0 && point2->longitude <= 0 &&
2012  point1->latitude > 0 && point2->latitude < 0) {
2013  return QUADRANT_1ST_TO_3RD;
2014  } else if (point1->longitude <= 0 && point2->longitude >= 0 &&
2015  point1->latitude < 0 && point2->latitude > 0) {
2016  return QUADRANT_3RD_TO_1ST;
2017  } else if (point1->longitude <= 0 && point2->longitude >= 0 &&
2018  point1->latitude > 0 && point2->latitude < 0) {
2019  return QUADRANT_2ND_TO_4TH;
2020  } else if (point1->longitude >= 0 && point2->longitude <= 0 &&
2021  point1->latitude < 0 && point2->latitude > 0) {
2022  return QUADRANT_4TH_TO_2ND;
2023  } else {
2024  /* FIXME */
2025  return QUADRANT_1ST;
2026  }
2027  }
2028 #undef QUADRANT_1ST_WITH_AXIS
2029 #undef QUADRANT_2ND_WITH_AXIS
2030 #undef QUADRANT_3RD_WITH_AXIS
2031 #undef QUADRANT_4TH_WITH_AXIS
2032 }
2033 
2034 static inline double
2035 geo_distance_rectangle_square_root(double start_longitude, double start_latitude,
2036  double end_longitude, double end_latitude)
2037 {
2038  double diff_longitude;
2039  double x, y;
2040 
2041  diff_longitude = end_longitude - start_longitude;
2042  x = diff_longitude * cos((start_latitude + end_latitude) * 0.5);
2043  y = end_latitude - start_latitude;
2044  return sqrt((x * x) + (y * y));
2045 }
2046 
2047 static inline double
2048 geo_distance_rectangle_short_dist_type(quadrant_type quad_type,
2049  double lng1, double lat1,
2050  double lng2, double lat2)
2051 {
2052  double distance;
2053  double longitude_delta, latitude_delta;
2054 
2055  if (quad_type == QUADRANT_1ST_TO_4TH ||
2056  quad_type == QUADRANT_4TH_TO_1ST ||
2057  quad_type == QUADRANT_2ND_TO_3RD ||
2058  quad_type == QUADRANT_3RD_TO_2ND) {
2059  longitude_delta = lng2 - lng1;
2060  if (longitude_delta > 0 || longitude_delta < 0) {
2061  if (lat2 > lat1) {
2062  distance = geo_distance_rectangle_square_root(lng1,
2063  lat1,
2064  lng2,
2065  lat2) * GRN_GEO_RADIUS;
2066  } else {
2067  distance = geo_distance_rectangle_square_root(lng2,
2068  lat2,
2069  lng1,
2070  lat1) * GRN_GEO_RADIUS;
2071  }
2072  } else {
2073  latitude_delta = fabs(lat1) + fabs(lat2);
2074  distance = sqrt(latitude_delta * latitude_delta) * GRN_GEO_RADIUS;
2075  }
2076  } else if (quad_type == QUADRANT_1ST_TO_3RD ||
2077  quad_type == QUADRANT_2ND_TO_4TH) {
2078  distance = geo_distance_rectangle_square_root(lng1,
2079  lat1,
2080  lng2,
2081  lat2) * GRN_GEO_RADIUS;
2082  } else if (quad_type == QUADRANT_3RD_TO_1ST ||
2083  quad_type == QUADRANT_4TH_TO_2ND) {
2084  distance = geo_distance_rectangle_square_root(lng2,
2085  lat2,
2086  lng1,
2087  lat1) * GRN_GEO_RADIUS;
2088  } else if (quad_type == QUADRANT_1ST_TO_2ND ||
2089  quad_type == QUADRANT_2ND_TO_1ST ||
2090  quad_type == QUADRANT_3RD_TO_4TH ||
2091  quad_type == QUADRANT_4TH_TO_3RD) {
2092  if (lat2 > lat1) {
2093  distance = geo_distance_rectangle_square_root(lng1,
2094  lat1,
2095  lng2,
2096  lat2) * GRN_GEO_RADIUS;
2097  } else if (lat2 < lat1) {
2098  distance = geo_distance_rectangle_square_root(lng2,
2099  lat2,
2100  lng1,
2101  lat1) * GRN_GEO_RADIUS;
2102  } else {
2103  longitude_delta = lng2 - lng1;
2104  distance = longitude_delta * cos(lat1);
2105  distance = sqrt(distance * distance) * GRN_GEO_RADIUS;
2106  }
2107  } else {
2108  distance = geo_distance_rectangle_square_root(lng1,
2109  lat1,
2110  lng2,
2111  lat2) * GRN_GEO_RADIUS;
2112  }
2113  return distance;
2114 }
2115 
2116 static inline double
2117 geo_distance_rectangle_long_dist_type(quadrant_type quad_type,
2118  double lng1, double lat1,
2119  double lng2, double lat2)
2120 {
2121 #define M_2PI 6.28318530717958647692
2122 
2123  double distance;
2124 
2125  if (quad_type == QUADRANT_1ST_TO_2ND ||
2126  quad_type == QUADRANT_4TH_TO_3RD) {
2127  if (lat1 > lat2) {
2128  distance = geo_distance_rectangle_square_root(lng2 + M_2PI,
2129  lat2,
2130  lng1,
2131  lat1) * GRN_GEO_RADIUS;
2132  } else {
2133  distance = geo_distance_rectangle_square_root(lng1,
2134  lat1,
2135  lng2 + M_2PI,
2136  lat2) * GRN_GEO_RADIUS;
2137  }
2138  } else if (quad_type == QUADRANT_2ND_TO_1ST ||
2139  quad_type == QUADRANT_3RD_TO_4TH) {
2140  if (lat1 > lat2) {
2141  distance = geo_distance_rectangle_square_root(lng2,
2142  lat2,
2143  lng1 + M_2PI,
2144  lat1) * GRN_GEO_RADIUS;
2145  } else {
2146  distance = geo_distance_rectangle_square_root(lng1 + M_2PI,
2147  lat1,
2148  lng2,
2149  lat2) * GRN_GEO_RADIUS;
2150  }
2151  } else if (quad_type == QUADRANT_1ST_TO_3RD) {
2152  distance = geo_distance_rectangle_square_root(lng2 + M_2PI,
2153  lat2,
2154  lng1,
2155  lat1) * GRN_GEO_RADIUS;
2156  } else if (quad_type == QUADRANT_3RD_TO_1ST) {
2157  distance = geo_distance_rectangle_square_root(lng1 + M_2PI,
2158  lat1,
2159  lng2,
2160  lat2) * GRN_GEO_RADIUS;
2161  } else if (quad_type == QUADRANT_2ND_TO_4TH) {
2162  distance = geo_distance_rectangle_square_root(lng2,
2163  lat2,
2164  lng1 + M_2PI,
2165  lat1) * GRN_GEO_RADIUS;
2166  } else if (quad_type == QUADRANT_4TH_TO_2ND) {
2167  distance = geo_distance_rectangle_square_root(lng1,
2168  lat1,
2169  lng2 + M_2PI,
2170  lat2) * GRN_GEO_RADIUS;
2171  } else {
2172  if (lng1 > lng2) {
2173  distance = geo_distance_rectangle_square_root(lng1,
2174  lat1,
2175  lng2 + M_2PI,
2176  lat2) * GRN_GEO_RADIUS;
2177  } else {
2178  distance = geo_distance_rectangle_square_root(lng2,
2179  lat2,
2180  lng1 + M_2PI,
2181  lat1) * GRN_GEO_RADIUS;
2182  }
2183  }
2184  return distance;
2185 #undef M_2PI
2186 }
2187 
2188 double
2189 grn_geo_distance_rectangle_raw(grn_ctx *ctx,
2190  grn_geo_point *point1, grn_geo_point *point2)
2191 {
2192 
2193  double lng1, lat1, lng2, lat2, distance;
2194  distance_type dist_type;
2195  quadrant_type quad_type;
2196 
2197  lat1 = GRN_GEO_INT2RAD(point1->latitude);
2198  lng1 = GRN_GEO_INT2RAD(point1->longitude);
2199  lat2 = GRN_GEO_INT2RAD(point2->latitude);
2200  lng2 = GRN_GEO_INT2RAD(point2->longitude);
2201  quad_type = geo_quadrant_type(point1, point2);
2202  if (quad_type <= QUADRANT_4TH) {
2203  distance = geo_distance_rectangle_square_root(lng1,
2204  lat1,
2205  lng2,
2206  lat2) * GRN_GEO_RADIUS;
2207  } else {
2208  dist_type = geo_longitude_distance_type(point1->longitude,
2209  point2->longitude);
2210  if (dist_type == LONGITUDE_SHORT) {
2211  distance = geo_distance_rectangle_short_dist_type(quad_type,
2212  lng1, lat1,
2213  lng2, lat2);
2214  } else {
2215  distance = geo_distance_rectangle_long_dist_type(quad_type,
2216  lng1, lat1,
2217  lng2, lat2);
2218  }
2219  }
2220  return distance;
2221 }
2222 
2223 double
2224 grn_geo_distance_sphere_raw(grn_ctx *ctx,
2225  grn_geo_point *point1, grn_geo_point *point2)
2226 {
2227  double lng1, lat1, lng2, lat2, x, y;
2228 
2229  lat1 = GRN_GEO_INT2RAD(point1->latitude);
2230  lng1 = GRN_GEO_INT2RAD(point1->longitude);
2231  lat2 = GRN_GEO_INT2RAD(point2->latitude);
2232  lng2 = GRN_GEO_INT2RAD(point2->longitude);
2233  x = sin(fabs(lng2 - lng1) * 0.5);
2234  y = sin(fabs(lat2 - lat1) * 0.5);
2235  return asin(sqrt((y * y) + cos(lat1) * cos(lat2) * x * x)) * 2 * GRN_GEO_RADIUS;
2236 }
2237 
2238 double
2239 grn_geo_distance_ellipsoid_raw(grn_ctx *ctx,
2240  grn_geo_point *point1, grn_geo_point *point2,
2241  int c1, int c2, double c3)
2242 {
2243  double lng1, lat1, lng2, lat2, p, q, r, m, n, x, y;
2244 
2245  lat1 = GRN_GEO_INT2RAD(point1->latitude);
2246  lng1 = GRN_GEO_INT2RAD(point1->longitude);
2247  lat2 = GRN_GEO_INT2RAD(point2->latitude);
2248  lng2 = GRN_GEO_INT2RAD(point2->longitude);
2249  p = (lat1 + lat2) * 0.5;
2250  q = (1 - c3 * sin(p) * sin(p));
2251  r = sqrt(q);
2252  m = c1 / (q * r);
2253  n = c2 / r;
2254  x = n * cos(p) * fabs(lng1 - lng2);
2255  y = m * fabs(lat1 - lat2);
2256  return sqrt((x * x) + (y * y));
2257 }
2258 
2259 double
2260 grn_geo_distance_ellipsoid_raw_tokyo(grn_ctx *ctx,
2261  grn_geo_point *point1,
2262  grn_geo_point *point2)
2263 {
2264  return grn_geo_distance_ellipsoid_raw(ctx, point1, point2,
2265  GRN_GEO_BES_C1,
2266  GRN_GEO_BES_C2,
2267  GRN_GEO_BES_C3);
2268 }
2269 
2270 double
2271 grn_geo_distance_ellipsoid_raw_wgs84(grn_ctx *ctx,
2272  grn_geo_point *point1,
2273  grn_geo_point *point2)
2274 {
2275  return grn_geo_distance_ellipsoid_raw(ctx, point1, point2,
2276  GRN_GEO_GRS_C1,
2277  GRN_GEO_GRS_C2,
2278  GRN_GEO_GRS_C3);
2279 }
2280 
2281 double
2282 grn_geo_distance(grn_ctx *ctx, grn_obj *point1, grn_obj *point2,
2283  grn_geo_approximate_type type)
2284 {
2285  double d = 0.0;
2286 
2287  switch (type) {
2288  case GRN_GEO_APPROXIMATE_RECTANGLE :
2289  d = grn_geo_distance_rectangle(ctx, point1, point2);
2290  break;
2291  case GRN_GEO_APPROXIMATE_SPHERE :
2292  d = grn_geo_distance_sphere(ctx, point1, point2);
2293  break;
2294  case GRN_GEO_APPROXIMATE_ELLIPSOID :
2295  d = grn_geo_distance_ellipsoid(ctx, point1, point2);
2296  break;
2297  default :
2298  ERR(GRN_INVALID_ARGUMENT, "unknown approximate type: <%d>", type);
2299  break;
2300  }
2301  return d;
2302 }
2303 
2304 double
2305 grn_geo_distance_rectangle(grn_ctx *ctx, grn_obj *point1, grn_obj *point2)
2306 {
2307  double d = 0;
2308  grn_bool point1_initialized = GRN_FALSE;
2309  grn_bool point2_initialized = GRN_FALSE;
2310  grn_obj point1_, point2_;
2311  grn_id domain1 = point1->header.domain;
2312  grn_id domain2 = point2->header.domain;
2313  if (domain1 == GRN_DB_TOKYO_GEO_POINT || domain1 == GRN_DB_WGS84_GEO_POINT) {
2314  if (domain1 != domain2) {
2315  GRN_OBJ_INIT(&point2_, GRN_BULK, 0, domain1);
2316  point2_initialized = GRN_TRUE;
2317  if (grn_obj_cast(ctx, point2, &point2_, GRN_FALSE)) { goto exit; }
2318  point2 = &point2_;
2319  }
2320  } else if (domain2 == GRN_DB_TOKYO_GEO_POINT ||
2321  domain2 == GRN_DB_WGS84_GEO_POINT) {
2322  GRN_OBJ_INIT(&point1_, GRN_BULK, 0, domain2);
2323  point1_initialized = GRN_TRUE;
2324  if (grn_obj_cast(ctx, point1, &point1_, GRN_FALSE)) { goto exit; }
2325  point1 = &point1_;
2326  } else if ((GRN_DB_SHORT_TEXT <= domain1 && domain1 <= GRN_DB_LONG_TEXT) &&
2327  (GRN_DB_SHORT_TEXT <= domain2 && domain2 <= GRN_DB_LONG_TEXT)) {
2328  GRN_OBJ_INIT(&point1_, GRN_BULK, 0, GRN_DB_WGS84_GEO_POINT);
2329  point1_initialized = GRN_TRUE;
2330  if (grn_obj_cast(ctx, point1, &point1_, GRN_FALSE)) { goto exit; }
2331  point1 = &point1_;
2332 
2333  GRN_OBJ_INIT(&point2_, GRN_BULK, 0, GRN_DB_WGS84_GEO_POINT);
2334  point2_initialized = GRN_TRUE;
2335  if (grn_obj_cast(ctx, point2, &point2_, GRN_FALSE)) { goto exit; }
2336  point2 = &point2_;
2337  } else {
2338  goto exit;
2339  }
2340  d = grn_geo_distance_rectangle_raw(ctx,
2341  GRN_GEO_POINT_VALUE_RAW(point1),
2342  GRN_GEO_POINT_VALUE_RAW(point2));
2343 exit :
2344  if (point1_initialized) {
2345  GRN_OBJ_FIN(ctx, &point1_);
2346  }
2347  if (point2_initialized) {
2348  GRN_OBJ_FIN(ctx, &point2_);
2349  }
2350  return d;
2351 }
2352 
2353 double
2354 grn_geo_distance_sphere(grn_ctx *ctx, grn_obj *point1, grn_obj *point2)
2355 {
2356  double d = 0;
2357  grn_bool point2_initialized = GRN_FALSE;
2358  grn_obj point2_;
2359  grn_id domain = point1->header.domain;
2360  if (domain == GRN_DB_TOKYO_GEO_POINT || domain == GRN_DB_WGS84_GEO_POINT) {
2361  if (point2->header.domain != domain) {
2362  GRN_OBJ_INIT(&point2_, GRN_BULK, 0, domain);
2363  point2_initialized = GRN_TRUE;
2364  if (grn_obj_cast(ctx, point2, &point2_, GRN_FALSE)) { goto exit; }
2365  point2 = &point2_;
2366  }
2367  d = grn_geo_distance_sphere_raw(ctx,
2368  GRN_GEO_POINT_VALUE_RAW(point1),
2369  GRN_GEO_POINT_VALUE_RAW(point2));
2370  } else {
2371  /* todo */
2372  }
2373 exit :
2374  if (point2_initialized) {
2375  GRN_OBJ_FIN(ctx, &point2_);
2376  }
2377  return d;
2378 }
2379 
2380 double
2381 grn_geo_distance_ellipsoid(grn_ctx *ctx, grn_obj *point1, grn_obj *point2)
2382 {
2383  double d = 0;
2384  grn_bool point2_initialized = GRN_FALSE;
2385  grn_obj point2_;
2386  grn_id domain = point1->header.domain;
2387  if (domain == GRN_DB_TOKYO_GEO_POINT || domain == GRN_DB_WGS84_GEO_POINT) {
2388  if (point2->header.domain != domain) {
2389  GRN_OBJ_INIT(&point2_, GRN_BULK, 0, domain);
2390  point2_initialized = GRN_TRUE;
2391  if (grn_obj_cast(ctx, point2, &point2_, GRN_FALSE)) { goto exit; }
2392  point2 = &point2_;
2393  }
2394  if (domain == GRN_DB_TOKYO_GEO_POINT) {
2395  d = grn_geo_distance_ellipsoid_raw_tokyo(ctx,
2396  GRN_GEO_POINT_VALUE_RAW(point1),
2397  GRN_GEO_POINT_VALUE_RAW(point2));
2398  } else {
2399  d = grn_geo_distance_ellipsoid_raw_wgs84(ctx,
2400  GRN_GEO_POINT_VALUE_RAW(point1),
2401  GRN_GEO_POINT_VALUE_RAW(point2));
2402  }
2403  } else {
2404  /* todo */
2405  }
2406 exit :
2407  if (point2_initialized) {
2408  GRN_OBJ_FIN(ctx, &point2_);
2409  }
2410  return d;
2411 }
Generated on Sun Nov 10 2013 09:49:02 for Groonga 3.0.9 Source Code Document by   doxygen 1.8.1.2