My name is Baraa, I recently started using Matlab (Bigginner), and I need your help in explaining the error I keep making in the following code. The code is supposed to measure the contact angle of a droplet based on defining the edges of the droplet, masking the capillary, dividing the pixels on the edges of the droplet by 2 (so we can have left and right), finding the reflection of each side to define the contact points between the surface and the droplet on each side, and finally connecting the contact points together. The code then fits the pixels to find the tangent and calculates the contact angle between the tangent and the line between the contact points.
The main problem is that the code doesn't work with all droplets and most of the time I have errors with droplets with high wettability and high volume (the droplets spread on the surface). The error message is as follows:
Index in position 1 exceeds array bounds. Index must not exceed 1113.
[x0L,y0L,indexL]=findreflection([edgeL.x,edgeL.y],70,0); % change the last entry (160) to zero and see what happens
The code always fails to read this image (Co-Eth-1, see attached) and it reads probably reads this image (Co-Eth-1b, see attached). Please refer to the two images.
clc
clear all
close all
addpath(genpath("C:\Users\Bara'a Al-khateeb\Desktop\Others\Students\Gideon Mensah\Laptop\codes\Contact Angle Fitting\Scripts\Subpixel Matlab v2.11"));
addpath(genpath("C:\Users\Bara'a Al-khateeb\Desktop\Others\Students\Gideon Mensah\Laptop\codes\Contact Angle Fitting\Scripts"));
main_path = "C:\Users\Bara'a Al-khateeb\Desktop\Others\Students\Gideon Mensah\Laptop\IMAGES\25.06.24";
output_folder = "C:\Users\Bara'a Al-khateeb\Desktop\Others\Students\Gideon Mensah\Laptop\sessileDrop_processed\";
%these are the coordinates of the polygon to mask the cannula
Xright=3100;%this the X value on the right side of the ploygon and it represents top and bottom if the right side of the polygon.
Xleft=2800;%this is the X value of the left side of the ploygon and it represents both top and bottom x cordinates
Y_top=0;% this is the y value on the top side of the polygon for both left and right side
Y_bottom=2050;% this is the y value of the bottom side of the ploygon for both left and right side
%******************************pixel intensity***************************
pixelValue=255;
%***********************threshold to be used*****************************
trsh_sub=4;
% Read a group of images from a folder
imageFiles = dir(fullfile(main_path, '*.png')); % Assuming images are in PNG format
numImages = input('Enter the number of images to work with: ');
% Create cell arrays to store data
imageNumbers = zeros(numImages, 1);
CAL_values = zeros(numImages, 1);
CAR_values = zeros(numImages, 1);
for k = 1:numImages
% Load image
filename = imageFiles(k).name;
im_load = fullfile(main_path, filename);
im = imread(im_load);
%-----------------------------------------------------------------
% Masking the cannula
%-----------------------------------------------------------------
% Define the coordinates of the four points
maskSize = size(im);
pointCoordinates = [Xleft, Y_bottom; Xright, Y_bottom; Xright, Y_top; Xleft, Y_top];
mask = zeros(maskSize(1), maskSize(2));
mask = poly2mask(pointCoordinates(:, 1), pointCoordinates(:, 2), maskSize(1), maskSize(2));
burnedImage = im;
burnedImage(mask) = pixelValue;
%---------------------------------------------------------------
% Step 2 Detect boundaries in image
%---------------------------------------------------------------
[edges, RI] = subpixelEdges(burnedImage, trsh_sub);
im_size=size(burnedImage);
im_size=size(burnedImage);
points = detectHarrisFeatures(burnedImage);
strongest = selectStrongest(points,10);
f2h=figure;
imshow(burnedImage)
hold on
plot(edges.x,edges.y,'r.','LineWidth',2)
% plot(edges.x(1), edges.y(1),'b*')
% plot(edges.x(im_size(1)),edges.y(im_size(2)),'g*')
plot(strongest)
set(f2h,'Units','normalized','Position',[0.31,0.6,0.3,0.3])
%---------------------------------------------------------------
% Step 3 Select longes boundary in image
%---------------------------------------------------------------
longestedge=findlongestedge(edges,size(burnedImage),5);
f3h=figure;
imshow(burnedImage)
hold on
plot(longestedge.x,longestedge.y,'r.','LineWidth',2)
set(f3h,'Units','normalized','Position',[0.62,0.6,0.3,0.3])
%---------------------------------------------------------------
% Step 4 Split edge into left and right and sort it
%---------------------------------------------------------------
[edgeL,edgeR]=leftrightedges(longestedge);
f4h=figure;
hold on
plot(edgeL.x,edgeL.y,'r','LineWidth',2)
plot(edgeR.x,edgeR.y,'b','LineWidth',2)
set(f4h,'Units','normalized','Position',[0.0,0.3,0.3,0.3])
%---------------------------------------------------------------
% Step 5 Find reflection
%---------------------------------------------------------------
[x0L,y0L,indexL]=findreflection([edgeL.x,edgeL.y],70,0); % change the last entry (160) to zero and see what happens
[x0R,y0R,indexR]=findreflection([edgeR.x,edgeR.y],60,0);
f5h=figure;
imshow(burnedImage)
hold on
plot(edgeL.x,edgeL.y,'r','LineWidth',2)
plot(edgeR.x,edgeR.y,'b','LineWidth',2)
plot(x0L,y0L,'yx','MarkerSize',10,'LineWidth',2)
plot(x0R,y0R,'yx','MarkerSize',10,'LineWidth',2)
t=linspace(-3,3);
plot((x0L-x0R)*t+x0R,(y0L-y0R)*t+y0R,'r--','LineWidth',2)
set(f5h,'Units','normalized','Position',[0.31,0.3,0.3,0.3])
%---------------------------------------------------------------
% Step 6 Fit data to polynomial
%---------------------------------------------------------------
PolyData=polynomialfit(edgeL,edgeR,[x0L,y0L],[x0R,y0R]);
f6h=figure;
imshow(burnedImage)
hold on
t=linspace(-3,3);
plot((x0L-x0R)*t+x0R,(y0L-y0R)*t+y0R,'r--','LineWidth',2)
plot(PolyData.EvalPolyL(:,1),PolyData.EvalPolyL(:,2),'r','LineWidth',2)
plot(PolyData.EvalPolyR(:,1),PolyData.EvalPolyR(:,2),'b','LineWidth',2)
radius=500;
tilt=atand((y0R-y0L)/(x0R-x0L));
plot([PolyData.TLL(1),PolyData.TLL(1)+radius*cosd(PolyData.CAL-tilt)],[PolyData.TLL(2),PolyData.TLL(2)-radius*sind(PolyData.CAL-tilt)],'LineWidth', 2,'color','g')
plot([PolyData.TLR(1),PolyData.TLR(1)-radius*cosd(PolyData.CAR+tilt)],[PolyData.TLR(2),PolyData.TLR(2)-radius*sind(PolyData.CAR+tilt)],'LineWidth', 2,'color','g')
legend(['contact angles, CA left= ',num2str(PolyData.CAL),' CA Right= ',num2str(PolyData.CAR)])
set(f6h,'Units','normalized','Position',[0.62,0.3,0.3,0.3])
display(['Polynomial fit return the contact angles, CA Left=',num2str(PolyData.CAL),' CA Right=',num2str(PolyData.CAR)])
% % Store data in the table
imageNumbers(k) = k;
CAL_values(k) = PolyData.CAL;
CAR_values(k) = PolyData.CAR;
% Save processed image to a new folder
if ~exist(output_folder, 'dir')
mkdir(output_folder);
end
% Save processed image in the specified formats
[~, name, ext] = fileparts(filename);
processed_filename = [name '_processed'];
% Check if the filename already exists, if so, append a unique identifier
index = 1;
while exist(fullfile(output_folder, [processed_filename, '_', num2str(index), '.jpg']), 'file')
index = index + 1;
end
% Save as JPG
imwrite(burnedImage, fullfile(output_folder, [processed_filename, '_', num2str(index), '.jpg']));
% Save as FIG
saveas(f6h, fullfile(output_folder, [processed_filename, '_', num2str(index), '.fig']));
% Save result to TXT
nsave = fullfile(output_folder, [name, '_', num2str(index), '.txt']);
result_vary = rand(10); % Example data, replace with your actual data
save(nsave, 'result_vary', '-ASCII');
disp(['Processed image saved as: ', processed_filename]);
end
% Create table from cell arrays
tableData = table(imageNumbers, CAL_values, CAR_values);
% Save table to a text file
output_filename = fullfile(output_folder, 'contactAngle_table_data.txt');
writetable(tableData, output_filename, 'Delimiter', '\t');
disp(['Table data saved as: ', output_filename]);
% Display the table
disp(tableData);
