How to evaluate the angle from an image?

I have several of the following images from an experimental campaign. I need to evaluate the teta angle for each of them.

My approach was to manually evaluate 3 points and obtain the angle between the 2 vectors (please see the code below). I was wondering how to introduce a methodology to automate the procedure.

%spray cone angle evaluation

[xa,ya,za]=improfile(a2,[x0 x1],[y0 y1]); ia=sub2ind(s,round(ya),round(xa));aa=a2;aa(ia)=255; [xa1,ya1,za1]=improfile(a2,[x0 x2],[y0 y2]); ia1=sub2ind(s,round(ya1),round(xa1));aa(ia1)=255; p1=[x1 y1];p0=[x0 y0]; p2=[x2 y2];v1=p1-p0;v2=p2-p0; costheta=dot(v1,v2)/(norm(v1)*norm(v2)); thetadeg=acosd(costheta); figure(1);imshow(aa);title(strcat('teta=',num2str(thetadeg)));

Answer 1

[Anton Semechko]

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);

Answer 2

[Davide Domenico Sciortino]

Thank you Anton Semechko. It seems working fine.