Branch data Line data Source code
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 : : #include "effectivearea.h"
26 : : #include "qgspoint.h"
27 : :
28 : 0 : static MINHEAP initiate_minheap( int npoints )
29 : : {
30 : 0 : MINHEAP tree;
31 : 0 : tree.key_array = ( areanode ** )lwalloc( npoints * sizeof( void * ) );
32 : 0 : tree.maxSize = npoints;
33 : 0 : tree.usedSize = 0;
34 : 0 : return tree;
35 : : }
36 : :
37 : 0 : static void destroy_minheap( MINHEAP tree )
38 : : {
39 : 0 : lwfree( tree.key_array );
40 : 0 : }
41 : :
42 : : /**
43 : : * Calculate the area of a triangle in 2d
44 : : */
45 : 0 : static double triarea2d( const QgsPoint &P1, const QgsPoint &P2, const QgsPoint &P3 )
46 : : {
47 : 0 : return std::fabs( 0.5 * ( ( P1.x() - P2.x() ) * ( P3.y() - P2.y() ) - ( P1.y() - P2.y() ) * ( P3.x() - P2.x() ) ) );
48 : : }
49 : :
50 : : /**
51 : : * Calculate the area of a triangle in 3d space
52 : : */
53 : 0 : static double triarea3d( const QgsPoint &P1, const QgsPoint &P2, const QgsPoint &P3 )
54 : : {
55 : : //LWDEBUG( 2, "Entered triarea3d" );
56 : : double ax, bx, ay, by, az, bz, cx, cy, cz, area;
57 : :
58 : 0 : ax = P1.x() - P2.x();
59 : 0 : bx = P3.x() - P2.x();
60 : 0 : ay = P1.y() - P2.y();
61 : 0 : by = P3.y() - P2.y();
62 : 0 : az = P1.z() - P2.z();
63 : 0 : bz = P3.z() - P2.z();
64 : :
65 : 0 : cx = ay * bz - az * by;
66 : 0 : cy = az * bx - ax * bz;
67 : 0 : cz = ax * by - ay * bx;
68 : :
69 : 0 : area = std::fabs( 0.5 * ( std::sqrt( cx * cx + cy * cy + cz * cz ) ) );
70 : 0 : return area;
71 : : }
72 : :
73 : : /**
74 : : * We create the minheap by ordering the minheap array by the areas in the areanode structs that the minheap keys refer to
75 : : */
76 : 0 : static int cmpfunc( const void *a, const void *b )
77 : : {
78 : 0 : double v1 = ( *( areanode ** )a )->area;
79 : 0 : double v2 = ( *( areanode ** )b )->area;
80 : :
81 : : /* qsort gives unpredictable results when comparing identical values.
82 : : * If two values are the same we force returning the last point in the point array.
83 : : * That way we get the same ordering on different machines and platforms
84 : : */
85 : 0 : if ( v1 == v2 )
86 : 0 : return ( *( areanode ** )a ) - ( *( areanode ** )b );
87 : : else
88 : 0 : return ( v1 > v2 ) ? 1 : -1;
89 : 0 : }
90 : :
91 : : /**
92 : : * Sift Down
93 : : */
94 : 0 : static void down( MINHEAP *tree, areanode *arealist, int parent )
95 : : {
96 : : //LWDEBUG( 2, "Entered down" );
97 : 0 : areanode **treearray = tree->key_array;
98 : 0 : int left = parent * 2 + 1;
99 : 0 : int right = left + 1;
100 : 0 : areanode *tmp = nullptr;
101 : 0 : int swap = parent;
102 : 0 : double leftarea = 0;
103 : 0 : double rightarea = 0;
104 : 0 : double parentarea = ( ( areanode * )treearray[parent] )->area;
105 : :
106 : 0 : if ( left < tree->usedSize )
107 : : {
108 : 0 : leftarea = ( ( areanode * )treearray[left] )->area;
109 : 0 : if ( parentarea > leftarea ) swap = left;
110 : 0 : }
111 : 0 : if ( right < tree->usedSize )
112 : : {
113 : 0 : rightarea = ( ( areanode * )treearray[right] )->area;
114 : 0 : if ( rightarea < parentarea && rightarea < leftarea ) swap = right;
115 : 0 : }
116 : 0 : if ( swap > parent )
117 : : {
118 : : // OK, we have to swap something
119 : 0 : tmp = treearray[parent];
120 : 0 : treearray[parent] = treearray[swap];
121 : : // Update reference
122 : 0 : ( ( areanode * )treearray[parent] )->treeindex = parent;
123 : 0 : treearray[swap] = tmp;
124 : : // Update reference
125 : 0 : ( ( areanode * )treearray[swap] )->treeindex = swap;
126 : 0 : if ( swap < tree->usedSize ) down( tree, arealist, swap );
127 : 0 : }
128 : 0 : }
129 : :
130 : : /**
131 : : * Sift Up
132 : : */
133 : 0 : static void up( MINHEAP *tree, areanode *arealist, int c )
134 : : {
135 : : //LWDEBUG( 2, "Entered up" );
136 : 0 : areanode *tmp = nullptr;
137 : :
138 : : Q_UNUSED( arealist )
139 : :
140 : 0 : areanode **treearray = tree->key_array;
141 : 0 : int parent = ( c - 1 ) / 2;
142 : :
143 : 0 : while ( ( ( areanode * )treearray[c] )->area < ( ( areanode * )treearray[parent] )->area )
144 : : {
145 : : // OK, we have to swap
146 : 0 : tmp = treearray[parent];
147 : 0 : treearray[parent] = treearray[c];
148 : : // Update reference
149 : 0 : ( ( areanode * )treearray[parent] )->treeindex = parent;
150 : 0 : treearray[c] = tmp;
151 : : // Update reference
152 : 0 : ( ( areanode * )treearray[c] )->treeindex = c;
153 : 0 : c = parent;
154 : 0 : parent = ( c - 1 ) / 2;
155 : : }
156 : 0 : }
157 : :
158 : : /**
159 : : * Gets a reference to the point with the smallest effective area from the root of the min heap
160 : : */
161 : 0 : static areanode *minheap_pop( MINHEAP *tree, areanode *arealist )
162 : : {
163 : : //LWDEBUG( 2, "Entered minheap_pop" );
164 : 0 : areanode *res = tree->key_array[0];
165 : :
166 : : // put last value first
167 : 0 : tree->key_array[0] = tree->key_array[( tree->usedSize ) - 1];
168 : 0 : ( ( areanode * )tree->key_array[0] )->treeindex = 0;
169 : :
170 : 0 : tree->usedSize--;
171 : 0 : down( tree, arealist, 0 );
172 : 0 : return res;
173 : : }
174 : :
175 : : /**
176 : : * The member of the minheap at index idx is changed. Update the tree and make restore the heap property
177 : : */
178 : 0 : static void minheap_update( MINHEAP *tree, areanode *arealist, int idx )
179 : : {
180 : 0 : areanode **treearray = tree->key_array;
181 : 0 : int parent = ( idx - 1 ) / 2;
182 : :
183 : 0 : if ( ( ( areanode * )treearray[idx] )->area < ( ( areanode * )treearray[parent] )->area )
184 : 0 : up( tree, arealist, idx );
185 : : else
186 : 0 : down( tree, arealist, idx );
187 : 0 : }
188 : :
189 : : /**
190 : : * To get the effective area, we have to check what area a point results in when all smaller areas are eliminated
191 : : */
192 : 0 : static void tune_areas( EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld )
193 : : {
194 : : //LWDEBUG( 2, "Entered tune_areas" );
195 : 0 : QgsPoint P1;
196 : 0 : QgsPoint P2;
197 : 0 : QgsPoint P3;
198 : : double area;
199 : 0 : int go_on = 1;
200 : 0 : double check_order_min_area = 0;
201 : :
202 : 0 : int npoints = ea->inpts.size();
203 : : int i;
204 : : int current, before_current, after_current;
205 : :
206 : 0 : MINHEAP tree = initiate_minheap( npoints );
207 : :
208 : : // Add all keys (index in initial_arealist) into minheap array
209 : 0 : for ( i = 0; i < npoints; i++ )
210 : : {
211 : 0 : tree.key_array[i] = ea->initial_arealist + i;
212 : : //LWDEBUGF( 2, "add nr %d, with area %lf, and %lf", i, ea->initial_arealist[i].area, tree.key_array[i]->area );
213 : 0 : }
214 : 0 : tree.usedSize = npoints;
215 : :
216 : : // order the keys by area, small to big
217 : 0 : qsort( tree.key_array, npoints, sizeof( void * ), cmpfunc );
218 : :
219 : : // We have to put references to our tree in our point-list
220 : 0 : for ( i = 0; i < npoints; i++ )
221 : : {
222 : 0 : ( ( areanode * )tree.key_array[i] )->treeindex = i;
223 : : //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 );
224 : 0 : }
225 : :
226 : : // OK, now we have a minHeap, just need to keep it
227 : 0 : i = 0;
228 : 0 : while ( go_on )
229 : : {
230 : : // Get a reference to the point with the currently smallest effective area
231 : 0 : current = minheap_pop( &tree, ea->initial_arealist ) - ea->initial_arealist;
232 : :
233 : : // We have found the smallest area. That is the resulting effective area for the "current" point
234 : 0 : if ( i < npoints - avoid_collaps )
235 : 0 : ea->res_arealist[current] = ea->initial_arealist[current].area;
236 : : else
237 : 0 : ea->res_arealist[current] = FLT_MAX;
238 : :
239 : 0 : if ( ea->res_arealist[current] < check_order_min_area )
240 : 0 : 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 mail at the mailing list. Returned area = %lf, and last area = %lf", ea->res_arealist[current], check_order_min_area );
241 : :
242 : 0 : check_order_min_area = ea->res_arealist[current];
243 : :
244 : : // 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
245 : :
246 : : // Find point before and after
247 : 0 : before_current = ea->initial_arealist[current].prev;
248 : 0 : after_current = ea->initial_arealist[current].next;
249 : :
250 : 0 : P2 = ea->inpts.at( before_current );
251 : 0 : P3 = ea->inpts.at( after_current );
252 : :
253 : : // Check if point before current point is the first in the point array.
254 : 0 : if ( before_current > 0 )
255 : : {
256 : 0 : P1 = ea->inpts.at( ea->initial_arealist[before_current].prev );
257 : :
258 : 0 : if ( ea->is3d )
259 : 0 : area = triarea3d( P1, P2, P3 );
260 : : else
261 : 0 : area = triarea2d( P1, P2, P3 );
262 : :
263 : 0 : ea->initial_arealist[before_current].area = FP_MAX( area, ea->res_arealist[current] );
264 : 0 : minheap_update( &tree, ea->initial_arealist, ea->initial_arealist[before_current].treeindex );
265 : 0 : }
266 : 0 : if ( after_current < npoints - 1 ) // Check if point after current point is the last in the point array.
267 : : {
268 : 0 : P1 = P2;
269 : 0 : P2 = P3;
270 : 0 : P3 = ea->inpts.at( ea->initial_arealist[after_current].next );
271 : :
272 : 0 : if ( ea->is3d )
273 : 0 : area = triarea3d( P1, P2, P3 );
274 : : else
275 : 0 : area = triarea2d( P1, P2, P3 );
276 : :
277 : 0 : ea->initial_arealist[after_current].area = FP_MAX( area, ea->res_arealist[current] );
278 : 0 : minheap_update( &tree, ea->initial_arealist, ea->initial_arealist[after_current].treeindex );
279 : 0 : }
280 : :
281 : : // rearrange the nodes so the eliminated point will be ignored on the next run
282 : 0 : ea->initial_arealist[before_current].next = ea->initial_arealist[current].next;
283 : 0 : ea->initial_arealist[after_current ].prev = ea->initial_arealist[current].prev;
284 : :
285 : : // Check if we are finished
286 : 0 : if ( ( !set_area && ea->res_arealist[current] > trshld ) || ( ea->initial_arealist[0].next == ( npoints - 1 ) ) )
287 : 0 : go_on = 0;
288 : :
289 : 0 : i++;
290 : : }
291 : 0 : destroy_minheap( tree );
292 : 0 : }
293 : :
294 : : /**
295 : : * We calculate the effective area for the first time
296 : : */
297 : 0 : void ptarray_calc_areas( EFFECTIVE_AREAS *ea, int avoid_collaps, int set_area, double trshld )
298 : : {
299 : : //LWDEBUG( 2, "Entered ptarray_calc_areas" );
300 : : int i;
301 : 0 : int npoints = ea->inpts.size();
302 : : double area;
303 : :
304 : 0 : QgsPoint P1;
305 : 0 : QgsPoint P2;
306 : 0 : QgsPoint P3;
307 : 0 : P1 = ea->inpts.at( 0 );
308 : 0 : P2 = ea->inpts.at( 1 );
309 : :
310 : : // The first and last point shall always have the maximum effective area. We use float max to not make trouble for bbox
311 : 0 : ea->initial_arealist[0].area = ea->initial_arealist[npoints - 1].area = FLT_MAX;
312 : 0 : ea->res_arealist[0] = ea->res_arealist[npoints - 1] = FLT_MAX;
313 : :
314 : 0 : ea->initial_arealist[0].next = 1;
315 : 0 : ea->initial_arealist[0].prev = 0;
316 : :
317 : 0 : for ( i = 1; i < ( npoints ) - 1; i++ )
318 : : {
319 : 0 : ea->initial_arealist[i].next = i + 1;
320 : 0 : ea->initial_arealist[i].prev = i - 1;
321 : 0 : P3 = ea->inpts.at( i + 1 );
322 : :
323 : 0 : if ( ea->is3d )
324 : 0 : area = triarea3d( P1, P2, P3 );
325 : : else
326 : 0 : area = triarea2d( P1, P2, P3 );
327 : :
328 : : //LWDEBUGF( 4, "Write area %lf to point %d on address %p", area, i, &( ea->initial_arealist[i].area ) );
329 : 0 : ea->initial_arealist[i].area = area;
330 : 0 : P1 = P2;
331 : 0 : P2 = P3;
332 : 0 : }
333 : 0 : ea->initial_arealist[npoints - 1].next = npoints - 1;
334 : 0 : ea->initial_arealist[npoints - 1].prev = npoints - 2;
335 : :
336 : 0 : for ( i = 1; i < ( npoints ) - 1; i++ )
337 : : {
338 : 0 : ea->res_arealist[i] = FLT_MAX;
339 : 0 : }
340 : 0 : tune_areas( ea, avoid_collaps, set_area, trshld );
341 : 0 : }
|
