initial commit; moved from imcf-toolbox

PSF_analyser_.ijm

+// PSF_analyser
+// wirtten by Kai Schleicher, IMCF Basel.
+//-------------------------------------------------------------------
+
+run("Bio-Formats Macro Extensions");
+//This will enable the bioformats macro functions:
+Ext.getVersionNumber(version);
+//Returns the version number of the currently installed version of Bio-Formats.
+
+
+print("Version of Bio-Formats required: 5.1.7-DEV or higher");
+print("Version of Bio-Formats installed: " + version);
+print(" ");
+print("If your version number is lower, please update via Help -> update FIJI and restart FIJI");
+print(" ");
+
+// --------------------set user environment -----------------
+// results are only correct no previouse image is open
+
+if(nImages>0){
+	exit("Please close all images before running PSF_analyser")
+}
+
+
+// batch open files ------------------------------------------------------------
+
+dir1 = getDirectory("Choose Source Directory");
+dir2 = getDirectory("Choose Destination Directory");
+// get an array that contains all image names in dir1
+list = getFileList(dir1);
+
+// open all files and work on them -----------------------------------------------------------------------
+
+//setBatchMode(true);
+for (i=0; i<list.length; i++) {
+    showProgress(i+1, list.length);
+    incoming= dir1+list[i];
+   
+    
+   	run("Bio-Formats Importer", "open=[" + incoming + "] color_mode=Default open_all_series split_channels view=Hyperstack stack_order=XYCZT");
+   	// open the image at position i as a hyperstack using the bio-formats
+	// opens all images of a container file (e.g. *.lif, *.sld)
+	// split channels
+	// get image IDs of all open images:
+
+	all=newArray(nImages);
+	//Array.print(all);
+	
+	for (k=0;k<nImages;k++) {
+        selectImage(k+1);
+        all[k]=getImageID;
+        title = getTitle;
+		// now do things ----------------------------------------------------------------------------------------------------------
+	
+		// select senter of the z-stack to better find the beads
+		Stack.getDimensions(width, height, channels, slices, frames);
+		Stack.setSlice(slices/2);  
+	
+		// adjust contrast to better see the beads
+		run("Brightness/Contrast...");
+		run("Enhance Contrast", "saturated=0.35");
+	
+		// userinput: draw an rectangular ROI around a bead
+		setTool("rectangle");
+		waitForUser("please select a rectangular ROI around a bead", " 1. Draw ROI \n 2. press [OK] when done");
+		roiManager("Add");
+	
+		// check if the user did his job properly
+		if(selectionType()!=0){
+			exit("please choose a rectangular ROI");
+		}
+		
+		// crop from the chose ROI
+		run("Duplicate...", "duplicate");
+		// reset B&C to original values, othervise the measurment is screwed up
+		run("Brightness/Contrast...");
+		resetMinAndMax();
+		XY = getImageID();
+		xytitle=getTitle();
+		print("working on image"); 
+		print(xytitle);
+		rename("XY_"+XY);
+		
+		// run orthogonal views on which to perform the measurments
+	
+		run("Orthogonal Views");
+		wait(3000); //Delays (sleeps) for n milliseconds, fiji vodoo
+		
+		//get a list (array) of all open windows, needs a delay after orthogonal views to work
+		imagetitles = getList("image.titles");
+		
+		// go through that list and rename the orthogonal images in a way they can be recalled later
+	
+		for (r=0; r<imagetitles.length; r++){
+			selectImage(imagetitles[r]); 
+			img=getTitle(); 
+			
+			if (substring(img,0,2)=="YZ") { 
+			        YZ = getImageID(); 
+			        rename("YZ_"+YZ); 
+			        yztitle=getTitle(); 
+			} 
+			if (substring(img,0,2)=="XZ") { 
+			        XZ = getImageID(); 
+			        rename("XZ_"+XZ); 
+			        xztitle=getTitle(); 
+			} 
+		} 
+		
+		// define line selection for YZ image
+	 	selectImage(YZ);
+	
+		getPixelSize(unit, pixelWidth, pixelHeight);
+		getDimensions(width, height, channels, slices, frames);
+		getRawStatistics(nPixels, mean, min, max, std, histogram);
+		dynrange=pow(2, bitDepth());
+	
+		//exit the macro if the ROI containes overexposed pixels
+		if(max==dynrange){
+			exit("your selection containes overexposed pixels");
+		}
+	
+		// define global variables
+		var maxX, maxY;
+	
+		// run the "findmaxima" function defined below to find the center of the PSF
+		// is a horizontal line through maximum as long as the image is wide
+		findmaxima();
+		selectImage(YZ);
+		makeLine(0, maxY, width, maxY);
+		roiManager("Add");
+		getSelectionBounds(a, b, c, d);
+		
+	
+		//name line accordingly
+		line="YZ-horizontal";
+		
+		// fit the YZ image with a gaussian fit along the horizontal line selection (results in FWHM_z)
+		gaussfit(a, b, c, d);
+	
+		// repeat for the other direction to get FWHM_y	
+		selectImage(YZ);
+		//vectical line through maximum as long as the image is high
+		makeLine(maxX, 0, maxX, height);
+		roiManager("Add");
+		getSelectionBounds(a, b, c, d);
+		//name line accordingly
+		line="YZ-vertical";
+	
+		// fit the YZ image with a gaussian fit along the line vertical selection (results in FWHM_y)
+		gaussfit(a, b, c, d);
+	
+		// save the psf images and plots as png
+		selectImage(YZ);
+		run("Remove Overlay");
+		//run("Spectrum");
+		//run("LUTforPSFs");
+		run("physics");
+		run("Enhance Contrast", "saturated=0.35");
+		saveAs("PNG", dir2+xytitle+"PSF_YZ.png");
+		rename("PSF");
+		
+		selectWindow("YZ-horizontal");
+		Plot.makeHighResolution("YZ-horizontal_HiRes",4.0);
+		saveAs("PNG", dir2+xytitle+"YZ-horizontal_HiRes.png");
+		rename("YZ-horizontal_HiRes");
+		
+		selectWindow("YZ-vertical");
+		Plot.makeHighResolution("YZ-vertical_HiRes",4.0);
+		saveAs("PNG", dir2+xytitle+"YZ-vertical_HiRes.png");
+		rename("YZ-vertical_HiRes");
+
+		wait(3000); //Delays (sleeps) for n milliseconds, fiji vodoo
+
+		// check if windows exist and if so close them
+		// needed because either Plot or orthogonal views is glitchy and does not always behave the same way
+		
+		if (isOpen("YZ-horizontal")==1) {		
+			selectWindow("YZ-horizontal");
+			close(); 
+		}	
+		
+		if (isOpen("YZ-vertical")==1) {		
+			selectWindow("YZ-vertical");
+			close();
+		}
+		
+		if (isOpen("Plot Values")==1) {		
+			selectWindow("Plot Values"); 
+			run("Close");
+		}
+
+		if (isOpen("YZ-horizontal_HiRes")==1) {		
+			selectWindow("YZ-horizontal_HiRes"); 
+			run("Close");
+		}
+		
+		if (isOpen("YZ-vertical_HiRes")==1) {		
+			selectWindow("YZ-vertical_HiRes"); 
+			run("Close");
+		}
+
+		if (isOpen("XY_"+XY)==1) {		
+			selectWindow("XY_"+XY);
+			close();
+		}
+
+		if (isOpen(yztitle + "-1")==1) {		
+		selectWindow(yztitle + "-1");
+		close();
+		}
+				
+		if (isOpen("YZ_"+YZ+"-1")==1) {		
+			selectWindow("YZ_"+YZ+"-1");
+			close();
+		}     
+	
+
+	} 
+
+	//close all images to free the memory
+	run("Close All");
+	
+	
+}
+
+// save the FWHM table to the destination directory
+selectWindow("FWHM");
+saveAs("Text", dir2+ "FWHMs.txt");
+
+print(" ");
+print("All done");
+
+// end of the batch ----------------------------------------------------------
+
+// Define functions -----------------------------------------------------------------
+
+/////////////////// find maxima /////////////////////////////////////////////////////////////////////////
+
+function findmaxima() {
+	run("Duplicate...", " ");
+	run("Gaussian Blur...", "sigma=5");
+	run("Find Maxima...", "noise=100000 output=List");
+	//returns a result table listing all maxima
+	
+	
+	cycle=0;
+	while (nResults() > 1) {
+		run("Gaussian Blur...", "sigma=5");
+		//sigma = bead diameter in px
+		cycle=cycle+1;
+		print("gaussian blurr cycle: " + cycle);
+	
+		run("Find Maxima...", "noise=100000 output=List");
+		
+		if (cycle == 5){
+			exit ("Cannot find sole maximum. Please change crop size");
+		}
+	}
+	// if several maxima are found, repeat run("Gaussian Blur...", "sigma=2"); until nResults=1,
+	// repeat no more than 5 times
+
+	//define coordinates of maximum from "find maxima" results table= middle of PSF
+	maxX = getResult("X", 0);
+	maxY = getResult("Y", 0);
+	
+	//close the results window
+	selectWindow("Results");
+	run("Close");
+	
+		
+}
+
+/////////////////// make plot axis /////////////////////////////////////////////////////////////////////
+
+function gaussfit(a, b, c, d) {
+	
+	if(c==1){
+		scale=pixelHeight;
+	} else {
+		scale=pixelWidth;
+	}
+	
+	if(channels>1){
+		Stack.getPosition(channel, dummy, dummy);
+		line="Ch"+channel+" "+a+" "+b+" "+c+" "+d;
+	}
+	
+	Y=getProfile();
+	len=Y.length;
+	
+	X=newArray(len);
+	for(i=0;i<len;i++){
+		X[i]=i*scale;
+	}
+	
+	//////////////// Gauss fit ///////////////////////////////////////////////////////////////////////////
+	
+	
+	Fit.doFit("Gaussian", X, Y);
+	r2=Fit.rSquared;
+	
+	//if(r2<0.9){
+	//	showMessage("Warning: Poor Fitting",r2);
+	//}
+	Fit.plot();
+	
+	if(isOpen("y = a + (b-a)*exp(-(x-c)*(x-c)/(2*d*d))")){
+		selectWindow("y = a + (b-a)*exp(-(x-c)*(x-c)/(2*d*d))");
+		rename(line);
+	}
+
+	Mean=Fit.p(2);
+	print(Mean);
+	sigma=Fit.p(3);
+	FWHM=abs(2*sqrt(2*log(2))*sigma);
+	// write results into a table using the myTable function as defined below
+	myTable(line,FWHM,unit,r2);
+	// selectImage(XY);
+	// selectWindow("FWHM");
+}
+  
+/////////////////// make table with fit results /////////////////////////////////////////////////////////////
+
+function myTable(a,b,c,d){
+	title1="FWHM";
+	title2="["+title1+"]";
+	if (isOpen(title1)){
+   		print(title2, xytitle+"\t"+a+"\t"+b+"\t"+c+"\t"+d+"\t"+maxX+"\t"+maxY);
+	}
+	else{
+   		run("Table...", "name="+title2+" width=600 height=200");
+   		print(title2, "\\Headings:Image\tLine\tFWHM\tUnit\tR-square\tCenter_z\tCenter_y");
+   		print(title2, xytitle+"\t"+a+"\t"+b+"\t"+c+"\t"+d+"\t"+maxX+"\t"+maxY);
+	}
+}
+
+// End of functions ----------------------------------------------------------------------------------------------------