/**
 * Copyright (c) 2011, The University of Southampton and the individual contributors.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 *
 *   *  Redistributions of source code must retain the above copyright notice,
 *      this list of conditions and the following disclaimer.
 *
 *   *  Redistributions in binary form must reproduce the above copyright notice,
 *      this list of conditions and the following disclaimer in the documentation
 *      and/or other materials provided with the distribution.
 *
 *   *  Neither the name of the University of Southampton nor the names of its
 *      contributors may be used to endorse or promote products derived from this
 *      software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
package org.openimaj.image.feature.dense.gradient.dsift;

import org.openimaj.feature.local.list.LocalFeatureList;
import org.openimaj.feature.local.list.MemoryLocalFeatureList;
import org.openimaj.image.FImage;
import org.openimaj.image.processing.convolution.FImageConvolveSeparable;
import org.openimaj.image.processing.convolution.FImageGradients;
import org.openimaj.math.geometry.shape.Rectangle;
import org.openimaj.util.array.ArrayUtils;

/**
 * Implementation of a dense SIFT feature extractor for {@link FImage}s.
 * Extracts upright SIFT features at a single scale on a grid.
 * <p>
 * Implementation inspired by the
 * <a href="http://www.vlfeat.org/api/dsift.html#dsift-usage">VLFeat
 * extractor</a>.
 * <p>
 * <b>Implementation Notes</b>. The analyser is not thread-safe, however, it is
 * safe to reuse the analyser. In multi-threaded environments, a separate
 * instance must be made for each thread. Internally, this implementation
 * allocates memory for the gradient images, and if possible re-uses these
 * between calls. Re-use requires that the input image is the same size between
 * calls to the analyser.
 *
 * @see "http://www.vlfeat.org/api/dsift.html#dsift-usage"
 *
 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
 *
 */
public class DenseSIFT extends AbstractDenseSIFT<FImage> {
        static class WorkingData {
                /**
                 * minimum X bound for sampling the descriptors (inclusive)
                 */
                protected int boundMinX;

                /**
                 * maximum X bound for sampling the descriptors (inclusive)
                 */
                protected int boundMaxX;

                /**
                 * minimum Y bound for sampling the descriptors (inclusive)
                 */
                protected int boundMinY;

                /**
                 * maximum X bound for sampling the descriptors (inclusive)
                 */
                protected int boundMaxY;

                /**
                 * Quantised orientation/gradient magnitude data. Each element
                 * corresponds to an angular bin of the orientations, and the individual
                 * images correspond to the gradient magnitudes for the respective
                 * orientations.
                 */
                protected FImage[] gradientMagnitudes;

                /**
                 * Setup the required space for holding gradient maps and descriptors
                 *
                 * @param image
                 *            the image being analysed
                 */
                protected void setupWorkingSpace(FImage image, DenseSIFT dsift) {
                        if (gradientMagnitudes == null) {
                                gradientMagnitudes = new FImage[dsift.numOriBins];
                        }

                        if (gradientMagnitudes[0] == null || gradientMagnitudes[0].width != image.width
                                        || gradientMagnitudes[0].height != image.height)
                        {
                                for (int i = 0; i < dsift.numOriBins; i++)
                                        gradientMagnitudes[i] = new FImage(image.width, image.height);
                        }

                        final int rangeX = boundMaxX - boundMinX - (dsift.numBinsX - 1) * dsift.binWidth;
                        final int rangeY = boundMaxY - boundMinY - (dsift.numBinsY - 1) * dsift.binHeight;

                        final int numWindowsX = (rangeX >= 0) ? rangeX / dsift.stepX + 1 : 0;
                        final int numWindowsY = (rangeY >= 0) ? rangeY / dsift.stepY + 1 : 0;

                        final int numFeatures = numWindowsX * numWindowsY;

                        dsift.descriptors = new float[numFeatures][dsift.numOriBins * dsift.numBinsX * dsift.numBinsY];
                        dsift.energies = new float[numFeatures];
                }
        }

        /**
         * Step size of sampling window in x-direction (in pixels)
         */
        protected int stepX = 5;

        /**
         * Step size of sampling window in y-direction (in pixels)
         */
        protected int stepY = 5;

