added Ilastik-based segmentation files

notebooks/ImageJ/Classify/Run_classification.ijm

@@ -0,0 +1,300 @@
+/*
+ * This script applies a previously trained classifier
+ * for HE images to all HE images that match a predefined
+ * directory structure.
+ * 
+ * The structure is expected to look like this:
+ * root	| Mouse 1	|
+ * 			...
+ * 		| Mouse N	|Histology	| Sample_1	|
+ * 									...
+ * 								| Sample_N	| HE filename
+ * 	
+ * 	The HE image is identified by the string "HE" in the filename
+ * 	The output image is then saved into the newly created directory 
+ * 	"..../Sample x/seg/"
+ * 	
+ * 	Steps that are taken:
+ * 		- Image is openend in max resolution
+ * 		- Image is saved as .hdf5 format for Ilastik to understand
+ * 		- Ilastik is called to segment image
+ * 		- hdf5 image is deleted again to save storage
+ * 		- Output image is moved to output directory
+ */
+
+// Config
+close("*");
+roiManager("reset");
+
+var usebatch =  false;
+var hasbeencropped = false;
+setBatchMode(usebatch);
+
+type = "IF";	// which type of image is segmented? "HE" or "IF"
+reseg = true; 	// if true, segmentation will be done all over again. If false, only unprocessed samples will be segmented
+
+// Ilastik files
+root = "E:/Promotion/Projects/2020_Radiomics/Data/";
+classifier = root + "../Scripts/Histology/Classify/"+type+"_classifier_AD_JM_2.ilp";
+Ilastik_exe = "C:/Program Files/ilastik-1.3.3post3/run-ilastik.bat";
+
+// Input dir
+mice = getsubdirs(root);
+
+// iterate over all mice
+for (i = 0; i < mice.length; i++) {
+
+	print(d2s(i+1, 0) + "/" + d2s(mice.length, 0) + " animals");
+	dir_sample = root + mice[i] + "Histology/";
+	samples = getsubdirs(dir_sample);
+
+	// iterate over samples (aka tissue sections
+	for (j = 0; j < samples.length; j++) {
+		fnames = findFileByID(dir_sample + samples[j], type);
+
+		// Check if all images are provided:
+		// IF HE or IF is missing, there's no need for segmentation
+		files = getFileList(dir_sample + samples[j]);
+		if (files.length <2) {
+			print("    Probably some data missing in " + dir_sample + samples[j]);
+			continue;
+		}
+
+		// check plausibility of filenames
+		if (fnames.length>1) {
+			intrpt = getBoolean("    Something's weird here. Check?");
+			if (intrpt) {
+				exit();
+			} else {
+				continue;
+			}
+		}
+
+		if (fnames.length == 0) {
+			print("    ---> No IF image found here");
+			continue;
+		}
+
+		fname = fnames[0];  // extract fname from single-element array <fnames>
+
+		// Make output directory ( if it doesn't already exist)
+		seg_dir = dir_sample + samples[j] + "1_seg/";
+		if (!File.exists(seg_dir)) {
+			File.makeDirectory(seg_dir);
+		}
+
+		// Check if seg has been done before
+		if (!(reseg || ! File.exists(seg_dir + type + "_seg_Simple_Segmentation.tif"))) {
+			print("    Has been segmented before. Skipping.");
+			continue;
+		}
+
+		// (Pre-) process
+		image = Load(dir_sample + samples[j] + fname);	// Load image
+
+		// Normalize if IF and identify black tiles if HE
+		if (type == "IF") {	image = Normalize(image); }
+		if (type == "HE") {	getTileBackground(image, seg_dir); }
+
+
+		// If a ROI has been provided
+		if ((roiManager("count") > 0) && (!hasbeencropped)) {
+			roiManager("save selected", dir_sample + samples[j] + File.getNameWithoutExtension(fname) + ".roi");
+			roiManager("select", 0);
+			run("Crop");
+		}		
+
+		image = getTitle();
+		f_hdf5 = ExportH5(image, seg_dir + type + "_seg");	// convert to H5 and save
+		RunIlastik(Ilastik_exe, classifier, f_hdf5);  //  Call Ilastik
+
+		// Clean up
+		close("*");
+		roiManager("reset");
+	}
+}
+
+function RunIlastik(IL_executable, classifier, filename){
+	// Runs Ilastik <executable> to apply <classifier> on <filename>
+	print("    Running pixel classification on " + filename);
+	outpath = File.getDirectory(filename);
+	exec(IL_executable, 
+			"--headless",
+			"--project=" + classifier,
+			filename, 
+			"--output_format=tif",
+			"--export_source=Simple Segmentation",
+			"--output_filename_format="+outpath+"/{nickname}_Simple_Segmentation");	
+	File.delete(filename);	// remove h5 image to save storage
+}
+
+function getTileBackground(image, dst_dir){
+	selectWindow(image);
+	setSlice(1);
+	run("Duplicate...", "duplicate title=TileBackground0 channels=1");
+	run("Duplicate...", "duplicate title=TileBackground1 channels=1");
+
+	// identify pixels with value 0
+	selectWindow("TileBackground0");
+	setThreshold(0, 0);
+	run("Convert to Mask");
+
+	// identify pixels with value 255
+	selectWindow("TileBackground1");
+	setThreshold(255, 255);
+	run("Convert to Mask");
+
+	imageCalculator("OR create", "TileBackground0","TileBackground1");
+	result = getTitle();
+	close("TileBackground0");
+	close("TileBackground1");
+
+	// extract ROI
+	selectWindow(result);
+	run("Erode");
+	setThreshold(128, 255);
+	run("Create Selection");
+
+	// smooth ROI outline
+	run("Enlarge...", "enlarge=-30 pixel");
+	run("Enlarge...", "enlarge=60 pixel");  // expand ROI by 20 pixels to include edges in segmented image
+	roiManager("add");
+	close(result);
+
+	roiManager("save", dst_dir + "BlackTiles.zip");
+
+}
+
+function ExportH5(img, f_output){
+	// saves input image <image> as <f_output>.h5
+	selectWindow(img);
+	run("Export HDF5", 	"select="+f_output+".h5 " + 
+						"exportpath=" + f_output + ".h5 " + 
+						"datasetname=data " +
+						"compressionlevel=0 " +
+						"input=["+img+"]");
+	return f_output +".h5";
+}
+
+function Normalize (image){
+	/*
+	 * Function for the preprocessing of IF image data.
+	 * Images are re-sorted into RGB stacks, converted to 16-bit
+	 * and normalized. Brightest pixels are clipped to eliminate
+	 * the influence of staining artifacts.
+	 */
+
+	selectWindow(image);
+	if (nSlices ==1) {
+		run("RGB Stack");
+	}
+	run("Make Composite");
+	run("16-bit");
+
+	// do BG subtraction and normalization
+	// check if sample comes from Slidescanner
+	N = nSlices;
+	if (N == 4){
+		C0 = 2;
+		C1 = 4;
+	} else {
+		C0 = 1;
+		C1 = 3;
+	}
+
+	// in every loop, copy channel, subtract background and normalize
+	index = 1;
+	for (i = C0; i <= C1; i++) {
+		selectWindow(image);
+		setSlice(i);
+		run("Duplicate...", "title=" + index);
+		selectWindow(index);
+		run("Subtract Background...", "rolling=50 stack");
+		run("Enhance Contrast...", "saturated=0.3 normalize");
+
+		index += 1;
+	}
+	exit();
+
+	// if it's slidescanner data (which is typically badly cropped)
+	// give user chance to draw his own ROI
+	if (N == 4) {
+		setBatchMode(false);
+		selectWindow("1");
+		selectWindow("2");
+		selectWindow("3");
+		waitForUser("Set a ROI if image is badly cropped.\n"+
+					"Or do nothing if it's fine");
+
+		if (selectionType() != -1) {
+			roiManager("add");
+			roiManager("select", 0);
+		}
+	}
+	setBatchMode(usebatch);
+
+	close(image);
+	run("Merge Channels...", "c1=1 c2=2 c3=3 create");
+	rename(image);
+	return image;
+}
+
+function Load(f_input){
+	// converts an input file <f_input> into an hdf5 file <f_ouput>.
+	// The ".h5" file extension will be added on the fly.
+	//print("    --->"+f_input);
+
+	hasbeencropped = false;
+
+	run("Bio-Formats Importer", "open="+f_input+" autoscale "+
+				"color_mode=Default "+
+				"rois_import=[ROI manager] "+
+				"view=Hyperstack "+
+				"stack_order=XYCZT "+
+				"series_1");
+
+	dir = File.getDirectory(f_input);
+	f  = File.getNameWithoutExtension(f_input);
+
+	if (File.exists(dir + f + ".roi")) {
+		roiManager("open", dir + f + ".roi");
+		roiManager("select", 0);
+		run("Crop");
+		hasbeencropped = true;
+	}
+
+	return getTitle();
+}
+
+function findFileByID(path, ID){
+	// browse directory <path> and check for files that match the given ID string <ID>
+	files = newArray();
+	filelist = getFileList(path);
+	for (i = 0; i < filelist.length; i++) {
+
+		// exclude ROIs from search
+		if (endsWith(filelist[i], "roi")) {
+			continue;
+		}
+
+		if (matches(filelist[i], ".*" + ID + ".*")) {
+			files = Array.concat(files, filelist[i]);
+		}
+	}
+
+	print("    Found " + d2s(files.length, 0) +" files matching " + ID + " in " + path);
+	Array.show(files);
+	return files;
+}
+
+function getsubdirs(path){
+ 	// Returns array with all subdirectories in <path>
+ 	list = getFileList(path);
+ 	dirlist = newArray();
+ 	for (i = 0; i < list.length; i++) {
+ 		if (File.isDirectory(path + list[i])) {
+ 			dirlist = Array.concat(dirlist, list[i]);
+ 		}
+ 	}
+ 	return dirlist;
+ }