/*
 *    GeoTools - The Open Source Java GIS Toolkit
 *    http://geotools.org
 *
 *    (C) 2011, Open Source Geospatial Foundation (OSGeo)
 *    (C) 2008-2011 TOPP - www.openplans.org.
 *
 *    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.process.vector;

import javax.measure.unit.NonSI;
import javax.measure.unit.SI;
import javax.measure.unit.Unit;

import org.geotools.coverage.CoverageFactoryFinder;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridCoverageFactory;
import org.geotools.data.Query;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.factory.GeoTools;
import org.geotools.factory.Hints;
import org.geotools.filter.text.cql2.CQLException;
import org.geotools.filter.text.ecql.ECQL;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.process.ProcessException;
import org.geotools.process.factory.DescribeParameter;
import org.geotools.process.factory.DescribeProcess;
import org.geotools.process.factory.DescribeResult;
import org.geotools.referencing.CRS;
import org.opengis.coverage.grid.GridCoverage;
import org.opengis.coverage.grid.GridGeometry;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.filter.Filter;
import org.opengis.filter.expression.Expression;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.util.ProgressListener;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;

/**
 * A Process that uses a {@link HeatmapSurface} to compute a heatmap surface over a set of
 * irregular data points as a {@link GridCoverage}.
 * Heatmaps are known more formally as <i>Multivariate Kernel Density Estimation</i>.
 * <p>
 * The appearance of the heatmap is controlled by the kernel radius, 
 * which determines the "radius of influence" of input points.
 * The radius is specified by the radiusPixels parameter,
 * which is in output pixels.
 * Using pixels allows easy estimation of a value which will give a visually
 * effective result, 
 * and ensures the heatmap appearance changes to match the zoom level.
 * <p>
 * By default each input point has weight 1.
 * Optionally the weights of points may be supplied by an attribute specified by the <code>weightAttr</code> parameter.
 * <p>
 * All geometry types are allowed as input. For non-point geometries the centroid is used.
 * <p>
 * To improve performance, the surface grid can be computed at a lower resolution than the requested
 * output image using the <code>pixelsPerCell</code> parameter. 
 * The grid is upsampled to match the required image size. Upsampling uses Bilinear
 * Interpolation to maintain visual quality. This gives a large improvement in performance, with
 * minimal impact on visual quality for small cell sizes (for instance, 10 pixels or less).
 * <p>
 * To ensure that the computed surface is stable (i.e. does not display obvious edge artifacts under
 * zooming and panning), the data extent is expanded to be larger than the specified output
 * extent. The expansion distance is equal to the size of <code>radiusPixels</code> in the input
 * CRS.
 * 
 * <h3>Parameters</h3> <i>M = mandatory, O = optional</i>
 * <ul>
 * <li><b>data</b> (M) - the FeatureCollection containing the point observations
 * <li><b>radiusPixels</b> (M)- the density kernel radius, in pixels
 * <li><b>weightAttr</b> (M)- the feature type attribute containing the observed surface value
 * <li><b>pixelsPerCell</b> (O) - The pixels-per-cell value determines the resolution of the
 * computed grid. Larger values improve performance, but degrade appearance. (Default = 1)
 * <li><b>outputBBOX</b> (M) - The georeferenced bounding box of the output area
 * <li><b>outputWidth</b> (M) - The width of the output raster
 * <li><b>outputHeight</b> (M) - The height of the output raster
 * </ul>
 * The output of the process is a {@linkplain GridCoverage2D} with a single band, with cell values
 * in the range [0, 1].
 * <p>
 * Computation of the surface takes places in the CRS of the output. 
 * If the data CRS is different to the output CRS, the input points are transformed into the output CRS.
 * 
 * <h3>Using the process as a Rendering Transformation</h3>
 * 
 * This process can be used as a RenderingTransformation, since it implements the
 * <tt>invertQuery(... Query, GridGeometry)</tt> method. In this case the <code>queryBuffer</code>
 * parameter should be specified to expand the query extent appropriately. The output raster
 * parameters may be provided from the request extents, using the following SLD environment
 * variables:
 * <ul>
 * <li><b>outputBBOX</b> - env var = <tt>wms_bbox</tt>
 * <li><b>outputWidth</b> - env var = <tt>wms_width</tt>
 * <li><b>outputHeight</b> - env var = <tt>wms_height</tt>
 * </ul>
 * When used as an Rendering Transformation the data query is rewritten to expand the query BBOX, to
 * ensure that enough data points are queried to make the computed surface stable under panning and
 * zooming.
 * 
 * <p>
 * 
 * @author Martin Davis - OpenGeo
 * 
 */
