GetDistanceUnstable.m

function [ distance  ] = GetDistanceUnstable( subPT, vertPT, deg_Increment )
% Code originally written by Skylar Gay
% and W. Scott Ingram at the University 
% of Texas M. D. Anderson Cancer 
% Center, Houston, Texas. All copies of  
% this code, even if modified, must 
% contain this message.

% No warranty or fitness for duty is
% expressed or implied by the authors. 
% All liabilities are assumed by user.    
% Usage of this code signifies that you 
% have read, understand, and agree to 
% these terms.

if numel(subPT)~=2
    A=('Cannot compute distances: subPT is not a single point!');
    disp(A)
else
if deg_Increment>=2
    %In some cases, two distances calculated by this function were
    %incorrect if deg_Increment=1. I believe it is due to rounding errors in the script where angle is calculated. Therefore, this does not allow
    %deg_Increments to equal 1. (This eliminates the need for "if
    %numel(distance)>180", but I have left it in case the problem described
    %above is eliminated.)
    if numel(vertPT)>=6
tic

theta=(0:deg_Increment:(360-deg_Increment));
%Sets the value of theta from 0 degrees to the largest angle that is not
%coterminal with 0 degrees and <360 degrees.

numVert=(numel(vertPT))/2; 

m_rays=zeros(1,length(theta));
b_rays=zeros(1,length(theta));

for a=1:length(theta)
    m_rays(a)=tand(theta(a));
    b_rays(a)=subPT(1,2)-m_rays(a)*(subPT(1,1));
end


vertPT_2=[vertPT;vertPT(1,:)];

for b=1:(numel(vertPT_2)/2)-1
    m_line(b)=(vertPT_2(b+1,2)-vertPT_2(b,2))/(vertPT_2(b+1,1)-vertPT_2(b,1));
    b_line(b)=vertPT_2(b+1,2)-(m_line(b)*vertPT_2(b+1,1));
end

for c=1:length(theta)
    clear angle
    clear x_Int y_Int
    
    for d=1:(numel(vertPT_2)/2)-1
        
        % If slope of ray is vertical (undefined, and therefore Inf)
        if abs(m_rays(c))==Inf && abs(m_rays(c))~=abs(m_line(d))
            x_Int(d)= subPT(1,1);
            y_Int(d)=(m_line(d)*x_Int(d))+b_line(d);
            angle(d) = atan2d(y_Int(d)-subPT(1,2),x_Int(d)-subPT(1,1));

        
        % If slope of line is vertical (undefined, and therefore Inf)    
        elseif abs(m_line(d))==Inf && abs(m_rays(c))~=abs(m_line(d))
            x_Int(d)=vertPT(d,1);
            y_Int(d)=(m_rays(c)*x_Int(d))+b_rays(c);
            angle(d) = atan2d(y_Int(d)-subPT(1,2),x_Int(d)-subPT(1,1));
        
        % If ray and line have equal slope (parallel)    
        elseif abs(m_rays(c))==abs(m_line(d))
            
           x_Int(d)=NaN;
           y_Int(d)=NaN;
           angle(d) = NaN;
       
        
        % Otherwise, use the normal calculation method    
        else
           x_Int(d)=(b_rays(c)-b_line(d))/(m_line(d)-m_rays(c));
           y_Int(d)=(m_line(d)*x_Int(d))+b_line(d);
           angle(d) = atan2d(y_Int(d)-subPT(1,2),x_Int(d)-subPT(1,1));
        end
        
    end
    
   angle = round(angle) + (round(angle)<0)*360;
   % Sets angle to a whole number to compensate for binary--decimal errors.
   correct_angle=find(round(angle)==theta(c));
   % Sets the correct angle to be equal to theta.
   testdistances = sqrt((x_Int(correct_angle)-repmat(subPT(1),1,length(correct_angle))).^2 + (y_Int(correct_angle)-repmat(subPT(2),1,length(correct_angle))).^2);
   % Calculates the distance along the ray whose angle corresponds to
   % theta.
   distance(c) = min(testdistances);
   % Selects the minimum distance as the correct distance.
 
    
   plotpoint(c)=min(find(testdistances==min(testdistances)));
   x_Int_plot(c)=x_Int(correct_angle(plotpoint(c)));
   y_Int_plot(c)=y_Int(correct_angle(plotpoint(c)));
   %Plots the points at which the rays intercept the perimeter.
    


end

%if numel(distance)>abs(180)
%    
%    X=('Cannot Export to Excel: numel(distance) too large to export!')
%    disp(X)
%    % I was unable to export distance to an Excel file if I used
%    % deg_Increment=1. This created 360 distances that were displayed on
%    % the Command Window but would not export them to an Excel file. I was
%    % unable to find any documentation on why this happened. deg_Increment
%    % of 2 or higher worked fine.
%    
%else
%    if exist('Distance.xls','file')==2
%    delete('Distance.xls')
%    xlswrite('Distance',distance)
%    open Distance.xls
%    else
%    xlswrite('Distance',distance)
%    open Distance.xls
%    end
%end

% Exports the distances (if possible) to a Microsoft Excel file.
hold on
grid on
 plot(vertPT_2(:,1),vertPT_2(:,2),'b-')
 plot(subPT(1),subPT(2),'ro')
 axis([min(vertPT(:,1))-3 max(vertPT(:,1))+3 min(vertPT(:,2))-3 max(vertPT(:,2))+3]) 
 axis square
 plot(x_Int_plot,y_Int_plot,'g^')
 
 for e=1:length(distance)
   plot([x_Int_plot(e), subPT(1)], [y_Int_plot(e), subPT(2)],'r--')
 end

 
 toc
    else
        Y=('Cannot compute distances: not enough vertices!');
        disp(Y)
    end
else
    Z=('Cannot compute distance: deg_Increment must be 2 or larger!');
    disp(Z)
end
end