The centroids of objects in an image have been found. It has then been cropped to a certain area. IS it possible to tell if more than one centroid is in a cropped image?

I have the centroids of all the objects of an image. I have than cropped a 288 by 288 pixel area around the images (the size of the area is to try and isolate one centroid uniformly). Is it possible to tell if I have more than one centroid within my cropped image.

Here is the code I have to find the centroid and crop it. I am using Otsu's method (graythresh the im2bw) to threshold because it is predetermined method I am suppose to use, even though I know there is better ways. Code at the beginning to take out all objects with less than 120 pixels connection is there as well just fyi.

The overall goal is to individually isolate the object that each centroid is on and to isolate it and find the area of it, after thresholding. Perhaps there is a better way.

if true 
tic;
load image=originalimage;
binaryimage1=imregionalmax(loadimage);
binaryimage2=bwareaopen(binaryimage1,120);
labeledimage=bwlabel(binaryimage2);
props=regionprops(binaryimage2,'Centroid);
allCentroids=[props.Centroid];
xcentroids=allCentroids(1:2:end);
ycentroids=allCentroids(2:2:end);
numberofobjects=size(props);
for m=1:numberofobjects
 xmin=xcentroids(m);
 ymin=ycentroids(m);
 cropthat=imcrop(loadimage, [xmin-144, ymin-144, 288, 288]);
 pixellevel=graythresh(cropthat);
 levelbw=im2bw(cropthat,pixellevel);
 regionprops2=regionprops(levelbw,'area');

end

Answer 1

[Image Analyst]

Upload your image so we can see what this strange algorithm is doing. I say strange because it depends on getting a binary image from regional max but only if the regional max is 120 pixels in area or more. This happens only for large perfectly flat plateaus.

[Image Analyst]

Heath, try this and get back to me. I'm still not sure what you want to measure, so I did the regional max thing and counted blobs in that, then did the same thing for the original image using threshold. Of course the cropped subimages from the original and regional max images will each give a different number of blobs since they arrived at the binary image via different algorithms. Not sure why you did that, but anyway... at least this is a more informative GUI I've developed for you.

See attached test3.m.

[Krispy Scripts]

Wow thank you so much. These are a bunch of neurons with there axons extending small areas around them. I am trying to find the area of each individual neuron. I have looked over your blob demo of image analysis and I think it would be better to do your way.

Do you think it would be better than what I've done if I did in the same way you did within your tutorial:

1) threshold with im2bw or logical operation 2)fill holes with imfill function 3) find objects with bwlabel or bwconncomp 4)then use region props for area, cell perimeter and other relative measurements?

