source: trunk/liblwgeom/lwalgorithm.c @ 5274

Last change on this file since 5274 was 5274, checked in by colivier, 7 years ago

make astyle session

  • Property svn:keywords set to Author Date Id Revision
File size: 23.5 KB
Line 
1/**********************************************************************
2 * $Id: lwalgorithm.c 5274 2010-02-21 12:36:27Z colivier $
3 *
4 * PostGIS - Spatial Types for PostgreSQL
5 * http://postgis.refractions.net
6 * Copyright 2008 Paul Ramsey
7 *
8 * This is free software; you can redistribute and/or modify it under
9 * the terms of the GNU General Public Licence. See the COPYING file.
10 *
11 **********************************************************************/
12
13#include "lwalgorithm.h"
14
15
16
17/**
18** lw_segment_side()
19**
20** Return < 0.0 if point Q is left of segment P
21** Return > 0.0 if point Q is right of segment P
22** Return = 0.0 if point Q in on segment P
23*/
24double lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
25{
26        double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );
27        if ( FP_IS_ZERO(side) )
28                return 0.0;
29        else
30                return side;
31}
32
33
34int lw_segment_envelope_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
35{
36        double minq=LW_MIN(q1->x,q2->x);
37        double maxq=LW_MAX(q1->x,q2->x);
38        double minp=LW_MIN(p1->x,p2->x);
39        double maxp=LW_MAX(p1->x,p2->x);
40
41        if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
42                return LW_FALSE;
43
44        minq=LW_MIN(q1->y,q2->y);
45        maxq=LW_MAX(q1->y,q2->y);
46        minp=LW_MIN(p1->y,p2->y);
47        maxp=LW_MAX(p1->y,p2->y);
48
49        if (FP_GT(minp,maxq) || FP_LT(maxp,minq))
50                return LW_FALSE;
51
52        return LW_TRUE;
53}
54
55/**
56** @brief returns the kind of #CG_SEGMENT_INTERSECTION_TYPE  behavior of lineseg 1 (constructed from p1 and p2) and lineseg 2 (constructed from q1 and q2)
57**      @param p1 start point of first straight linesegment
58**      @param p2 end point of first straight linesegment
59**      @param q1 start point of second line segment
60**      @param q2 end point of second line segment
61**      @return a #CG_SEGMENT_INTERSECTION_TYPE
62**      Returns one of
63**              SEG_ERROR = -1,
64**              SEG_NO_INTERSECTION = 0,
65**              SEG_COLINEAR = 1,
66**              SEG_CROSS_LEFT = 2,
67**              SEG_CROSS_RIGHT = 3,
68*/
69int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)
70{
71
72        double pq1, pq2, qp1, qp2;
73
74        /* No envelope interaction => we are done. */
75        if (!lw_segment_envelope_intersects(p1, p2, q1, p2))
76        {
77                return SEG_NO_INTERSECTION;
78        }
79
80        /* Are the start and end points of q on the same side of p? */
81        pq1=lw_segment_side(p1,p2,q1);
82        pq2=lw_segment_side(p1,p2,q2);
83        if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0))
84        {
85                return SEG_NO_INTERSECTION;
86        }
87
88        /* Are the start and end points of p on the same side of q? */
89        qp1=lw_segment_side(q1,q2,p1);
90        qp2=lw_segment_side(q1,q2,p2);
91        if ( (qp1 > 0.0 && qp2 > 0.0) || (qp1 < 0.0 && qp2 < 0.0) )
92        {
93                return SEG_NO_INTERSECTION;
94        }
95
96        /* Nobody is on one side or another? Must be colinear. */
97        if ( pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0 )
98        {
99                return SEG_COLINEAR;
100        }
101
102        /*
103        ** When one end-point touches, the sidedness is determined by the
104        ** location of the other end-point. Only touches by the first point
105        ** will be considered "real" to avoid double counting.
106        */
107        LWDEBUGF(4, "pq1=%.15g pq2=%.15g", pq1, pq2);
108        LWDEBUGF(4, "qp1=%.15g qp2=%.15g", qp1, qp2);
109
110        /* Second point of p or q touches, it's not a crossing. */
111        if ( pq2 == 0.0 || qp2 == 0.0 )
112        {
113                return SEG_NO_INTERSECTION;
114        }
115
116        /* First point of p touches, it's a "crossing". */
117        if ( pq1 == 0.0 )
118        {
119                if ( FP_GT(pq2,0.0) )
120                        return SEG_CROSS_RIGHT;
121                else
122                        return SEG_CROSS_LEFT;
123        }
124
125        /* First point of q touches, it's a crossing. */
126        if ( qp1 == 0.0 )
127        {
128                if ( FP_LT(pq1,pq2) )
129                        return SEG_CROSS_RIGHT;
130                else
131                        return SEG_CROSS_LEFT;
132        }
133
134        /* The segments cross, what direction is the crossing? */
135        if ( FP_LT(pq1,pq2) )
136                return SEG_CROSS_RIGHT;
137        else
138                return SEG_CROSS_LEFT;
139
140        /* This should never happen! */
141        return SEG_ERROR;
142}
143
144/**
145** @brief lwline_crossing_direction: returns the kind of #CG_LINE_CROSS_TYPE behavior  of 2 linestrings
146** @param l1 first line string
147** @param l2 second line string
148** @return a #CG_LINE_CROSS_TYPE
149**   LINE_NO_CROSS = 0
150**   LINE_CROSS_LEFT = -1
151**   LINE_CROSS_RIGHT = 1
152**   LINE_MULTICROSS_END_LEFT = -2
153**   LINE_MULTICROSS_END_RIGHT = 2
154**   LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3
155**   LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3
156**
157*/
158int lwline_crossing_direction(LWLINE *l1, LWLINE *l2)
159{
160        int i = 0, j = 0, rv = 0;
161        POINT2D p1, p2, q1, q2;
162        POINTARRAY *pa1 = NULL, *pa2 = NULL;
163        int cross_left = 0;
164        int cross_right = 0;
165        int first_cross = 0;
166        int this_cross = 0;
167
168        pa1 = (POINTARRAY*)l1->points;
169        pa2 = (POINTARRAY*)l2->points;
170
171        /* One-point lines can't intersect (and shouldn't exist). */
172        if ( pa1->npoints < 2 || pa2->npoints < 2 )
173                return LINE_NO_CROSS;
174
175        LWDEBUGF(4, "l1 = %s", lwgeom_to_ewkt((LWGEOM*)l1,0));
176        LWDEBUGF(4, "l2 = %s", lwgeom_to_ewkt((LWGEOM*)l2,0));
177
178        /* Initialize first point of q */
179        rv = getPoint2d_p(pa2, 0, &q1);
180
181        for ( i = 1; i < pa2->npoints; i++ )
182        {
183
184                /* Update second point of q to next value */
185                rv = getPoint2d_p(pa2, i, &q2);
186
187                /* Initialize first point of p */
188                rv = getPoint2d_p(pa1, 0, &p1);
189
190                for ( j = 1; j < pa1->npoints; j++ )
191                {
192
193                        /* Update second point of p to next value */
194                        rv = getPoint2d_p(pa1, j, &p2);
195
196                        this_cross = lw_segment_intersects(&p1, &p2, &q1, &q2);
197
198                        LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", this_cross, i, j, p1.x, p1.y, p2.x, p2.y);
199
200                        if ( this_cross == SEG_CROSS_LEFT )
201                        {
202                                LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");
203                                cross_left++;
204                                if ( ! first_cross )
205                                        first_cross = SEG_CROSS_LEFT;
206                        }
207
208                        if ( this_cross == SEG_CROSS_RIGHT )
209                        {
210                                LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");
211                                cross_right++;
212                                if ( ! first_cross )
213                                        first_cross = SEG_CROSS_LEFT;
214                        }
215
216                        /*
217                        ** Crossing at a co-linearity can be turned handled by extending
218                        ** segment to next vertext and seeing if the end points straddle
219                        ** the co-linear segment.
220                        */
221                        if ( this_cross == SEG_COLINEAR )
222                        {
223                                LWDEBUG(4,"this_cross == SEG_COLINEAR");
224                                /* TODO: Add logic here and in segment_intersects()
225                                continue;
226                                */
227                        }
228
229                        LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");
230
231                        /* Turn second point of p into first point */
232                        p1 = p2;
233
234                }
235
236                /* Turn second point of q into first point */
237                q1 = q2;
238
239        }
240
241        LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);
242
243        if ( !cross_left && !cross_right )
244                return LINE_NO_CROSS;
245
246        if ( !cross_left && cross_right == 1 )
247                return LINE_CROSS_RIGHT;
248
249        if ( !cross_right && cross_left == 1 )
250                return LINE_CROSS_LEFT;
251
252        if ( cross_left - cross_right == 1 )
253                return LINE_MULTICROSS_END_LEFT;
254
255        if ( cross_left - cross_right == -1 )
256                return LINE_MULTICROSS_END_RIGHT;
257
258        if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )
259                return LINE_MULTICROSS_END_SAME_FIRST_LEFT;
260
261        if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )
262                return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;
263
264        return LINE_NO_CROSS;
265
266}
267
268
269
270/*
271** lwpoint_get_ordinate(point, ordinate) => double
272*/
273double lwpoint_get_ordinate(const POINT4D *p, int ordinate)
274{
275        if ( ! p )
276        {
277                lwerror("Null input geometry.");
278                return 0.0;
279        }
280
281        if ( ordinate > 3 || ordinate < 0 )
282        {
283                lwerror("Cannot extract ordinate %d.", ordinate);
284                return 0.0;
285        }
286
287        if ( ordinate == 3 )
288                return p->m;
289        if ( ordinate == 2 )
290                return p->z;
291        if ( ordinate == 1 )
292                return p->y;
293
294        return p->x;
295
296}
297void lwpoint_set_ordinate(POINT4D *p, int ordinate, double value)
298{
299        if ( ! p )
300        {
301                lwerror("Null input geometry.");
302                return;
303        }
304
305        if ( ordinate > 3 || ordinate < 0 )
306        {
307                lwerror("Cannot extract ordinate %d.", ordinate);
308                return;
309        }
310
311        LWDEBUGF(4, "    setting ordinate %d to %g", ordinate, value);
312
313        switch ( ordinate )
314        {
315        case 3:
316                p->m = value;
317                return;
318        case 2:
319                p->z = value;
320                return;
321        case 1:
322                p->y = value;
323                return;
324        case 0:
325                p->x = value;
326                return;
327        }
328}
329
330
331int lwpoint_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int ndims, int ordinate, double interpolation_value)
332{
333        double p1_value = lwpoint_get_ordinate(p1, ordinate);
334        double p2_value = lwpoint_get_ordinate(p2, ordinate);
335        double proportion;
336        int i = 0;
337
338        if ( ordinate < 0 || ordinate >= ndims )
339        {
340                lwerror("Ordinate (%d) is not within ndims (%d).", ordinate, ndims);
341                return 0;
342        }
343
344        if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
345                FP_MAX(p1_value, p2_value) < interpolation_value )
346        {
347                lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
348                return 0;
349        }
350
351        proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
352
353        for ( i = 0; i < ndims; i++ )
354        {
355                double newordinate = 0.0;
356                p1_value = lwpoint_get_ordinate(p1, i);
357                p2_value = lwpoint_get_ordinate(p2, i);
358                newordinate = p1_value + proportion * (p2_value - p1_value);
359                lwpoint_set_ordinate(p, i, newordinate);
360                LWDEBUGF(4, "   clip ordinate(%d) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", i, p1_value, p2_value, proportion, newordinate );
361        }
362
363        return 1;
364}
365
366LWCOLLECTION *lwmline_clip_to_ordinate_range(LWMLINE *mline, int ordinate, double from, double to)
367{
368        LWCOLLECTION *lwgeom_out = NULL;
369
370        if ( ! mline )
371        {
372                lwerror("Null input geometry.");
373                return NULL;
374        }
375
376        if ( mline->ngeoms == 1)
377        {
378                lwgeom_out = lwline_clip_to_ordinate_range(mline->geoms[0], ordinate, from, to);
379        }
380        else
381        {
382                LWCOLLECTION *col;
383                char hasz = TYPE_HASZ(mline->type);
384                char hasm = TYPE_HASM(mline->type);
385                char hassrid = TYPE_HASSRID(mline->type);
386                int i, j;
387                char homogeneous = 1;
388                size_t geoms_size = 0;
389                lwgeom_out = lwcollection_construct_empty(mline->SRID, hasz, hasm);
390                lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, MULTILINETYPE);
391                for ( i = 0; i < mline->ngeoms; i ++ )
392                {
393                        col = lwline_clip_to_ordinate_range(mline->geoms[i], ordinate, from, to);
394                        if ( col )
395                        {
396                                /* Something was left after the clip. */
397                                if ( lwgeom_out->ngeoms + col->ngeoms > geoms_size )
398                                {
399                                        geoms_size += 16;
400                                        if ( lwgeom_out->geoms )
401                                        {
402                                                lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, geoms_size * sizeof(LWGEOM*));
403                                        }
404                                        else
405                                        {
406                                                lwgeom_out->geoms = lwalloc(geoms_size * sizeof(LWGEOM*));
407                                        }
408                                }
409                                for ( j = 0; j < col->ngeoms; j++ )
410                                {
411                                        lwgeom_out->geoms[lwgeom_out->ngeoms] = col->geoms[j];
412                                        lwgeom_out->ngeoms++;
413                                }
414                                if ( TYPE_GETTYPE(col->type) != TYPE_GETTYPE(mline->type) )
415                                {
416                                        homogeneous = 0;
417                                }
418                                /* Shallow free the struct, leaving the geoms behind. */
419                                if ( col->bbox ) lwfree(col->bbox);
420                                lwfree(col);
421                        }
422                }
423                lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
424                lwgeom_add_bbox((LWGEOM*)lwgeom_out);
425                if ( ! homogeneous )
426                {
427                        lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
428                }
429        }
430
431        if ( ! lwgeom_out || lwgeom_out->ngeoms == 0 ) /* Nothing left after clip. */
432        {
433                return NULL;
434        }
435
436        return lwgeom_out;
437
438}
439
440
441/*
442** lwline_clip_to_ordinate_range(line, ordinate, from, to) => lwmline
443**
444** Take in a LINESTRING and return a MULTILINESTRING of those portions of the
445** LINESTRING between the from/to range for the specified ordinate (XYZM)
446*/
447LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double from, double to)
448{
449
450        POINTARRAY *pa_in = NULL;
451        LWCOLLECTION *lwgeom_out = NULL;
452        POINTARRAY *pa_out = NULL;
453        DYNPTARRAY *dp = NULL;
454        int i, rv;
455        int added_last_point = 0;
456        POINT4D *p = NULL, *q = NULL, *r = NULL;
457        double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
458        char hasz = TYPE_HASZ(line->type);
459        char hasm = TYPE_HASM(line->type);
460        char dims = TYPE_NDIMS(line->type);
461        char hassrid = TYPE_HASSRID(line->type);
462
463        LWDEBUGF(4, "hassrid = %d", hassrid);
464
465        /* Null input, nothing we can do. */
466        if ( ! line )
467        {
468                lwerror("Null input geometry.");
469                return NULL;
470        }
471
472        /* Ensure 'from' is less than 'to'. */
473        if ( to < from )
474        {
475                double t = from;
476                from = to;
477                to = t;
478        }
479
480        LWDEBUGF(4, "from = %g, to = %g, ordinate = %d", from, to, ordinate);
481        LWDEBUGF(4, "%s", lwgeom_to_ewkt((LWGEOM*)line, PARSER_CHECK_NONE));
482
483        /* Asking for an ordinate we don't have. Error. */
484        if ( ordinate >= dims )
485        {
486                lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
487                return NULL;
488        }
489
490        /* Prepare our working point objects. */
491        p = lwalloc(sizeof(POINT4D));
492        q = lwalloc(sizeof(POINT4D));
493        r = lwalloc(sizeof(POINT4D));
494
495        /* Construct a collection to hold our outputs. */
496        lwgeom_out = lwalloc(sizeof(LWCOLLECTION));
497        lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, MULTILINETYPE);
498        if (hassrid)
499                lwgeom_out->SRID = line->SRID;
500        else
501                lwgeom_out->SRID = -1;
502        lwgeom_out->bbox = NULL;
503        lwgeom_out->ngeoms = 0;
504        lwgeom_out->geoms = NULL;
505
506        pa_in = (POINTARRAY*)line->points;
507
508        for ( i = 0; i < pa_in->npoints; i++ )
509        {
510                LWDEBUGF(4, "Point #%d", i);
511                if ( i > 0 )
512                {
513                        q->x = p->x;
514                        q->y = p->y;
515                        q->z = p->z;
516                        q->m = p->m;
517                        ordinate_value_q = ordinate_value_p;
518                }
519                rv = getPoint4d_p(pa_in, i, p);
520                ordinate_value_p = lwpoint_get_ordinate(p, ordinate);
521                LWDEBUGF(4, " ordinate_value_p %g (current)", ordinate_value_p);
522                LWDEBUGF(4, " ordinate_value_q %g (previous)", ordinate_value_q);
523
524                /* Is this point inside the ordinate range? Yes. */
525                if ( ordinate_value_p >= from && ordinate_value_p <= to )
526                {
527                        LWDEBUGF(4, " inside ordinate range (%g, %g)", from, to);
528
529                        if ( ! added_last_point )
530                        {
531                                LWDEBUG(4,"  new ptarray required");
532                                /* We didn't add the previous point, so this is a new segment.
533                                *  Make a new point array. */
534                                if ( dp ) lwfree(dp);
535                                dp = dynptarray_create(64, line->type);
536
537                                /* We're transiting into the range so add an interpolated
538                                *  point at the range boundary.
539                                *  If we're on a boundary and crossing from the far side,
540                                *  we also need an interpolated point. */
541                                if ( i > 0 && ( /* Don't try to interpolate if this is the first point */
542                                            ( ordinate_value_p > from && ordinate_value_p < to ) || /* Inside */
543                                            ( ordinate_value_p == from && ordinate_value_q > to ) || /* Hopping from above */
544                                            ( ordinate_value_p == to && ordinate_value_q < from ) ) ) /* Hopping from below */
545                                {
546                                        double interpolation_value;
547                                        (ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
548                                        rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
549                                        rv = dynptarray_addPoint4d(dp, r, 0);
550                                        LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
551                                }
552                        }
553                        /* Add the current vertex to the point array. */
554                        rv = dynptarray_addPoint4d(dp, p, 0);
555                        if ( ordinate_value_p == from || ordinate_value_p == to )
556                        {
557                                added_last_point = 2; /* Added on boundary. */
558                        }
559                        else
560                        {
561                                added_last_point = 1; /* Added inside range. */
562                        }
563                }
564                /* Is this point inside the ordinate range? No. */
565                else
566                {
567                        LWDEBUGF(4, "  added_last_point (%d)", added_last_point);
568                        if ( added_last_point == 1 )
569                        {
570                                /* We're transiting out of the range, so add an interpolated point
571                                *  to the point array at the range boundary. */
572                                double interpolation_value;
573                                (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
574                                rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
575                                rv = dynptarray_addPoint4d(dp, r, 0);
576                                LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
577                        }
578                        else if ( added_last_point == 2 )
579                        {
580                                /* We're out and the last point was on the boundary.
581                                *  If the last point was the near boundary, nothing to do.
582                                *  If it was the far boundary, we need an interpolated point. */
583                                if ( from != to && (
584                                            (ordinate_value_q == from && ordinate_value_p > from) ||
585                                            (ordinate_value_q == to && ordinate_value_p < to) ) )
586                                {
587                                        double interpolation_value;
588                                        (ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
589                                        rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
590                                        rv = dynptarray_addPoint4d(dp, r, 0);
591                                        LWDEBUGF(4, " interpolating between (%g, %g) with interpolation point (%g)", ordinate_value_q, ordinate_value_p, interpolation_value);
592                                }
593                        }
594                        else if ( i && ordinate_value_q < from && ordinate_value_p > to )
595                        {
596                                /* We just hopped over the whole range, from bottom to top,
597                                *  so we need to add *two* interpolated points! */
598                                pa_out = ptarray_construct(hasz, hasm, 2);
599                                /* Interpolate lower point. */
600                                rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
601                                setPoint4d(pa_out, 0, r);
602                                /* Interpolate upper point. */
603                                rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
604                                setPoint4d(pa_out, 1, r);
605                        }
606                        else if ( i && ordinate_value_q > to && ordinate_value_p < from )
607                        {
608                                /* We just hopped over the whole range, from top to bottom,
609                                *  so we need to add *two* interpolated points! */
610                                pa_out = ptarray_construct(hasz, hasm, 2);
611                                /* Interpolate upper point. */
612                                rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
613                                setPoint4d(pa_out, 0, r);
614                                /* Interpolate lower point. */
615                                rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
616                                setPoint4d(pa_out, 1, r);
617                        }
618                        /* We have an extant point-array, save it out to a multi-line. */
619                        if ( dp || pa_out )
620                        {
621                                LWGEOM *oline;
622                                LWDEBUG(4, "saving pointarray to multi-line (1)");
623                                if ( dp )
624                                {
625                                        /* Only one point, so we have to make an lwpoint to hold this
626                                        *  and set the overall output type to a generic collection. */
627                                        if ( dp->pa->npoints == 1 )
628                                        {
629                                                oline = (LWGEOM*)lwpoint_construct(line->SRID, NULL, dp->pa);
630                                                oline->type = lwgeom_makeType(hasz, hasm, hassrid, POINTTYPE);
631                                                lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
632                                        }
633                                        else
634                                        {
635                                                oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
636                                                oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
637                                        }
638                                }
639                                else
640                                {
641                                        oline = (LWGEOM*)lwline_construct(line->SRID, NULL, pa_out);
642                                }
643                                lwgeom_out->ngeoms++;
644                                if ( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
645                                {
646                                        lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
647                                }
648                                else
649                                {
650                                        lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
651                                }
652                                lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
653                                lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
654                                lwgeom_add_bbox((LWGEOM*)lwgeom_out);
655                                if ( dp ) lwfree(dp);
656                                dp = NULL;
657                                if ( pa_out ) pa_out = NULL;
658                        }
659                        added_last_point = 0;
660
661                }
662        }
663
664        /* Still some points left to be saved out. */
665        if ( dp && dp->pa->npoints > 0 )
666        {
667                LWGEOM *oline;
668                LWDEBUG(4, "saving pointarray to multi-line (2)");
669                LWDEBUGF(4, "dp->pa->npoints == %d", dp->pa->npoints);
670                LWDEBUGF(4, "lwgeom_out->ngeoms == %d", lwgeom_out->ngeoms);
671
672                if ( dp->pa->npoints == 1 )
673                {
674                        oline = (LWGEOM*)lwpoint_construct(line->SRID, NULL, dp->pa);
675                        oline->type = lwgeom_makeType(hasz, hasm, hassrid, POINTTYPE);
676                        lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
677                }
678                else
679                {
680                        oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
681                        oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
682                }
683
684                lwgeom_out->ngeoms++;
685                if ( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
686                {
687                        lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
688                }
689                else
690                {
691                        lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
692                }
693                lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
694                lwgeom_drop_bbox((LWGEOM*)lwgeom_out);
695                lwgeom_add_bbox((LWGEOM*)lwgeom_out);
696                if ( dp ) lwfree(dp);
697                dp = NULL;
698        }
699
700        lwfree(p);
701        lwfree(q);
702        lwfree(r);
703
704        if ( lwgeom_out->ngeoms > 0 )
705                return lwgeom_out;
706
707        return NULL;
708
709}
710
711static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";
712
713/*
714** Calculate the geohash, iterating downwards and gaining precision.
715** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>
716** Released under the MIT License.
717*/
718char *geohash_point(double longitude, double latitude, int precision)
719{
720        int is_even=1, i=0;
721        double lat[2], lon[2], mid;
722        char bits[] = {16,8,4,2,1};
723        int bit=0, ch=0;
724        char *geohash = NULL;
725
726        geohash = lwalloc(precision + 1);
727
728        lat[0] = -90.0;
729        lat[1] = 90.0;
730        lon[0] = -180.0;
731        lon[1] = 180.0;
732
733        while (i < precision)
734        {
735                if (is_even)
736                {
737                        mid = (lon[0] + lon[1]) / 2;
738                        if (longitude > mid)
739                        {
740                                ch |= bits[bit];
741                                lon[0] = mid;
742                        }
743                        else
744                        {
745                                lon[1] = mid;
746                        }
747                }
748                else
749                {
750                        mid = (lat[0] + lat[1]) / 2;
751                        if (latitude > mid)
752                        {
753                                ch |= bits[bit];
754                                lat[0] = mid;
755                        }
756                        else
757                        {
758                                lat[1] = mid;
759                        }
760                }
761
762                is_even = !is_even;
763                if (bit < 4)
764                {
765                        bit++;
766                }
767                else
768                {
769                        geohash[i++] = base32[ch];
770                        bit = 0;
771                        ch = 0;
772                }
773        }
774        geohash[i] = 0;
775        return geohash;
776}
777
778int lwgeom_geohash_precision(BOX3D bbox, BOX3D *bounds)
779{
780        double minx, miny, maxx, maxy;
781        double latmax, latmin, lonmax, lonmin;
782        double lonwidth, latwidth;
783        double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;
784        int precision = 0;
785
786        /* Get the bounding box, return error if things don't work out. */
787        minx = bbox.xmin;
788        miny = bbox.ymin;
789        maxx = bbox.xmax;
790        maxy = bbox.ymax;
791
792        if ( minx == maxx && miny == maxy )
793        {
794                /* It's a point. Doubles have 51 bits of precision.
795                ** 2 * 51 / 5 == 20 */
796                return 20;
797        }
798
799        lonmin = -180.0;
800        latmin = -90.0;
801        lonmax = 180.0;
802        latmax = 90.0;
803
804        /* Shrink a world bounding box until one of the edges interferes with the
805        ** bounds of our rectangle. */
806        while ( 1 )
807        {
808                lonwidth = lonmax - lonmin;
809                latwidth = latmax - latmin;
810                latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;
811
812                if ( minx > lonmin + lonwidth / 2.0 )
813                {
814                        lonminadjust = lonwidth / 2.0;
815                }
816                else if ( maxx < lonmax - lonwidth / 2.0 )
817                {
818                        lonmaxadjust = -1 * lonwidth / 2.0;
819                }
820                if ( miny > latmin + latwidth / 2.0 )
821                {
822                        latminadjust = latwidth / 2.0;
823                }
824                else if (maxy < latmax - latwidth / 2.0 )
825                {
826                        latmaxadjust = -1 * latwidth / 2.0;
827                }
828                /* Only adjust if adjustments are legal (we haven't crossed any edges). */
829                if ( (lonminadjust || lonmaxadjust) && (latminadjust || latmaxadjust ) )
830                {
831                        latmin += latminadjust;
832                        lonmin += lonminadjust;
833                        latmax += latmaxadjust;
834                        lonmax += lonmaxadjust;
835                        /* Each adjustment cycle corresponds to 2 bits of storage in the
836                        ** geohash.     */
837                        precision += 2;
838                }
839                else
840                {
841                        break;
842                }
843        }
844
845        /* Save the edges of our bounds, in case someone cares later. */
846        bounds->xmin = lonmin;
847        bounds->xmax = lonmax;
848        bounds->ymin = latmin;
849        bounds->ymax = latmax;
850
851        /* Each geohash character (base32) can contain 5 bits of information.
852        ** We are returning the precision in characters, so here we divide. */
853        return precision / 5;
854}
855
856
857/*
858** Return a geohash string for the geometry. <http://geohash.org>
859** Where the precision is non-positive, calculate a precision based on the
860** bounds of the feature. Big features have loose precision.
861** Small features have tight precision.
862*/
863char *lwgeom_geohash(const LWGEOM *lwgeom, int precision)
864{
865        BOX3D *bbox = NULL;
866        BOX3D precision_bounds;
867        double lat, lon;
868
869        bbox = lwgeom_compute_box3d(lwgeom);
870        if ( ! bbox ) return NULL;
871
872        /* Return error if we are being fed something outside our working bounds */
873        if ( bbox->xmin < -180 || bbox->ymin < -90 || bbox->xmax > 180 || bbox->ymax > 90 )
874        {
875                lwerror("Geohash requires inputs in decimal degrees.");
876                lwfree(bbox);
877                return NULL;
878        }
879
880        /* What is the center of our geometry bounds? We'll use that to
881        ** approximate location. */
882        lon = bbox->xmin + (bbox->xmax - bbox->xmin) / 2;
883        lat = bbox->ymin + (bbox->ymax - bbox->ymin) / 2;
884
885        if ( precision <= 0 )
886        {
887                precision = lwgeom_geohash_precision(*bbox, &precision_bounds);
888        }
889
890        lwfree(bbox);
891
892        /*
893        ** Return the geohash of the center, with a precision determined by the
894        ** extent of the bounds.
895        ** Possible change: return the point at the center of the precision bounds?
896        */
897        return geohash_point(lon, lat, precision);
898}
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
Note: See TracBrowser for help on using the repository browser.