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

/** digital darkfield Tableau
  *
  * Converts a real imaginary fourier stack into a 
  * tiled array of darkfield (and one brightfield) images.
  *
  * This example adapted by pf from tutorial plugins named
  * ColorInverter & StackAverage found with the tutorial at
  * http://www.fh-hagenberg.at/mtd/depot/imaging/imagej
  * 
  */
public class ddfTableaux_ implements PlugInFilter { // , PlugIn

	protected ImagePlus imp; // defines an instance variable, protected means what?
	static double wrad = 32;	
	static boolean normalized = true;

	public int setup(String arg, ImagePlus imp) {
		if (arg.equals("about"))
			{showAbout(); return DONE;}
		this.imp=imp;
                return DOES_32+STACK_REQUIRED; // checks for stacked floats as input
	}

	public void run(ImageProcessor ip) { //line30

		// if (!IJ.showMessageWithCancel("Warning","Are you ready to over write the selected image?")) return;
		if (!showDialog()) return;
		
		ImageStack stack=imp.getStack(); // sets up the image stack for access
		// get width, height and the region of interest
		int w = ip.getWidth();     
		int h = ip.getHeight();    
		// int dimension = ip.getWidth()*ip.getHeight();	
		float[] pixels01 = (float[]) stack.getPixels(1); // gets the first slice pixel array
		float[] pixels02 = (float[]) stack.getPixels(2); // gets the 2nd slice pixel array //line40
	
		// create a scratch image with the same size
		ImagePlus oImage = NewImage.createFloatImage("Digital Darkfield Tableau",w,h,2,NewImage.FILL_BLACK);
		ImageStack oStack = oImage.getStack();
		// oStack.ConvertToGray32();
		float[] pixels01o = (float[]) oStack.getPixels(1);
		float[] pixels02o = (float[]) oStack.getPixels(2); 

		// now calculate offsets from dialog inputs        		
		// int xoff = (int) (roff*Math.cos(aoff*Math.PI/180.)); 
        		// int yoff = (int) (roff*Math.sin(aoff*Math.PI/180.));

		// 
		// Rectangle roi01;
		int xnew, ynew; 
		int oldj, newj;
		int roix=0, roiy=0, roiw=(int) wrad, roih=(int) wrad; // later prompt for roiw/roih
		int nxwin, nywin;
		nxwin = w/roiw; nywin = h/roih;
		// nxwin = 1; nywin = 1; // debug statement
		String s = imp.getTitle();

        		for (int y = 0; y < h; y++){ // load old image into the scratch image space
            		for (int x = 0; x < w; x++){
				oldj = x + w*y;
				pixels01o[oldj]=pixels01[oldj];
				pixels02o[oldj]=pixels02[oldj];
			}
		}

		int k, kroi, newx, newy;
		float scale=1;
		double mag;
                	Rectangle fitRect = new Rectangle();
                	fitRect.x = 0; fitRect.y=0; fitRect.width = roiw; fitRect.height = roih;
		ImageStack stackin = new ImageStack(roiw,roih);
		float[] rein = new float[roiw*roih];
		float[] imin = new float[roiw*roih];
		stackin.addSlice("Real", rein);
		stackin.addSlice("Imaginary", imin);
		float[] inrepix = (float[]) stackin.getPixels(1);
		float[] inimpix = (float[]) stackin.getPixels(2);
		ImageStack stackout = new ImageStack(roiw,roih);
		float[] reout = new float[roiw*roih];
		float[] imout = new float[roiw*roih];
		stackout.addSlice("Real", reout);
		stackout.addSlice("Imaginary", imout);
		float[] outrepix = (float[]) stackout.getPixels(1);
		float[] outimpix = (float[]) stackout.getPixels(2);
		int offset = roiw/2;
		for (int iw = 0; iw<nxwin; iw++){ // begin cycle through all input tiles iw,jw
			roix = iw*roiw;
			// fitRect.x = roix;
			for (int jw = 0; jw<nywin; jw++) {
				roiy = jw*roih;
				// load the window location in oStack (or stack) into our working tile
				for (int i = 0; i<roiw; i++) { // first load input tile from complex array
					// newx = roix+i;
					if(roix+i+offset<w) {newx=roix+i+offset;} else {newx=i+offset;}
					for (int j = 0; j<roih; j++) {
						// newy=roiy+j;
						if(roiy+j+offset<h) {newy=roiy+j+offset;} else {newy=j+offset;}
						kroi = i + roiw*j;
						k = newx + w*newy;
						inrepix[kroi]= pixels01o[k];
						inimpix[kroi] = pixels02o[k];
					}
				}
				// set up and execute inverse transform of the working tile
				swapQuadrants(stackin); // couldn't figure out how to call external
				rein = (float[])stackin.getPixels(1);
				imin = (float[])stackin.getPixels(2);
				c2c2DFFT(rein, imin, roiw, reout, imout); // likewise
				stackout.setPixels(reout,1);
				stackout.setPixels(imout,2);
				// return the output tile to current location in oStack
				// initialize contrast normalization if checked
				if(normalized) {
					scale = 0;
					for (int i = 0; i<roiw; i++) {
						for (int j=0; j<roih; j++) {
							kroi = i + roiw*j;
							mag = Math.sqrt(outrepix[kroi]*outrepix[kroi]+outimpix[kroi]*outimpix[kroi]);
							if (scale<mag) scale = (float) mag;
						}
					}
				} else {
					scale = 1;
				}
				for (int i = 0; i<roiw; i++) { // now return output tile to original array
					for (int j = 0; j<roih; j++) {
						kroi = i + roiw*j;
						k = (roix+i) + w*(roiy+j);
						// add contrast normalization if checked
						pixels01o[k] = outrepix[kroi]/scale;
						pixels02o[k] = outimpix[kroi]/scale;
					}
				}
			}
		}

     		oImage.getProcessor().resetMinAndMax();
		oImage.show();
		oImage.updateAndDraw();
	}

