/**
    * 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.global;
  
  import org.openimaj.citation.annotation.Reference;
  import org.openimaj.citation.annotation.ReferenceType;
  import org.openimaj.feature.FloatFV;
  import org.openimaj.image.FImage;
  import org.openimaj.image.Image;
  import org.openimaj.image.MBFImage;
  import org.openimaj.image.analyser.ImageAnalyser;
  import org.openimaj.image.colour.ColourMap;
  import org.openimaj.image.colour.ColourSpace;
  import org.openimaj.image.colour.RGBColour;
  import org.openimaj.image.processing.algorithm.FourierTransform;
  import org.openimaj.image.processing.convolution.FourierConvolve;
  import org.openimaj.image.processing.convolution.GaborFilters;
  import org.openimaj.image.processing.resize.ResizeProcessor;
  import org.openimaj.image.processor.SinglebandImageProcessor;
  
  import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D;
  
  /**
   * Implementation of the "Gist" spatial envelope feature. Based on (and tested
   * against) the original Matlab implementation from <a
   * href="http://people.csail.mit.edu/torralba/code/spatialenvelope/"
   * >http://people.csail.mit.edu/torralba/code/spatialenvelope/</a>, and designed
   * to produce comparable features.
   * <p>
   * The Gist or Spatial Envelope is a very low dimensional representation of the
   * scene. Gist encodes a set of perceptual dimensions (naturalness, openness,
   * roughness, expansion, ruggedness) that represents the dominant spatial
   * structure of a scene. These perceptual dimensions are reliably estimated
   * using coarsely localized spectral information from the image.
   * <p>
   * <b>Implementation notes:</b> This class supports both fixed and variable size
   * Gist. In the fixed size mode, the image is resized and centre-cropped to a
   * pre-defined size (allowing the Gabor filters to be pre-computed and reused).
   * In the variable size mode the Gabor filters will be re-computed as necessary
   * if the image size changes. For computing features to compare images, fixed
   * size mode makes more sense.
   * <p>
   * <b>Example usage for image comparison:</b><br>
   * </br> <code>
   * <pre>
   * FImage i1, i2 = ...
   * Gist<FImage> fsg = new Gist<FImage>(256, 256);
   *
   * fsg.analyseImage(i1);
   * FloatFV f1 = fsg.getResponse();
   *
   * fsg.analyseImage(i2);
   * FloatFV f2 = fsg.getResponse();
   *
   * double d = FloatFVComparison.SUM_SQUARE.compare(f1, f2);
   * </pre>
   * </code>
   * <p>
   * The ability to generate a visualisation of the current descriptor in the same
   * way as provided by the original Matlab code is given by the
   * {@link #visualiseDescriptor(int)} method.
   *
   * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
   *
   * @param <IMAGE>
   *            The type of image. {@link MBFImage} and {@link FImage} are
   *            supported; other types might not work.
   */
  @Reference(
  		type = ReferenceType.Article,
  		author = { "Oliva, Aude", "Torralba, Antonio" },
  		title = "Modeling the Shape of the Scene: A Holistic Representation of the Spatial Envelope",
 		year = "2001",
 		journal = "Int. J. Comput. Vision",
 		pages = { "145", "", "175" },
 		url = "http://dx.doi.org/10.1023/A:1011139631724",
 		month = "May",
 		number = "3",
 		publisher = "Kluwer Academic Publishers",
 		volume = "42",
 		customData = {
 				"issn", "0920-5691",
 				"numpages", "31",
 				"doi", "10.1023/A:1011139631724",
 				"acmid", "598462",
 				"address", "Hingham, MA, USA",
 				"keywords", "energy spectrum, natural images, principal components, scene recognition, spatial layout"
 		})
 public class Gist<IMAGE extends Image<?, IMAGE> & SinglebandImageProcessor.Processable<Float, FImage, IMAGE>>
 implements
 ImageAnalyser<IMAGE>
 {
 	/**
 	 * The default number of filter orientations per scale (from HF to LF)
 	 */
 	public static final int[] DEFAULT_ORIENTATIONS_PER_SCALE = { 8, 8, 8, 8 };
 
 	/**
 	 * The default number spatial blocks
 	 */
 	public static final int DEFAULT_NUMBER_OF_BLOCKS = 4;
 
 	/**
 	 * The default number of cycles per image for the pre-filter Gaussian
 	 */
 	public static final int DEFAULT_PREFILTER_FC = 4;
 
 	/**
 	 * The default amount of padding to apply before convolving with the Gabor
 	 * filters
 	 */
 	public static final int DEFAULT_BOUNDARY_EXTENSION = 32;
 
 	/**
 	 * The number of blocks in each direction
 	 */
 	public int numberOfBlocks = DEFAULT_NUMBER_OF_BLOCKS;
 
 	/**
 	 * The number of cycles per image for the pre-filter Gaussian
 	 */
 	public int prefilterFC = DEFAULT_PREFILTER_FC;
 
 	/**
 	 * The amount of padding to add before convolving with the Gabor functions
 	 */
 	public int boundaryExtension = DEFAULT_BOUNDARY_EXTENSION;
 
 	/**
 	 * Default image size (both height and width) for fixed size Gist
 	 */
 	public static final int DEFAULT_SIZE = 128;
 
 	protected FImage[] gaborFilters;
 	protected FloatFV response;
 	protected int[] orientationsPerScale;
 	protected boolean fixedSize;
 	protected int imageWidth;
 	protected int imageHeight;
 
 	/**
 	 * Construct a fixed size Gist extractor using the default values.
 	 */
 	public Gist() {
 		this(true);
 	}
 
 	/**
 	 * Construct a variable or fixed size Gist extractor using the default
 	 * values.
 	 *
 	 * @param fixedSize
 	 *            if true the images should be resized and cropped
 	 */
 	public Gist(boolean fixedSize) {
 		this(DEFAULT_SIZE, DEFAULT_SIZE, DEFAULT_ORIENTATIONS_PER_SCALE, fixedSize);
 	}
 
 	/**
 	 * Construct a fixed-size Gist extractor, with the given image size and the
 	 * default number of orientations per scale for the Gabor jet.
 	 *
 	 * @param width
 	 *            the image width
 	 * @param height
 	 *            the image height
 	 */
 	public Gist(int width, int height) {
 		this(width, height, DEFAULT_ORIENTATIONS_PER_SCALE);
 	}
 
 	/**
 	 * Construct a fixed-size Gist extractor, with the given image size and
 	 * number of orientations per scale for the Gabor jet.
 	 *
 	 * @param width
 	 *            the image width
 	 * @param height
 	 *            the image height
 	 * @param orientationsPerScale
 	 *            the number of Gabor orientations per scale (from HF to LF)
 	 */
 	public Gist(int width, int height, int[] orientationsPerScale) {
 		this(width, height, orientationsPerScale, true);
 	}
 
 	/**
 	 * Construct a variable or fixed size Gist extractor with the given number
 	 * of orientations per scale for the Gabor jet.
 	 *
 	 * @param orientationsPerScale
 	 *            the number of Gabor orientations per scale (from HF to LF)
 	 * @param fixedSize
 	 *            if true the images should be resized and cropped
 	 */
 	public Gist(int[] orientationsPerScale, boolean fixedSize) {
 		this(DEFAULT_SIZE, DEFAULT_SIZE, orientationsPerScale, fixedSize);
 	}
 
 	/**
 	 * Construct a variable or fixed size Gist extractor with the given number
 	 * of orientations per scale for the Gabor jet. The given height and width
 	 * are ignored in variable size mode.
 	 *
 	 * @param width
 	 *            the image width
 	 * @param height
 	 *            the image height
 	 * @param orientationsPerScale
 	 *            the number of Gabor orientations per scale (from HF to LF)
 	 * @param fixedSize
 	 *            if true the images should be resized and cropped
 	 */
 	public Gist(int width, int height, int[] orientationsPerScale, boolean fixedSize) {
 		this.fixedSize = fixedSize;
 		this.orientationsPerScale = orientationsPerScale;
 		this.imageWidth = width;
 		this.imageHeight = height;
 
 		if (fixedSize)
 			this.gaborFilters = GaborFilters.createGaborJets(width + 2 * this.boundaryExtension, height + 2
 					* this.boundaryExtension, orientationsPerScale);
 	}
 
 	@Override
 	public void analyseImage(IMAGE image) {
 		if (fixedSize) {
 			final double sc = Math.max((double) imageWidth / (double) image.getWidth(), (double) imageHeight
 					/ (double) image.getHeight());
 
 			final IMAGE resized = image.process(new ResizeProcessor((float) sc));
 			final IMAGE roi = resized.extractCenter(imageWidth, imageHeight);
 			extractGist(roi);
 		} else {
 			if (gaborFilters == null || gaborFilters[0].width != image.getWidth()
 					|| gaborFilters[0].height != image.getHeight())
 			{
 				gaborFilters = GaborFilters.createGaborJets(image.getWidth() + 2 * this.boundaryExtension,
 						image.getHeight() + 2 * this.boundaryExtension, orientationsPerScale);
 			}
 
 			extractGist(image.clone()); // clone to stop side effects from
 			// normalisation further down
 		}
 	}
 
 	protected Gist(int[] orientationsPerScale) {
 		this.orientationsPerScale = orientationsPerScale;
 	}
 
 	protected void extractGist(IMAGE image) {
 		MBFImage mbfimage;
 		if (image instanceof FImage) {
 			mbfimage = new MBFImage((FImage) image);
 		} else if (image instanceof MBFImage) {
 			mbfimage = (MBFImage) image;
 		} else {
 			throw new UnsupportedOperationException("Image type " + image.getClass()
 					+ " is not currently supported. Please file a bug report.");
 		}
 
 		final MBFImage o = prefilter(mbfimage.normalise());
 		this.response = gistGabor(o);
 	}
 
 	private MBFImage prefilter(MBFImage img) {
 		final int w = 5;
 		final double s1 = this.prefilterFC / Math.sqrt(Math.log(2));
 
 		final int sw = img.getWidth() + 2 * w;
 		final int sh = img.getHeight() + 2 * w;
 		int n = Math.max(sw, sh);
 		n = n + n % 2;
 		img = img.paddingSymmetric(w, w, w + n - sw, w + n - sh);
 
 		final FImage filter = new FImage(2 * n, n);
 		for (int j = 0; j < n; j++) {
 			final int fy = j - n / 2;
 
 			for (int i = 0; i < n * 2; i += 2) {
 				final int fx = (i / 2) - n / 2;
 
 				filter.pixels[j][i] = (float) Math.exp(-(fx * fx + fy * fy) / (s1 * s1));
 			}
 		}
 
 		final MBFImage output = new MBFImage();
 		for (int b = 0; b < img.numBands(); b++) {
 			final FImage band = img.getBand(b);
 			for (int y = 0; y < band.height; y++) {
 				for (int x = 0; x < band.width; x++) {
 					band.pixels[y][x] = (float) Math.log(1 + band.pixels[y][x] * 255);
 				}
 			}
 
 			output.bands.add(band.subtractInplace(FourierConvolve.convolvePrepared(band, filter, true)));
 		}
 		final FImage mean = output.flatten();
 		final FImage meansq = mean.multiply(mean);
 		final FImage localstd = FourierConvolve.convolvePrepared(meansq, filter, true);
 
 		for (int b = 0; b < img.numBands(); b++) {
 			final FImage band = output.getBand(b);
 			for (int y = 0; y < localstd.height; y++)
 				for (int x = 0; x < localstd.width; x++)
 					band.pixels[y][x] = (float) (band.pixels[y][x] / (0.2 + Math
 							.sqrt(Math.abs(localstd.pixels[y][x]))));
 		}
 
 		return output.extractROI(w, w, sw - w - w, sh - w - w);
 	}
 
 	private FloatFV gistGabor(MBFImage img) {
 		final int blocksPerFilter = computeNumberOfSamplingBlocks();
 		final int nFeaturesPerBand = gaborFilters.length * blocksPerFilter;
 		final int nFilters = this.gaborFilters.length;
 
 		// pad the image
 		img = img.paddingSymmetric(boundaryExtension, boundaryExtension, boundaryExtension, boundaryExtension);
 
 		final int cols = img.getCols();
 		final int rows = img.getRows();
 		final FloatFFT_2D fft = new FloatFFT_2D(rows, cols);
 
 		final float[][] workingSpace = new float[rows][cols * 2];
 		final FloatFV fv = new FloatFV(nFeaturesPerBand * img.numBands());
 
 		for (int b = 0; b < img.numBands(); b++) {
 			final FImage band = img.bands.get(b);
 
 			final float[][] preparedImage =
 					FourierTransform.prepareData(band.pixels, rows, cols, true);
 			fft.complexForward(preparedImage);
 
 			for (int i = 0; i < nFilters; i++) {
 				// convolve with the filter
 				FImage ig = performConv(fft, preparedImage, workingSpace, this.gaborFilters[i], rows, cols);
 
 				// remove padding
 				ig = ig.extractROI(boundaryExtension, boundaryExtension, band.width - 2 * boundaryExtension, band.height
 						- 2
 						* boundaryExtension);
 
 				sampleResponses(ig, fv.values, b * nFeaturesPerBand + i * blocksPerFilter);
 			}
 		}
 
 		return fv;
 	}
 
 	/**
 	 * Compute the number of sampling blocks that are used for every filter. The
 	 * default implementation returns {@link #numberOfBlocks}*
 	 * {@link #numberOfBlocks}, but can be overridden in combination with
 	 * {@link #sampleResponses(FImage, float[], int)} in a subclass to support
 	 * different spatial sampling strategies.
 	 *
 	 * @return the number of sampling blocks per filter.
 	 */
 	protected int computeNumberOfSamplingBlocks() {
 		return numberOfBlocks * numberOfBlocks;
 	}
 
 	/**
 	 * Sample the average response from each of the blocks in the image and
 	 * insert into the vector. This method could be overridden to support
 	 * different spatial aggregation strategies (in which case
 	 * {@link #computeNumberOfSamplingBlocks()} should also be overridden).
 	 *
 	 * @param image
 	 *            the image to sample
 	 * @param v
 	 *            the vector to write into
 	 * @param offset
 	 *            the offset from which to sample
 	 */
 	protected void sampleResponses(FImage image, float[] v, int offset) {
 		final int gridWidth = image.width / this.numberOfBlocks;
 		final int gridHeight = image.height / this.numberOfBlocks;
 
 		for (int iy = 0; iy < this.numberOfBlocks; iy++) {
 			final int starty = gridHeight * iy;
 			final int stopy = Math.min(starty + gridHeight, image.height);
 
 			for (int ix = 0; ix < this.numberOfBlocks; ix++) {
 				final int startx = gridWidth * ix;
 				final int stopx = Math.min(startx + gridWidth, image.width);
 
 				float avg = 0;
 				for (int y = starty; y < stopy; y++) {
 					for (int x = startx; x < stopx; x++) {
 						avg += image.pixels[y][x];
 					}
 				}
 				avg /= ((stopx - startx) * (stopy - starty));
 
 				// note y and x transposed to conform to the matlab
 				// implementation
 				v[offset + iy + ix * this.numberOfBlocks] = avg;
 			}
 		}
 	}
 
 	/*
 	 * Perform convolution in the frequency domain an reconstruct the resultant
 	 * image as the magnitudes of the complex components from the ifft.
 	 */
 	private FImage performConv(FloatFFT_2D fft, float[][] preparedImage, float[][] workingSpace, FImage filterfft,
 			int rows, int cols)
 	{
 		final float[][] preparedKernel = filterfft.pixels;
 
 		for (int y = 0; y < rows; y++) {
 			for (int x = 0; x < cols; x++) {
 				final float reImage = preparedImage[y][x * 2];
 				final float imImage = preparedImage[y][1 + x * 2];
 
 				final float reKernel = preparedKernel[y][x * 2];
 				final float imKernel = preparedKernel[y][1 + x * 2];
 
 				final float re = reImage * reKernel - imImage * imKernel;
 				final float im = reImage * imKernel + imImage * reKernel;
 
 				workingSpace[y][x * 2] = re;
 				workingSpace[y][1 + x * 2] = im;
 			}
 		}
 
 		fft.complexInverse(workingSpace, true);
 
 		final FImage out = new FImage(cols, rows);
 		for (int r = 0; r < rows; r++) {
 			for (int c = 0; c < cols; c++) {
 				out.pixels[r][c] = (float) Math.sqrt(workingSpace[r][c * 2] * workingSpace[r][c * 2]
 						+ workingSpace[r][1 + c * 2]
 								* workingSpace[r][1 + c * 2]);
 			}
 		}
 		return out;
 	}
 
 	/**
 	 * Compute the descriptor visualisation in the same form as the original
 	 * matlab code. The resultant image illustrates the dominant filter
 	 * responses for each spatial bin. The filter scales are drawn with
 	 * differing hues, and brightness is proportional to response value.
 	 *
 	 * @param height
 	 *            the desired height of the produced image
 	 * @return the visualisation TODO: colour descriptors
 	 */
 	public MBFImage visualiseDescriptor(int height) {
 		final Float[][] C = ColourMap.HSV.generateColours(orientationsPerScale.length);
 
 		final FImage[] G = new FImage[this.gaborFilters.length];
 		for (int i = 0; i < gaborFilters.length; i++) {
 			G[i] = new FImage(this.gaborFilters[i].width / 2, this.gaborFilters[i].height);
 			FourierTransform.unprepareData(gaborFilters[i].pixels, G[i], false);
 			G[i] = ResizeProcessor.halfSize(G[i]);
 			G[i].addInplace(G[i].clone().flipY().flipX());
 		}
 
 		final int blocksPerFilter = computeNumberOfSamplingBlocks();
 		final int nFeaturesPerBand = gaborFilters.length * blocksPerFilter;
 		final int nBands = this.response.values.length / nFeaturesPerBand;
 
 		final MBFImage output = new MBFImage(nBands * height, height, ColourSpace.RGB);
 		output.fill(RGBColour.WHITE);
 
 		for (int b = 0; b < nBands; b++) {
 			float max = 0;
 			final MBFImage[] blockImages = new MBFImage[numberOfBlocks * numberOfBlocks];
 			for (int y = 0, k = 0; y < numberOfBlocks; y++) {
 				for (int x = 0; x < numberOfBlocks; x++, k++) {
 					blockImages[k] = new MBFImage(G[0].width, G[0].height, 3);
 
 					for (int s = 0, j = 0; s < orientationsPerScale.length; s++) {
 						for (int i = 0; i < orientationsPerScale[s]; i++, j++) {
 							final MBFImage col = new MBFImage(G[0].width, G[0].height, 3);
 							col.fill(C[s]).multiplyInplace(
 									this.response.values[y + x * numberOfBlocks + j * numberOfBlocks * numberOfBlocks]);
 
 							blockImages[k].addInplace(G[j].toRGB().multiply(col));
 						}
 					}
 
 					for (int i = 0; i < blockImages[k].numBands(); i++)
 						max = Math.max(max, blockImages[k].bands.get(i).max());
 				}
 			}
 
 			final int ts = (height / 4);
 			final ResizeProcessor rp = new ResizeProcessor(ts, ts, true);
 			for (int y = 0, k = 0; y < numberOfBlocks; y++) {
 				for (int x = 0; x < numberOfBlocks; x++, k++) {
 					blockImages[k].divideInplace(max);
 					final MBFImage tmp = blockImages[k].process(rp);
 					tmp.drawLine(0, 0, 0, ts - 1, 1, RGBColour.WHITE);
 					tmp.drawLine(0, 0, ts - 1, 0, 1, RGBColour.WHITE);
 					tmp.drawLine(ts - 1, 0, ts - 1, ts - 1, 1, RGBColour.WHITE);
 					tmp.drawLine(0, ts - 1, ts - 1, ts - 1, 1, RGBColour.WHITE);
 					output.drawImage(tmp, x * ts + b * height, y * ts);
 				}
 			}
 		}
 
 		return output;
 	}
 
 	/**
 	 * Get the response vector from the previous call to
 	 * {@link #analyseImage(Image)}.
 	 *
 	 * @return the response
 	 */
 	public FloatFV getResponse() {
 		return response;
 	}
 }