@DescribeProcess(title = "Heatmap", description = "Computes a heatmap surface over a set of data points and outputs as a single-band raster.")
public class HeatmapProcess implements VectorProcess {

    @DescribeResult(name = "result", description = "Output raster")
    public GridCoverage2D execute(

            // process data
            @DescribeParameter(name = "data", description = "Input features") SimpleFeatureCollection obsFeatures,

            // process parameters
            @DescribeParameter(name = "radiusPixels", description = "Radius of the density kernel in pixels") Integer argRadiusPixels,
            //CR#temp change the type of valueAttr to Expression
            //@DescribeParameter(name = "weightAttr", description = "Name of the attribute to use for data point weight", min = 0, max = 1) String valueAttr,
            @DescribeParameter(name = "weightAttr", description = "Name of the attribute to use for data point weight", min = 0, max = 1) Expression valueAttr,
            @DescribeParameter(name = "pixelsPerCell", description = "Resolution at which to compute the heatmap (in pixels). Default = 1", defaultValue="1", min = 0, max = 1) Integer argPixelsPerCell,

            // output image parameters
            @DescribeParameter(name = "outputBBOX", description = "Bounding box of the output") ReferencedEnvelope argOutputEnv,
            @DescribeParameter(name = "outputWidth", description = "Width of output raster in pixels") Integer argOutputWidth,
            @DescribeParameter(name = "outputHeight", description = "Height of output raster in pixels") Integer argOutputHeight,

            ProgressListener monitor) throws ProcessException {

        /**
         * -------- Extract required information from process arguments -------------
         */
        int pixelsPerCell = 1;
        if (argPixelsPerCell != null && argPixelsPerCell > 1) {
            pixelsPerCell = argPixelsPerCell;
        }
        int outputWidth = argOutputWidth;
        int outputHeight = argOutputHeight;
        int gridWidth = outputWidth;
        int gridHeight = outputHeight;
        if (pixelsPerCell > 1) {
            gridWidth = outputWidth / pixelsPerCell;
            gridHeight = outputHeight / pixelsPerCell;
        }

        /**
         * Compute transform to convert input coords into output CRS
         */
        CoordinateReferenceSystem srcCRS = obsFeatures.getSchema().getCoordinateReferenceSystem();
        CoordinateReferenceSystem dstCRS = argOutputEnv.getCoordinateReferenceSystem();
        MathTransform trans = null;
        try {
            trans = CRS.findMathTransform(srcCRS, dstCRS);
        } catch (FactoryException e) {
            throw new ProcessException(e);
        }

        //------------ Kernel Radius
        /*
         * // not used for now - only pixel radius values are supported double
         * distanceConversionFactor = distanceConversionFactor(srcCRS, dstCRS); double dstRadius =
         * argRadius * distanceConversionFactor;
         */
        int radiusCells = 100;
        if (argRadiusPixels != null)
            radiusCells = argRadiusPixels;
        if (pixelsPerCell > 1) {
            radiusCells /= pixelsPerCell;
        }


        /**
         * -------------- Extract the input observation points -----------
         */
        HeatmapSurface heatMap = new HeatmapSurface(radiusCells, argOutputEnv, gridWidth,
                gridHeight);
        try {
            extractPoints(obsFeatures, valueAttr, trans, heatMap);
        } catch (CQLException e) {
            throw new ProcessException(e);
        }

        /**
         * --------------- Do the processing ------------------------------
         */
        // Stopwatch sw = new Stopwatch();
        // compute the heatmap at the specified resolution
        float[][] heatMapGrid = heatMap.computeSurface();

        // flip now, since grid size may be smaller
        heatMapGrid = flipXY(heatMapGrid);

        // upsample to output resolution if necessary
        float[][] outGrid = heatMapGrid;
        if (pixelsPerCell > 1)
            outGrid = upsample(heatMapGrid, -999, outputWidth, outputHeight);

        // convert to the GridCoverage2D required for output
        GridCoverageFactory gcf = CoverageFactoryFinder.getGridCoverageFactory(GeoTools.getDefaultHints());
        GridCoverage2D gridCov = gcf.create("Process Results", outGrid, argOutputEnv);

        // System.out.println("**************  Heatmap computed in " + sw.getTimeString());

        return gridCov;
    }

