Sludge settling (SV45) level measurement using image processing or computer vision

I am new to MATLAB. I am trying to use an external cam to monitor a graduated cylinder (0.36 m high with 1000 mL capacity) from its front for sludge settling (SV45). I am also trying to make it plot the height of sludge as it settles over time (m vs min). The time intervals are 0, 0.5, 1, 2, 3, 4, 5, 10, 15, 20, 30, 45.

Any help on this would be much appreciated. Happy holidays to everyone.

This is currently what I am working on:

1. Load Sample Image (Serves as background reference)
This example uses the reference image example.png. 
refImage = imread("example.png");
This image is of an empty glass with a resolution of 640 x 480. 

Alternatively, you can capture this image using the Acquire Webcam Image task.

2. Set up Output Window (Plot GUI)
Set up the figure to plot the sludge level changes in real time.  
fig = figure("Name","Sludge level tracker","Visible","off");
t = tiledlayout(2,2,"Parent",fig);
t.Padding = 'none';
t.TileSpacing = 'none';
ax1 = nexttile(t);
hold (ax1,'on');
title("Reference Image","Parent",ax1);
imshow(refImage,"Parent",ax1); 
ax2 = nexttile(t,[2 1]);
title('Sludge level vs Time',"Parent",ax2);
xlabel('Time (in s)',"Parent",ax2);
ylabel('Sludge level (normalized)',"Parent",ax2);
ax3 = nexttile(t);
title("Sludge Level","Parent",ax3);
grid on;
h = animatedline(ax2,"Color","red","LineWidth",3);
webcamLETFig = figure("Name","Acquire Webcam Image Output");

3. Plot Sludge Level Change 
To determine the change in sludge level, use image subtraction. Subtract each sample image from the reference image. To capture the sample images,use the Acquire Webcam Image task.
Note: 
Ensure the following in the Live Task - 
Select a webcam from the webcams connected to your computer. 
Clear the Show image check box to prevent displaying the image in each loop iteration.
Turn off the Autorun feature if you want the live task to execute only when the entire script is executed. To disable autorun, click the autorun   icon.
Select Live preview to view the live video from the webcam while adjusting webcam properties.
Do the following - 
Capture images using the Acquire Webcam Image task in a loop. You can vary the time for which the loop executes by entering the time in seconds in the textbox. 
Subtract each sampleImage from the refImage to get the diffImage. 
Extract the the RGB components from diffImage and convert the image to grayscale to reduce the noise in the image.
Find the mean pixel intensity value across each row of the image. The mean pixel intensity is used to find the image pixel corresponding to the top of the sludge level.
Define a window size and a threshold for the intensity by moving the slider. 
Find the difference in the instensities between the n+10th pixel and nth pixel in the windowSize. 'n' is the row index and 10 is the window size used in this example. If the intensity of the pixel value is greater than the threshold intensity, store the pixel value as the sludge level.  
Plot the sludge level against time.

% Start time-stamp timer for plotting sludge level against time
t = tic;

while(toc(t) < 2700)
   % Read the image captured from the webcam
  
   % Make the "Acquire Webcam Image Output" figure the active figure
  figure(webcamLETFig);
   

% Connect to selected webcam
cam = webcam("TOSHIBA Web Camera - FHD",...
	"Resolution","640x480",...
	"Brightness",50,...
	"Contrast",32,...
	"Saturation",64,...
	"Hue",0,...
	"Gamma",300,...
	"Sharpness",50,...
	"WhiteBalanceMode","auto",...
	"BacklightCompensation",0);

% Acquire webcam image
img = snapshot(cam);

% Show acquired image
imshow(img);

% Clear existing webcam connection
clear cam
    sampleImage = img;
    
    % Store the relative timestamp
    tStamp = toc(t);
    
    % Subtract the sample image from the reference image
    % 'imsubtract' from the Image Processing Toolbox can be used here
    diffImage = refImage - sampleImage;
    dims = size(diffImage);
    
    % Extract RGB components from the difference image
    R = diffImage(:,:,1);
    G = diffImage(:,:,2);
    B = diffImage(:,:,3);
    
    % Convert the difference image to grayscale
    greyDiffImage = zeros(dims(1),dims(2),'uint8');
    for row = 1:dims(1)
        for col = 1:dims(2)
            greyDiffImage(row,col) = (0.2989 * R(row,col)) + (0.5870 * G(row,col)) + (0.1140 * B(row,col));
        end
    end
    
    % Find mean pixel intensity value across each row of the image
    rowMean = zeros(1,dims(1));
    for row = 1: dims(1)
        rowMean(row) = mean(greyDiffImage(row,:));
    end
    
    % Find the image pixel corresponding to the top of sludge level
    % Find the difference instesities between the n+10th and nth pixel.
    % Store this pixel value as the sludge level, if this intesnity is greater than a threshold
    windowSize = 10;
    intensityThreshold = 15;
    sludgeLevelPixel = dims(1);
    for rowIndex = 1:length(rowMean)-windowSize
        if rowMean(rowIndex+windowSize) - rowMean(rowIndex) > intensityThreshold
            sludgeLevelPixel = rowIndex;
            break;
        end
    end
        
    % Display images
    
    % Make the "Sludge Level Tracker" figure the active figure
    figure(fig);
    
    % Display the reference image
    
    % Display the current sludge level
    hold(ax3, 'all');  
    imshow(sampleImage,"Parent",ax3); 
    
    quiver(dims(2)/2,sludgeLevelPixel,0,dims(1)-sludgeLevelPixel,...
        "Color","yellow","LineWidth",3,'Parent',ax3);
    
    hold off;
    
    % Plot sludge level vs time
    addpoints(h,tStamp,abs(1 - (sludgeLevelPixel/dims(1))));
    drawnow;
    
    fig.Visible = "on"';