        /**
         * Width of a single bin of the sampling window (in pixels). Sampling window
         * width is this multiplied by #numBinX.
         */
        protected int binWidth = 5;

        /**
         * Height of a single bin of the sampling window (in pixels). Sampling
         * window height is this multiplied by #numBinY.
         */
        protected int binHeight = 5;

        /**
         * Number of spatial bins in the X direction
         */
        protected int numBinsX = 4;

        /**
         * Number of spatial bins in the Y direction
         */
        protected int numBinsY = 4;

        /** The number of orientation bins */
        protected int numOriBins = 8;

        /**
         * Size of the Gaussian window (in relative to of the size of a bin)
         */
        protected float gaussianWindowSize = 2f;

        /**
         * Threshold for clipping the SIFT features
         */
        protected float valueThreshold = 0.2f;

        protected volatile WorkingData data = new WorkingData();

        /**
         * Extracted descriptors
         */
        protected volatile float[][] descriptors;

        /**
         * Descriptor energies
         */
        protected volatile float[] energies;

        /**
         * Construct with the default configuration: standard SIFT geometry (4x4x8),
         * 5px x 5px spatial bins, 5px step size, gaussian window size of 2 and
         * value threshold of 0.2.
         */
        public DenseSIFT() {
        }

        /**
         * Construct with the given step size (for both x and y) and binSize. All
         * other values are the defaults.
         *
         * @param step
         *            the step size
         * @param binSize
         *            the spatial bin size
         */
        public DenseSIFT(int step, int binSize) {
                this.binWidth = binSize;
                this.binHeight = binSize;
                this.stepX = step;
                this.stepY = step;
        }

        /**
         * Construct with the given configuration. The gaussian window size is set
         * to 2, and value threshold to 0.2.
         *
         * @param stepX
         *            step size in x direction
         * @param stepY
         *            step size in y direction
         * @param binWidth
         *            width of spatial bins
         * @param binHeight
         *            height of spatial bins
         * @param numBinsX
         *            number of bins in x direction for each descriptor
         * @param numBinsY
         *            number of bins in y direction for each descriptor
         * @param numOriBins
         *            number of orientation bins for each descriptor
         */
        public DenseSIFT(int stepX, int stepY, int binWidth, int binHeight, int numBinsX, int numBinsY, int numOriBins) {
                this(stepX, stepY, binWidth, binHeight, numBinsX, numBinsY, numOriBins, 2f);
        }

        /**
         * Construct with the given configuration. The value threshold is set to
         * 0.2.
         *
         * @param stepX
         *            step size in x direction
         * @param stepY
         *            step size in y direction
         * @param binWidth
         *            width of spatial bins
         * @param binHeight
         *            height of spatial bins
         * @param numBinsX
         *            number of bins in x direction for each descriptor
         * @param numBinsY
         *            number of bins in y direction for each descriptor
         * @param numOriBins
         *            number of orientation bins for each descriptor
         * @param gaussianWindowSize
         *            the size of the gaussian weighting window
         */
        public DenseSIFT(int stepX, int stepY, int binWidth, int binHeight, int numBinsX, int numBinsY, int numOriBins,
                        float gaussianWindowSize)
        {
                this(stepX, stepY, binWidth, binHeight, numBinsX, numBinsY, numOriBins, gaussianWindowSize, 0.2f);
        }

        /**
         * Construct with the given configuration.
         *
         * @param stepX
         *            step size in x direction
         * @param stepY
         *            step size in y direction
         * @param binWidth
         *            width of spatial bins
         * @param binHeight
         *            height of spatial bins
         * @param numBinsX
         *            number of bins in x direction for each descriptor
         * @param numBinsY
         *            number of bins in y direction for each descriptor
         * @param numOriBins
         *            number of orientation bins for each descriptor
         * @param gaussianWindowSize
         *            the size of the gaussian weighting window
         * @param valueThreshold
         *            the threshold for clipping features
         */
        public DenseSIFT(int stepX, int stepY, int binWidth, int binHeight, int numBinsX, int numBinsY, int numOriBins,
                        float gaussianWindowSize, float valueThreshold)
        {
                this.binWidth = binWidth;
                this.binHeight = binHeight;
                this.stepX = stepX;
                this.stepY = stepY;
                this.numBinsX = numBinsX;
                this.numBinsY = numBinsY;
                this.numOriBins = numOriBins;
                this.gaussianWindowSize = gaussianWindowSize;
                this.valueThreshold = valueThreshold;
        }

