1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // shapes_points //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // points_thinning.cpp //
14 // //
15 // Copyright (C) 2011 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "points_thinning.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
CPoints_Thinning(void)59 CPoints_Thinning::CPoints_Thinning(void)
60 {
61 Set_Name (_TL("Point Thinning"));
62
63 Set_Author ("O.Conrad (c) 2011");
64
65 Set_Description (_TW(
66 "The Points Thinning tool aggregates points at a level "
67 "that fits the specified resolution. The information of "
68 "those points that become aggregated is based on basic "
69 "statistics, i.e. mean values for coordinates and mean, "
70 "minimum, maximum, standard deviation for the selected "
71 "attribute. Due to the underlying spatial structure "
72 "the quadtree and the raster method lead to differing, "
73 "though comparable results. "
74
75 ));
76
77 Parameters.Add_Shapes("",
78 "POINTS" , _TL("Points"),
79 _TL(""),
80 PARAMETER_INPUT, SHAPE_TYPE_Point
81 );
82
83 Parameters.Add_Table_Field("POINTS",
84 "FIELD" , _TL("Attribute"),
85 _TL("")
86 );
87
88 Parameters.Add_Bool("",
89 "OUTPUT_PC" , _TL("Output to Point Cloud"),
90 _TL(""),
91 false
92 );
93
94 Parameters.Add_Shapes("",
95 "THINNED" , _TL("Thinned Points"),
96 _TL(""),
97 PARAMETER_OUTPUT, SHAPE_TYPE_Point
98 );
99
100 Parameters.Add_PointCloud("",
101 "THINNED_PC", _TL("Thinned Points"),
102 _TL(""),
103 PARAMETER_OUTPUT
104 );
105
106 Parameters.Add_Double("",
107 "RESOLUTION", _TL("Resolution"),
108 _TL(""),
109 1., 0., true
110 );
111
112 Parameters.Add_Choice("",
113 "METHOD" , _TL("Method"),
114 _TL(""),
115 CSG_String::Format("%s|%s",
116 _TL("quadtree"),
117 _TL("raster")
118 ), 1
119 );
120 }
121
122
123 ///////////////////////////////////////////////////////////
124 // //
125 ///////////////////////////////////////////////////////////
126
127 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)128 int CPoints_Thinning::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
129 {
130 if( pParameter->Cmp_Identifier("OUTPUT_PC") )
131 {
132 pParameters->Set_Enabled("THINNED" , pParameter->asBool() == false);
133 pParameters->Set_Enabled("THINNED_PC", pParameter->asBool() == true);
134 }
135
136 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
137 }
138
139
140 ///////////////////////////////////////////////////////////
141 // //
142 ///////////////////////////////////////////////////////////
143
144 //---------------------------------------------------------
On_Execute(void)145 bool CPoints_Thinning::On_Execute(void)
146 {
147 m_pPoints = Parameters("POINTS")->asShapes();
148
149 if( !m_pPoints->is_Valid() )
150 {
151 Error_Set(_TL("invalid points layer"));
152
153 return( false );
154 }
155
156 if( m_pPoints->Get_Count() < 2 )
157 {
158 Error_Set(_TL("nothing to do, there are less than two points in layer"));
159
160 return( false );
161 }
162
163 //-----------------------------------------------------
164 m_Resolution = Parameters("RESOLUTION")->asDouble();
165
166 if( m_Resolution <= 0. )
167 {
168 Error_Set(_TL("resolution has to be greater than zero"));
169
170 return( false );
171 }
172
173 if( m_Resolution >= m_pPoints->Get_Extent().Get_XRange()
174 && m_Resolution >= m_pPoints->Get_Extent().Get_YRange() )
175 {
176 Error_Set(_TL("nothing to do, resolution needs to be set smaller than the points' extent"));
177
178 return( false );
179 }
180
181 //-----------------------------------------------------
182 m_pPoints->Select();
183
184 if( Parameters("OUTPUT_PC")->asBool() )
185 {
186 m_pThinned = Parameters("THINNED_PC")->asShapes();
187 m_pThinned->asPointCloud()->Create();
188 }
189 else
190 {
191 m_pThinned = Parameters("THINNED" )->asShapes();
192 m_pThinned->asShapes ()->Create(SHAPE_TYPE_Point);
193 }
194
195 m_Field = Parameters("FIELD")->asInt();
196
197 m_pThinned->Fmt_Name("%s [%s]", m_pPoints->Get_Name(), m_pPoints->Get_Field_Name(m_Field));
198
199 m_pThinned->Add_Field("Count" , SG_DATATYPE_Int );
200 m_pThinned->Add_Field("Mean" , SG_DATATYPE_Double);
201 m_pThinned->Add_Field("Minimun", SG_DATATYPE_Double);
202 m_pThinned->Add_Field("Maximun", SG_DATATYPE_Double);
203 m_pThinned->Add_Field("StdDev" , SG_DATATYPE_Double);
204
205 //-----------------------------------------------------
206 switch( Parameters("METHOD")->asInt() )
207 {
208 default: if( !QuadTree_Execute(Get_Extent(0)) ) { return( false ); } break;
209 case 1: if( ! Raster_Execute(Get_Extent(1)) ) { return( false ); } break;
210 }
211
212 //-----------------------------------------------------
213 if( m_pThinned->Get_Count() == m_pPoints->Get_Count() )
214 {
215 Message_Add(_TL("no points removed"));
216 }
217 else
218 {
219 Message_Fmt("\n%d %s", m_pPoints->Get_Count() - m_pThinned->Get_Count(), _TL("points removed"));
220 }
221
222 return( true );
223 }
224
225
226 ///////////////////////////////////////////////////////////
227 // //
228 ///////////////////////////////////////////////////////////
229
230 //---------------------------------------------------------
Get_Extent(int Method)231 CSG_Rect CPoints_Thinning::Get_Extent(int Method)
232 {
233 if( Method == 0 )
234 {
235 CSG_Rect Extent(m_pPoints->Get_Extent());
236
237 Extent.Assign(
238 Extent.Get_XCenter() - 0.5 * m_Resolution, Extent.Get_YCenter() - 0.5 * m_Resolution,
239 Extent.Get_XCenter() + 0.5 * m_Resolution, Extent.Get_YCenter() + 0.5 * m_Resolution
240 );
241
242 while( Extent.Intersects(m_pPoints->Get_Extent()) != INTERSECTION_Contains )
243 {
244 Extent.Inflate(200.);
245 }
246
247 return( Extent );
248 }
249
250 //-----------------------------------------------------
251 CSG_Rect Extent(m_pPoints->Get_Extent());
252
253 double dx = 0.5 * m_Resolution * (1 + (int)(Extent.Get_XRange() / m_Resolution));
254 double dy = 0.5 * m_Resolution * (1 + (int)(Extent.Get_YRange() / m_Resolution));
255
256 Extent.Assign(
257 Extent.Get_XCenter() - dx, Extent.Get_YCenter() - dy,
258 Extent.Get_XCenter() + dx, Extent.Get_YCenter() + dy
259 );
260
261 return( Extent );
262 }
263
264
265 ///////////////////////////////////////////////////////////
266 // //
267 ///////////////////////////////////////////////////////////
268
269 //---------------------------------------------------------
Add_Point(double x,double y,int Count,double Mean,double Minimum,double Maximum,double StdDev)270 inline void CPoints_Thinning::Add_Point(double x, double y, int Count, double Mean, double Minimum, double Maximum, double StdDev)
271 {
272 if( m_pThinned->asPointCloud() )
273 {
274 m_pThinned->asPointCloud()->Add_Point (x, y, Mean);
275 m_pThinned->asPointCloud()->Set_Attribute(0, Count );
276 m_pThinned->asPointCloud()->Set_Attribute(1, Mean );
277 m_pThinned->asPointCloud()->Set_Attribute(2, Minimum);
278 m_pThinned->asPointCloud()->Set_Attribute(3, Maximum);
279 m_pThinned->asPointCloud()->Set_Attribute(4, StdDev );
280 }
281 else if( m_pThinned->asShapes() )
282 {
283 CSG_Shape *pPoint = m_pThinned->Add_Shape();
284
285 pPoint->Add_Point(x, y);
286
287 pPoint->Set_Value(0, Count );
288 pPoint->Set_Value(1, Mean );
289 pPoint->Set_Value(2, Minimum);
290 pPoint->Set_Value(3, Maximum);
291 pPoint->Set_Value(4, StdDev );
292 }
293 }
294
295
296 ///////////////////////////////////////////////////////////
297 // //
298 ///////////////////////////////////////////////////////////
299
300 //---------------------------------------------------------
QuadTree_Execute(const CSG_Rect & Extent)301 bool CPoints_Thinning::QuadTree_Execute(const CSG_Rect &Extent)
302 {
303 Process_Set_Text(_TL("initializing..."));
304
305 if( !m_QuadTree.Create(Extent, true) )
306 {
307 Error_Set(_TL("initialization failed"));
308
309 return( false );
310 }
311
312 for(int i=0; i<m_pPoints->Get_Count() && Set_Progress(i, m_pPoints->Get_Count()); i++)
313 {
314 CSG_Shape *pPoint = m_pPoints->Get_Shape(i);
315
316 m_QuadTree.Add_Point(pPoint->Get_Point(0), pPoint->asDouble(m_Field));
317 }
318
319 //-----------------------------------------------------
320 Process_Set_Text(_TL("evaluating..."));
321
322 QuadTree_Get_Points(m_QuadTree.Get_Root_Pointer());
323
324 m_QuadTree.Destroy();
325
326 return( true );
327 }
328
329 //---------------------------------------------------------
QuadTree_Get_Points(CSG_PRQuadTree_Item * pItem)330 void CPoints_Thinning::QuadTree_Get_Points(CSG_PRQuadTree_Item *pItem)
331 {
332 if( pItem )
333 {
334 if( pItem->is_Leaf() )
335 {
336 QuadTree_Add_Point(pItem->asLeaf());
337 }
338 else if( pItem->Get_Size() <= m_Resolution )
339 {
340 QuadTree_Add_Point((CSG_PRQuadTree_Node_Statistics *)pItem);
341 }
342 else for(int i=0; i<4; i++)
343 {
344 QuadTree_Get_Points(pItem->asNode()->Get_Child(i));
345 }
346 }
347 }
348
349 //---------------------------------------------------------
QuadTree_Add_Point(CSG_PRQuadTree_Leaf * pLeaf)350 inline void CPoints_Thinning::QuadTree_Add_Point(CSG_PRQuadTree_Leaf *pLeaf)
351 {
352 if( pLeaf->has_Statistics() )
353 {
354 CSG_PRQuadTree_Leaf_List *pList = (CSG_PRQuadTree_Leaf_List *)pLeaf;
355
356 Add_Point(pLeaf->Get_X(), pLeaf->Get_Y(),
357 (int)pList->Get_Count (), // Count
358 pList->Get_Mean (), // Mean
359 pList->Get_Minimum(), // Minimun
360 pList->Get_Maximum(), // Maximun
361 pList->Get_StdDev () // StdDev
362 );
363 }
364 else
365 {
366 Add_Point(pLeaf->Get_X(), pLeaf->Get_Y(),
367 1 , // Count
368 pLeaf->Get_Z(), // Mean
369 pLeaf->Get_Z(), // Minimun
370 pLeaf->Get_Z(), // Maximun
371 0. // StdDev
372 );
373 }
374 }
375
376 //---------------------------------------------------------
QuadTree_Add_Point(CSG_PRQuadTree_Node_Statistics * pNode)377 inline void CPoints_Thinning::QuadTree_Add_Point(CSG_PRQuadTree_Node_Statistics *pNode)
378 {
379 Add_Point(pNode->Get_X()->Get_Mean(), pNode->Get_Y()->Get_Mean(),
380 (int)pNode->Get_Z()->Get_Count (), // Count
381 pNode->Get_Z()->Get_Mean (), // Mean
382 pNode->Get_Z()->Get_Minimum(), // Minimun
383 pNode->Get_Z()->Get_Maximum(), // Maximun
384 pNode->Get_Z()->Get_StdDev () // StdDev
385 );
386 }
387
388
389 ///////////////////////////////////////////////////////////
390 // //
391 ///////////////////////////////////////////////////////////
392
393 //---------------------------------------------------------
Raster_Execute(const CSG_Rect & Extent)394 bool CPoints_Thinning::Raster_Execute(const CSG_Rect &Extent)
395 {
396 Process_Set_Text(_TL("initializing..."));
397
398 CSG_Grid_System System(m_Resolution, Extent);
399
400 CSG_Grid X (System, SG_DATATYPE_Double );
401 CSG_Grid Y (System, SG_DATATYPE_Double );
402 CSG_Grid N (System, SG_DATATYPE_Word );
403 CSG_Grid Sum (System, SG_DATATYPE_Double );
404 CSG_Grid Sum2(System, SG_DATATYPE_Double );
405 CSG_Grid Min (System, m_pPoints->Get_Field_Type(m_Field));
406 CSG_Grid Max (System, m_pPoints->Get_Field_Type(m_Field));
407
408 if( !X.is_Valid() || !Y.is_Valid() || !N.is_Valid() || !Sum.is_Valid() || !Sum2.is_Valid() || !Min.is_Valid() || !Max.is_Valid() )
409 {
410 Error_Set(_TL("initialization failed"));
411
412 return( false );
413 }
414
415 //---------------------------------------------------------
416 for(int i=0; i<m_pPoints->Get_Count() && Set_Progress(i, m_pPoints->Get_Count()); i++)
417 {
418 int x, y; CSG_Shape *pPoint = m_pPoints->Get_Shape(i);
419
420 if( System.Get_World_to_Grid(x, y, pPoint->Get_Point(0)) )
421 {
422 double Value = pPoint->asDouble(m_Field);
423
424 X .Add_Value(x, y, pPoint->Get_Point(0).x);
425 Y .Add_Value(x, y, pPoint->Get_Point(0).y);
426 N .Add_Value(x, y, 1.);
427 Sum .Add_Value(x, y, Value);
428 Sum2.Add_Value(x, y, Value * Value);
429
430 if( N.asInt(x, y) <= 1 )
431 {
432 Max.Set_Value(x, y, Value);
433 Min.Set_Value(x, y, Value);
434 }
435 else
436 {
437 if( Max.asDouble(x, y) < Value ) { Max.Set_Value(x, y, Value); }
438 if( Min.asDouble(x, y) > Value ) { Min.Set_Value(x, y, Value); }
439 }
440 }
441 }
442
443 //---------------------------------------------------------
444 Process_Set_Text(_TL("evaluating..."));
445
446 for(int y=0; y<System.Get_NY() && Set_Progress(y, System.Get_NY()); y++)
447 {
448 for(int x=0; x<System.Get_NX(); x++)
449 {
450 int n = N.asInt(x, y);
451
452 if( n > 0 )
453 {
454 double Mean = Sum.asDouble(x, y) / n;
455
456 Add_Point(X.asDouble(x, y) / n, Y.asDouble(x, y) / n,
457 n,
458 Mean,
459 Min.asDouble(x, y),
460 Max.asDouble(x, y),
461 sqrt(Sum2.asDouble(x, y) / n - Mean*Mean)
462 );
463 }
464 }
465 }
466
467 //-----------------------------------------------------
468 return( true );
469 }
470
471
472 ///////////////////////////////////////////////////////////
473 // //
474 // //
475 // //
476 ///////////////////////////////////////////////////////////
477
478 //---------------------------------------------------------
479