Count number of blobs in rows only

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.

Answer 1

[Image Analyst]

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:

https://www.olympus-lifescience.com/en/microscope-resource/primer/techniques/confocal/resolutionintro/

"... 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...."

https://www.olympus-lifescience.com/en/microscope-resource/primer/digitalimaging/deconvolution/deconresolution/

https://en.wikipedia.org/wiki/Optical_resolution#Measuring_optical_resolution

[Marcus Glover]

Thanks! That is very helpful and I think your approach makes more sense.

The goal is to automate the task of a person scoring these, and the viewer is allowed to window and level the image to the best viewing, which tends to be the smallest window. So that's whay I was binarizing and looping through various threholds. Perhaps that was important information I should have mentioned.

I feel like setting a criteria for the greyscale will be equivalent, but thought I would mention it.

The bounding box idea is a good one- that will alleviate any positioning concerns. For that I think I may have to loop through various thresholds until I arrive at the 'correct' size.

Thanks!

[Image Analyst]

I don't see how a person interactively thresholding

https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle

an image would tell you whether two points are resolvable or not.

[Marcus Glover]

Perhaps I am not using the terms correctly- the originals are DICOM and the person can adjust the window and level interactively. It is easiest to see the points when the window is very small (so the displayed image is practically binary), and then move the level until the points are distinguishable.

The rules of the test make it so four local maximum detected at once in a row will result in a 'resolvable' or 'distinguishable' pass, regardless of the local minimums near them, so that may not meet the formal definition of resolution either.

Here is how I was attempting to emulate that. You can resolve the points in the first image, but I think it's a bit more subjctibve by eye.

load im.mat

bim=imbinarize(rescale(im),.77); 

subplot(1,2,1),imshow(im,[],'InitialMagnification','fit');
subplot(1,2,2),imshow(bim,'InitialMagnification','fit');

[Image Analyst]

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]

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');