import ij.*;
import ij.process.*;
import ij.gui.*;
import java.awt.*;
import ij.plugin.PlugIn;
import ij.WindowManager;
//import edu.emory.mathcs.jtransforms.fft.*;
//import edu.emory.mathcs.utils.*;
import org.jtransforms.fft.*;
import org.jtransforms.utils.*;
import flanagan.complex.*;
import flanagan.math.*;
import ij.plugin.frame.RoiManager;
import ij.gui.Roi;


public class Align_ComputeShifts implements PlugIn {
	protected ImagePlus imp_r, imp_i;
	protected int reference_slide;
	protected Roi roi;
	protected int usfac;
	protected boolean debug = true;

	public void run(String arg) {
		int[] wList = WindowManager.getIDList();

		if (null == wList) { //there are no images
			IJ.error("There must be at least one open image");
			return;
		}

		//get all the image titles so they can be shown in the dialog
		String[] titles = new String[wList.length+1];
		ImagePlus[] limps = new ImagePlus[wList.length+1];

		titles[0] = "None";
		limps[0] = null;
		for (int i=0, k=1; i<wList.length; i++) {
			ImagePlus limp = WindowManager.getImage(wList[i]);
			if (limp != null) {
				titles[k] = limp.getTitle();
				limps[k] = limp;
				k++;
			}
		}

		GenericDialog gd = new GenericDialog("Alignement Procedure");
		String[] tasks = {"Calculate direct shifts", "Calculate all shifts"};
		gd.addChoice("Task:", tasks, tasks[0]);
		gd.addChoice("Real part image:", titles, titles[1]);
		gd.addChoice("Imag part image (optional):", titles, titles[0]);
		gd.addNumericField("Upscale factor (1 < integer < 100)", 100, 0);
		
		gd.showDialog();
      		if (gd.wasCanceled()) return;
      		
      		this.usfac = (int)Math.round(gd.getNextNumber());
		int task =  gd.getNextChoiceIndex();
		int i1Index = gd.getNextChoiceIndex();
		int i2Index = gd.getNextChoiceIndex();

		//get the images stack(s)
		if (i1Index == 0) {
			IJ.error("Real part image must exist!");
			return;
		}

		imp_r = limps[i1Index];
		if (i2Index == 0) {
			imp_i = null;
		} else {
			imp_i = limps[i2Index];
		}

		//get the reference slide
		reference_slide = imp_r.getCurrentSlice();

		// get the ROI
		RoiManager manager = RoiManager.getInstance();
		if (manager == null) {
			IJ.error("No ROI manager opened");
			return;
		}

		Roi rois[] = manager.getSelectedRoisAsArray();
		if (rois.length != 1) {
			IJ.error("There must be one and only one ROI selected");
			return;
		}
		
		roi = rois[0];
		if (roi.getType()!=Roi.RECTANGLE) {
			IJ.error("The selected ROI must be a rectangle");
			return;
		}
		
		Rectangle box = roi.getBounds();
		if (!ConcurrencyUtils.isPowerOf2(box.height)
			|| !ConcurrencyUtils.isPowerOf2(box.width)) {
			IJ.error("The selected ROI height and with must be a power of 2");
			return;
		}
		
      		if (task == 0) {
      			calculateShiftsRun();
      		} else {
      			calculateAllShiftsRun();
      		}
      		
      		return;
	}

	private void calculateShiftsRun() {
		// perform the FFT of each slice
        
		IJ.showStatus("1/2 Perform FFT of each slice");
		ComplexMatrix[] ffts = computeFFT();
        
		// calculate shifts
		IJ.showStatus("2/2 Calculate shifts between slices");
		double[][] shifts = calculateShifts(ffts);
		
		// save shifts
		ShiftsIO sio = new ShiftsIO();
		sio.save(shifts, "directshifts");
       
	}

	private void calculateAllShiftsRun() {
		// perform the FFT of each slice
		IJ.showStatus("1/2 Perform FFT of each slice");
		ComplexMatrix[] ffts = computeFFT();

		// calculate shifts
		IJ.showStatus("2/2 Calculate shifts between slices");
		double[][] shifts = calculateAllShifts(ffts);
		// save shifts
		ShiftsIO sio = new ShiftsIO();
		sio.save(shifts, "allshifts");
	}

