What modification is required in the given code on image processing to accurately capture the dimensions of a spreading blob from a series of images?

I have a series of images of a drop of liquid spreading on a solid surface, the google drive link for the images are attached in this query mail. The images are taken with a help of a high speed camera and I need to find the spreading length of the liquid blob. I have constructed a code which reads all the images in a folder and processes them and returns the spreading width and height. The code correctly thresholds some of the images but fails to do so for some of the image due to reflection of light on the spreading blob. What modification is required in the code so as to get all the blob captured correctly and returns the almost accurate spreading datas? link for the image folder on Google Drive

The code is given in the following :

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format short;
format compact;
fontSize = 18;

% Ask user to enter the conversion factor
%conversionFactor = input('Enter the conversion factor to convert pixel dimensions to real dimensions: ');
conversionFactor = 0.021636;

% Folder containing images
folder = 'H:\Reduced no. files for analysis\at 1.64 m_sec\on 0.25 cmc\water on 0.25 cmc\4 mm off';
% Get a list of all files in the folder
fileList = dir(fullfile(folder, '*.tif'));

% Initialize arrays to store spreading width and droplet height
spreadingWidthArray = [];
dropletHeightArray = [];

% Process each image
for i = 1:numel(fileList)
    baseFileName = fileList(i).name;
    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
        fullFileName = fullFileNameOnSearchPath;
    end
    % Read the image
    [grayImage, map] = imread(fullFileName);
    grayImage = imflatfield(grayImage, 6); % Correct uneven illumination *50 for 1 cmc on 1 cmc (1 m_sec)
    grayImage = rgb2gray(grayImage);
      
    % Display preprocessed image
    figure;
    imshow(grayImage, map);
    title('Preprocessed Image');
   
    % Thresholding
    binaryImage = ~imbinarize(grayImage);
    % Erase from line 758 down:
    binaryImage(773:end, :) = false; % 758 % 725 for 1 m/sec surface tension contrast % 759 for a little bit bright *743 for 1 on 1 (1 m/sec)
    % Fill holes
    binaryImage = imfill(binaryImage, 'holes');
    % Get rid of small noise blobs
    binaryImage = bwareafilt(binaryImage, 1); % Take largest blob only
    % Compute properties of the binary image
    props = regionprops(binaryImage, 'BoundingBox');
    spreadingWidth = props.BoundingBox(3);
    dropletHeight = props.BoundingBox(4);
    % Convert pixel dimensions to real dimensions
    spreadingWidthReal = spreadingWidth * conversionFactor;
    dropletHeightReal = dropletHeight * conversionFactor;
    % Round to 4 decimal places
    spreadingWidthReal = round(spreadingWidthReal, 4);
    dropletHeightReal = round(dropletHeightReal, 4);
    % Store spreading width and droplet height in real dimensions
    spreadingWidthArray = [spreadingWidthArray, spreadingWidthReal];
    dropletHeightArray = [dropletHeightArray, dropletHeightReal];
    % Display binary image with bounding box
    figure;
    imshow(binaryImage);
    hold on;
    rectangle('Position', props.BoundingBox, 'Edgecolor', 'g', 'LineWidth', 0.25);
    hold off;
    title('Binary Image with Bounding Box');
    impixelinfo;
    fprintf('Done running %s.m ...\n', mfilename);
end

% Display the arrays of real dimension values
disp('Spreading Widths in Real Dimensions:');
fprintf('%.4f\n', spreadingWidthArray);
disp('Droplet Heights in Real Dimensions:');
fprintf('%.4f\n', dropletHeightArray);

link for the folder of images on Google Drive

Answer 1

[Arnav]

Hi @Pragyan Kumar Sarma,

To deal with parts of the water blob that got cutoff due to reflection of light, we can use the following methods:

• For enclosing the parts of the blob that got left out of the bounding box, we can use edge detection and combine the results that we got from dark areas. These are combined using a weighted sum.
• I found having the weight of edges a bit smaller than dark areas works better. The resultant image is converted to a binary image using a threshold value of 0.5.
• Decreasing the contrast before edge detection will help leave out the edges introduced by the background noise. This is done using imadjust function.
• Other than these, we close gaps and remove small objects using imclose and bwareaopen functions.

Here is the modified code:

clc;    % Clear the command window.
close all;  % Close all figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format short;
format compact;
fontSize = 18;

% Conversion factor for pixel to real dimensions
conversionFactor = 0.021636;

% Folder containing images
folder = '.\droplets\';
fileList = dir(fullfile(folder, '*.tif'));

% Initialize arrays to store spreading width and droplet height
spreadingWidthArray = [];
dropletHeightArray = [];

