Try this:
spray_angle_demo
Output:
Spray angle: 45.98 degrees
function spray_angle_demo
% Sample binary image
im=imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/122806/im1.jpg');
bw=max(im,[],3)>100;
clear im
% Get the largest foreground object (in case there is more than one)
L=bwlabel(bw);
N=max(L(:));
if N>1
RP=regionprops(L,'Area');
A=zeros(N,1);
for i=1:N, A(i)=RP(i).Area; end
[~,id]=max(A);
bw=L==id;
end
clear L
figure
imshow(bw)
axis on
hold on
% Get boundary coordinates
XY=bwboundaries(bw,8,'noholes');
XY=XY{1}(:,[2 1]); % make x coordinate the 1st one
% Identify left and right edges of the triangular spray profile
% -------------------------------------------------------------------------
% Boundary vertices comprising the convex hull (CH)
CH=convhull(XY,'simplify',true);
CH=flipud(CH); % counter-clocwise order when superiposed on the image
B=XY(CH,:);
plot(B(:,1),B(:,2),'--r','LineWidth',2)
% Top-most vertex
B(end,:)=[]; % remove duplicate vertex
[~,id_top]=min(B(:,2));
B=circshift(B ,[1-id_top 0]);
% Left edge (vertices ordered from top to bottom as viewed on the image)
[~,id_lft]=min(B(:,1)); % left-most vertex
B_lft=B(1:id_lft,:);
%plot(B_lft(:,1),B_lft(:,2),'-g','LineWidth',2)
%plot(B_lft(:,1),B_lft(:,2),'.g','MarkerSize',20)
% Right edge (vertices ordered from top to bottom as viewed on the image)
[~,id_rgt]=max(B(:,1)); % right-most vertex
B_rgt=B;
B_rgt(2:(id_rgt-1),:)=[];
B_rgt=circshift(B_rgt,[-1 0]);
B_rgt=flipud(B_rgt);
%plot(B_rgt(:,1),B_rgt(:,2),'-b','LineWidth',2)
%plot(B_rgt(:,1),B_rgt(:,2),'.b','MarkerSize',20)
% Identify straight portions of the left and right edges
% -------------------------------------------------------------------------
E_lft=straight_edge_segment(B_lft);
E_rgt=straight_edge_segment(B_rgt);
% Visualize solution
% -------------------------------------------------------------------------
P=LineIntersection(E_lft(1,:),E_lft(2,:),E_rgt(1,:),E_rgt(2,:)); % point of intersection
E_lft(1,:)=P;
E_rgt(1,:)=P;
plot(E_lft(:,1),E_lft(:,2),'-g','LineWidth',3)
plot(E_rgt(:,1),E_rgt(:,2),'-g','LineWidth',3)
plot(P(1),P(2),'*y','MarkerSize',15,'LineWidth',3)
d_lft=E_lft(2,:)-E_lft(1,:);
d_lft=d_lft/norm(d_lft);
d_rgt=E_rgt(2,:)-E_rgt(1,:);
d_rgt=d_rgt/norm(d_rgt);
t=(180/pi)*acos(dot(d_lft,d_rgt));
fprintf('Spray angle: %.2f degrees\n',t)
% =========================================================================
function E=straight_edge_segment(B)
N=size(B,1);
E=B(1:2,:);
if N==2, return; end
% Lenghts of boundary segments comprising B
D=B(2:N,:)-B(1:(N-1),:);
L=sqrt(sum(D.^2,2));
% Accumulate vertices
L_net=sum(L);
L=L(1);
i=2;
f=L/L_net;
df=0;
t_thr=5;
while f<.99 && df<(2/3)
i=i+1;
L=L + norm(B(i,:)-E(end,:));
fi=L/L_net;
df=fi-f;
f=fi;
De=E(end,:)-E(1,:);
De=De/norm(De);
Di=B(i,:)-E(end,:);
Di=Di/norm(Di);
t=dot(De,Di);
t(t>1)=1;
t=acos(t)*(180/pi);
%disp([i t f df])
if t<t_thr
E=cat(1,E,B(i,:));
elseif size(E,1)>2
De=E(end,:)-E(2,:);
De=De/norm(De);
Di=B(i,:)-E(end,:);
Di=Di/norm(Di);
t=dot(De,Di);
t(t>1)=1;
t=acos(t)*(180/pi);
%disp([i t f])
if t<t_thr
E=cat(1,E,B(i,:));
E(1,:)=[];
else
break
end
elseif f<0.99
E=cat(1,E,B(i,:));
E(1,:)=[];
else
break
end
end
% End-points
E=cat(1,E(1,:),E(end,:));
% =========================================================================
function P=LineIntersection(A1,A2,B1,B2)
% Find points of intersection between line pairs rdefined by point tuples.
%
% INPUT:
% - A1,A2 : N-by-d arrays containing end point coordinates of N line
% segments in set A
% - B1,B2 : N-by-d arrays containing end point coordinates of N line
% segments in set B
%
% OUTPUT:
% - Pa,Pb : N-by-d arrays containing coordinates of points of
% intersection between corresponding lines is sets A and B.
Daa=sum((A2-A1).^2,2);
Dbb=sum((B2-B1).^2,2);
Dab=sum((A2-A1).*(B2-B1),2);
N=size(A1,1);
S=zeros(2,2,N);
S(1,1,:)= Daa(:)+eps;
S(1,2,:)=-Dab(:);
S(2,1,:)=-Dab(:);
S(2,2,:)= Dbb(:)+eps;
g=[sum((A2-A1).*(A1-B1),2) -sum((B2-B1).*(A1-B1),2)];
g=permute(g,[2 3 1]);
t=Cramer_2by2_Inverse(S,-g);
t=permute(t,[3 1 2]);
%t(t<0)=0; %t(t>1)=1; % apply constraits to get closest points between line segments
%Pa=bsxfun(@times,t(:,1),A2-A1) + A1;
%Pb=bsxfun(@times,t(:,2),B2-B1) + B1;
P=bsxfun(@times,t(:,1),A2-A1) + A1;
% =========================================================================
function uv=Cramer_2by2_Inverse(A,b)
% Get solution to a 2-by-2 linear system using Cramer's rule
%
% - A : 2-by-2-by-N stack of N matrices
% - b : 2-by-1-by-N stack of N column vectors
% - uv : 2-by-1-by-N stack of N column vectors, so that uv(:,:,i)=A(:,:,i)\b(:,:,i)
% detA = A(1,1)*A(2,2) - A(1,2)*A(2,1)
% detU = b(1)*A(2,2) - b(2)*A(1,2)
% detV = b(2)*A(1,1) - b(1)*A(2,1)
detA=A(1,1,:).*A(2,2,:)-A(1,2,:).*A(2,1,:);
% u=detU/det(A)
u=(b(1,:,:).*A(2,2)-b(2,:,:).*A(1,2))./detA;
% v=detV/det(A)
v=(b(2,:,:).*A(1,1)-b(1,:,:).*A(2,1))./detA;
% uv=A\b
uv=cat(1,u,v);
