1 /**********************************************************
2 * Version $Id: SurfaceSpecificPoints.cpp 1921 2014-01-09 10:24:11Z oconrad $
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // ta_morphometry //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // SurfaceSpecificPoints.cpp //
17 // //
18 // Copyright (C) 2003 by //
19 // Olaf Conrad //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: oconrad@saga-gis.org //
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54
55
56 ///////////////////////////////////////////////////////////
57 // //
58 // //
59 // //
60 ///////////////////////////////////////////////////////////
61
62 //---------------------------------------------------------
63 #include "SurfaceSpecificPoints.h"
64
65
66 ///////////////////////////////////////////////////////////
67 // //
68 // //
69 // //
70 ///////////////////////////////////////////////////////////
71
72 //---------------------------------------------------------
CSurfaceSpecificPoints(void)73 CSurfaceSpecificPoints::CSurfaceSpecificPoints(void)
74 {
75 CSG_Parameter *pNode;
76
77 Set_Name (_TL("Surface Specific Points"));
78
79 Set_Author (SG_T("(c) 2001 by O.Conrad"));
80
81 Set_Description (_TW(
82 "References:\n"
83 "Peucker, T.K. and Douglas, D.H., 1975:\n"
84 "'Detection of surface-specific points by local parallel processing of discrete terrain elevation data',\n"
85 "Computer Graphics and Image Processing, 4, 375-387\n"
86 ));
87
88 Parameters.Add_Grid(
89 NULL , "ELEVATION" , _TL("Elevation"),
90 _TL(""),
91 PARAMETER_INPUT
92 );
93
94 Parameters.Add_Grid(
95 NULL , "RESULT" , _TL("Result"),
96 _TL(""),
97 PARAMETER_OUTPUT
98 );
99
100 pNode = Parameters.Add_Choice(
101 NULL , "METHOD" , _TL("Method"),
102 _TL("Algorithm for the detection of Surface Specific Points"),
103 CSG_String::Format(SG_T("%s|%s|%s|%s|%s|"),
104 _TL("Mark Highest Neighbour"),
105 _TL("Opposite Neighbours"),
106 _TL("Flow Direction"),
107 _TL("Flow Direction (up and down)"),
108 _TL("Peucker & Douglas")
109 ), 1
110 );
111
112 Parameters.Add_Value(
113 pNode , "THRESHOLD" , _TL("Threshold"),
114 _TL("Threshold for Peucker & Douglas Algorithm"),
115 PARAMETER_TYPE_Double , 2
116 );
117 }
118
119
120 ///////////////////////////////////////////////////////////
121 // //
122 // //
123 // //
124 ///////////////////////////////////////////////////////////
125
126 //---------------------------------------------------------
On_Execute(void)127 bool CSurfaceSpecificPoints::On_Execute(void)
128 {
129 CSG_Grid *pGrid, *pResult;
130
131 pGrid = Parameters("ELEVATION") ->asGrid();
132 pResult = Parameters("RESULT") ->asGrid();
133
134 switch( Parameters("METHOD")->asInt() )
135 {
136 case 0:
137 Do_MarkHighestNB (pGrid, pResult);
138 break;
139
140 case 1:
141 Do_OppositeNB (pGrid, pResult);
142 break;
143
144 case 2:
145 Do_FlowDirection (pGrid, pResult);
146 break;
147
148 case 3:
149 Do_FlowDirection2 (pGrid, pResult);
150 break;
151
152 case 4:
153 Do_PeuckerDouglas (pGrid, pResult, Parameters("THRESHOLD")->asDouble());
154 break;
155 }
156
157 return( true );
158 }
159
160
161 ///////////////////////////////////////////////////////////
162 // //
163 // //
164 // //
165 ///////////////////////////////////////////////////////////
166
167 //---------------------------------------------------------
Do_MarkHighestNB(CSG_Grid * pGrid,CSG_Grid * pResult)168 void CSurfaceSpecificPoints::Do_MarkHighestNB(CSG_Grid *pGrid, CSG_Grid *pResult) // Band & Lammers...
169 {
170 int i, x, y, ix, iy, xlo, ylo, xhi, yhi;
171
172 double lo, hi, z;
173
174 CSG_Grid *clo, *chi;
175
176 clo = SG_Create_Grid(pGrid, SG_DATATYPE_Char);
177 chi = SG_Create_Grid(pGrid, SG_DATATYPE_Char);
178
179 // Pass 1: Auszaehlen...
180 for(y=0; y<Get_NY() && Set_Progress(y); y++)
181 {
182 for(x=0; x<Get_NX(); x++)
183 {
184 lo = hi = pGrid->asDouble(x,y);
185 xhi = xlo = x;
186 yhi = ylo = y;
187
188 for(i=0; i<4; i++)
189 {
190 ix = Get_xTo(i,x);
191 iy = Get_yTo(i,y);
192
193 if( is_InGrid(ix,iy) )
194 {
195 z = pGrid->asDouble(ix,iy);
196
197 if( z > hi )
198 {
199 hi = z;
200 xhi = ix;
201 yhi = iy;
202 }
203 else if( z < lo )
204 {
205 lo = z;
206 xlo = ix;
207 ylo = iy;
208 }
209 }
210 }
211
212 clo->Add_Value(xlo,ylo,1);
213 chi->Add_Value(xhi,yhi,1);
214 }
215 }
216
217 // Pass 2: Setzen...
218 for(y=0; y<Get_NY() && Set_Progress(y); y++)
219 {
220 for(x=0; x<Get_NX(); x++)
221 {
222 if( !chi->asChar(x,y) )
223 {
224 if( !clo->asChar(x,y) )
225 pResult->Set_Value(x,y, 2); // Sattel
226 else
227 pResult->Set_Value(x,y, 1); // Tiefenlinie
228 }
229 else if( !clo->asChar(x,y) )
230 pResult->Set_Value(x,y, -1); // Wasserscheide
231 else
232 pResult->Set_Value(x,y, 0); // Nichts...
233 }
234 }
235
236 delete(clo);
237 delete(chi);
238 }
239
240
241 ///////////////////////////////////////////////////////////
242 // //
243 // //
244 // //
245 ///////////////////////////////////////////////////////////
246
247 //---------------------------------------------------------
Do_OppositeNB(CSG_Grid * pGrid,CSG_Grid * pResult)248 void CSurfaceSpecificPoints::Do_OppositeNB(CSG_Grid *pGrid, CSG_Grid *pResult)
249 {
250 int i, x, y, ix, iy, jx, jy;
251
252 double z, iz, jz;
253
254 CSG_Grid *clo, *chi;
255
256 clo = SG_Create_Grid(pGrid, SG_DATATYPE_Char);
257 chi = SG_Create_Grid(pGrid, SG_DATATYPE_Char);
258
259 // Pass 1: Auszaehlen...
260 for(y=0; y<Get_NY() && Set_Progress(y); y++)
261 {
262 for(x=0; x<Get_NX(); x++)
263 {
264 z = pGrid->asDouble(x,y);
265
266 for(i=0; i<4; i++)
267 {
268 ix = Get_xTo(i,x);
269 iy = Get_yTo(i,y);
270
271 if( is_InGrid(ix,iy) )
272 {
273 jx = Get_xFrom(i,x);
274 jy = Get_yFrom(i,y);
275
276 if( is_InGrid(jx,jy) )
277 {
278 iz = pGrid->asDouble(ix,iy);
279 jz = pGrid->asDouble(jx,jy);
280
281 if( iz>z && jz>z )
282 chi->Add_Value(x,y,1);
283
284 else if( iz<z && jz<z )
285 clo->Add_Value(x,y,1);
286 }
287 }
288 }
289 }
290 }
291
292 // Pass 2: Setzen...
293 for(y=0; y<Get_NY() && Set_Progress(y); y++)
294 {
295 for(x=0; x<Get_NX(); x++)
296 {
297 if( chi->asChar(x,y) )
298 {
299 if( clo->asChar(x,y) )
300 pResult->Set_Value(x,y, 5); // Sattel
301 else
302 pResult->Set_Value(x,y, chi->asChar(x,y) ); // Tiefenlinie
303 }
304 else if( clo->asChar(x,y) )
305 pResult->Set_Value(x,y, - clo->asChar(x,y) ); // Wasserscheide
306 else
307 pResult->Set_Value(x,y, 0); // Nichts...
308 }
309 }
310
311 delete(clo);
312 delete(chi);
313 }
314
315
316 ///////////////////////////////////////////////////////////
317 // //
318 // //
319 // //
320 ///////////////////////////////////////////////////////////
321
322 //---------------------------------------------------------
Do_FlowDirection(CSG_Grid * pGrid,CSG_Grid * pResult)323 void CSurfaceSpecificPoints::Do_FlowDirection(CSG_Grid *pGrid, CSG_Grid *pResult)
324 {
325 bool bLower;
326
327 int x, y, i, ix, iy, xLow, yLow;
328
329 double z, iz, zLow;
330
331 pResult->Assign();
332
333 for(y=0; y<Get_NY() && Set_Progress(y); y++)
334 {
335 for(x=0; x<Get_NX(); x++)
336 {
337 z = pGrid->asDouble(x,y);
338 bLower = false;
339
340 for(i=0; i<8; i++)
341 {
342 ix = Get_xTo(i,x);
343 iy = Get_yTo(i,y);
344
345 if( is_InGrid(ix,iy) )
346 {
347 iz = pGrid->asDouble(ix,iy);
348
349 if(iz<z)
350 {
351 if(!bLower)
352 {
353 bLower = true;
354 zLow = iz;
355 xLow = ix;
356 yLow = iy;
357 }
358 else if(iz<zLow)
359 {
360 zLow = iz;
361 xLow = ix;
362 yLow = iy;
363 }
364 }
365 }
366 }
367
368 if(bLower)
369 {
370 pResult->Add_Value(xLow, yLow, 1);
371 }
372 }
373 }
374 }
375
376
377 ///////////////////////////////////////////////////////////
378 // //
379 // //
380 // //
381 ///////////////////////////////////////////////////////////
382
383 //---------------------------------------------------------
Do_FlowDirection2(CSG_Grid * pGrid,CSG_Grid * pResult)384 void CSurfaceSpecificPoints::Do_FlowDirection2(CSG_Grid *pGrid, CSG_Grid *pResult)
385 {
386 CSG_Grid Grid(*pGrid), Result(*pResult);
387
388 Do_FlowDirection(&Grid, &Result);
389
390 Grid.Invert();
391
392 Do_FlowDirection(&Grid, pResult);
393
394 for(sLong n=0; n<Get_NCells(); n++)
395 {
396 pResult->Add_Value(n, -Result.asInt(n));
397 }
398 }
399
400
401 ///////////////////////////////////////////////////////////
402 // //
403 // //
404 // //
405 ///////////////////////////////////////////////////////////
406
407 //---------------------------------------------------------
Do_PeuckerDouglas(CSG_Grid * pGrid,CSG_Grid * pResult,double Threshold)408 void CSurfaceSpecificPoints::Do_PeuckerDouglas(CSG_Grid *pGrid, CSG_Grid *pResult, double Threshold)
409 {
410 bool wasPlus;
411
412 int x, y, i, ix, iy, nSgn;
413
414 double d, dPlus, dMinus, z, alt[8];
415
416 for(y=0; y<Get_NY() && Set_Progress(y); y++)
417 {
418 for(x=0; x<Get_NX(); x++)
419 {
420 z = pGrid->asDouble(x,y);
421
422 for(i=0; i<8; i++)
423 {
424 ix = Get_xTo(i,x);
425 iy = Get_yTo(i,y);
426
427 if( is_InGrid(ix,iy) )
428 alt[i] = pGrid->asDouble(ix,iy);
429 else
430 alt[i] = z;
431 }
432
433 dPlus = dMinus = 0;
434 nSgn = 0;
435 wasPlus = (alt[7] - z > 0) ? true : false;
436
437 for(i=0; i<8; i++)
438 {
439 d = alt[i] - z;
440
441 if(d>0)
442 {
443 dPlus += d;
444 if(!wasPlus)
445 {
446 nSgn++;
447 wasPlus = true;
448 }
449 }
450 else if(d<0)
451 {
452 dMinus -= d;
453 if(wasPlus)
454 {
455 nSgn++;
456 wasPlus = false;
457 }
458 }
459 }
460
461 i = 0;
462 if(!dPlus) // Peak...
463 i = 9;
464 else if(!dMinus) // Pit
465 i = -9;
466 else if(nSgn==4) // Pass
467 i = 1;
468 else if(nSgn==2)
469 {
470 i = nSgn = 0;
471
472 if(alt[7]>z)
473 {
474 while(alt[i++]>z);
475 do nSgn++; while(alt[i++]<z);
476 }
477 else
478 {
479 while(alt[i++]<z);
480 do nSgn++; while(alt[i++]>z);
481 }
482
483 i = 0;
484
485 if(nSgn==4)
486 {
487 if(dMinus-dPlus > Threshold) // convex break...
488 i = 2;
489 else if(dPlus-dMinus > Threshold) // concave break...
490 i = -2;
491 }
492 else // lines:
493 {
494 if(dMinus-dPlus>0) // Ridge
495 i = 7;
496 else // Channel
497 i = -7;
498 }
499 }
500
501 pResult->Set_Value(x,y,i);
502 }
503 }
504 }
505
506
507 ///////////////////////////////////////////////////////////
508 // //
509 // //
510 // //
511 ///////////////////////////////////////////////////////////
512
513 //---------------------------------------------------------
514