import ij.*;
import ij.plugin.*;
import ij.gui.*;
import ij.process.*;
import ij.measure.Calibration;

/** This plugin does MRI t2 calculations (s2 = k/ln(s1/s2)) on two images or stacks.
 *
 *   ERK (Eric Kischell, keesh@ieee.org) changes 11/27/99 
 *       8-bit enhancement 
 *   ERK changes 04/13/02
 *       Aded option to output to a new, third image/stack, rather
 *       than overwriting the second (S2) stack.
*/ 

public class MRI_t2_Calculator implements PlugIn {

    private static String title = "MRI_t2_Calculator";
    private static double k = 35.0;

    public void run(String arg) {
        int[] wList = WindowManager.getIDList();
        if (wList==null || wList.length<2) {
            IJ.showMessage(title, "There must be at least two windows open.");
            return;
        }
        String[] titles = new String[wList.length];
        for (int i=0; i<wList.length; i++) {
            ImagePlus imp = WindowManager.getImage(wList[i]);
            if (imp!=null)
                titles[i] = imp.getTitle();
            else
                titles[i] = "";
        }
        GenericDialog gd = new GenericDialog("MRI t2 Calculator");
        gd.addChoice("S1:", titles, titles[0]);
        gd.addChoice("S2:", titles, titles[1]);
        gd.addNumericField("K:", k, 1);
        gd.addMessage("S2 = K/ln(S1/S2)");
        gd.showDialog();
        if (gd.wasCanceled())
            return;
        int s1Index = gd.getNextChoiceIndex();
        int s2Index = gd.getNextChoiceIndex();
        k = gd.getNextNumber();
        ImagePlus i1 = WindowManager.getImage(wList[s1Index]);
        ImagePlus i2 = WindowManager.getImage(wList[s2Index]);
        int width  = i1.getWidth();
        int height = i1.getHeight();
                
        if (width!=i2.getWidth() || height!=i2.getHeight()) {
            IJ.showMessage(title, "Images must have the same dimensions");
            return;
        }

        int size1 = i1.getStackSize();
        int size2 = i2.getStackSize();
        if (size1!=size2) {
            IJ.showMessage(title, "Stacks must have the same number of slices");
            return;
        }
        ImagePlus i3 = IJ.createImage((int)k+"/ln(s1/s2)", "32-bit", width, height, size1);
        i3.setCalibration(i1.getCalibration());
        calculate(i1, i2, i3, k);
        i3.show();
     }
    
    /** Performs MRI t2 calculations (s2 = k/ln(i1/i2)) on two images or stacks. */
    public void calculate(ImagePlus i1, ImagePlus i2, ImagePlus i3, double k) {
        int width  = i1.getWidth();
        int height = i1.getHeight();
        int slices = i1.getStackSize();
        ImageStack stack1 = i1.getStack();
        ImageStack stack2 = i2.getStack();
        ImageStack stack3 = i3.getStack();
	 Calibration cal1 = i1.getCalibration();
	 Calibration cal2 = i2.getCalibration();
        double v1, v2;
        for (int n=1; n<=slices; n++) {
            ImageProcessor ip1 = stack1.getProcessor(n);
            ip1.setCalibrationTable(cal1.getCTable()); // not needed in  ImageJ 1.40a or later
            ImageProcessor ip2 = stack2.getProcessor(n);
            ip2.setCalibrationTable(cal2.getCTable()); // not needed in  ImageJ 1.40a or later
            ImageProcessor ip3 = stack3.getProcessor(n);
            for (int y=0; y<height; y++) {
                for (int x=0; x<width; x++) {
                    v1 = ip1.getPixelValue(x, y);
                    v2 = ip2.getPixelValue(x, y);
                    ip3.putPixelValue(x, y, k/Math.log(v1/v2));
               }
            }          
            IJ.showProgress((double)n/slices);
        }
    }

}