/*--- formatted by Jindent 2.1, (www.c-lab.de/~jindent) ---*/

/*
 * JIGL--Java Imaging and Graphics Library
 * Copyright (C)1999 Brigham Young University
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 * A copy of the GNU Library General Public Licence is contained in
 * /jigl/licence.txt
 */

package jm.jigl; // Modifizierte Paketstruktur [JM]
import jm.jigl.*;
import ij.*;


/**
 * Class declaration
 * Modifikation der urspruenglichen jigl Paketstruktur durch Jens Martin
 * Aenderungen im Sourcecode sind mit [JM] gekennzeichnet.
 * @version 1.3, Modifikationen der Methode doFFT(): Schnellere Ausfuehrung unter SUN JRE 1.1.x.
 * @see jigl.utils.signal.FFT
 */
public class FFT {

	static final int		MAX_POW = 31;
	static final int[]	pow2 = new int[MAX_POW + 1];

	// static initializer
	static {
		pow2[0] = 1;

		for (int i = 1; i <= MAX_POW; i++) {
			pow2[i] = 2 * pow2[i - 1];
		}
	}


	/**
	 * Method declaration
	 *
	 *
	 * @param i
	 * @param k
	 *
	 * @return
	 *
	 * @see
	 */
	static int BitReverse(int i, int k) {
		int rev = 0;

		for (int j = 0; j < k; ++j) {
			if ((i & pow2[k - j - 1]) > 0) {
				rev |= pow2[j];
			}
		}

		return rev;
	}


	/**
	 * Method declaration
	 *
	 *
	 * @param d
	 *
	 * @see
	 */
	static void bitDump(int d) {
		System.out.println(d + ":");
		for (int i = 31; i >= 0; i--) {
			if ((d & pow2[i]) != 0) {
				System.out.print("1");
			} else {
				System.out.print("0");
			}
		}

		System.out.println();
	}


	/**
	 * Dummy Constructor
	 */
	public FFT() {}


	/**
	 * Do Forward FFT calculation
	 */
	static public ComplexImage forward(Image im) {
		return doFFT(im, true);
	}


	/**
	 * Do Inverse FFT calculation
	 */
	static public ComplexImage inverse(Image im) {
		return doFFT(im, false);
	}


