# source:branches/1.4/liblwgeom/lwalgorithm.c@4787

Last change on this file since 4787 was 4787, checked in by pramsey, 7 years ago

Simplify code and improve consistency of linecrossing results (#272)

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