Hough Circular Transform analysis script to count the number of cells present in bright-field image.

hough_circle.ijm

//ImageJ macro written by Dominic Waithe for Caroline Scott. 2017.
//This macro uses the Hough circlular transform to count the number of cells present.
//Works very well, but artefacts on sample will lower accuracy.


///Parameters
//Radii to focus the hough-transform on
r_start = 6; //radii minimum size
r_stop = 8; //radii maximum size.
edge_mag = 70; //Threshold for edge magnitude, the lower this value the more of the edges will be captured.
noise_tol = 7; //This is the final noise tolerance for peaks. The lower this value the more peaks will be found.


////// Start of script.
gtitle =getTitle();
//Batch mode for stopping windows from opening.
setBatchMode(true);
run("Select None");
run("Duplicate...", " ");

//Reduce size of image for increased speed.
run("Bin...", "x=4 y=4 bin=Average");
run("RGB to Luminance");
run("Find Edges");

run("Clear Results");
title = getTitle();

//run through the radius
for(r=6;r<=8;r++){
selectWindow(title);
C = 2 * PI * r;
width = getWidth();
height = getHeight();

//I set the number of points to be sampled to be the same as the circumference
//Works with replicates, which we later remove.
m = floor(C)+1;
angles = newArray(0,2*PI);
nangles = Array.resample(angles,m);
cos_ang = newArray(m);
sin_ang = newArray(m);

//We calculate the position of the points

for (i=0;i<m;i++){
	ca = round(cos(nangles[i])*r);
	sa = round(sin(nangles[i])*r);
	if (i>0){
	//Check for replicates. If at least one is different, proceed.
	if(cos_ang[ct-1] != ca || sin_ang[ct-1] != sa ){
	cos_ang[ct] = ca;
	sin_ang[ct] = sa;
	ct ++;}}else{
		cos_ang[0] = ca;
		sin_ang[0] = sa;
		ct = 1;
}}
//Removes repetitions
cos_ang = Array.trim(cos_ang,ct);
sin_ang = Array.trim(sin_ang,ct);
//redefines the length.
m = cos_ang.length;

//Defines the array.
l = newArray(height*width);
c = 0;
for (y=0;y<height;y++){
	for (x=0;x<width;x++){
		if(getPixel(x,y) >edge_mag)
		{l[c]=1;}
		c++;
}}
//
newImage("radius_"+d2s(r,0), "32-bit black", width, height, 1);
c = 0;
for (y=0;y<height;y++){
	for (x=0;x<width;x++){
		if(l[c] ==1)
		{
		for (i=0;i<m;i++){
			px = getPixel(cos_ang[i]+x,sin_ang[i]+y);
			setPixel(cos_ang[i]+x,sin_ang[i]+y,1+px);
		}}
		c++;
}}}

run("Images to Stack", "name=stack title=radius use");
run("Z Project...", "projection=[Max Intensity]");
//run("Enhance Contrast", "saturated=0.35");
run("8-bit");
//run("FFT");
//run("Specify...", "width=50 height=50 x=256 y=256 oval centered");
//run("Clear");

//run("Inverse FFT");
run("Find Maxima...", "noise=45 output=[List]");





selectWindow(gtitle);
count = 0
for(i=0;i<nResults;i++){
xPt = getResult("X",i)*4;
yPt = getResult("Y",i)*4;
//Excludes the point which is routinely found on the scale bar.
if(xPt > 1710 && xPt < 1730 && yPt > 1230 && yPt < 1250){}else{
count +=1;
makePoint(xPt, yPt);
setKeyDown("Shift");
}
	}
//Print output count.
print(gtitle,"\tcount\t",count);
selectWindow(gtitle);
//setBatchMode(false);
//run("Restore Selection");