The reason I did it the way I did was to try and find each neuron(find centroids), then use the original image with those locations, threshold it(Otsu's method), and then measure the area, perimeter and spread of neuron extensions of each neuron(imcrop/regionprops).

With neuron staining there tends to be background and I know whatever I choose will not exactly get the area and other measurements right exactly, but I think there is a better way then I had in mind.

[Image Analyst]

No, sorry, but actually I'm not a big fan of your algorithm. Or maybe I just don't know what you want to do. The main issues I have with it are using imregionalmax() to find centroids of bright blobs. imregionalmax() finds peaks and the problem is that each connected bright blob may have several peaks in it. And each peak location does not necessarily occur at the weighted centroid (center of gravity). Then you call bwareaopen() on it which requires that the blob be 120 pixels or bigger. For imregionalmax() to give a blob that big, you'd need a flat plateau that's 120 pixels big and I don't think that would be very likely, though it does seem to find a lot of them. However I think you could be missing a lot of other blobs just because they didn't have a 120 pixels or larger flat region.

But then you use a regular Otsu-style threshold on the original image in a 288x288 square region around the peak centroids. Well I don't see any reason for the 288 size, and the centroid from imregionalmax() is not necessarily at the same centroid are you'd find from using an Otsu threshold. So you have a set or imregionalmax centroids and a set of centroids from graythresh. So you have two counts instead of one count. Which is "right"? Are they both "right"? What do you do with these two counts? It looks like you're ignoring the imregionalmax centroids and just using the graythresh centroids and then you said you want to count the centroids (graythresh ones I supposed) that are in this square 288x288 window around the regionalmax centroids. For one thing, I don't know what theoretical meaning this count has, and why use a shape of that size? Why not a circle or hexagon or something? I have no idea.

Now you know that graythresh() will give different thresholds as you move around the image don't you? Now perhaps you want that and perhaps you don't. If the background level changes as you move around, and you're interested in the stuff brighter than the local background then that would be good. However if the background is flat, like you've already corrected/flattened the background, then maybe you can just use a global threshold. Or maybe you can use adapthisteq() to flatten your image and then use a global threshold after that.

In general I don't like graythresh() or the Otsu threshold. Please please read Steve's three part blog on thresholding:

http://blogs.mathworks.com/steve/2016/05/16/image-binarization-new-r2016a-functions/

http://blogs.mathworks.com/steve/2016/06/14/image-binarization-otsus-method/

http://blogs.mathworks.com/steve/2016/07/25/adaptive-thresholding-for-binarization/

I have two methods I prefer much more than the Otsu method. The first is a fixed threshold. Why use that when virtually the whole world of novice image analysis people think you should use automatic thresholding? Well, because often they're wrong and they don't know how to design imaging systems or software. If you have good control of your lighting (it doesn't change brightness from image to image), and you have good background correction (so the brightness of an object doesn't depend on where it is in the field of view like it will for virtually all optical lens systems), then you can use a fixed threshold. Why fixed? What advantage is there to that? Well what if you want to find ANY percentage of "stuff" in the foreground - anything from 0% to 50% to 100%. Well if you have around 50% and have a nice bimodal histogram with well separated humps, then Otsu is fine. But what about if it's a bimodal histogram with a single hump? Otsu will be forced to find a threshold in the hump and it will probably be wrong since it's likely you have only foreground or only background and not some pixels of each. The advantage of a fixed threshold is that it can robustly detect any amount of foreground and background from 0% to 100%. You will not get that with Otsu. Otsu will fail and tell you there is something there when there is NOT.

The other threshold that is great is a triangle threshold. Steve is considering adding this to the Toolbox I think and that would be great. It's claim to fame is finding the threshold for unimodal but skewed histograms (a histogram with a long tail). Very very often distributions of almost anything (gray levels, particle areas or diameters, perimeters, etc.) will have a skewed distribution that looks like a log-normal shape. Triangle thresholding is great in that it finds the "corner" of the histogram on the skewed/long-tail side. And what it finds correlates very well with what people pick when they do it interactively, like you can do with my interactive thresholding application (give it a try).

Anyway, that's how I assess things at the moment. I still don't know why you're doing all that and I think that to get the area and extensions, etc. there are better ways. Like see my attached demo to find the farthest points in a blob (the feret diameters). There are a lot of ways to characterize an image, like getting distributions of areas, feret diameters, perimeters, solidities, and from those the means, medians, min, 25%, 50% and 75% points, and max of the distribution. If you care about spatial uniformity or non-uniformity or texture, there are a whole bunch of other things you can measure like distribution of inter-blob distances, Morishita plots, Fry Plots, Ripley's K function, etc. Sometimes you don't know what's important so you can just measure a whole bunch of things and get help from your statistician to tell you what correlates with some other trusted measure, like a human visual grade or the known concentration of a chemical in your sample or whatever.

[Krispy Scripts]

Thank you, thank you, thank you. I was going off previous methods used in the scientific literature, but should have come to image analysis' for analyzing images. Seriously cannot say thank you enough.

[Krispy Scripts]

I have been working on thresholding with your application. So the method I am trying to use is based off of this paper

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3294422/

It is open to the public, but could you look at figure 2 (Don't feel then the need to read any of it). I am trying to do this, most importantly cell soma size. Do you have any advice for this?

[Image Analyst]

First of all, what you did is not what they did. Secondly what they did is not that efficient and has some useless things in it.

They did not take a 120x120 pixel image around each local max. They took a 120 micron by 120 micron region. How many pixels this is depends on the magnification, but it's not going to be 120 pixels unless you have exactly 1 pixel per micron.

Then they take the region around each max, however one of their cells may have more than one local max in it (like I already said) so they may analyze all or a portion of a cell more than once. This will require you to have an accumulator binary image because you may only get part of the cell in each window. If you did get the whole cell then you may analyze it twice. Let's say your blob was a dumbbell with a bright max spot on each end. If the 120x120 field of view (FOV) encompassed the whole dumbell, you'd get the mask there twice because it does the 120x120 local region once at each end. If the 120x120 region was only big enough to get one end of the dumbbell, then one box would get one end and the other box would get the other end and you'd have to OR those two together into your accumulator image to record/save masks from both box locations.

Then they say that they do Otsu on a box and iterate it until they get a foreground area that is as close as possible to some "target cell size" that they defined in advance. Well, that's inefficient. No need to do Otsu and then iterate to "home in on" the final threshold. If you require that the foreground size (# pixels above threshold) is some number (I think they said 500 pixels or something), then you can simply use cumsum() on the histogram and find out where it first reaches or exceeds the target cell size. And you're not doing that (adjusting the original Ostu threshold guess). They say they did that, but like I said, iterating to find the final threshold (that will produce a foreground area that equals the TCS) will get you there but it's a complicated indirect way of getting the threshold, and probably slower as well.

[Krispy Scripts]

OK, thank you. They did 50 microns which converted to 120 pixels, but that is what I was confused on is finding regional max did not make sense because I thought there was most likely multiple regional max. I also did not think they were thresholding correctly. Thank you for looking at this.

Answer 2

[Walter Roberson]

http://www.mathworks.com/help/matlab/ref/ismembertol.html

[Krispy Scripts]

Would it be more appropriate and efficient to just use regionprops after thresholding to find the area of the objects?

[Walter Roberson]

That will not tell you whether you have more than one centroid in your cropped region. You are cropping relative to each centroid, and the only reason that might generate less than (288 x 288 x number of centroids) total area would be if one of the centroids happens to be less than 144 from an edge.

With your current code you are extracting a constant size around each centroid and that constant size might happen to include other centroids as well, but those centroids will also have their own centered extracted images.

If you want to crop in a way that if two centroids are closer together than 288 that the crops are merged, then that is a different task.

Answer 3

[Krispy Scripts]

Here is the attached image. I am basically trying to individually find the area, perimeter, circularity (ECD in image analyst's tutorials). So perhaps I am going about this all wrong with this code.

Any feedback or advice on how to analyze this would be so appreciated!