package ij.process;
import ij.*;
import ij.plugin.FFT;
import ij.plugin.ContrastEnhancer;
import java.awt.image.ColorModel; 

/**
This class contains a Java implementation of the Fast Hartley
Transform. It is based on Pascal code in NIH Image contributed 
by Arlo Reeves (http://imagej.nih.gov/ij/docs/ImageFFT/).
The Fast Hartley Transform was restricted by U.S. Patent No. 4,646,256, 
but was placed in the public domain by Stanford University in 1995 
and is now freely available.
*/
public class FHT extends FloatProcessor {
	private boolean isFrequencyDomain;
	private int maxN;
	private float[] C;
	private float[] S;
	private int[] bitrev;
	private float[] tempArr;
	private boolean showProgress = true;
	
	/** Used by the FFT class. */
	public boolean quadrantSwapNeeded;
	/** Used by the FFT class. */
	public ColorProcessor rgb;
	/** Used by the FFT class. */
	public int originalWidth;
	/** Used by the FFT class. */
	public int originalHeight;
	/** Used by the FFT class. */
	public int originalBitDepth;
	/** Used by the FFT class. */
	public ColorModel originalColorModel;

	/** Constructs a FHT object from an ImageProcessor. Byte, short and RGB images 
		are converted to float. Float images are duplicated. */
	public FHT(ImageProcessor ip) {
		this(ip, false);
	}

	public FHT(ImageProcessor ip, boolean isFrequencyDomain) {
		super(ip.getWidth(), ip.getHeight(), (float[])((ip instanceof FloatProcessor)?ip.duplicate().getPixels():ip.convertToFloat().getPixels()), null);
		this.isFrequencyDomain = isFrequencyDomain;
		maxN = getWidth();
		resetRoi();
	}

	public FHT() {
		super(8,8); //create dummy FloatProcessor
	}

	/** Returns true of this FHT contains a square image with a width that is a power of two. */
	public boolean powerOf2Size() {
		int i=2;
		while(i0, use (x[i]+x[maxN-i])/(2*maxN).
	 *  The imaginary part of the complex FFT, with i>0, is given by (x[i]-x[maxN-i])/(2*maxN)
	 *  The coefficients of cosine and sine are like the real and imaginary values above,
	 *  but you have to divide by maxN instead of 2*maxN.
	 */
	public void transform1D(float[] x) {
		int n = x.length;
		if (S==null || n!=maxN) {
			if (!isPowerOf2(n))
				throw new IllegalArgumentException("Not power of 2 length: "+n);
			initializeTables(n);
		}
		dfht3(x, 0, false, n);
	}

    /** Performs an inverse 1D Fast Hartley Transform (FHT) of an array */
	public void inverseTransform1D(float[] fht) {
		int n = fht.length;
		if (S==null || n!=maxN) {
			if (!isPowerOf2(n))
				throw new IllegalArgumentException("Not power of 2 length: "+n);
			initializeTables(n);
		}
		dfht3(fht, 0, true, n);
	}

	void transform(boolean inverse) {
		//IJ.log("transform: "+maxN+" "+inverse);
		if (!powerOf2Size())
			throw new  IllegalArgumentException("Image not power of 2 size or not square: "+width+"x"+height);
		maxN = width;
		if (S==null)
			initializeTables(maxN);
		float[] fht = (float[])getPixels();
	 	rc2DFHT(fht, inverse, maxN);
		isFrequencyDomain = !inverse;
	}
	
	void initializeTables(int maxN) {
	    if (maxN>0x40000000)
	        throw new  IllegalArgumentException("Too large for FHT:  "+maxN+" >2^30");
		makeSinCosTables(maxN);
		makeBitReverseTable(maxN);
		tempArr = new float[maxN];
	}

	void makeSinCosTables(int maxN) {
		int n = maxN/4;
		C = new float[n];
		S = new float[n];
		double theta = 0.0;
		double dTheta = 2.0 * Math.PI/maxN;
		for (int i=0; i 2) {
			 // third + stages computed here
			gpSize = 4;
			numBfs = 2;
			numGps = numGps / 2;
			for (stage=2; stage 2 */

		if (inverse)  {
			for (i=0; imax) max = r;
			}
		}

		max = (float)Math.log(max);
		min = (float)Math.log(min);
		if (Float.isNaN(min) || max-min>50)
			min = max - 50; //display range not more than approx e^50
		scale = (float)(253.999/(max-min));

		for (int row=0; row
		    2 1
		    3 4