Form Index Image Processing

Hello Guys,

I am struggling to get the Form Index from a bw image. The form index (FI) isdescribed as the sum of the changes in radius: and it is represented by the attached model. where R is the radius of the particle in different directions, and θ is the directional angle. In Equation 1, the form index of a circle is zero, because there are no changes in radii. I would apreciate if there is someone that can help me with it

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 = [];
baseFileName = 'eight.tif';
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
grayImage = imread(fullFileName);
% 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(grayImage)

rows =    242

columns =    308

numberOfColorChannels =      1

%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
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 = grayImage(:, :, 3);
end
% Update 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(grayImage)

rows =    242

columns =    308

numberOfColorChannels =      1

% Maximize window.
g = gcf;
g.WindowState = 'maximized';
drawnow;
%--------------------------------------------------------------------------------------------------------
% Show the histogram
subplot(2, 2, 2);
imhist(grayImage);
title('Histogram of Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
grid on;
%--------------------------------------------------------------------------------------------------------
% Interactively and visually set a threshold on a gray scale image.
% https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle
lowThreshold = 0;
highThreshold = 200;
% [lowThreshold, highThreshold] = threshold(lowThreshold, highThreshold, grayImage)
%--------------------------------------------------------------------------------------------------------
% Create a mask
mask = grayImage >= lowThreshold & grayImage <= highThreshold;
% Fill blobs
mask = imfill(mask,"holes");
% Take blobs only 1000 pixels or larger in area
mask = bwareaopen(mask, 1000);
subplot(2, 2, 3);
imshow(mask, []);
impixelinfo;
axis('on', 'image');
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Measure blob sizes
props = regionprops(mask, 'Area', 'Centroid');
allAreas = [sort([props.Area], 'descend')]

allAreas = 1×4

4786 4773 4770 4769

subplot(2, 2, 4);
histogram(allAreas);
title('Histogram of blob areas', 'FontSize', fontSize, 'Interpreter', 'None');
grid on;
% 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.
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.
    % Compute the Form Index
    xCenter = props(k).Centroid(1);
    yCenter = props(k).Centroid(2);
    % Get all the radii as we move around the boundary.
    radii = sqrt((x - xCenter).^2 + (y - yCenter).^2);
    % Get the sum of differences to get the Form Index
    FI(k) = sum(abs(diff(radii)));
end
FI % Show in command window.

FI = 1×4

72.8062472085245 71.5919083047267 69.0889051455168 65.4482079330172

hold off;
caption = sprintf('%d Outlines, from bwboundaries()', numberOfBoundaries); 
fontSize = 15;
title(caption, 'FontSize', fontSize); 
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.

[Nikolau]

Thank you so much for your help. I could see that the radii from your answer is consistent, which is good. However, it seems that the FI is not right because the FI is the variation in the figure radius and since the coins are almos circular the FI was supposed to be zero or close to zero. I think there maybe something wrong the the lFI formula that you used.

Sincerely

[Image Analyst]

Yes, it was zero but I thought it might be that for almost any shape because it always ends up at the same place after you've gone all the way around the clock. What goes up must come down. The sum of differences, if you don't take the absolute value will be close to zero. But I think the want the displacement from the prior radius so that's why I took the absolute value. That way it might be around 70 for a coin but higher, likr 150 for a gear shape. Can you provide an image where there is a smooth shape and a very irregular shape and try it both with and without the abs()?

[Nikolau]

Thank you so much for the help.

[Image Analyst]

Look, here is one with irregularly shaped blobs, some smooth and convex and some not. I compute the FI both taking abs and without taking abs.

% 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 = [];
baseFileName = 'toyobjects.png';
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
grayImage = imread(fullFileName);
% 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(grayImage)
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
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 = grayImage(:, :, 3);
end
% Update 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(grayImage)
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
drawnow;
%--------------------------------------------------------------------------------------------------------
% Show the histogram
subplot(2, 2, 2);
imhist(grayImage);
title('Histogram of Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
grid on;
%--------------------------------------------------------------------------------------------------------
% Interactively and visually set a threshold on a gray scale image.
% https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle
lowThreshold = 0;
highThreshold = 200;
% [lowThreshold, highThreshold] = threshold(lowThreshold, highThreshold, grayImage)
%--------------------------------------------------------------------------------------------------------
% Create a mask
% mask = grayImage >= lowThreshold & grayImage <= highThreshold;
mask = grayImage ~= 127;
% Fill blobs
mask = imfill(mask,"holes");
% Take blobs only 1000 pixels or larger in area
mask = bwareaopen(mask, 1000);
subplot(2, 2, 3);
imshow(mask, []);
impixelinfo;
axis('on', 'image');
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Measure blob sizes
props = regionprops(mask, 'Area', 'Centroid');
allAreas = [sort([props.Area], 'descend')]
subplot(2, 2, 4);
histogram(allAreas);
title('Histogram of blob areas', 'FontSize', fontSize, 'Interpreter', 'None');
grid on;
% 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.
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.
    % Compute the Form Index
    xCenter = props(k).Centroid(1);
    yCenter = props(k).Centroid(2);
    % Get all the radii as we move around the boundary.
    radii = sqrt((x - xCenter).^2 + (y - yCenter).^2);
    % Get the sum of differences to get the Form Index
    FI(k) = sum(diff(radii));
    FIabs(k) = sum(abs(diff(radii)));
end
FI % Show in command window.
FIabs
hold off;
caption = sprintf('%d Outlines, from bwboundaries()', numberOfBoundaries); 
fontSize = 15;
title(caption, 'FontSize', fontSize); 
axis('on', 'image'); % Make sure image is not artificially stretched because of screen's aspect ratio.

You see the FI without taking the abs are all zero as expected:

FI =

0 0 0 0

FIabs =

149.325027041666 218.54116484628 245.790271754549 84.1212568797374