	// perform the FFT of each slice
	private ComplexMatrix[] computeFFT() {
		
		int slices = imp_r.getStackSize();
		
		ComplexMatrix[] ffts = new ComplexMatrix[slices];
		for (int i=1; i <= slices; i++) {
			if (imp_i == null) {
				ImageProcessor ip = imp_r.getStack().getProcessor(i);
				ip.setRoi(roi);
				ImageProcessor curr = ip.crop().convertToFloat();
				double[][] data = ImageProcessor_to_FFTArray2D(curr);	
				ffts[i-1] = fft2(data);
			} else {
				ImageProcessor ip1, ip2;
				ip1 = imp_r.getStack().getProcessor(i);
				ip1.setRoi(roi);
				ImageProcessor curr_r = ip1.crop().convertToFloat();
				ip2 = imp_i.getStack().getProcessor(i);
				ip2.setRoi(roi);
				ImageProcessor curr_i = ip2.crop().convertToFloat();
				double[][] data = ImageProcessor_to_FFTComplexArray2D(curr_r, curr_i);
				ffts[i-1] = cfft2(data);
			}
			IJ.showProgress(i, slices);
		}

		return ffts;
	}

	//calculate the shifts between ffts
	private double[][] calculateShifts(ComplexMatrix[] ffts) {
        
		double[][] shifts = new double[ffts.length][6];        
		for (int i = 0; i < ffts.length; i++) {
			shifts[i][0] = reference_slide; shifts[i][1] = i+1;
			double[] temp = DFTRegistration(ffts[reference_slide - 1], ffts[i]);
			shifts[i][2] = temp[2]; shifts[i][3] = temp[3];
			shifts[i][4] = temp[0]; shifts[i][5] = temp[1];            
			IJ.showProgress(i + 1, ffts.length);
		}
		return shifts; // [ref, drifted, dr, dc, error, diffphase]
	
	}

	//calculate all the shifts between ffts
	private double[][] calculateAllShifts(ComplexMatrix[] ffts) {

		double[][] shifts = new double[ffts.length*(ffts.length-1)/2][6];
		int id = 0;
		for (int i = 0; i < ffts.length-1; i++) {
			for (int j = i+1; j < ffts.length; j++) {
				shifts[id][0] = i+1; shifts[id][1] = j+1;
				double[] temp = DFTRegistration(ffts[i], ffts[j]);
				shifts[id][2] = temp[2]; shifts[id][3] = temp[3];
				shifts[id][4] = temp[0]; shifts[id][5] = temp[1];
				id = id + 1;
				IJ.showProgress(id + 1, ffts.length*(ffts.length-1)/2);
			}
		}
		
		return shifts; // [ref,drifted,dr,dc,error, diffphase]
	}
	
	// compute 2D fft from an image
	private ComplexMatrix fft2(double[][] data) {

		int h = data.length;
		int w = data[0].length;
		DoubleFFT_2D fft = new DoubleFFT_2D(h, w);
		
		fft.realForward(data);
		ComplexMatrix m = FFTArray2D_to_ComplexMatrix(data, h, w);

		return m;
	}

	// compute complex 2D fft from an image
	private ComplexMatrix cfft2(double[][] data) {

		int h = data.length;
		int w = data[0].length;
		DoubleFFT_2D fft = new DoubleFFT_2D(h, w/2);
		
		fft.complexForward(data);
		ComplexMatrix m = FFTComplexArray2D_to_ComplexMatrix(data, h, w/2);

		return m;
	}
	
	// compute inverse 2D fft from a complex matrix
	private double[][] ifft2(ComplexMatrix m) {
		int w = m.getNcol();
		int h = m.getNrow();
		
		DoubleFFT_2D fft = new DoubleFFT_2D(h, w);
		
		double[][] data = ComplexMatrix_to_FFTArray2D(m);	
		fft.realInverse(data, true);

		return data;
	}
	
	// compute complex inverse 2D fft from a complex matrix
	private ComplexMatrix cifft2(ComplexMatrix m) {
		int w = m.getNcol();
		int h = m.getNrow();
		
		DoubleFFT_2D fft = new DoubleFFT_2D(h, w);
		
		double[][] data = new double[h][2*w];
		for (int j=0; j<h; j++) {
			for (int i=0; i<w; i++) {
				data[j][2*i] = m.getElementReference(j, i).getReal();
				data[j][2*i+1] = m.getElementReference(j, i).getImag();				
			}
		}

		fft.complexInverse(data, true);

		ComplexMatrix out = new ComplexMatrix(h, w);
		for (int j=0; j<h; j++) {
			for (int i=0; i<w; i++) {
				out.setElement(j, i, data[j][2*i], data[j][2*i+1]);			
			}
		}

		return out;
	}
	
