How to Get The Coordinates of The Point on a Line That Has The Smallest Distance from Another Point

How to get the coordinates of the point on a line that has the smallest distance from another point

In a mathematical point of view you can :

• first find the equation of the line :

      y1 = a1x1+b1 
      a1 = (y2-y1) / (x2-x1)
      b1 = y1-a1*x1
  

• Then calculate the gradient of the second line knowing :

      a1 * a2 = -1 <-> 
      a2 = -1/a1
  

• with a2 you can find the value of b for the second equation :

      y3 = a2*x3 + b2 <-> 
      b2 = y3 - a2*x3
  

• Finally calculate the intersection of the 2 lines :

      xi = (b2-b1) / (a1-a2)
      y = a1*xi + b1
  

Then it's quite straightforward to bring that to swift :

typealias Line = (gradient:CGFloat, intercept:CGFloat)

func getLineEquation(point1:CGPoint, point2:CGPoint) -> Line {
    guard point1.x != point2.x else {
        if(point1.y != point2.y)
        {
            print("Vertical line : x = \(point1.x)")
        }
        return (gradient: .nan, intercept: .nan)
    }
    let gradient = (point2.y - point1.y)/(point2.x-point1.x)
    let intercept = point1.y - gradient*point1.x
    return (gradient: gradient, intercept: intercept)
}

func getPerpendicularGradient(gradient:CGFloat) -> CGFloat
{
    guard gradient != 0 else {
        print("horizontal line, the perpendicilar line is vertical")
        return .nan
    }
    return -1/gradient
}

func getIntercept(forPoint point:CGPoint, withGradient gradient:CGFloat) -> CGFloat
{
    return point.y - gradient * point.x
}

func getIntersectionPoint(line1:Line, line2:Line)-> CGPoint
{
    guard line1.gradient != line2.gradient else {return CGPoint(x: CGFloat.nan, y: CGFloat.nan)}
    let x = (line2.intercept - line1.intercept)/(line1.gradient-line2.gradient)
    return CGPoint(x:x, y: line1.gradient*x + line1.intercept)
}

func getClosestIntersectionPoint(forLine line:Line, point:CGPoint) -> CGPoint
{
    let line2Gradient = getPerpendicularGradient(gradient:line.gradient)
    let line2 = (
        gradient: line2Gradient,
        intercept: getIntercept(forPoint: point, withGradient: line2Gradient))
    
    return getIntersectionPoint(line1:line, line2:line2)
}

func getClosestIntersectionPoint(forLinePoint1 linePoint1:CGPoint, linePoint2:CGPoint, point:CGPoint) -> CGPoint
{
    return getClosestIntersectionPoint(
        forLine:getLineEquation(point1: linePoint1, point2: linePoint2),
        point:point)
}

Shortest distance between a point and a line segment

Eli, the code you've settled on is incorrect. A point near the line on which the segment lies but far off one end of the segment would be incorrectly judged near the segment. Update: The incorrect answer mentioned is no longer the accepted one.

Here's some correct code, in C++. It presumes a class 2D-vector class vec2 {float x,y;}, essentially, with operators to add, subract, scale, etc, and a distance and dot product function (i.e. x1 x2 + y1 y2).

float minimum_distance(vec2 v, vec2 w, vec2 p) {
  // Return minimum distance between line segment vw and point p
  const float l2 = length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  if (l2 == 0.0) return distance(p, v);   // v == w case
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line. 
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  // We clamp t from [0,1] to handle points outside the segment vw.
  const float t = max(0, min(1, dot(p - v, w - v) / l2));
  const vec2 projection = v + t * (w - v);  // Projection falls on the segment
  return distance(p, projection);
}

