setBatchMode(false);

// Select the slice with highest intensity and Set display;
Maxstack=0; Minstack=65500; OptSlice=0;
for (i=1; i<nSlices+1; i++) {
	setSlice(i);
	getStatistics(area, mean, min, max, std, histogram);
	if (max>Maxstack) {
		Maxstack=max;
		Minstack=min;
		OptSlice=i;
		}
}
setSlice(OptSlice);
setMinAndMax(Minstack, Maxstack);

// Get info on the setup;

getDateAndTime(year, month, dayOfWeek, dayOfMonth, hour, minute, second, msec);
Dayadjust = "";
if (dayOfMonth<10) {Dayadjust = "0";}
month=month+1;
Monthadjust = "";
if (month<10) {Monthadjust = "0";}
date = ""+year+"-"+Monthadjust+month+"-"+Dayadjust+dayOfMonth;
Dialog.create("Image and microscope information");
  Dialog.addString("Microscope Name:", "");
  Dialog.addNumber("Magnification:", 100);
  Dialog.addString("NA:", "1.4");
  Dialog.addNumber("Pixel size (nm):", 60);
  Dialog.addString("Date:", date);
Dialog.show();
  microscope = Dialog.getString();
  MA = Dialog.getNumber();
  NA = Dialog.getString();
  xyVoxel = Dialog.getNumber();
  date = Dialog.getString();

MainName=date+"_"+microscope+"_"+MA+"x_"+NA;
zVoxel=200;
setVoxelSize(xyVoxel, xyVoxel, zVoxel, "nm");

rename("Stack");
ROIsize=round(15000/xyVoxel);
halfROIsize = round(ROIsize/2);
      boucle=1;
      while (boucle!=0) {
          getCursorLoc(x, y, z, flags);
          if (flags==4) boucle=0;
          }

OptSlice=z;
zProfileX = newArray(100); zProfileY = newArray(100);
for (i=0; i<nSlices; i++) {
	setSlice(i+1);
	zProfileX[i] = i;
	zProfileY[i] = getPixel(x, y);
	}

xROI=x-halfROIsize; yROI=y-halfROIsize;
makeRectangle(xROI, yROI, ROIsize, ROIsize);
run("Crop"); run("32-bit");
run("Duplicate...", "title=[Stack FWHM] duplicate");



// _______________________Projections_______________________ ;

selectWindow("Stack");
makeLine(0, halfROIsize, ROIsize, halfROIsize);
run("Reslice [/]...", "input="+zVoxel+" output="+zVoxel+" slice=1");
rename ("xProj");
H=getHeight();
selectWindow("Stack");
makeLine(halfROIsize, 0, halfROIsize, ROIsize);
run("Reslice [/]...", "input="+zVoxel+" output="+zVoxel+" slice=1 rotate");
rename ("yProj");
selectWindow("Stack");
SliceNum=getSliceNumber();
run("Z Project...", "start=1 stop="+SliceNum+" projection=[Max Intensity]");
rename("Project");
ProjectWidth=ROIsize+H;
run("Canvas Size...", "width="+ProjectWidth+" height="+ProjectWidth+" position=Top-Left zero");
selectWindow("yProj");
run("Canvas Size...", "width="+ProjectWidth+" height="+ProjectWidth+" position=Top-Right zero");
selectWindow("xProj");
run("Canvas Size...", "width="+ProjectWidth+" height="+ProjectWidth+" position=Bottom-Left zero");
run("Image Calculator...", "image1=Project operation=Add image2=xProj");
selectWindow("Project");
run("Image Calculator...", "image1=Project operation=Add image2=yProj");
run("Size...", "width=550 height=550 constrain interpolate");
selectWindow("Stack");
close();
selectWindow("xProj");
close();
selectWindow("yProj");
close();
selectWindow("Project");
run("Square Root");
ROIsizeBG = round(ROIsize/10);
makeRectangle(ROIsizeBG, ROIsizeBG, ROIsizeBG, ROIsizeBG);
getStatistics(area, mean, min, max, std, histogram);
minSelection=min;
run("Select None");
getStatistics(area, mean, min, max, std, histogram);
setMinAndMax(minSelection, max);
run("Fire");
run("8-bit");
run("RGB Color");



// _______________________FWHM axial_______________________ ;

selectWindow("Stack FWHM");

Fit.doFit("Gaussian", zProfileX, zProfileY); a=Fit.p(0); b=Fit.p(1); c=Fit.p(2); d=Fit.p(3);

