Thresholding of images gone wrong

Hello,

I have attached a small file of 3 images...

When I run the following code of the function 'graythresh', I end up with a strange binary image that also changes part of the background to 1, I am not sure how to use graythresh to select the threshold level efficiently...

%%
clc;
clear;
close all;
load('doubleBeadSample');

A=doubleBeadSample(:,:,1);           

%% Size of Matrix A
[r,c]=size(A);

%% Otsu threshold
B=im2gray(A);
%figure;
 temp=graythresh(B);
 bin_im=im2bw(B,temp);

%%
figure; 
subplot(1,2,1)
imagesc(doubleBeadSample(:,:,1)); %
title(' frame');
subplot(1,2,2)
imagesc(bin_im); %  image after thresholding
title(' frame with threshold');

[Gucci]

@Image Analyst @Image.Processor is it possible to see where I went wrong here? This is the output image I am getting on the right..

Answer 1

[Image Analyst]

Try this:

% Demo by Image Analyst
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 = 22;
markerSize = 40;
%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'doubleBeadSample.mat';
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
s = load(fullFileName)
rgbImage = s.doubleBeadSample;
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
if numberOfColorChannels > 1
    % It's not really gray scale like we expected - it's color.
    fprintf('It is not really gray scale like we expected - it is color\n');
    % Extract the blue channel.
    grayImage = rgbImage(:, :, 3);
else
    grayImage = rgbImage;
