Segmentation and analysis of spherical caps on spherical particles

fiji
imagej
segmentation

#1

Hello, I have (many) SEM images of spherical glass nanoparticles with metal spherical caps on them (see example images attached). In a nutshell I would like to determine through image analysis:

  1. the fraction of all spherical glass particles possessing a cap
  2. for each spherical glass particle having a cap the cap coverage

I am very familiar with the basic fiji/imagej functions of filtering, thresholding and particle analysis (we use them a lot to analyse “bare” particles like the pure spherical glass particles shown. However, I have quickly established that these functions will not help me with the partially coated particles and a more sophisticated segmentation approach is needed (currently out of my field of expertise).

I wondered if these (inorganic) structures to be analyzed are analogous to biological structures for which an analysis recipe already exists? I would be grateful for any general or specific pointers as to how I can approach this.

Many thanks in advance!

Robin

!


#2

Good day!

for each spherical glass particle having a cap the cap coverage

How do you expect this to happen?
These are 3D structures and most often the caps are only partially visible.

(currently out of my field of expertise)

Mine too …

Regards

Herbie


#3

Many thanks for the feedback Herbie,

You are right that the caps are only partially visible. However, I believe the coverage fraction of even those caps can be obtained provided core sphere diameter and the outline of the visible part of the cap are known. Here we will make the assumption that the caps are perfect spherical caps. This is very reasonable because of the way these caps form (by crystal growth from a central nucleus at a constant rate in all directions). Thus, if we know the dimensions (perimeter coordinates) of the visible part of the cap we can maybe get back to the dimensions of the complete cap via geometrical analysis. Or failing that we can do a brute force fitting routine of many “test caps” to the real one. So I think that is not the problem. The challenge is (hence my posting here) to get a good hold on the perimeter coordinates of the caps / core particles through a workflow that can be applied to many images of these particles.

Robin


#4

the fraction of all spherical glass particles possessing a cap

I’d use a two step process for doing this but before you start you need to get rid of the uneven “illumination” that is highly disturbing.

  1. Get the position of all totally visible spheres
  2. Determine the fraction of spheres with caps

This two step process may be helpful because you may need different thresholds for 1. and 2.

Good luck

Herbie


#5

So I think that is not the problem

I think it is, at least with the spatial resolution of your second sample image.
You don’t explain which resolution is given, that of image one or two.

In any case you may get a statistical estimate of the coverage. I doubt that you will get a reasonable result “via geometrical analysis” from images having a spatial resolution of your second sample image.

Even your point one will be quite difficult to achieve but it seems possible with a certain accuracy.
What are your expectations and goals for the latter?

Regards

Herbie


#6

Hello @Robin_KT !

Interesting task.
You have several difficulties:

  • resolution, but this is not the most important one (and possibly you have a high resolution setting on your SEM for new images)
  • uneven illumination,
  • the fact that spheres overlap.

For uneven illumination, try to grab an non focused image for each image you take, and use it for background subtraction.
In the following image I was able to segment some of the spheres by first using a black-top hat operation (since spheres are separated by dark layers), then treshold, then watershed to separate touching elements.


This works well for spheres that lie in a single layer. For the other ones (bottom right), that does not work because the spheres are separated by “less dark” features (actually, the spheres that lie undernath). I did it very fast, so finer adjustment of the parameters may provide better results.
But I would suggest you try to take pictures of spheres in single layers, or get rid of “difficult parts” in the image. I know the latter is not satisfying, but anyway, you’re not going to process hundreds of images automatically, I guess. Since your objects should have the same diameter, you can also choose to exclude the ones that are too large.

After the spheres have been segmented, it think becomes easier, but tedious. Possibly you would like to evaluate the gray level distribution in each sphere: spheres with a cap should show high maximum gray values. If you’re lucky the histogram for each individual sphere will allow you to segment it into cap/non cap. At that point, a tool like WEKA segmentation can help.

Then you have the issue of statistics and representativity: what is the fraction of spheres that have a cap which cannot be observed ? And what about multiple caps? Good luck with that! Anyway, the point is to provide some kind of uncertainty associated with your measurement.

Regards,

Nicolas


#7

Nicolas,

most of what you write has already been proposed by me.

This is what I would start from:

Regards

Herbie


#8

Sorry, I read too fast. I didn’t realize you already mentionned the top-hat operation, the watershed, the flat-field correction of SEM images, and the need of statistical correction for caps occlusion, in your answers.
Regards,

Nicolas


#9

Robin,
I know you will rectify the uneven illumination, but since you mentioned that the spheres were a glass and the caps were metallic, could you try using a different color or type of lighting?
I am very curious as to a way to solve this problem.
Bob


#10

Hi Robin,

I think the WEKA Trainable Segmentation might be a good tool for this task. First apply pseudo flatfield correction to even out your image, and then train WEKA with typical background, sphere and cap elements. I’ve attached the result form a few minutes work.

yempski

pseudo-flat


#11

Many thanks @Herbie, @Nicolas and @yempski for your suggestions and examples of what can be done with these images. I’ve now had a chance to return to this problem and consider your comments.

First a comment on spatial resolution. The images I provided were just examples to give you an idea of the problem. I am new to this forum and I had not expected that someone would start actually processing these images (thanks very much that you are prepared to do that!). So we have images covering a range of magnifications and we have the possibility to make new images, lots of them by automated acquisition if necessary. The lower magnification image I provided is obviously too low resolution and we would never use that for analysis. I just wanted to show the kind of distribution that we get. The higher magnification image has an acceptable resolution but is complicated by the particles lying on top of each other. We can usually find sample regions with isolated particles or single layers of touching particles. Indeed, if we can use image analysis to obtain the information we need then it will not be a problem to improve our sample preparation and imaging to get suitable images. This also regards the issue of background removal mentioned by several of you (@Nicolas – I like this suggestion for getting an equivalent out-of-focus image for each image in order obtain the background – I will try that).

Regarding the results of the analysis I am fully aware that I am going to get statistical data with some uncertainties (e.g. due to incorrectly identified bare-particles having a hidden patch). Fortunately image analysis is not the only tool we have for analysing these particles. They have optical spectra that are highly sensitive to the cap number and size. Also we have analytical centrifugation that is sensitive to the cap-to-core mass ratio. These techniques also have their shortcomings. To a certain extent, the analysis of these micrographs is in order to validate these “indirect” techniques.

I have tried to put into practice some of the new tricks I have learnt from your responses. Firstly I have taken a more reasonable image in terms of resolution and no overlapping particles. Here is the workflow that I applied with the images following.

  1. Apply kuwahara filter (Sampling Window 5) – to remove noise

  2. Apply white Top Hat filter (Element = Disk, Radius = 70px) [@Nicolas – thanks for introducing me to this family of filters which I did not know. I chose the filter conditions very empirically to give a good contrast between the particles and background. Is this OK or is there a more systematic approach?]

  3. Apply threshold and watershed algorithm

  4. Analyze particles and obtain the Mean, Standard Deviation of the level and the Maximum of the level for each particle

From visual inspection, 21 of the 46 found particles have caps. I looked to see if I could get this number from the approach suggested by @Nicolas (Thanks for pointing out that the presence or absence of a cap should be indicated by the maximum level in the particle. It seems seems so obvious now)

Here are plots of the mean level, the standard deviation the max level and the standard deviation multiplied by the max level for all particles.

Clearly the mean level is pretty useless. The Standard Deviation and Max Level are better, each revealing 16 particles with caps. In fact, by multiplying these together (see bottom right plot) I could even identify 17 particles with caps. Considering that 3 of the particles with caps that are not detected have uncharacteristically small caps (and thus can be neglected) this is really an excellent result. Possibly this could be further improved by levelling the background which I did not do.

The next task is to find the caps. I have tried a second thresholding step (suggested by @herbie and @nicolas) locally on each found capped particle. This works well for some, but not all patches, presumably because the histogram doesn’t contain two clear peaks.

@Nicolas and @yempski have introduced me to another new tool: WEKA Trainable Segmentation. I have tried to do this for the image processed above. Firstly I used the ROI data to completely wipe the background

Here the WEKA processing:

The WEKA routine seems to find most of the caps. However, it also finds parts of the bare spheres including a thin line around their circumference. What I am doing wrong?

Any other tips on how to get a hold on these caps (any experience with edge detection?) much appreciated.

Robin


#12

@smithrobertj The images were made in an SEM (secondary electron imaging mode). Hence I cannot easily change the color or type of illumination. I could possibly take a backscattered electron image in parallel to the SE one. The metal caps would then appear brighter. However, that imaging mode has lower resolution and less topological detail.


#13

Good day Robin,

let’s assume you have a segmetation like this:

Then the following ImageJ-macro gives you the number of spheres with caps:

// imagej-macro "spheres&caps" (Herbie G., 13. Nov. 2018)
requires( "1.52h" );
setBackgroundColor(0, 0, 0);
setOption("BlackBackground", true);
orig=getImageID();
run("To ROI Manager");
roiManager("Combine");
roiManager("Add");
run("Clear Outside");
n=roiManager("count")-1;
roiManager("Select", n);
roiManager("Delete");
setBatchMode(true);
for ( i=0; i<n; i++ ) { 
   selectImage(orig);
   roiManager("select", i)
   run("Duplicate...", "title="+(i+1));
   run("Variance...", "radius=2");
   run("Enlarge...", "enlarge=-3");
   run("Clear Outside");
   setAutoThreshold("Yen dark no-reset");
   run("Convert to Mask");
   run("Analyze Particles...", "size=150-Infinity exclude summarize");
   close();
}
run("Select None");
selectWindow("ROI Manager");
run("Close");
m=Table.getColumn("Count","Summary");
caps=0;
for ( i=0; i<m.length; i++ ) if (m[i]>0) caps++;
Table.deleteColumn("Median");
Table.deleteColumn("Mode");
Table.deleteColumn("Mean");
Table.deleteColumn("%Area");
Table.deleteColumn("Average Size");
Table.deleteColumn("Total Area");
Table.set("Total Caps", 0, caps);
Table.update("Summary");
setBatchMode(false);
exit();
// imagej-macro "spheres&caps" (Herbie G., 13. Nov. 2018)

Paste the above macro code to an empty macro window (Plugins >> New >> Macro) and run it with the following ZIP-compressed image open in ImageJ:
Ag@SiO2_AP53_11_02_WTH_disc_70px_OL.tif.zip (626.6 KB)
Please use this image (after un-zipping it), not the above posted PNG-version!

There are three critical numerical parameter values that must be adapted according to the magnification (sphere size):

  1. run(“Variance…”, “radius=2”);
  2. run(“Enlarge…”, “enlarge=-3”);
  3. run(“Analyze Particles…”, “size=150-Infinity exclude summarize”);

HTH

Herbie


#14

@Herbie, that is a very neat thing you did there. I was thinking up all sorts of tricks that I could perform on the individual particle histograms in order to get at the last remaining uncounted caps. It would have been a lot of code. Your method, to convert the edge of the cap within the particle into a particle to be found is much more elegant. I tried your code on another of my images. It was quite successful but not 100%. Some cups were not detected because of weak contrast between the cup and particle background. I was going to start tweaking the parameters you mentioned. Then I noticed that you seem to have performed an operation on my image before applying your code. The contrast appears better on your image than on mine. Was this the “enhance constrast” function? I tried that function but did not get the same image that you had.

Thanks again for this very cool piece of code!

Robin


#15

Hi Robin,
I’m trying to figure why you are getting such a preferential direction of the images. Now I know that you are acquiring via SEM I will look into this from a different direction. Also if you send in a backscattered image, it may prove more useful for what you are wanting.
Bob


#16

Hello Robin,
Please take a look at the attached images and files. They are all based on the 1st Tiff image you submited earlier. The files are all summarized at the end of each file and I’m sure you can handle all of the math required (if any is needed).
I’m running a little late so I will write up the whole precedure tomorrow. It uses no special workings and can be applied to any of the future images you may take for Analysis.
Type at you later,
Bob

Just Caps Results.csv (11.1 KB)
Whole Results.csv (7.8 KB)


#17

Good day Bob,

AFAIK your result is not up to what the OP and I have achieved already.

There are 46 spheres that are not touching the image border and 21 of them have caps. My macro gives the correct number and identifies these spheres.
According to your first image, your approach does not.

Regards

Herbie


#18

Good day Robin,

here is the ImageJ-macro that I used for the pre-processing:

// imagej-macro "segmentSpheres" (Herbie G., 13. Nov. 2018)
requires( "1.52h" );
setOption("BlackBackground", true);
setBatchMode(true);
run("32-bit");
run("Median...", "radius=2");
run("Subtract Background...", "rolling=100 sliding");
run("Duplicate...", "title=temp");
setAutoThreshold("Default dark no-reset");
run("Convert to Mask");
run("Watershed");
run("Analyze Particles...", "size=100-Infinity circularity=0.75-1.00 exclude add");
close("temp");
run("From ROI Manager");
setBatchMode(false);
exit();
// imagej-macro "segmentSpheres" (Herbie G., 13. Nov. 2018)

It could be further improved.

Regards

Herbie


#19

First an update: the first challenge in my original post has been met - I now have a very robust macro which can take a raw image and count the particles with caps!

Now the detail: Many thanks @Herbie for your upstream code. With just the 2px median filter and 100px background subtraction I still could not get the same result as you. The median filter makes the edges of the caps blurry and then the variance filter applied later on doesn’t do its job so well. I worked around this issue by adding an unsharp mask step after the median filter and background subtraction. There is probably a more elegant solution, but it worked satisfactorily for me, and better than with the processing methods I used earlier.

I have packaged all this into a longer code which enables an open raw image to be processed and analyzed with a few clicks and all intermediate images (including the cutouts of the capped and non-capped particles) and data saved in a dedicated directory. This automates a lot of things, like the determination of the pixel size of the minimum cap during the cap search routine, which obviously will depend on on the image resolution. This code makes it easy to override such automated settings though. This is necessary because, while I have (by trial and error) found the settings that find most caps on most images at medium magnification (about 25-75 particles in the image), if I go to really low or high magnification some settings have to be tweaked. To get quick visual feedback on the result I have got the macro to produce images with just the capped or just the bare particles. Then the miscounted particles can be immediately seen.
I also found it useful, in the case that too many bare particles were being counted as capped, that the results from the first search are refined by looking at the average area% of found “caps”. This slightly raises the bar and after one or two attempts one can find a good condition which leads to just the capped particles being found. This is particularly useful for very low magnification images (>300 particles per image) and thus valuable for getting good statistics without analysing hundreds of images!

But first the code (in case it is of any use to anyone):

//Macro to determine cap yield (Robin KT (with significant input from Herbie G), November 2018)
requires( "1.52h" );
//run("Input/Output...", "jpeg=90 gif=-1 file=.xls save_column");
setOption("BlackBackground", false);
// close any open windows
if (isOpen("ROI Manager")) { 
       selectWindow("ROI Manager"); 
       run("Close"); 
   }
if (isOpen("Results")) { 
       selectWindow("Results"); 
       run("Close"); 
   } 
if (isOpen("Summary")) { 
       selectWindow("Results"); 
       run("Close"); 
   } 

// get file information from image
olddir = getDirectory("image");
namese = getInfo("image.filename");
if (namese=="") exit ("name not available");
oldpath = olddir + namese;
imagetit = replace(namese, ".tif", "");
newdir = olddir + imagetit + File.separator;
dircount=d2s(0,0);

// create directory for data
i=1;
while (File.exists(newdir) ==1){
dircount = d2s(i, 0);
newdir = olddir + imagetit + "_" + dircount + File.separator;
i = i+1;
}
newpath = newdir + imagetit;
File.makeDirectory(newdir);

//Dialog to set information on image scale and pre-processing
Dialog.create("Set image scale and define preprocessing");
items = newArray("Raw SEM image from ZEISS microscope (Tifftags plugin required!)", "Image already scaled");
Dialog.setInsets(5, 20, 3)
Dialog.addRadioButtonGroup("Choose scale option", items, 2, 1, "Image already scaled");
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addMessage("Choose processing options");
Dialog.addCheckbox("Median filter", true);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Radius in pixels", 2);
Dialog.setInsets(5, 30, 3)
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addCheckbox("Kuwahara filter", false);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Sampling window width (must be odd)", 5);
Dialog.setInsets(5, 30, 3)
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addCheckbox("Subtract background (auto radius is 1.8 times average particle radius)", true);
Dialog.setInsets(5, 30, 3)
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addCheckbox("Unsharp mask", true);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Radius", 1);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("mask", 0.6);
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addCheckbox("Enhance Contrast (Global)", false);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Saturated pixels (%)", 0.3);
Dialog.addMessage("-------------------------------------------------------------------------------------");
Dialog.addCheckbox("Enhance Contrast (Local)", false);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Blocksize", 127);
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Histogram", 256); 
Dialog.setInsets(5, 30, 3)
Dialog.addNumber("Maximum", 2.5);
Dialog.show;
scaledyet = charCodeAt(Dialog.getRadioButton, 1);
procopts = newArray(14);
procopts[0] = Dialog.getCheckbox(); //median filter
procopts[1] = Dialog.getNumber();//median filter radius
procopts[2] = Dialog.getCheckbox(); //kuwahara filter
procopts[3] = Dialog.getNumber(); //kuwahara filter sampling window width
procopts[4] = Dialog.getCheckbox(); //Subtract background
procopts[5] = Dialog.getCheckbox(); //Unsharp mask
procopts[6] = Dialog.getNumber(); //Unsharp mask - radius
procopts[7] = Dialog.getNumber(); //Unsharp mask - mask
procopts[8] = Dialog.getCheckbox(); //Enhance contrast (global)
procopts[9] = Dialog.getNumber(); //Enhance contrast (global) - saturated pixels
procopts[10] = Dialog.getCheckbox(); //Enhance constrast (local)
procopts[11] = Dialog.getNumber(); //Enhance constrast (local) blocksize
procopts[12] = Dialog.getNumber(); //Enhance constrast (local) histogram
procopts[13] = Dialog.getNumber(); //Enhance constrast (local) maximum

//Dialog about binary processing and particle segmentation
Dialog.create("Binary processing and particle search");
items = newArray("Yes", "No");
Dialog.addRadioButtonGroup("Do you want to execute the watershed algorithm?", items, 2, 1, "Yes");
Dialog.addRadioButtonGroup("Do you want to fill holes inside the particles?", items, 2, 1, "No");
Dialog.addNumber("Core Particle Diameter in nm?", 500);
Dialog.addNumber("Minimum Circularity:", 0.75);
Dialog.show()
ws_query= charCodeAt(Dialog.getRadioButton, 1);
fill_query=charCodeAt(Dialog.getRadioButton, 1);
corediameter = Dialog.getNumber();
mincirc = Dialog.getNumber();

//Scaling the image if necessary
if (scaledyet < 100) // Not yet scaled - Set from tiff header
{
tagnum=34118;
tag = call("TIFF_Tags.getTag", oldpath, tagnum);
pixsizestr = substring(tag,9,17);
exp1str = substring(tag,19,20);
exp1 = parseInt(exp1str);
exp2str = substring(tag,20,21);
exp2 = parseInt(exp2str);
exp3str = substring(tag,21,22);
exp3 = parseInt(exp3str);
expl = 100*exp1 + 10*exp2 + exp3;
pixsizenum = parseFloat(pixsizestr);
pixsize = pixsizenum*pow(10, -1*expl);
pixsizenm = pixsize*1E9;
pixsizereal = pixsizenm;
//print(pixsizereal);
setVoxelSize(pixsizereal, pixsizereal, 1, "nm");
save(newpath + "_00_raw.tif");

//crop off ZEISS scalebar
makeRectangle(0, 0, 1024, 708);
run("Crop");
}

else getPixelSize(unit, pixsizereal, pixsizereal); // Already scaled - get from image
save(newpath + "_00_raw.tif");

//calculate some useful stuff

pixdiameter = corediameter/pixsizereal; //diameter of the core particle in pixels
pixarea = PI*pixdiameter*pixdiameter/4; //Area of the core particle in pixels
smallest = corediameter - 0.25*corediameter; //smallest expected particle in nm
largest = corediameter + 0.25*corediameter; //largest expected particles in nm
smallestarea = round(PI*smallest*smallest/4); //smallest expected particle in nm^2
largestarea = round(PI*largest*largest/4); //largest expected particle in nm^2
bkgdradius = round(0.9*pixdiameter); // determine the background removal moving ball radius according to particle size

//Do selected processing and save processed image

run("32-bit");
if (procopts[0] == 1) run("Median...", "radius=" + procopts[1]);
if (procopts[2] == 1) run("Kuwahara Filter", "sampling=" + procopts[3]);
if (procopts[4] == 1) run("Subtract Background...", "rolling="+ bkgdradius +" sliding");
if (procopts[5] == 1) run("Unsharp Mask...", "radius=" + procopts[6] + " mask=" + procopts[7]);
if (procopts[8] == 1) run("Enhance Contrast...", "saturated=" + procopts[9]);
if (procopts[10] == 1) run( "Enhance Local Contrast (CLAHE)", "blocksize=" + procopts[11] + " histogram=" + procopts[12] + " maximum=" + procopts[13]);
newname = getTitle();
save(newpath + "_01_proc.tif");

//Do auto thresholding
setAutoThreshold("MinError dark no-reset");
getThreshold(lower, upper);
if (lower == 0) thss = d2s(upper, 0);
else thss = d2s(lower, 0);
run("Convert to Mask");

//carry out watershed algorithm if selected (default = Yes)
//carry out fill holes algorithm if selected (default = No)
if (ws_query < 105) 
{
run("Watershed");
if (fill_query < 105) run("Fill Holes");
//save watershed image
save(newpath + "_02_ws_threshold=" + thss + "_" + dircount + ".tif");
}

else 
{
if (fill_query < 105) run("Fill Holes");
save(newpath + "_threshold=" + thss + "_" + dircount + ".tif");
}


//Find particles

run("Set Measurements...", "area centroid fit area_fraction redirect=None decimal=5");

run("Analyze Particles...", "size="+ smallestarea +"-"+ largestarea +" circularity="+ mincirc +"-1.00 show=Ellipses display exclude add in_situ");

//close thresholded image and reopen original and set scale
selectWindow(newname);
close();
open(newpath + "_01_proc.tif");

//import particle outlines to original processed image
run("From ROI Manager");

//request manual removal of unwanted particles
waitForUser("Remove unwanted regions of interest (ROIs) by clicking them with the mouse and hitting 'backspace' on the keyboard \n When you are done click OK in this pop-up.");

//Update the ROI manager and save the image showing ROIs
run("To ROI Manager");
save(newpath + "_03_roi_" + dircount +".tif");

if (isOpen("Results")) { 
       selectWindow("Results"); 
       run("Close"); 
   } 

//Measure the segmented particles and find the circle equivalent area diameter xc
roiManager("Deselect");
roiManager("Measure");
foundcount=getValue("results.count");
run("Summarize");
foundmean=getResult("Area", foundcount);
foundmeanpix = round(foundmean/(pixsizereal*pixsizereal));
xc= sqrt(4*foundmean/PI);
patchcutoff = round(0.0075*foundmeanpix); //Determines suggested value for minimum cut off for cap size to be used below

selectWindow("Results"); 
       run("Close"); 

//close all images
while (nImages>0) { 
          selectImage(nImages); 
          close(); 
      }

//dialog to get information on cap search process      
Dialog.create("Cap Search Parameters");
	Dialog.addMessage(foundcount + " particles found");	
	Dialog.addMessage("Mean circle-equivalent diameter is " + xc + " nm or " + xc/pixsizereal + " pixels.");
	Dialog.addMessage("Mean particle area is " + foundmeanpix + " pixels");	
	Dialog.addNumber("Variance filter radius in pixels (increase for higher magnification images): ",2);
	Dialog.addNumber("Selection boundary contraction amount in pixels: ",4);
	Dialog.addNumber("Minimum number of pixels for identification of cap (increase for lower magnification images):", patchcutoff);
	Dialog.addNumber("Minimum area fraction (%) for found object to be considered a cap (if some bare particles are counted as being capped this can be raised slightly):", 0.71);
Dialog.show();
varradius = Dialog.getNumber();
contract = Dialog.getNumber();
cutoff = Dialog.getNumber();
minareafrac = Dialog.getNumber();

//open the segmented image again and remove the scale (because the cap search is all done in pixels!)
open(newpath + "_03_roi_" + dircount +".tif");
run("Set Scale...", "distance=0 known=0 pixel=1 unit=pixel");

// Here starts the imagej-macro "spheres&caps" (Modified from macro created by Herbie G., 13. Nov. 2018)
setBackgroundColor(0, 0, 0);
setOption("BlackBackground", true);


//select all segmented particles and remove background
orig=getImageID();
run("To ROI Manager");
roiManager("Combine");
roiManager("Add");
run("Clear Outside");
n=roiManager("count")-1;
roiManager("Select", n);
roiManager("Delete");

//Define directories for cutout images
cutoutdir = newdir + imagetit + "_" + dircount + "_extracted" + File.separator;
capdir = cutoutdir + imagetit + "_" + dircount + "_caps" + File.separator;
nocapdir = cutoutdir + imagetit + "_" + dircount + "_no_caps" + File.separator;
File.makeDirectory(cutoutdir);
File.makeDirectory(capdir);
File.makeDirectory(nocapdir);

//The search for caps (by Herbie G., 13. Nov. 2018)
setBatchMode(true);
for ( i=0; i<n; i++ ) { 
   selectImage(orig);
   roiManager("select", i)
   run("Duplicate...", "title="+(i+1));
   save(cutoutdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif"); //Save cut out image for future inspection
   run("Variance...", "radius=" + varradius);
   run("Enlarge...", "enlarge=-" + contract); //was 3 in original
   run("Clear Outside");
   setAutoThreshold("Yen dark no-reset");
   run("Convert to Mask");
   save(cutoutdir + imagetit +"_ROI=" + i+1 + "_Proc_" + dircount + ".tif"); //Saves binary image for future inspection
   run("Analyze Particles...", "size=" + cutoff + "-Infinity exclude summarize");
   close();
						}
//Close the ROI Manager
run("Select None");
selectWindow("ROI Manager");
run("Close");

//Get useful parameters from cap search which might help filter results further
m=Table.getColumn("Count","Summary");
p=Table.getColumn("%Area","Summary");
q=Table.getColumn("Major","Summary");
r=Table.getColumn("Minor","Summary");
caps=0;

//Carry out secondary filtering of results based on %area (cutoff defined in above dialog)
for ( i=0; i<m.length; i++ ) { 
if (m[i]>0 && p[i]/m[i]>minareafrac ) {
	caps++; 
	Table.set("Cap", i, "1");
	Table.update("Summary");
							}
//else if (m[i]>0 && p[i]/m[i]>2 && p[i]/m[i]<4.2 && q[i]/r[i] < 1.5) { 	//possible third filtering by allowing some particles with %area between 2 and 5 to be counted if their bright region has low aspect ratio (i.e. probably a cap)
//caps++;
	//Table.set("Cap", i, "1");
	//Table.update("Summary");						
		//				}
else {
	Table.set("Cap", i, "0");
Table.update("Summary");
		}
								}

//Set up arrays with the indexes of particles with and without caps and move cutout images into two different directories for easy future analysis
truecaps=newArray(caps);
falsecaps=newArray(m.length-caps);
c=Table.getColumn("Cap","Summary");
i_true=0;
i_false=0;
for ( i=0; i<m.length; i++ ) { 
if (c[i]==1) {
	truecaps[i_true]=i;
	i_true++;
	File.copy(cutoutdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif",capdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif");
	File.delete(cutoutdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif");
	File.copy(cutoutdir + imagetit +"_ROI=" + i+1 + "_Proc_" + dircount + ".tif",capdir + imagetit +"_ROI=" + i+1 + "_Proc_" + dircount + ".tif");
	File.delete(cutoutdir + imagetit + "_ROI=" + i+1 + "_Proc_" + dircount + ".tif");
			}
else {
	falsecaps[i_false]=i;
	i_false++;
	File.copy(cutoutdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif",nocapdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif");
	File.delete(cutoutdir + imagetit + "_ROI=" + i+1 + "_" + dircount + ".tif");
	File.copy(cutoutdir + imagetit +"_ROI=" + i+1 + "_Proc_" + dircount + ".tif",nocapdir + imagetit +"_ROI=" + i+1 + "_Proc_" + dircount + ".tif");
	File.delete(cutoutdir + imagetit + "_ROI=" + i+1 + "_Proc_" + dircount + ".tif");
	}
								}
close("Log");
//Save images showing just the capped and non-capped particles (for visual appraisal of the success of the cap search algorithm)
open(newpath + "_03_roi_" + dircount +".tif");
run("To ROI Manager");
roiManager("select",truecaps);
roiManager("Combine");
roiManager("Add");
run("Clear Outside");
n=roiManager("count")-1;
roiManager("Select", n);
roiManager("Delete");
save(newpath + "_04_roi_caps_" + dircount +".tif");
run("Remove Overlay");
save(newpath + "_04_found_caps_" + dircount +".tif");

open(newpath + "_03_roi_" + dircount +".tif");
run("To ROI Manager");
roiManager("select",falsecaps);
roiManager("Combine");
roiManager("Add");
run("Clear Outside");
n=roiManager("count")-1;
roiManager("Select", n);
roiManager("Delete");
save(newpath + "_05_roi_no_caps_" + dircount +".tif");
run("Remove Overlay");
save(newpath + "_05_found_no_caps_" + dircount +".tif");
while (nImages>0) { 
          selectImage(nImages); 
          close(); 
      }
setBatchMode(false);

//Save summary of processing and particle search to file
f = File.open(newpath + "_07_Summary_" + dircount +".txt");
print(f,"Image processing:");
if (procopts[0] == 1) print(f,"Median filter used with radius = " + procopts[1] + " pixels");
if (procopts[2] == 1) print(f,"Kuwahara filter used with sampling window width = " + procopts[3]);
if (procopts[4] == 1) print(f,"Background substract used with rolling ball radius = " + bkgdradius + " pixels or " + pixsizereal*bkgdradius + " nm");
if (procopts[5] == 1) print(f,"Unsharp mask used with radius = " + procopts[6] + ", mask = " + procopts[7]);
if (procopts[8] == 1) print(f,"Enhance contrast (global) used with saturated pixels (%) = " + procopts[9]);
if (procopts[10] == 1) print(f,"Enhance constrast (local) used with blocksize = " + procopts[11] + ", histogram = " + procopts[12] + ", maximum = " + procopts[13]);
print(f," ");
print(f,"Particle Search:");
print(f,"Total number of particles found: " + foundcount);
print(f,"Mean diameter of particles found: " + xc + " nm");
print(f," ");
print(f,"Cap Search:");
print(f,"Variance filter applied with radius: " + varradius + " pixels");
print(f,"ROI contracted by: " + contract + " pixels");
print(f,"Minimum number of pixels for identification of cap: " + cutoff);
print(f,"Minimum area fraction (%) for found object to be considered a cap: " + minareafrac);
print(f," ");
print(f,"Results:");
print(f,"Total number of particles found: " + foundcount);
print(f,"Total number of caps found: " + caps);
print(f,"Cap yield: " + 100*caps/foundcount + "%");
File.close(f);
//Save the raw statistics on the particles with found caps (can help to refine the secondary filtering as necessary)
Table.rename("Summary","Results");
saveAs("Results", newpath + "_06_Rawstats_" + dircount +".csv");
close("Results");
exit();
// N.B. includes parts of imagej-macro "spheres&caps" (Herbie G., 13. Nov. 2018)

This code will take the image I shared in my second big post above and with all the default settings it will return the 21 capped particles. Here is the raw image (this time as a zipped tif and including the scale which is needed for the code to work):
Ag@SiO2_AP53_11.zip (593.8 KB)
As mentioned above the code will output lots of things but particularly useful are 2 images which should show only the capped and bare particles respectively. For the above image these then look like:

Original:

Capped particles only (75 particles found)

Non-capped particles only:

Another example is this image which contains 249 particles:
Ag_SiO2_AP31_03.zip (569.6 KB)

Again the code can be run with all the default settings apart from the second cap refinement (Minimum area fraction (%) for found object to be considered a cap) which needs to be raised to 2.5. Then the following images are returned:
Original:

Capped particles only:

Non-capped particles only:

Clearly there are a few very small caps that were counted as bare particles. But as noted before, these don’t bother me too much.

So overall, this is a great code and I am very thankful for the hints and contributions that helped get it set up so quickly! Regarding my 2nd question I will write a separate post after this one.

Robin


#20

The second task I want to do with these capped particle images is to obtain the cap coverage fraction for each capped particle. The code I posted above outputs, amongst other things, images of the individual capped particles.
Let’s take the following raw image:

By applying the cap-search code we get these “cutouts” (tifs in this zip: Ag_SiO2_AP31_05_cutouts.zip (378.3 KB)):

Ag_SiO2_AP31_05_ROI%3D5_1 Ag_SiO2_AP31_05_ROI%3D8_1 Ag_SiO2_AP31_05_ROI%3D11_1 Ag_SiO2_AP31_05_ROI%3D13_1 Ag_SiO2_AP31_05_ROI%3D18_1 Ag_SiO2_AP31_05_ROI%3D20_1 Ag_SiO2_AP31_05_ROI%3D23_1 Ag_SiO2_AP31_05_ROI%3D24_1 Ag_SiO2_AP31_05_ROI%3D25_1 Ag_SiO2_AP31_05_ROI%3D26_1 Ag_SiO2_AP31_05_ROI%3D30_1 Ag_SiO2_AP31_05_ROI%3D32_1

Simply put, I need the perimeter coordinates of each cap. As I have noted in an earlier post, thresholding works well for only some of the found caps. As @herbie has shown, a simple edge detection (variance filter) could distinguish the caps from everything else sufficiently enough to count the caps. However, it often misses some parts of the caps so I couldn’t find a way to use that to get the complete outline of the cap. So I am wondering if there is an even more sophisticated edge detection algorithm that could give me these perimeters? As ever, many thanks in advance for any hints and tips!

Robin