import ij.*;
import ij.plugin.filter.PlugInFilter;
import ij.process.*;
import java.awt.*;

/* Compute the 3D Fast Harley Transform of a 32-bit image stack.
	The width and height of the images must be powers of 2.  The need not match.
	The number of slices is not restricted, but a slow transform
	technique is used if the number of slices is not a power of 2.
	The normalization is such that the forward and inverse transforms are
	performed by the same plugin.

	Bob Dougherty
	Version 0 4/24/2005.
	Version 0.1 4/25/2005  Explicitly redraw the image upon completion.

	Much is the code is from FFT_Filter by Joachim Walter.  Some items are from
	the related FHT by Wayne Rasband.  The following notice is from the FHT source code in ImageJ:

		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://rsb.info.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_3D implements  PlugInFilter {
	private ImagePlus imp;

	public int setup(String arg, ImagePlus imp) {
		IJ.register(FHT_3D.class);
 		this.imp = imp;
		return DOES_32;
	}

	public void run(ImageProcessor ip) {
		ImageStack stack = imp.getStack();
		int w = stack.getWidth();
		int h = stack.getHeight();
		if(!(powerOf2Size(w)&&powerOf2Size(h))){
			IJ.error("The width and height must be powers of 2.");
			return;
		}
		int d = imp.getStackSize();
		float[][] data = new float[d][];
		for (int i = 0; i < d; i++){
 			data[i] = (float[])stack.getProcessor(i+1).getPixels();
		}
		FHT3D(data,w,h,d);
		imp.updateAndDraw();
	}
    boolean powerOf2Size(int w) {
        int i=2;
        while(i<w) i *= 2;
        return i==w;
    }
	public void FHT3D(float[][] data,int w, int h, int d) {
		float[] sw = new float[w/4];
		float[] cw = new float[w/4];
		float[] sh = new float[h/4];
		float[] ch = new float[h/4];
		makeSinCosTables(w,sw,cw);
		makeSinCosTables(h,sh,ch);
		for (int i = 0; i < d; i++){
	 		rc2DFHT(data[i], w, h, sw, cw, sh, ch);
		}
		float[] u = new float[d];
		//if(IJ.getNumber("0 for fast, 1 for slow",0)==0){
		if(powerOf2Size(d)){
			float[] s = new float[d/4];
			float[] c = new float[d/4];
			makeSinCosTables(d,s,c);
			for(int k2 = 0; k2 < h; k2++){
				for(int k1 = 0; k1 < w; k1++){
					int ind = k1 + k2*w;
					for(int k3 = 0; k3 < d; k3++){
						u[k3] = data[k3][ind];
					}
					dfht3(u, 0, d, s, c);
					for(int k3 = 0; k3 < d; k3++){
						data[k3][ind] = u[k3];
					}
				}
			}
		}else{
			float[] cas = hartleyCoefs(d);
			float[] work = new float[d];
			for(int k2 = 0; k2 < h; k2++){
				for(int k1 = 0; k1 < w; k1++){
					int ind = k1 + k2*w;
					for(int k3 = 0; k3 < d; k3++){
						u[k3] = data[k3][ind];
					}
					slowHT(u,cas,d,work);
					for(int k3 = 0; k3 < d; k3++){
						data[k3][ind] = u[k3];
					}
				}
			}
		}
		//Convert to actual Hartley transform
		float A,B,C,D,E,F,G,H;
		int k1C,k2C,k3C;
		for(int k3 = 0; k3 <= d/2; k3++){
			k3C = (d - k3) % d;
			for(int k2 = 0; k2 <= h/2; k2++){
				k2C = (h - k2) % h;
				for (int k1 = 0; k1 <= w/2; k1++){
					k1C = (w - k1) % w;
					A = data[k3][k1 + w*k2C];
					B = data[k3][k1C + w*k2];
					C = data[k3C][k1 + w*k2];
					D = data[k3C][k1C + w*k2C];
					E = data[k3C][k1 + w*k2C];
					F = data[k3C][k1C + w*k2];
					G = data[k3][k1 + w*k2];
					H = data[k3][k1C + w*k2C];
					data[k3][k1 + w*k2] = (A+B+C-D)/2;
					data[k3C][k1 + w*k2] = (E+F+G-H)/2;
					data[k3][k1 + w*k2C] = (G+H+E-F)/2;
					data[k3C][k1 + w*k2C] = (C+D+A-B)/2;
					data[k3][k1C + w*k2] = (H+G+F-E)/2;
					data[k3C][k1C + w*k2] = (D+C+B-A)/2;
					data[k3][k1C + w*k2C] = (B+A+D-C)/2;
					data[k3C][k1C + w*k2C] = (F+E+H-G)/2;
				}
			}
		}
		//normalize
		float norm = (float)Math.sqrt(d*h*w);
		for(int k3 = 0; k3 < d; k3++){
			for(int k2 = 0; k2 < h; k2++){
				for (int k1 = 0; k1 < w; k1++){
					data[k3][k1 + w*k2] /= norm;
				}
			}
		}
	}
	float[] hartleyCoefs(int max){
		float[] cas = new float[max*max];
		int ind = 0;
		for(int n = 0; n < max; n++){
			for (int k = 0; k < max; k++){
				double arg = (2*Math.PI*k*n)/max;
				cas[ind++] = (float)(Math.cos(arg) + Math.sin(arg));
			}
		}
		return cas;
	}
	void slowHT(float[] u, float[] cas, int max, float[] work){
		int ind = 0;
		for(int k = 0; k < max; k++){
			float sum = 0;
			for(int n = 0; n < max; n++){
				sum += u[n]*cas[ind++];
			}
			work[k] = sum;
		}
		for (int k = 0; k < max; k++){
			u[k] = work[k];
		}
	}
	void makeSinCosTables(int maxN, float[] s, float[] c) {
		int n = maxN/4;
		double theta = 0.0;
		double dTheta = 2.0 * Math.PI/maxN;
		for (int i=0; i<n; i++) {
			c[i] = (float)Math.cos(theta);
			s[i] = (float)Math.sin(theta);
			theta += dTheta;
		}
	}
	/** Row-column Fast Hartley Transform */
	void rc2DFHT(float[] x, int w, int h, float[] sw, float[] cw, float[] sh, float[] ch) {
		for (int row=0; row<h; row++)
			dfht3(x, row*w, w, sw, cw);
		float[] temp = new float[h];
		for(int col = 0; col < w; col++){
			for (int row = 0; row < h; row++){
				temp[row] = x[col + w*row];
			}
			dfht3(temp, 0, h, sh, ch);
			for (int row = 0; row < h; row++){
				x[col + w*row] = temp[row];
			}
		}
	}
	/* An optimized real FHT */
	void dfht3 (float[] x, int base, int maxN, float[] s, float[] c) {
		int i, stage, gpNum, gpIndex, gpSize, numGps, Nlog2;
		int bfNum, numBfs;
		int Ad0, Ad1, Ad2, Ad3, Ad4, CSAd;
		float rt1, rt2, rt3, rt4;

		Nlog2 = log2(maxN);
		BitRevRArr(x, base, Nlog2, maxN);	//bitReverse the input array
		gpSize = 2;     //first & second stages - do radix 4 butterflies once thru
		numGps = maxN / 4;
		for (gpNum=0; gpNum<numGps; gpNum++)  {
			Ad1 = gpNum * 4;
			Ad2 = Ad1 + 1;
			Ad3 = Ad1 + gpSize;
			Ad4 = Ad2 + gpSize;
			rt1 = x[base+Ad1] + x[base+Ad2];   // a + b
			rt2 = x[base+Ad1] - x[base+Ad2];   // a - b
			rt3 = x[base+Ad3] + x[base+Ad4];   // c + d
			rt4 = x[base+Ad3] - x[base+Ad4];   // c - d
			x[base+Ad1] = rt1 + rt3;      // a + b + (c + d)
			x[base+Ad2] = rt2 + rt4;      // a - b + (c - d)
			x[base+Ad3] = rt1 - rt3;      // a + b - (c + d)
			x[base+Ad4] = rt2 - rt4;      // a - b - (c - d)
		}
		if (Nlog2 > 2) {
			 // third + stages computed here
			gpSize = 4;
			numBfs = 2;
			numGps = numGps / 2;
			//IJ.write("FFT: dfht3 "+Nlog2+" "+numGps+" "+numBfs);
			for (stage=2; stage<Nlog2; stage++) {
				for (gpNum=0; gpNum<numGps; gpNum++) {
					Ad0 = gpNum * gpSize * 2;
					Ad1 = Ad0;     // 1st butterfly is different from others - no mults needed
					Ad2 = Ad1 + gpSize;
					Ad3 = Ad1 + gpSize / 2;
					Ad4 = Ad3 + gpSize;
					rt1 = x[base+Ad1];
					x[base+Ad1] = x[base+Ad1] + x[base+Ad2];
					x[base+Ad2] = rt1 - x[base+Ad2];
					rt1 = x[base+Ad3];
					x[base+Ad3] = x[base+Ad3] + x[base+Ad4];
					x[base+Ad4] = rt1 - x[base+Ad4];
					for (bfNum=1; bfNum<numBfs; bfNum++) {
					// subsequent BF's dealt with together
						Ad1 = bfNum + Ad0;
						Ad2 = Ad1 + gpSize;
						Ad3 = gpSize - bfNum + Ad0;
						Ad4 = Ad3 + gpSize;

						CSAd = bfNum * numGps;
						rt1 = x[base+Ad2] * c[CSAd] + x[base+Ad4] * s[CSAd];
						rt2 = x[base+Ad4] * c[CSAd] - x[base+Ad2] * s[CSAd];

						x[base+Ad2] = x[base+Ad1] - rt1;
						x[base+Ad1] = x[base+Ad1] + rt1;
						x[base+Ad4] = x[base+Ad3] + rt2;
						x[base+Ad3] = x[base+Ad3] - rt2;

					} /* end bfNum loop */
				} /* end gpNum loop */
				gpSize *= 2;
				numBfs *= 2;
				numGps = numGps / 2;
			} /* end for all stages */
		} /* end if Nlog2 > 2 */
	}
	int log2 (int x) {
		int count = 15;
		while (!btst(x, count))
			count--;
		return count;
	}
	private boolean btst (int  x, int bit) {
		//int mask = 1;
		return ((x & (1<<bit)) != 0);
	}
	void BitRevRArr (float[] x, int base, int bitlen, int maxN) {
		int    l;
		float[] tempArr = new float[maxN];
		for (int i=0; i<maxN; i++)  {
			l = BitRevX (i, bitlen);  //i=1, l=32767, bitlen=15
			tempArr[i] = x[base+l];
		}
		for (int i=0; i<maxN; i++)
			x[base+i] = tempArr[i];
	}
	private int BitRevX (int  x, int bitlen) {
		int  temp = 0;
		for (int i=0; i<=bitlen; i++)
			if ((x & (1<<i)) !=0)
				temp  |= (1<<(bitlen-i-1));
		return temp & 0x0000ffff;
	}
	private int bset (int x, int bit) {
		x |= (1<<bit);
		return x;
	}
}