Improving Droplet Shape Detection and Tracking from Video

I am trying to write code to track a water droplet over time. This includes identifying the center of the droplet and measuring changes in its circumference. I then calculate the droplet’s displacement in the x and y directions, perform FFT analysis, and plot the 2D trajectory of the droplet's movement. Finally, I generate a 3D volume video showing how the droplet’s shape changes over time.

I used AI-assisted coding to help structure the code. I am not good at creating complicated code; however, the code failed to track the droplet accurately. The binary image is significantly different from the actual droplet. I am attaching the code, along with a video link of a ball as an example (valid for 3 days: https://we.tl/t-5hMniU7T6o), and an image from the first frame of the video. Maybe additional algorithms are required.

Thank you very much!

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

% ======= Parameters =======
videoFile = 'Q1M 1156-N139.avi';
ballDiameterMM = 1;
ballDiameterPixels = 324;
scale = ballDiameterMM / ballDiameterPixels; % mm/pixel
fprintf('Using scale: %.6f mm/pixel\n', scale);
% ======= Read video =======
vidObj = VideoReader(videoFile);
numFrames = floor(vidObj.FrameRate * vidObj.Duration);
fprintf('Total frames: %d\n', numFrames);
% ======= Manual center selection =======
vidObj.CurrentTime = 0;
firstFrame = readFrame(vidObj);
if size(firstFrame, 3) == 3
    firstFrameGray = rgb2gray(firstFrame);
else
    firstFrameGray = firstFrame;
end
figure('Name', 'Select Ball Center');
imshow(firstFrameGray);
title('Click the center of the ball, then press Enter');
[xManual, yManual] = ginput(1);
hold on; plot(xManual, yManual, 'r+', 'MarkerSize', 15, 'LineWidth', 2); hold off;
manualCenter = [xManual, yManual];
close;
% ======= Manual ellipse reference selection =======
figure('Name', 'Define Ellipse Diameters');
imshow(firstFrameGray);
title('Click two points along the **major axis**, then press Enter');
[majX, majY] = ginput(2);
hold on; plot(majX, majY, 'g-', 'LineWidth', 2);
title('Click two points along the **minor axis**, then press Enter');
[minX, minY] = ginput(2);
plot(minX, minY, 'c-', 'LineWidth', 2); hold off; close;
% Compute reference ellipse parameters
Xc_ref = mean(majX);
Yc_ref = mean(majY);
a_ref = sqrt((majX(2) - majX(1))^2 + (majY(2) - majY(1))^2) / 2;
b_ref = sqrt((minX(2) - minX(1))^2 + (minY(2) - minY(1))^2) / 2;
phi_ref = atan2(majY(2) - majY(1), majX(2) - majX(1));  % radians
fprintf('Reference Ellipse: Center=(%.2f, %.2f), a=%.2f px, b=%.2f px, phi=%.2f°\n', ...
    Xc_ref, Yc_ref, a_ref, b_ref, rad2deg(phi_ref));
% ======= Manual ROI selection for enhancement =======
figure('Name', 'Select ROI for Enhancement');
imshow(firstFrameGray);
title('Draw a polygon around the ball for enhancement, then double-click to finish');
roiMask = roipoly();
close;
% Ensure ROI has sufficient size
if sum(roiMask(:)) < 4
    warning('ROI is too small. Using full image for enhancement.');
    roiMask = true(size(firstFrameGray));
end
% ======= Preallocate =======
centroids = nan(numFrames, 3); % Includes x, y, z coordinates
majorAxisLengths = nan(numFrames, 1);
minorAxisLengths = nan(numFrames, 1);
orientations = nan(numFrames, 1);
% ======= Output video setup for combined 2D and 3D visualization =======
outputVideoFile = 'output_combined.avi';
outputVid = VideoWriter(outputVideoFile, 'Uncompressed AVI');
outputVid.FrameRate = vidObj.FrameRate;
open(outputVid);
% ======= Set figure positions =======
figWidth = 700; % Width for each figure
figHeight = 420;
% Position figure 1 (2D) on the left
figure(1);
set(gcf, 'Position', [100, 100, figWidth, figHeight]);
% Position figure 2 (3D) on the right, next to figure 1
figure(2);
set(gcf, 'Position', [100 + figWidth + 10, 100, figWidth, figHeight]);
% ======= Process each frame =======
vidObj.CurrentTime = 0;
for k = 1:numFrames
    frame = readFrame(vidObj);
    if size(frame, 3) == 3
        grayFrame = rgb2gray(frame);
    else
        grayFrame = frame;
    end
    % Dynamically resize and validate ROI mask for current frame
    currentMask = imresize(roiMask, size(grayFrame), 'nearest');
    if sum(currentMask(:)) < 4
        currentMask = true(size(grayFrame)); % Fallback to full image
        warning('ROI too small for frame %d. Using full image.', k);
    end
    % Apply enhancement and brightening with robust handling
    enhancedFrame = grayFrame;
    enhancedFrame(~currentMask) = imadjust(enhancedFrame(~currentMask), stretchlim(enhancedFrame(~currentMask), [0.2 0.8]), [0 1]);
    % Use imadjust for ROI if adapthisteq would fail
    enhancedFrame(currentMask) = imadjust(enhancedFrame(currentMask), stretchlim(enhancedFrame(currentMask), [0.2 0.8]), [0 1]);
    brightBoosted = imadjust(enhancedFrame, stretchlim(enhancedFrame, [0.2 0.8]), [0 1]);
    filteredFrame = medfilt2(brightBoosted, [3 3]);
    bw = imbinarize(filteredFrame, 'adaptive', 'ForegroundPolarity', 'bright', 'Sensitivity', 0.6);
    bw = bwareaopen(bw, 50);
    % Create reference mask based on pink dashed ellipse
    [x, y] = meshgrid(1:size(grayFrame, 2), 1:size(grayFrame, 1));
    theta = atan2(y - Yc_ref, x - Xc_ref) - phi_ref;
    r = ((x - Xc_ref) * cos(phi_ref) + (y - Yc_ref) * sin(phi_ref)).^2 / a_ref^2 + ...
        (-(x - Xc_ref) * sin(phi_ref) + (y - Yc_ref) * cos(phi_ref)).^2 / b_ref^2;
    refMask = r <= 1;
    bw = bw & refMask; % Apply mask to limit to reference boundary
    bw = imdilate(bw, strel('disk', 1));
    cc = bwconncomp(bw);
    stats = regionprops(cc, 'Centroid', 'Area', 'Orientation', 'MajorAxisLength', 'MinorAxisLength');
    if isempty(stats)
        if k > 1 && ~isnan(centroids(k-1,1))
            centroids(k,:) = centroids(k-1,:); % Hold last known position
        else
            centroids(k,:) = [NaN, NaN, NaN];
        end
        majorAxisLengths(k) = NaN;
        minorAxisLengths(k) = NaN;
        orientations(k) = NaN;
        continue;
    end
    [~, ballIdx] = max([stats.Area]); % Use largest component within mask
    selected = stats(ballIdx);
    % Ensure ellipse is slightly larger
    selected.MajorAxisLength = selected.MajorAxisLength * 1.1;
    selected.MinorAxisLength = selected.MinorAxisLength * 1.1;
    centroids(k,1:2) = selected.Centroid;
    majorAxisLengths(k) = selected.MajorAxisLength * scale;
    minorAxisLengths(k) = selected.MinorAxisLength * scale;
    orientations(k) = rad2deg(-selected.Orientation);
    % Estimate z-coordinate
    if ~isnan(minorAxisLengths(k)) && minorAxisLengths(k) > 0
        z = (b_ref / (selected.MinorAxisLength * scale / b_ref)) * ballDiameterMM;
        z = min(max(z, 0), ballDiameterMM * 10);
    else
        z = ballDiameterMM;
    end
    centroids(k, 3) = z;
    % ======= Display results for 2D visualization (Figure 1) =======
    figure(1); clf;
    subplot(1,3,1);
    imshow(grayFrame); title('Original Grayscale');
    subplot(1,3,2);
    imshow(brightBoosted); title('Enhanced + Brightened');
    subplot(1,3,3);
    imshow(bw); hold on;
    plot(centroids(k,1), centroids(k,2), 'ro', 'MarkerSize', 10, 'LineWidth', 2);
    % Draw detected ellipse (yellow)
    theta = linspace(0, 2*pi, 100);
    a = selected.MajorAxisLength / 2;
    b = selected.MinorAxisLength / 2;
    Xc = centroids(k,1);
    Yc = centroids(k,2);
    phi = deg2rad(-orientations(k));
    x_ellipse = Xc + a * cos(theta) * cos(phi) - b * sin(theta) * sin(phi);
    y_ellipse = Yc + a * cos(theta) * sin(phi) + b * sin(theta) * cos(phi);
    plot(x_ellipse, y_ellipse, 'y-', 'LineWidth', 2);
    % Draw reference ellipse (magenta dashed)
    x_ref = Xc_ref + a_ref * cos(theta) * cos(phi_ref) - b_ref * sin(theta) * sin(phi_ref);
    y_ref = Yc_ref + a_ref * cos(theta) * sin(phi_ref) + b_ref * sin(theta) * cos(phi_ref);
    plot(x_ref, y_ref, 'm--', 'LineWidth', 1.5);
    title(sprintf('Binary + Ellipses (Frame %d/%d)', k, numFrames));
    hold off; drawnow;
    % Ensure figure is rendered before capturing
    pause(0.1); % Increased delay
    frame2D = getframe(figure(1));
    if isempty(frame2D.cdata)
        warning('No frame captured for 2D at frame %d. Skipping write.', k);
        continue;
    end
    % ======= 3D Visualization (Figure 2) =======
    figure(2); clf;
    [X, Y, Z] = sphere(20);
    a_radius = (selected.MajorAxisLength * scale) / 2;
    b_radius = (selected.MinorAxisLength * scale) / 2;
    c_radius = ballDiameterMM / 2;
    X = X * a_radius;
    Y = Y * b_radius;
    Z = Z * c_radius;
    phi = deg2rad(-orientations(k));
    X_rot = X * cos(phi) - Y * sin(phi);
    Y_rot = X * sin(phi) + Y * cos(phi);
    X_rot = X_rot + centroids(k,1) * scale;
    Y_rot = Y_rot + centroids(k,2) * scale;
    Z = Z + z;
    surf(X_rot, Y_rot, Z, 'FaceColor', 'b', 'EdgeColor', 'none', 'FaceAlpha', 0.8);
    hold on;
    plot3(centroids(k,1) * scale, centroids(k,2) * scale, z, ...
        'ro', 'MarkerSize', 10, 'LineWidth', 2);
    axis equal;
    xlabel('X (mm)'); ylabel('Y (mm)'); zlabel('Z (mm)');
    title(sprintf('3D Ball Tracking (Frame %d/%d)', k, numFrames));
    grid on;
    currentCentroid = centroids(k, :) * scale;
    if ~isnan(currentCentroid(1))
        xRange = currentCentroid(1) + [-ballDiameterMM*5, ballDiameterMM*5];
        yRange = currentCentroid(2) + [-ballDiameterMM*5, ballDiameterMM*5];
        zRange = [0, z + ballDiameterMM*5];
    else
        xRange = [0, vidObj.Width * scale];
        yRange = [0, vidObj.Height * scale];
        zRange = [0, ballDiameterMM * 10];
    end
    axis([xRange, yRange, zRange]);
    lighting gouraud;
    camlight;
    view(45, 30);
    hold off; drawnow;
    pause(0.1); % Increased delay
    frame3D = getframe(figure(2));
    if isempty(frame3D.cdata)
        warning('No frame captured for 3D at frame %d. Skipping write.', k);
        continue;
    end
    targetHeight = min(size(frame2D.cdata, 1), size(frame3D.cdata, 1));
    frame2D_resized = imresize(frame2D.cdata, [targetHeight, NaN]);
    frame3D_resized = imresize(frame3D.cdata, [targetHeight, NaN]);
    combinedFrame = cat(2, frame2D_resized, frame3D_resized);
    writeVideo(outputVid, combinedFrame);
end
% ======= Close video writer =======
close(outputVid);
fprintf('Combined video saved to %s\n', outputVideoFile);
% ======= Save video to SavedVideos directory =======
outputDir = 'SavedVideos';
if ~exist(outputDir, 'dir')
    mkdir(outputDir);
end
try
    movefile(outputVideoFile, fullfile(outputDir, outputVideoFile));
    fprintf('Moved combined video to %s\n', fullfile(outputDir, outputVideoFile));
catch e
    fprintf('Error moving video to %s: %s\n', outputDir, e.message);
end
% ======= Displacement Analysis =======
validIdx = ~isnan(centroids(:,1));
timeVec = (0:numFrames-1)' / vidObj.FrameRate;
startIdx = find(validIdx, 1, 'first');
x0 = centroids(startIdx,1);
y0 = centroids(startIdx,2);
displacementX = centroids(:,1) - x0;
displacementY = centroids(:,2) - y0;
displacementX_mm = displacementX * scale;
displacementY_mm = displacementY * scale;
figure;
subplot(2,1,1)
plot(timeVec(validIdx), displacementX_mm(validIdx), 'b-');
xlabel('Time (s)'); ylabel('Displacement X (mm)');
title('X Displacement of Styrofoam Ball');
subplot(2,1,2)
plot(timeVec(validIdx), displacementY_mm(validIdx), 'r-');
xlabel('Time (s)'); ylabel('Displacement Y (mm)');
title('Y Displacement of Styrofoam Ball');
figure;
plot(displacementX_mm(validIdx), displacementY_mm(validIdx), 'ko-', 'MarkerFaceColor', 'k');
xlabel('X Displacement (mm)'); ylabel('Y Displacement (mm)');
title('2D Trajectory of Styrofoam Ball');
axis equal; grid on;
% ======= FFT Analysis =======
Fs = vidObj.FrameRate;
X_valid = displacementX_mm(validIdx);
Y_valid = displacementY_mm(validIdx);
N = length(X_valid);
f = (0:N-1)*(Fs/N);
X_centered = X_valid - mean(X_valid);
Y_centered = Y_valid - mean(Y_valid);
X_FFT = fft(X_centered); Y_FFT = fft(Y_centered);
magX = abs(X_FFT)/N; magY = abs(Y_FFT)/N;
ylimMax = max([magX; magY]);
figure;
subplot(2,1,1); plot(f, magX); xlabel('Frequency (Hz)');
ylabel('Magnitude'); title('FFT of X Displacement'); ylim([0 ylimMax]);
subplot(2,1,2); plot(f, magY); xlabel('Frequency (Hz)');
ylabel('Magnitude'); title('FFT of Y Displacement'); ylim([0 ylimMax]);
figure;
subplot(2, 1, 1);
plot(timeVec(validIdx), majorAxisLengths(validIdx), 'k-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Major Axis Length (mm)');
title('Major Axis Length Over Time');
subplot(2, 1, 2);
plot(timeVec(validIdx), minorAxisLengths(validIdx), 'k-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Minor Axis Length (mm)');
title('Minor Axis Length Over Time');

[Image Analyst]

Does the ball look very similar to your droplet? If not, getting something to work with this ball video may fail completely for your actual droplet video.

Note that your binarization says 'ForegroundPolarity', 'bright' while your forground is actually dark.

[Zahra]

Please check the video link showing the droplet (Valid for 3 days : https://we.tl/t-MsW6Jn9h8X ). Thank you.

[Image Analyst]

It's already expired.

[Zahra]

Please check the new link (Valid for 3 days : https://we.tl/t-AmAIvcYUYi ).