/* Wipe_Background.ijm
 * IJ BAR: https://github.com/tferr/Scripts#scripts
 *
 * ImageJ macro that clears (sets to zero) clusters of thresholded pixels of defined
 * circularity & size within the active area ROI (or the whole image if no ROI exists).
 *
 * It tries to be as unobtrusive as possible and accepts any grayscale image including
 * multidimensional hyperstacks. Currently, filtering is restricted to the Z-dimension,
 * but it could be easily extended to other dimensions. It could be greatly simplified
 * if the goal was to deal exclusively with 2D/3D images
 */

// Set requirements
requires("1.49t");
getThreshold(lower, upper);
if (lower==-1 && upper==-1)
	exit("A thresholded image is required.");

// Gather information about image
imgID = getImageID();
isHyper = Stack.isHyperstack;
Stack.getDimensions(null, null, channels, depth, frames);
Stack.getPosition(channel, activeSlice, frame);
roi = selectionType();

// Define options
optnsA = newArray("Active ROI", "Whole image");
optnsB = newArray("Current slice only", "All slices", "Preceding slices", "Subsequent slices");
scope = optnsB[0];
wholeImg = false;
start = activeSlice;
end = activeSlice;
areaRoi = roi > -1 && roi < 5 || roi == 9;

// Prompt for range(s) and scope of filtering
Dialog.create("Wipe Background - Filtering Particles");
	Dialog.addString("Size (pixels):", "1-10000", 12);
	Dialog.addString("Circularity:", "0.00-1.00", 12);
	if (areaRoi)
		Dialog.addChoice("Scope:", optnsA);
	if (depth>1)
		Dialog.addChoice("Apply to:", optnsB);
	Dialog.addMessage("Clusters of thresholded pixels\n"
			+ "fulfilling both conditions will\n"
			+ "be cleared (set to zero).");
	Dialog.addHelp(helpMsg());
Dialog.show();

size = Dialog.getString();
circ = Dialog.getString();
if (areaRoi)
	if (Dialog.getChoice==optnsA[1])
		whole = true;
if (depth>1) {
	scope = Dialog.getChoice();
	if (scope==optnsB[1]) {
		start = 1; end = depth;
	} else if (scope==optnsB[2]) {
		start = 1; end = activeSlice;
	} else if (scope==optnsB[3]) {
		start = activeSlice; end = depth;
	} else {
		start = activeSlice; end = activeSlice;
	}
}

// Do not display any intermediate images
setBatchMode(true);

// Keep track of detected particles
nCount = nResults();

// Retrieve mask
analyzerArg = "size="+ size +" pixel circularity="+ circ +" show=Masks ";
if (scope!=optnsB[0])
	analyzerArg += " stack";
if (wholeImg)
	run("Select None");
run("Analyze Particles...", analyzerArg);

// Active image should now be a mask generated by ParticleAnalyzer
maskID = getImageID();
nCount = nResults() - nCount - 1;

if (nCount>0 && maskID!=imgID) {

	// The ParticleAnalyzer's mask is never a multidimensional image.
	// We'll ensure it has the same dimensions of the original image
	run("Invert", "stack");
	if (isHyper && scope!=optnsB[0]) {
		run("Stack to Hyperstack...",
			"channels="+ channels +" slices="+ depth +" frames="+ frames);
		maskID = getImageID();
	}

	// The ParticleAnalyzer's mask has background = 255, foreground = 0
	// (IJ 1.49u). We'll set foreground particles= 0, remaining pixels to 1
	for (c=1; c<=channels; c++) {
		for (t=1; t<=frames; t++) {
			for (z=1; z<=depth; z++) {
				Stack.setPosition(c, z, t);
				if (z>=start && z<=end)
					run("Divide...", "value=255 slice");
				else
					run("Set...", "value=1 slice");
			}
		}
	}

	// Perform the filtering
	if (scope==optnsB[0])
		imageCalculator("Multiply stack", imgID, maskID);
	else
		imageCalculator("Multiply", imgID, maskID);

	// Ensure mask is not displayed (for whatever reason hyperstacks
	// seem to escape batchMode (at least with 1J 1.49t, Java 1.6)
	selectImage(maskID);
	close();

	// Report progress
	showStatus("Filtered out "+ nCount +" particles...");
	if (isOpen("Log"))
		print(">>> "+ getTitle() +": Filtered out "+ nCount +" particles...");

} else {

	showMessage("No particles were detected with\nthe specified settings:\n \n"
			+ "\n    Size: "+ size
			+ "\n    Circ: "+ circ
			+ "\n    Threshold: "+ lower +"-"+ upper);

}

// Restore original state
Stack.setPosition(channel, activeSlice, frame);
if (areaRoi && (wholeImg || isHyper))
	run("Restore Selection");
setThreshold(lower, upper);
updateDisplay();


function helpMsg() {
	msg = '<html>'
	+ '<div WIDTH=420>'
	+ 'Since bitmap objects are discretised into a regular lattice of pixels, '
	+ 'circularity is highly affected by the particle size and may not be '
	+ 'valid for very small sizes. Here are some calculations using IJ 1.49t: '
	+ '<table align="center">'
	+ '  <tr>'
	+ '    <th>Particle size (pixels)</th>'
	+ '    <th>Circularity range</th>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>1</td>'
	+ '    <td>Always 1</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>2</td>'
	+ '    <td>0.79&mdash;1</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>3</td>'
	+ '    <td>0.52&mdash;1</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>4</td>'
	+ '    <td>0.39&mdash;1</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>&hellip;</td>'
	+ '    <td>&hellip;</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>10</td>'
	+ '    <td>0.31&mdash;1</td>'
	+ '  </tr>'
	+ '  <tr align="center">'
	+ '    <td>64</td>'
	+ '    <td>0.02&mdash;1</td>'
	+ '  </tr>'
	+ '</table>'
	+ 'This means that e.g., if you set <i>Size</i> to <tt>&le;4</tt> and '
	+ '<i>Circularity</i> to <tt>0.0-0.3</tt> no filtering operation will '
	+ 'be performed because it is not possible to fullfill both conditions. '
	+ 'You can read more about shape descriptors on the '
	+ '<a href="http://imagej.nih.gov/ij/docs/guide/">IJ User guide</a> and '
	+ '<a href="http://fiji.sc/Category:Particle_analysis">Fiji website</a>.'
	+ '</div>'
	+ '</html>';
	return msg;
}