rollingBallTestCLIJ.ijm

// Rolling Ball Test via Binary 3D Landscape
// Mark Kittisopikul
// December 2019

// References
// 1. SR Sternberg. Biomedical Image Processing. IEEE Computer 1983.
//    https://doi.ieeecomputersociety.org/10.1109/MC.1983.1654163
// 2. SR Sternberg. Grayscale Morphology. CVGIP 1986
//    https://doi.org/10.1016/0734-189X(86)90004-6
// 3. R Adams. “Radial Decomposition of Discs and Spheres”.
//    CVGIP: Graphical Models and Image Processing,
//    vol. 55, no. 5, September 1993, pp. 325-332.
//    https://doi.org/10.1006/cgip.1993.1024
// 

// Setup

// Parameters
var radius = 10; // default rolling ball radius
var smooth = 3;  // default Gaussian sigma for smoothing
var output8Bit = 1; // default parameter

// Input Image Properties
var imageBitDepth = 0; // bit depth of input image
var inputID = 0; // id of input image
var min = 0;     // min pixel value for building binary
var max = 0;     // max pixel value for building binary
var title = "";  // title of input image

// Technical
var i = 0;       // for loop iterator
var startTime = 0; // time in millisecds for profiling

// Output
var output = "CLIJ Rolling Ball BG Subtraction of "; // output prefix
var bgOutput = "CLIJ Rolling Ball BG of "; // background output prefix


// If no image is open, open blobs
list = getList("image.titles");
if (list.length == 0) {
    run("Blobs (25K)");
}

// Check bit depth and force 8-bit
imageBitDepth = bitDepth();
if (imageBitDepth != 8) {
	run("8-bit");
}

// Acquire image title and setup output
title = getTitle();
inputID = getImageID();
output = output + title;
bgOutput = bgOutput + title;


// Ask for parameters
Dialog.create("Rolling Ball via CLIJ Binary 3D");
Dialog.addNumber("Rolling Ball Radius",radius);
Dialog.addNumber("Gaussian Pre-Smoothing Sigma",smooth);
Dialog.addCheckbox("Output Unsigned 8-bit",output8Bit)
Dialog.show();

// Retrieve parameters
radius = Dialog.getNumber();
smooth = Dialog.getNumber();
output8Bit = Dialog.getCheckbox();

// Start Clock
startTime = getTime();

// Get Min and Max
getMinAndMax(min, max);

// Start CLIJ
run("CLIJ Macro Extensions", "cl_device=");
Ext.CLIJ_clear();
Ext.CLIJ_push(title);

// Smooth if requested
if(smooth > 0) {
	Ext.CLIJ_blur2D(title, "blurred", smooth, smooth);
} else {
	Ext.CLIJ_copy(title, "blurred");
}

print("Blurring Done");
print(getTime()-startTime);

// Create Binary Volume for Rolling Ball
width  = getWidth();
height = getHeight();

stackDepth = 256;


// Replicate blurred in a stack
Ext.CLIJ_create3D("blurredStack", width, height, stackDepth,8);
Ext.CLIJx_resample("blurred", "blurredStack", 1, 1, stackDepth, false);

print('Resampling to make blurredStack Done');
print(getTime()-startTime);

// Create seqStack where each level is equal to z-coordinate
// Consider use of multiplyImageAndCoordinate when it becomes available
// Create a zStack to interpolate 0 to 255
Ext.CLIJ_create3D("zStack", width, height, 2, 32);
// Create a max-valued stack at the top of zStack
Ext.CLIJ_create2D("max", width, height, 8);
Ext.CLIJ_set("max",255);
Ext.CLIJ_copySlice("max", "zStack", 1);
// Interpolate zStack using resampling
// Initialize as 8-bit before resampling
Ext.CLIJ_create3D("bigStack", width, height, stackDepth*2,8);
Ext.CLIJx_resample("zStack", "bigStack", 1, 1, stackDepth, true);
Ext.CLIJ_create3D("cropStack", width, height, stackDepth,8);
Ext.CLIJ_crop3D("bigStack", "seqStack", 0, 0, stackDepth/2, width, height, stackDepth);

print('Seq Stack done');
print(getTime()-startTime);

// Output binary3D
Ext.CLIJx_greaterOrEqual("blurredStack", "seqStack", "binary3D");

print('Binary Volume Done');
print(getTime()-startTime);

//Clean up
Ext.CLIJ_release("blurredStack");
Ext.CLIJ_release("zStack");
Ext.CLIJ_release("max");
Ext.CLIJ_release("bigStack");
Ext.CLIJ_release("cropStack");

//Preview the Binary 3D Volume
//Ext.CLIJ_copy("binary3D","beforeOpening");
//Ext.CLIJ_pullBinary("beforeOpening");

print('Clean Up Done');
print(getTime()-startTime);

// Opening part 1: Erosion
for(var i=0; i < radius/2; i++) {
	Ext.CLIJ_erodeSphere("binary3D", "binary3D2");
	Ext.CLIJ_erodeBox("binary3D2", "binary3D");
}
// If radius is odd do one more sphere erosion
if(radius%2 == 1) {
	Ext.CLIJ_erodeSphere("binary3D", "binary3D2");
	Ext.CLIJ_copy("binary3D2","binary3D");
}


print('Erosion Done');
print(getTime()-startTime);

// Opening part 2: Dilation
for(var i=0; i < radius/2; i++) {
	Ext.CLIJ_dilateSphere("binary3D", "binary3D2");
	Ext.CLIJ_dilateBox("binary3D2", "binary3D");
}
// If radius is odd do one more sphere dilation
if(radius%2 == 1) {
	Ext.CLIJ_dilateSphere("binary3D", "binary3D2");
	Ext.CLIJ_copy("binary3D2","binary3D");
}

print('Dilation Done');
print(getTime()-startTime);

//Examine Binary 3D volume after opening
//Ext.CLIJ_pullBinary("binary3D");

// Counting the remaining heights
Ext.CLIJ_sumZProjection("binary3D", "zSum");

print('Projection Done');
print(getTime()-startTime);

//Make zSum 8-bit for Output
Ext.CLIJ_convertUInt8("zSum",bgOutput);

//Convert original to float for subtraction
Ext.CLIJ_convertFloat(title,"origFloat");
if(output8Bit) {
	Ext.CLIJ_subtractImages("origFloat", "zSum", "floatOutput")
	Ext.CLIJ_convertUInt8("floatOutput",output);
} else {
	Ext.CLIJ_subtractImages("origFloat", "zSum", output);
}


// Get background and background subtracted image
Ext.CLIJ_pull(bgOutput);
Ext.CLIJ_pull(output);

print('All Done');
print(getTime()-startTime);