/*
* GeoTools - The Open Source Java GIS Toolkit
* http://geotools.org
*
* (C) 2001-2008, Open Source Geospatial Foundation (OSGeo)
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation;
* version 2.1 of the License.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*/
package org.geotools.coverage.processing.operation;
// J2SE, JAI and extensions
import java.awt.Color;
import java.awt.geom.AffineTransform;
import javax.media.jai.KernelJAI;
import javax.media.jai.ParameterList;
import javax.media.jai.ParameterListDescriptor;
import javax.media.jai.ParameterListDescriptorImpl;
import javax.media.jai.operator.GradientMagnitudeDescriptor; // For javadoc
import javax.measure.unit.Unit;
// OpenGIS dependencies
import org.opengis.util.InternationalString;
import org.opengis.referencing.cs.CoordinateSystem;
import org.opengis.referencing.operation.MathTransform1D;
import org.opengis.referencing.operation.MathTransform2D;
// Geotools dependencies
import org.geotools.coverage.Category;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.processing.OperationJAI;
import org.geotools.referencing.operation.matrix.XAffineTransform;
import org.geotools.util.NumberRange;
/**
* Edge detector which computes the magnitude of the image gradient vector in two orthogonal
* directions. The result of the {@code "GradientMagnitude"} operation may be defined as:
* <p>
* <BLOCKQUOTE><CODE>
* dst[<var>x</var>][<var>y</var>][<var>b</var>] = {@linkplain Math#sqrt sqrt}(
* <strong>SH</strong>(<var>x</var>,<var>y</var>,<var>b</var>)<sup>2</sup> +
* <strong>SV</strong>(<var>x</var>,<var>y</var>,<var>b</var>)<sup>2</sup>)
* × <var>scale</var>
* </CODE></BLOCKQUOTE>
* <p>
* where {@code SH(x,y,b)} and {@code SV(x,y,b)} are the horizontal and vertical gradient images
* generated from band <var>b</var> of the source image by correlating it with the supplied
* orthogonal (horizontal and vertical) gradient masks.
* <p>
* Before to compute the gradients, the kernels are tested against artificials horizontal and
* vertical gradients of one <cite>unit of sample / unit of localization</cite>. For example:
* <p>
* <UL>
* <LI>For an image of elevation (meters) using a geographic coordinate system (degrees
* of latitude and longitude), the units are <strong>meters/degree</strong>.</LI>
* <LI>For an image of temperature (°C) using a projected coordinate system (kilometers),
* the units are <strong>°C/km</strong>.</LI>
* </UL>
* <p>
* Kernels are normalized by dividing all their coefficients by the result of this test. In other
* words, kernels are normalized in such a way that applying the {@code "GradientMagnitude"}
* operation on a horizontal or vertical gradient of 1 such "geophysical" units will give a result
* of 1. This is an attempt to give geophysical meaning to the numbers produced by the
* {@code "GradientMagnitude"} operation. This normalization depends of the coverage's
* {@linkplain GridCoverage2D#getGridGeometry grid geometry}.
* <p>
* <strong>NOTE:</strong> When the masks are symetric (e.g. Sobel, Prewitt (or Smoothed),
* isotropic, <cite>etc.</cite>), then the above-cited algorithm produces the same result
* than the "<cite>normalization factor</cite>" × "<cite>spatial factor</cite>" published by:
*
* <blockquote>
* Simpson, J.J. (1990), "<u>On the accurate detection and enhancement of oceanic features
* observed in satellite data</u>" <i>in</i> Remote sensing environment, <b>33:17-33</b>.
* </blockquote>
*
* However, for non-symetric masks (e.g. Kirsch), then a difference is found.
*
* <P><STRONG>Name:</STRONG> <CODE>"GradientMagnitude"</CODE><BR>
* <STRONG>JAI operator:</STRONG> <CODE>"{@linkplain GradientMagnitudeDescriptor GradientMagnitude}"</CODE><BR>
* <STRONG>Parameters:</STRONG></P>
* <table border='3' cellpadding='6' bgcolor='F4F8FF'>
* <tr bgcolor='#B9DCFF'>
* <th>Name</th>
* <th>Class</th>
* <th>Default value</th>
* <th>Minimum value</th>
* <th>Maximum value</th>
* </tr>
* <tr>
* <td>{@code "Source"}</td>
* <td>{@link org.geotools.coverage.grid.GridCoverage2D}</td>
* <td align="center">N/A</td>
* <td align="center">N/A</td>
* <td align="center">N/A</td>
* </tr>
* <tr>
* <td>{@code "Mask1"}</td>
* <td>{@link KernelJAI}</td>
* <td>{@link KernelJAI#GRADIENT_MASK_SOBEL_HORIZONTAL}</td>
* <td align="center">N/A</td>
* <td align="center">N/A</td>
* </tr>
* <tr>
* <td>{@code "Mask2"}</td>
* <td>{@link KernelJAI}</td>
* <td>{@link KernelJAI#GRADIENT_MASK_SOBEL_VERTICAL}</td>
* <td align="center">N/A</td>
* <td align="center">N/A</td>
* </tr>
* </table>
*
* @since 2.2
*
* @source $URL$
* @version $Id$
* @author Martin Desruisseaux (IRD)
*
* @see org.geotools.coverage.processing.Operations#gradientMagnitude
* @see GradientMagnitudeDescriptor
*/
public class GradientMagnitude extends OperationJAI {
/**
* Serial number for interoperability with different versions.
*/
private static final long serialVersionUID = -1514713427236924048L;
/**
* Set to {@code true} for enabling some tracing code.
*/
private static final boolean DEBUG = false;
/**
* Set to {@code true} to enable automatic kernel normalization. Normalization modifies
* kernel coefficients according the "grid to coordinate system" transform in order to get
* some meaningful engineering units (e.g. °C/km). The normalization factor is computed by
* testing the original kernels against synthetic horizontal and vertical gradients of
* 1 sampleUnit/csUnit.
*/
private static final boolean NORMALIZE = true;
/**
* The default scale factor to apply on the range computed by {@link #deriveCategory}. For
* example a value of 0.25 means that only values from 0 to 25% of the expected maximum will
* appears in different colors.
*/
private static final double DEFAULT_RANGE_SCALE = 1.0;
/**
* The default color palette for the gradients.
*/
private static final Color[] DEFAULT_COLOR_PALETTE = new Color[] {
new Color( 16, 32, 64),
new Color(192,224,255)
};
/**
* A flag indicating that {@link #getNormalizationFactorSquared}
* should tests the horizontal gradient computed by the supplied kernel.
*/
private static final int HORIZONTAL = 1;
/**
* A flag indicating that {@link #getNormalizationFactorSquared}
* should tests the vertical gradient computed by the supplied kernel.
*/
private static final int VERTICAL = 2;
/**
* Constructs a default gradient magnitude operation.
*/
public GradientMagnitude() {
super("GradientMagnitude");
}
/**
* Returns a scale factor for the supplied kernel. If {@code kernel} computes horizontal
* grandient, this method returns {@code scaleX}. Otherwise, if {@code kernel} computes
* vertical gradient, then this method returns {@code scaleY}. Otherwise, returns a geometric
* combinaison of both.
*/
private static double getScaleFactor(final KernelJAI kernel, double scaleX, double scaleY) {
scaleX *= scaleX;
scaleY *= scaleY;
double factorX = getNormalizationFactorSquared(kernel, HORIZONTAL);
double factorY = getNormalizationFactorSquared(kernel, VERTICAL);
double factor2 = (factorX*scaleX + factorY*scaleY) / (factorX+factorY);
return Math.sqrt(factor2);
}
/**
* Returns the square of a normalization factor for the supplied kernel.
* The kernel can be normalized by invoking {@link #divide(KernelJAI,double)}
* with the square root of this value.
*
* @param kernel The kernel for which to compute normalization factor.
* @param type Any combinaison of {@link #HORIZONTAL} and {@link #VERTICAL}.
* @return The square of a normalization factor that could be applied on the kernel.
*/
private static double getNormalizationFactorSquared(final KernelJAI kernel, final int type) {
double sumH = 0;
double sumV = 0;
final int width = kernel.getWidth();
final int height = kernel.getHeight();
/*
* Tests the kernel with a horizontal gradient [ -1 0 1 ]
* of 1/pixel. For example, we get sumH=8 with [ -2 0 2 ]
* the horizontal Sobel kernel show on right: [ -1 0 1 ]
*/
if ((type & HORIZONTAL) != 0) {
int value = kernel.getYOrigin();
for (int y=height; --y>=0;) {
for (int x=width; --x>=0;) {
sumH += value * kernel.getElement(x,y);
}
value--;
}
}
/*
* Tests the kernel with a vertical gradient of [ -1 -2 -1 ]
* 1/pixel. For example, we get sumV=8 with the [ 0 0 0 ]
* vertical Sobel kernel show on right: [ 1 2 1 ]
*/
if ((type & VERTICAL) != 0) {
int value = kernel.getXOrigin();
for (int x=width; --x>=0;) {
for (int y=height; --y>=0;) {
sumV += value * kernel.getElement(x,y);
}
value--;
}
}
return (sumH*sumH) + (sumV*sumV);
}
/**
* Returns the normalization factor for the supplied kernel. The kernel can
* be normalized by invoking {@link #divide(KernelJAI,double)} with this factor.
*
* @param mask1 The first kernel for which to compute a normalization factor.
* @param mask2 The second kernel for which to compute a normalization factor.
* @return The normalization factor that could be applied on both kernels.
*
* @todo When the masks are symetric (e.g. Sobel, Prewitt (or Smoothed), isotropic, etc.),
* then this algorithm matches the "normalization factor" times "spatial factor"
* provided by
*
* J.J. Simpson (1990), "On the accurate detection and enhancement of oceanic features
* observed in satellite data" in Remote sensing environment, 33:17-33.
*
* However, for non-symetric masks (e.g. Kirsch), then a difference is found.
* We should provides a way to disable normalization when the user did it himself
* according some other rules than the one used here.
*/
private static double getNormalizationFactor(final KernelJAI mask1, final KernelJAI mask2) {
double factor;
factor = getNormalizationFactorSquared(mask1, HORIZONTAL|VERTICAL);
factor += getNormalizationFactorSquared(mask2, HORIZONTAL|VERTICAL);
factor = Math.sqrt(factor/2);
return factor;
}
/**
* Divides a kernel by some number.
*
* @param kernel The kernel to divide.
* @param denominator The factor to divide by.
* @return The resulting kernel.
*/
private static KernelJAI divide(KernelJAI kernel, final double denominator) {
if (denominator != 1) {
final float[] data = kernel.getKernelData();
final int length = data.length;
for (int i=0; i<length; i++) {
data[i] /= denominator;
}
kernel = new KernelJAI(kernel.getWidth(), kernel.getHeight(),
kernel.getXOrigin(), kernel.getYOrigin(), data);
}
return kernel;
}
/**
* Applies the operation on grid coverage. The default implementation looks for kernels
* specified in the {@link org.geotools.coverage.processing.OperationJAI.Parameters#parameters
* parameter block} and divide them by the distance between pixels, in the grid coverage's
* coordinate reference system.
*/
protected GridCoverage2D deriveGridCoverage(final GridCoverage2D[] sources,
final Parameters parameters)
{
if (NORMALIZE) {
final ParameterList block = parameters.parameters;
KernelJAI mask1 = (KernelJAI) block.getObjectParameter("Mask1");
KernelJAI mask2 = (KernelJAI) block.getObjectParameter("Mask2");
/*
* Normalizes the kernel in such a way that pixel values likes
* [-2 -1 0 +1 +2] will give a gradient of about 1 unit/pixel.
*/
double factor = getNormalizationFactor(mask1, mask2);
if (!(factor > 0)) {
// Do not transform if factor is 0 or NaN.
factor = 1;
}
/*
* Computes a scale factor taking in account the transformation from
* grid to coordinate system. This scale will convert gradient from
* 1 unit/pixel to 1 unit/meters or 1 unit/degrees, depending the
* coordinate reference systems axis unit.
*/
double scaleMask1 = 1;
double scaleMask2 = 1;
if (sources.length != 0) {
final MathTransform2D mtr;
mtr = ((GridGeometry2D) sources[0].getGridGeometry()).getGridToCRS2D();
if (mtr instanceof AffineTransform) {
final AffineTransform tr = (AffineTransform) mtr;
final double scaleX = XAffineTransform.getScaleX0(tr);
final double scaleY = XAffineTransform.getScaleY0(tr);
scaleMask1 = getScaleFactor(mask1, scaleX, scaleY);
scaleMask2 = getScaleFactor(mask2, scaleX, scaleY);
if (!(scaleMask1>0 && scaleMask2>0)) {
// Do not rescale if scale is 0 or NaN.
scaleMask1 = 1;
scaleMask2 = 1;
}
if (DEBUG) {
System.out.print("factor= "); System.out.println(factor);
System.out.print("scaleMask1= "); System.out.println(scaleMask1);
System.out.print("scaleMask1= "); System.out.println(scaleMask2);
}
}
}
block.setParameter("Mask1", divide(mask1, factor*scaleMask1));
block.setParameter("Mask2", divide(mask2, factor*scaleMask2));
}
return super.deriveGridCoverage(sources, parameters);
}
/**
* Derives the quantitative category for a band in the destination image.
* This implementation computes the expected gradient range from the two
* masks and the value range in the source grid coverage.
*/
protected Category deriveCategory(final Category[] categories, final Parameters parameters) {
NumberRange range = null;
Category category = categories[0];
final NumberRange samples = category.geophysics(false).getRange();
final boolean isGeophysics = (category == category.geophysics(true));
/*
* Computes a default range of output values one from the normalized kernels.
* The normalization has been done by 'deriveGridCoverage' before this method
* is invoked. The algorithm is as below:
*
* - Computes the value produced by the kernels for an artificial gradient of 1 unit/pixel.
* - Transforms into a lower gradient of 1 unit/(kernel size).
* - Transforms into a gradient of (maximal range)/(kernel size).
* - Applies an arbitrary correction factor for more convenient range in most cases.
*/
final ParameterList block = parameters.parameters;
final KernelJAI mask1 = (KernelJAI) block.getObjectParameter("Mask1");
final KernelJAI mask2 = (KernelJAI) block.getObjectParameter("Mask2");
final double size = (mask1.getWidth() + mask1.getHeight() +
mask2.getWidth() + mask2.getHeight()) / 4.0;
double factor = getNormalizationFactor(mask1, mask2) / (size-1);
if (factor>0 && !Double.isInfinite(factor)) {
range = category.geophysics(true).getRange();
final double minimum = range.getMinimum();
final double maximum = range.getMaximum();
factor *= (maximum - minimum) * DEFAULT_RANGE_SCALE;
range = NumberRange.create(0, factor);
}
if (range != null) {
category = new Category(category.getName(), DEFAULT_COLOR_PALETTE, samples, range);
return category.geophysics(isGeophysics);
}
return super.deriveCategory(categories, parameters);
}
/**
* Derives the unit of data for a band in the destination image.
* This method compute {@code sample}/{@code axis} where:
*
* <ul>
* <li>{@code sample} is the sample unit in source image.</li>
* <li>{@code axis} is the coordinate reference system axis unit.</li>
* </ul>
*/
protected Unit deriveUnit(final Unit[] units, final Parameters parameters) {
final CoordinateSystem cs = parameters.crs.getCoordinateSystem();
if (!DEBUG && units.length==1 && units[0]!=null) {
final Unit spatialUnit = cs.getAxis(0).getUnit();
for (int i=Math.min(cs.getDimension(), 2); --i>=0;) {
if (!spatialUnit.equals(cs.getAxis(i).getUnit())) {
return super.deriveUnit(units, parameters);
}
}
try {
return units[0].divide(spatialUnit);
} catch (RuntimeException exception) {
// Can't compute units... We will compute image data
// anyway, but the result will have no know unit.
// TODO: Catch a more specific exception.
}
}
return super.deriveUnit(units, parameters);
}
}