end
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
drawnow;
% Display histogram.
subplot(2, 2, 2);
imhist(grayImage);
grid on;
drawnow;
title('Histogram of Image', 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Binarize the image to get a mask.
lowThreshold = 4000;
highThreshold = intmax("int16")
% threshold(lowThreshold, highThreshold, grayImage);
mask = grayImage >= lowThreshold & grayImage <= highThreshold;
% Put red threshold line on histogram so they know where it was thresholded at.
xline(lowThreshold, 'Color', 'r', 'LineWidth', 2)
% Delete blobs touching border.
mask = imclearborder(mask);
mask = imfill(mask, 'holes');
% Take at most 2 blobs.
mask = bwareafilt(mask, 2);
% Display mask image.
subplot(2, 2, 3);
imshow(mask);
hold on;
impixelinfo;
axis('on', 'image');
drawnow;
title('Mask, a Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Plot the borders of all the blobs in the overlay above the original grayscale image 
% using the coordinates returned by bwboundaries().
% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
subplot(2,2,4);
imshow(grayImage, []); % Optional : show the original image again.  Or you can leave the binary image showing if you want.
% Here is where we actually get the boundaries for each blob.
boundaries = bwboundaries(mask);
% boundaries is a cell array - one cell for each blob.
% In each cell is an N-by-2 list of coordinates in a (row, column) format.  Note: NOT (x,y).
% Column 1 is rows, or y.    Column 2 is columns, or x.
numberOfBoundaries = size(boundaries, 1); % Count the boundaries so we can use it in our for loop
% Here is where we actually plot the boundaries of each blob in the overlay.
hold on; % Don't let boundaries blow away the displayed image.
for k = 1 : numberOfBoundaries
    thisBoundary = boundaries{k}; % Get boundary for this specific blob.
    x = thisBoundary(:,2); % Column 2 is the columns, which is x.
    y = thisBoundary(:,1); % Column 1 is the rows, which is y.
    plot(x, y, 'r-', 'LineWidth', 2); % Plot boundary in red.
end
hold off;
caption = sprintf('%d Outlines, from bwboundaries()', numberOfBoundaries); 
title(caption, 'FontSize', fontSize); 
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.
% Get all the blob properties.  Can only pass in originalImage in version R2008a and later.
props = regionprops(mask, grayImage, 'Area', 'Centroid', 'WeightedCentroid');
allAreas = [props.Area]
centroids = vertcat(props.Centroid)
weightedCentroid = vertcat(props.WeightedCentroid)
% Plot centroids with a red dot
hold on;
plot(centroids(:, 1), centroids(:, 2), 'r.', 'MarkerSize', 50);
% Plot a yellow line between them.
plot(centroids(:, 1), centroids(:, 2), 'y-', 'LineWidth', 4);
hold off;
deltas = diff(centroids, 1)
distance = sqrt(deltas(1) .^ 2 + deltas(2) .^ 2)
% Tell the user
message = sprintf('The distance between centroids is %.1f pixels', distance);
title(message, 'FontSize', fontSize); 
uiwait(helpdlg(message))

[Gucci]

Thanks @Image Analyst, this works, i want to ask you what type of automatic threshold is this? Is this a cluster based method? Also I see you have set the low threshold to 4000, however is there a way to automate this? As in to let the algorithm determine this value with graythresh or otsuthresh?

And how can I add a loop around it for multiple set of images?

[Image Analyst]

It was a manual, fixed threshold which is best if you have control over your image capture situation like exposure time and lighting level. If you don't have control over your image capture conditions, and the histogram shifts around from image to image, you'll have to have an automatic, variable threshold but this is not ideal.

Since the threshold I picked seems to be at the valley between two peaks, you can negate the histogram counts and find the first peak, which would be the first valley between the uninverted histogram. You can use findpeaks for that. You might want to smooth the counts with sgolayfilt() first. You can probably do it but if you just can't figure it out, you can come back for help. But again, it's best if you don't have to do this and can just use a fixed, permanent threshold. Either way, the centroid locations, and thus the distance between them, shouldn't vary much with the threshold since it's a pretty symmetric circle we're looking at.

[Gucci]

Okay thank you @Image Analyst, I understand this now.

I am now trying to loop this around roughly 10000 images, and will then get a matrice of all the centroid values.

I had one more question, I am trying to compare 2 particle tracking techniques, by finding the displacement between the two beads, using firstly this technique (Regionprops center of mass) and now I am attempting to do Cross-correlation. Would you have any experience in this?

I have gotten so far as to cross correlate the first 2 images, and obtain a matrix, however how can I obtain the centroid between the two peaks on either side of the matrix, while negating the one peak in the middle? i fail to understand what the peak in the middle is: here is my code and my result..

I understand I use the centroid method right after the cross correlation technique on the new matrix 'X', however is there another way to determine the displacement between them? (Of course after using cross correlation)

Thank you...

clc;
clear;
close all;
%%
load('doubleBeadSample.mat');
    
%% Center Container
center=zeros(2,2*size(ims,3));   % preallocation

%%
%for k = 2 : size(ims, 3)
 %   template = ims(:, :, k-1); % Template is prior image.
 %   full = ims(:, :, k); % This is the "current" image.

template=ims(:,:,1);
full=ims(:,:,2);

S_temp = size(template);
S_full = size(full);
%%
X=normxcorr2(template, full);
%r=max(X(:));
%[i,j]=find(X==r);

%% Otsu Thresholding of cross-correlated matrix
B=im2gray(X);
%figure;
 temp=graythresh(B);
 bin_im=im2bw(B,temp);
 %imshow(bin_im);
 %title('Thresholded image of cross-correlation matrix');

%% centroid of cross-correlation matrix
%tot_mass = sum(bin_im(:));
%[ii,jj] = ndgrid(1:size(bin_im,1),1:size(bin_im,2));
%R(k) = sum(ii(:).*bin_im(:))/tot_mass;
%C(k) = sum(jj(:).*bin_im(:))/tot_mass;


props = regionprops(bin_im, B, 'WeightedCentroid', 'Centroid');
% Extract from  structure into matrices
% xy = vertcat(props.Centroid); % Centroid of binary blobs.
xyWeighted = vertcat(props.WeightedCentroid); % Centroid of binary blobs but weighted by gray level.

%center(:,k) = [xyWeighted(1);xyWeighted(2)];
%end


%% figures
figure;
subplot(2,2,2), imagesc(full), title('Subsequent image');
subplot(2,2,1), imagesc(template), title('Image used as template');
subplot(2,2,3), imagesc(X), title('Cross-Correlation matrix');

R = zeros(S_temp);
shift_a = [0 0];
shift_b = [i j] - S_temp;
%R((1:S_full(1))+shift_a(1), (1:S_full(2))+shift_a(2)) = full;
%R((1:S_temp(1))+shift_b(1), (1:S_temp(2))+shift_b(2)) = template;
subplot(2,2,4), imagesc(bin_im), title('Otsu-Binary mask on Cross correlation matrix');

figure;
mesh(X), title('Mesh graph of cross-correlation');

This is my figure results I get...

[Image Analyst]

I'm attaching a normxcorr2() demo, but I don't see how it's applicable here. Since you have two similar regions you will get 3 peaks in the cross correlation (as you've shown). But you'd then just have to use regionprops() still, same as before, but now you'd do it on the cross correlation image. Then you'd have two sets of blobs to compute the centroids and distance for instead of 1 set. And you'd get two shifts (distances) which should be almost identical to what you got using my method. So what's the point? How is that better?

[Gucci]

I understand the reasoning, and am sure the two algorithms will be almost identical, but I am trying to get the standard deviation of the distance and compare the two different methods.

Whichever will attain the lowest Std will be the most optimal algorithm in further experiments.

[Gucci]

@Image Analyst I am a bit confused in my correlation matrix, I understand there will be 3 peaks, however I am not sure why there are 3 instead of 2?

And how is the central peak calculated? Can you shed some light on this if possible

[Image Analyst]

Let's say the blobs are called A and B. And say that the fixed array is on the top and the moving array is below. So when the shift is most of the way to the left you'll have this

A B

A B

So that will give the first peak. Now when the shift is 0 the moving image aligns with the fixed image

A B

A B

So that gives you the second (middle) peak.

Now when the shift is most of the way to the right, you'll have this:

A B

A B

That will give you the third peak.

[Gucci]

Thank you very much for this explanation @Image Analyst

To clarify, my initial matrix size for both images is [m X n], Now my cross correlation matrix size is

[(m*2)-1 X (n*2)-1].

Let us say centroid of Blob 1 is A (left particle), centroid of Blob 2 is B (Right particle) and centroid of blob 3 is C, (Central peak position):

How can I now find the true position of Blob A and Blob B in all matrices?

I took mean(A,C) to get the position of the Left particle, and mean(B,C) to get the position of the Right particle, but after your explanation I fear this is incorrect.

[Image Analyst]

@Gucci, using correlation you don't get the true position. You can do that like I showed you at first. If you want to get it, you have to indirectly get it with some complicated bookkeeping. It might help to look at a simplified example where you have known parameters, like this example. There is a known peak at location 4 and 7 in a vector of 11 elements. The first peak is at 4 and the second is at 7 and they are separated by 3. If it's the separation you want through this complicated way of arriving at it, see the example below:

v = [0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0]

v = 1×11

0 0 0 1 0 0 1 0 0 0 0

trueLocations = find(v)

trueLocations = 1×2

4 7

c = xcorr(v, v);
% Get rid of numbers smaller than 1e-10
c(c<1e-10) = 0

c = 1×21

0 0 0 0 0 0 0 1 0 0 2 0 0 1 0 0 0 0 0 0 0

% Find the indexes of the peaks
peakIndexes = find(c>0)

peakIndexes = 1×3

8 11 14

meanSeparation = mean(diff(peakIndexes))

meanSeparation = 3

[Gucci]

@Image Analyst Thank you, I understand your method above for simple matrices to get the separation,

however the complicated bookeeping that you said is unfortunately my task with this, to perform the same operation as the normal code so that I can measure the error between the 2..

I have tried but keep failing to convert the values I get for the 3 blobs, into 2 positions and then find the distance between object 1 and object 2, which should be equivalent to that of the normal code..

[Image Analyst]

I showed you how to do it with a simple 1-D example, and I showed you how to do it with images (attached again) so you should be able to figure it out. I will be traveling for the next 4 days and probably won't be able to answer any but the simplest of questions from people (not so easy to do code on my iPad which is the only computer I'm taking with me). Plus I'll be busy all day. But it's not hard and I'm sure by Saturday you'll have it all figured out.

Answer 2

[Image Analyst]

What parts of the image would you want to be foreground and background?

I think it's best to look at the image in grayscale.

See my interactive threshold in my File Exchange:

Image Analyst's File Exchange