        private float[] buildKernel(int binSize, int numBins, int binIndex, float windowSize) {
                final int kernelSize = 2 * binSize - 1;
                final float[] kernel = new float[kernelSize];
                final float delta = binSize * (binIndex - 0.5F * (numBins - 1));

                final float sigma = binSize * windowSize;

                for (int x = -binSize + 1, i = 0; x <= +binSize - 1; x++, i++) {
                        final float z = (x - delta) / sigma;

                        kernel[i] = (1.0F - Math.abs(x) / binSize) * ((binIndex >= 0) ? (float) Math.exp(-0.5F * z * z) : 1.0F);
                }

                return kernel;
        }

        /**
         * Extract the DSIFT features
         */
        protected void extractFeatures() {
                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;

                for (int biny = 0; biny < numBinsY; biny++) {
                        final float[] yker = buildKernel(binHeight, numBinsY, biny, gaussianWindowSize);

                        for (int binx = 0; binx < numBinsX; binx++) {
                                final float[] xker = buildKernel(binWidth, numBinsX, binx, gaussianWindowSize);

                                for (int bint = 0; bint < numOriBins; bint++) {
                                        final FImage conv = data.gradientMagnitudes[bint].process(new FImageConvolveSeparable(xker, yker));
                                        final float[][] src = conv.pixels;

                                        final int descriptorOffset = bint + binx * numOriBins + biny * (numBinsX * numOriBins);
                                        int descriptorIndex = 0;

                                        for (int framey = data.boundMinY; framey <= data.boundMaxY - frameSizeY + 1; framey += stepY) {
                                                for (int framex = data.boundMinX; framex <= data.boundMaxX - frameSizeX + 1; framex += stepX) {
                                                        descriptors[descriptorIndex][descriptorOffset] = src[framey + biny * binHeight][framex
                                                                        + binx * binWidth];
                                                        descriptorIndex++;
                                                }
                                        }
                                }
                        }
                }
        }

        @Override
        public void analyseImage(FImage image, Rectangle bounds) {
                if (data == null)
                        data = new WorkingData();

                data.boundMinX = (int) bounds.x;
                data.boundMaxX = (int) (bounds.width - 1);
                data.boundMinY = (int) bounds.y;
                data.boundMaxY = (int) (bounds.height - 1);

                data.setupWorkingSpace(image, this);

                FImageGradients.gradientMagnitudesAndQuantisedOrientations(image, data.gradientMagnitudes);

                extractFeatures();

                normaliseDescriptors();
        }

        private void normaliseDescriptors() {
                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;
                final float energyNorm = frameSizeX * frameSizeY;

                for (int j = 0; j < descriptors.length; j++) {
                        final float[] arr = descriptors[j];

                        energies[j] = ArrayUtils.sumValues(arr) / energyNorm;

                        ArrayUtils.normalise(arr);

                        boolean changed = false;
                        for (int i = 0; i < arr.length; i++) {
                                if (arr[i] > valueThreshold) {
                                        arr[i] = valueThreshold;
                                        changed = true;
                                }
                        }

                        if (changed)
                                ArrayUtils.normalise(arr);
                }
        }

        @Override
        public LocalFeatureList<FloatDSIFTKeypoint> getFloatKeypoints() {
                final MemoryLocalFeatureList<FloatDSIFTKeypoint> keys = new MemoryLocalFeatureList<FloatDSIFTKeypoint>(numOriBins
                                * numBinsX * numBinsY, descriptors.length);

                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;

                final float deltaCenterX = 0.5F * binWidth * (numBinsX - 1);
                final float deltaCenterY = 0.5F * binHeight * (numBinsY - 1);

                for (int framey = data.boundMinY, i = 0; framey <= data.boundMaxY - frameSizeY + 1; framey += stepY) {
                        for (int framex = data.boundMinX; framex <= data.boundMaxX - frameSizeX + 1; framex += stepX, i++) {
                                keys.add(new FloatDSIFTKeypoint(framex + deltaCenterX, framey + deltaCenterY, descriptors[i],
                                                energies[i]));
                        }
                }

                return keys;
        }

