001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, Geotools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation; either
010: * version 2.1 of the License, or (at your option) any later version.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: */
017: package org.geotools.coverage.processing.operation;
018:
019: // J2SE, JAI and extensions
020: import java.awt.Color;
021: import java.awt.geom.AffineTransform;
022: import javax.media.jai.KernelJAI;
023: import javax.media.jai.ParameterList;
024: import javax.media.jai.ParameterListDescriptor;
025: import javax.media.jai.ParameterListDescriptorImpl;
026: import javax.media.jai.operator.GradientMagnitudeDescriptor; // For javadoc
027: import javax.units.Unit;
028:
029: // OpenGIS dependencies
030: import org.opengis.util.InternationalString;
031: import org.opengis.referencing.cs.CoordinateSystem;
032: import org.opengis.referencing.operation.MathTransform1D;
033: import org.opengis.referencing.operation.MathTransform2D;
034:
035: // Geotools dependencies
036: import org.geotools.coverage.Category;
037: import org.geotools.coverage.grid.GridGeometry2D;
038: import org.geotools.coverage.grid.GridCoverage2D;
039: import org.geotools.coverage.processing.OperationJAI;
040: import org.geotools.referencing.operation.matrix.XAffineTransform;
041: import org.geotools.util.NumberRange;
042:
043: /**
044: * Edge detector which computes the magnitude of the image gradient vector in two orthogonal
045: * directions. The result of the {@code "GradientMagnitude"} operation may be defined as:
046: * <p>
047: * <BLOCKQUOTE><CODE>
048: * dst[<var>x</var>][<var>y</var>][<var>b</var>] = {@linkplain Math#sqrt sqrt}(
049: * <strong>SH</strong>(<var>x</var>,<var>y</var>,<var>b</var>)<sup>2</sup> +
050: * <strong>SV</strong>(<var>x</var>,<var>y</var>,<var>b</var>)<sup>2</sup>)
051: * × <var>scale</var>
052: * </CODE></BLOCKQUOTE>
053: * <p>
054: * where {@code SH(x,y,b)} and {@code SV(x,y,b)} are the horizontal and vertical gradient images
055: * generated from band <var>b</var> of the source image by correlating it with the supplied
056: * orthogonal (horizontal and vertical) gradient masks.
057: * <p>
058: * Before to compute the gradients, the kernels are tested against artificials horizontal and
059: * vertical gradients of one <cite>unit of sample / unit of localization</cite>. For example:
060: * <p>
061: * <UL>
062: * <LI>For an image of elevation (meters) using a geographic coordinate system (degrees
063: * of latitude and longitude), the units are <strong>meters/degree</strong>.</LI>
064: * <LI>For an image of temperature (°C) using a projected coordinate system (kilometers),
065: * the units are <strong>°C/km</strong>.</LI>
066: * </UL>
067: * <p>
068: * Kernels are normalized by dividing all their coefficients by the result of this test. In other
069: * words, kernels are normalized in such a way that applying the {@code "GradientMagnitude"}
070: * operation on a horizontal or vertical gradient of 1 such "geophysical" units will give a result
071: * of 1. This is an attempt to give geophysical meaning to the numbers produced by the
072: * {@code "GradientMagnitude"} operation. This normalization depends of the coverage's
073: * {@linkplain GridCoverage2D#getGridGeometry grid geometry}.
074: * <p>
075: * <strong>NOTE:</strong> When the masks are symetric (e.g. Sobel, Prewitt (or Smoothed),
076: * isotropic, <cite>etc.</cite>), then the above-cited algorithm produces the same result
077: * than the "<cite>normalization factor</cite>" × "<cite>spatial factor</cite>" published by:
078: *
079: * <blockquote>
080: * Simpson, J.J. (1990), "<u>On the accurate detection and enhancement of oceanic features
081: * observed in satellite data</u>" <i>in</i> Remote sensing environment, <b>33:17-33</b>.
082: * </blockquote>
083: *
084: * However, for non-symetric masks (e.g. Kirsch), then a difference is found.
085: *
086: * <P><STRONG>Name:</STRONG> <CODE>"GradientMagnitude"</CODE><BR>
087: * <STRONG>JAI operator:</STRONG> <CODE>"{@linkplain GradientMagnitudeDescriptor GradientMagnitude}"</CODE><BR>
088: * <STRONG>Parameters:</STRONG></P>
089: * <table border='3' cellpadding='6' bgcolor='F4F8FF'>
090: * <tr bgcolor='#B9DCFF'>
091: * <th>Name</th>
092: * <th>Class</th>
093: * <th>Default value</th>
094: * <th>Minimum value</th>
095: * <th>Maximum value</th>
096: * </tr>
097: * <tr>
098: * <td>{@code "Source"}</td>
099: * <td>{@link org.geotools.coverage.grid.GridCoverage2D}</td>
100: * <td align="center">N/A</td>
101: * <td align="center">N/A</td>
102: * <td align="center">N/A</td>
103: * </tr>
104: * <tr>
105: * <td>{@code "Mask1"}</td>
106: * <td>{@link KernelJAI}</td>
107: * <td>{@link KernelJAI#GRADIENT_MASK_SOBEL_HORIZONTAL}</td>
108: * <td align="center">N/A</td>
109: * <td align="center">N/A</td>
110: * </tr>
111: * <tr>
112: * <td>{@code "Mask2"}</td>
113: * <td>{@link KernelJAI}</td>
114: * <td>{@link KernelJAI#GRADIENT_MASK_SOBEL_VERTICAL}</td>
115: * <td align="center">N/A</td>
116: * <td align="center">N/A</td>
117: * </tr>
118: * </table>
119: *
120: * @since 2.2
121: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/coverage/src/main/java/org/geotools/coverage/processing/operation/GradientMagnitude.java $
122: * @version $Id: GradientMagnitude.java 23157 2006-12-01 01:29:53Z desruisseaux $
123: * @author Martin Desruisseaux
124: *
125: * @see org.geotools.coverage.processing.Operations#gradientMagnitude
126: * @see GradientMagnitudeDescriptor
127: */
128: public class GradientMagnitude extends OperationJAI {
129: /**
130: * Serial number for interoperability with different versions.
131: */
132: private static final long serialVersionUID = -1514713427236924048L;
133:
134: /**
135: * Set to {@code true} for enabling some tracing code.
136: */
137: private static final boolean DEBUG = false;
138:
139: /**
140: * Set to {@code true} to enable automatic kernel normalization. Normalization modifies
141: * kernel coefficients according the "grid to coordinate system" transform in order to get
142: * some meaningful engineering units (e.g. °C/km). The normalization factor is computed by
143: * testing the original kernels against synthetic horizontal and vertical gradients of
144: * 1 sampleUnit/csUnit.
145: */
146: private static final boolean NORMALIZE = true;
147:
148: /**
149: * The default scale factor to apply on the range computed by {@link #deriveCategory}. For
150: * example a value of 0.25 means that only values from 0 to 25% of the expected maximum will
151: * appears in different colors.
152: */
153: private static final double DEFAULT_RANGE_SCALE = 1.0;
154:
155: /**
156: * The default color palette for the gradients.
157: */
158: private static final Color[] DEFAULT_COLOR_PALETTE = new Color[] {
159: new Color(16, 32, 64), new Color(192, 224, 255) };
160:
161: /**
162: * A flag indicating that {@link #getNormalizationFactorSquared}
163: * should tests the horizontal gradient computed by the supplied kernel.
164: */
165: private static final int HORIZONTAL = 1;
166:
167: /**
168: * A flag indicating that {@link #getNormalizationFactorSquared}
169: * should tests the vertical gradient computed by the supplied kernel.
170: */
171: private static final int VERTICAL = 2;
172:
173: /**
174: * Constructs a default gradient magnitude operation.
175: */
176: public GradientMagnitude() {
177: super ("GradientMagnitude");
178: }
179:
180: /**
181: * Returns a scale factor for the supplied kernel. If {@code kernel} computes horizontal
182: * grandient, this method returns {@code scaleX}. Otherwise, if {@code kernel} computes
183: * vertical gradient, then this method returns {@code scaleY}. Otherwise, returns a geometric
184: * combinaison of both.
185: */
186: private static double getScaleFactor(final KernelJAI kernel,
187: double scaleX, double scaleY) {
188: scaleX *= scaleX;
189: scaleY *= scaleY;
190: double factorX = getNormalizationFactorSquared(kernel,
191: HORIZONTAL);
192: double factorY = getNormalizationFactorSquared(kernel, VERTICAL);
193: double factor2 = (factorX * scaleX + factorY * scaleY)
194: / (factorX + factorY);
195: return Math.sqrt(factor2);
196: }
197:
198: /**
199: * Returns the square of a normalization factor for the supplied kernel.
200: * The kernel can be normalized by invoking {@link #divide(KernelJAI,double)}
201: * with the square root of this value.
202: *
203: * @param kernel The kernel for which to compute normalization factor.
204: * @param type Any combinaison of {@link #HORIZONTAL} and {@link #VERTICAL}.
205: * @return The square of a normalization factor that could be applied on the kernel.
206: */
207: private static double getNormalizationFactorSquared(
208: final KernelJAI kernel, final int type) {
209: double sumH = 0;
210: double sumV = 0;
211: final int width = kernel.getWidth();
212: final int height = kernel.getHeight();
213: /*
214: * Tests the kernel with a horizontal gradient [ -1 0 1 ]
215: * of 1/pixel. For example, we get sumH=8 with [ -2 0 2 ]
216: * the horizontal Sobel kernel show on right: [ -1 0 1 ]
217: */
218: if ((type & HORIZONTAL) != 0) {
219: int value = kernel.getYOrigin();
220: for (int y = height; --y >= 0;) {
221: for (int x = width; --x >= 0;) {
222: sumH += value * kernel.getElement(x, y);
223: }
224: value--;
225: }
226: }
227: /*
228: * Tests the kernel with a vertical gradient of [ -1 -2 -1 ]
229: * 1/pixel. For example, we get sumV=8 with the [ 0 0 0 ]
230: * vertical Sobel kernel show on right: [ 1 2 1 ]
231: */
232: if ((type & VERTICAL) != 0) {
233: int value = kernel.getXOrigin();
234: for (int x = width; --x >= 0;) {
235: for (int y = height; --y >= 0;) {
236: sumV += value * kernel.getElement(x, y);
237: }
238: value--;
239: }
240: }
241: return (sumH * sumH) + (sumV * sumV);
242: }
243:
244: /**
245: * Returns the normalization factor for the supplied kernel. The kernel can
246: * be normalized by invoking {@link #divide(KernelJAI,double)} with this factor.
247: *
248: * @param mask1 The first kernel for which to compute a normalization factor.
249: * @param mask2 The second kernel for which to compute a normalization factor.
250: * @return The normalization factor that could be applied on both kernels.
251: *
252: * @todo When the masks are symetric (e.g. Sobel, Prewitt (or Smoothed), isotropic, etc.),
253: * then this algorithm matches the "normalization factor" times "spatial factor"
254: * provided by
255: *
256: * J.J. Simpson (1990), "On the accurate detection and enhancement of oceanic features
257: * observed in satellite data" in Remote sensing environment, 33:17-33.
258: *
259: * However, for non-symetric masks (e.g. Kirsch), then a difference is found.
260: * We should provides a way to disable normalization when the user did it himself
261: * according some other rules than the one used here.
262: */
263: private static double getNormalizationFactor(final KernelJAI mask1,
264: final KernelJAI mask2) {
265: double factor;
266: factor = getNormalizationFactorSquared(mask1, HORIZONTAL
267: | VERTICAL);
268: factor += getNormalizationFactorSquared(mask2, HORIZONTAL
269: | VERTICAL);
270: factor = Math.sqrt(factor / 2);
271: return factor;
272: }
273:
274: /**
275: * Divides a kernel by some number.
276: *
277: * @param kernel The kernel to divide.
278: * @param denominator The factor to divide by.
279: * @return The resulting kernel.
280: */
281: private static KernelJAI divide(KernelJAI kernel,
282: final double denominator) {
283: if (denominator != 1) {
284: final float[] data = kernel.getKernelData();
285: final int length = data.length;
286: for (int i = 0; i < length; i++) {
287: data[i] /= denominator;
288: }
289: kernel = new KernelJAI(kernel.getWidth(), kernel
290: .getHeight(), kernel.getXOrigin(), kernel
291: .getYOrigin(), data);
292: }
293: return kernel;
294: }
295:
296: /**
297: * Applies the operation on grid coverage. The default implementation looks for kernels
298: * specified in the {@link org.geotools.coverage.processing.OperationJAI.Parameters#parameters
299: * parameter block} and divide them by the distance between pixels, in the grid coverage's
300: * coordinate reference system.
301: */
302: protected GridCoverage2D deriveGridCoverage(
303: final GridCoverage2D[] sources, final Parameters parameters) {
304: if (NORMALIZE) {
305: final ParameterList block = parameters.parameters;
306: KernelJAI mask1 = (KernelJAI) block
307: .getObjectParameter("Mask1");
308: KernelJAI mask2 = (KernelJAI) block
309: .getObjectParameter("Mask2");
310: /*
311: * Normalizes the kernel in such a way that pixel values likes
312: * [-2 -1 0 +1 +2] will give a gradient of about 1 unit/pixel.
313: */
314: double factor = getNormalizationFactor(mask1, mask2);
315: if (!(factor > 0)) {
316: // Do not transform if factor is 0 or NaN.
317: factor = 1;
318: }
319: /*
320: * Computes a scale factor taking in account the transformation from
321: * grid to coordinate system. This scale will convert gradient from
322: * 1 unit/pixel to 1 unit/meters or 1 unit/degrees, depending the
323: * coordinate reference systems axis unit.
324: */
325: double scaleMask1 = 1;
326: double scaleMask2 = 1;
327: if (sources.length != 0) {
328: final MathTransform2D mtr;
329: mtr = ((GridGeometry2D) sources[0].getGridGeometry())
330: .getGridToCRS2D();
331: if (mtr instanceof AffineTransform) {
332: final AffineTransform tr = (AffineTransform) mtr;
333: final double scaleX = XAffineTransform
334: .getScaleX0(tr);
335: final double scaleY = XAffineTransform
336: .getScaleY0(tr);
337: scaleMask1 = getScaleFactor(mask1, scaleX, scaleY);
338: scaleMask2 = getScaleFactor(mask2, scaleX, scaleY);
339: if (!(scaleMask1 > 0 && scaleMask2 > 0)) {
340: // Do not rescale if scale is 0 or NaN.
341: scaleMask1 = 1;
342: scaleMask2 = 1;
343: }
344: if (DEBUG) {
345: System.out.print("factor= ");
346: System.out.println(factor);
347: System.out.print("scaleMask1= ");
348: System.out.println(scaleMask1);
349: System.out.print("scaleMask1= ");
350: System.out.println(scaleMask2);
351: }
352: }
353: }
354: block.setParameter("Mask1", divide(mask1, factor
355: * scaleMask1));
356: block.setParameter("Mask2", divide(mask2, factor
357: * scaleMask2));
358: }
359: return super .deriveGridCoverage(sources, parameters);
360: }
361:
362: /**
363: * Derives the quantitative category for a band in the destination image.
364: * This implementation computes the expected gradient range from the two
365: * masks and the value range in the source grid coverage.
366: */
367: protected Category deriveCategory(final Category[] categories,
368: final Parameters parameters) {
369: NumberRange range = null;
370: Category category = categories[0];
371: final NumberRange samples = category.geophysics(false)
372: .getRange();
373: final boolean isGeophysics = (category == category
374: .geophysics(true));
375: /*
376: * Computes a default range of output values one from the normalized kernels.
377: * The normalization has been done by 'deriveGridCoverage' before this method
378: * is invoked. The algorithm is as below:
379: *
380: * - Computes the value produced by the kernels for an artificial gradient of 1 unit/pixel.
381: * - Transforms into a lower gradient of 1 unit/(kernel size).
382: * - Transforms into a gradient of (maximal range)/(kernel size).
383: * - Applies an arbitrary correction factor for more convenient range in most cases.
384: */
385: final ParameterList block = parameters.parameters;
386: final KernelJAI mask1 = (KernelJAI) block
387: .getObjectParameter("Mask1");
388: final KernelJAI mask2 = (KernelJAI) block
389: .getObjectParameter("Mask2");
390: final double size = (mask1.getWidth() + mask1.getHeight()
391: + mask2.getWidth() + mask2.getHeight()) / 4.0;
392: double factor = getNormalizationFactor(mask1, mask2)
393: / (size - 1);
394: if (factor > 0 && !Double.isInfinite(factor)) {
395: range = category.geophysics(true).getRange();
396: final double minimum = range.getMinimum();
397: final double maximum = range.getMaximum();
398: factor *= (maximum - minimum) * DEFAULT_RANGE_SCALE;
399: range = new NumberRange(0, factor);
400: }
401: if (range != null) {
402: category = new Category(category.getName(),
403: DEFAULT_COLOR_PALETTE, samples, range);
404: return category.geophysics(isGeophysics);
405: }
406: return super .deriveCategory(categories, parameters);
407: }
408:
409: /**
410: * Derives the unit of data for a band in the destination image.
411: * This method compute {@code sample}/{@code axis} where:
412: *
413: * <ul>
414: * <li>{@code sample} is the sample unit in source image.</li>
415: * <li>{@code axis} is the coordinate reference system axis unit.</li>
416: * </ul>
417: */
418: protected Unit deriveUnit(final Unit[] units,
419: final Parameters parameters) {
420: final CoordinateSystem cs = parameters.crs
421: .getCoordinateSystem();
422: if (!DEBUG && units.length == 1 && units[0] != null) {
423: final Unit spatialUnit = cs.getAxis(0).getUnit();
424: for (int i = Math.min(cs.getDimension(), 2); --i >= 0;) {
425: if (!spatialUnit.equals(cs.getAxis(i).getUnit())) {
426: return super .deriveUnit(units, parameters);
427: }
428: }
429: try {
430: return units[0].divide(spatialUnit);
431: } catch (RuntimeException exception) {
432: // Can't compute units... We will compute image data
433: // anyway, but the result will have no know unit.
434: // TODO: Catch a more specific exception.
435: }
436: }
437: return super.deriveUnit(units, parameters);
438: }
439: }
|