Create boundaries in an image based on porosity

Hi I am trying to create a boundary in my image between 3 regions of different porosities. I have an algorithm which tells me the porosity of the whole image and when I give it different parts of this image it gives me different porosity values, but I am not sure how to apply that to create the boundaries I want. Please help. Here is the binary image and roughly what I'm trying to do.

Answer 1

[Abhaya]

Hi @Nyeli Kratz,

As per my understanding, you have an algorithm which provides the porosity values given a range of an image, and you want to divide the regions in the image based on that porosity.

To achieve this you can follow the steps given below.

• Define porosity thresholds to find specific regions.
• Create masks for different porosity regions.
• Apply the boundaries to the image to distinguish different regions.

Here is a sample code where regions with different greyscale values are divided into different regions. You can use porosity values instead of pixel values to divide the regions.

% Read the grayscale image
grayImage = imread('image.jpg'); 

if size(grayImage, 3) == 3
    grayImage = rgb2gray(grayImage);  % Convert to grayscale if it's RGB
end

% Define the intensity thresholds for different intensity regions
threshold_1 = 85;  
threshold_2 = 170; 

% Create binary masks for different regions based on intensity
lowPorosityMask = grayImage < threshold_1;            % Low porosity (dark pixels)
mediumPorosityMask = grayImage >= threshold_1 & grayImage < threshold_2; % Medium porosity
highPorosityMask = grayImage >= threshold_2;          % High porosity (bright pixels)

combinedMask = mediumPorosityMask;

boundaryMask = bwperim(combinedMask);
SE = strel('line', 10,0);  % Create a linear structuring element (line shape)
boundaryMask = imdilate(boundaryMask, SE);

% Create a 3D matrix for RGB image, starting with the grayscale image
coloredImage = repmat(grayImage, [1, 1, 3]);  % Copy grayscale into the 3 RGB channels

% Set the boundary pixels to red 
coloredImage = double(coloredImage);
coloredImage(boundaryMask) = 255;  

figure;
subplot(1, 2, 1), imshow(grayImage), title('Original Grayscale Image');
subplot(1, 2, 2), imshow(coloredImage), title('Image with boundaries separating different intensity regions');

The result can be visualized as the figure shown below.

For more information, please refer to MATLAB documentation for 'bwperim' and 'imdilate' function.

• https://www.mathworks.com/help/images/ref/bwperim.html
• https://www.mathworks.com/help/images/ref/imdilate.html

[Image Analyst]

You can't define porosity as the pixel value at each point, especially since this image is practically binary. You need to look at the average value in a local region.

Here is your code using his image:

% Read the grayscale image
grayImage = imread('porosity layers.png'); 

if size(grayImage, 3) == 3
    grayImage = rgb2gray(grayImage);  % Convert to grayscale if it's RGB
end

% Define the intensity thresholds for different intensity regions
threshold_1 = 85;  
threshold_2 = 170; 

% Create binary masks for different regions based on intensity
lowPorosityMask = grayImage < threshold_1;            % Low porosity (dark pixels)
mediumPorosityMask = grayImage >= threshold_1 & grayImage < threshold_2; % Medium porosity
highPorosityMask = grayImage >= threshold_2;          % High porosity (bright pixels)
figure
subplot(3, 1, 1);
imshow(lowPorosityMask);
title('Low Porosity Mask')
subplot(3, 1, 2);
imshow(mediumPorosityMask);
title('Medium Porosity Mask')
subplot(3, 1, 3);
imshow(highPorosityMask);
title('High Porosity Mask')

combinedMask = mediumPorosityMask;

boundaryMask = bwperim(combinedMask);
SE = strel('line', 10,0);  % Create a linear structuring element (line shape)
boundaryMask = imdilate(boundaryMask, SE);

% Create a 3D matrix for RGB image, starting with the grayscale image
coloredImage = repmat(grayImage, [1, 1, 3]);  % Copy grayscale into the 3 RGB channels

