1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // Grid_Filter //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // geomrec.cpp //
17 // //
18 // Copyright (C) 2013 by //
19 // HfT Stuttgart //
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: johannes.engels@hft-stuttgart.de //
43 // //
44 // contact: Johannes Engels //
45 // Hochschule fuer Technik Stuttgart //
46 // Schellingstr. 24 //
47 // 70174 Stuttgart //
48 // Germany //
49 // //
50 ///////////////////////////////////////////////////////////
51
52 //---------------------------------------------------------
53
54
55 ///////////////////////////////////////////////////////////
56 // //
57 // //
58 // //
59 ///////////////////////////////////////////////////////////
60
61 //---------------------------------------------------------
62 #include "geomrec.h"
63
64 extern "C" {
65 #include "geodesic_morph_rec/storeorg.h"
66 #include "geodesic_morph_rec/geovinc.h"
67 }
68
69 ///////////////////////////////////////////////////////////
70 // //
71 // //
72 // //
73 ///////////////////////////////////////////////////////////
74
75 //---------------------------------------------------------
CGeomrec(void)76 CGeomrec::CGeomrec(void)
77 {
78 Set_Name (_TL("Geodesic Morphological Reconstruction"));
79
80 Set_Author (SG_T("HfT Stuttgart (c) 2013"));
81
82 Set_Description (_TW(
83 "Geodesic morphological reconstruction according to \n"
84 "L. Vincent (1993): Morphological Grayscale Reconstruction in Image Analysis: "
85 "Applications and Efficient Algorithms. "
86 "IEEE Transactions on Image Processing, Vol. 2, No 2\n"
87 "Here we use the algorithm on p. 194: Computing of Regional Maxima and Breadth-first Scanning.\n\n"
88 "A marker is derived from the input image INPUT_GRID by subtracting a constant SHIFT_VALUE. Optionally "
89 "the SHIFT_VALUE can be set to zero at the border of the grid (\"Preserve 1px border Yes/No\"). "
90 "OUTPUT_GRID is the difference between the input image and the morphological reconstruction of "
91 "the marker under the input image as mask. If the Option \"Create a binary mask\" is selected, "
92 "the OUTPUT_GRID is thresholded with THRESHOLD, creating a binary image of maxima regions.\n"
93 ));
94
95 CSG_Parameter *pNode;
96
97 Parameters.Add_Grid (NULL, "INPUT_GRID", _TL ("Input Grid"), _TL ("Input for the morphological reconstruction"), PARAMETER_INPUT);
98 Parameters.Add_Grid (NULL, "OBJECT_GRID", _TL("Object Grid"), _TL("Binary object mask"), PARAMETER_OUTPUT, true, SG_DATATYPE_Char);
99 Parameters.Add_Grid (NULL, "DIFFERENCE_GRID", _TL ("Difference Input - Reconstruction"), _TL ("Difference Input - Reconstruction"), PARAMETER_OUTPUT);
100 Parameters.Add_Value (Parameters ("SHIFT"), "SHIFT_VALUE", _TL ("Shift value"), _TL ("Shift value"), PARAMETER_TYPE_Double, 5);
101 Parameters.Add_Value (NULL, "BORDER_YES_NO", _TL ("Preserve 1px border Yes/No"), _TL ("Preserve 1px border Yes/No"), PARAMETER_TYPE_Bool, true);
102 pNode = Parameters.Add_Value (NULL, "BIN_YES_NO", _TL ("Create a binary mask Yes/No"), _TL ("Create a binary mask Yes/No"), PARAMETER_TYPE_Bool, true);
103 Parameters.Add_Value (pNode, "THRESHOLD", _TL ("Threshold"), _TL ("Threshold"), PARAMETER_TYPE_Double, 1);
104 }
105
106 ///////////////////////////////////////////////////////////
107 // //
108 // //
109 // //
110 ///////////////////////////////////////////////////////////
111
On_Execute(void)112 bool CGeomrec::On_Execute(void)
113 {
114 CSG_Grid *pinpgrid, *bingrid, *poutgrid;
115
116
117 unsigned short numrows;
118 unsigned short numcols;
119 double **mask;
120 double **marker;
121
122
123 double h, t;
124 bool pborder, bin;
125
126 pinpgrid = Parameters ("INPUT_GRID")->asGrid ();
127 bingrid = Parameters ("OBJECT_GRID")->asGrid();
128 poutgrid = Parameters ("DIFFERENCE_GRID")->asGrid ();
129
130 h = Parameters ("SHIFT_VALUE")->asDouble ();
131 t = Parameters ("THRESHOLD")->asDouble ();
132 pborder = Parameters ("BORDER_YES_NO")->asBool ();
133 bin = Parameters ("BIN_YES_NO")->asBool ();
134
135 numcols = Get_NY();
136 numrows = Get_NX();
137
138 mask = (double **) matrix_all_alloc (numrows, numcols, 'D', 0);
139 marker = (double **) matrix_all_alloc (numrows, numcols, 'D', 0);
140
141 for (int y = 0; y < Get_NY () && Set_Progress(y, Get_NY()); y++)
142 {
143 #pragma omp parallel for
144 for (int x = 0; x < Get_NX (); x++)
145 {
146 if (pinpgrid->is_NoData(x,y)) // check if there are no_data in input datasets
147 {
148 mask [x][y] = -999999.9;
149 marker [x][y] = -999999.9;
150 }
151 else if ((pborder) && ((x==0)||(y==0)||(x==Get_NX()-1)||(y==Get_NY()-1)))
152 {
153 mask [x][y] = pinpgrid->asDouble(x,y);
154 marker [x][y] = pinpgrid->asDouble(x,y);
155 }
156 else
157 {
158 mask [x][y] = pinpgrid->asDouble(x,y);
159 marker [x][y] = pinpgrid->asDouble(x,y) - h;
160 }
161 }
162 }
163
164 geodesic_morphological_reconstruction (numrows, numcols, mask, marker);
165
166 for (int y = 0; y < Get_NY () && Set_Progress(y, Get_NY()); y++)
167 {
168 #pragma omp parallel for
169 for (int x = 0; x < Get_NX (); x++)
170 {
171 if (pinpgrid->is_NoData(x,y))
172 poutgrid->Set_NoData(x,y);
173 else
174 poutgrid->Set_Value(x,y, mask[x][y]-marker[x][y]);
175 }
176 }
177
178 if (bin)
179 {
180 for (int y = 0; y < Get_NY () && Set_Progress(y, Get_NY()); y++)
181 #pragma omp parallel for
182 for (int x = 0; x < Get_NX (); x++)
183 if ((mask[x][y]-marker[x][y])>t)
184 bingrid->Set_Value(x, y, 1);
185 else
186 bingrid->Set_Value(x, y, 0);
187 }
188
189 matrix_all_free ((void **) mask);
190 matrix_all_free ((void **) marker);
191
192 return( true );
193 }
194
195
196 ///////////////////////////////////////////////////////////
197 // //
198 // //
199 // //
200 ///////////////////////////////////////////////////////////
201
202 //---------------------------------------------------------
203