        @Override
        public LocalFeatureList<ByteDSIFTKeypoint> getByteKeypoints() {
                final MemoryLocalFeatureList<ByteDSIFTKeypoint> keys = new MemoryLocalFeatureList<ByteDSIFTKeypoint>(numOriBins
                                * numBinsX * numBinsY, descriptors.length);

                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;

                final float deltaCenterX = 0.5F * binWidth * (numBinsX - 1);
                final float deltaCenterY = 0.5F * binHeight * (numBinsY - 1);

                for (int framey = data.boundMinY, i = 0; framey <= data.boundMaxY - frameSizeY + 1; framey += stepY) {
                        for (int framex = data.boundMinX; framex <= data.boundMaxX - frameSizeX + 1; framex += stepX, i++) {
                                keys.add(
                                                new ByteDSIFTKeypoint(framex + deltaCenterX, framey + deltaCenterY, descriptors[i], energies[i]));
                        }
                }

                return keys;
        }

        @Override
        public LocalFeatureList<FloatDSIFTKeypoint> getFloatKeypoints(float energyThreshold) {
                final MemoryLocalFeatureList<FloatDSIFTKeypoint> keys = new MemoryLocalFeatureList<FloatDSIFTKeypoint>(numOriBins
                                * numBinsX * numBinsY);

                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;

                final float deltaCenterX = 0.5F * binWidth * (numBinsX - 1);
                final float deltaCenterY = 0.5F * binHeight * (numBinsY - 1);

                for (int framey = data.boundMinY, i = 0; framey <= data.boundMaxY - frameSizeY + 1; framey += stepY) {
                        for (int framex = data.boundMinX; framex <= data.boundMaxX - frameSizeX + 1; framex += stepX, i++) {
                                if (energies[i] >= energyThreshold)
                                        keys.add(new FloatDSIFTKeypoint(framex + deltaCenterX, framey + deltaCenterY, descriptors[i],
                                                        energies[i]));
                        }
                }

                return keys;
        }

        @Override
        public LocalFeatureList<ByteDSIFTKeypoint> getByteKeypoints(float energyThreshold) {
                final MemoryLocalFeatureList<ByteDSIFTKeypoint> keys = new MemoryLocalFeatureList<ByteDSIFTKeypoint>(numOriBins
                                * numBinsX * numBinsY);

                final int frameSizeX = binWidth * (numBinsX - 1) + 1;
                final int frameSizeY = binHeight * (numBinsY - 1) + 1;

                final float deltaCenterX = 0.5F * binWidth * (numBinsX - 1);
                final float deltaCenterY = 0.5F * binHeight * (numBinsY - 1);

                for (int framey = data.boundMinY, i = 0; framey <= data.boundMaxY - frameSizeY + 1; framey += stepY) {
                        for (int framex = data.boundMinX; framex <= data.boundMaxX - frameSizeX + 1; framex += stepX, i++) {
                                if (energies[i] >= energyThreshold)
                                        keys.add(new ByteDSIFTKeypoint(framex + deltaCenterX, framey + deltaCenterY, descriptors[i],
                                                        energies[i]));
                        }
                }

                return keys;
        }

        /**
         * Get the computed raw dense SIFT descriptors from the previous call to
         * {@link #analyseImage(FImage)} or {@link #analyseImage(FImage, Rectangle)}
         * .
         *
         * @return the descriptors.
         */
        @Override
        public float[][] getDescriptors() {
                return descriptors;
        }

        @Override
        public DenseSIFT clone() {
                final DenseSIFT clone = (DenseSIFT) super.clone();

                clone.descriptors = null;
                clone.energies = null;
                clone.data = null;

                return clone;
        }

        @Override
        public void setBinWidth(int size) {
                this.binWidth = size;
        }

        @Override
        public void setBinHeight(int size) {
                this.binHeight = size;
        }

        @Override
        public int getBinWidth() {
                return binWidth;
        }

        @Override
        public int getBinHeight() {
                return binHeight;
        }

        @Override
        public int getNumBinsX() {
                return numBinsX;
        }

        @Override
        public int getNumBinsY() {
                return numBinsY;
        }

        @Override
        public int getNumOriBins() {
                return numOriBins;
        }
}