root / org.gvsig.toolbox / trunk / org.gvsig.toolbox / org.gvsig.toolbox.algorithm / src / main / java / es / unex / sextante / morphometry / curvatures / CurvaturesAlgorithm.java @ 59
History | View | Annotate | Download (11.4 KB)
1 | 59 | nbrodin | package es.unex.sextante.morphometry.curvatures; |
---|---|---|---|
2 | |||
3 | import es.unex.sextante.core.GeoAlgorithm; |
||
4 | import es.unex.sextante.core.AnalysisExtent; |
||
5 | import es.unex.sextante.core.Sextante; |
||
6 | import es.unex.sextante.dataObjects.IRasterLayer; |
||
7 | import es.unex.sextante.exceptions.GeoAlgorithmExecutionException; |
||
8 | import es.unex.sextante.exceptions.RepeatedParameterNameException; |
||
9 | |||
10 | public class CurvaturesAlgorithm |
||
11 | extends
|
||
12 | GeoAlgorithm { |
||
13 | |||
14 | private final static double PLAN_THRESHOLD = 0.00001; |
||
15 | private final static int m_iOffsetX[] = { 0, 1, 1, 1, 0, -1, -1, -1 }; |
||
16 | private final static int m_iOffsetY[] = { 1, 1, 0, -1, -1, -1, 0, 1 }; |
||
17 | |||
18 | public static final String METHOD = "METHOD"; |
||
19 | public static final String DEM = "DEM"; |
||
20 | public static final String HORZ = "HORZ"; |
||
21 | public static final String VERT = "VERT"; |
||
22 | public static final String GLOBAL = "GLOBAL"; |
||
23 | public static final String CLASS = "CLASS"; |
||
24 | |||
25 | IRasterLayer m_DEM = null;
|
||
26 | IRasterLayer m_Curv; |
||
27 | IRasterLayer m_hCurv; |
||
28 | IRasterLayer m_vCurv; |
||
29 | IRasterLayer m_CurvClass; |
||
30 | |||
31 | private double _6DX; |
||
32 | private double _2DX; |
||
33 | private double _4_DX2; |
||
34 | private double _DX2; |
||
35 | |||
36 | |||
37 | @Override
|
||
38 | public boolean processAlgorithm() throws GeoAlgorithmExecutionException { |
||
39 | |||
40 | int x, y;
|
||
41 | int iNX, iNY;
|
||
42 | |||
43 | final int iMethod = m_Parameters.getParameterValueAsInt(METHOD); |
||
44 | m_DEM = m_Parameters.getParameterValueAsRasterLayer(DEM); |
||
45 | |||
46 | m_hCurv = getNewRasterLayer(HORZ, Sextante.getText("Horizontal_curvature"), IRasterLayer.RASTER_DATA_TYPE_FLOAT);
|
||
47 | m_vCurv = getNewRasterLayer(VERT, Sextante.getText("Vertical_curvature"), IRasterLayer.RASTER_DATA_TYPE_FLOAT);
|
||
48 | m_Curv = getNewRasterLayer(GLOBAL, Sextante.getText("Curvature"), IRasterLayer.RASTER_DATA_TYPE_FLOAT);
|
||
49 | m_CurvClass = getNewRasterLayer(CLASS, Sextante.getText("Classification"), IRasterLayer.RASTER_DATA_TYPE_SHORT);
|
||
50 | final AnalysisExtent extent = m_hCurv.getWindowGridExtent();
|
||
51 | |||
52 | m_DEM.setWindowExtent(extent); |
||
53 | |||
54 | iNX = m_DEM.getNX(); |
||
55 | iNY = m_DEM.getNY(); |
||
56 | |||
57 | _2DX = m_DEM.getWindowCellSize() * 2.;
|
||
58 | _DX2 = m_DEM.getWindowCellSize() * m_DEM.getWindowCellSize(); |
||
59 | _4_DX2 = 4 * m_DEM.getWindowCellSize() * m_DEM.getWindowCellSize();
|
||
60 | _6DX = m_DEM.getWindowCellSize() * 6.;
|
||
61 | |||
62 | for (y = 0; (y < iNY) && setProgress(y, iNY); y++) { |
||
63 | for (x = 0; x < iNX; x++) { |
||
64 | switch (iMethod) {
|
||
65 | case 0: |
||
66 | Do_FD_BRM(x, y); |
||
67 | break;
|
||
68 | |||
69 | case 1: |
||
70 | Do_FD_Heerdegen(x, y); |
||
71 | break;
|
||
72 | |||
73 | case 2: |
||
74 | Do_FD_Zevenbergen(x, y); |
||
75 | break;
|
||
76 | |||
77 | case 3: |
||
78 | Do_FD_Haralick(x, y); |
||
79 | break;
|
||
80 | } |
||
81 | } |
||
82 | } |
||
83 | |||
84 | return !m_Task.isCanceled();
|
||
85 | |||
86 | } |
||
87 | |||
88 | |||
89 | @Override
|
||
90 | public void defineCharacteristics() { |
||
91 | |||
92 | final String[] sMethod = { Sextante.getText("Fit_2_Degree_Polynom__Bauer_Rohdenburg_Bork_1985"), |
||
93 | Sextante.getText("Fit_2_Degree_Polynom__Heerdegen_&_Beran_1982"),
|
||
94 | Sextante.getText("Fit_2_Degree_Polynom__Zevenbergen_&_Thorne_1987"),
|
||
95 | Sextante.getText("Fit_3_Degree_Polynom__Haralick_1983") };
|
||
96 | |||
97 | setName(Sextante.getText("Curvatures"));
|
||
98 | setGroup(Sextante.getText("Geomorphometry_and_terrain_analysis"));
|
||
99 | setUserCanDefineAnalysisExtent(true);
|
||
100 | |||
101 | try {
|
||
102 | m_Parameters.addInputRasterLayer(DEM, Sextante.getText("Elevation"), true); |
||
103 | m_Parameters.addSelection(METHOD, Sextante.getText("Method"), sMethod);
|
||
104 | addOutputRasterLayer(HORZ, Sextante.getText("Horizontal_curvature"));
|
||
105 | addOutputRasterLayer(VERT, Sextante.getText("Vertical_curvature"));
|
||
106 | addOutputRasterLayer(GLOBAL, Sextante.getText("Curvature"));
|
||
107 | addOutputRasterLayer(CLASS, Sextante.getText("Classification"));
|
||
108 | } |
||
109 | catch (final RepeatedParameterNameException e) { |
||
110 | Sextante.addErrorToLog(e); |
||
111 | } |
||
112 | |||
113 | } |
||
114 | |||
115 | |||
116 | private void Do_FD_BRM(final int x, |
||
117 | final int y) { |
||
118 | |||
119 | double zm[], D, E, F, G, H; |
||
120 | |||
121 | zm = new double[9]; |
||
122 | |||
123 | if (Get_SubMatrix3x3(x, y, zm)) {
|
||
124 | D = ((zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8]) - 2 * (zm[1] + zm[4] + zm[7])) / _DX2; |
||
125 | E = ((zm[0] + zm[6] + zm[1] + zm[7] + zm[2] + zm[8]) - 2 * (zm[3] + zm[4] + zm[5])) / _DX2; |
||
126 | F = (zm[8] + zm[0] - zm[7]) / _4_DX2; |
||
127 | G = ((zm[2] - zm[0]) + (zm[5] - zm[3]) + (zm[8] - zm[6])) / _6DX; |
||
128 | H = ((zm[6] - zm[0]) + (zm[7] - zm[1]) + (zm[8] - zm[2])) / _6DX; |
||
129 | |||
130 | Set_Parameters_Derive(x, y, D, E, F, G, H); |
||
131 | } |
||
132 | } |
||
133 | |||
134 | |||
135 | private void Do_FD_Heerdegen(final int x, |
||
136 | final int y) { |
||
137 | |||
138 | double zm[], a, b, D, E, F, G, H; |
||
139 | |||
140 | zm = new double[9]; |
||
141 | |||
142 | if (Get_SubMatrix3x3(x, y, zm)) {
|
||
143 | a = zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8]; |
||
144 | b = zm[0] + zm[1] + zm[2] + zm[6] + zm[7] + zm[8]; |
||
145 | D = (0.3 * a - 0.2 * b) / _DX2; |
||
146 | E = (0.3 * b - 0.2 * a) / _DX2; |
||
147 | F = (zm[0] - zm[2] - zm[6] + zm[8]) / _4_DX2; |
||
148 | G = (-zm[0] + zm[2] - zm[3] + zm[5] - zm[6] + zm[8]) / _6DX; |
||
149 | H = (-zm[0] - zm[1] - zm[2] + zm[6] + zm[7] + zm[8]) / _6DX; |
||
150 | |||
151 | Set_Parameters_Derive(x, y, D, E, F, G, H); |
||
152 | } |
||
153 | } |
||
154 | |||
155 | |||
156 | private void Do_FD_Zevenbergen(final int x, |
||
157 | final int y) { |
||
158 | |||
159 | double zm[], D, E, F, G, H; |
||
160 | |||
161 | zm = new double[9]; |
||
162 | |||
163 | if (Get_SubMatrix3x3(x, y, zm)) {
|
||
164 | |||
165 | D = ((zm[3] + zm[5]) / 2.0 - zm[4]) / _DX2; |
||
166 | E = ((zm[1] + zm[7]) / 2.0 - zm[4]) / _DX2; |
||
167 | F = (zm[0] - zm[2] - zm[6] + zm[8]) / _4_DX2; |
||
168 | G = (zm[5] - zm[3]) / _2DX; |
||
169 | H = (zm[7] - zm[1]) / _2DX; |
||
170 | |||
171 | Set_Parameters_Derive(x, y, D, E, F, G, H); |
||
172 | } |
||
173 | } |
||
174 | |||
175 | |||
176 | private void Do_FD_Haralick(final int x, |
||
177 | final int y) { |
||
178 | |||
179 | final int Mtrx[][][] = { |
||
180 | { { 31, -5, -17, -5, 31 }, { -44, -62, -68, -62, -44 }, { 0, 0, 0, 0, 0 }, { 44, 62, 68, 62, 44 }, |
||
181 | { -31, 5, 17, 5, -31 } }, |
||
182 | { { 31, -44, 0, 44, -31 }, { -5, -62, 0, 62, 5 }, { -17, -68, 0, 68, 17 }, { -5, -62, 0, 62, 5 }, |
||
183 | { 31, -44, 0, 44, -31 } }, |
||
184 | { { 2, 2, 2, 2, 2 }, { -1, -1, -1, -1, -1 }, { -2, -2, -2, -2, -2 }, { -1, -1, -1, -1, -1 }, { 2, 2, 2, 2, 2 } }, |
||
185 | { { 4, 2, 0, -2, -4 }, { 2, 1, 0, -1, -2 }, { 0, 0, 0, 0, 0 }, { -2, -1, 0, 1, 2 }, { -4, -2, 0, 2, 4 } }, |
||
186 | { { 2, -1, -2, -1, 2 }, { 2, -1, -2, -1, 2 }, { 2, -1, -2, -1, 2 }, { 2, -1, -2, -1, 2 }, { 2, -1, -2, -1, 2 } } }; |
||
187 | |||
188 | final int QMtrx[] = { 4200, 4200, 700, 1000, 700 }; |
||
189 | |||
190 | int i, ix, iy, n;
|
||
191 | double Sum, zm[], k[]; |
||
192 | |||
193 | zm = new double[25]; |
||
194 | k = new double[5]; |
||
195 | |||
196 | if (Get_SubMatrix5x5(x, y, zm)) {
|
||
197 | for (i = 0; i < 5; i++) { |
||
198 | for (n = 0, Sum = 0.0, iy = 0; iy < 5; iy++) { |
||
199 | for (ix = 0; ix < 5; ix++, n++) { |
||
200 | Sum += zm[n] * Mtrx[i][ix][iy]; |
||
201 | } |
||
202 | } |
||
203 | |||
204 | k[i] = Sum / QMtrx[i]; |
||
205 | } |
||
206 | |||
207 | Set_Parameters_Derive(x, y, k[4], k[2], k[3], k[1], k[0]); |
||
208 | } |
||
209 | } |
||
210 | |||
211 | |||
212 | // additional methods
|
||
213 | |||
214 | private void Set_Parameters(final int x, |
||
215 | final int y, |
||
216 | final double dCurv, |
||
217 | final double dHCurv, |
||
218 | final double dVCurv) { |
||
219 | |||
220 | int iClass;
|
||
221 | |||
222 | m_Curv.setCellValue(x, y, dCurv); |
||
223 | m_hCurv.setCellValue(x, y, dHCurv); |
||
224 | m_vCurv.setCellValue(x, y, dVCurv); |
||
225 | |||
226 | |||
227 | iClass = dHCurv < -PLAN_THRESHOLD ? 0 : (dHCurv <= PLAN_THRESHOLD ? 3 : 6); |
||
228 | iClass += dVCurv < -PLAN_THRESHOLD ? 0 : (dVCurv <= PLAN_THRESHOLD ? 1 : 2); |
||
229 | |||
230 | m_CurvClass.setCellValue(x, y, iClass); |
||
231 | |||
232 | } |
||
233 | |||
234 | |||
235 | private void Set_Parameters_Derive(final int x, |
||
236 | final int y, |
||
237 | final double D, |
||
238 | final double E, |
||
239 | final double F, |
||
240 | final double G, |
||
241 | final double H) { |
||
242 | |||
243 | double k1, k2, Curv, vCurv, hCurv;
|
||
244 | |||
245 | k1 = F * G * H; |
||
246 | k2 = G * G + H * H; |
||
247 | |||
248 | if (k2 != 0) { |
||
249 | Curv = -2.0 * (E + D);
|
||
250 | vCurv = -2.0 * (D * G * G + E * H * H + k1) / k2;
|
||
251 | hCurv = -2.0 * (D * H * H + E * G * G - k1) / k2;
|
||
252 | |||
253 | } |
||
254 | else {
|
||
255 | Curv = vCurv = hCurv = 0.0;
|
||
256 | } |
||
257 | |||
258 | Set_Parameters(x, y, Curv, hCurv, vCurv); |
||
259 | } |
||
260 | |||
261 | |||
262 | private void Set_Parameters_NoData(final int x, |
||
263 | final int y) { |
||
264 | |||
265 | m_Curv.setNoData(x, y); |
||
266 | m_hCurv.setNoData(x, y); |
||
267 | m_vCurv.setNoData(x, y); |
||
268 | |||
269 | } |
||
270 | |||
271 | |||
272 | // ---------------------------------------------------------
|
||
273 | // Indexing of the Submatrix:
|
||
274 | //
|
||
275 | // +-------+ +-------+
|
||
276 | // | 7 0 1 | | 2 5 8 |
|
||
277 | // | 6 * 2 | => | 1 4 7 |
|
||
278 | // | 5 4 3 | | 0 3 6 |
|
||
279 | // +-------+ +-------+
|
||
280 | //
|
||
281 | // ---------------------------------------------------------
|
||
282 | private boolean Get_SubMatrix3x3(final int x, |
||
283 | final int y, |
||
284 | final double SubMatrix[]) { |
||
285 | |||
286 | final int iSub[] = { 5, 8, 7, 6, 3, 0, 1, 2 }; |
||
287 | |||
288 | int i;
|
||
289 | double z, z2;
|
||
290 | |||
291 | z = m_DEM.getCellValueAsDouble(x, y); |
||
292 | |||
293 | if (m_DEM.isNoDataValue(z)) {
|
||
294 | Set_Parameters_NoData(x, y); |
||
295 | } |
||
296 | else {
|
||
297 | SubMatrix[4] = 0.0; |
||
298 | for (i = 0; i < 8; i++) { |
||
299 | z2 = m_DEM.getCellValueAsDouble(x + m_iOffsetX[i], y + m_iOffsetY[i]); |
||
300 | if (!m_DEM.isNoDataValue(z2)) {
|
||
301 | SubMatrix[iSub[i]] = z2 - z; |
||
302 | } |
||
303 | else {
|
||
304 | z2 = m_DEM.getCellValueAsDouble(x + m_iOffsetY[(i + 4) % 8], y + m_iOffsetY[(i + 4) % 8]); |
||
305 | if (!m_DEM.isNoDataValue(z2)) {
|
||
306 | SubMatrix[iSub[i]] = z - z2; |
||
307 | } |
||
308 | else {
|
||
309 | SubMatrix[iSub[i]] = 0.0;
|
||
310 | } |
||
311 | } |
||
312 | } |
||
313 | |||
314 | return (true); |
||
315 | } |
||
316 | |||
317 | return (false); |
||
318 | } |
||
319 | |||
320 | |||
321 | private boolean Get_SubMatrix5x5(final int x, |
||
322 | final int y, |
||
323 | final double SubMatrix[]) { |
||
324 | int i, ix, iy, jx, jy;
|
||
325 | double z, z2;
|
||
326 | |||
327 | z = m_DEM.getCellValueAsDouble(x, y); |
||
328 | |||
329 | if (!m_DEM.isNoDataValue(z)) {
|
||
330 | for (i = 0, iy = y - 2; iy <= y + 2; iy++) { |
||
331 | jy = iy < 0 ? 0 : (iy >= m_DEM.getNY() ? m_DEM.getNY() - 1 : iy); |
||
332 | for (ix = x - 2; ix <= x + 2; ix++, i++) { |
||
333 | jx = ix < 0 ? 0 : (ix >= m_DEM.getNX() ? m_DEM.getNY() - 1 : ix); |
||
334 | z2 = m_DEM.getCellValueAsDouble(jx, jy); |
||
335 | if (!m_DEM.isNoDataValue(z2)) {
|
||
336 | SubMatrix[i] = z2 - z; |
||
337 | } |
||
338 | else {
|
||
339 | SubMatrix[i] = 0.0;
|
||
340 | } |
||
341 | } |
||
342 | } |
||
343 | |||
344 | return (true); |
||
345 | } |
||
346 | |||
347 | return (false); |
||
348 | } |
||
349 | |||
350 | |||
351 | } |