end

Answer 1

[Image Analyst]

Try this. Adapt thresholds 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 = 14;
markerSize = 40;
%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = [];
baseFileName = 'cylinder dirty.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)
if numberOfColorChannels > 1
    % It's not really gray scale like we expected - it's color.
    % Extract the blue channel.
    grayImage = grayImage(:, :, 3);
end
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 2, 1);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
title('Original Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
hold on
% Maximize window.
g = gcf;
g.WindowState = 'maximized'
drawnow;
%--------------------------------------------------------------------------------------------------------
% Get a horizontal profile
horizontalProfile = mean(grayImage, 1);
% Plot it
subplot(2, 2, 3);
plot(1:columns, horizontalProfile, 'b-')
grid on;
title('Horizontal Profile', 'FontSize',fontSize)
xlabel('columns', 'FontSize',fontSize)
ylabel('Intensity', 'FontSize',fontSize)
xlim([1, columns])
% Find out where it dips below 190 so we can avoid edges of glass.
margin = 10; % Pixels inward, so we can avoid glass edge artifact.
col1 = find(horizontalProfile < 190, 1, "first") + margin;
col2 = find(horizontalProfile < 190, 1, "last") - margin;
% Plot edges on plot.
xline(col1, 'Color', 'r', 'LineWidth', 2);
xline(col2, 'Color', 'r', 'LineWidth', 2);
% Plot edges on image.
subplot(2, 2, 1);
xline(col1, 'Color', 'r', 'LineWidth', 2);
xline(col2, 'Color', 'r', 'LineWidth', 2);
%--------------------------------------------------------------------------------------------------------
% Get a vertical profile
verticalProfile = mean(grayImage(:, col1:col2), 2);
% It has spikes, like at the top of the cylinder so get rid of those.
verticalProfile = movmean(verticalProfile, 19);
% Plot it
subplot(2, 2, 2);
plot(verticalProfile, 1:rows, 'b-')
grid on;
axis ij % Flip so that origin is at top like the image.
title('Vertical Profile', 'FontSize',fontSize)
xlabel('columns', 'FontSize',fontSize)
ylim([1, rows])
%--------------------------------------------------------------------------------------------------------
% Find the height of the dirty fluid.
thresholdValue = 160
row1 = find(verticalProfile < thresholdValue, 1, "first");
row2 = find(verticalProfile < thresholdValue, 1, "last");
% Display original image.
subplot(2, 2, 4);
imshow(grayImage);
impixelinfo;
axis('on', 'image');
yline(row1, 'Color', 'g', 'LineWidth', 2);
yline(row2, 'Color', 'g', 'LineWidth', 2);
drawnow;
fluidHeight = row2 - row1
caption = sprintf('Fluid Height = %d pixels', fluidHeight);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

Multiply the fluidHeight by the calibration factor in mm per pixel to get the height in mm.

[Lime Gatnil]

Thank you so much for this Doctor. Could you clarify what the calibration factor is?

[Image Analyst]

Try this:

img = imshow('cylinder empty.jpeg')
uiwait(helpdlg('Draw a line fropm 1000 to 100, adjust endpoints, then hit enter'))
hLine = drawline("Color",'r')
position = hLine.Position
deltax = position(1,1) - position(2,1)
deltay = position(1,2) - position(2,2)
distanceInPixels = abs(sqrt(deltax^2 + deltay^2))
% Define the distance in milliliters.
distanceInML = abs(1000 - 100)
mlPerPixel = distanceInML / distanceInPixels
% Yay!!!  You're now calibrated!

% Now to get any volume in ml, multiply the distance
% in pixels by mlPerPixel to get ml.
% For example
distanceInPixels = 321; % Whatever
volumeInMl = distanceInPixels * mlPerPixel

[Lime Gatnil]

Thank you again Doctor. You have been a big help!

Answer 2

[Image Analyst]

I have a program to identify suds level in a cylinder. But I can't try it on your image because you forgot to post your image.

[Lime Gatnil]

Good day Doctor. Would this be okay?

[Image Analyst]

Yes, for calibration. But we'd also need one with a colored liquid in it so that we can measure its height.

[Lime Gatnil]

Hello again Doctor. Would it be possible to use a different image of graduated cylinder with the colored liquid?

[Lime Gatnil]

[Lime Gatnil]