//Needs 3D ImageJ suite, ImageScience and IJPB added to Fiji
Dialog.create("Set cutoff");

Dialog.addNumber("Threshold: ", 330.0);
Dialog.addNumber("Padding (pixels): ", 0);
Dialog.addNumber("Minimum object size (pixels): ", 1000);
Dialog.addCheckbox("Use uncalibrated pixel value?",false);
Dialog.addCheckbox("Crop? ", true);
Dialog.addDirectory("Batch folder (blank for open images): ", "");
Dialog.show();

threshold = Dialog.getNumber();
buffer = Dialog.getNumber();
minSize = Dialog.getNumber();
calibrated = !Dialog.getCheckbox();
crop = Dialog.getCheckbox();
dir = Dialog.getString();

if(dir == ""){
	list = getList("image.titles");
	setBatchMode(true);
	for (i = 0; i < list.length; i++) {
		selectImage(list[i]);
		processImage(getTitle());		
	}
	setBatchMode(false);
}

else{
	filelist = getFileList(dir);
    outputFolder = dir + File.separator + "Output";
    File.makeDirectory(outputFolder);
	for (i = 0; i < lengthOf(filelist); i++) {
	    if (endsWith(filelist[i], ".tif")) { 
	        open(dir + File.separator + filelist[i]);
	        imgTitle = getTitle();
	        processImage(imgTitle);
	        selectImage(imgTitle);
	        close();
	        list = getList("image.titles");
	        for(j=0; j < lengthOf(list); j++){
	        	selectImage(list[j]);
	        	saveAs("tiff", outputFolder + File.separator + getTitle());
	        }	        
	        close("*");
	    }	    	    
	}
	setBatchMode(false);
}

function processImage(title){
	
	dims = newArray(6);

	Stack.getStatistics(voxelCount, mean, min, tMax, stdDev);
	if(calibrated){			
		setThreshold(threshold, tMax);
		run("Convert to Mask", "method=Default background=Dark black create");
	}
	else{
		setThreshold(threshold, Math.ceil(((tMax-calibrate(0))/(calibrate(1)-calibrate(0)))), "raw");
		run("Convert to Mask", "method=Default background=Dark black create");
	}
	
	for(k=0; k<buffer; ++k){
		run("Dilate (3D)", "iso=255");	
	}	
	
	run("3D Simple Segmentation", "seeds=None low_threshold=128 min_size=" + minSize +" max_size=-1");
	close("Bin");
	selectImage("Seg");
	run("Calibrate...", "function=None unit=mgHA/ccm");
	Stack.getStatistics(voxelCount, mean, min, labelMax, stdDev);
	for (j = 1; j <= labelMax; j++) {
		selectImage("Seg");
		run("Select Label(s)", "label(s)=" + j);
		
		run("Divide...", "value=" + j + " stack");			
		
		imageCalculator("Multiply create stack", title,"Seg-keepLabels");			
		rename(title + " - " + j);
		objTitle = getTitle();;	
						
		if(calibrated){			
			setThreshold(threshold, tMax);
			run("Convert to Mask", "method=Default background=Dark black create");
		}
		else{
			setThreshold(threshold, Math.ceil(((tMax-calibrate(0))/(calibrate(1)-calibrate(0)))), "raw");
			run("Convert to Mask", "method=Default background=Dark black create");
		}
		
		setVoxelSize(1, 1, 1, "pixel");		
		if(crop)
			dims = getBoundingBox();
		
		run("Subtract...", "value=254 stack");
		imageCalculator("Multiply stack", objTitle,"MASK_" + objTitle);
		
		if(!crop) {
			setBatchMode("show");
		}
		if(crop){
			run("TransformJ Crop", "x-range=" + dims[0] + "," + dims[1]  + " y-range=" + dims[2] + "," + dims[3]  + " z-range=" + dims[4] + "," + dims[5]);
			setBatchMode("show");
			close(objTitle);
		}
		
		close("MASK_" + objTitle);
		
		close("Seg-keepLabels");
	}
	close("Seg");
}		


function getBoundingBox(){
	run("Set Measurements...", "bounding stack redirect=None decimal=3");
	run("Analyze Particles...", "display clear stack");
	Array.getStatistics(Table.getColumn("BX"), xMin, max, mean, stdDev);
	Array.getStatistics(Table.getColumn("BY"), yMin, max, mean, stdDev);
	Array.getStatistics(Table.getColumn("Slice"), zMin, zMax, mean, stdDev);
	
	run("Rotate 90 Degrees Right");
	run("Rotate 90 Degrees Right");
	
	run("Analyze Particles...", "display clear stack");
	getDimensions(width, height, channels, slices, frames);
	Array.getStatistics(Table.getColumn("BX"), min, max, mean, stdDev);
	xMax = width-min-1;
	Array.getStatistics(Table.getColumn("BY"), min, max, mean, stdDev);
	yMax = height-min-1;	
	
	run("Rotate 90 Degrees Right");
	run("Rotate 90 Degrees Right");
		
	selectWindow("Results");
	run("Close");

	return newArray(xMin,xMax,yMin,yMax,zMin,zMax);
}