	private double[] DFTRegistration(ComplexMatrix ref, ComplexMatrix drifted) {
		int m = ref.getNrow();
		int n = ref.getNcol();
		double[] output = new double[4];
		
		// First upsample by a factor of 2 to obtain initial estimate
    		// Embed Fourier data in a 2x larger array
    		int mlarge = m*2;
    		int nlarge = n*2;
    		ComplexMatrix large = new ComplexMatrix(mlarge, nlarge);
    		ComplexMatrix c = fftshift(ElementProduct(ref, drifted.conjugate()));

    		for (int j = 0; j < m; j++) {
    			for (int i = 0; i < n; i++) {
    				large.setElement((int)(j + m - Math.floor(m/2.0)), (int)(i + n - Math.floor(n/2.0)), c.getElementReference(j, i));
    			}
    		}
  
    		// Compute crosscorrelation and locate the peak
    		ComplexMatrix CC = cifft2(ifftshift(large));
    		double[] peak = cFindPeak(CC); //max, r, c, max_r, max_c
    		    
    		// Obtain shift in original pixel grid from the position of the
    		// crosscorrelation peak 
    		if (peak[1] > m) { 
        		peak[1] = peak[1] - mlarge;
    		}
		if (peak[2] > n) {
			peak[2] = peak[2] - nlarge;
		}

    		//% If upsampling > 2, then refine estimate with matrix multiply DFT
    		if (this.usfac > 2) {
     			//   %%% DFT computation %%%
     			//   % Initial shift estimate in upsampled grid
     			double row_shift = Math.round(peak[1]/2.0*this.usfac)/this.usfac; 
     			double col_shift = Math.round(peak[2]/2.0*this.usfac)/this.usfac;     
     			int dftshift = (int)Math.floor(Math.ceil(this.usfac*1.5)/2); // Center of output array at dftshift+1
     			//   % Matrix multiply DFT around the current shift estimate
    			ComplexMatrix in = ElementProduct(drifted, ref.conjugate());
     			ComplexMatrix nCC = dftups(in, (int)Math.ceil(this.usfac*1.5), (int)Math.ceil(this.usfac*1.5),
     					dftshift-row_shift*this.usfac, dftshift-col_shift*this.usfac);
     		     	nCC = nCC.times(1.0/(m*n*this.usfac*this.usfac)).conjugate();
     			//   % Locate maximum and map back to original pixel grid 
     			double[] npeak = cFindPeak(nCC); //max_r, max_i, r, c
    
 			ComplexMatrix mrg00 = dftups(ElementProduct(ref, ref.conjugate()),1,1,0,0);
 			double rg00 = mrg00.getElementReference(0, 0).abs()/(m*n*this.usfac*this.usfac);
     			ComplexMatrix mrf00 = dftups(ElementProduct(drifted, drifted.conjugate()),1,1,0,0);
     			double rf00 = mrf00.getElementReference(0, 0).abs()/(m*n*this.usfac*this.usfac);
     			
     			npeak[1] = npeak[1] - dftshift;
     			npeak[2] = npeak[2] - dftshift;
     			output[0] = Math.sqrt(Math.abs(1.0 - npeak[0]*npeak[0]/(rg00*rf00)));	//error
			output[1] = Math.atan2(npeak[4], npeak[3]);				//diffphase		
     			output[2] = row_shift + npeak[1]/this.usfac;			//delta row
     			output[3] = col_shift + npeak[2]/this.usfac;			//delta col
     			
    		} else {
   			// % If upsampling = 2, no additional pixel shift refinement
   			double rg00 = SumSquareAbs(ref)/(mlarge*nlarge);
			double rf00 = SumSquareAbs(drifted)/(mlarge*nlarge);
			
			output[0] = Math.sqrt(Math.abs(1.0 - peak[0]*peak[0]/(rg00*rf00)));	//error
			output[1] = Math.atan2(peak[4], peak[3]);				//diffphase		
			output[2] = peak[1]/2.0;					//delta row
    			output[3] = peak[2]/2.0;					//delta col
    		}
    		
		return output;	
	}

	private double SumSquareAbs(ComplexMatrix m) {
		double sum = 0.0;

		for (int j = 0; j < m.getNrow(); j ++){
			for (int i = 0; i < m.getNcol(); i++) {
				sum += m.getElementReference(j, i).squareAbs();
			}
		}

		return sum;
	}
	