    	void swapQuadrants(ImageStack stack) {
        		FHT fht = new FHT(new FloatProcessor(1, 1));
        		for (int i=1; i<=stack.getSize(); i++) {
            		fht.swapQuadrants(stack.getProcessor(i));
		}
    	}

    	/** Complex to Complex Inverse Fourier Transform
    	*   @author Joachim Wesner
    	*/
   	void c2c2DFFT(float[] rein, float[] imin, int maxN, float[] reout, float[] imout) {
            	FHT fht = new FHT(new FloatProcessor(maxN,maxN));
            	float[] fhtpixels = (float[])fht.getPixels();
            	// Real part of inverse transform
            	for (int iy = 0; iy < maxN; iy++)
                  		cplxFHT(iy, maxN, rein, imin, false, fhtpixels);
            	fht.inverseTransform();
            	// Save intermediate result, so we can do a "in-place" transform
            	float[] hlp = new float[maxN*maxN];
            	System.arraycopy(fhtpixels, 0, hlp, 0, maxN*maxN);
            	// Imaginary part of inverse transform
            	for (int iy = 0; iy < maxN; iy++)
                  		cplxFHT(iy, maxN, rein, imin, true, fhtpixels);
            	fht.inverseTransform();
            	System.arraycopy(hlp, 0, reout, 0, maxN*maxN);
            	System.arraycopy(fhtpixels, 0, imout, 0, maxN*maxN);
      	}

    	/** Build FHT input for equivalent inverse FFT
    	*   @author Joachim Wesner
    	*/
    	void cplxFHT(int row, int maxN, float[] re, float[] im, boolean reim, float[] fht) {
            	int base = row*maxN;
            	int offs = ((maxN-row)%maxN) * maxN;
            	if (!reim) {
                  		for (int c=0; c<maxN; c++) {
                        		int l =  offs + (maxN-c)%maxN;
                        		fht[base+c] = ((re[base+c]+re[l]) - (im[base+c]-im[l]))*0.5f;
                  		}
            	} else {
                  		for (int c=0; c<maxN; c++) {
                        		int l = offs + (maxN-c)%maxN;
                        		fht[base+c] = ((im[base+c]+im[l]) + (re[base+c]-re[l]))*0.5f;
                  		}
            	}
      	}

	void showAbout() {
		// called by setup if the string argument is "about"
		IJ.showMessage("offsetXYstack",
			"offsets pixels in an XY stack by integer amounts"
		);
	}

	public boolean showDialog() {
		GenericDialog gd = new GenericDialog("Parameters");
		// gd.addNumericField("horizontal offset frequency in fpixels:", xoff, 64);
		// gd.addNumericField("vertical offset frequency in fpixels:", yoff, 64);
		gd.addNumericField("tile width in fpixels:", wrad, 32);
		gd.addCheckbox("Normalize Tile Contrast?", normalized);
		gd.showDialog();
		if (gd.wasCanceled()) return false;
		// xoff = gd.getNextNumber();
		// yoff = gd.getNextNumber();
		wrad = gd.getNextNumber();
		normalized = gd.getNextBoolean();
        		return true;
	}

}