content: package edu.stanford.rsl.conrad.utils;
  
  import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;
  import edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D;
  import edu.stanford.rsl.conrad.data.numeric.Grid2D;
  import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators;
  import edu.stanford.rsl.conrad.filtering.rampfilters.RampFilter;
  
  
  /**
   * This class is a wrapper for the FFT as implemented in JTransforms.
   * 
   * 
   * @author Andreas Maier
   * 
   */
  
  public abstract class FFTUtil {
  
  	public static void init1DFFT(int width){
  		DoubleFFT_1D fft = new DoubleFFT_1D(FFTUtil.getNextPowerOfTwo(width));
  		double [] test = new double [FFTUtil.getNextPowerOfTwo(width) * 2];
  		fft.complexForward(test);
  		fft.complexInverse(test, true);
  	}
  	
  	public static double []  fft(double []  array){
  		DoubleFFT_1D fft = new DoubleFFT_1D(FFTUtil.getNextPowerOfTwo(array.length));
  		double [] test = new double [FFTUtil.getNextPowerOfTwo(array.length) * 2];
  		for (int i = 0; i < array.length; i++){
  			test[i*2] = array[i];
  		}
  		fft.complexForward(test);
  		return test;
  	}
  	
  	public static double []  fft(double []  array, int padding){
  		DoubleFFT_1D fft = new DoubleFFT_1D(FFTUtil.getNextPowerOfTwo(array.length + padding));
  		double [] test = new double [FFTUtil.getNextPowerOfTwo(array.length+padding) * 2];
  		for (int i = 0; i < array.length; i++){
  			test[i*2] = array[i];
  		}
  		fft.complexForward(test);
  		return test;
  	}
  
  	/**
  	 * Computes the complex signum for the given array. 
  	 * @param array the array
  	 */
  	public static double [] complexSignum (double [] array){
  		double [] revan = new double [array.length];
  		for (int i=0;i<array.length/2;i++){
  			double abs = abs(i, array);
  			if (abs != 0) {
  				revan[i*2] = array[i*2] / abs;
  				revan[(i*2)+1] = array[(i*2)+1] / abs;
  			} else {
  				revan[i*2] = 0;
  				revan[(i*2)+1] = 0;
  			}
  		}
  		return revan;
  	}
  
  	/**
  	 * Computes the signum for the given complex array. 
  	 * @param array the array
  	 */
  	public static double [] cSignum (double [] array){
  		double [] revan = new double [array.length];
  		for (int i=0;i<array.length/2;i++){
  			double abs = abs(i, array);
  			if (abs != 0) {
  				double val = 1;
  				if (array[i*2] < 0 ||(array[i*2]==0 && array[(i*2)+1]<0)) val =-1;
  				revan[i*2] = val;
  				revan[(i*2)+1] = 0;
  			} else {
  				revan[i*2] = 0;
  				revan[(i*2)+1] = 0;
  			}
  		}
  		return revan;
  	}
  
  
  	/**
  	 * Multiplies the array with the complex number defined as two double values;
  	 * @param array the array
  	 * @param realOne the real part of the complex number
  	 * @param imagOne the imaginary part of the complex number
  	 */
  	public static void timesComplexNumber(double [] array, double realOne, double imagOne){
  		for (int i=0;i<array.length/2;i++){
  			double [] newValue = multiplyComplex(array[i*2], array[(i*2)+1], realOne, imagOne);
  			array[i*2] = newValue[0];
  			array[(i*2)+1] = newValue[1];
  		}
  	}
  
  	public static double[] complexOneOverPiT(int width){
  		double [] revan = new double [width];
  		for(int i=0; i< width/2; i++){
  			revan[i*2] = 1/((i-(width/4))*Math.PI);
  			if ((i-(width/4))==0){
  				revan[i*2] =10000;
  			}
  		}
  		return revan;
  	}
  
  	public static double[] hilbertKernel(int width){
  		double [] revan = new double [width];
  		for(int i=0; i< width/2; i++){
  			revan[i*2] = Math.signum(i-(width/4));
  		}
  		return revan;
  	}
  
  	/**
  	 * Computes the 1-D Hilbert Transform of the given Array. Uses the Fast Fourier transform to compute the transform
  	 * @param array the array
  	 * @param nTimes the size of the zero padding:<pre>
  	 *            n = 1 zero pad to next power of 2
  	 *            n = 2 zero pad to 2 times the next power of 2</pre>
  	 * @return the Hilbert transformed array
  	 */
  	public static double [] hilbertTransform(double [] array, int nTimes){
  		int newDimension = FFTUtil.getNextPowerOfTwo(array.length);
  		// Write the data to JTransforms complex format
  		double [] fftData = new double[newDimension*2*nTimes];
  		for (int i = 0 ; i < array.length; i++){
  			fftData[i*2] = array[i];
  		}
  		// FFT
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension*nTimes);
  		fft.complexForward(fftData);
  		// compute Hilbert Transform
  		for (int i=0;i<fftData.length/2;i++){
  			double sign = Math.signum((fftData.length/4)-i);
  			if (sign == 0) {
  			} else {
  				if (i==0) sign =0;
  				if(i == (fftData.length/2) -1) sign =0;
  				double real = sign * -1.0 * fftData[(2*i)+1]; 
  				double imag = sign * 1.0 * fftData[(2*i)];
  				fftData[2*i] = real;
  				fftData[(2*i)+1] = imag;
  			}
  		}
  		// iFFT
  		fft.complexInverse(fftData, true);
  		// rewrite as real array.
  		double [] revan = new double[array.length];
  		for (int i = 0 ; i < array.length; i++){
  			revan[i] = 
  				//abs( i, fftData);
  				fftData[2*i];
  		}
  		return revan;
  	}
  
  	public static double [] discreteHilbertTransform(double [] array, int nTimes){
  		double [] revan = new double [array.length];
  		for (int i =0; i<array.length; i++){
  			double sum = 0;
  			if(i % 2 == 0){ // i even
  				for (int j = 0; j < array.length; j++){
  					if (j % 2 == 1){ // j odd
  						sum += array[j] / (i - j);
  					}
  				}
  				sum *= 2.0 / Math.PI;
  			} else { // i odd
  				for (int j = 0; j < array.length; j++){
  					if (j % 2 == 0){ // j even
  						sum += array[j] / (i - j);
  					}
  				}
  				sum *= 2.0 / Math.PI;
  			}
  			revan[i] = sum;
  		}
  		return revan;
  	}
  	/**
  	 * Exact finite discrete Hilbert transformation after Kak SC. Hilbert Transformation for discrete data. Int J Electronics 34(2):177-83. 1973. Eq. 14
  	 * @param array the Array
  	 * @param nTimes Approximation for infinity
  	 * @return the array with the Hilbert transform
  	 */
  	public static double [] exactFiniteHilbertTransform(double [] array, int nTimes){
  		double [] revan = new double [array.length];
  		for (int n =0; n<array.length; n++){
  			double sum1 = 0;
  			for (int k = 0; k< array.length; k++){
  				if (n-k !=0){
  					double cos = Math.cos(Math.PI*(n-k));
  					double sum2 = (1.0 - cos) / (Math.PI*(n-k));
  					for (int p = 1; p < nTimes; p++){
  						sum2 += (2.0 / Math.PI) *
  						(((n-k)*(1.0-(Math.pow(-1.0, p) * cos))) / 
  								(Math.pow(n-k,2) - (Math.pow(p,2)*Math.pow(array.length, 2)))
  						);
  					}
  					sum1 += array[k] * sum2;
  					//System.out.println(sum1);
  				}
  			}
  			revan[n] = sum1;
  		}
  		return revan;
  	}
  
  
  	/**
  	 * Estimates an inverse of a Ramp Filter. Works in place. The array in the argument list will be changed!
  	 * @param filter the filter
  	 * @param max the maximal scaling (1.0 is a good number)
  	 * @return the inverted RampFilter
  	 */
  	public static double [] invertRampFilter(double [] filter, double max){
  		double maximum = -Double.MAX_VALUE;
  		for (int i=0; i< filter.length/2; i++){
  			// Compute Inverse
  			double value = 1.0 / abs(i, filter);
  			// set to real part
  			filter[2*i] = value;
  			// set complex part to 0
  			filter[(2*i) + 1] = 0;
  			// update maximum if needed
  			if (value > maximum) maximum = value;		
  		}
  		// compute scaling factor
  		double scale = max / maximum;
  		// scale the inverted filter
  		boolean scaling = true;
  		if (scaling) {
  			for (int i=0; i< filter.length/2; i++){
  				if (Double.isInfinite(filter[2*i])) {
  					filter [2*i] = 0;
  				} else {
  					filter [2*i] *= scale;
  				}
  			}
  		}
  		return filter;
  	}
  
  	/**
  	 * Removes a ramp filter from the ImageProcessor. Note that frequencies which were multiplied with 0 cannot
  	 * be restored. This has to be estimated differently. (If you apply the same filter later on again, it should
  	 * not matter anyway, as the same frequencies will be multiplied with 0 again.)
  	 * 
  	 * @param imp the filtered ImageProcessor
  	 * @param ramp the ramp
  	 * @param max the maximum for the inverted ramp filter. (try 1.0)
  	 * @return the unfiltered ImageProcessor
  	 */
  	public static Grid2D removeRampFilter(Grid2D imp, RampFilter ramp, double max){
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = imp.getWidth();
  		originalHeight = imp.getHeight();
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		// Write the data to JTransforms format
  		double [] fftData = new double[newDimension*2];
  		double [] filter = ramp.getRampFilter1D(newDimension);
  		invertRampFilter(filter, max);
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		Grid2D revan = new Grid2D(originalWidth, originalHeight);
  		revan.setOrigin(imp.getOrigin());
  		revan.setSpacing(imp.getSpacing());
  		for (int j = 0;  j < imp.getHeight(); j++){
  			for (int i = 0 ; i < imp.getWidth(); i++){
  				fftData[i*2] = imp.getPixelValue(i, j);
  			}
  			// Perform forward transform
  			fft.complexForward(fftData);
  			// Filter
  			fftData = FFTUtil.multiplyAbsolute(fftData, filter);
  			// Perform backward transform
  			fft.complexInverse(fftData, true);
  			for (int i = 0 ; i< revan.getWidth(); i++){
  				double value = fftData[i*2];
  				revan.putPixelValue(i, j, value);
  			}
  		}
  		boolean adjust = false;
  		if (adjust){
  			double [] scale = new double [revan.getWidth()];
  			double [] stats = new double [2];
  			stats[0] = NumericPointwiseOperators.mean(revan);
  			stats[1] = NumericPointwiseOperators.stddev(revan, stats[0]);
  			System.out.println("Mean :"  + stats[0]);
  			double avgBorder = 0;
  			for (int i = 0; i < revan.getHeight(); i++){
  				avgBorder += (revan.getPixelValue(0, i) + revan.getPixelValue(revan.getWidth()-1, i)) / (2* revan.getHeight());
  			}
  			double factor = ((stats[0] - Math.abs(stats[0] - avgBorder)) / avgBorder) - 1;
  			System.out.println("Factor :"  + factor);
  			for (int i = 0; i < scale.length; i++){
  				scale[i] =  (1 + (factor * Math.cos(i * 2 * Math.PI / revan.getWidth())));
  			}
  			for (int j = 0;  j < imp.getHeight(); j++){
  				for (int i = 0 ; i< revan.getWidth(); i++){
  					double value = revan.getPixelValue(i, j) * scale[i];
  					revan.putPixelValue(i, j, value);
  				}
  			}
  		}
  		return revan;
  	}	
  
  
  	/**
  	 * low pass filters the given array. All indices with a distance  greater or equal to cutOffIndex from the highest frequency are set to 0;
  	 * @param array the real double array
  	 * @param cutOffIndex the index to start cutting
  	 * @return the low pass filtered array.
  	 */
  	public static double [] lowPassFilterRealDoubleArray (double [] array, int cutOffIndex){
  		int newDimension = FFTUtil.getNextPowerOfTwo(array.length);
  		// Write the data to JTransforms complex format
  		double [] fftData = new double[newDimension*2];
  		for (int i = 0 ; i < array.length; i++){
  			fftData[i*2] = array[i];
  		}
  		// FFT
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		fft.complexForward(fftData);
  		// low pass data
  		for (int i = 0 ; i < cutOffIndex; i++){
  			int center = newDimension;
  			fftData[center + (i*2)] = 0;
  			fftData[center + (i*2) + 1] = 0;
  			fftData[center - (i*2)] = 0;
  			fftData[center - (i*2) + 1] = 0;
  		}
  		// iFFT
  		fft.complexInverse(fftData, true);
  		// rewrite as real array.
  		double [] revan = new double[newDimension];
  		for (int i = 0 ; i < array.length; i++){
  			revan[i] = fftData[2*i];
  		}
  		return revan;
  	}
  	
  
  
  	/**
  	 * Applies a ramp filter to the an ImageProcessor
  	 * @param imp the ImageProcessor to be filtered
  	 * @param ramp the ramp
  	 * @return the filtered ImageProcessor
  	 */
  	public static Grid2D applyRampFilter_ECC(Grid2D imp, RampFilter ramp){
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		boolean error = false;
  		originalWidth = imp.getWidth();
  		originalHeight = imp.getHeight();
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		// Write the data to JTransforms format
  		double [] fftData = new double[newDimension*2];
  		double [] filter = ramp.getRampFilter1D(newDimension);
  		if (DoubleArrayUtil.isNaN(filter)){
  			System.out.println("NaN found in Filter!");
  		}
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		Grid2D revan = new Grid2D(originalWidth, originalHeight);
  		revan.setOrigin(imp.getOrigin());
  		revan.setSpacing(imp.getSpacing());
  		for (int j = 0;  j < imp.getHeight(); j++){
  			for (int i = 0 ; i < newDimension *2; i++){
  				fftData[i] = 0;
  			}
  			if (DoubleArrayUtil.isNaN(fftData)){
  				System.out.println("Never happens!");
  			}
  			for (int i = 0 ; i < imp.getWidth(); i++){
  				fftData[i*2] = imp.getPixelValue(i, j);
  			}
  			if (DoubleArrayUtil.isNaN(fftData)){
  				System.out.println("NaN found in input Data! line " + j);
  				for (int i = 0 ; i < imp.getWidth(); i++){
  					if (Double.isNaN(fftData[i*2])){
  						System.out.println("index " + i);
  					}
  				}
  				error = true;
  			}
  			// Perform forward transform
  			fft.complexForward(fftData);
  			if (DoubleArrayUtil.isNaN(fftData)){
  				System.out.println("NaN found after FFT!");
  			}
  			// Filter
  			fftData = FFTUtil.multiplyAbsolute(fftData, filter);
  			if (DoubleArrayUtil.isNaN(fftData)){
  				System.out.println("NaN found after filter application!");
  			}
  			// Perform backward transform
  			fft.complexInverse(fftData, true);
  			if (DoubleArrayUtil.isNaN(fftData)){
  				System.out.println("NaN found after iFFT!");
  			}
  			for (int i = 0 ; i< revan.getWidth(); i++){
  				double value = fftData[i*2];
  				revan.putPixelValue(i, j, value);
  			}
  		}
  		if (error) {
  			//ImagePlus gi = ImageUtil.wrapGrid3D(imp, "");
  			//gi.setTitle("Errors");
  			//gi.show();
  		}
  		return revan;
  	}
  
  	/**
  	 * Applies a ramp filter to the an ImageProcessor
  	 * @param imp the ImageProcessor to be filtered
  	 * @param ramp the ramp
  	 * @return the filtered ImageProcessor
  	 */
  	public static Grid2D applyRampFilter(Grid2D imp, RampFilter ramp){
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = imp.getWidth();
  		originalHeight = imp.getHeight();
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN); // Satisfy Nyquist?
  		// Write the data to JTransforms format
  		double [] fftData = new double[newDimension*2];	// Zero padding?
  		double [] filter = ramp.getRampFilter1D(newDimension);		
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		Grid2D revan = new Grid2D(originalWidth, originalHeight);
  		revan.setOrigin(imp.getOrigin());
  		revan.setSpacing(imp.getSpacing());
  		for (int j = 0;  j < imp.getHeight(); j++){
  			for (int i = 0 ; i < newDimension *2; i++){
  				fftData[i] = 0;
  			}
  			for (int i = 0 ; i < imp.getWidth(); i++){
  				fftData[i*2] = imp.getPixelValue(i, j);
  			}
  			// Perform forward transform
  			fft.complexForward(fftData);
  			// Filter
  			fftData = FFTUtil.multiplyAbsolute(fftData, filter);
  			// Perform backward transform
  			fft.complexInverse(fftData, true);
  			for (int i = 0 ; i< revan.getWidth(); i++){
  				double value = fftData[i*2];
  				revan.putPixelValue(i, j, value);
  			}
  		}
  		return revan;
  	}
  
  	/**
  	 * Applies a ramp filter to the a detector row
  	 * @param detectorRow the row
  	 * @param ramp the ramp
  	 * @return the filtered row
  	 */
  	public static double [] applyRampFilter(double [] detectorRow, RampFilter ramp){
  		int originalWidth;
  		int maxN;
  		originalWidth = detectorRow.length;
  		maxN = originalWidth;
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		
  		// Write the data to JTransforms format
  		double [] fftData = new double[newDimension*2];
  		double [] filter = ramp.getRampFilter1D(newDimension);		
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		double [] revan = new double[originalWidth];
  		for (int i = 0 ; i < newDimension *2; i++){
  			fftData[i] = 0;
  		}
  		for (int i = 0 ; i < detectorRow.length; i++){
  			fftData[i*2] = detectorRow[i];
  		}
  		// Perform forward transform
  		fft.complexForward(fftData);
  		// Filter
  		fftData = FFTUtil.multiplyAbsolute(fftData, filter);
  		// Perform backward transform
  		fft.complexInverse(fftData, true);
  		for (int i = 0 ; i< revan.length; i++){
  			double value = fftData[i*2];
  			revan[i] = value;
  		}
  		return revan;
  	}
  
  	/**
  	 * Applies a 2D filter to the an ImageProcessor
  	 * @param imp the ImageProcessor to be filtered
  	 * @param filter the filter
  	 * @return the filtered ImageProcessor
  	 */
  	public static Grid2D apply2DFilter(Grid2D imp, Grid2D filter){
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = imp.getWidth();
  		originalHeight = imp.getHeight();
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		// Write the data to JTransforms format
  		double [][] fftData = new double[newDimension][newDimension*2];
  		DoubleFFT_2D fft = new DoubleFFT_2D(newDimension, newDimension);
  		Grid2D revan = new Grid2D(originalWidth, originalHeight);
  		revan.setOrigin(imp.getOrigin());
  		revan.setSpacing(imp.getSpacing());
  		for (int j = 0;  j < imp.getHeight(); j++){
  			for (int i = 0 ; i < imp.getWidth(); i++){
  				fftData[i][j*2] = imp.getPixelValue(i, j);
  			}
  		}
  		// Perform forward transform
  		fft.complexForward(fftData);
  		// Filter
  		for (int j = 0;  j < filter.getHeight(); j++){
  			for (int i = 0 ; i < filter.getWidth(); i++){
  				fftData[i][j*2] *= filter.getPixelValue(i, j);
  				fftData[i][(j*2)+1] *= filter.getPixelValue(i, j);
  			}
  		}
  		// Perform backward transform
  		fft.complexInverse(fftData, true);
  		for (int j = 0;  j < imp.getHeight(); j++){
  			for (int i = 0 ; i< revan.getWidth(); i++){
  				double value = fftData[i][j*2];
  				revan.putPixelValue(i, j, value);
  			}
  		}
  		return revan;
  	}
  
  
  	/**
  	 * Divides two complex values
  	 * @param realOne real part one
  	 * @param imagOne imaginary part one
  	 * @param realTwo real part two
  	 * @param imagTwo imaginary part two
  	 * @return an array of two values: first entry is real, second imaginary
  	 */
  	public static double [] divideComplex(double realOne, double imagOne, double realTwo, double imagTwo){
  		double [] revan = new double [2];
  		double denominator = Math.pow(realTwo, 2) + Math.pow(imagTwo, 2);
  		revan [0] = ((realOne * realTwo) + (imagOne * imagTwo)) / denominator;
  		revan [1] = ((imagOne * realTwo) - (realOne * imagTwo)) / denominator;
  		return revan;
  	}
  
  	/**
  	 * Divides two arrays of complex numbers. The absolute value of the second array is used to divide the complex first array.
  	 * @param one the first array
  	 * @param two the second array
  	 * @return the array of null if the lengths of the arrays don't match
  	 */
  	public static double [] divideAbsolute(double [] one, double [] two){
  		double [] revan = null;
  		if (one.length==two.length) {
  			revan = new double [one.length];
  			for (int i = 0; i < one.length/2; i++){
  				revan[2*i] = (one[2*i] / abs(i, two));
  				revan[(2*i)+1] = (one[(2*i)+1] / abs(i, two));
  			}
  		}
  		return revan;
  	}
  
  	/**
  	 * Multiplies two complex values
  	 * @param realOne real part one
  	 * @param imagOne imaginary part one
  	 * @param realTwo real part two
  	 * @param imagTwo imaginary part two
  	 * @return an array of two values: first entry is real, second imaginary
  	 */
  	public static double [] multiplyComplex(double realOne, double imagOne, double realTwo, double imagTwo){
  		double [] revan = new double [2];
  		revan [0] = (realOne * realTwo) - (imagOne * imagTwo);
  		revan [1] = (imagOne * realTwo) + (realOne * imagTwo);
  		return revan;
  	}
  
  	/**
  	 * Multiplies two arrays of complex numbers pairwise
  	 * @param one the first array
  	 * @param two the second array
  	 * @return the array of null if the lengths of the arrays don't match
  	 */
  	public static double [] multiplyComplex(double [] one, double [] two){
  		double [] revan = null;
  		if (one.length==two.length) {
  			revan = new double [one.length];
  			for (int i = 0; i < one.length/2; i++){
  				revan[2*i] = (one[2*i] * two[2*i]) - (one[(2*i)+1] * two[(2*i)+1]);
  				revan[(2*i)+1] = (one[(2*i)+1] * two[2*i]) + (one[2*i] * two[(2*i)+1]);
  			}
  		}
  		return revan;
  	}
  
  	/**
  	 * Multiplies two arrays of complex numbers. The absolute value of the second array is multiplied to the complex first array.
  	 * @param one the first array
  	 * @param two the second array
  	 * @return the array of null if the lengths of the arrays don't match
  	 */
  	public static double [] multiplyAbsolute(double [] one, double [] two){
  		double [] revan = null;
  		if (one.length==two.length) {
  			revan = new double [one.length];
  			for (int i = 0; i < one.length/2; i++){
  				revan[2*i] = (one[2*i] * abs(i, two));
  				revan[(2*i)+1] = (one[(2*i)+1] * abs(i, two));
  			}
  		}
  		return revan;
  	}
  	/**
  	 * shift zero frequency to center, or vice verse, 1D.
  	 * @param data the double data to be shifted
  	 * @param bComplex true: complex; false: real
  	 * @param bSign true: fftshift; false: ifftshift
  	 * @return the fft shifted array
  	 */
  
  	public static double []  fftshift(double [] data, boolean bComplex, boolean bSign)
  	{
  		double [] revan = new double [data.length];
  
  		int step = 1;
  		if (bComplex) step = 2;
  		int len = data.length/step;
  		int p = 0;
  		if(bSign) 
  			p = (int) Math.ceil(len/2.0);
  		else
  			p = (int) Math.floor(len/2.0);
  
  		int i=0;
  		if (step==1){
  			for (i=p;i<len;i++){
  				revan[i-p] = data[i];
  			}
  			for (i=0;i<p;i++){
  				revan[i+len-p] = data[i];
  			}
  		}
  		else{
  			for (i=p;i<len;i++){
  				revan[(i-p)*2] = data[i*2];
  				revan[(i-p)*2+1] = data[i*2+1];
  			}
  			for (i=0;i<p;i++){
  				revan[(i+len-p)*2] = data[i*2];
  				revan[(i+len-p)*2+1] = data[i*2+1];
  			}			
  		}			
  		return revan;
  	}
  	
  	/**
  	 * shift zero frequency to center, or vice verse, 1D.
  	 * @param data the float data to be shifted
  	 * @param bComplex true: complex; false: real
  	 * @param bSign true: fftshift; false: ifftshift
  	 * @return the fft shifted array
  	 */
  
  	public static float []  fftshift(float [] data, boolean bComplex, boolean bSign)
  	{
  		float [] revan = new float [data.length];
  
  		int step = 1;
  		if (bComplex) step = 2;
  		int len = data.length/step;
  		int p = 0;
  		if(bSign) 
  			p = (int) Math.ceil(len/2.0);
  		else
  			p = (int) Math.floor(len/2.0);
  
  		int i=0;
  		if (step==1){
  			for (i=p;i<len;i++){
  				revan[i-p] = data[i];
  			}
  			for (i=0;i<p;i++){
  				revan[i+len-p] = data[i];
  			}
  		}
  		else{
  			for (i=p;i<len;i++){
  				revan[(i-p)*2] = data[i*2];
  				revan[(i-p)*2+1] = data[i*2+1];
  			}
  			for (i=0;i<p;i++){
  				revan[(i+len-p)*2] = data[i*2];
  				revan[(i+len-p)*2+1] = data[i*2+1];
  			}			
  		}			
  		return revan;
  	}
  	
  	/**
  	 * shift zero frequency to center, or vice verse, 2D.
  	 * @param data the double data to be shifted
  	 * @param bComplex true: complex; false: real
  	 * @param bSign true: fftshift; false: ifftshift
  	 * @return the fft shifted array
  	 */
  
  	public static double[][]  fftshift(double [][] data, boolean bComplex, boolean bSign)
  	{
  		int step = 1;
  		if (bComplex) step = 2;
  		int height = data.length;
  	
  		int width = data[0].length/step;
  		
  		double [][] revan = new double [data.length][data[0].length];
  
  		int pH = 0;
  		int pW = 0;
  		if(bSign) {
  			pH = (int) Math.ceil(height/2.0);
  			pW = (int) Math.ceil(width/2.0);
  		}
  		else{
  			pH = (int) Math.floor(height/2.0);
  			pW = (int) Math.floor(width/2.0);
  		}
  		int i=0;
  		int j=0;
  		if (step==1){
  			for(j=pH;j<height;j++){
  				for (i=pW;i<width;i++){				
  					revan[j-pH][i-pW] = data[j][i];
  				}
  				for (i=0;i<pW;i++){
  					revan[j-pH][i+width-pW] = data[j][i];
  				}				
  			}
  			for(j=0;j<pH;j++){
  				for (i=pW;i<width;i++){				
  					revan[j+height-pH][i-pW] = data[j][i];
  				}
  				for (i=0;i<pW;i++){
  					revan[j+height-pH][i+width-pW] = data[j][i];
  				}				
  			}
  		}
  		else{
  			for(j=pH;j<height;j++){
  				for (i=pW;i<width;i++){				
  					revan[j-pH][(i-pW)*2] = data[j][i*2];
  					revan[j-pH][(i-pW)*2+1] = data[j][i*2+1];
  				}
  				for (i=0;i<pW;i++){
  					revan[j-pH][(i+width-pW)*2] = data[j][i*2];
  					revan[j-pH][(i+width-pW)*2+1] = data[j][i*2+1];
  				}				
  			}
  			for(j=0;j<pH;j++){
  				for (i=pW;i<width;i++){				
  					revan[j+height-pH][(i-pW)*2] = data[j][i*2];
  					revan[j+height-pH][(i-pW)*2+1] = data[j][i*2+1];					
  				}
  				for (i=0;i<pW;i++){
  					revan[j+height-pH][(i+width-pW)*2] = data[j][i*2];
  					revan[j+height-pH][(i+width-pW)*2+1] = data[j][i*2+1];
  				}				
  			}		
  		}			
  		return revan;
  	}
  	
  	/**
  	 * shift zero frequency to center, or vice verse, 2D.
  	 * @param data the double data to be shifted
  	 * @param bComplex true: complex; false: real
  	 * @param bSign true: fftshift; false: ifftshift
  	 * @return the fft shifted array
  	 */
  
  	public static float[][]  fftshift(float [][] data, boolean bComplex, boolean bSign)
  	{
  		int step = 1;
  		if (bComplex) step = 2;
  		int height = data.length;
  	
  		int width = data[0].length/step;
  		
  		float [][] revan = new float [data.length][data[0].length];
  
  		int pH = 0;
  		int pW = 0;
  		if(bSign) {
  			pH = (int) Math.ceil(height/2.0);
  			pW = (int) Math.ceil(width/2.0);
  		}
  		else{
  			pH = (int) Math.floor(height/2.0);
  			pW = (int) Math.floor(width/2.0);
  		}
  		int i=0;
  		int j=0;
  		if (step==1){
  			for(j=pH;j<height;j++){
  				for (i=pW;i<width;i++){				
  					revan[j-pH][i-pW] = data[j][i];
  				}
  				for (i=0;i<pW;i++){
  					revan[j-pH][i+width-pW] = data[j][i];
  				}				
  			}
  			for(j=0;j<pH;j++){
  				for (i=pW;i<width;i++){				
  					revan[j+height-pH][i-pW] = data[j][i];
  				}
  				for (i=0;i<pW;i++){
  					revan[j+height-pH][i+width-pW] = data[j][i];
  				}				
  			}
  		}
  		else{
  			for(j=pH;j<height;j++){
  				for (i=pW;i<width;i++){				
  					revan[j-pH][(i-pW)*2] = data[j][i*2];
  					revan[j-pH][(i-pW)*2+1] = data[j][i*2+1];
  				}
  				for (i=0;i<pW;i++){
  					revan[j-pH][(i+width-pW)*2] = data[j][i*2];
  					revan[j-pH][(i+width-pW)*2+1] = data[j][i*2+1];
  				}				
  			}
  			for(j=0;j<pH;j++){
  				for (i=pW;i<width;i++){				
  					revan[j+height-pH][(i-pW)*2] = data[j][i*2];
  					revan[j+height-pH][(i-pW)*2+1] = data[j][i*2+1];					
  				}
  				for (i=0;i<pW;i++){
  					revan[j+height-pH][(i+width-pW)*2] = data[j][i*2];
  					revan[j+height-pH][(i+width-pW)*2+1] = data[j][i*2+1];
  				}				
  			}		
  		}			
  		return revan;
  	}
  	/**
  	 * Get Lowpass and HighPass Images for "real" image using "real" filters
  	 * @param data image in real
  	 * @param LowPassFilter Low pass filter in real
  	 * @param HighPassFilter High pass filter in real
  	 * @param LowPassImage Low pass image in real
  	 * @param HighPassImage High Pass image in real
  	 */
  	public static void  GetLowandHighPassImage(float [][] data, 
  			float [][] LowPassFilter, 
  			float [][] HighPassFilter, 
  			float [][] LowPassImage, 
  			float [][] HighPassImage)
  	{
  		int nHeight = data.length;
  		int nWidth = data[0].length;
  		double [][] temp1 = new double[nHeight][nWidth*2]; // for high pass 
  		double [][] temp2 = new double[nHeight][nWidth*2]; // for low pass
  		
  		for (int i=0;i<nHeight;i++){
  			for(int j=0;j<nWidth;j++){
  				temp1[i][2*j] = data[i][j];
  			}
  		}
  		//FloatFFT_2D fft2 = new FloatFFT_2D(nHeight,nWidth);
  		DoubleFFT_2D fft2 = new DoubleFFT_2D(nHeight,nWidth);
  		fft2.complexForward(temp1);
  		
  		for (int i=0;i<nHeight;i++){
  			for(int j=0;j<nWidth;j++){
  				temp2[i][2*j] = temp1[i][2*j] * LowPassFilter[i][j];
  				temp1[i][2*j] = temp1[i][2*j] * HighPassFilter[i][j];
  				temp2[i][2*j+1] = temp1[i][2*j+1] * LowPassFilter[i][j];	
  				temp1[i][2*j+1] = temp1[i][2*j+1] * HighPassFilter[i][j];
  			}
  		}
  	
  		temp1 = fftshift(temp1,true,false);
  		fft2.complexInverse(temp1,true);
  		fft2.complexInverse(temp2,true);
  		for (int i=0;i<nHeight;i++){
  			for(int j=0;j<nWidth;j++){
  				LowPassImage[i][j] = (float)temp2[i][2*j];// Math.sqrt(temp2[i][2*j]*temp2[i][2*j]+temp2[i][2*j+1]*temp2[i][2*j+1]);
  				HighPassImage[i][j] = (float)temp1[i][2*j];// Math.sqrt(temp1[i][2*j]*temp1[i][2*j]+temp1[i][2*j+1]*temp1[i][2*j+1]);
  			}
  		}
  		
  	}
  
  	/**
  	 * Returns the 2D power spectrum of a given image processor.
  	 * 
  	 * @param imp the ImageProcessor
  	 * @return the power spectrum as ImageProcessor 
  	 */
  	public static Grid2D getPowerSpectrum(Grid2D imp){
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = imp.getWidth();
  		originalHeight = imp.getHeight();
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		// Write the data to JTransforms format
  		double [][] fftData = new double[newDimension][newDimension*2];
  		DoubleFFT_2D fft = new DoubleFFT_2D(newDimension, newDimension);
  		for (int i = 0 ; i< imp.getWidth(); i++){
  			for (int j = 0;  j< imp.getHeight(); j++){
  				fftData[i][j*2] = imp.getPixelValue(i, j);
  			}
  		}
  		// Perform forward transform
  		fft.complexForward(fftData);
  		// Filter
  		Grid2D revan = new Grid2D(newDimension, newDimension);
  		for (int i = 0 ; i< revan.getWidth(); i++){
  			for (int j = 0;  j< revan.getHeight(); j++){
  				double value = FFTUtil.abs(i, j, fftData);
  				revan.putPixelValue(i, j, value);
  			}
  		}
  		return revan;
  	}
  
  	/**
  	 * Returns true if the number is a power of two
  	 * 
  	 * @param value the input number.
  	 */
  	public static boolean isPowerOfTwo(int value){
  		return Integer.bitCount(value)==1;
  	}
  	
  	/**
  	 * Returns the NEXT power of 2 given a certain integer value 
  	 * For power of twos it also returns the next power of 2, e.g. for 512 -> 1024
  	 * 
  	 * Code was partially taken from ij.plugin.FFT.java::pad().
  	 * Thanks for the inspiration!
  	 * 
  	 * @param value the input number.
  	 * @return the next power of two.
  	 */
  	public static int getNextPowerOfTwo(int value){
  		if (isPowerOfTwo(value)) 
  			return value*2;
  		else
  		{
  			int i = 2;
  			while (i <= value) {
  				i *= 2;
  			}
  			return i*2;
  		}
  	}
  
  	/**
  	 * Estimates the applied filter given an input and an output image.
  	 * 
  	 * @param before the input image
  	 * @param after the output image
  	 * @param threshold the maximal value which may appear in the estimate. (to avoid outliers.)
  	 * @return the estimate of the applied filter
  	 * @throws Exception may occur.
  	 */
  	public static Grid2D estimateFilter2D(Grid2D before, Grid2D after, double threshold) throws Exception{
  		if ((before.getWidth() != after.getWidth()|| (after.getHeight() != before.getHeight()))){
  			throw new Exception ("Image dimensions do not match!");
  		}
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = before.getWidth();
  		originalHeight = before.getHeight();
  		// Padding required? Pad if necessary
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		Grid2D filterEstimate = new Grid2D(newDimension, newDimension);
  		double [][] fftBefore = new double [newDimension][newDimension*2];
  		double [][] fftAfter = new double [newDimension][newDimension*2];
  		DoubleFFT_2D fft = new DoubleFFT_2D(newDimension, newDimension);
  		double outlier = 0;
  		for (int j = 0;  j < before.getHeight(); j++){ // for all image rows
  			// Convert to JTransforms format
  			for (int i = 0 ; i < newDimension; i++){
  				fftBefore[i][j*2] = before.getPixelValue(i, j);
  				fftBefore[i][(j*2)+1] = 0;
  				fftAfter[i][j*2] = after.getPixelValue(i, j);
  				fftAfter[i][(j*2)+1] = 0;
  			}
  		}
  		// Perform the FFTs
  		fft.complexForward(fftBefore);
  		fft.complexForward(fftAfter);
  		// estimate the filter
  		for (int j = 0;  j < newDimension; j++){ // for all image rows
  			for (int i = 0 ; i < newDimension; i++){
  				double beforeValue = abs(i, j, fftBefore);
  				double add = 0;
  				if (beforeValue != 0) {
  					add = abs(0, divideComplex(fftAfter[i][2*j], fftAfter[i][(2*j)+1], fftBefore[i][2*j], fftBefore[i][(2*j)+1]));
  					if (add < threshold) {
  						filterEstimate.putPixelValue(i, j, add);
  					} else {
  						outlier++;
  					}
  				} else {
  					outlier++;
  				}
  			}
  		}
  		//System.out.println("Outliers per frame: " + outlier);
  		return filterEstimate;
  	}
  
  
  	/**
  	 * Estimates the applied ramp filter given an input and an output image.
  	 * 
  	 * @param after the output image
  	 * @param before the input image
  	 * @return the estimate of the applied filter
  	 * @throws Exception may occur.
  	 */
  	public static double [] estimateFilter(Grid2D after, Grid2D before, double threshold, boolean meanSquare) throws Exception{
  		if ((before.getWidth() != after.getWidth()|| (after.getHeight() != before.getHeight()))){
  			throw new Exception ("Image dimensions do not match!");
  		}
  		int originalWidth;
  		int originalHeight;
  		int maxN;
  		originalWidth = before.getWidth();
  		originalHeight = before.getHeight();
  		// Padding required? Pad if necessary
  		maxN = Math.max(originalWidth, originalHeight);
  		int newDimension = FFTUtil.getNextPowerOfTwo(maxN);
  		double [] filterEstimate = new double [newDimension];
  		double [] fftBefore = new double [newDimension*2];
  		double [] fftAfter = new double [newDimension*2];
  		DoubleFFT_1D fft = new DoubleFFT_1D(newDimension);
  		double outlier = 0;
  		for (int j = 0;  j < before.getHeight(); j++){ // for all image rows
  			// Perform the fft transform on both images
  			for (int i = 0 ; i < newDimension; i++){
  				fftBefore[i*2] = before.getPixelValue(i, j);
  				fftBefore[(i*2)+1] = 0;
  				fftAfter[i*2] = after.getPixelValue(i, j);
  				fftAfter[(i*2)+1] = 0;
  			}
  			fft.complexForward(fftBefore);
  			fft.complexForward(fftAfter);
  			// estimate the filter
  			if (meanSquare) {
  				double [] input = new double [newDimension];
  				double [] output = new double [newDimension];
  				for (int i = 0 ; i < newDimension; i++){
  					double in = abs(i,fftBefore);
  					double out = abs(i, fftAfter);
  					// Just consider observations which are not 0;
  					if ((Math.abs(in) > 0.01)){
  						// Least Square estimate
  						input[i] += in * in;
  						output[i] += in * out;
  					}
  				}
  				for (int i = 0 ; i < newDimension; i++){
  					filterEstimate[i] = output[i] / input[i];
  					if ((filterEstimate[i] > threshold)||(Double.isNaN(filterEstimate[i]))) filterEstimate[i] = -1;
  				}
  			} else {
  				for (int i = 0 ; i < newDimension; i++){
  					double beforeValue = abs(i, fftBefore);
  					double add = 0;
  					if (beforeValue != 0) {
  						add = abs(0, divideComplex(fftAfter[2*i], fftAfter[(2*i)+1], fftBefore[2*i], fftBefore[(2*i)+1]));
  						if (add < threshold) {
  							filterEstimate[i] += add / before.getHeight();
  						} else {
  							outlier++;
  						}
  					} else {
  						outlier++;
  					}
  				}	
  			}
  		}
  		outlier /= before.getHeight(); // Outliers per row
  		if (!meanSquare) System.out.println("Outliers per row: " + outlier);
  		return filterEstimate;
  	}
  
  	/**
  	 * Computes the absolute value of the complex number at position pos in the array
  	 * @param pos the position
  	 * @param array the array which contains the values
  	 * @return the absolute value
  	 */
  	public static double abs (int pos, double[] array){
  		return Math.sqrt(Math.pow(array[pos *2], 2) + Math.pow(array[(2*pos)+1], 2));
  	}
  
  	/**
  	 * Computes the absolute values of the complex number at posx, posy in the 2D array array.
  	 * 
  	 * @param posx x position
  	 * @param posy y position
  	 * @param array the array
  	 * @return the absolute value of the complex number
  	 */
  	public static double abs (int posx, int posy, double[][] array){
  		return Math.sqrt(Math.pow(array[posx][posy *2], 2) + Math.pow(array[posx][(2*posy)+1], 2));
  	}
  
  	/**
  	 * Prints the absolute values in the given array
  	 * @param array the array of complex values
  	 */
  	public static void printAbsolute(double [] array){
  		System.out.println("Begin Absolute Out");
  		for (int i = 0; i < array.length / 2; i++){
  			System.out.println(FFTUtil.abs(i, array));
  		}
  		System.out.println("End Absolute Out");
  	}
  
  	/**
  	 * Prints the complex values in the given array.
  	 * @param array the array
  	 */
  	public static void printComplex(double [] array){
  		System.out.println("Begin Compex Out");
  		for (int i = 0; i < array.length / 2; i++){
  			System.out.println(array[i*2] + " " + array[(i*2)+1]);
  		}
  		System.out.println("End Compex Out");
  	}
  	
  
  
  }
  /*
   * Copyright (C) 2010-2014  Andreas Maier
   * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
  */