	private double[] cFindPeak(ComplexMatrix m) {
		double max = 0.0;
		double realmax = 0.0;
		double imagmax = 0.0;
		int cmax = 0, rmax = 0;
		
		for (int j = 0; j < m.getNrow(); j ++){
			for (int i = 0; i < m.getNcol(); i++) {
				if (m.getElementReference(j, i).abs() > max) {
					max = m.getElementReference(j, i).abs();
					realmax = m.getElementReference(j, i).getReal();
					imagmax = m.getElementReference(j, i).getImag();
					rmax = j;
					cmax = i;
				}
			}
		}

		double[] res = new double[5];
		res[0] = Math.sqrt(realmax*realmax+imagmax*imagmax); res[1] = rmax; res[2] = cmax;
		res[3] = realmax; res[4] = imagmax;
		return res;
	}
	
	private ComplexMatrix fftshift(ComplexMatrix in) {
		int nc = in.getNcol();
		int nr = in.getNrow();

		ComplexMatrix out = new ComplexMatrix (nr, nc);
		
		int midi = (int)Math.floor(nc/2.0);
		int offi = (int)Math.ceil(nc/2.0);
		int midj = (int)Math.floor(nr/2.0);
		int offj = (int)Math.ceil(nr/2.0);
		
		for (int j = 0; j < nr; j ++){
			for (int i = 0; i < nc; i++) {
				if (j < midj) {
					if (i < midi) {
						out.setElement(j, i, in.getElementReference(j+offj, i+offi));
					} else {
						out.setElement(j, i, in.getElementReference(j+offj, i-midi));		
					}
				} else {
					if (i < midi) {
						out.setElement(j, i, in.getElementReference(j-midj, i+offi));
					} else {
						out.setElement(j, i, in.getElementReference(j-midj, i-midi));		
					}
				}
			}
		}

		return out;
	}

	private ComplexMatrix ifftshift(ComplexMatrix in) {
		int nc = in.getNcol();
		int nr = in.getNrow();

		ComplexMatrix out = new ComplexMatrix (nr, nc);
		
		int midi = (int)Math.ceil(nc/2.0);
		int offi = (int)Math.floor(nc/2.0);
		int midj = (int)Math.ceil(nr/2.0);
		int offj = (int)Math.floor(nr/2.0);
		
		for (int j = 0; j < nr; j ++){
			for (int i = 0; i < nc; i++) {
				if (j < midj) {
					if (i < midi) {
						out.setElement(j, i, in.getElementReference(j+offj, i+offi));
					} else {
						out.setElement(j, i, in.getElementReference(j+offj, i-midi));		
					}
				} else {
					if (i < midi) {
						out.setElement(j, i, in.getElementReference(j-midj, i+offi));
					} else {
						out.setElement(j, i, in.getElementReference(j-midj, i-midi));		
					}
				}
			}
		}

		return out;
	}
	
	private Matrix ifftshift(Matrix in) {
		int nc = in.getNcol();
		int nr = in.getNrow();

		Matrix out = new Matrix (nr, nc);
		
		int midi = (int)Math.ceil(nc/2.0);
		int offi = (int)Math.floor(nc/2.0);
		int midj = (int)Math.ceil(nr/2.0);
		int offj = (int)Math.floor(nr/2.0);
		
		for (int j = 0; j < nr; j ++){
			for (int i = 0; i < nc; i++) {
				if (j < midj) {
					if (i < midi) {
						out.setElement(j, i, in.getElement(j+offj, i+offi));
					} else {
						out.setElement(j, i, in.getElement(j+offj, i-midi));		
					}
				} else {
					if (i < midi) {
						out.setElement(j, i, in.getElement(j-midj, i+offi));
					} else {
						out.setElement(j, i, in.getElement(j-midj, i-midi));		
					}
				}
			}
		}

		return out;
	}
	