% Set the boundary pixels to red 
coloredImage = double(coloredImage);
coloredImage(boundaryMask) = 255;  

figure;
subplot(2, 1, 1), imshow(grayImage), title('Original Grayscale Image');
subplot(2, 1, 2), imshow(coloredImage), title('Image with boundaries separating different intensity regions');

Answer 2

[Image Analyst]

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 = "porosity layers.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.
grayImage = imread(fullFileName);
% Get size
[rows, columns, numberOfColorChannels] = size(grayImage)

rows =    647

columns =         1159

numberOfColorChannels =      3

% Get gray scale version of it.
if numberOfColorChannels == 3
    grayImage = grayImage(:, :, 1); % Take red channel.
end
% Display the image.
subplot(3, 2, 1);
imshow(grayImage);
axis('on', 'image');
impixelinfo;
caption = sprintf('Original image : %s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Maximize window.
g = gcf;
g.WindowState = 'maximized';
g.Name = 'Demo by Image Analyst';
g.NumberTitle = 'off';
drawnow;

%--------------------------------------------------------------------------------------------------------
% FIRST FIND THE FAR RIGHT EDGE BETWEEN THE TISSUE AND THE AIR
air = grayImage == 255;
% Extract only the blob touching the right border
mask = imclearborder(air, 4, "Borders","right");
air(mask) = false;
% Fill holes and take the largest blob.
air = bwareafilt(imfill(air, 'holes'), 1);
% Display the masked image.
subplot(3, 2, 2);
imshow(air, []);
axis('on', 'image');
impixelinfo;
title('Air Mask Image', 'FontSize', fontSize);
drawnow;
% Get boundary of the
airBoundaries = bwboundaries(air);
hold on;
visboundaries(airBoundaries, 'LineWidth',2);

%--------------------------------------------------------------------------------------------------------
% GET THE AVERAGE VALUE IN A BOX.  THIS WILL BE THE LOCAL POROSITY.
% Experiment around with different window sizes until you get what you like.
kernel = ones(71, 25); % A window of 71 rows and 25 columns.
kernel = kernel / sum(kernel, "all");
% localSumImage = conv2(grayImage, kernel, "same");
localSumImage = imfilter(grayImage, kernel, 'replicate');
% Display the masked image.
subplot(3, 2, 2);
imshow(localSumImage, []);
axis('on', 'image');
impixelinfo;
title('Local Porosity Image', 'FontSize', fontSize);
drawnow;

%--------------------------------------------------------------------------------------------------------
% GET AVERAGE VERTICAL PROFILE.  (NOT USED - JUST FOR YOUR INFO)
verticalProfile = sum(localSumImage, 1);
subplot(3, 2, 3);
plot(verticalProfile, 'b-', 'LineWidth', 2);
grid on;
xlabel('Column Number', 'FontSize',fontSize)
ylabel('Sum within column', 'FontSize',fontSize)
title('Average Vertical Profile', 'FontSize', fontSize);
drawnow;

%--------------------------------------------------------------------------------------------------------
% Use imsegkmeans to divide it into regions.
numberOfClasses = 3; % Adjust as necessary
[labeledImage, Centers] = imsegkmeans(localSumImage, numberOfClasses);
B = labeloverlay(localSumImage, labeledImage);
subplot(3, 2, 4);
imshow(B)
axis('on', 'image');
title("Initial K-Means Labeled Image", 'FontSize', fontSize)

%--------------------------------------------------------------------------------------------------------
% Take the largest blob for each label and fill in islands of other labels -- fill holes.
finalClasses = zeros(size(labeledImage), class(labeledImage)); % Make output classification image.
for k = 1 : numberOfClasses
    thisClass = ismember(labeledImage, k);
    % Take largest blob
    thisClass = bwareafilt(thisClass, 1);
    
    % There are some holes at the top and bottom that won't get filled by
    % imfill because they are not completely enclosed.  We can completely
    % enclose those by drawing a line along the top and bottom and then
    % calling imfill.
    % Find first and last pixel in top line of this class
    col1 = find(thisClass(1, :), 1, 'first');
    col2 = find(thisClass(1, :), 1, 'last');
    if ~isempty(col1) && ~isempty(col2)
        thisClass(1, col1:col2) = true; % Write a white line across.
    end
    % Find first and last pixel in last line of this class
    col1 = find(thisClass(end, :), 1, 'first');
    col2 = find(thisClass(end, :), 1, 'last');
    if ~isempty(col1) && ~isempty(col2)
        thisClass(end, col1:col2) = true; % Write a white line across.
    end
    % Fill holes
    thisClass = imfill(thisClass, 'holes');
    finalClasses(thisClass) = k;
    % Show this class in the lower right of the figure
    imagePosition = min(3*6, k+15);
    subplot(3, 6, imagePosition)
    imshow(thisClass);
    axis('on', 'image');
    caption = sprintf('Class #%d', k);
    title(caption, 'FontSize', fontSize);
end

% Right now the class numbers are assigned randomly so that the left most
% porisity zone might be class 1 for one image and class 2 for some other
% image.  Let's sort them so that the zones are labeled from 1 to 3 going
% from left to right, so that the leftmost label is always labeled as 1.
% Get the centroids of each label
props = regionprops(finalClasses, 'Centroid');
xy = vertcat(props.Centroid);
xCentroid = xy(:, 1)

xCentroid = 3×1

438.73480887561 835.390401814233 147.547561186356

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Get the initial order of the classes.
[sortedValues, initialClassOrder] = sort(xCentroid, 'ascend')

sortedValues = 3×1

147.547561186356 438.73480887561 835.390401814233

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

initialClassOrder = 3×1

3 1 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Make a new labeled image where they are sorted along the x axis.
finalClasses2 = zeros(size(finalClasses), class(finalClasses));
for k = 1 : numberOfClasses
    % Extract original class
    thisClass = ismember(finalClasses, initialClassOrder(k));
    % Reassign its label number to be it's x position order.
    fprintf('Assigning original class %d (at x = %.1f) to new class label number %d.\n', k, xCentroid(k), initialClassOrder(k))
    finalClasses2(thisClass) = k;
end

Assigning original class 1 (at x = 438.7) to new class label number 3.
Assigning original class 2 (at x = 835.4) to new class label number 1.
Assigning original class 3 (at x = 147.5) to new class label number 2.

% Put sorted labeled image back into finalClasses.
finalClasses = finalClasses2;

% To verify they were sorted correctly, let's display then in class order
for k = 1 : numberOfClasses
    thisClass = ismember(finalClasses, k);
    % Show this class in the lower right of the figure
    imagePosition = min(3*6, k+15);
    subplot(3, 6, imagePosition)
    imshow(thisClass);
    axis('on', 'image');
    caption = sprintf('Class #%d', k);
    title(caption, 'FontSize', fontSize);
end

●

% Display the final classes, sorted by x.
h5 = subplot(3, 2, 5);
imshow(finalClasses, [])
colormap(h5, 'turbo')
axis('on', 'image');
title("Final Labeled Image", 'FontSize', fontSize)

% If you want, we can also get the boundary of each class
% Here is where we actually get the boundaries for each blob.
for k = 1 : numberOfClasses
    % Extract this class number.
    thisClass = ismember(finalClasses2, k);
    % Get its boundaries.
    boundaries = bwboundaries(mask);
    % Extract and save the boundary for this class.  There will be only 1 boundary.
    classBoundaries{k} = boundaries{k};
end

So now you have a binary image mask for each porosity zone as well as the boundaries of each zone. You shouldn't really need both but anyway with them you should now be able to do whatever you want to do, such as getting the mean porosity of each zone or whatever.