/*
 * GridGaussLowPassOperator.java
 *
 */
package edu.oregonstate.cartography.grid.operators;

import edu.oregonstate.cartography.grid.Grid;
import java.text.DecimalFormat;

/**
 * Gaussian blur or low pass filter. Uses the fact that a 2D Gaussian
 * convolution can be replaced by a horizontal and a vertical 1D convolution. A
 * horizontal 1D convolution is applied, then the grid is transposed, the
 * horizontal 1D convolution applied again, and the grid transposed again. This
 * is much faster than a combination of a horizontal and a vertical convolution,
 * due to the fact that horizontal rows are stored in contiguous locations. The
 * horizontal 1D convolution and the transposition operators are merged into one
 * multi-threaded operator. See http://en.wikipedia.org/wiki/Gaussian_blur
 * August 26, 2010, and April 14, 2011.
 *
 * @author Bernhard Jenny, Institute of Cartography, ETH Zurich.
 */
public class GridGaussLowPassOperator implements GridOperator {

    /**
     * Standard deviation of the Gaussian distribution. Higher values produce
     * stronger smoothing.
     */
    private double std = 0.8;

    /**
     * The size of the kernel relative to the standard deviation. The kernel's
     * dimension in pixels in one direction is: relativeFilterSize * std A value
     * of 6 is sufficient for most applications. If gradients are computed after
     * filtering with a kernel size of 6, however, artifacts (vertical and
     * horizontal stripes) are likely to become visible.
     */
    private int relativeFilterSize = 8;

    /**
     * Applies horizontal Gaussian convolution and stores results in a
     * transposed grid.
     */
    private class HorizontalTransposedConvolution extends ThreadedGridOperator {

        /**
         * Create a transposed grid
         *
         * @param srcGrid
         * @return
         */
        @Override
        protected Grid initDestinationGrid(Grid srcGrid) {
            final int nrows = srcGrid.getRows();
            final int ncols = srcGrid.getCols();
            Grid tempTransposedGrid = new Grid(nrows, ncols, srcGrid.getCellSize());
            tempTransposedGrid.setWest(srcGrid.getWest());
            tempTransposedGrid.setSouth(srcGrid.getSouth());
            return tempTransposedGrid;
        }

        @Override
        public Grid operate(Grid src, Grid dst) {
            if (src.getCols() != dst.getRows() || src.getRows() != dst.getCols()) {
                throw new IllegalStateException("destination grid has wrong size");
            }
            return super.operate(src, dst);
        }

        private float convolveWithNaN(float[] srcRow, int col, float[] kernel) {
            final int nCols = srcRow.length;
            final int halfFilterSize = kernelSize() / 2;
            float sum = 0;
            float coefSum = 0;
            for (int f = -halfFilterSize; f <= halfFilterSize; f++) {
                int c = f + col;
                if (c < 0 || c >= nCols) {
                    continue;
                }
                float v = srcRow[c];
                if (Float.isNaN(v)) {
                    continue;
                }
                float w = kernel[f + halfFilterSize];
                sum += v * w;
                coefSum += w;
            }
            return sum / coefSum;
        }

        @Override
        public void operate(Grid src, Grid dst, int startRow, int endRow) {

            final int ncols = src.getCols();
            final int halfFilterSize = kernelSize() / 2;
            final float[] kernel = kernel();

            for (int row = startRow; row < endRow; row++) {
                final float[][] dstGrid = dst.getGrid();
                final float[] srcRow = src.getGrid()[row];

                // convolve left border area
                final int maxCol = Math.min(halfFilterSize, ncols);
                for (int col = 0; col < maxCol; col++) {
                    float sum = 0;
                    float coefSum = 0;
                    for (int f = -col; f <= halfFilterSize; f++) {
                        if (col + f < ncols) {
                            final float s = kernel[f + halfFilterSize];
                            sum += srcRow[col + f] * s;
                            coefSum += s;
                        }
                    }
                    float v = sum / coefSum;
                    if (Float.isNaN(v)) {
                        v = convolveWithNaN(srcRow, col, kernel);
                    }
                    dstGrid[col][row] = v; // transposed destination
                }

                // convolve center area
                for (int col = halfFilterSize; col < ncols - halfFilterSize; col++) {
                    float sum = 0;
                    for (int c = col - halfFilterSize, f = 0; c <= col + halfFilterSize; c++, f++) {
                        sum += srcRow[c] * kernel[f];
                    }
                    if (Float.isNaN(sum)) {
                        sum = convolveWithNaN(srcRow, col, kernel);
                    }
                    dstGrid[col][row] = sum; // transposed destination
                }

                // convolve right border area
                final int minCol = Math.max(0, ncols - halfFilterSize);
                for (int col = minCol; col < ncols; col++) {
                    float sum = 0;
                    float coefSum = 0;
                    for (int f = -halfFilterSize; f < ncols - col; f++) {
                        if (col + f >= 0) {
                            final float s = kernel[f + halfFilterSize];
                            sum += srcRow[col + f] * s;
                            coefSum += s;
                        }
                    }
                    float v = sum / coefSum;
                    if (Float.isNaN(v)) {
                        v = convolveWithNaN(srcRow, col, kernel);
                    }
                    dstGrid[col][row] = v; // transposed destination
                }
            }
        }