    /*
     * An approximate value for the length of a degree at the equator in meters. This doesn't have
     * to be precise, since it is only used to convert values which are themselves rough
     * approximations.
     */
    private static final double METRES_PER_DEGREE = 111320;

    private static double distanceConversionFactor(CoordinateReferenceSystem srcCRS,
            CoordinateReferenceSystem dstCRS) {
        Unit<?> srcUnit = srcCRS.getCoordinateSystem().getAxis(0).getUnit();
        Unit<?> dstUnit = dstCRS.getCoordinateSystem().getAxis(0).getUnit();
        if (srcUnit == dstUnit) {
            return 1;
        } else if (srcUnit == NonSI.DEGREE_ANGLE && dstUnit == SI.METER) {
            return METRES_PER_DEGREE;
        } else if (srcUnit == SI.METER && dstUnit == NonSI.DEGREE_ANGLE) {
            return 1.0 / METRES_PER_DEGREE;
        }
        throw new IllegalStateException("Unable to convert distances from " + srcUnit + " to "
                + dstUnit);
    }

    /**
     * Flips an XY matrix along the X=Y axis, and inverts the Y axis. Used to convert from
     * "map orientation" into the "image orientation" used by GridCoverageFactory. The surface
     * interpolation is done on an XY grid, with Y=0 being the bottom of the space. GridCoverages
     * are stored in an image format, in a YX grid with Y=0 being the top.
     * 
     * @param grid the grid to flip
     * @return the flipped grid
     */
    private float[][] flipXY(float[][] grid) {
        int xsize = grid.length;
        int ysize = grid[0].length;

        float[][] grid2 = new float[ysize][xsize];
        for (int ix = 0; ix < xsize; ix++) {
            for (int iy = 0; iy < ysize; iy++) {
                int iy2 = ysize - iy - 1;
                grid2[iy2][ix] = grid[ix][iy];
            }
        }
        return grid2;
    }

    private float[][] upsample(float[][] grid, float noDataValue, int width, int height) {
        BilinearInterpolator bi = new BilinearInterpolator(grid, noDataValue);
        float[][] outGrid = bi.interpolate(width, height, false);
        return outGrid;
    }

    /**
     * Given a target query and a target grid geometry returns the query to be used to read the
     * input data of the process involved in rendering. In this process this method is used to:
     * <ul>
     * <li>determine the extent & CRS of the output grid
     * <li>expand the query envelope to ensure stable surface generation
     * <li>modify the query hints to ensure point features are returned
     * </ul>
     * Note that in order to pass validation, all parameters named here must also appear in the
     * parameter list of the <tt>execute</tt> method, even if they are not used there.
     * 
     * @param argRadiusPixels the feature type attribute that contains the observed surface value
     * @param targetQuery the query used against the data source
     * @param targetGridGeometry the grid geometry of the destination image
     * @return The transformed query
     */
    public Query invertQuery(
            @DescribeParameter(name = "radiusPixels", description = "Radius to use for the kernel", min = 0, max = 1) Integer argRadiusPixels,
            // output image parameters
            @DescribeParameter(name = "outputBBOX", description = "Georeferenced bounding box of the output") ReferencedEnvelope argOutputEnv,
            @DescribeParameter(name = "outputWidth", description = "Width of the output raster") Integer argOutputWidth,
            @DescribeParameter(name = "outputHeight", description = "Height of the output raster") Integer argOutputHeight,

            Query targetQuery, GridGeometry targetGridGeometry) throws ProcessException {

        // TODO: handle different CRSes in input and output
        
        int radiusPixels = argRadiusPixels > 0 ? argRadiusPixels : 0;
        // input parameters are required, so should be non-null
        double queryBuffer = radiusPixels / pixelSize(argOutputEnv, argOutputWidth, argOutputHeight);
        /*
         * if (argQueryBuffer != null) { queryBuffer = argQueryBuffer; }
         */
        targetQuery.setFilter(expandBBox(targetQuery.getFilter(), queryBuffer));

        // clear properties to force all attributes to be read
        // (required because the SLD processor cannot see the value attribute specified in the
        // transformation)
        // TODO: set the properties to read only the specified value attribute
        targetQuery.setProperties(null);

        // set the decimation hint to ensure points are read
        Hints hints = targetQuery.getHints();
        hints.put(Hints.GEOMETRY_DISTANCE, 0.0);

        return targetQuery;
    }

