How to remove black (no perfectly black) background and calculate area of cracks?

Dear all,

I'm trying to remove the background (not perfectly black) from the attached image. I have tried many methods, including setting thresholds and visualizing the data with the attached code. However, I realized that this might not be the best approach for this image because it also removes information within the circle (tomography of a cylindrical sample). I believe the best way would be to define a mask with the radius of the circle, but I don't know how to determine it. Can anyone help fix the issues in my code or define a new one for the circular mask? The final aim is to calculate the area of the cracks.

I really appreciate any help!

clear all
close all
clc
I= imread('image_1567.png');
FG = fliplr(imadjust(I,[0.05 1])); 
sout = size(I);
squaresize = [10 10];
xx = mod(0:(sout(2)-1),squaresize(2)*2)<squaresize(2);
yy = mod(0:(sout(1)-1),squaresize(1)*2)<squaresize(1);
BG = im2uint8(0.3 + bsxfun(@xor,xx,yy')*0.4);
mask = FG>10;
outpict = BG;
outpict(mask) = FG(mask);
outpict = uint8(double(FG).*mask + double(BG).*(1-mask));
subplot(3,2,1); 
imshow(I,[]);  
title('Original Grayscale Image', 'FontSize', 15);
subplot(3,2,2);
imhist(I);
title('Histogram of original image');
subplot(3,2,3);
imhist(outpict);

[counts, grayLevels] = imhist(outpict);
bar(grayLevels, counts, 'EdgeColor', 'r', 'FaceColor', 'b', 'BarWidth', 1);
xlim([0, max(outpict(:))]);
title('Histogram of output image with threshold 60');
subplot(3,2,4);
thresholdValue1 = 60;
f.InvertHardcopy = 'off';
binaryImage1 = outpict > thresholdValue1;
imshow(binaryImage1, []);
title('Binary Image threshold 60');

f.InvertHardcopy = 'off';
figure (1)

[Cris LaPierre]

For an interactive approach, try out the Image Segmenter app.

What color would you like it to be instead? Black indicates a value of 0. If you apply a mask to your image to capture the circle, everything outside the mask will be set to 0, which will appear black again.

img = imread('image_1567.png');

% Adjust data to span data range.
X = imadjust(img);

% Threshold image with global threshold
BW = imbinarize(im2gray(X));

% Create masked image.
maskedImage = X;
maskedImage(~BW) = 0;

figure
imshowpair(BW,maskedImage,'montage')

[Elisa]

dear @Cris LaPierre thank you so much for your kind response!! I would like to set NaN to the black background of the original image, and then to extract only the cracks within the circle

[Cris LaPierre]

That can be accomplished with a small change to the shared code, but it will still appear black.

img = imread('image_1567.png');

% Adjust data to span data range.
X = imadjust(img);

% Threshold image with global threshold
BW = imbinarize(im2gray(X));

% Create masked image.
maskedImage = X;
maskedImage(~BW) = nan;

figure
imshowpair(BW,maskedImage,'montage')

[Elisa]

@Cris LaPierre that for sure better than my tries!! however if you see the final result, we have lost the cracks (the dark gray to black lines in the sample). Can you suggest a way to extract them?

[Cris LaPierre]

Mask is shown on the left (0s and 1s), while the masked image is shown on the right (background removed). You would continue working with the image on the right to identify and compute properties of the cracks.

[Cris LaPierre]

You might be interested in the concepts covered in out Image Processing for Engineering and Science specialization on Coursera. Its free to enroll. One of the examples used is cracks in concrete (Automate Image Processing course)

Here's a sample video: Using the Image Batch Processor App

Answer 1

[Image Analyst]

@Elisa try this:

% Demo by Image Analyst
% 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 = 16;
%--------------------------------------------------------------------------------------------------------
%    READ IN TEST IMAGE
folder = pwd;
baseFileName = "image_1567.png";
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~isfile(fullFileName)
    % 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
% Read in image file.
rgbImage = imread(fullFileName);
% Get size
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Get gray scale version of it.
if numberOfColorChannels == 3
    grayImage = rgbImage(:, :, 2); % Take green channel.
else
    grayImage = rgbImage;
end
% Display the image.
subplot(2, 3, 1);
imshow(grayImage);
axis('on', 'image');
impixelinfo;
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
g.Name = 'Demo by Image Analyst';
g.NumberTitle = 'off';
drawnow;
%--------------------------------------------------------------------------------------------------------
% Show histogram
subplot(2, 3, 2);
imhist(grayImage);
grid on;
title('Gray Scale Histogram', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
%--------------------------------------------------------------------------------------------------------
% CREATE MASK.
% Threshold.
thresholdValue = 60;
xline(thresholdValue, 'Color', 'r');
% Create mask
circleMask = grayImage >= thresholdValue;
circleMask = bwconvhull(circleMask);
% Fill any holes.
circleMask = imfill(circleMask, 'holes');
% Take largest blob.
circleMask = bwareafilt(circleMask, 1);
% Shrink the mask a bit to avoid edge artifacts later on.
se = strel('disk', 2, 0);
circleMask = imerode(circleMask, se);
subplot(2, 3, 3);
imshow(circleMask);
axis('on', 'image');
impixelinfo;
title('Final Mask Image', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
%--------------------------------------------------------------------------------------------------------
% CREATE MASK.
sigma = 20;
Iflatfield = imflatfield(grayImage,sigma);
% Mask out background:
Iflatfield(~circleMask) = 0;
subplot(2, 3, 4);
imshow(Iflatfield) 
impixelinfo;
title('Flat Field Image', 'FontSize', fontSize, 'Interpreter', 'None');
%--------------------------------------------------------------------------------------------------------
% Show histogram
subplot(2, 3, 5);
histogram(Iflatfield(circleMask));
grid on;
title('Gray Scale Histogram of Flat Field Inside Circle', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
subplot(2, 3, 6);
crackMask = (Iflatfield < 72) & circleMask;
% Get rid of small blobs.
crackMask = bwareaopen(crackMask, 400);
imshow(crackMask)
% Label each blob with 8-connectivity, so we can color each blob individually
[labeledImage, numberOfBlobs] = bwlabel(crackMask, 8);
% Apply a variety of pseudo-colors to the regions.
coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
% Display the pseudo-colored image.
imshow(coloredLabelsImage);
%--------------------------------------------------------------------------------------------------------
% GET AREA(S)
blobMeasurements = regionprops(crackMask, 'Area');
numberOfBlobs = size(blobMeasurements, 1)
allAreas = sort([blobMeasurements.Area])
caption = sprintf('Final Mask Image with %d Blobs', numberOfBlobs);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
figure;
histogram(allAreas, 100)

Answer 2

[Image Analyst]

Defining the mask is easy. Assuming it's pure zero you can just do something like

circleMask = grayImage > 0;

It looks like the gray background inside the circle is not uniform. And you need to flatten it before you threshold to find the cracks. So you can flatten it in several ways. You can try adapthisteq or imflatfield. Hopefully that will make your histogram more bimodal, but I'm sure there will still be a judgment call as to exactly where the crack ends since they fade away into the background. One thing you could try is entropyfilt or stdfilt. These will give a high signal wherever the pixels in the mask vary a lot (like over the cracks) and give a low signal where the image is smooth (like the background). So something like this:

circleMask = grayImage > 0;
grayImage = imflatfield(.......);
grayImage = stdfilt(.......);
crackMask = grayImage < someThreshold;
crackMask = crackMask & circleMask;

If there are some edge effects out near the edge of the circle, you can call imerode to shrink the circle mask a little.

You might also look to the File Exchange for filters that find veins in the retina, like

https://www.mathworks.com/matlabcentral/fileexchange/49172-trainable-cosfire-filters-for-curvilinear-structure-delineation-in-images

or look for ridge-finding filters like Hessian in the File Exchange.

Also search this forum for cracks because it's been asked before for cracks in concrete or tiles.

By the way, you don't need to call imadjust or set the background to nan. Everything works fine without doing those things.

[Elisa]

dear @Image Analyst thank you for precious help and suggestions!!

I tried following you in that way (code below). I tried to change the threshold but as you can see not all the cracks are included. Some of them are too thin. Can you please suggest an improvement of my code? at the end, how can i calculate the area of the cracks?

thank you a lot!

clear all
close all
clc
A= imread('image_1567.png');
[rows, columns, numberOfColorChannels] = size(A);
[centers,radii] = imfindcircles(A,[500 650],Sensitivity=0.95); 
mask = circles2mask(centers,radii,size(A));
new_image = A;
new_image(~mask) = nan;
figure(1)
sigma = 20;
Iflatfield = imflatfield(new_image,sigma);
imshow(Iflatfield) 
figure(2)
imhist(Iflatfield);
figure(3)
crackMask = Iflatfield < 50;
imshow(crackMask)