	private ComplexMatrix dftups(ComplexMatrix in, int nor, int noc, double roff, double coff) {
		// function out=dftups(in,nor,noc,usfac,roff,coff);
		// Upsampled DFT by matrix multiplies, can compute an upsampled DFT in just
		// a small region.
		// usfac         Upsampling factor (default usfac = 1)
		// [nor,noc]     Number of pixels in the output upsampled DFT, in
		//               units of upsampled pixels (default = size(in))
		// roff, coff    Row and column offsets, allow to shift the output array to
		//               a region of interest on the DFT (default = 0)
		// Recieves DC in upper left corner, image center must be in (1,1)
		// Loïc Le Guyader - Jun 11, 2011 Java version for ImageJ plugin 
		// Manuel Guizar - Dec 13, 2007
		// Modified from dftus, by J.R. Fienup 7/31/06

		// This code is intended to provide the same result as if the following
		// operations were performed
		//   - Embed the array "in" in an array that is usfac times larger in each
		//     dimension. ifftshift to bring the center of the image to (1,1).
		//   - Take the FFT of the larger array
		//   - Extract an [nor, noc] region of the result. Starting with the 
		//     [roff+1 coff+1] element.

		// It achieves this result by computing the DFT in the output array without
		// the need to zeropad. Much faster and memory efficient than the
		// zero-padded FFT approach if [nor noc] are much smaller than [nr*usfac nc*usfac]

		
		int nr = in.getNrow();
		int nc = in.getNcol();
	
		// Compute kernels and obtain DFT by matrix products
		double amplitude = -2.0*Math.PI/(nc*usfac);
		
		Matrix u = new Matrix(nc, 1);
		for (int i = 0; i < nc; i++) {
			u.setElement(i, 0, i - Math.floor(nc/2.0));
		}
		u = ifftshift(u);

		Matrix v = new Matrix(1, noc);
		for (int i = 0; i < noc; i++) {
			v.setElement(0, i, i-coff);
		}

		Matrix phase = u.times(v);
		ComplexMatrix kernc = new ComplexMatrix(nc, noc);
		for (int j = 0; j < nc; j++) {
			for (int i = 0; i < noc; i++) {
				Complex t = new Complex();
				t.polar(1.0, amplitude*phase.getElement(j, i));
				kernc.setElement(j, i, t);
			}		
		}
		//ComplexMatrixPrint(kernc);
		
		amplitude = -2.0*Math.PI/(nr*usfac);
		
		Matrix w = new Matrix(nor, 1); 
		for (int i = 0; i < nor; i++) {
			w.setElement(i, 0, i - roff);
		}
		
		Matrix x = new Matrix(1, nr);
		for (int i = 0; i < nr; i++) {
			x.setElement(0, i, i - Math.floor(nr/2.0));
		}
		x = ifftshift(x);

		Matrix nphase = w.times(x);
		ComplexMatrix kernr = new ComplexMatrix(nor, nr);
		for (int j = 0; j < nor; j++) {
			for (int i = 0; i < nr; i++) {
				Complex t = new Complex();
				t.polar(1.0, amplitude*nphase.getElement(j, i));
				kernr.setElement(j, i, t);
			}		
		}
		//ComplexMatrixPrint(kernr);
		
		ComplexMatrix out = kernr.times(in.times(kernc));

		return out;
	}
	
	private double[][] CrossCorrelation(ComplexMatrix ref, ComplexMatrix drifted) {
		int h = ref.getNrow();
		int w = ref.getNcol();
		ComplexMatrix b = drifted.conjugate();
		ComplexMatrix res = ElementProduct(ref, b);
		
		double[][] data = ComplexMatrix_to_FFTArray2D(res);
		DoubleFFT_2D fft = new DoubleFFT_2D(h, w);
		fft.realInverse(data, true);
		
		return data;
	}

	private ComplexMatrix ElementProduct(ComplexMatrix a, ComplexMatrix b) {
		int nr = a.getNrow();
		int nc = a.getNcol();
		
		ComplexMatrix res = new ComplexMatrix(nr, nc);

		for(int j = 0; j < nr; j++) {
			for(int i = 0; i < nc; i++) {
				res.setElement(j, i, a.getElementReference(j, i).times(b.getElementReference(j, i)));
			}
		}

		return res;
	}
	
	private double[][] ImageProcessor_to_FFTArray2D(ImageProcessor ip) {

		float[] pixels = (float[])ip.getPixels();
		int w = ip.getWidth();
		int h = ip.getHeight();
		double[][] data = new double[h][w];
		
		for (int j = 0; j < h; j++) { 
			for (int i = 0; i < w; i++) {
				data[j][i] = (double)pixels[j*w + i];
    			}
		}

		return data;
	}
	
	private double[][] ImageProcessor_to_FFTComplexArray2D(ImageProcessor ip_r, ImageProcessor ip_i) {

		float[] pixels_r = (float[])ip_r.getPixels();
		float[] pixels_i = (float[])ip_i.getPixels();
		int w = ip_r.getWidth();
		int h = ip_r.getHeight();
		double[][] data = new double[h][2*w];
		
		for (int j = 0; j < h; j++) { 
			for (int i = 0; i < w; i++) {
				data[j][2*i] = (double)pixels_r[j*w + i];
				data[j][2*i+1] = (double)pixels_i[j*w + i];
    			}
		}

		return data;
	}
	
