Solved question video analysis

Hey all, so i have this script that analyzes the contact angle of a sessile droplet right now i do it with the code attached below and on the video attached to this. The issue is mainly that the aquired contact angles are quite jaggery and inconsitent since it is depended on the surface selection could anyone help me solve this?(i edited the question to make my issue and the answer more logical for readers)

        contact_angle_left = NaN;
        contact_angle_right = NaN;
        contact_angle_avg = NaN;

        if is_impact
            dropletMask = false(size(filledImage));
            dropletMask(dropletStats.PixelIdxList) = true;

            % Find droplet boundary
            dropletBoundary = bwboundaries(dropletMask);
            if ~isempty(dropletBoundary)
                outline = dropletBoundary{1};
                x_outline = outline(:,2); % Column indices (x-coordinates)
                y_outline = outline(:,1); % Row indices (y-coordinates)

                if numel(x_outline) >= 20 % Need sufficient points
                    % Convert to physical coordinates (meters)
                    x_phys = x_outline * scale_Image * 1e-6;
                    y_phys = y_outline * scale_Image * 1e-6;

                    % Find contact points (points near the surface)
                    surface_y_phys = minYSurface * scale_Image * 1e-6;
                    contact_threshold = equiv_radius * 0.3; % 30% of droplet radius

                    % Find all points near the surface
                    near_surface_indices = find(abs(y_phys - surface_y_phys) < contact_threshold);

                    if length(near_surface_indices) >= 4
                        % Find the leftmost and rightmost contact points
                        x_contact = x_phys(near_surface_indices);
                        y_contact = y_phys(near_surface_indices);

                        [x_min, left_idx] = min(x_contact);
                        [x_max, right_idx] = max(x_contact);

                        left_contact_point = [x_min, y_contact(left_idx)];
                        right_contact_point = [x_max, y_contact(right_idx)];

                        % Calculate contact angles using tangent method
                        % For each contact point, find nearby points to calculate tangent
                        tangent_range = equiv_radius * 0.2; % Distance for tangent calculation

                        % LEFT CONTACT ANGLE
                        left_contact_idx = near_surface_indices(left_idx);
                        % Find points within tangent range from left contact point
                        distances_from_left = sqrt((x_phys - left_contact_point(1)).^2 + ...
                            (y_phys - left_contact_point(2)).^2);
                        tangent_indices_left = find(distances_from_left <= tangent_range & ...
                            distances_from_left > 0);

                        if length(tangent_indices_left) >= 3

                            % Filter points that are above the contact point (droplet side)
                            above_contact = tangent_indices_left(y_phys(tangent_indices_left) < left_contact_point(2));

                            if length(above_contact) >= 3
                                % Fit line to these points
                                x_fit = x_phys(above_contact);
                                y_fit = y_phys(above_contact);

                                % Remove outliers and sort by x
                                [x_fit, sort_idx] = sort(x_fit);
                                y_fit = y_fit(sort_idx);

                                % Linear fit with error handling
                                if length(unique(x_fit)) >= 2 && length(x_fit) >= 2
                                    try
                                        p_left = polyfit(x_fit, y_fit, 1);
                                        slope_left = p_left(1);
                                    catch
                                        slope_left = NaN;
                                    end
                                else
                                    slope_left = NaN;
                                end

                                % Contact angle calculations with error handling
                                if ~isnan(slope_left)
                                    angle_rad_left = atan(slope_left);
                                    contact_angle_left = 180 - rad2deg(angle_rad_left);

                                    % Ensure angle is in valid range
                                    if contact_angle_left < 0
                                        contact_angle_left = contact_angle_left + 180;
                                    elseif contact_angle_left > 180
                                        contact_angle_left = contact_angle_left - 180;
                                    end
                                end
                            end
                        end

                        % RIGHT CONTACT ANGLE
                        right_contact_idx = near_surface_indices(right_idx);
                        % Find points within tangent range from right contact point
                        distances_from_right = sqrt((x_phys - right_contact_point(1)).^2 + ...
                            (y_phys - right_contact_point(2)).^2);
                        tangent_indices_right = find(distances_from_right <= tangent_range & ...
                            distances_from_right > 0);

                        if length(tangent_indices_right) >= 3
                            % Filter points that are above the contact point (droplet side)
                            above_contact = tangent_indices_right(y_phys(tangent_indices_right) < right_contact_point(2));

                            if length(above_contact) >= 3
                                % Fit line to these points
                                x_fit = x_phys(above_contact);
                                y_fit = y_phys(above_contact);


                                % Remove outliers and sort by x
                                [x_fit, sort_idx] = sort(x_fit);
                                y_fit = y_fit(sort_idx);

                                % Linear fit with error handling
                                if length(unique(x_fit)) >= 2 && length(x_fit) >= 2
                                    try
                                        p_right = polyfit(x_fit, y_fit, 1);
                                        slope_right = p_right(1);
                                    catch
                                        slope_right = NaN;
                                    end
                                else
                                    slope_right = NaN;
                                end

                                % Contact angle calculations with error handling
                                if ~isnan(slope_right)
                                    angle_rad_right = atan(slope_right);
                                    contact_angle_right = rad2deg(angle_rad_right);

                                    % Ensure angle is in valid range
                                    if contact_angle_right < 0
                                        contact_angle_right = contact_angle_right + 180;
                                    elseif contact_angle_right > 180
                                        contact_angle_right = contact_angle_right - 180;
                                    end
                                end
                            end
                        end

