/*Radial_And_Line_profile.txt

N. Vischer
17.08.19, 20:51

Demo to measure FWHM of spots in two different ways:

1. A noisy spot with sigma = sigA and amplidute = ampA
	is used as test image
2. center of noisy spot is detected 
2. Radial profile (from area) is fitted to Gauss
3. Line profile (from line) is fitted to Gauss
4. optionaly repeat with next random seed
*/

var seed = 0, //for reproducible noise
	noiseA = 10,  //=stdev
	sigA = 2,  //sigma of peak
	ampA = 50,  //amplidude of peak
	kRad = 7, //evaluation radius in pixels
	backgnd = 1000,
	expandFactor = 3, //tmp expansion to find sub-pixel max position
	sigM = 8, //apply tmp gaussian smoothing before find maxima
	;

macro "Create Image [F1]"{
	while(true){
		random("seed", ++seed);
		run("Close All");
		createNoisySpot();
		findSubResMaximum(expandFactor);//
		plotRadialProfile();
		lineProfile();
		close("Plot Values");
		waitForUser("Press 'OK' to continue\nPress 'Esc' to stop");
 	}
}
	
	
function createNoisySpot(){
	close("Noisy*");
	ww=80;
	cc = ww/2;
	newImage("NoisySpot", "32-bit black", ww, ww, 1);
	setMetadata("label", "sigmaA=" + d2s(sigA, 1) + " ampA=" + ampA);
	pxSize = 1;
	setVoxelSize(pxSize, pxSize, 1, "um");
	run("Set... ", "zoom=100 x=&cc y=&cc");
	setLocation(4, 100, 400, 400);
	run("Set... ", "zoom=2000 x=&cc y=&cc");
	setPixel(cc, cc, 1000);//will be smoothed
	run("Gaussian Blur...", "sigma=&sigA scaled");
	mul = ampA / getPixel(cc, cc);
	run("Multiply...", "value=&mul");
	run("Add...", "value=&backgnd");
	run("Enhance Contrast", "saturated=0.35");
	run("Add Specified Noise...", "standard=&noiseA");
}

function findSubResMaximum(expandF){
		selectImage("NoisySpot");
		run("Scale...", "x=&expandF y=&expandF interpolation=Bilinear create title=Big");
		run("Gaussian Blur...", "sigma=sigM");
		run("Find Maxima...", "prominence=1 output=[Point Selection]");
		getSelectionCoordinates(xx, yy);
		close;
		selectImage("NoisySpot");
		x = xx[0]/expandF;
		y = yy[0]/expandF;
		makePoint(x, y);
}


//Expects an image with Point selection
//Creates a radial profile with gauss fit