	private ComplexMatrix FFTArray2D_to_ComplexMatrix(double[][] data, int h, int w) {
	
		ComplexMatrix m = new ComplexMatrix(h,w);

		for (int j = 0; j < h; j++) {
			for (int i = 0; i <= w/2; i++) {
				if (j > 0 && i > 0 && i < w/2) { 
      					m.setElement(j,i, new Complex(data[j][2*i], data[j][2*i+1]));
      					m.setElement(h-j, w-i, new Complex(data[j][2*i], -data[j][2*i+1]));
				}
				if (j == 0 && i > 0 && i < w/2) {
					m.setElement(0, i, new Complex(data[0][2*i], data[0][2*i+1]));
      					m.setElement(0, w-i, new Complex(data[0][2*i], -data[0][2*i+1]));
				}
				if (i == 0 && j > 0 && j < h/2) {
					m.setElement(j,0, new Complex(data[j][0], data[j][1]));
      					m.setElement(h-j, 0, new Complex(data[j][0], -data[j][1]));
      					m.setElement(j, w/2, new Complex(data[h-j][1], -data[h-j][0]));
      					m.setElement(h-j, w/2, new Complex(data[h-j][1], data[h-j][0]));	
				}
				if (j == 0 && i == 0) {
					m.setElement(0, 0, new Complex(data[0][0], 0));
				}
				if (j == 0 && i == w/2) {
					m.setElement(0, w/2, new Complex(data[0][1], 0));
				}
				if (j == h/2 && i == 0) {
					m.setElement(h/2, 0, new Complex(data[h/2][0], 0));
				}
				if (j == h/2 && i == w/2) {
					m.setElement(h/2, w/2, new Complex(data[h/2][1], 0));
				}
			}
		}

		return m;
	}

	private ComplexMatrix FFTComplexArray2D_to_ComplexMatrix(double[][] data, int h, int w) {
	
		ComplexMatrix m = new ComplexMatrix(h,w);

		for (int j = 0; j < h; j++) {
			for (int i = 0; i < w; i++) {
				m.setElement(j,i, new Complex(data[j][2*i], data[j][2*i+1]));
			}
		}

		return m;
	}

	private double[][] ComplexMatrix_to_FFTArray2D(ComplexMatrix m) {
		int w = m.getNcol();
		int h = m.getNrow();
		double[][] data = new double[h][w];

		for (int j = 0; j < h; j++) {
			for (int i = 0; i <= w/2; i++) {
				if (j > 0 && i > 0 && i < w/2) {
      					data[j][2*i] = m.getElementReference(j,i).getReal();
      					data[j][2*i+1] = m.getElementReference(j,i).getImag();
				}
				if (j == 0 && i > 0 && i < w/2) {
					data[0][2*i] = m.getElementReference(0,i).getReal();
      					data[0][2*i+1] = m.getElementReference(0,i).getImag();
				}
				if (i == 0 && j > 0 && j < h/2) {
					data[j][0] = m.getElementReference(j,0).getReal();
      					data[j][1] = m.getElementReference(j,0).getImag();
      					data[h-j][1] = m.getElementReference(j,w/2).getReal();
      					data[h-j][0] = m.getElementReference(h-j,w/2).getImag();	
				}
				if (j == 0 && i == 0) {
					data[0][0] = m.getElementReference(0,0).getReal();
				}
				if (j == 0 && i == w/2) {
					data[0][1] = m.getElementReference(0,w/2).getReal();
				}
				if (j == h/2 && i == 0) {
					data[h/2][0] = m.getElementReference(h/2,0).getReal();
				}
				if (j == h/2 && i == w/2) {
					data[h/2][1] = m.getElementReference(h/2,w/2).getReal();
				}
			}
		}

		return data;
	}

	// convert a Complex Matrix into an 2d real part array data[0][][]
	// and 2d imaginary part data[1][][]
	private double[][][] ComplexMatrix_to_RealArray2D(ComplexMatrix m) {
		int w = m.getNcol();
		int h = m.getNrow();
		double[][][] data = new double[2][h][w];

		for (int j = 0; j < h; j++) {
			for (int i = 0; i < w; i++) {
				data[0][j][i] = m.getElementReference(j,i).getReal();
      				data[1][j][i] = m.getElementReference(j,i).getImag();
			}
		}
		
		return data;
	}
}