Answer 1

[Mathieu NOE]

hello

maybe this ? NB, this is a no Image Processing Tbx solution (as I don't have it) but I wanted to give it a try

for the small video (8 frames if I'm correct) that you posted here are my results (in degrees ) :

angle_left = 81.9921

82.0263

82.4190

81.5316

82.5122

82.5768

81.7864

81.7500

angle_right = -80.5377

-80.8377

-80.6901

-81.5111

-80.6901

-80.9807

-81.3844

-80.8377

the code basically do a gradient of the image , smooth it a bit (I have used this Fex submission : smooth2a - File Exchange - MATLAB Central) , then I do top (red dots) and bottom (yellow dots) boundary of the gradient data (as well as the "mean" curve => green dots)

personnaly I opted to use the top boundary points and taking the first (left) 5 samples (same on the right side) you can use polyfit to estimate the tangent slope - from where you can get the angles , of course. The tangent vector is displayed as the magenta segment

zoom on the LHS of the plot :

zoom on the RHS of the plot :

NB that the results may vary according to the threshold used to define the boundary points

hope it helps !

code :

% Create a VideoReader object for the AVI file
vidObj = VideoReader('4dropletshorter.avi');

figure(1)
axis equal

% Read video frames until available
k = 0;
  while hasFrame(vidObj)
      k = k +1;
      vidFrame = readFrame(vidObj);
      frame = double(rgb2gray(vidFrame))/255;
      [angle_left(k,1),angle_right(k,1)] = do_the_stuff(frame); % main function 
      pause(1);

  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [angle_left,angle_right] = do_the_stuff(frame)

      A = flipud(frame); % to have image displayed with correct y direction
      B = abs(gradient(A)); % abs gradient
      % smooth the gradient
      nx = 10;
      ny = 1;
      B = smooth2a(B,nx,ny); % Fex : https://fr.mathworks.com/matlabcentral/fileexchange/23287-smooth2a?s_tid=srchtitle
      B = B -min(B,[],'all'); % min is now 0
      B = B./max(B,[],"all"); % normalize 0 to 1


    [y,x] = find(B>0.25); % adjust the threshold !!
    [xtop,ytop,xbottom,ybottom,xm,ym] = top_bottom_boundary(x,y);
    
    % figure(2)
    imagesc(B);
    colorbar('vert');
    set(gca,'YDir','normal');
    colormap('gray');
    hold on
    plot(xtop, ytop, '.r');         % TOP boundary
    plot(xm,ym,'.g');               % mean curve
    plot(xbottom,ybottom, '.y');    % BOTTOM boundary
    
    % tangent vector left side
    samples = 5;
    p = polyfit(xtop(1:samples), ytop(1:samples),1);
    f = polyval(p,xtop(1:samples));
    plot(xtop(1:samples),f,'m-','LineWidth',2)
    angle_left = atan(p(1))*180/pi;
    
    % tangent vector right side
    p = polyfit(xtop(end-samples+1:end), ytop(end-samples+1:end),1);
    f = polyval(p,xtop(end-samples+1:end));
    plot(xtop(end-samples+1:end),f,'m-','LineWidth',2)
    angle_right = atan(p(1))*180/pi;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x3,y3,x4,y4,x5,y5] = top_bottom_boundary(x,y)
% based on FEX : https://fr.mathworks.com/matlabcentral/answers/299796-tight-boundary-around-a-set-of-points
    %split data into classes and find max/min for each class
    class_label = unique(x);
    upper_boundary = zeros(size(class_label));
    lower_boundary = zeros(size(class_label));
    for idx = 1:numel(class_label)
        class = y(x == class_label(idx));
        upper_boundary(idx) = max(class);
        lower_boundary(idx) = min(class);

        % mean curve computation : 
        indy = diff(class)<2; % indice of contiguous y data 
        [begin,ends] = find_start_end_group(indy);
        duration = ends - begin; % in samples
        [md,ii] = max(duration);

        % keep only mean values with y dist below threshold
        if md>2 % let's find contiguous block of pixel that are at least 3 samples long
            mean_curve(idx) = mean(class(begin(ii):ends(ii)));
        else
            mean_curve(idx) = NaN;
        end
        

    end
%     left_boundary = y(x == class_label(1));
%     right_boundary = y(x == class_label(end));
% 
%     % left_boundary
%     x1 = class_label(1)*ones(size(left_boundary));
%     y1 = left_boundary;
% 
%     % right_boundary
%     x2 = class_label(end)*ones(size(right_boundary));
%     y2 = right_boundary;

    % top boundary
    x3 = class_label;
    y3 = upper_boundary;
    % bottom boundary
    x4 = class_label;
    y4 = lower_boundary;
    % mean y curve
    x5= class_label;
    y5 = mean_curve;    
end
 
 

[Tom]

Thank you so much for the answer, I think it worked and will try to implement it in my code since it seemed to work much better than before. If i have any question i will let you know

[Mathieu NOE]

attached the function "find_start_end_group.m" I forgot to share previously

all the best