	/**
	 * Do FFT calculation @param im image to do FFT on.
	 * @param forward true=forward false=inverse
	 */
	static public ComplexImage doFFT(Image im, boolean forward) {

		int				x = 0;
		int				lgn = 0;
		int				n = 0;
		int				rev = 0;
		int				size = im.X() * im.Y();
		float[][] a = null;

		while (x <= MAX_POW) {

			// find the next-highest power of two after the total size of im
			if (pow2[x] >= size) {
				lgn = x;
				n = pow2[lgn];

				a = new float[2][n];

				// load the values from im into the new ComplexImage,
				// in reverse bit order, padding with 0's.
				for (int y = 0; y < n; y++) {
					rev = BitReverse(y, lgn);
					try {
						if (forward == true) {
							if (im instanceof GrayImage) {
								GrayImage ub = (GrayImage) im;

								if (y < (ub.X() * ub.Y())) {
									a[0][rev] = (float) ub.get(y % ub.X(), y / ub.X());
									a[1][rev] = 0f;
								} else {
									a[0][rev] = 0f;
								}
								a[1][rev] = 0f;
							} else if (im instanceof ComplexImage) {
								ComplexImage	ci = (ComplexImage) im;

								if (y < (ci.X() * ci.Y())) {
									a[0][rev] = (float) ci.getReal(y % im.X(), y / im.X());
									a[1][rev] = (float) ci.getImag(y % im.X(), y / im.X());
								} else {
									a[0][rev] = 0f;
									a[1][rev] = 0f;
								}
							} else if (im instanceof RealGrayImage) {
								RealGrayImage ub = (RealGrayImage) im;

								if (y < (ub.X() * ub.Y())) {
									a[0][rev] = ub.get(y % ub.X(), y / ub.X());
									a[1][rev] = 0f;
								} else {
									a[0][rev] = 0f;
									a[1][rev] = 0f;
								}
							} else {

								// unsupported type - throw exception?
							}
						} else {
							if (im instanceof GrayImage) {
								GrayImage ub = (GrayImage) im;

								if (y < (ub.X() * ub.Y())) {
									a[0][y] = (float) ub.get(y % ub.X(), y / ub.X()) / (float) n;
								} else {
									a[0][rev] = (float) 0f;
								}
								a[1][rev] = (float) 0f;
							} else if (im instanceof ComplexImage) {
								ComplexImage	ci = (ComplexImage) im;

								if (y < (ci.X() * ci.Y())) {
									a[0][rev] = (float) ci.getReal(y % ci.X(), y / ci.X()) / (float) n;
									a[1][rev] = (float) ci.getImag(y % ci.X(), y / ci.X()) / (float) n;
								} else {
									a[0][rev] = (float) 0f;
									a[1][rev] = (float) 0f;
								}
							} else if (im instanceof RealGrayImage) {
								RealGrayImage ub = (RealGrayImage) im;

								if (y < (ub.X() * ub.Y())) {
									a[0][y] = ub.get(y % ub.X(), y / ub.X()) / (float) n;
								} else {
									a[0][rev] = (float) 0f;
								}
								a[1][rev] = (float) 0f;
							} else {

								// unsupported type - throw exception?
							}
						}
					}		// try
					 catch (ArrayIndexOutOfBoundsException e) {
						a[0][rev] = (float) 0f;
						a[1][rev] = (float) 0f;
					} catch (Exception e) {
						System.out.println("FFT Error:  " + e.toString());
					}
				}			// for
				 break;
			}				// if
			 x++;
		}

		// float[] t = new float[2]; // dieses Array wird durch 2 Variablen substituiert [JM]
		// float[] u = new float[2]; // dieses Array wird durch 2 Variablen substituiert [JM]
		float			t0, t1, u0, u1; // neue Variablen [JM]
		float			temp;
		float			t2;
		int			    m = 0;
		int			    mdiv2 = 0;

		// float[] omega_m = null; // dieses Array wird durch 2 Variablen substituiert [JM]
		float			omega_m0 = 0f, omega_m1 = 0f; // neue Variablen [JM]
		float			omega0, omega1; // neue Hilfsvariablen [JM]
		float[][] b = new float[2][n]; // Initialisierung von b ausserhalb der Hauptschleife [JM]
		int				i, j, k, s; // Initialisierung ausserhalb der Hauptschleife [JM]
		int				w, h, kplusmdiv2; // neue Hilfsvariablen [JM]

		w = im.X(); // neue Hilfsvariable: Bildbreite [JM]
		h = im.Y(); // neue Hilfsvariable: Bildhoehe [JM]
		double	minustwotimesPIdivm; // neue Hilfsvariable: spaetere Zuweisung von (-2.0) * Math.PI / m [JM]
		double	minustwotimesPI = (-2.0) * Math.PI; // neue Hilfsvariable [JM]
		double	twotimesPIdivm; // neue Hilfsvariable: spaetere Zuweisung von (2.0) * Math.PI / m [JM]
		double	twotimesPI = (2.0) * Math.PI; // neue Hilfsvariable [JM]

		// Hilfsvariable zur Messung der Laufzeit [JM]
        // Wenn zu Debug-Zwecken benoetigt, bitte wieder einkommentieren
        // long    starttime = System.currentTimeMillis();

		// do the FFT
		for (s = 1; s <= lgn; s++) {
			m = pow2[s];

			// mdiv2 = m / 2; // auskommentierter Originalcode [JM]
			mdiv2 = m >> 1;	// neue Substitution [JM]

			// omega_m = new float[2]; // auskommentierter Originalcode [JM]

			if (forward == true) {
                // omega_m[0] = (float)(Math.cos((-2.0 * Math.PI) / (float) m)); // auskommentierter Originalcode [JM]
		        // omega_m[1] = (float)(Math.sin((-2.0 * Math.PI) / (float) m)); // auskommentierter Originalcode [JM]
                minustwotimesPIdivm = minustwotimesPI / (double) m;  // neu [JM]
				omega_m0 = (float) (Math.cos(minustwotimesPIdivm));  // neu [JM]
				omega_m1 = (float) (Math.sin(minustwotimesPIdivm));  // neu [JM]
			} else {
                // omega_m[0] = (float)(Math.cos(2.0 * Math.PI / (double) m)); // auskommentierter Originalcode [JM]
	            // omega_m[1] = (float)(Math.sin(2.0 * Math.PI / (double) m)); // auskommentierter Originalcode [JM]
                twotimesPIdivm = twotimesPI / (double) m; // neu [JM]
				omega_m0 = (float) (Math.cos(twotimesPIdivm)); // neu [JM]
				omega_m1 = (float) (Math.sin(twotimesPIdivm)); // neu [JM]
			}

			// float[] omega = {(float)1f, (float)0f};  // auskommentierter Originalcode [JM]
			omega0 = 1f; // neu [JM]
			omega1 = 0f; // neu [JM]

			// for (int j = 0; j < mdiv2; ++j) { // auskommentierter Originalcode [JM]
	        for (j = 0; j < mdiv2; ++j) { // neu [JM]
				for (k = j; k < n; k = k + m) {

	                kplusmdiv2 = k + mdiv2; // neu [JM]
					//t[0] = omega[0]; t[1] = omega[1]; // auskommentierter Originalcode [JM]
                    //u[0] = a[0][k]; u[1] = a[1][k]; // auskommentierter Originalcode [JM]
                    t0 = omega0; // neu [JM]
					t1 = omega1; // neu [JM]
					u0 = a[0][k]; // neu [JM]
					u1 = a[1][k]; // neu [JM]

					//temp = (t[0]*a[0][k+mdiv2]) - (t[1]*a[1][k+mdiv2]);  // auskommentierter Originalcode [JM]
	                temp = (t0 * a[0][kplusmdiv2]) - (t1 * a[1][kplusmdiv2]); // neu [JM]
					//t[1] = (t[0]*a[1][k+mdiv2]) + (t[1]*a[0][k+mdiv2]);  // auskommentierter Originalcode [JM]
	                t1 = (t0 * a[1][kplusmdiv2]) + (t1 * a[0][kplusmdiv2]); // neu [JM]
					//t[0] = temp;  // auskommentierter Originalcode [JM]
	                t0 = temp; // neu [JM]

					//a[k] = u + t; better way to do this?
	                //a[0][k] = u[0]+t[0]; a[1][k] = u[1]+t[1];   // auskommentierter Originalcode [JM]
	                a[0][k] = u0 + t0; // neu [JM]
					a[1][k] = u1 + t1; // neu [JM]
					//a[0][k+mdiv2] = u[0]-t[0]; a[1][k+mdiv2] = u[1]-t[1];   // auskommentierter Originalcode [JM]
                    a[0][kplusmdiv2] = u0 - t0; // neu [JM]
					a[1][kplusmdiv2] = u1 - t1; // neu [JM]
				}


				//t2 = (omega[0]*omega_m[0]) - (omega[1]*omega_m[1]); // auskommentierter Originalcode [JM]
	            t2 = (omega0 * omega_m0) - (omega1 * omega_m1); // neu [JM]
				//omega[1] = (omega[0]*omega_m[1]) + (omega[1]*omega_m[0]); // auskommentierter Originalcode [JM]
	            omega1 = (omega0 * omega_m1) + (omega1 * omega_m0); // neu [JM]
				//omega[0] = t2; // auskommentierter Originalcode [JM]
                omega0 = t2; // neu [JM]
			}
		}

		// ComplexImage ar = new ComplexImage(im.X(),im.Y());  // auskommentierter Originalcode [JM]
		ComplexImage ar = new ComplexImage(w, h);  // neu [JM]

		//for (int i=0; i<n; i++) {  // auskommentierter Originalcode [JM]
        for (i = 0; i < n; i++) {  // neu [JM]
			// ar.set(i % w, i / w, (float)a[0][i], (float)a[1][i]);  // auskommentierter Originalcode [JM]
			ar.set(i % w, i / w, a[0][i], a[1][i]);  // neu [JM]
		}

        // Hilfsvariable zur Messung der Laufzeit [JM]
        // Wenn zu Debug-Zwecken benoetigt, bitte wieder einkommentieren
        // long endtime = System.currentTimeMillis();

        // Ausgabe der Laufzeitmessung im ImageJ-Hauptfenster [JM]
        // Wenn zu Debug-Zwecken benoetigt, bitte wieder einkommentieren
        // IJ.write("Dauer der FFT: " + (endtime - starttime) + "ms");
		return ar;
	}
}



/*--- formatting done in "My Own Convention" style on 08-08-2000 ---*/

