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