[3D Shell Intensities] Segments nucleus and creates concentric shells to measure the average intensity of a channel of interest #Fiji #ImageSC #Macro

3D Shell Measurement.groovy

// Make 3D Shell and measure intensity at given thicknesses
// By Olivier Burri, EPFL - SV - PTECH - BIOP
// As a response to an ImageSC forum post: 
// https://forum.image.sc/t/help-with-shell-analysis-of-nuclei/54342
// Last update: 2021/07/23

raw = getTitle();

close("\\Others");

// Get voxel size to make isotropic Z
getVoxelSize(width, height, depth, unit);
ratio = depth/width;


setBatchMode(true);

// Get nuclei channel
run("Duplicate...", "title=nuclei duplicate channels=4");

// Make isotropic, blur and threshold
run("Scale...", "x=1.0 y=1.0 z="+ratio+" interpolation=Bilinear average process create");
run("Gaussian Blur...", "sigma=2 stack");
run("Median...", "radius=2 stack");
setAutoThreshold("Huang dark stack");
setOption("BlackBackground", true);
run("Convert to Mask", "method=Huang background=Dark black");

// Create distance map (in pixels)
run("Exact Euclidean Distance Transform (3D)");
run("mpl-inferno");

distance = getTitle();

// Get max projection for getting maximum distance
run("Z Project...", "projection=[Max Intensity]");
mProj = getTitle();

getStatistics(area, mean, min, max_edt);

// Prepare to create shells around the nucleus and measure the area of each shell and the intensity of the signal in that area
channel = 1;
// Duplicate channel
selectImage(raw);
run("Duplicate...", "title=measure duplicate channels="+channel);
run("Scale...", "x=1.0 y=1.0 z="+ratio+" interpolation=Bilinear average process");

measure = getTitle();

// Make shells that are 2px thick
thickness = 2
for (i = 1; i <= max_edt/thickness; i++) {
	lower = i*thickness;
	upper = (i+1)*thickness;
	makeShellOfImage(distance, lower, upper, measure); // custom function, see below
	
	
	// Get all statistics from 3D images
	// Project the masked intensities we just got
	run("Z Project...", "projection=[Sum Slices]");
	getStatistics(areas, means);
	r = nResults;
	setResult("Label", r, raw);
	setResult("Channel", r, channel);
	setResult("Shell", r, lower);
	// Get the total intensity in that shell
	setResult("Total intensity", r, (areas*means));

	// get the volume of that shell in pixels
	selectImage("Shell "+lower+"-"+upper+" px");
	run("Z Project...", "projection=[Sum Slices]");
	getStatistics(aream, meanm);
	setResult("Total volume of shell", r, (aream*meanm));
	setResult("Average Signal in Shell", r, ((areas*means) / (aream*meanm)));	
}

// Combine all max projections into a 2 channel image (channel of interest + shell)
run("Images to Stack", "name=Stack title=SUM use");
getDimensions(width, height, channels, slices, frames);
run("Stack to Hyperstack...", "order=xyczt(default) channels=2 slices="+(slices/2)+" frames=1 display=Color");
setBatchMode(false);

// END of script

// Get the intensity of the channel for each shell 
function makeShellOfImage(distImage, thrLow, thrHigh, imageToMask) {
	selectImage(distImage);
	run("Duplicate...", "duplicate");
	setThreshold(thrLow, thrHigh);
	run("Convert to Mask", "background=Dark black");
	rename("Shell "+thrLow+"-"+thrHigh+" px");
	run("Divide...", "value=255 stack");
	imageCalculator("Multiply create stack", imageToMask,"Shell "+thrLow+"-"+thrHigh+" px");
	rename("Shell "+thrLow);
}