    private double pixelSize(ReferencedEnvelope outputEnv, int outputWidth, int outputHeight)
    {
        // error-proofing
        if (outputEnv.getWidth() <= 0) return 0;
        // assume view is isotropic
        return outputWidth / outputEnv.getWidth();
    }
    
    protected Filter expandBBox(Filter filter, double distance) {
        return (Filter) filter.accept(new BBOXExpandingFilterVisitor(distance, distance, distance,
                distance), null);
    }

    /**
     * Extract points from a feature collection, and stores them in the heatmap
     *
     * @param obsPoints features to extract
     * @param attrName expression or property name used to evaluate the geometry from a feature
     * @param trans transform for extracted points
     * @param heatMap heatmap to add points to
     * @throws CQLException if attrName can't be parsed
     */
    protected void extractPoints(SimpleFeatureCollection obsPoints, Expression attrName,         //CR#temp change the type of AttrName to Expression
            MathTransform trans, HeatmapSurface heatMap) throws CQLException {
        Expression attrExpr = attrName;  //CR#temp
        /* //CR#temp
        Expression attrExpr = null;
        if (attrName != null) {
            attrExpr = ECQL.toExpression(attrName);
        }
        */
        SimpleFeatureIterator obsIt = obsPoints.features();

        double[] srcPt = new double[2];
        double[] dstPt = new double[2];

        int i = 0;
        try {
            while (obsIt.hasNext()) {
                SimpleFeature feature = obsIt.next();

                try {
                    // get the weight value, if any
                    double val = 1;
                    if (attrExpr != null) {
                        val = getPointValue(feature, attrExpr);
                    }

                    // get the point location from the geometry
                    Geometry geom = (Geometry) feature.getDefaultGeometry();
                    Coordinate p = getPoint(geom);
                    srcPt[0] = p.x;
                    srcPt[1] = p.y;
                    trans.transform(srcPt, 0, dstPt, 0, 1);
                    Coordinate pobs = new Coordinate(dstPt[0], dstPt[1], val);

                    heatMap.addPoint(pobs.x, pobs.y, val);
                } catch (Exception e) {
                    // just carry on for now (debugging)
                    // throw new ProcessException("Expression " + attrExpr +
                    // " failed to evaluate to a numeric value", e);
                }
            }
        } finally {
            obsIt.close();
        }
    }

    /**
     * Gets a point to represent the Geometry. If the Geometry is a point, this is returned.
     * Otherwise, the centroid is used.
     * 
     * @param g the geometry to find a point for
     * @return a point representing the Geometry
     */
    private static Coordinate getPoint(Geometry g) {
        if (g.getNumPoints() == 1)
            return g.getCoordinate();
        return g.getCentroid().getCoordinate();
    }

    /**
     * Gets the value for a point from the supplied attribute.
     * The value is checked for validity,
     * and a default of 1 is used if necessary.
     * 
     * @param feature the feature to extract the value from
     * @param attrExpr the expression specifying the attribute to read
     * @return the value for the point
     */
    private static double getPointValue(SimpleFeature feature, Expression attrExpr) {
        Double valObj = attrExpr.evaluate(feature, Double.class);
        if (valObj != null) {
            return valObj;
        }
        return 1;
    }
}