Count number of blobs in rows only

1 view (last 30 days)
The object I am imagng has 4 rows of 4 dots. I am interested in knowing the maximum number of blobs I can resolve in one row. I am only interested in the 4 rows of 4 in the top left of the object (I will similarly repeat for the four columns in the bottom right). The original object is much larger, and I am masking out a 24x24 portion with the dots which is always in the same relative position in the image.
I am not sure I am takling the best approach and would appreciate any insight. Here is the original object:
Here is the original image (mat file attached):
My process is to binarize the uint16 image, and then iterate through threshold levels. I then use bwconncomp to count the objects in each row, and take the maximum.
If any row in the object has all 4 objects resolved at any threshold level, the image is a scored as a 'pass', regardless of the resolution in the column direction. These 3 examples of the same image at different thresholds all pass as there are rows on the top half that have 4 objects in aat least one row in them.
Note that in the Thresh 0.48 image on the left there is a row that would be scored as having 5 components. I do apply a mask slightly larger then the top left quadrant for counting, but the images can move somewhat so this does not always work, and sometimes backfires.
I also get the situation where certain threshold levels add noise that will reult in a row having a score of 4 making an apparent pass, like this at threshold 0.05:
My question is how best to avoid these situations? Is there a better method to prepare the image for resolution assesment than my binarize/threshold? I do have some rules that help eliminate the odd ones (such as not counting any score where there are >x number of white pixels), but with a lot of images many things can happen so I am reluctant to just choose an acceptable threshold range.
I also thought of registering the original image to a template and then masking the top left quadrant for scoring, but not sure if that is worth the effort.
Thanks.

Accepted Answer

Image Analyst
Image Analyst on 15 Jan 2022
I don't think that's the correct approach at all. What I'd do is to first identify the rows and columns of the centers of the points. That can be done in a variety of ways. For example take the bounding box and just divide it by the number of expected rows and columns, and add an offset. Ther are other ways and I'm sure you can get them So, once you have the center rows and columns, then you can get the horizontal profile
horizontalProfile = grayImage(row, :);
plot(horizontalProfile, 'b-');
Now you will have peaks and valleys and you have to decide what is your criteria for two peaks being "resolved". You can use the well known Rayleigh criteria:
"... the well-known Rayleigh criterion for resolution states that two points are resolved when the first minimum (zero crossing) of one Airy disk is aligned with the central maximum of the second Airy disk. Under optimum imaging conditions, the Rayleigh criterion separation distance corresponds to a contrast value of 26.4 percent. Although any contrast value greater than zero can be specified in defining resolution, the 26-percent contrast of the Rayleigh criterion is considered reasonable in typical fluorescence microscopy applications, and is the basis for the common expression defining lateral resolution...."
  5 Comments
Image Analyst
Image Analyst on 16 Jan 2022
Edited: Image Analyst on 16 Jan 2022
If the ideal situation is to have single pixels/squares, and not rectangles/bars, then you can get the areas of all of them
props = regionprops(mask, 'Area');
allAreas = [props.Area];
numSquares = sum(allAreas == areaOfSingleSquare)
numRectangles = length(allAreas) - numSquares
You'll need to determine areaOfSingleSquare in advance. Of course these values will change for every threshold value.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
s = load('im.mat')
grayImage = s.im;
subplot(2, 2, 1);
imshow(grayImage, [], 'InitialMagnification', 1000);
axis('on', 'image')
title('Original Image')
ugl = unique(grayImage(:))
maxGL = max(ugl)
numSquares = zeros(1, maxGL+1);
numRectangles = zeros(1, maxGL+1);
areaOfSingleSquare = 1;
for k = 1 : length(ugl)
thisGrayLevel = ugl(k);
mask = grayImage >= thisGrayLevel;
subplot(2, 2, 2);
imshow(mask, 'InitialMagnification', 1000);
caption = sprintf('Thresholded from %d and brighter', thisGrayLevel);
title(caption);
drawnow;
props = regionprops(mask, 'Area');
allAreas = [props.Area];
v = allAreas == areaOfSingleSquare;
if ~isempty(v)
numSquares(thisGrayLevel+1 : end) = sum(v);
else
numSquares(thisGrayLevel+1 : end) = 0;
end
numRectangles(thisGrayLevel+1 : end) = length(allAreas) - numSquares(thisGrayLevel + 1);
end
gls = 0:maxGL;
% Plot number of squares as a function of gray level threshold.
subplot(2, 2, 3);
plot(gls, numSquares, 'b-', 'LineWidth', 2)
grid on;
xlabel('Threshold');
ylabel('Number of Squares');
title('Number of Squares');
% Plot number of rectangles as a function of gray level threshold.
subplot(2, 2, 4);
plot(gls, numRectangles, 'b-', 'LineWidth', 2)
grid on;
xlabel('Threshold');
ylabel('Number of Other Shaped Blobs');
title('Number of Other Shaped Blobs');
Marcus Glover
Marcus Glover on 16 Jan 2022
Edited: Marcus Glover on 16 Jan 2022
Thanks. Unfortunately since the resolution in the column direction does not matter for this test, I can run into cases like this where the area is large but four points in a row (and it only has to be one of the four rows to qualify) can be seen. Would work for masking each idividual row, similar to my bwconncomp idea.
The column resolution is tested using the lower right quadrant- where rowwise resolution does not matter- so this one does indeed pass both.
No test or solution is perfect of course, and you've given me some great tips. I think I will maybe try a few and see which ends up being faster. Thanks again!
load im.mat
bim=imbinarize(rescale(im),.57);
imshow(bim,'InitialMagnification','fit');

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!