package ij.plugin;
import ij.*;
import ij.process.*;
import ij.gui.*;
import ij.measure.*;
import ij.plugin.ContrastEnhancer;
import ij.measure.Calibration;
import ij.util.Tools;
import java.awt.*;
import java.util.*;

/** 
This class implements the FFT, Inverse FFT and Redisplay Power Spectrum commands 
in the Process/FFT submenu. It is based on Arlo Reeves'  
Pascal implementation of the Fast Hartley Transform from NIH Image 
(http://imagej.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.

Version 2008-08-25 inverse transform: mask is always symmetrized
*/
public class FFT implements  PlugIn, Measurements {

    static boolean displayFFT = true;
    public static boolean displayRawPS;
    public static boolean displayFHT;
    public static boolean displayComplex;
    public static String fileName;
    
    private ImagePlus imp;
    private boolean padded;
    private int originalWidth;
    private int originalHeight;
    private int stackSize = 1;
    private int slice = 1;
    private boolean doFFT;

    public void run(String arg) {
        if (arg.equals("options")) {
            showDialog();
            if (doFFT)
            	arg="fft";
            else
            	return;
        }
        imp = IJ.getImage();
        if (arg.equals("fft") && imp.isComposite()) {
        	if (!GUI.showCompositeAdvisory(imp,"FFT"))
        		return;
        }
        if (arg.equals("redisplay"))
            {redisplayPowerSpectrum(); return;}
        if (arg.equals("swap"))
            {swapQuadrants(imp.getStack()); imp.updateAndDraw(); return;}
       if (arg.equals("inverse")) {
            if (imp.getTitle().startsWith("FHT of"))
                {doFHTInverseTransform(); return;}
            if (imp.getStackSize()==2)
                {doComplexInverseTransform(); return;}
        }
        ImageProcessor ip = imp.getProcessor();
        Object obj = imp.getProperty("FHT");
        FHT fht = (obj instanceof FHT)?(FHT)obj:null;
        stackSize = imp.getStackSize();
        boolean inverse;
        if (fht==null && arg.equals("inverse")) {
            IJ.error("FFT", "Frequency domain image required");
            return;
        }
        if (fht!=null) {
            inverse = true;
            imp.deleteRoi();
        } else {
            if (imp.getRoi()!=null)
                ip = ip.crop();
            fht = newFHT(ip);
            inverse = false;
        }
        if (inverse)
            doInverseTransform(fht);
        else {
            fileName = imp.getTitle();
            doForwardTransform(fht);   
        }    
        IJ.showProgress(1.0);
    }
    
    void doInverseTransform(FHT fht) {
        fht = fht.getCopy();
        doMasking(fht);
        showStatus("Inverse transform");
        fht.inverseTransform();
        if (fht.quadrantSwapNeeded)
            fht.swapQuadrants();
        fht.resetMinAndMax();
        ImageProcessor ip2 = fht;
        if (fht.originalWidth>0) {
            fht.setRoi(0, 0, fht.originalWidth, fht.originalHeight);
            ip2 = fht.crop();
        }
        int bitDepth = fht.originalBitDepth>0?fht.originalBitDepth:imp.getBitDepth();
        switch (bitDepth) {
            case 8: ip2 = ip2.convertToByte(false); break;
            case 16: ip2 = ip2.convertToShort(false); break;
            case 24:
                showStatus("Setting brightness");
                if (fht.rgb==null || ip2==null) {
                    IJ.error("FFT", "Unable to set brightness");
                    return;
                }
                ColorProcessor rgb = (ColorProcessor)fht.rgb.duplicate();
                rgb.setBrightness((FloatProcessor)ip2);
                ip2 = rgb; 
                fht.rgb = null;
                break;
            case 32: break;
        }
        if (bitDepth!=24 && fht.originalColorModel!=null)
            ip2.setColorModel(fht.originalColorModel);
        String title = imp.getTitle();
        if (title.startsWith("FFT of "))
            title = title.substring(7, title.length());
        ImagePlus imp2 = new ImagePlus("Inverse FFT of "+title, ip2);
        imp2.setCalibration(imp.getCalibration());
        imp2.show();
    }

    void doForwardTransform(FHT fht) {
        showStatus("Forward transform");
        fht.transform();
        showStatus("Calculating power spectrum");
        long t0 = System.currentTimeMillis();
        ImageProcessor ps = fht.getPowerSpectrum();
        if (!(displayFHT||displayComplex||displayRawPS))
        	displayFFT = true;
        if (displayFFT) {
            ImagePlus imp2 = new ImagePlus("FFT of "+imp.getTitle(), ps);
            imp2.show((System.currentTimeMillis()-t0)+" ms");
            imp2.setProperty("FHT", fht);
            imp2.setCalibration(imp.getCalibration());
            String properties = "Fast Hartley Transform\n";
            properties += "width: "+fht.originalWidth + "\n";
            properties += "height: "+fht.originalHeight + "\n";
            properties += "bitdepth: "+fht.originalBitDepth + "\n";
            imp2.setProperty("Info", properties);
        }
    }
    
    FHT newFHT(ImageProcessor ip) {
        FHT fht;
        if (ip instanceof ColorProcessor) {
            showStatus("Extracting brightness");
            ImageProcessor ip2 = ((ColorProcessor)ip).getBrightness();
            fht = new FHT(pad(ip2));
            fht.rgb = (ColorProcessor)ip.duplicate(); // save so we can later update the brightness
        } else
            fht = new FHT(pad(ip));
        if (padded) {
            fht.originalWidth = originalWidth;
            fht.originalHeight = originalHeight;
        }
        int bitDepth = imp.getBitDepth();
        fht.originalBitDepth = bitDepth;
        if (bitDepth!=24)
        	fht.originalColorModel = ip.getColorModel();
        return fht;
    }
    
    ImageProcessor pad(ImageProcessor ip) {
        originalWidth = ip.getWidth();
        originalHeight = ip.getHeight();
        int maxN = Math.max(originalWidth, originalHeight);
        int i = 2;
        while(i1)
            IJ.showStatus("FFT: " + slice+"/"+stackSize);
        else
            IJ.showStatus(msg);
    }
    
    void doMasking(FHT ip) {
        if (stackSize>1)
            return;
        float[] fht = (float[])ip.getPixels();
        ImageProcessor mask = imp.getProcessor();
        mask = mask.convertToByte(false);
        if (mask.getWidth()!=ip.getWidth() || mask.getHeight()!=ip.getHeight())
            return;
        ImageStatistics stats = ImageStatistics.getStatistics(mask, MIN_MAX, null);
        if (stats.histogram[0]==0 && stats.histogram[255]==0)
            return;
        boolean passMode = stats.histogram[255]!=0;
        IJ.showStatus("Masking: "+(passMode?"pass":"filter"));
        mask = mask.duplicate();
        if (passMode)
            changeValuesAndSymmetrize(mask, (byte)255, (byte)0); //0-254 become 0
        else
            changeValuesAndSymmetrize(mask, (byte)0, (byte)255); //1-255 become 255
        //long t0=System.currentTimeMillis();
        for (int i=0; i<3; i++)
            smooth(mask);
        //IJ.log("smoothing time:"+(System.currentTimeMillis()-t0));
        if (IJ.debugMode || IJ.altKeyDown())
        	new ImagePlus("mask", mask.duplicate()).show();
        ip.swapQuadrants(mask);
        byte[] maskPixels = (byte[])mask.getPixels();
        for (int i=0; i0) pixels[n*n-i] = v1;
                } else if (i