Otsu Thresholding using graythresh

Hello, I have a mat file of roughly 1000 images, uint16 format, which I convert to double format to then threshold each image and subsequently find the center of mass of each image.

Currently, I have to manually look at some images and determine a suitable Threshold value myself to then input into the loop and determine the center of mass, however I need to find a way to perform an automated way of obtaining this 'Th' value, for example using the function 'graythresh' as seen commented below following the file load statement.

When I use this, however it always returns a Th value of 0 with no errors present. Please help in how to properly implement 'graythresh'

clc;
clear;
close all;
load('C:\Users\Mitre\Desktop\trappedbeaddata\new_data\2022_02_17\30pc_laser.mat');

%% Global image threshold using Otsu's method (graythresh)
Th = graythresh(ims);
ImsThresholded = double(ims).*(ims>(Th*max(ims,[],'all')));

%% Simple threshold input
% Th=input('Threshold limit: ');

%% Center Container
center=zeros(2,size(ims,3));   % zeros array

%% Image loop
for x=1:1000
    A=ims(:,:,x);           % suppress output with semicolon
    AA = double(A);         % Convert to double 
%% Size of Matrix A
[r,c]=size(AA);

%% Compare each element of matrix with threshold value
for i=1:r
    for j=1:c
        if AA(i,j)<Th
           AA(i,j)=0;        % suppress output with semicolon
        else
         %  AA(i,j)>Th;      % Remove for textured center
         %  AA(i,j)=1;       % Remove for textured center % suppress output with semicolon

            AA(i,j)=AA(i,j)-Th;   % Value - Threshold value = More complex threshold
        end
    end
end

%% Center of mass method 1: using mean
% https://uk.mathworks.com/matlabcentral/answers/363181-center-of-mass-and-total-mass-of-a-matrix
tot_mass = sum(AA(:));
[ii,jj] = ndgrid(1:size(AA,1),1:size(AA,2));
R = sum(ii(:).*AA(:))/tot_mass;
C = sum(jj(:).*AA(:))/tot_mass;

%% display / store values of x, (C,R) and FrameTime
% display 'x'
%disp("x=" + num2str(x));
% Columns = x-axis // Rows = y-axis
center(:,x) = [C;R];          

end

Answer 1

[Image Analyst]

How many blobs do you expect in the image?

You don't need to convert to double before thresholding.

To compute the center of mass for each blob you can do

props = regionprops(ImsThresholded, ims, 'WeightedCentroid');

[Gucci]

Thanks, I've already calculated the center of mass using the mean method, but at the moment I'm using the simple threshold which I have commented out, but will include below:

%% Simple threshold input
Th=input('Threshold limit: ');

This allows me to select a value from just observing some images as a rough threshold Th value

However I need to get graythresh to work to automatically obtain an accurate Th value as below: (But at the moment this just gives me a Th value = 0)

%% Global image threshold using Otsu's method (graythresh)
Th = graythresh(ims);
ImsThresholded = double(ims).*(ims>(Th*max(ims,[],'all')));

[Image Analyst]

Attach ims in a .mat file with the paperclip icon.

save('answers.mat', 'ims');

[Gucci]

ims is already in the .mat file, 30pc_laser.mat...

[Image Analyst]

I figured that. OK, Let's see if it's a standard demo file that ships with MATLAB and therefore I would have it:

load('30pc_laser.mat')

Error using load
Unable to find file or directory '30pc_laser.mat'.

OK, it's not. But good luck with it though.

[Gucci]

Sorry, forgot to mention that it is a local file I have containing 'ims' which is a 130x296x50000 uint16 struct of pixel values, and another meta data struct, however I am unable to upload it as it is much larger than 5mb even when compressed

[Image Analyst]

OK. Good luck then, unless you can crop down the data to a smaller subset of it.

[Gucci]

Hi @Image Analyst, how could I crop the 130x296x50000 struct to maybe a 130x296x3 struct, so that the file is smaller and runs through 3 iterations rather than 50000?

[Image Analyst]

Maybe

str = str(:,:,1:3);

[Gucci]

Okay thanks @Image Analyst, uploading the new 130x296x3 struct here, please see if you can get the graythresh to work on it.

[Gucci]

Here is the edited code again,

clc;
clear;
close all;
load('30pc_laser.mat');
%% Global image threshold using Otsu's method (graythresh)
Th = graythresh(ims);
ImsThresholded = double(ims).*(ims>(Th*max(ims,[],'all')));
%% Simple threshold input
% Th=input('Threshold limit: ');
%% Center Container
center=zeros(2,size(ims,3));   % zeros array
%% Image loop
for x=1:3
    A=ims(:,:,x);           
    AA = double(A);         % Convert to double 
%% Size of Matrix A
[r,c]=size(AA);
%% Compare each element of matrix with threshold value
for i=1:r
    for j=1:c
        if AA(i,j)<Th
           AA(i,j)=0;        % suppress output with semicolon
        else
         %  AA(i,j)>Th;      % Remove for textured center
         %  AA(i,j)=1;       % Remove for textured center % suppress output with semicolon
            AA(i,j)=AA(i,j)-Th;   % Value - Threshold value = More complex threshold
        end
    end
end
%% Center of mass method 1: using mean

tot_mass = sum(AA(:));
[ii,jj] = ndgrid(1:size(AA,1),1:size(AA,2));
R = sum(ii(:).*AA(:))/tot_mass;
C = sum(jj(:).*AA(:))/tot_mass;
%% display / store values of x, (C,R) and FrameTime
% display 'x'
%disp("x=" + num2str(x));
% Columns = x-axis // Rows = y-axis
center(:,x) = [C;R];          
end

[Gucci]

hi @Image Analyst, Im still a bit confused, been trying this graythresh method, but getting no progress

[Image Analyst]

You can mask your image to erase/blacken pixels below the threshold then if you don't want to use regionprops() you can do it manually like this:

grayImage = imread('moon.tif');
[rows, columns, numberOfColorChannels] = size(grayImage);
imshow(grayImage);
axis('on', 'image');
sumgx = 0;
sumgy = 0;
sumg = 0;
sumg = sum(grayImage(:))
for col = 1 : columns
    for row = 1 : rows
        gl = double(grayImage(row, col));
        sumgx = sumgx + col * gl;
        sumgy = sumgy + row * gl;
    end
end
xCOG = sumgx / sumg
yCOG = sumgy / sumg
hold on
xline(xCOG, 'LineWidth', 2, 'Color','r');
yline(yCOG, 'LineWidth', 2, 'Color','r');
caption = sprintf('Center of Gravity: (%.2f, %.2f)', xCOG, yCOG);
title(caption, 'FontSize', 18)