1 /**********************************************************************
2 *
3 * PostGIS - Spatial Types for PostgreSQL
4 * http://postgis.net
5 *
6 * PostGIS 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 * PostGIS 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 PostGIS. If not, see <http://www.gnu.org/licenses/>.
18 *
19 **********************************************************************
20 *
21 * Copyright 2014 Nicklas Avén
22 *
23 **********************************************************************/
24
25
26 #include "effectivearea.h"
27
28
29 EFFECTIVE_AREAS*
initiate_effectivearea(const POINTARRAY * inpts)30 initiate_effectivearea(const POINTARRAY *inpts)
31 {
32 LWDEBUG(2, "Entered initiate_effectivearea");
33 EFFECTIVE_AREAS *ea;
34 ea=lwalloc(sizeof(EFFECTIVE_AREAS));
35 ea->initial_arealist = lwalloc(inpts->npoints*sizeof(areanode));
36 ea->res_arealist = lwalloc(inpts->npoints*sizeof(double));
37 ea->inpts=inpts;
38 return ea;
39 }
40
41
destroy_effectivearea(EFFECTIVE_AREAS * ea)42 void destroy_effectivearea(EFFECTIVE_AREAS *ea)
43 {
44 lwfree(ea->initial_arealist);
45 lwfree(ea->res_arealist);
46 lwfree(ea);
47 }
48
49
50 static MINHEAP
initiate_minheap(int npoints)51 initiate_minheap(int npoints)
52 {
53 MINHEAP tree;
54 tree.key_array = lwalloc(npoints*sizeof(void*));
55 tree.maxSize=npoints;
56 tree.usedSize=0;
57 return tree;
58 }
59
60
61 static void
destroy_minheap(MINHEAP tree)62 destroy_minheap(MINHEAP tree)
63 {
64 lwfree(tree.key_array);
65 }
66
67
68 /**
69
70 Calculate the area of a triangle in 2d
71 */
triarea2d(const double * P1,const double * P2,const double * P3)72 static double triarea2d(const double *P1, const double *P2, const double *P3)
73 {
74 return fabs(0.5*((P1[0]-P2[0])*(P3[1]-P2[1])-(P1[1]-P2[1])*(P3[0]-P2[0])));
75 }
76
77 /**
78
79 Calculate the area of a triangle in 3d space
80 */
triarea3d(const double * P1,const double * P2,const double * P3)81 static double triarea3d(const double *P1, const double *P2, const double *P3)
82 {
83 LWDEBUG(2, "Entered triarea3d");
84 double ax,bx,ay,by,az,bz,cx,cy,cz, area;
85
86 ax=P1[0]-P2[0];
87 bx=P3[0]-P2[0];
88 ay=P1[1]-P2[1];
89 by=P3[1]-P2[1];
90 az=P1[2]-P2[2];
91 bz=P3[2]-P2[2];
92
93 cx = ay*bz - az*by;
94 cy = az*bx - ax*bz;
95 cz = ax*by - ay*bx;
96
97 area = fabs(0.5*(sqrt(cx*cx+cy*cy+cz*cz)));
98 return area;
99 }
100
101 /**
102
103 We create the minheap by ordering the minheap array by the areas in the areanode structs that the minheap keys refer to
104 */
cmpfunc(const void * a,const void * b)105 static int cmpfunc (const void * a, const void * b)
106 {
107 double v1 = (*(areanode**)a)->area;
108 double v2 = (*(areanode**)b)->area;
109 /*qsort gives unpredictable results when comparing identical values.
110 If two values is the same we force returning the last point in the point array.
111 That way we get the same ordering on different machines and platforms*/
112 if (v1==v2)
113 return (*(areanode**)a)-(*(areanode**)b);
114 else
115 return (v1 > v2) ? 1 : ((v1 < v2) ? -1 : 0);
116 }
117
118
119 /**
120
121 Sift Down
122 */
down(MINHEAP * tree,areanode * arealist,int parent)123 static void down(MINHEAP *tree,areanode *arealist,int parent)
124 {
125 LWDEBUG(2, "Entered down");
126 areanode **treearray=tree->key_array;
127 int left=parent*2+1;
128 int right = left +1;
129 void *tmp;
130 int swap=parent;
131 double leftarea=0;
132 double rightarea=0;
133
134 double parentarea=((areanode*) treearray[parent])->area;
135
136 if(left<tree->usedSize)
137 {
138 leftarea=((areanode*) treearray[left])->area;
139 if(parentarea>leftarea)
140 swap=left;
141 }
142 if(right<tree->usedSize)
143 {
144 rightarea=((areanode*) treearray[right])->area;
145 if(rightarea<parentarea&&rightarea<leftarea)
146 swap=right;
147 }
148 if(swap>parent)
149 {
150 /*ok, we have to swap something*/
151 tmp=treearray[parent];
152 treearray[parent]=treearray[swap];
153 /*Update reference*/
154 ((areanode*) treearray[parent])->treeindex=parent;
155 treearray[swap]=tmp;
156 /*Update reference*/
157 ((areanode*) treearray[swap])->treeindex=swap;
158 if(swap<tree->usedSize)
159 down(tree,arealist,swap);
160 }
161 return;
162 }
163
164
165 /**
166
167 Sift Up
168 */
up(MINHEAP * tree,areanode * e,int c)169 static void up(MINHEAP *tree, __attribute__((__unused__)) areanode *e,int c)
170 {
171 LWDEBUG(2, "Entered up");
172 void *tmp;
173
174 areanode **treearray=tree->key_array;
175
176 int parent=floor((c-1)/2);
177
178 while(((areanode*) treearray[c])->area<((areanode*) treearray[parent])->area)
179 {
180 /*ok, we have to swap*/
181 tmp=treearray[parent];
182 treearray[parent]=treearray[c];
183 /*Update reference*/
184 ((areanode*) treearray[parent])->treeindex=parent;
185 treearray[c]=tmp;
186 /*Update reference*/
187 ((areanode*) treearray[c])->treeindex=c;
188 c=parent;
189 parent=floor((c-1)/2);
190 }
191 return;
192 }
193
194
195 /**
196
197 Get a reference to the point with the smallest effective area from the root of the min heap
198 */
minheap_pop(MINHEAP * tree,areanode * arealist)199 static areanode* minheap_pop(MINHEAP *tree,areanode *arealist )
200 {
201 LWDEBUG(2, "Entered minheap_pop");
202 areanode *res = tree->key_array[0];
203
204 /*put last value first*/
205 tree->key_array[0]=tree->key_array[(tree->usedSize)-1];
206 ((areanode*) tree->key_array[0])->treeindex=0;
207
208 tree->usedSize--;
209 down(tree,arealist,0);
210 return res;
211 }
212
213
214 /**
215
216 The member of the minheap at index idx is changed. Update the tree and make restore the heap property
217 */
minheap_update(MINHEAP * tree,areanode * arealist,int idx)218 static void minheap_update(MINHEAP *tree,areanode *arealist , int idx)
219 {
220 areanode **treearray=tree->key_array;
221 int parent=floor((idx-1)/2);
222
223 if(((areanode*) treearray[idx])->area<((areanode*) treearray[parent])->area)
224 up(tree,arealist,idx);
225 else
226 down(tree,arealist,idx);
227 return;
228 }
229
230 /**
231
232 To get the effective area, we have to check what area a point results in when all smaller areas are eliminated
233 */
tune_areas(EFFECTIVE_AREAS * ea,int avoid_collaps,int set_area,double trshld)234 static void tune_areas(EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld)
235 {
236 LWDEBUG(2, "Entered tune_areas");
237 const double *P1;
238 const double *P2;
239 const double *P3;
240 double area;
241 int go_on=1;
242 double check_order_min_area = 0;
243
244 int npoints=ea->inpts->npoints;
245 int i;
246 int current, before_current, after_current;
247
248 MINHEAP tree = initiate_minheap(npoints);
249
250 int is3d = FLAGS_GET_Z(ea->inpts->flags);
251
252
253 /*Add all keys (index in initial_arealist) into minheap array*/
254 for (i=0;i<npoints;i++)
255 {
256 tree.key_array[i]=ea->initial_arealist+i;
257 LWDEBUGF(2, "add nr %d, with area %lf, and %lf",i,ea->initial_arealist[i].area, tree.key_array[i]->area );
258 }
259 tree.usedSize=npoints;
260
261 /*order the keys by area, small to big*/
262 qsort(tree.key_array, npoints, sizeof(void*), cmpfunc);
263
264 /*We have to put references to our tree in our point-list*/
265 for (i=0;i<npoints;i++)
266 {
267 ((areanode*) tree.key_array[i])->treeindex=i;
268 LWDEBUGF(4,"Check ordering qsort gives, area=%lf and belong to point %d",((areanode*) tree.key_array[i])->area, tree.key_array[i]-ea->initial_arealist);
269 }
270 /*Ok, now we have a minHeap, just need to keep it*/
271
272 /*for (i=0;i<npoints-1;i++)*/
273 i=0;
274 while (go_on)
275 {
276 /*Get a reference to the point with the currently smallest effective area*/
277 current=minheap_pop(&tree, ea->initial_arealist)-ea->initial_arealist;
278
279 /*We have found the smallest area. That is the resulting effective area for the "current" point*/
280 if (i<npoints-avoid_collaps)
281 ea->res_arealist[current]=ea->initial_arealist[current].area;
282 else
283 ea->res_arealist[current]=FLT_MAX;
284
285 if(ea->res_arealist[current]<check_order_min_area)
286 lwerror("Oh no, this is a bug. For some reason the minHeap returned our points in the wrong order. Please file a ticket in PostGIS ticket system, or send a mial at the mailing list.Returned area = %lf, and last area = %lf",ea->res_arealist[current],check_order_min_area);
287
288 check_order_min_area=ea->res_arealist[current];
289
290 /*The found smallest area point is now regarded as eliminated and we have to recalculate the area the adjacent (ignoring earlier eliminated points) points gives*/
291
292 /*FInd point before and after*/
293 before_current=ea->initial_arealist[current].prev;
294 after_current=ea->initial_arealist[current].next;
295
296 P2= (double*)getPoint_internal(ea->inpts, before_current);
297 P3= (double*)getPoint_internal(ea->inpts, after_current);
298
299 /*Check if point before current point is the first in the point array. */
300 if(before_current>0)
301 {
302
303 P1= (double*)getPoint_internal(ea->inpts, ea->initial_arealist[before_current].prev);
304 if(is3d)
305 area=triarea3d(P1, P2, P3);
306 else
307 area=triarea2d(P1, P2, P3);
308
309 ea->initial_arealist[before_current].area = FP_MAX(area,ea->res_arealist[current]);
310 minheap_update(&tree, ea->initial_arealist, ea->initial_arealist[before_current].treeindex);
311 }
312 if(after_current<npoints-1)/*Check if point after current point is the last in the point array. */
313 {
314 P1=P2;
315 P2=P3;
316
317 P3= (double*)getPoint_internal(ea->inpts, ea->initial_arealist[after_current].next);
318
319
320 if(is3d)
321 area=triarea3d(P1, P2, P3);
322 else
323 area=triarea2d(P1, P2, P3);
324
325
326 ea->initial_arealist[after_current].area = FP_MAX(area,ea->res_arealist[current]);
327 minheap_update(&tree, ea->initial_arealist, ea->initial_arealist[after_current].treeindex);
328 }
329
330 /*rearrange the nodes so the eliminated point will be ingored on the next run*/
331 ea->initial_arealist[before_current].next = ea->initial_arealist[current].next;
332 ea->initial_arealist[after_current].prev = ea->initial_arealist[current].prev;
333
334 /*Check if we are finished*/
335 if((!set_area && ea->res_arealist[current]>=trshld) || (ea->initial_arealist[0].next==(npoints-1)))
336 go_on=0;
337
338 i++;
339 };
340 destroy_minheap(tree);
341 return;
342 }
343
344
345 /**
346
347 We calculate the effective area for the first time
348 */
ptarray_calc_areas(EFFECTIVE_AREAS * ea,int avoid_collaps,int set_area,double trshld)349 void ptarray_calc_areas(EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld)
350 {
351 LWDEBUG(2, "Entered ptarray_calc_areas");
352 int i;
353 int npoints=ea->inpts->npoints;
354 int is3d = FLAGS_GET_Z(ea->inpts->flags);
355 double area;
356
357 const double *P1;
358 const double *P2;
359 const double *P3;
360
361 P1 = (double*)getPoint_internal(ea->inpts, 0);
362 P2 = (double*)getPoint_internal(ea->inpts, 1);
363
364 /*The first and last point shall always have the maximum effective area. We use float max to not make trouble for bbox*/
365 ea->initial_arealist[0].area=ea->initial_arealist[npoints-1].area=FLT_MAX;
366 ea->res_arealist[0]=ea->res_arealist[npoints-1]=FLT_MAX;
367
368 ea->initial_arealist[0].next=1;
369 ea->initial_arealist[0].prev=0;
370
371 for (i=1;i<(npoints)-1;i++)
372 {
373 ea->initial_arealist[i].next=i+1;
374 ea->initial_arealist[i].prev=i-1;
375 P3 = (double*)getPoint_internal(ea->inpts, i+1);
376
377 if(is3d)
378 area=triarea3d(P1, P2, P3);
379 else
380 area=triarea2d(P1, P2, P3);
381
382 LWDEBUGF(4,"Write area %lf to point %d on address %p",area,i,&(ea->initial_arealist[i].area));
383 ea->initial_arealist[i].area=area;
384 P1=P2;
385 P2=P3;
386
387 }
388 ea->initial_arealist[npoints-1].next=npoints-1;
389 ea->initial_arealist[npoints-1].prev=npoints-2;
390
391 for (i=1;i<(npoints)-1;i++)
392 {
393 ea->res_arealist[i]=FLT_MAX;
394 }
395
396 tune_areas(ea,avoid_collaps,set_area, trshld);
397 return ;
398 }
399
400
401
ptarray_set_effective_area(POINTARRAY * inpts,int avoid_collaps,int set_area,double trshld)402 static POINTARRAY * ptarray_set_effective_area(POINTARRAY *inpts,int avoid_collaps,int set_area, double trshld)
403 {
404 LWDEBUG(2, "Entered ptarray_set_effective_area");
405 uint32_t p;
406 POINT4D pt;
407 EFFECTIVE_AREAS *ea;
408 POINTARRAY *opts;
409 int set_m;
410 if(set_area)
411 set_m=1;
412 else
413 set_m=FLAGS_GET_M(inpts->flags);
414 ea=initiate_effectivearea(inpts);
415
416 opts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), set_m, inpts->npoints);
417
418 ptarray_calc_areas(ea,avoid_collaps,set_area,trshld);
419
420 if(set_area)
421 {
422 /*Only return points with an effective area above the threshold*/
423 for (p=0;p<ea->inpts->npoints;p++)
424 {
425 if(ea->res_arealist[p]>=trshld)
426 {
427 pt=getPoint4d(ea->inpts, p);
428 pt.m=ea->res_arealist[p];
429 ptarray_append_point(opts, &pt, LW_TRUE);
430 }
431 }
432 }
433 else
434 {
435 /*Only return points with an effective area above the threshold*/
436 for (p=0;p<ea->inpts->npoints;p++)
437 {
438 if(ea->res_arealist[p]>=trshld)
439 {
440 pt=getPoint4d(ea->inpts, p);
441 ptarray_append_point(opts, &pt, LW_TRUE);
442 }
443 }
444 }
445 destroy_effectivearea(ea);
446
447 return opts;
448
449 }
450
lwline_set_effective_area(const LWLINE * iline,int set_area,double trshld)451 static LWLINE* lwline_set_effective_area(const LWLINE *iline,int set_area, double trshld)
452 {
453 LWDEBUG(2, "Entered lwline_set_effective_area");
454
455 /* Skip empty case or too small to simplify */
456 if( lwline_is_empty(iline) || iline->points->npoints<3)
457 return lwline_clone(iline);
458
459 int set_m;
460 if(set_area)
461 set_m=1;
462 else
463 set_m=FLAGS_GET_M(iline->flags);
464
465 LWLINE *oline = lwline_construct_empty(iline->srid, FLAGS_GET_Z(iline->flags), set_m);
466
467
468
469 oline = lwline_construct(iline->srid, NULL, ptarray_set_effective_area(iline->points,2,set_area,trshld));
470
471 oline->type = iline->type;
472 return oline;
473
474 }
475
476
lwpoly_set_effective_area(const LWPOLY * ipoly,int set_area,double trshld)477 static LWPOLY* lwpoly_set_effective_area(const LWPOLY *ipoly,int set_area, double trshld)
478 {
479 LWDEBUG(2, "Entered lwpoly_set_effective_area");
480 uint32_t i;
481 int set_m;
482 int avoid_collapse=4;
483 if(set_area)
484 set_m=1;
485 else
486 set_m=FLAGS_GET_M(ipoly->flags);
487 LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), set_m);
488
489 if( lwpoly_is_empty(ipoly) )
490 return opoly; /* should we return NULL instead ? */
491
492 for (i = 0; i < ipoly->nrings; i++)
493 {
494 POINTARRAY *pa = ptarray_set_effective_area(ipoly->rings[i],avoid_collapse,set_area,trshld);
495 /* Add ring to simplified polygon */
496 if(pa->npoints>=4)
497 {
498 if( lwpoly_add_ring(opoly,pa ) == LW_FAILURE )
499 return NULL;
500 }
501 /*Inner rings we allow to collapse and then we remove them*/
502 avoid_collapse=0;
503 }
504
505
506 opoly->type = ipoly->type;
507
508 if( lwpoly_is_empty(opoly))
509 return NULL;
510
511 return opoly;
512
513 }
514
515
lwcollection_set_effective_area(const LWCOLLECTION * igeom,int set_area,double trshld)516 static LWCOLLECTION* lwcollection_set_effective_area(const LWCOLLECTION *igeom,int set_area, double trshld)
517 {
518 LWDEBUG(2, "Entered lwcollection_set_effective_area");
519 uint32_t i;
520 int set_m;
521 if(set_area)
522 set_m=1;
523 else
524 set_m=FLAGS_GET_M(igeom->flags);
525 LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), set_m);
526
527 if( lwcollection_is_empty(igeom) )
528 return out; /* should we return NULL instead ? */
529
530 for( i = 0; i < igeom->ngeoms; i++ )
531 {
532 LWGEOM *ngeom = lwgeom_set_effective_area(igeom->geoms[i],set_area,trshld);
533 if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
534 }
535
536 return out;
537 }
538
539
lwgeom_set_effective_area(const LWGEOM * igeom,int set_area,double trshld)540 LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom,int set_area, double trshld)
541 {
542 LWDEBUG(2, "Entered lwgeom_set_effective_area");
543 switch (igeom->type)
544 {
545 case POINTTYPE:
546 case MULTIPOINTTYPE:
547 return lwgeom_clone(igeom);
548 case LINETYPE:
549 return (LWGEOM*)lwline_set_effective_area((LWLINE*)igeom,set_area, trshld);
550 case POLYGONTYPE:
551 return (LWGEOM*)lwpoly_set_effective_area((LWPOLY*)igeom,set_area, trshld);
552 case MULTILINETYPE:
553 case MULTIPOLYGONTYPE:
554 case COLLECTIONTYPE:
555 return (LWGEOM*)lwcollection_set_effective_area((LWCOLLECTION *)igeom,set_area, trshld);
556 default:
557 lwerror("lwgeom_simplify: unsupported geometry type: %s",lwtype_name(igeom->type));
558 }
559 return NULL;
560 }
561
562