[Concentric Circles Around Donut] Macro to create concentric circles based on relative distance to an outer edge for a "donut-like" shape #fiji #imagej #imagesc

concentric-circles.ijm

/**
 * Macro to create concentric circles and sectors based on relative distance to an outer edge for a "donut-like" shape
 * 
 * Author: Olivier Burri, EPFL - SV - PTECH - BIOP
 * 
 * For Kailie Batsche, from a request on the Image.sc forum:
 * https://forum.image.sc/t/cross-sectional-analysis-cut-into-8-equal-sectors-with-three-concentric-rings/50302?u=oburri
 * 
 * Provided under GPL3
 * https://www.gnu.org/licenses/gpl-3.0.en.html
 * 
 * Please cite this Gist if it was useful to you!
 */

#@ImagePlus rawImage
#@Integer nRings (label="Number of Rings", value=3)
#@String ringNames =  (label="Ring Names (From inner to outer)", value="Inner, Mid, Outer")
#@Integer nSectors (label="Number of Sectors", value=8)
#@Integer downsample (label="Downsampling for Detection", value=4)
setBatchMode(true);

close("\\Others");
roiManager("reset");

run("Select None");

run("Scale...", "x="+(1/downsample)+" y="+(1/downsample)+" create");

// Remove calibration
setVoxelSize(1.0, 1.0, 1.0, "pixel");

//run("Duplicate...", " ");
run("Gaussian Blur...", "sigma=2");

//setThreshold(0, 211);
setOption("BlackBackground", true);
run("Median...", "radius=15");

setAutoThreshold("Huang");
run("Create Mask");
//run("Threshold...");
run("Set Measurements...", "centroid center redirect=None decimal=3");
run("Analyze Particles...", "size=100000-Infinity clear display pixel include add");

// Get Centroid to get center. Useful for determining inner area and for segment origin
xm = getResult("XM", nResults-1);
ym = getResult("YM", nResults-1);
selectImage("mask");
doWand(xm, ym);
roiManager("Add");


roiManager("select", 0);
// Need to do some math to get the normalzied distances.

// Set Distance Map to 32-bit
run("Options...", "iterations=1 count=1 black edm=32-bit do=Nothing");

// Set colors to draw necessary masks
setForegroundColor(255, 255, 255);
setBackgroundColor(0,0,0);

// Get the dimensions of the downsampled image to get the max radius available and to be able to create the mask images.
getDimensions(width, height, channels, slices, frames);

// Compute distance to outer. That is, how far inside we are from the outermost perimeter
newImage("Outer", "8-bit black", width, height, 1);

// Select the outer region
roiManager("select", 0);

// Fill it up
run("Fill", "slice");

// Get the Distance Map
run("Distance Map");

// Set the outside and inside to 0 so it ends up as Not a Number (NaN) for subsequent calculations
roiManager("select", 0);
run("Clear Outside");
roiManager("select", 1);
run("Clear");

// Get the distances to the inner perimeter. This is the inverse image of the inner region
newImage("Inner", "8-bit black", width, height, 1);
roiManager("select", 1);
run("Fill", "slice");
run("Select None");

// After the next step, everything is white, except for the inner region, which is black
run("Invert");
run("Distance Map");

// Again clean up the region by setting it to 0 outside of the outer region. The inner region is already at 0
roiManager("select", 0);
run("Clear Outside");

// Sum the two distances
imageCalculator("Add create 32-bit", "EDM of Inner","EDM of Outer");
rename("Sum");

// Divide the inner distance by the sum. This gives us a normalized map where 0 is when it touches the inner area and 1.0 is when it touches the outer area.
imageCalculator("Divide create 32-bit", "EDM of Inner","Sum");

roiManager("reset");

// Draw a line pointing outwards from the center in the image to create the sectors
setColor(0);

for( i=0; i< nSectors; i++ ) {
	xt = xm + getWidth() * cos(2*PI*i/ nSectors);
	yt = ym + getWidth() * sin(2*PI*i/ nSectors);
	setLineWidth(2);
	drawLine(xm, ym, xt, yt );
}

// Now cut everything up by ring and sector and save it with the designated names in the ROI Manager
allNames = split(ringNames, ",");

// This loop creates the rings, based on the relative distance map we created.
for(r=0; r< nRings; r++) {
	// Set a threshold based on the number of rings
	setThreshold(r/nRings+0.0001, ( (r+1) / nRings ));

	// This loop handles creating the N sectors around each ring.
	for( s=0; s< nSectors; s++ ) {

		// Get the name we want for all the sectors, and add a label to them (1, 2, ..., N)
		name = String.trim(allNames[r]+"_S_"+(s+1));

		// This custom method will create a selection at the right place. See below
		select(r,s);
		
        // Set the name
		Roi.setName(name);

		// Add it to the manager
		roiManager("add");
	}
}

// We are almost done, we need to set the regions back to the original size
nR = roiManager("count");
for(r = 0; r<nR; r++) {
	roiManager("Select", r);
	run("Scale... ", "x="+downsample+" y="+downsample);
	roiManager("Update");
}

// Place them onto the original image
selectImage(rawImage);

roiManager("Show All with labels");
setBatchMode(false);

// end of script

// Helper methods below

// This method computes the coordinate of a point that is sure to be inside and in the middle of the sectors
function select( ring, sector ) {
		// Make a line that will span the whole length
		xt = xm + getWidth() * cos(2*PI*sector / nSectors + 2*PI / nSectors / 2 );
		yt = ym + getWidth() * sin(2*PI*sector / nSectors + 2*PI / nSectors / 2 );
		
		makeLine(xm, ym, xt, yt );
		
		// Get the intensity (Values of the normalized distance) that this line intersects
		profile = getProfile();

		// Get the middle position along the ring
		rPosition = ring/nRings + ( 1 / nRings  / 2 );

		// find the radius of the circle that crosses that distance value
		radius = findRadiusAt( profile, rPosition );

		// This point is then inside the sector we want in the ring we want
		xf = xm + radius * cos(2*PI*sector / nSectors + 2*PI / nSectors / 2 );
		yf = ym + radius * sin(2*PI*sector / nSectors + 2*PI / nSectors / 2 );

		// This allows us to create the selection
		doWand(xf, yf);
		
}

// We assume that the distances from inner to outer increase gradually (it is a distance)
// so this return the radius right after the value along the profile has exceeded the given rPosition
function findRadiusAt( profile, rPosition ) {
	//Array.show(profile);
	for( i=0; i< profile.length; i++ ) {
		if( profile[i] > rPosition ) return i;
	}
}