Image processing : can we remove objects attached to others ?

Ok I finally ask a question because I can't figure it out despite my researchs.

I want to process footprints photos to calculate the Arch Index (Cavanagh and al. 1987) :

To calculate the arch index we use the easy equation above, however we have to remove toes before calculation.

Therefore, I took a photo of my footprint ('testPab.png') , there is two footprints on paper (left and right foot, right foot has less ink) then I tried a way to process the data to get a binary image with footprints.

data = imread("testPab.png");         

% Trying to smooth data
H=fspecial('gaussian',30,10);
AdjSmooth=imfilter(data,H,'replicate');
grayIMG=im2gray(AdjSmooth);

% improve contrast
edgeThreshold = 0.8;
amount = 0.4;
AdjGrayIMG = localcontrast(grayIMG, edgeThreshold, amount);
contrastIM = imadjust(AdjGrayIMG);

% Closing the image
SE=strel('rectangle',[10000 10000]); %[height width]
BW = imbothat(contrastIM,SE);


BW=imlocalbrighten(BW,0.6);
filteredIM = medfilt2(BW,[100 100]);

bwData = imbinarize(filteredIM,"adaptive","ForegroundPolarity","bright",'Sensitivity',0.6);

bigObject = bwareafilt(bwData,2);   % keep the 2 biggest forms
binaryFoot = imfill(bigObject,'holes'); % fill the holes
montage({data,binaryFoot})

How can I detect and remove the toes? Is ther a way to interpolate the form without the toes ?

I thought about finding circles (imfindcircles) or the number of [0 1] I have for each rows to cut the toes but I still can't figure out the good way to do it.

PS: On R2022a I get this result with the same code, so I have to improve robustness if you have some advices :D

Answer 1

[Image Analyst]

Here's a start. Adapt as needed

% 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 = 'testPab.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 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;
%--------------------------------------------------------------------------------------------------------
% Display the histogram.
subplot(2, 2, 2);
histogram(grayImage(grayImage>0), 256);
grid on;
title('Histogram of Image Gray Levels', 'FontSize', fontSize, 'Interpreter', 'None');
% Threshold the image to get the dark stuff.
% https://www.mathworks.com/matlabcentral/fileexchange/29372-thresholding-an-image?s_tid=srchtitle
lowThreshold = 0;
highThreshold = 201;
% [lowThreshold, highThreshold, lastThresholdedBand] = threshold(lowThreshold, highThreshold, grayImage)
mask = grayImage >= lowThreshold & grayImage <= highThreshold;
% Do an opening to break lines.
windowSize = 11;
se = zeros(windowSize,windowSize);
se(ceil(windowSize/2), :) = 1;
se(:, ceil(windowSize/2)) = 1;
% se = ones(windowSize);
% mask = imopen(mask, se);
mask = imerode(mask, se);
windowSize = 15;
se = zeros(windowSize,windowSize);
se(ceil(windowSize/2), :) = 1;
se(:, ceil(windowSize/2)) = 1;
mask = imdilate(mask, se);
% Fill in the mask.
mask = imfill(mask, 'holes');
% Smooth it some more.
windowSize = 21;
kernel = ones(windowSize)/ windowSize^2;
mask = conv2(double(mask), kernel, 'same');
mask = mask > 0.5;
% Take 2 largest blobs
mask = bwareafilt(mask, 2);
subplot(2, 2, 3);
imshow(mask)
title('Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
impixelinfo
drawnow;
subplot(2, 2, 4);
imshow(grayImage)
title('Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
% 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.
% 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.

[Pablo Rozier-Delgado]

Thank you very much for your help! I'm looking forward to implement this to my original code.

So if i'm correct, I did filter too much and your solution resides in the fact that I have to filter less to keep the toes quite separated from the footprint and then keep the two biggest forms.

[Image Analyst]

It's kind of tricky because some of the lighter stuff on the left of the right foot you might want to include but it's the same darkness as the region connecting the toes to the ball, so if you raise the threshold to include that you'll also connect the toes. If you have enough images, like over a hundred, then you might try using SegNet to automatically mask the region.

[Pablo Rozier-Delgado]

Yes I have hundreds of footprints. Thank you I'm going to see the SegNet thing

Answer 2

[Pablo Rozier-Delgado]

After @Image Analyst comments and suggestions, for those who are interested, here it is the final way I did it.

The erosion was really a huge part of what I was missing.

I also added the watershed transform to remove toes, because on many others pictures the preprocessing wasn't enough to remove them.

Finally we can calculate the Arch Index :D

It's almost working on all footprints, but I still have some fails when there is not that much ink, or when picture luminosity is poor.