Distance from reference line at right angle

Question

A few weeks back I received a script that, I think, calculates the distance of a point from a reference line (at right angle) and returns the distance along the reference line as well as the distance from the reference line to that point. To be honest, this is what I think it does, however I do not fully understand the script even though it is only 26 lines long.

Since I think the script is the problem why another script, which depends on it, gets unstable I figured it is actually handy to really understand what is going on. I had hoped someone could help me with this.

The script (written in Matlab but I've also got it in python if necessary):

function [sp,np]=locate(s,x,y,xp,yp)
ipp=0;
ns=length(s);
for ip=1:length(xp);
    for i=1:ns-2
        cosa=(x(i+2)-x(i))/(s(i+2)-s(i));
        sina=(y(i+2)-y(i))/(s(i+2)-s(i));
        sproj=s(i)+(xp(ip)-x(i))*cosa+(yp(ip)-y(i))*sina;
        if sproj<s(1)
            ipp=ipp+1;
            sp(ipp)=sproj;
            np(ipp)=-(xp(ip)-x(1))*sina+(yp(ip)-y(1))*cosa;
            break
        elseif sproj>=s(i)&sproj<=s(i+2)
            ipp=ipp+1;
            sp(ipp)=sproj;
            np(ipp)=-(xp(ip)-x(i))*sina+(yp(ip)-y(i))*cosa;
            break
        elseif sproj>s(ns)
            ipp=ipp+1;
            sp(ipp)=sproj;
            np(ipp)=-(xp(ip)-x(ns))*sina+(yp(ip)-y(ns))*cosa;
            break
        end
    end
end

Here s is the distance along the reference line, x and y are the x and y points on the reference line and xp and yp are the coordinates of the points of which the distance to the line (at right angle) as well as the distance along the reference line need to be calculated. The script is called as follows:

dist(1)=0;
for i=2:length(xref);
    dist(i)=dist(i-1)+sqrt((xref(i)-xref(i-1))^2+(yref(i)-yref(i-1))^2);
end

%% Create computational grid
ds=(dist(end)-dist(1))/(ns-1);  % stepsize
s=0:ds:dist(end);               % distance
xr=spline(dist,xref,s);            % x of line grid points
yr=spline(dist,yref,s);            % y of line grid points

%% Compute locations of initial line
[si,ni]=locate(s,xr,yr,xi,yi);
n=interp1(si,ni,s,'linear','extrap');    % distance change at right angle

Is there anyone that can help me by explaining what exactly happens in the first script, cause I can't really get my head around it.

Answer 1

My opinion on the first function depends on whether the points described by arrays x,y lie on a line (as you say) or not (as the second function suggests; in the second script they are obtained from spline interpolation.

If x,y are consecutive points on the line

Then the computation of sina, cosa within the loop is redundant since the angle is always the same. For each point in the arrays xp,yp the function calculates its orthogonal projection onto the line. It returns the position along the line, measured from the first point x(1),y(1), and signed normal distance. All of this happens in a rather inefficient way for this task.

If x,y do not lie on a line

Then the computation of sina, cosa is incorrect. These are meant to be components of a unit vector in the direction from (x(i),y(i)) to (x(i+2),y(i+2)). But the denominator is s(i+2)-s(i), which [looking at the second function] will generally be greater than the straight-line distance from (x(i),y(i)) to (x(i+2),y(i+2)).

The following would make more sense:

cosa = (x(i+2)-x(i))/sqrt((x(i+2)-x(i))^2 + (y(i+2)-y(i))^2);
sina = (y(i+2)-y(i))/sqrt((x(i+2)-x(i))^2 + (y(i+2)-y(i))^2);

The confusion of straight-line distances and along-the-curve distances persists through the rest of the function: the comparisons of sproj with s(i+2) and especially with s(ns) are not geometrically motivated.

Anyway, the function attempts to locate a point near the orthogonal projection of xp,yp onto the curve (which is not unique to begin with); depending on how wiggly the curve is, it may have moderate success or fail badly.

Answer 2

If I understand what you're looking for, this is it.

This function takes in three points, the first two as a vector of x values and a vector of y values (I'll call them p1 and p2) and a third point (p3). The function returns the minimum distance from p3 to the line intersecting p1&p2 and the minimum of the distance from (p1 to the intersection point) or (p2 to the intersection point).

Example:

The two points to define the line are (1,2) and (5,10). The line defining them would be y = 2x. If the third point were (1,0). The perpendicular line is y = -(x-1)/2. The intersection point of the two lines would be (2/5, 1/5). The closest point defining the reference line would be (1,2) and the distance from (2/5, 1/5) to (1,2) is 1.7 and the distance from (2/5, 1/5) to (1,0) is 0.89. The function would return [1.7, 0.89].

Please let me know if this is what you want.

function [ dRef, dInt ] = StackExchange( x, y, xP, yP )
%StackExchange: Calculate line intersecting (x(1), y(1)) & (x(2), y(2))
%finds perpendicular line which intersects (xP, yP) which occurs as
%(IP(1), IP(2)) then calculated the minimum distance from the reference points
%(x(1), y(1)) or (x(2), y(2)) to (IP(1), IP(2))
%and from (IP(1), IP(2)) to (xP, yP)

%Calculate  slope of reference line, the slope of the perpendicular line
%will be -1/m
m = diff(y) ./ diff(x) ;

%IP(1) is the x-coordinate of the intersection point
%This formula is solving the equation
% m * (x - x(1)) + y(1) = -(1/m) * (x - xP) + yP
%for x
IP(1) = (m^2 * x(1) + xP + m * (yP - y(1))) ./ (m^2 + 1);
%IP(2) is the y-coordinate of the intersection point
IP(2) = m * (IP(1) - x(1)) + y(1);

%Minimum distance from the reference points to the intersection point
dRef = sqrt(min((IP(1) - x).^2 + (IP(2) - y).^2));

%Distance from the intersection point to  point of interest
dInt = sqrt(sum(([xP, yP] - IP).^2));

end