EDIT: I needed a Javascript implementation, so here it is, with no dependencies (or comments, but it's a direct port of the above). Points are represented as objects with x and y attributes.

function sqr(x) { return x * x }
function dist2(v, w) { return sqr(v.x - w.x) + sqr(v.y - w.y) }
function distToSegmentSquared(p, v, w) {
  var l2 = dist2(v, w);
  if (l2 == 0) return dist2(p, v);
  var t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
  t = Math.max(0, Math.min(1, t));
  return dist2(p, { x: v.x + t * (w.x - v.x),
                    y: v.y + t * (w.y - v.y) });
}
function distToSegment(p, v, w) { return Math.sqrt(distToSegmentSquared(p, v, w)); }

EDIT 2: I needed a Java version, but more important, I needed it in 3d instead of 2d.

float dist_to_segment_squared(float px, float py, float pz, float lx1, float ly1, float lz1, float lx2, float ly2, float lz2) {
  float line_dist = dist_sq(lx1, ly1, lz1, lx2, ly2, lz2);
  if (line_dist == 0) return dist_sq(px, py, pz, lx1, ly1, lz1);
  float t = ((px - lx1) * (lx2 - lx1) + (py - ly1) * (ly2 - ly1) + (pz - lz1) * (lz2 - lz1)) / line_dist;
  t = constrain(t, 0, 1);
  return dist_sq(px, py, pz, lx1 + t * (lx2 - lx1), ly1 + t * (ly2 - ly1), lz1 + t * (lz2 - lz1));
}

Here, in the function parameters, <px,py,pz> is the point in question and the line segment has the endpoints <lx1,ly1,lz1> and <lx2,ly2,lz2>. The function dist_sq (which is assumed to exist) finds the square of the distance between two points.

Find the shortest distance between a point and line segments (not line)

Here is the answer. This code belongs to Malcolm Kesson, the source is here. I provided it before with just link itself and it was deleted by the moderator. I assume that the reason for that is because of not providing the code (as an answer).

import math

def dot(v,w):
    x,y,z = v
    X,Y,Z = w
    return x*X + y*Y + z*Z

def length(v):
    x,y,z = v
    return math.sqrt(x*x + y*y + z*z)

def vector(b,e):
    x,y,z = b
    X,Y,Z = e
    return (X-x, Y-y, Z-z)

def unit(v):
    x,y,z = v
    mag = length(v)
    return (x/mag, y/mag, z/mag)

def distance(p0,p1):
    return length(vector(p0,p1))

def scale(v,sc):
    x,y,z = v
    return (x * sc, y * sc, z * sc)

def add(v,w):
    x,y,z = v
    X,Y,Z = w
    return (x+X, y+Y, z+Z)

# Given a line with coordinates 'start' and 'end' and the
# coordinates of a point 'pnt' the proc returns the shortest 
# distance from pnt to the line and the coordinates of the 
# nearest point on the line.
#
# 1  Convert the line segment to a vector ('line_vec').
# 2  Create a vector connecting start to pnt ('pnt_vec').
# 3  Find the length of the line vector ('line_len').
# 4  Convert line_vec to a unit vector ('line_unitvec').
# 5  Scale pnt_vec by line_len ('pnt_vec_scaled').
# 6  Get the dot product of line_unitvec and pnt_vec_scaled ('t').
# 7  Ensure t is in the range 0 to 1.
# 8  Use t to get the nearest location on the line to the end
#    of vector pnt_vec_scaled ('nearest').
# 9  Calculate the distance from nearest to pnt_vec_scaled.
# 10 Translate nearest back to the start/end line. 
# Malcolm Kesson 16 Dec 2012

def pnt2line(pnt, start, end):
    line_vec = vector(start, end)
    pnt_vec = vector(start, pnt)
    line_len = length(line_vec)
    line_unitvec = unit(line_vec)
    pnt_vec_scaled = scale(pnt_vec, 1.0/line_len)
    t = dot(line_unitvec, pnt_vec_scaled)    
    if t < 0.0:
        t = 0.0
    elif t > 1.0:
        t = 1.0
    nearest = scale(line_vec, t)
    dist = distance(nearest, pnt_vec)
    nearest = add(nearest, start)
    return (dist, nearest)

How to calculate the (smallest) distance between a set of lines and a point?

I am assuming you are working with the spherical model of the Earth. Then I suppose what you call "lines" are in fact arc-segments of great circles (the straight lines of spherical geometry). In other words, you have one point p1 on the surface of the sphere with lat-lon coordinates [lat1, lon1] and another point p2 on the sphere with lat-lon coordinates [lat2, lon2]. What is considered "the straight line" on the sphere that passes through p1 and p2 is the unique circle (called great circle) obtained by the intersection of the sphere with the plane that passes through the center of the sphere and the two points p1 and p2. Then, what you call a "line", is the smaller of the two circular arcs from this great circle, bounded by the two points p1 and p2.

The distance you would like to calculate (in radians) could be the distance from a third point p with lat-lon coordinates [lat, lon] to the arc-segment determined by the two points p1 and p2. The said distance should be the arc-length of the arc from the great circle passing through p and perpendicular to the great circle of p1 and p2. This perpendicular great circle is the one determined by the intersection of the sphere with the plane perpendicular to the plane of the great circle of p1 and p2 and passing through the point p and the center of the sphere. If the intersection of the perpendicular great circle with the great-circle arc p1 p2 is a point h inside the arc-segment p1 p2, then the length of the great circle arc p h is the sought distance. If, however, h is outside the arc p1 p2 the the sought distance is either p p1 or p p2 whichever is smaller.

Here is some Matlab code that calculates the shortest distance between a point and an arc-interval:

lat_lon = [lat, lon];
lat_lon1 = [lat1, lon1];
lat_lon2 = [lat2, lon2];

function dist = dist_point_2_road(lat_lon, lat_lon1, lat_lon2)

   % you may have to convert lat-long angles from degrees to radians 

   % First, convert from lat-long coordinates to 3D coordinates of points on the unit    
   % sphere. Since Earth's radius cancels out in our computations, we simply assume it 
   % is R = 1 
   lat = lat_lon(1);
   lon = lat_lon(2);

   lat1 = lat_lon1(1);
   lon1 = lat_lon1(2);

   lat2 = lat_lon2(1);
   lon2 = lat_lon2(2);

   p1 = [ cosd(lat1)*cosd(lon1),  cosd(lat1)*sind(lon1),  sind(lat1) ]; %cosd = cos(degrees)
   p2 = [ cosd(lat2)*cosd(lon2),  cosd(lat2)*sind(lon2),  sind(lat2) ]; %sind = sin(degrees)
   p = [ cosd(lat)*cosd(lon),  cosd(lat)*sind(lon),  sind(lat) ];

   % n12 is the unit vector perpendicular to the plane of the great circle 
   % determined by the points p1 and p2  
   n12 = cross(p1, p2);
   n12 = n12 / sqrt(dot(n12, n12));
   sin_of_dist = dot(p, n12); % sine of the angle that equals arc-distance 
                              % from point p to the great arc p1 p2  

   dist = pi/2 - acos(abs(sin_of_dist)); % acos = arccos, abs() = absolute value
          % dist is the shortest distance in radians from p to the 
          % great circle determined by the points p1 and p2

   n1 = cross(p1, p); 
   n1 = n1 / sqrt(dot(n1, n1));
   % unit normal vector perpendicular to the great-arc determined by p and p1

   n2 = cross(p, p2);
   n2 = n1 / sqrt(dot(n2, n2));
   % unit normal vector perpendicular to the great-arc determined by p and p2

   if dot(n12, n1) < 0 % if the angle of spherical triangle p p1 p2 at vertex p1 is obtuse 
      dist = acos(dot(p, p1)); % the shortest distance is p p1 
   elseif dot(n12, n2) < 0 % if the angle of spherical triangle p p1 p2 at vertex p2 is obtuse 
      dist = acos(dot(p, p2)); % the shortest distance is p p2 
   end

   % the function returns the appropriate dist as output 

end

You can iterate this for the sequence of arc-intervals that form a road and select the smallest distance to an arc-interval.

According to this computation, the distance of the point to the first "vector 1" is
0.0000615970599633145 radians and the distance to the second "vector 2" is
0.00162840939265068 radians. The point is closest to a point inside vector 1, but for vector 2, it is closest to to one of the ends of the arc-interval.

Edit. Now, if one wants to use Euculdiean (flat) approximation, ignoring Earths curvature, one may needs to convert lat-lon coordinates into Euclidean flat approximation coordinates. To avoid any map specific coordinates, one may be tempted to plot latitude against longitude coordinates. That may be ok around the equator, but the closer to the poles one gets, the more innaccurate these coordinates gets at representing distance data. This is because closer to the poles, distance along a fixed latitude is much shorter than distance along a fixed longitude. That is why we need to correct this discrepancy. This is done by using the Riemannian metric on the sphere in lat-long coordinates, or simply by looking at the 3D geometry of latitude and longitude circles near a given point on the sphere.

lat_lon = [lat, lon];
lat_lon1 = [lat1, lon1];
lat_lon2 = [lat2, lon2];

%center of approximate Euclidean coordinate system is point p 
% with lat_long coordinates and the scaling coefficient of longitude, 
% which equalizes longitude and latitude distance at point p, is

a = cosd(lat_long(1));  

function  x = convert_2_Eucl(lat_long1, lat_long, a)
   x = [lat_long1(1) - lat_long(1),  a*(lat_long1(2) - lat_long(2))];
end

% convert from lat-long to approximate Euclidean coordinates
x1 = convert_2_Eucl(lat_long1, lat_long, a);
x2 = convert_2_Eucl(lat_long2, lat_long, a);

function dist = dist_point_2_road(x1, x2)

   dist = dot(x1, x1) * dot(x2 - x1, x2 - x1) - dot(x1, x2 - x1)^2 ;
   dist = sqrt( dist / ( dot(x2 - x1, x2 - x1)^2) );
   % dist is the distance from the point p, which has Eucl coordinates [0,0] 
   % to the straight Euclidean interval x1 x2 representing the interval p1 p2

   if dot(x1, x2 - x1) > 0
      dist = sqrt( dot(x1, x1) );
   elseif dot(x2, x1 - x2) > 0
      dist = sqrt( dot(x2, x2) );
   end

end

Remark: the latter function calculates distance, but it might be equally convenient to just calculate dist^2, avoiding the calculation of the square root sqrt in order to speed up the performance. Measuring with respect to dist^2 should work just as well.

You choose which function you want, the spherical one or approximately Euclidean. The latter is probably faster. You can choose to remove the square root and calculate distance squared to make things even faster.

I wrote this in a hurry, so there might be some inaccuracies.

Shortest distance between points algorithm

http://en.wikipedia.org/wiki/Closest_pair_of_points

The problem can be solved in O(n log n) time using the recursive divide and conquer approach, e.g., as follows:

• Sort points along the x-coordinate
• Split the set of points into two equal-sized subsets by a vertical line x = xmid
• Solve the problem recursively in the left and right subsets. This will give the left-side and right-side minimal distances dLmin and dRmin respectively.
• Find the minimal distance dLRmin among the pair of points in which one point lies on the left of the dividing vertical and the second point lies to the right.
• The final answer is the minimum among dLmin, dRmin, and dLRmin.

Finding the smallest distance in a set of points from the origin

I'll work off of your original code. Firstly, you want to compute the norm of all of the points and store them as individual elements in an array. Your current code isn't doing that and is overwriting the variable L which is a single value at each iteration of the loop.

You'll want to make L an array and store the norms at each iteration of the loop. Once you do this, you'll want to find the location as well as the minimum distance itself. That can be done with one call to min where the first output gives you the minimum distance and the second output gives you the location of the minimum. You can use the second output to slice into your S array to retrieve the actual point.

Last but not least, you need to define S first before calling this function. You are defining S inside the function and that will probably give you unintended results if you want to change the input into this function at each invocation. Therefore, define S first, then call the function:

 S= [-6.8667,  -44.7967; 
    -38.0136, -35.5284; 
     14.4552, -27.1413; 
     8.4996,   31.7294; 
    -17.2183,  28.4815; 
    -37.5100,  14.1941; 
    -4.2664,  -24.4428; 
    -18.6655,  26.9427; 
    -15.8828,  18.0170; 
     17.8440, -22.9164];

function [dist,koor] = bonus4(S)             
    %// New - Create an array to store the distances
    L = zeros(size(S,1), 1);

    %// Change to iterate over number of rows
    for i=1:size(S,1) 
        L(i)=norm(S(i, :)); %// Change
    end

    [dist,ind] = min(L); %// Find the minimum distance
    koor = S(ind,:); %// Get the actual point
end

Or, make sure you save the above function in a file called bonus4.m, then do this in the Octave command prompt:

octave:1> S= [-6.8667,  -44.7967; 
>         -38.0136, -35.5284; 
>          14.4552, -27.1413; 
>          8.4996,   31.7294; 
>         -17.2183,  28.4815; 
>         -37.5100,  14.1941; 
>         -4.2664,  -24.4428; 
>         -18.6655,  26.9427; 
>         -15.8828,  18.0170; 
>          17.8440, -22.9164];
octave:2> [dist,koor] = bonus4(S);

Though this code works, I'll debate that it's slow as you're using a for loop. A faster way would be to do this completely vectorized. Because using norm for matrices is different than with vectors, you'll have to compute the distance yourself. Because you are measuring the distance from the origin, you can simply square each of the columns individually then add the columns of each row.

Therefore, you can just do this:

S= [-6.8667,  -44.7967; 
    -38.0136, -35.5284; 
     14.4552, -27.1413; 
     8.4996,   31.7294; 
    -17.2183,  28.4815; 
    -37.5100,  14.1941; 
    -4.2664,  -24.4428; 
    -18.6655,  26.9427; 
    -15.8828,  18.0170; 
     17.8440, -22.9164]; 

function [dist,koor] = bonus4(S)     
    %// New - Computes the norm of each point
    L = sqrt(sum(S.^2, 2));

    [dist,ind] = min(L); %// Find the minimum distance
    koor = S(ind,:); %// Get the actual point

end

The function sum can be used to sum over a dimension independently. As such, by doing S.^2, you are squaring each term in the points matrix, then by using sum with the second parameter as 2, you are summing over all of the columns for each row. Taking the square root of this result computes the distance of each point to the origin, exactly the way the for loop functions. However, this (at least to me) is more readable and I daresay faster for larger sizes of points.

---

Related Topics