% Process each image
for i = 1:numel(fileList)
    baseFileName = fileList(i).name;
    fullFileName = fullfile(folder, baseFileName);
    
    % Check if file exists.
    if ~isfile(fullFileName)
        errorMessage = sprintf('Error: %s does not exist.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
    
    % Read and preprocess the image
    rgbImage = imread(fullFileName);
    grayImage = rgb2gray(rgbImage);
    grayImage = imflatfield(grayImage, 6);
    
    gammaValue = 1.5; % Increase this value for lower contrast
    reducedContrastImage = imadjust(grayImage, [], [], gammaValue);
    enhancedImage = adapthisteq(reducedContrastImage, 'clipLimit', 0.001, 'Distribution', 'rayleigh');
    
    
    % Use morphological operations to enhance dark regions
    se = strel('disk', 12);
    tophatFiltered = imtophat(enhancedImage, se);
    bottomhatFiltered = imbothat(enhancedImage, se);
    contrastEnhanced = imsubtract(imadd(enhancedImage, tophatFiltered), bottomhatFiltered);

    % Adaptive thresholding
    binaryImage = ~imbinarize(contrastEnhanced, 'adaptive', 'ForegroundPolarity', 'dark', 'Sensitivity', 0.4);
    
    % Morphological operations
    binaryImage = bwareaopen(binaryImage, 50); % Remove small objects
    se = strel('disk', 5); % Structuring element
    binaryImage = imclose(binaryImage, se); % Close gaps
    
    % More aggressive edge detection
    edges = edge(contrastEnhanced, 'Canny', [0.9, 0.98]); % Adjusted thresholds
    
    % Apply a lower weight to edges
    edgeWeight = 0.4;  % Adjust this weight to control the influence of edges
    combinedImage = binaryImage * (1 - edgeWeight) + edges * edgeWeight;
    
    % Threshold the combined image to create a binary image
    thresholdValue = 0.5;  % Adjust this threshold as needed
    finalBinaryImage = combinedImage > thresholdValue;
    
    % Fill holes
    finalBinaryImage = imfill(finalBinaryImage, 'holes');
    
    finalBinaryImage(775:end, :) = false;

    % Keep only the largest blob
    finalBinaryImage = bwareafilt(finalBinaryImage, 1);
    
    % Compute properties of the binary image
    props = regionprops(finalBinaryImage, 'BoundingBox');
    if ~isempty(props)
        spreadingWidth = props.BoundingBox(3);
        dropletHeight = props.BoundingBox(4);
        
        % Convert pixel dimensions to real dimensions
        spreadingWidthReal = spreadingWidth * conversionFactor;
        dropletHeightReal = dropletHeight * conversionFactor;
        
        % Round to 4 decimal places
        spreadingWidthReal = round(spreadingWidthReal, 4);
        dropletHeightReal = round(dropletHeightReal, 4);
        
        % Store spreading width and droplet height
        spreadingWidthArray = [spreadingWidthArray, spreadingWidthReal];
        dropletHeightArray = [dropletHeightArray, dropletHeightReal];
    end
    
    % Display the original and processed images side by side
    figure;
    imshow(rgbImage);
    rectangle('Position', props.BoundingBox, 'EdgeColor', 'r', 'LineWidth', 0.5);
    title('Original with Bounding Box');
    fprintf('Processed %s\n', baseFileName);
end

% Display the arrays of real dimension values
disp('Spreading Widths in Real Dimensions:');
fprintf('%.4f\n', spreadingWidthArray);
disp('Droplet Heights in Real Dimensions:');
fprintf('%.4f\n', dropletHeightArray);

Outputs:

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0064.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0074.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0084.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0114.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0144.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0254.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0294.tif

Processed water on 0.25 cmc 1.64 m_sec 4 mm off0374.tif

Spreading Widths in Real Dimensions:

6.8802

9.8444

10.3420

12.8301

13.1331

12.5056

12.2676

12.1378

Droplet Heights in Real Dimensions:

2.7910

2.1852

2.7478

1.1467

1.2333

0.8654

0.7573

0.7356

You can learn more about the functions used below:

• imadjust: https://www.mathworks.com/help/releases/R2022b/images/ref/imadjust.html
• adapthisteq: https://www.mathworks.com/help/releases/R2022b/images/ref/adapthisteq.html
• bwareaopen: https://www.mathworks.com/help/releases/R2022b/images/ref/bwareaopen.html
• imclose: https://www.mathworks.com/help/releases/R2022b/images/ref/imclose.html
• edge: https://www.mathworks.com/help/releases/R2022b/images/ref/edge.html

Hope it helped you.