1 /**********************************************************************
2 *
3 * rttopo - topology library
4 * http://git.osgeo.org/gitea/rttopo/librttopo
5 *
6 * rttopo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * rttopo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with rttopo. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright (C) 2012 Sandro Santilli <strk@kbt.io>
22 * Copyright (C) 2001-2006 Refractions Research Inc.
23 *
24 **********************************************************************/
25
26
27
28 /* basic RTLINE functions */
29
30 #include "rttopo_config.h"
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string.h>
34 #include "librttopo_geom_internal.h"
35 #include "rtgeom_log.h"
36
37
38
39 /*
40 * Construct a new RTLINE. points will *NOT* be copied
41 * use SRID=SRID_UNKNOWN for unknown SRID (will have 8bit type's S = 0)
42 */
43 RTLINE *
rtline_construct(const RTCTX * ctx,int srid,RTGBOX * bbox,RTPOINTARRAY * points)44 rtline_construct(const RTCTX *ctx, int srid, RTGBOX *bbox, RTPOINTARRAY *points)
45 {
46 RTLINE *result;
47 result = (RTLINE*) rtalloc(ctx, sizeof(RTLINE));
48
49 RTDEBUG(ctx, 2, "rtline_construct called.");
50
51 result->type = RTLINETYPE;
52
53 result->flags = points->flags;
54 RTFLAGS_SET_BBOX(result->flags, bbox?1:0);
55
56 RTDEBUGF(ctx, 3, "rtline_construct type=%d", result->type);
57
58 result->srid = srid;
59 result->points = points;
60 result->bbox = bbox;
61
62 return result;
63 }
64
65 RTLINE *
rtline_construct_empty(const RTCTX * ctx,int srid,char hasz,char hasm)66 rtline_construct_empty(const RTCTX *ctx, int srid, char hasz, char hasm)
67 {
68 RTLINE *result = rtalloc(ctx, sizeof(RTLINE));
69 result->type = RTLINETYPE;
70 result->flags = gflags(ctx, hasz,hasm,0);
71 result->srid = srid;
72 result->points = ptarray_construct_empty(ctx, hasz, hasm, 1);
73 result->bbox = NULL;
74 return result;
75 }
76
77
rtline_free(const RTCTX * ctx,RTLINE * line)78 void rtline_free(const RTCTX *ctx, RTLINE *line)
79 {
80 if ( ! line ) return;
81
82 if ( line->bbox )
83 rtfree(ctx, line->bbox);
84 if ( line->points )
85 ptarray_free(ctx, line->points);
86 rtfree(ctx, line);
87 }
88
89
printRTLINE(const RTCTX * ctx,RTLINE * line)90 void printRTLINE(const RTCTX *ctx, RTLINE *line)
91 {
92 rtnotice(ctx, "RTLINE {");
93 rtnotice(ctx, " ndims = %i", (int)RTFLAGS_NDIMS(line->flags));
94 rtnotice(ctx, " srid = %i", (int)line->srid);
95 printPA(ctx, line->points);
96 rtnotice(ctx, "}");
97 }
98
99 /* @brief Clone RTLINE object. Serialized point lists are not copied.
100 *
101 * @see ptarray_clone
102 */
103 RTLINE *
rtline_clone(const RTCTX * ctx,const RTLINE * g)104 rtline_clone(const RTCTX *ctx, const RTLINE *g)
105 {
106 RTLINE *ret = rtalloc(ctx, sizeof(RTLINE));
107
108 RTDEBUGF(ctx, 2, "rtline_clone called with %p", g);
109
110 memcpy(ret, g, sizeof(RTLINE));
111
112 ret->points = ptarray_clone(ctx, g->points);
113
114 if ( g->bbox ) ret->bbox = gbox_copy(ctx, g->bbox);
115 return ret;
116 }
117
118 /* Deep clone RTLINE object. RTPOINTARRAY *is* copied. */
119 RTLINE *
rtline_clone_deep(const RTCTX * ctx,const RTLINE * g)120 rtline_clone_deep(const RTCTX *ctx, const RTLINE *g)
121 {
122 RTLINE *ret = rtalloc(ctx, sizeof(RTLINE));
123
124 RTDEBUGF(ctx, 2, "rtline_clone_deep called with %p", g);
125 memcpy(ret, g, sizeof(RTLINE));
126
127 if ( g->bbox ) ret->bbox = gbox_copy(ctx, g->bbox);
128 if ( g->points ) ret->points = ptarray_clone_deep(ctx, g->points);
129 RTFLAGS_SET_READONLY(ret->flags,0);
130
131 return ret;
132 }
133
134
135 void
rtline_release(const RTCTX * ctx,RTLINE * rtline)136 rtline_release(const RTCTX *ctx, RTLINE *rtline)
137 {
138 rtgeom_release(ctx, rtline_as_rtgeom(ctx, rtline));
139 }
140
141 void
rtline_reverse(const RTCTX * ctx,RTLINE * line)142 rtline_reverse(const RTCTX *ctx, RTLINE *line)
143 {
144 if ( rtline_is_empty(ctx, line) ) return;
145 ptarray_reverse(ctx, line->points);
146 }
147
148 RTLINE *
rtline_segmentize2d(const RTCTX * ctx,RTLINE * line,double dist)149 rtline_segmentize2d(const RTCTX *ctx, RTLINE *line, double dist)
150 {
151 RTPOINTARRAY *segmentized = ptarray_segmentize2d(ctx, line->points, dist);
152 if ( ! segmentized ) return NULL;
153 return rtline_construct(ctx, line->srid, NULL, segmentized);
154 }
155
156 /* check coordinate equality */
157 char
rtline_same(const RTCTX * ctx,const RTLINE * l1,const RTLINE * l2)158 rtline_same(const RTCTX *ctx, const RTLINE *l1, const RTLINE *l2)
159 {
160 return ptarray_same(ctx, l1->points, l2->points);
161 }
162
163 /*
164 * Construct a RTLINE from an array of point and line geometries
165 * RTLINE dimensions are large enough to host all input dimensions.
166 */
167 RTLINE *
rtline_from_rtgeom_array(const RTCTX * ctx,int srid,uint32_t ngeoms,RTGEOM ** geoms)168 rtline_from_rtgeom_array(const RTCTX *ctx, int srid, uint32_t ngeoms, RTGEOM **geoms)
169 {
170 int i;
171 int hasz = RT_FALSE;
172 int hasm = RT_FALSE;
173 RTPOINTARRAY *pa;
174 RTLINE *line;
175 RTPOINT4D pt;
176
177 /*
178 * Find output dimensions, check integrity
179 */
180 for (i=0; i<ngeoms; i++)
181 {
182 if ( RTFLAGS_GET_Z(geoms[i]->flags) ) hasz = RT_TRUE;
183 if ( RTFLAGS_GET_M(geoms[i]->flags) ) hasm = RT_TRUE;
184 if ( hasz && hasm ) break; /* Nothing more to learn! */
185 }
186
187 /* ngeoms should be a guess about how many points we have in input */
188 pa = ptarray_construct_empty(ctx, hasz, hasm, ngeoms);
189
190 for ( i=0; i < ngeoms; i++ )
191 {
192 RTGEOM *g = geoms[i];
193
194 if ( rtgeom_is_empty(ctx, g) ) continue;
195
196 if ( g->type == RTPOINTTYPE )
197 {
198 rtpoint_getPoint4d_p(ctx, (RTPOINT*)g, &pt);
199 ptarray_append_point(ctx, pa, &pt, RT_TRUE);
200 }
201 else if ( g->type == RTLINETYPE )
202 {
203 ptarray_append_ptarray(ctx, pa, ((RTLINE*)g)->points, -1);
204 }
205 else
206 {
207 ptarray_free(ctx, pa);
208 rterror(ctx, "rtline_from_ptarray: invalid input type: %s", rttype_name(ctx, g->type));
209 return NULL;
210 }
211 }
212
213 if ( pa->npoints > 0 )
214 line = rtline_construct(ctx, srid, NULL, pa);
215 else {
216 /* Is this really any different from the above ? */
217 ptarray_free(ctx, pa);
218 line = rtline_construct_empty(ctx, srid, hasz, hasm);
219 }
220
221 return line;
222 }
223
224 /*
225 * Construct a RTLINE from an array of RTPOINTs
226 * RTLINE dimensions are large enough to host all input dimensions.
227 */
228 RTLINE *
rtline_from_ptarray(const RTCTX * ctx,int srid,uint32_t npoints,RTPOINT ** points)229 rtline_from_ptarray(const RTCTX *ctx, int srid, uint32_t npoints, RTPOINT **points)
230 {
231 int i;
232 int hasz = RT_FALSE;
233 int hasm = RT_FALSE;
234 RTPOINTARRAY *pa;
235 RTLINE *line;
236 RTPOINT4D pt;
237
238 /*
239 * Find output dimensions, check integrity
240 */
241 for (i=0; i<npoints; i++)
242 {
243 if ( points[i]->type != RTPOINTTYPE )
244 {
245 rterror(ctx, "rtline_from_ptarray: invalid input type: %s", rttype_name(ctx, points[i]->type));
246 return NULL;
247 }
248 if ( RTFLAGS_GET_Z(points[i]->flags) ) hasz = RT_TRUE;
249 if ( RTFLAGS_GET_M(points[i]->flags) ) hasm = RT_TRUE;
250 if ( hasz && hasm ) break; /* Nothing more to learn! */
251 }
252
253 pa = ptarray_construct_empty(ctx, hasz, hasm, npoints);
254
255 for ( i=0; i < npoints; i++ )
256 {
257 if ( ! rtpoint_is_empty(ctx, points[i]) )
258 {
259 rtpoint_getPoint4d_p(ctx, points[i], &pt);
260 ptarray_append_point(ctx, pa, &pt, RT_TRUE);
261 }
262 }
263
264 if ( pa->npoints > 0 )
265 line = rtline_construct(ctx, srid, NULL, pa);
266 else
267 line = rtline_construct_empty(ctx, srid, hasz, hasm);
268
269 return line;
270 }
271
272 /*
273 * Construct a RTLINE from a RTMPOINT
274 */
275 RTLINE *
rtline_from_rtmpoint(const RTCTX * ctx,int srid,const RTMPOINT * mpoint)276 rtline_from_rtmpoint(const RTCTX *ctx, int srid, const RTMPOINT *mpoint)
277 {
278 uint32_t i;
279 RTPOINTARRAY *pa = NULL;
280 RTGEOM *rtgeom = (RTGEOM*)mpoint;
281 RTPOINT4D pt;
282
283 char hasz = rtgeom_has_z(ctx, rtgeom);
284 char hasm = rtgeom_has_m(ctx, rtgeom);
285 uint32_t npoints = mpoint->ngeoms;
286
287 if ( rtgeom_is_empty(ctx, rtgeom) )
288 {
289 return rtline_construct_empty(ctx, srid, hasz, hasm);
290 }
291
292 pa = ptarray_construct(ctx, hasz, hasm, npoints);
293
294 for (i=0; i < npoints; i++)
295 {
296 rt_getPoint4d_p(ctx, mpoint->geoms[i]->point, 0, &pt);
297 ptarray_set_point4d(ctx, pa, i, &pt);
298 }
299
300 RTDEBUGF(ctx, 3, "rtline_from_rtmpoint: constructed pointarray for %d points", mpoint->ngeoms);
301
302 return rtline_construct(ctx, srid, NULL, pa);
303 }
304
305 /**
306 * Returns freshly allocated #RTPOINT that corresponds to the index where.
307 * Returns NULL if the geometry is empty or the index invalid.
308 */
309 RTPOINT*
rtline_get_rtpoint(const RTCTX * ctx,const RTLINE * line,int where)310 rtline_get_rtpoint(const RTCTX *ctx, const RTLINE *line, int where)
311 {
312 RTPOINT4D pt;
313 RTPOINT *rtpoint;
314 RTPOINTARRAY *pa;
315
316 if ( rtline_is_empty(ctx, line) || where < 0 || where >= line->points->npoints )
317 return NULL;
318
319 pa = ptarray_construct_empty(ctx, RTFLAGS_GET_Z(line->flags), RTFLAGS_GET_M(line->flags), 1);
320 pt = rt_getPoint4d(ctx, line->points, where);
321 ptarray_append_point(ctx, pa, &pt, RT_TRUE);
322 rtpoint = rtpoint_construct(ctx, line->srid, NULL, pa);
323 return rtpoint;
324 }
325
326
327 int
rtline_add_rtpoint(const RTCTX * ctx,RTLINE * line,RTPOINT * point,int where)328 rtline_add_rtpoint(const RTCTX *ctx, RTLINE *line, RTPOINT *point, int where)
329 {
330 RTPOINT4D pt;
331 rt_getPoint4d_p(ctx, point->point, 0, &pt);
332
333 if ( ptarray_insert_point(ctx, line->points, &pt, where) != RT_SUCCESS )
334 return RT_FAILURE;
335
336 /* Update the bounding box */
337 if ( line->bbox )
338 {
339 rtgeom_drop_bbox(ctx, rtline_as_rtgeom(ctx, line));
340 rtgeom_add_bbox(ctx, rtline_as_rtgeom(ctx, line));
341 }
342
343 return RT_SUCCESS;
344 }
345
346
347
348 RTLINE *
rtline_removepoint(const RTCTX * ctx,RTLINE * line,uint32_t index)349 rtline_removepoint(const RTCTX *ctx, RTLINE *line, uint32_t index)
350 {
351 RTPOINTARRAY *newpa;
352 RTLINE *ret;
353
354 newpa = ptarray_removePoint(ctx, line->points, index);
355
356 ret = rtline_construct(ctx, line->srid, NULL, newpa);
357 rtgeom_add_bbox(ctx, (RTGEOM *) ret);
358
359 return ret;
360 }
361
362 /*
363 * Note: input will be changed, make sure you have permissions for this.
364 */
365 void
rtline_setPoint4d(const RTCTX * ctx,RTLINE * line,uint32_t index,RTPOINT4D * newpoint)366 rtline_setPoint4d(const RTCTX *ctx, RTLINE *line, uint32_t index, RTPOINT4D *newpoint)
367 {
368 ptarray_set_point4d(ctx, line->points, index, newpoint);
369 /* Update the box, if there is one to update */
370 if ( line->bbox )
371 {
372 rtgeom_drop_bbox(ctx, (RTGEOM*)line);
373 rtgeom_add_bbox(ctx, (RTGEOM*)line);
374 }
375 }
376
377 /**
378 * Re-write the measure ordinate (or add one, if it isn't already there) interpolating
379 * the measure between the supplied start and end values.
380 */
381 RTLINE*
rtline_measured_from_rtline(const RTCTX * ctx,const RTLINE * rtline,double m_start,double m_end)382 rtline_measured_from_rtline(const RTCTX *ctx, const RTLINE *rtline, double m_start, double m_end)
383 {
384 int i = 0;
385 int hasm = 0, hasz = 0;
386 int npoints = 0;
387 double length = 0.0;
388 double length_so_far = 0.0;
389 double m_range = m_end - m_start;
390 double m;
391 RTPOINTARRAY *pa = NULL;
392 RTPOINT3DZ p1, p2;
393
394 if ( rtline->type != RTLINETYPE )
395 {
396 rterror(ctx, "rtline_construct_from_rtline: only line types supported");
397 return NULL;
398 }
399
400 hasz = RTFLAGS_GET_Z(rtline->flags);
401 hasm = 1;
402
403 /* Null points or npoints == 0 will result in empty return geometry */
404 if ( rtline->points )
405 {
406 npoints = rtline->points->npoints;
407 length = ptarray_length_2d(ctx, rtline->points);
408 rt_getPoint3dz_p(ctx, rtline->points, 0, &p1);
409 }
410
411 pa = ptarray_construct(ctx, hasz, hasm, npoints);
412
413 for ( i = 0; i < npoints; i++ )
414 {
415 RTPOINT4D q;
416 RTPOINT2D a, b;
417 rt_getPoint3dz_p(ctx, rtline->points, i, &p2);
418 a.x = p1.x;
419 a.y = p1.y;
420 b.x = p2.x;
421 b.y = p2.y;
422 length_so_far += distance2d_pt_pt(ctx, &a, &b);
423 if ( length > 0.0 )
424 m = m_start + m_range * length_so_far / length;
425 /* #3172, support (valid) zero-length inputs */
426 else if ( length == 0.0 && npoints > 1 )
427 m = m_start + m_range * i / (npoints-1);
428 else
429 m = 0.0;
430 q.x = p2.x;
431 q.y = p2.y;
432 q.z = p2.z;
433 q.m = m;
434 ptarray_set_point4d(ctx, pa, i, &q);
435 p1 = p2;
436 }
437
438 return rtline_construct(ctx, rtline->srid, NULL, pa);
439 }
440
441 RTGEOM*
rtline_remove_repeated_points(const RTCTX * ctx,const RTLINE * rtline,double tolerance)442 rtline_remove_repeated_points(const RTCTX *ctx, const RTLINE *rtline, double tolerance)
443 {
444 RTPOINTARRAY* npts = ptarray_remove_repeated_points_minpoints(ctx, rtline->points, tolerance, 2);
445
446 RTDEBUGF(ctx, 3, "%s: npts %p", __func__, npts);
447
448 return (RTGEOM*)rtline_construct(ctx, rtline->srid,
449 rtline->bbox ? gbox_copy(ctx, rtline->bbox) : 0,
450 npts);
451 }
452
453 int
rtline_is_closed(const RTCTX * ctx,const RTLINE * line)454 rtline_is_closed(const RTCTX *ctx, const RTLINE *line)
455 {
456 if (RTFLAGS_GET_Z(line->flags))
457 return ptarray_is_closed_3d(ctx, line->points);
458
459 return ptarray_is_closed_2d(ctx, line->points);
460 }
461
462 int
rtline_is_trajectory(const RTCTX * ctx,const RTLINE * line)463 rtline_is_trajectory(const RTCTX *ctx, const RTLINE *line)
464 {
465 RTPOINT3DM p;
466 int i, n;
467 double m = -1 * FLT_MAX;
468
469 if ( ! RTFLAGS_GET_M(line->flags) ) {
470 rtnotice(ctx, "Line does not have M dimension");
471 return RT_FALSE;
472 }
473
474 n = line->points->npoints;
475 if ( n < 2 ) return RT_TRUE; /* empty or single-point are "good" */
476
477 for (i=0; i<n; ++i) {
478 rt_getPoint3dm_p(ctx, line->points, i, &p);
479 if ( p.m <= m ) {
480 rtnotice(ctx, "Measure of vertex %d (%g) not bigger than measure of vertex %d (%g)",
481 i, p.m, i-1, m);
482 return RT_FALSE;
483 }
484 m = p.m;
485 }
486
487 return RT_TRUE;
488 }
489
490
491 RTLINE*
rtline_force_dims(const RTCTX * ctx,const RTLINE * line,int hasz,int hasm)492 rtline_force_dims(const RTCTX *ctx, const RTLINE *line, int hasz, int hasm)
493 {
494 RTPOINTARRAY *pdims = NULL;
495 RTLINE *lineout;
496
497 /* Return 2D empty */
498 if( rtline_is_empty(ctx, line) )
499 {
500 lineout = rtline_construct_empty(ctx, line->srid, hasz, hasm);
501 }
502 else
503 {
504 pdims = ptarray_force_dims(ctx, line->points, hasz, hasm);
505 lineout = rtline_construct(ctx, line->srid, NULL, pdims);
506 }
507 lineout->type = line->type;
508 return lineout;
509 }
510
rtline_is_empty(const RTCTX * ctx,const RTLINE * line)511 int rtline_is_empty(const RTCTX *ctx, const RTLINE *line)
512 {
513 if ( !line->points || line->points->npoints < 1 )
514 return RT_TRUE;
515 return RT_FALSE;
516 }
517
518
rtline_count_vertices(const RTCTX * ctx,RTLINE * line)519 int rtline_count_vertices(const RTCTX *ctx, RTLINE *line)
520 {
521 assert(line);
522 if ( ! line->points )
523 return 0;
524 return line->points->npoints;
525 }
526
rtline_simplify(const RTCTX * ctx,const RTLINE * iline,double dist,int preserve_collapsed)527 RTLINE* rtline_simplify(const RTCTX *ctx, const RTLINE *iline, double dist, int preserve_collapsed)
528 {
529 static const int minvertices = 2; /* TODO: allow setting this */
530 RTLINE *oline;
531 RTPOINTARRAY *pa;
532
533 RTDEBUG(ctx, 2, "function called");
534
535 /* Skip empty case */
536 if( rtline_is_empty(ctx, iline) )
537 return NULL;
538
539 pa = ptarray_simplify(ctx, iline->points, dist, minvertices);
540 if ( ! pa ) return NULL;
541
542 /* Make sure single-point collapses have two points */
543 if ( pa->npoints == 1 )
544 {
545 /* Make sure single-point collapses have two points */
546 if ( preserve_collapsed )
547 {
548 RTPOINT4D pt;
549 rt_getPoint4d_p(ctx, pa, 0, &pt);
550 ptarray_append_point(ctx, pa, &pt, RT_TRUE);
551 }
552 /* Return null for collapse */
553 else
554 {
555 ptarray_free(ctx, pa);
556 return NULL;
557 }
558 }
559
560 oline = rtline_construct(ctx, iline->srid, NULL, pa);
561 oline->type = iline->type;
562 return oline;
563 }
564
rtline_length(const RTCTX * ctx,const RTLINE * line)565 double rtline_length(const RTCTX *ctx, const RTLINE *line)
566 {
567 if ( rtline_is_empty(ctx, line) )
568 return 0.0;
569 return ptarray_length(ctx, line->points);
570 }
571
rtline_length_2d(const RTCTX * ctx,const RTLINE * line)572 double rtline_length_2d(const RTCTX *ctx, const RTLINE *line)
573 {
574 if ( rtline_is_empty(ctx, line) )
575 return 0.0;
576 return ptarray_length_2d(ctx, line->points);
577 }
578
579
580
rtline_grid(const RTCTX * ctx,const RTLINE * line,const gridspec * grid)581 RTLINE* rtline_grid(const RTCTX *ctx, const RTLINE *line, const gridspec *grid)
582 {
583 RTLINE *oline;
584 RTPOINTARRAY *opa;
585
586 opa = ptarray_grid(ctx, line->points, grid);
587
588 /* Skip line3d with less then 2 points */
589 if ( opa->npoints < 2 ) return NULL;
590
591 /* TODO: grid bounding box... */
592 oline = rtline_construct(ctx, line->srid, NULL, opa);
593
594 return oline;
595 }
596
597