function plotRadialProfile(){
	
 	selectWindow("NoisySpot");
	if (selectionType!= 10)
		exit("Point selection expected");
	getSelectionCoordinates(xx, yy);
	
	//nPoints = xx.length;
	//nPoints = 1;


	Plot.create("Radial Profile", "Radius [um]", "Intensity");
	Plot.setFrameSize(450, 300 );
	Plot.setFontSize(18); 
	
	
	xC = xx[0];
	yC = yy[0];
	close("Peak Profile_*");
	getPixelSize(dummy, pxSize, dummy);
	rr1 = newArray(kRad * kRad * 9);//will be trimmed
	val1 = newArray(kRad * kRad * 9);
	cnt = 0;
	for(y = -kRad; y< kRad; y++){
		for(x = - kRad;  x< kRad; x++){
			r= (sqrt(x*x + y*y));
			if(r <= kRad + 0.1){
				r *= pxSize;
				val = getPixel(xC + x,  yC + y);
				val1[cnt] = val;
				rr1[cnt++] = r;
				val1[cnt] = val;
				rr1[cnt++] = -r;
			}
		}
	}
	val1 = Array.trim(val1, cnt);
	rr1 = Array.trim(rr1, cnt);
	
	Fit.doFit("Gaussian", rr1, val1);
	for(jj = 0; jj < cnt; jj++){//don't paint left circles
		if(rr1[jj] < 0){
			rr1[jj] = NaN;
			val1[jj] = NaN;
		}
	}
	rr1 = Array.deleteValue(rr1, NaN);
	val1 = Array.deleteValue(val1, NaN);
	rSq = Fit.rSquared;
	yoff = Fit.p(0);//backgnd
	ampB = Fit.p(1) - Fit.p(0);//amp
	
	
	
	
	sigB = Fit.p(3);//sigma
	fwhm = sigB * 2.355;	
	intFl1 = 2*PI*ampB*sigB * sigB; 

	
	n = 100;//plot gauss points
	xFit = newArray(-8*pxSize, 8*pxSize);//range in um
	xFit = Array.resample(xFit, n);
	yFit = newArray(n);
	for (i=0; i<n; i++)
		yFit[i] = Fit.f(xFit[i]);

	Plot.setLineWidth(1);
	Plot.setColor("red", "red");
	Plot.add("circle", rr1, val1);
	Plot.setColor("blue");
	Plot.setLineWidth(2);
	Plot.add("line", xFit, yFit);
	Plot.setFontSize(14);
	msg = "";
	msg += "\nSigmaB = " + d2s(sigB, 3);
	msg += "\nAmpl = " + d2s(ampB, 1);
	msg += "\nyOff = " + d2s(yoff,3);

	msg += "\n \nFWHM = " + d2s(fwhm, 3);
	msg += "\nIntFl = " + d2s(intFl1, 2);
	msg += "\nrSqr = " + d2s(rSq,3);
	msg += "\nseed = " + seed;
	
	Plot.addText(msg, 0.05, 0.12);
	Plot.setLimits(NaN, NaN, NaN,NaN );
	Plot.show;
	setSlice(1);
	setLocation(360, 100);
	selectWindow("ImageJ");
	selectWindow("Radial Profile");
	setLocation(360, 100);//work around
}

function lineProfile(){
	boxLen = kRad *2 + 1;
 	selectWindow("NoisySpot");
	if (selectionType!= 10)
		exit("Point selection expected");
	getPixelSize(dummy, pxSize, dummy);
	getSelectionCoordinates(xx, yy);
	x = round(xx[0]);
	y = round(yy[0]);
	rr1 = newArray(-kRad, kRad);
	rr1 = Array.resample(rr1, boxLen);
	makeRectangle(x - kRad, y,  boxLen, 1);
	val1 = getProfile;
	
	
	Fit.doFit("Gaussian", rr1, val1);
	rSq = Fit.rSquared;
	yoff = Fit.p(0);//backgnd
	amp = Fit.p(1) - Fit.p(0);//amp
	sigC = Fit.p(3);//sigma
	fwhm = sigC * 2.355;	
	intFl1 = 2*PI*amp*sigC * sigC; 
	
	n = 100;//plot gauss points
	xFit = newArray(-kRad*pxSize, kRad*pxSize);//range in um
	xFit = Array.resample(xFit, n);
	yFit = newArray(n);
	for (i=0; i<n; i++)
		yFit[i] = Fit.f(xFit[i]);

	Plot.create("Line Profile", "Radius [um]", "Intensity");
	Plot.setFrameSize(450, 300 );
	Plot.setFontSize(18);  
	
	
	Plot.setLineWidth(1);
	Plot.setColor("red", "red");
	Plot.add("circle", rr1, val1);
	Plot.setColor("blue");
	Plot.setLineWidth(2);
	Plot.add("line", xFit, yFit);
	Plot.setFontSize(14);
	pt = 0;
	msg = "";
	msg += "\nSigmaC = " + d2s(sigC, 3);
	msg += "\nAmpl = " + d2s(amp, 1);
	msg += "\nyOff = " + d2s(yoff,3);

	msg += "\n \nFWHM = " + d2s(fwhm, 3);
	msg += "\nIntFl = " + d2s(intFl1, 2);
	msg += "\nrSqr = " + d2s(rSq,3);
	msg += "\nseed = " + seed;

	Plot.addText(msg, 0.05, 0.12);
	Plot.setLimits(NaN, NaN, NaN,NaN );

	FWHM = fwhm;
	RSqr = rSq;
	Ampl = amp;
	YOff = yoff;
	IntFl1 = intFl1;
	Img = getTitle;		

	Plot.show;
	//rename("Noisy Peak Profile");
	setSlice(1);
	setLocation(920,100);
	selectWindow("ImageJ");
	selectWindow("Line Profile");
	setLocation(920,100);//work around

}