contact_angle_left = NaN;
contact_angle_right = NaN;
dropletMask = false(size(filledImage));
dropletMask(dropletStats.PixelIdxList) = true;
dropletBoundary = bwboundaries(dropletMask);
if ~isempty(dropletBoundary)
outline = dropletBoundary{1};
x_outline = outline(:,2);
y_outline = outline(:,1);
if numel(x_outline) >= 20
x_phys = x_outline * scale_Image * 1e-6;
y_phys = y_outline * scale_Image * 1e-6;
surface_y_phys = minYSurface * scale_Image * 1e-6;
contact_threshold = equiv_radius * 0.3;
near_surface_indices = find(abs(y_phys - surface_y_phys) < contact_threshold);
if length(near_surface_indices) >= 4
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)];
tangent_range = equiv_radius * 0.2;
left_contact_idx = near_surface_indices(left_idx);
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
above_contact = tangent_indices_left(y_phys(tangent_indices_left) < left_contact_point(2));
if length(above_contact) >= 3
x_fit = x_phys(above_contact);
y_fit = y_phys(above_contact);
[x_fit, sort_idx] = sort(x_fit);
if length(unique(x_fit)) >= 2 && length(x_fit) >= 2
p_left = polyfit(x_fit, y_fit, 1);
angle_rad_left = atan(slope_left);
contact_angle_left = 180 - rad2deg(angle_rad_left);
if contact_angle_left < 0
contact_angle_left = contact_angle_left + 180;
elseif contact_angle_left > 180
contact_angle_left = contact_angle_left - 180;
right_contact_idx = near_surface_indices(right_idx);
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
above_contact = tangent_indices_right(y_phys(tangent_indices_right) < right_contact_point(2));
if length(above_contact) >= 3
x_fit = x_phys(above_contact);
y_fit = y_phys(above_contact);
[x_fit, sort_idx] = sort(x_fit);
if length(unique(x_fit)) >= 2 && length(x_fit) >= 2
p_right = polyfit(x_fit, y_fit, 1);
slope_right = p_right(1);
angle_rad_right = atan(slope_right);
contact_angle_right = rad2deg(angle_rad_right);
if contact_angle_right < 0
contact_angle_right = contact_angle_right + 180;
elseif contact_angle_right > 180
contact_angle_right = contact_angle_right - 180;