    001    //$HeadURL: https://svn.wald.intevation.org/svn/deegree/base/branches/2.4_testing/src/org/deegree/processing/raster/filter/Convolve.java $
    002    /*----------------------------------------------------------------------------
    003     This file is part of deegree, http://deegree.org/
    004     Copyright (C) 2001-2009 by:
    005       Department of Geography, University of Bonn
    006     and
    007       lat/lon GmbH
    008
    009     This library is free software; you can redistribute it and/or modify it under
    010     the terms of the GNU Lesser General Public License as published by the Free
    011     Software Foundation; either version 2.1 of the License, or (at your option)
    012     any later version.
    013     This library is distributed in the hope that it will be useful, but WITHOUT
    014     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
    015     FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
    016     details.
    017     You should have received a copy of the GNU Lesser General Public License
    018     along with this library; if not, write to the Free Software Foundation, Inc.,
    019     59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
    020
    021     Contact information:
    022
    023     lat/lon GmbH
    024     Aennchenstr. 19, 53177 Bonn
    025     Germany
    026     http://lat-lon.de/
    027
    028     Department of Geography, University of Bonn
    029     Prof. Dr. Klaus Greve
    030     Postfach 1147, 53001 Bonn
    031     Germany
    032     http://www.geographie.uni-bonn.de/deegree/
    033
    034     e-mail: info@deegree.org
    035    ----------------------------------------------------------------------------*/
    036    package org.deegree.processing.raster.filter;
    037
    038    /**
    039     * This class offeres some convenience methods for convolving a data matrix. <br>
    040     * Convolution filtering is often used to reduce the effects of noise in images or to sharpen the
    041     * detail in rasters. Convolution filtering is a form of spatial filtering that computes each output
    042     * sample by multiplying elements of a kernel with the samples surrounding a particular source
    043     * sample. <br>
    044     * Convolution filtering operates on a group of input values surrounding a center pixel. The
    045     * adjoining values provide important information about values trends in the area of the values
    046     * being processed. <br>
    047     * Convolution filtering moves across the source raster, cell by cell, placing resulting values into
    048     * the destination raster. The resulting value of each source cell depends on the group of values
    049     * surrounding the source cell. Using the cell values of the source cell's neighbors, the
    050     * convolution process calculates the spatial frequency activity in the area, making it possible to
    051     * filter the values based on the spatial frequency of the area. <br>
    052     * Convolution filtering uses a convolve kernel, containing an array of convolution coefficient
    053     * values, called key elements. The array is not restricted to any particular size, and does not
    054     * even have to be square. The kernel can be 1 x 1, 3 x 3, 5 x 5, M x N, and so on. A larger kernel
    055     * size affords a more precise filtering operation by increasing the number of neighboring cells
    056     * used in the calculation. However, the kernel cannot be bigger in any dimension than the source
    057     * data. Also, the larger the kernel, the more computations that are required to be performed. For
    058     * example, given a 640 x 480 raster and a 3 x 3 kernel, the convolve operation requires over five
    059     * million total multiplications and additions. <br>
    060     * The convolution filtering operation computes each output sample by multiplying the key elements
    061     * of the kernel with the samples surrounding a particular source cell. For each destination cell,
    062     * the kernel is rotated 180 degrees and its key element is placed over the source pixel
    063     * corresponding with the destination pixel. The key elements are multiplied with the source cell
    064     * value under them, and the resulting products are summed together to produce the destination
    065     * sample value. <br>
    066     * The selection of the weights for the key elements determines the nature of the filtering action,
    067     * such as high-pass or low-pass. If the values of the key elements are the reciprocal of the number
    068     * of key elements in the kernel (for example, 1/9 for a 3 x 3 kernel), the result is a conventional
    069     * low-pass averaging process. If the weights are altered, certain cells in the kernel will have an
    070     * increased or decreased influence in the average. <br>
    071     * (modified from JAI documentation)
    072     *
    073     * @version $Revision: 18195 $
    074     * @author <a href="mailto:poth@lat-lon.de">Andreas Poth</a>
    075     * @author last edited by: $Author: mschneider $
    076     *
    077     * @version 1.0. $Revision: 18195 $, $Date: 2009-06-18 17:55:39 +0200 (Do, 18. Jun 2009) $
    078     *
    079     * @since 2.0
    080     */
    081    public class Convolve {
    082
    083        /**
    084         * convolves the passed float matrix by applying the passed kernel. A kernel for lowpass
    085         * filtering may looks like:
    086         *
    087         * <pre>
    088         *  float[][] kernel = new float[5][5];
    089         *  for ( int i = 0; i &lt; kernel.length; i++ ) {
    090         *       for ( int j = 0; j &lt; kernel[i].length; j++ ) {
    091         *           kernel[i][j] = 1;
    092         *       }
    093         *   }
    094         * &lt;pre&gt;
    095         * which results in a strong smoothing of the raster.
    096         * <BR>
    097         * A kernel for highpass filtering may looks like:
    098         * &lt;pre&gt;
    099         *  float[][] kernel = new float[3][3];
    100         *  for ( int i = 0; i &lt; kernel.length; i++ ) {
    101         *       for ( int j = 0; j &lt; kernel[i].length; j++ ) {
    102         *           kernel[i][j] = -1;
    103         *       }
    104         *   }
    105         *   kernel[kernel.length/2][kernel[0].length/2] = 2;
    106         * &lt;/pre&gt;
    107         * Notice that the sum of all kernel values does not must be
    108         * == 1! Correct weighting will be done by this method.
    109         *
    110         * @param source to apply the kernel to.
    111         * @param kernel to apply
    112         * @return the raster after applying the given kernel.
    113         * @throws RasterFilterException if the kernel does not have an odd number of cells.
    114         *
    115         */
    116        public static float[][] perform( float[][] source, float[][] kernel )
    117                                throws RasterFilterException {
    118
    119            if ( kernel.length % 2 == 0 || kernel[0].length % 2 == 0 ) {
    120                throw new RasterFilterException( "A Kernel must have an odd number of cells in x- and y-direction" );
    121            }
    122
    123            kernel = rotateKernel( kernel );
    124
    125            float[][] dest = new float[source.length][source[0].length];
    126
    127            int ww = kernel.length;
    128            int hh = kernel[0].length;
    129            int xOrigin = ww / 2;
    130            int yOrigin = hh / 2;
    131            for ( int y = 0; y < source.length; y++ ) {
    132                for ( int x = 0; x < source[y].length; x++ ) {
    133                    float g = 0;
    134                    float v = 0;
    135                    for ( int i = -xOrigin; i < -xOrigin + ww; i++ ) {
    136                        for ( int j = -yOrigin; j < -yOrigin + ww; j++ ) {
    137                            if ( y + i >= 0 && y + i < source.length && x + j >= 0 && x + j < source[y].length ) {
    138                                v += source[y + i][x + j] * kernel[xOrigin + i][yOrigin + j];
    139                                g += kernel[xOrigin + i][yOrigin + j];
    140                            }
    141                        }
    142                    }
    143                    dest[y][x] = v / g;
    144                }
    145            }
    146
    147            return dest;
    148        }
    149
    150        /**
    151         * rotates a convolution kernel by 180°
    152         *
    153         * @param kernel
    154         * @return the rotated kernel.
    155         */
    156        private static float[][] rotateKernel( float[][] kernel ) {
    157            for ( int i = 0; i < kernel.length / 2; i++ ) {
    158                for ( int j = 0; j < kernel[i].length / 2; j++ ) {
    159                    float v = kernel[i][j];
    160                    kernel[i][j] = kernel[kernel.length - i - 1][kernel[i].length - j - 1];
    161                    kernel[kernel.length - i - 1][kernel[i].length - j - 1] = v;
    162                }
    163            }
    164            return kernel;
    165        }
    166
    167    }