// Plot ;
Amplitude = 40;
XplotLatReal = newArray(Amplitude); YplotLatReal = newArray(Amplitude);
MaxGraph = 0;
for (i=0; i<Amplitude; i++) {
	XplotLatReal[i] = (i-Amplitude/2)*zVoxel;
	YplotLatReal[i] = zProfileY[OptSlice-Amplitude/2+i];
	if (YplotLatReal[i]>=MaxGraph) {MaxGraph = YplotLatReal[i];}
	}
XplotLatFit = newArray(Amplitude*4);
YplotLatFit = newArray(Amplitude*4);
Ymin=66000;
Ymax=0;
for (i=0; i<Amplitude*4; i++) {
	XplotLatFit[i] = (i/4-Amplitude/2)*zVoxel;
	X = OptSlice-Amplitude/2+i/4;
	YplotLatFit[i] =  a + (b-a)*exp(-(X-c)*(X-c)/(2*d*d));
	if (YplotLatFit[i]>=MaxGraph) {MaxGraph = YplotLatFit[i];}
	if (Ymin>YplotLatFit[i]) {Ymin=YplotLatFit[i];}
	if (Ymax<YplotLatFit[i]) {Ymax=YplotLatFit[i];}
	}

HM=(Ymax-Ymin)/2; k=-2*d*d*log((HM-a)/(b-a)); FWHMa=2*zVoxel*sqrt(k);

Plot.create("FWHM axial", "Z", "Intensity", XplotLatFit, YplotLatFit);
Plot.setLimits(-4000, 4000, 0, MaxGraph*1.1);
Plot.add("circles", XplotLatReal, YplotLatReal);
Text="FWHM axial ="+d2s(FWHMa,0)+"nm";
Plot.addText(Text, 0, 0);
Plot.show();
run("Canvas Size...", "width=528 height=510 position=Bottom-Center zero");



// _______________________FWHM lateral_______________________ ;

x = newArray(-8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8);
selectWindow("Stack FWHM");
setSlice(OptSlice);
makeLine(halfROIsize-8, halfROIsize, halfROIsize+8, halfROIsize);
y = getProfile();

Fit.doFit("Gaussian", x, y); a=Fit.p(0); b=Fit.p(1); c=Fit.p(2); d=Fit.p(3);

// Plot ;

XplotLatReal = newArray(17); YplotLatReal = newArray(17);
Ymin=66000;
Ymax=0;
MaxGraph = 0;
for (i=0; i<17; i++) {
	XplotLatReal[i] = (i-8)*xyVoxel;
	YplotLatReal[i] = y[i];
	if (y[i]>=MaxGraph) {MaxGraph = y[i];}
	}
XplotLatFit = newArray(65); YplotLatFit = newArray(65);
for (i=0; i<65; i++) {
	XplotLatFit[i] = (i/4-8)*xyVoxel;
	X = i/4-8;
	YplotLatFit[i] =  a + (b-a)*exp(-(X-c)*(X-c)/(2*d*d));
	if (YplotLatFit[i]>=MaxGraph) {MaxGraph = YplotLatFit[i];}
	if (Ymin>YplotLatFit[i]) {Ymin=YplotLatFit[i];}
	if (Ymax<YplotLatFit[i]) {Ymax=YplotLatFit[i];}
	}

HM=(Ymax-Ymin)/2; k=-2*d*d*log((HM-a)/(b-a)); FWHMl=2*xyVoxel*sqrt(k);

Plot.create("FWHM lateral", "X", "Intensity", XplotLatFit, YplotLatFit);
Plot.setLimits(-8*xyVoxel, 8*xyVoxel, 0, MaxGraph*1.1);
Plot.add("circles", XplotLatReal, YplotLatReal);
//Text="FWHM lateral ="+d2s(FWHMl,0)+"nm";
Plot.addText("FWHM lateral ="+d2s(FWHMl,0)+"nm", 0, 0);
Plot.show();
run("Canvas Size...", "width=528 height=510 position=Top-Center zero");
run("Image Calculator...", "image1=[FWHM lateral] operation=Add image2=[FWHM axial]");
run("Canvas Size...", "width=550 height=550 position=Center");
rename("FWHM");
run("RGB Color");
selectWindow("FWHM axial");
close();
selectWindow("Stack FWHM");
close();

MainName=MainName+" FWHMa "+d2s(FWHMa,0)+"nm - FWHMl "+d2s(FWHMl,0)+"nm";
run("Images to Stack");
rename(MainName);