1 /**********************************************************
2 * Version $Id: wombling.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 // Grid_Filter //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // wombling.cpp //
17 // //
18 // Copyright (C) 2015 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 Hamburg //
47 // Germany //
48 // //
49 ///////////////////////////////////////////////////////////
50
51 //---------------------------------------------------------
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
61 #include "wombling.h"
62
63
64 ///////////////////////////////////////////////////////////
65 // //
66 // //
67 // //
68 ///////////////////////////////////////////////////////////
69
70 //---------------------------------------------------------
CWombling_Base(void)71 CWombling_Base::CWombling_Base(void)
72 {
73 //-----------------------------------------------------
74 Parameters.Add_Double("",
75 "TMAGNITUDE" , _TL("Minimum Magnitude"),
76 _TL("Minimum magnitude as percentile."),
77 90.0, 0.0, true, 100.0, true
78 );
79
80 Parameters.Add_Double("",
81 "TDIRECTION" , _TL("Maximum Angle"),
82 _TL("Maximum angular difference as degree between adjacent segment points."),
83 30.0, 0.0, true, 180.0, true
84 );
85
86 Parameters.Add_Int("",
87 "TNEIGHBOUR" , _TL("Minimum Neighbours"),
88 _TL("Minimum number of neighbouring potential edge cells with similar direction."),
89 1, 0, true
90 );
91
92 Parameters.Add_Choice("",
93 "ALIGNMENT" , _TL("Alignment"),
94 _TL(""),
95 CSG_String::Format(SG_T("%s|%s|"),
96 _TL("between cells"),
97 _TL("on cell")
98 ), 1
99 );
100
101 Parameters.Add_Choice("",
102 "NEIGHBOUR" , _TL("Edge Connectivity"),
103 _TL(""),
104 CSG_String::Format("%s|%s|",
105 _TL("Rooke's case"),
106 _TL("Queen's case")
107 ), 1
108 );
109 }
110
111
112 ///////////////////////////////////////////////////////////
113 // //
114 ///////////////////////////////////////////////////////////
115
116 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)117 int CWombling_Base::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
118 {
119 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
120 }
121
122
123 ///////////////////////////////////////////////////////////
124 // //
125 ///////////////////////////////////////////////////////////
126
127 //---------------------------------------------------------
Initialize(CSG_Grid Gradient[2],CSG_Grid * pEdges)128 bool CWombling_Base::Initialize(CSG_Grid Gradient[2], CSG_Grid *pEdges)
129 {
130 m_Neighbour = Parameters( "NEIGHBOUR")->asInt() == 0 ? 2 : 1;
131 m_minNeighbours = Parameters("TNEIGHBOUR")->asInt();
132 m_maxAngle = Parameters("TDIRECTION")->asDouble() * M_DEG_TO_RAD;
133
134 //-----------------------------------------------------
135 CSG_Grid_System System;
136
137 if( Parameters("ALIGNMENT")->asInt() == 1 )
138 {
139 System = Get_System();
140 }
141 else
142 {
143 System.Assign(Get_Cellsize(),
144 Get_XMin() + 0.5 * Get_Cellsize(),
145 Get_YMin() + 0.5 * Get_Cellsize(),
146 Get_NX() - 1, Get_NY() - 1
147 );
148 }
149
150 //-----------------------------------------------------
151 Gradient[0].Create(System); Gradient[0].Set_NoData_Value(-1.0);
152 Gradient[1].Create(System); Gradient[1].Set_NoData_Value(-1.0);
153
154 pEdges->Create(System, SG_DATATYPE_Char); pEdges->Set_NoData_Value(0.0);
155
156 //-----------------------------------------------------
157 return( true );
158 }
159
160
161 ///////////////////////////////////////////////////////////
162 // //
163 ///////////////////////////////////////////////////////////
164
165 //---------------------------------------------------------
Get_Gradient(CSG_Grid Gradient[2],CSG_Grid * pFeature,bool bOrientation)166 bool CWombling_Base::Get_Gradient(CSG_Grid Gradient[2], CSG_Grid *pFeature, bool bOrientation)
167 {
168 int Alignment = Parameters("ALIGNMENT")->asInt() == 1;
169
170 Gradient[0].Fmt_Name("%s [%s]", pFeature->Get_Name(), _TL("Magnitude"));
171 Gradient[1].Fmt_Name("%s [%s]", pFeature->Get_Name(), _TL("Direction"));
172
173 //-----------------------------------------------------
174 for(int y=0; y<Gradient[0].Get_NY() && Set_Progress(y, Gradient[0].Get_NY()); y++)
175 {
176 #pragma omp parallel for
177 for(int x=0; x<Gradient[0].Get_NX(); x++)
178 {
179 double Slope, Aspect, z[4];
180
181 if( Alignment == 1 )
182 {
183 if( !pFeature->Get_Gradient(x, y, Slope, Aspect) )
184 {
185 Slope = -1.0;
186 }
187 }
188 else if( pFeature->is_NoData(x , y )
189 || pFeature->is_NoData(x + 1, y )
190 || pFeature->is_NoData(x , y + 1)
191 || pFeature->is_NoData(x + 1, y + 1) )
192 {
193 Slope = -1.0;
194 }
195 else
196 {
197 z[0] = pFeature->asDouble(x , y );
198 z[1] = pFeature->asDouble(x + 1, y );
199 z[2] = pFeature->asDouble(x , y + 1);
200 z[3] = pFeature->asDouble(x + 1, y + 1);
201
202 _Get_Gradient_2x2(z, Slope, Aspect);
203 }
204
205 //---------------------------------------------
206 if( Slope < 0.0 )
207 {
208 Gradient[0].Set_NoData(x, y);
209 }
210 else
211 {
212 Gradient[0].Set_Value(x, y, Slope);
213 }
214
215 if( Aspect < 0.0 )
216 {
217 Gradient[1].Set_NoData(x, y);
218 }
219 else
220 {
221 Gradient[1].Set_Value(x, y, bOrientation && Aspect > M_PI_180 ? Aspect - M_PI_180 : Aspect);
222 }
223 }
224 }
225
226 return( true );
227 }
228
229 //---------------------------------------------------------
_Get_Gradient_2x2(double z[4],double & Slope,double & Aspect)230 inline void CWombling_Base::_Get_Gradient_2x2(double z[4], double &Slope, double &Aspect)
231 {
232 double a = (-z[2] + z[3] + z[0] - z[1]) * 0.5 / Get_Cellsize();
233 double b = ( z[2] + z[3] - z[0] - z[1]) * 0.5 / Get_Cellsize();
234
235 Slope = sqrt(a*a + b*b);
236 Aspect = a != 0.0 ? M_PI_180 + atan2(b, a)
237 : b > 0.0 ? M_PI_270
238 : b < 0.0 ? M_PI_090 : -1;
239 }
240
241
242 ///////////////////////////////////////////////////////////
243 // //
244 ///////////////////////////////////////////////////////////
245
246 //---------------------------------------------------------
Get_Edge_Cells(CSG_Grid * Gradient,CSG_Grid * pEdges)247 bool CWombling_Base::Get_Edge_Cells(CSG_Grid *Gradient, CSG_Grid *pEdges)
248 {
249 CSG_Shapes *pPoints = Parameters("EDGE_POINTS") ? Parameters("EDGE_POINTS")->asShapes() : NULL;
250
251 if( pPoints )
252 {
253 pPoints->Create(SHAPE_TYPE_Point, CSG_String::Format("%s.%s", Parameters("FEATURE")->asGrid()->Get_Name(), _TL("Edges")));
254
255 pPoints->Add_Field("ID" , SG_DATATYPE_Int );
256 pPoints->Add_Field("MAGNITUDE" , SG_DATATYPE_Double);
257 pPoints->Add_Field("DIRECTION" , SG_DATATYPE_Double);
258 pPoints->Add_Field("NEIGHBOURS", SG_DATATYPE_Int );
259 }
260
261 //-----------------------------------------------------
262 int y;
263
264 Lock_Create();
265
266 //-----------------------------------------------------
267 // 1. magnitude
268
269 double Threshold = Gradient[0].Get_Percentile(Parameters("TMAGNITUDE")->asDouble());
270
271 for(y=0; y<Gradient[0].Get_NY() && Set_Progress(y, Gradient[0].Get_NY()); y++)
272 {
273 #pragma omp parallel for
274 for(int x=0; x<Gradient[0].Get_NX(); x++)
275 {
276 if( !Gradient[1].is_NoData(x, y) && !Gradient[0].is_NoData(x, y) && Gradient[0].asDouble(x, y) >= Threshold )
277 {
278 Lock_Set(x, y);
279 }
280 }
281 }
282
283 //-----------------------------------------------------
284 // 2. neighbours and direction
285
286 for(y=0; y<Gradient[0].Get_NY() && Set_Progress(y, Gradient[0].Get_NY()); y++)
287 {
288 for(int x=0; x<Gradient[0].Get_NX(); x++)
289 {
290 int n = _is_Edge_Cell(Gradient, x, y);
291
292 if( n >= m_minNeighbours )
293 {
294 pEdges->Set_Value(x, y, n);
295
296 if( pPoints )
297 {
298 CSG_Shape *pPoint = pPoints->Add_Shape();
299
300 pPoint->Set_Point(Gradient[0].Get_System().Get_Grid_to_World(x, y), 0);
301 pPoint->Set_Value(0, pPoints->Get_Count());
302 pPoint->Set_Value(1, Gradient[0].asDouble(x, y));
303 pPoint->Set_Value(2, Gradient[1].asDouble(x, y) * M_RAD_TO_DEG);
304 pPoint->Set_Value(3, n);
305 }
306 }
307 else
308 {
309 pEdges->Set_NoData(x, y);
310 }
311 }
312 }
313
314 //-----------------------------------------------------
315 Lock_Destroy();
316
317 return( true );
318 }
319
320 //---------------------------------------------------------
_is_Edge_Cell(CSG_Grid Gradient[2],int x,int y)321 int CWombling_Base::_is_Edge_Cell(CSG_Grid Gradient[2], int x, int y)
322 {
323 int n = 0;
324
325 if( Lock_Get(x, y) )
326 {
327 for(int i=0; i<8; i+=m_Neighbour)
328 {
329 int ix = Get_xTo(i, x);
330 int iy = Get_yTo(i, y);
331
332 if( Gradient[0].is_InGrid(ix, iy) && Lock_Get(ix, iy) && SG_Get_Angle_Difference(Gradient[1].asDouble(x, y), Gradient[1].asDouble(ix, iy)) <= m_maxAngle )
333 {
334 n++;
335 }
336 }
337 }
338
339 return( n );
340 }
341
342
343 ///////////////////////////////////////////////////////////
344 // //
345 ///////////////////////////////////////////////////////////
346
347 //---------------------------------------------------------
Get_Edge_Lines(CSG_Grid Gradient[2],CSG_Grid * pEdge)348 bool CWombling_Base::Get_Edge_Lines(CSG_Grid Gradient[2], CSG_Grid *pEdge)
349 {
350 CSG_Shapes *pLines = Parameters("EDGE_LINES") ? Parameters("EDGE_LINES")->asShapes() : NULL;
351
352 if( !pLines )
353 {
354 return( false );
355 }
356
357 pLines->Create(SHAPE_TYPE_Line, CSG_String::Format("%s %s", Parameters("FEATURE")->asGrid()->Get_Name(), _TL("Edges")));
358
359 pLines->Add_Field("ID" , SG_DATATYPE_Int);
360 pLines->Add_Field("ANGLE", SG_DATATYPE_Double);
361
362 for(int y=0; y<Gradient[0].Get_NY() && Set_Progress(y); y++)
363 {
364 for(int x=0; x<Gradient[0].Get_NX(); x++)
365 {
366 if( !pEdge->is_NoData(x, y) )
367 {
368 for(int i=0; i<4; i+=2)
369 {
370 int ix = Get_xTo(i, x);
371 int iy = Get_yTo(i, y);
372
373 if( pEdge->is_InGrid(ix, iy) )
374 {
375 double diff = SG_Get_Angle_Difference(Gradient[1].asDouble(x, y), Gradient[1].asDouble(ix, iy));
376
377 if( diff <= m_maxAngle )
378 {
379 CSG_Shape *pLine = pLines->Add_Shape();
380
381 pLine->Add_Point(Gradient[0].Get_System().Get_Grid_to_World( x, y));
382 pLine->Add_Point(Gradient[0].Get_System().Get_Grid_to_World(ix, iy));
383 pLine->Set_Value(0, pLines->Get_Count());
384 pLine->Set_Value(1, diff * M_RAD_TO_DEG);
385 }
386 }
387 }
388 }
389 }
390 }
391
392 return( true );
393 }
394
395
396 ///////////////////////////////////////////////////////////
397 // //
398 // //
399 // //
400 ///////////////////////////////////////////////////////////
401
402 //---------------------------------------------------------
CWombling(void)403 CWombling::CWombling(void)
404 {
405 //-----------------------------------------------------
406 Set_Name (_TL("Wombling (Edge Detection)"));
407
408 Set_Author ("O.Conrad (c) 2015");
409
410 Set_Description (_TW(
411 "Continuous Wombling for edge detection. Uses magnitude of gradient "
412 "to detect edges between adjacent cells. Edge segments connect such "
413 "edges, when the difference of their gradient directions is below given threshold."
414 ));
415
416 Add_Reference(
417 "Fitzpatrick, M.C., Preisser, E.L., Porter, A., Elkinton, J., Waller, L.A., Carlin, B.P., Ellison, A.M.", "2010",
418 "Ecological boundary detection using Bayesian areal wombling", "Ecology 91(12): 3448-3455. doi:10.1890/10-0807.1."
419 );
420
421 Add_Reference(
422 "Fortin, M.-J. and Dale, M.R.T", "2005",
423 "Spatial Analysis - A Guide for Ecologists", "Cambridge University Press."
424 );
425
426 //-----------------------------------------------------
427 Parameters.Add_Grid("",
428 "FEATURE" , _TL("Feature"),
429 _TL(""),
430 PARAMETER_INPUT
431 );
432
433 Parameters.Add_Shapes("",
434 "EDGE_POINTS" , _TL("Edge Points"),
435 _TL(""),
436 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Point
437 );
438
439 Parameters.Add_Shapes("",
440 "EDGE_LINES" , _TL("Edge Segments"),
441 _TL(""),
442 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Line
443 );
444
445 Parameters.Add_Bool("",
446 "GRADIENTS_OUT" , _TL("Output of Gradients"),
447 _TL(""),
448 false
449 );
450
451 Parameters.Add_Grid_List("",
452 "GRADIENTS" , _TL("Gradients"),
453 _TL(""),
454 PARAMETER_OUTPUT_OPTIONAL, false
455 );
456 }
457
458
459 ///////////////////////////////////////////////////////////
460 // //
461 ///////////////////////////////////////////////////////////
462
463 //---------------------------------------------------------
On_Execute(void)464 bool CWombling::On_Execute(void)
465 {
466 //-----------------------------------------------------
467 CSG_Grid Gradient[2], Edges;
468
469 if( !Initialize(Gradient, &Edges) )
470 {
471 return( false );
472 }
473
474 //-----------------------------------------------------
475 CSG_Grid *pFeature = Parameters("FEATURE")->asGrid();
476
477 Edges.Fmt_Name("%s [%s]", pFeature->Get_Name(), _TL("Edges"));
478
479 Get_Gradient (Gradient, pFeature, false);
480
481 Get_Edge_Cells (Gradient, &Edges);
482 Get_Edge_Lines (Gradient, &Edges);
483
484 //-----------------------------------------------------
485 if( Parameters("GRADIENTS_OUT")->asBool() )
486 {
487 CSG_Parameter_Grid_List *pGrids = Parameters("GRADIENTS")->asGridList();
488
489 if( pGrids->Get_Grid(0) && pGrids->Get_Grid(0)->Get_System().is_Equal(Gradient[0].Get_System())
490 && pGrids->Get_Grid(1) && pGrids->Get_Grid(1)->Get_System().is_Equal(Gradient[1].Get_System()) )
491 {
492 pGrids->Get_Grid(0)->Assign(&Gradient[0]);
493 pGrids->Get_Grid(1)->Assign(&Gradient[1]);
494 }
495 else
496 {
497 pGrids->Del_Items();
498
499 pGrids->Add_Item(SG_Create_Grid(Gradient[0]));
500 pGrids->Add_Item(SG_Create_Grid(Gradient[1]));
501 }
502 }
503
504 //-----------------------------------------------------
505 return( true );
506 }
507
508
509 ///////////////////////////////////////////////////////////
510 // //
511 // //
512 // //
513 ///////////////////////////////////////////////////////////
514
515 //---------------------------------------------------------
CWombling_MultiFeature(void)516 CWombling_MultiFeature::CWombling_MultiFeature(void)
517 {
518 //-----------------------------------------------------
519 Set_Name (_TL("Wombling for Multiple Features (Edge Detection)"));
520
521 Set_Author ("O.Conrad (c) 2015");
522
523 Set_Description (_TW(
524 "Continuous Wombling for edge detection. Uses magnitude of gradient "
525 "to detect edges between adjacent cells. Edge segments connect such "
526 "edges, when the difference of their gradient directions is below given threshold."
527 ));
528
529 Add_Reference(
530 "Fitzpatrick, M.C., Preisser, E.L., Porter, A., Elkinton, J., Waller, L.A., Carlin, B.P., Ellison, A.M.", "2010",
531 "Ecological boundary detection using Bayesian areal wombling", "Ecology 91(12): 3448-3455. doi:10.1890/10-0807.1."
532 );
533
534 Add_Reference(
535 "Fortin, M.-J. and Dale, M.R.T", "2005",
536 "Spatial Analysis - A Guide for Ecologists", "Cambridge University Press."
537 );
538
539 //-----------------------------------------------------
540 Parameters.Add_Grid_List("",
541 "FEATURES" , _TL("Features"),
542 _TL(""),
543 PARAMETER_INPUT
544 );
545
546 Parameters.Add_Grid_List("",
547 "EDGE_CELLS" , _TL("Edges"),
548 _TL(""),
549 PARAMETER_OUTPUT, false
550 );
551
552 //Parameters.Add_Choice("",
553 // "ASPECT_CMP" , _TL("Direction Difference"),
554 // _TL(""),
555 // CSG_String::Format("%s|%s|",
556 // _TL("direction"),
557 // _TL("orientation")
558 // ), 0
559 //);
560
561 Parameters.Add_Choice("",
562 "OUTPUT_ADD" , _TL("Additional Output"),
563 _TL(""),
564 CSG_String::Format("%s|%s|%s|",
565 _TL("none"),
566 _TL("gradients"),
567 _TL("edge cells")
568 ), 0
569 );
570
571 Parameters.Add_Grid_List("",
572 "OUTPUT" , _TL("Output"),
573 _TL(""),
574 PARAMETER_OUTPUT_OPTIONAL, false
575 );
576
577 Parameters.Add_Bool("",
578 "ZERO_AS_NODATA" , _TL("Zero as No-Data"),
579 _TL(""),
580 true
581 );
582 }
583
584
585 ///////////////////////////////////////////////////////////
586 // //
587 ///////////////////////////////////////////////////////////
588
589 //---------------------------------------------------------
On_Execute(void)590 bool CWombling_MultiFeature::On_Execute(void)
591 {
592 CSG_Parameter_Grid_List *pFeatures = Parameters("FEATURES")->asGridList(), *pOutput = NULL;
593
594 //-----------------------------------------------------
595 CSG_Grid Gradient[2], Edges;
596
597 if( !Initialize(Gradient, &Edges) )
598 {
599 return( false );
600 }
601
602 //-----------------------------------------------------
603 if( Parameters("OUTPUT_ADD")->asInt() )
604 {
605 pOutput = Parameters("OUTPUT")->asGridList();
606
607 for(int i=pOutput->Get_Grid_Count()-1; i>=0; i--)
608 {
609 if( !pOutput->Get_Grid(i)->Get_System().is_Equal(Gradient[0].Get_System()) )
610 {
611 pOutput->Del_Item(i);
612 }
613 }
614 }
615
616 //-----------------------------------------------------
617 CSG_Grid Count, *pCount = Parameters("EDGE_CELLS")->asGridList()->Get_Grid(0);
618
619 if( !pCount || !pCount->Get_System().is_Equal(Gradient[0].Get_System()) )
620 {
621 Parameters("EDGE_CELLS")->asGridList()->Del_Items();
622 Parameters("EDGE_CELLS")->asGridList()->Add_Item(pCount = SG_Create_Grid(Gradient[0].Get_System(), SG_DATATYPE_Char));
623 }
624
625 pCount->Set_Name(_TL("Edge Cells"));
626 pCount->Assign(0.0);
627 pCount->Set_NoData_Value(-1.0);
628
629 //-----------------------------------------------------
630 for(int i=0; i<pFeatures->Get_Grid_Count() && Process_Get_Okay(); i++)
631 {
632 Edges.Fmt_Name("%s [%s]", pFeatures->Get_Grid(i)->Get_Name(), _TL("Edges"));
633
634 Get_Gradient(Gradient, pFeatures->Get_Grid(i), false);
635
636 Get_Edge_Cells(Gradient, &Edges);
637
638 pCount->Add(Edges);
639
640 if( pOutput )
641 {
642 if( !pOutput->Get_Grid(i) )
643 {
644 pOutput->Add_Item(SG_Create_Grid());
645 }
646
647 pOutput->Get_Grid(i)->Create(Parameters("OUTPUT_ADD")->asInt() == 1 ? Gradient[0] : Edges);
648 }
649 }
650
651 if( Parameters("ZERO_AS_NODATA")->asBool() )
652 {
653 pCount->Set_NoData_Value(0.0);
654 }
655
656 //-----------------------------------------------------
657 return( true );
658 }
659
660
661 ///////////////////////////////////////////////////////////
662 // //
663 // //
664 // //
665 ///////////////////////////////////////////////////////////
666
667 //---------------------------------------------------------
668