import ij.*;
import ij.plugin.PlugIn;
import ij.gui.*;
import ij.process.*;
/* Bob Dougherty, OptiNav, Inc.  Plugin to compute the 2D point spread function of a diffraction limited
microscope.  See http://yakko.bme.virginia.edu/phy506/.
Version 0	May25, 2004
*/
/*	License:
	Copyright (c) 2004, 2005, OptiNav, Inc.
	All rights reserved.

	Redistribution and use in source and binary forms, with or without
	modification, are permitted provided that the following conditions
	are met:

		Redistributions of source code must retain the above copyright
	notice, this list of conditions and the following disclaimer.
		Redistributions in binary form must reproduce the above copyright
	notice, this list of conditions and the following disclaimer in the
	documentation and/or other materials provided with the distribution.
		Neither the name of OptiNav, Inc. nor the names of its contributors
	may be used to endorse or promote products derived from this software
	without specific prior written permission.

	THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
	"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
	LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
	A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
	CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
	EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
	PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
	PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
	LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
	NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
	SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
public class Diffraction_Limit_PSF implements PlugIn {
	//Choices for normalization
    private static final int PEAK1=0, PEAK255=1, AREA1=2;
   	private static String[] norm = {"Peak = 1", "Peak = 255", "Sum of pixel values = 1" };
	//Constants for Bessel function approximation.
	private static double[] p = new double[]{
		-2.35619449,
		0.12499612,
		0.00005650,
		-0.00637879,
		0.00074348,
		0.00079824,
		-0.00029166};
	private static double[] f = new double[]{
		0.79788456,
		0.00000156,
		0.01659667,
		0.00017105,
		-0.00249511,
		0.00113653,
		-0.00020033};
	private static double[] t = new double[]{
		0.5,
		-0.56249985,
		0.21093573,
		-0.03954289,
		0.00443319,
		-0.00031761,
		0.00001109};

	public void run(String arg) {
		if (IJ.versionLessThan("1.32c"))
			return;

		double lambda = Prefs.get("difflimitpsf.lambda", 510);
		double pixelSpacing = Prefs.get("difflimitpsf.pixelSpacing", 30);
		double na = Prefs.get("difflimitpsf.na", 0.6);
		int w = (int)Prefs.get("difflimitpsf.w", 256);
		int h = (int)Prefs.get("difflimitpsf.h", 256);
		int normalization = (int)Prefs.get("difflimitpsf.normalization", 2);
		String title = Prefs.get("difflimitpsf.title", "PSF");
		boolean dB = Prefs.get("difflimitpsf.dB", false);

		GenericDialog gd = new GenericDialog("Specify psf", IJ.getInstance());
		gd.addMessage("Rayleigh resolution: 0.6*lambda/NA");
		gd.addNumericField("Wavelength (perhaps in nm)", lambda, 1);
		gd.addNumericField("Image pixel spacing, same units (ccd cell spacing / magnification)", pixelSpacing, 0);
		gd.addNumericField("Numerical Aperture, n*sin(theta)", na, 2);
        gd.addNumericField("PSF width, pixels",w,0);
        gd.addNumericField("PSF height, pixels",h,0);
 		gd.addChoice("Normalization", norm, norm[normalization]);
		gd.addStringField("Title", title, 12);
		gd.addCheckbox("PSF in dB",dB);

        gd.showDialog();
        if (gd.wasCanceled())
            return;

  		lambda = gd.getNextNumber();
		pixelSpacing = gd.getNextNumber();
		na = gd.getNextNumber();
		w = (int)gd.getNextNumber();
		h = (int)gd.getNextNumber();
		normalization = gd.getNextChoiceIndex();
		title = gd.getNextString();
		dB = gd.getNextBoolean();

		Prefs.set("difflimitpsf.lambda", lambda);
		Prefs.set("difflimitpsf.pixelSpacing", pixelSpacing);
		Prefs.set("difflimitpsf.na", na);
		Prefs.set("difflimitpsf.w", w);
		Prefs.set("difflimitpsf.h", h);
		Prefs.set("difflimitpsf.normalization", normalization);
		Prefs.set("difflimitpsf.title", title);
		Prefs.set("difflimitpsf.dB", dB);

		FloatProcessor fp = new FloatProcessor(w,h);
		float[] pixels = (float[])fp.getPixels();
		int ic = w/2;
		int jc = h/2;
		float a = (float)(2*Math.PI*na/lambda);

		double d = 0.6*lambda/(pixelSpacing*na);
		if (!IJ.showMessageWithCancel("PSF: peak to first dark ring","Rayleigh resolution = "+IJ.d2s(d)+" pixels"))
			return;
		IJ.showStatus("Computing psf");
		for (int j = 0; j < h; j++){
			IJ.showProgress((float)j/h);
			for (int i = 0; i < w; i++){
				float r = (float)Math.sqrt((i - ic)*(i-ic) + (j-jc)*(j - jc));
				r *= pixelSpacing;
				pixels[i + w*j] = psf(r,a);
			}
		}
		int n = w*h;
		float peak = pixels[ic + w*jc];
		if(normalization == PEAK1){
			float f = 1/peak;
			for (int ind = 0; ind < n; ind++){
				pixels[ind] *= f;
				if(pixels[ind] > 255)pixels[ind] = 255;
			}
		}else if(normalization == PEAK255){
			float f = 255/peak;
			for (int ind = 0; ind < n; ind++){
				pixels[ind] *= f;
				if(pixels[ind] > 255)pixels[ind] = 255;
			}
		}else if(normalization == AREA1){
			float area = 0;
			for (int ind = 0; ind < n; ind++){
				area += pixels[ind];
			}
			for (int ind = 0; ind < n; ind++){
				pixels[ind] /= area;
			}
		}
		if(dB){
			double SCALE = 10/Math.log(10);
			for (int ind = 0; ind < n; ind++){
				if(pixels[ind] > 0.000000001)
					pixels[ind] = (float)(SCALE*Math.log(pixels[ind]));
				else
					pixels[ind] = -90;
			}
		}
		fp.setMinAndMax(0,0);
		ImagePlus imp = new ImagePlus(title,fp);
		imp.show();
	}

	//Diffraction limit for the point spread function of a microscope with a circular aperture
	//input a = 2*pi*NA/lambda, r = radius from the point source.
	//See http://yakko.bme.virginia.edu/phy506/
	float psf(float r, float a){
		float b = 2*a*J1OverX(r*a);
		return b*b;
	}

	//Bessel function J1(x)/x.  Uses the polynomial approximations on p. 370 of Abramowitz & Stegun
	//The error in J1 is supposed to be less than or equal to 4 x 10^-8.
	float J1OverX(float xIn){
		double x = xIn;
		if (x < 0) x *= -1;
		double r;
		if (x <= 3){
			double y = x*x/9;
			r = t[0] + y*(t[1] + y*(t[2] + y*(t[3] + y*(t[4] + y*(t[5] + y*t[6])))));
		}else{
			double y = 3/x;
			double theta1 = x + p[0] + y*(p[1] + y*(p[2] + y*(p[3] + y*(p[4] + y*(p[5] + y*p[6])))));
			double f1 = f[0] + y*(f[1] + y*(f[2] + y*(f[3] + y*(f[4] + y*(f[5] + y*f[6])))));
			r = Math.sqrt(1/x)*f1*Math.cos(theta1);
			r /= x;
		}
		return (float)r;
	}
}