        @Override
        public String getName() {
            return "Horizontal Transposed 1D Convolution";
        }
    }

    /**
     * Creates a new instance of GridGaussLowPassOperator
     */
    public GridGaussLowPassOperator() {
    }

    /**
     * Creates a new instance of GridGaussLowPassOperator
     *
     * @param std
     */
    public GridGaussLowPassOperator(double std) {
        setStandardDeviation(std);
    }

    @Override
    public String getName() {
        return "Gauss Low Pass";
    }

    public String kernelToString() {
        float[] kernel = kernel();
        StringBuilder sb = new StringBuilder();
        DecimalFormat df = new DecimalFormat("0.####");
        sb.append("Kernel ");
        sb.append(kernel.length);
        sb.append("x");
        sb.append(kernel.length);
        sb.append(System.getProperty("line.separator"));
        for (int i = 0; i < kernel.length; i++) {
            sb.append(df.format(kernel[i])).append(" ");
        }
        return sb.toString();
    }

    /**
     * Evaluates the Gaussian function.
     *
     * @param x Evaluate at distance x from 0.
     * @return The vertical distance at position x.
     */
    private double gaussian(double x) {
        final double stdSqr2 = 2 * std * std;
        return Math.exp(-x * x / stdSqr2) / Math.sqrt(stdSqr2 * Math.PI);
    }

    /**
     * Returns the size of the kernel.
     *
     * @return
     */
    private int kernelSize() {
        // compute odd filter size and even half size
        // the size is relativeFilterSize * the standard deviation on both 
        // sides of the origin. Kernel values outside this range are neglected.
        int filterSize = (int) Math.ceil(relativeFilterSize * std);
        if (filterSize % 2 == 0) {
            filterSize += 1;
        }
        return filterSize;
    }

    /**
     * Computes the coefficients for the Gaussian kernel
     *
     * @param filterSize
     * @return
     */
    private float[] kernel() {
        int size = kernelSize();
        float coef[] = new float[size];
        for (int i = 0; i <= size / 2; i++) {
            coef[i] = (float) gaussian(size / 2 - i);
            coef[size - i - 1] = coef[i];
        }

        // normalize sum of coefficients
        float coefSum = 0;
        for (int i = 0; i < size; i++) {
            coefSum += coef[i];
        }
        for (int i = 0; i < size; i++) {
            coef[i] /= coefSum;
        }
        return coef;
    }

    @Override
    public Grid operate(Grid grid) {
        Grid dst = new Grid(grid.getCols(), grid.getRows(), grid.getCellSize());
        dst.setWest(grid.getWest());
        dst.setSouth(grid.getSouth());
        return operate(grid, dst);
    }

    public Grid operate(Grid src, Grid dst, double std) {
        setStandardDeviation(std);
        return operate(src, dst);
    }

    public Grid operate(Grid src, Grid dst, double std, int relativeFilterSize) {
        setStandardDeviation(std);
        setRelativeFilterSize(relativeFilterSize);
        return operate(src, dst);
    }

    public Grid operate(Grid src, Grid dst) {
        if (std == 0) {
            if (src == dst) {
                return dst;
            } else {
                return new GridCopyOperator().operate(src, dst);
            }
        }

        HorizontalTransposedConvolution hop = new HorizontalTransposedConvolution();
        Grid transposedGrid = hop.operate(src);
        return hop.operate(transposedGrid, dst);
    }

    /**
     * Get the standard deviation of the Gaussian distribution.
     *
     * @return the standard deviation
     */
    public double getStandardDeviation() {
        return std;
    }

    /**
     * Set the standard deviation of the Gaussian distribution. If std is 0, no
     * filter is applied.
     *
     * @param std the standard deviation to set
     */
    public final void setStandardDeviation(double std) {
        if (std < 0) {
            throw new IllegalArgumentException("negative standard deviation");
        }
        this.std = std;
    }

    /**
     * Get the size of the kernel, relative to the standard deviation. The
     * kernel size in pixels in one dimension is relativeFilterSize * std
     *
     * @return the relativeFilterSize
     */
    public int getRelativeFilterSize() {
        return relativeFilterSize;
    }

    /**
     * Set the size of the kernel, relative to the standard deviation. The
     * kernel size in pixels in one dimension is relativeFilterSize * std
     *
     * @param relativeFilterSize the relativeFilterSize to set
     */
    public void setRelativeFilterSize(int relativeFilterSize) {
        this.relativeFilterSize = relativeFilterSize;
    }
}
