Perpendicular on a line segment from a given point
Asked Answered
T

12

36

I want to calculate a point on a given line that is perpendicular from a given point.

I have a line segment AB and have a point C outside line segment. I want to calculate a point D on AB such that CD is perpendicular to AB.

Find point D

I have to find point D.

It quite similar to this, but I want to consider to Z coordinate also as it does not show up correctly in 3D space.

Tirewoman answered 24/4, 2012 at 15:23 Comment(3)
In the absence of a programming question, this would be better on Mathematics (where it is almost certainly already duplicated)Uniliteral
#1812049Campball
It would be good to specify what language you want this in.Verger
C
38

Proof: Point D is on a line CD perpendicular to AB, and of course D belongs to AB. Write down the Dot product of the two vectors CD.AB = 0, and express the fact D belongs to AB as D=A+t(B-A).

We end up with 3 equations:

 Dx=Ax+t(Bx-Ax)
 Dy=Ay+t(By-Ay)
(Dx-Cx)(Bx-Ax)+(Dy-Cy)(By-Ay)=0

Subtitute the first two equations in the third one gives:

(Ax+t(Bx-Ax)-Cx)(Bx-Ax)+(Ay+t(By-Ay)-Cy)(By-Ay)=0

Distributing to solve for t gives:

(Ax-Cx)(Bx-Ax)+t(Bx-Ax)(Bx-Ax)+(Ay-Cy)(By-Ay)+t(By-Ay)(By-Ay)=0

which gives:

t= -[(Ax-Cx)(Bx-Ax)+(Ay-Cy)(By-Ay)]/[(Bx-Ax)^2+(By-Ay)^2]

getting rid of the negative signs:

t=[(Cx-Ax)(Bx-Ax)+(Cy-Ay)(By-Ay)]/[(Bx-Ax)^2+(By-Ay)^2]

Once you have t, you can figure out the coordinates for D from the first two equations.

 Dx=Ax+t(Bx-Ax)
 Dy=Ay+t(By-Ay)
Campball answered 3/3, 2013 at 15:56 Comment(3)
This ignores the Z component of the original question.Antibiosis
If you want to limit the found point D to the line between A and B limit t for 0 >= t <= 1.Klemens
If this not work for u.. Try this solution #1812049Extortion
H
27
function getSpPoint(A,B,C){
    var x1=A.x, y1=A.y, x2=B.x, y2=B.y, x3=C.x, y3=C.y;
    var px = x2-x1, py = y2-y1, dAB = px*px + py*py;
    var u = ((x3 - x1) * px + (y3 - y1) * py) / dAB;
    var x = x1 + u * px, y = y1 + u * py;
    return {x:x, y:y}; //this is D
}

question

Headstall answered 19/9, 2012 at 17:1 Comment(5)
Would be nice to see a little explanation of what you did.Monsour
Hi. I have point a,b and d and i want to calculate x and y coordinates of point c. Can i use your function? If no please tell me the modifications.Anschluss
@MeanCoder point c has countless possibilities if only know point a,b,d.Headstall
What if A,B, C are GPS coordinates (lat,lon)?Fulfillment
It doesn't matter if they're GPS or not. They're still coordinates. You will have to account for cases where A and B lie on different sides of the 0 coordinate.Maniemanifest
W
9

There is a simple closed form solution for this (requiring no loops or approximations) using the vector dot product.

Imagine your points as vectors where point A is at the origin (0,0) and all other points are referenced from it (you can easily transform your points to this reference frame by subtracting point A from every point).

In this reference frame point D is simply the vector projection of point C on the vector B which is expressed as:

// Per wikipedia this is more efficient than the standard (A . Bhat) * Bhat
Vector projection = Vector.DotProduct(A, B) / Vector.DotProduct(B, B) * B

The result vector can be transformed back to the original coordinate system by adding point A to it.

Westerman answered 24/4, 2012 at 16:1 Comment(0)
B
6

A point on line AB can be parametrized by:

M(x)=A+x*(B-A), for x real.

You want D=M(x) such that DC and AB are orthogonal:

dot(B-A,C-M(x))=0.

That is: dot(B-A,C-A-x*(B-A))=0, or dot(B-A,C-A)=x*dot(B-A,B-A), giving:

x=dot(B-A,C-A)/dot(B-A,B-A) which is defined unless A=B.

Byssus answered 24/4, 2012 at 18:1 Comment(0)
R
2

What you are trying to do is called vector projection

Rosy answered 24/4, 2012 at 15:57 Comment(0)
T
1

Here i have converted answered code from "cuixiping" to matlab code.

function Pr=getSpPoint(Line,Point)
% getSpPoint(): find Perpendicular on a line segment from a given point
x1=Line(1,1);
y1=Line(1,2);
x2=Line(2,1);
y2=Line(2,1);
x3=Point(1,1);
y3=Point(1,2);

px = x2-x1;
py = y2-y1;
dAB = px*px + py*py;

u = ((x3 - x1) * px + (y3 - y1) * py) / dAB;
x = x1 + u * px;
y = y1 + u * py;

Pr=[x,y];

end
Tristis answered 9/3, 2016 at 11:50 Comment(0)
T
0

I didn't see this answer offered, but Ron Warholic had a great suggestion with the Vector Projection. ACD is merely a right triangle.

  1. Create the vector AC i.e (Cx - Ax, Cy - Ay)
  2. Create the Vector AB i.e (Bx - Ax, By - Ay)
  3. Dot product of AC and AB is equal to the cosine of the angle between the vectors. i.e cos(theta) = ACx*ABx + ACy*ABy.
  4. Length of a vector is sqrt(x*x + y*y)
  5. Length of AD = cos(theta)*length(AC)
  6. Normalize AB i.e (ABx/length(AB), ABy/length(AB))
  7. D = A + NAB*length(AD)
Tuner answered 21/12, 2016 at 5:12 Comment(0)
E
0

For anyone who might need this in C# I'll save you some time:

double Ax = ;
double Ay = ;
double Az = ;
    
double Bx = ;
double By = ;
double Bz = ;
    
double Cx = ;
double Cy = ;
double Cz = ; 
    
double t = ((Cx - Ax) * (Bx - Ax) + (Cy - Ay) * (By - Ay)) / (Math.Pow(Bx - Ax, 2) + Math.Pow(By - Ay, 2));
    
double Dx = Ax + t*(Bx - Ax);
double Dy = Ay + t*(By - Ay);
Enesco answered 19/3, 2021 at 7:38 Comment(0)
C
0

Here is another python implementation without using a for loop. It works for any number of points and any number of line segments. Given p_array as a set of points, and x_array , y_array as continues line segments or a polyline.

This uses the equation Y = mX + n and considering that the m factor for a perpendicular line segment is -1/m.


      import numpy as np

      def ortoSegmentPoint(self, p_array, x_array, y_array):
            """
    
            :param p_array: np.array([[ 718898.941  9677612.901 ], [ 718888.8227 9677718.305 ], [ 719033.0528 9677770.692 ]])
            :param y_array: np.array([9677656.39934991 9677720.27550726 9677754.79])
            :param x_array: np.array([718895.88881594 718938.61392781 718961.46])
            :return: [POINT, LINE] indexes where point is orthogonal to line segment
            """
            # PENDIENTE "m" de la recta, y = mx + n
            m_array = np.divide(y_array[1:] - y_array[:-1], x_array[1:] - x_array[:-1])
            # PENDIENTE INVERTIDA, 1/m
            inv_m_array = np.divide(1, m_array)
            # VALOR "n", y = mx + n
            n_array = y_array[:-1] - x_array[:-1] * m_array
            # VALOR "n_orto" PARA LA RECTA PERPENDICULAR
            n_orto_array = np.array(p_array[:, 1]).reshape(len(p_array), 1) + inv_m_array * np.array(p_array[:, 0]).reshape(len(p_array), 1)
            # PUNTOS DONDE SE INTERSECTAN DE FORMA PERPENDICULAR
            x_intersec_array = np.divide(n_orto_array - n_array, m_array + inv_m_array)
            y_intersec_array = m_array * x_intersec_array + n_array
            # LISTAR COORDENADAS EN PARES
            x_coord = np.array([x_array[:-1], x_array[1:]]).T
            y_coord = np.array([y_array[:-1], y_array[1:]]).T
            # FILAS: NUMERO DE PUNTOS, COLUMNAS: NUMERO DE TRAMOS
            maskX = np.where(np.logical_and(x_intersec_array < np.max(x_coord, axis=1), x_intersec_array > np.min(x_coord, axis=1)), True, False)
            maskY = np.where(np.logical_and(y_intersec_array < np.max(y_coord, axis=1), y_intersec_array > np.min(y_coord, axis=1)), True, False)
            mask = maskY * maskX
            return np.argwhere(mask == True)

Creaky answered 30/12, 2021 at 3:33 Comment(0)
G
0

As Ron Warholic and Nicolas Repiquet answered, this can be solved using vector projection. For completeness I'll add a python/numpy implementation of this here in case it saves anyone else some time:

import numpy as np

# Define some test data that you can solve for directly.
first_point = np.array([4, 4])
second_point = np.array([8, 4])
target_point = np.array([6, 6])

# Expected answer
expected_point = np.array([6, 4])

# Create vector for first point on line to perpendicular point.
point_vector = target_point - first_point
# Create vector for first point and second point on line.
line_vector = second_point - first_point

# Create the projection vector that will define the position of the resultant point with respect to the first point.
projection_vector = (np.dot(point_vector, line_vector) / np.dot(line_vector, line_vector)) * line_vector

# Alternative method proposed in another answer if for whatever reason you prefer to use this.
_projection_vector = (np.dot(point_vector, line_vector) / np.linalg.norm(line_vector)**2) * line_vector

# Add the projection vector to the first point
projected_point = first_point + projection_vector

# Test
(projected_point == expected_point).all()
Gadget answered 3/8, 2022 at 16:41 Comment(0)
T
-1

Since you're not stating which language you're using, I'll give you a generic answer:

Just have a loop passing through all the points in your AB segment, "draw a segment" to C from them, get the distance from C to D and from A to D, and apply pithagoras theorem. If AD^2 + CD^2 = AC^2, then you've found your point.

Also, you can optimize your code by starting the loop by the shortest side (considering AD and BD sides), since you'll find that point earlier.

Thoracoplasty answered 24/4, 2012 at 15:31 Comment(0)
H
-1

Here is a python implementation based on Corey Ogburn's answer from this thread.
It projects the point q onto the line segment defined by p1 and p2 resulting in the point r.
It will return null if r falls outside of line segment:

def is_point_on_line(p1, p2, q):

    if (p1[0] == p2[0]) and (p1[1] == p2[1]):
        p1[0] -= 0.00001

    U = ((q[0] - p1[0]) * (p2[0] - p1[0])) + ((q[1] - p1[1]) * (p2[1] - p1[1]))
    Udenom = math.pow(p2[0] - p1[0], 2) + math.pow(p2[1] - p1[1], 2)
    U /= Udenom

    r = [0, 0]
    r[0] = p1[0] + (U * (p2[0] - p1[0]))
    r[1] = p1[1] + (U * (p2[1] - p1[1]))

    minx = min(p1[0], p2[0])
    maxx = max(p1[0], p2[0])
    miny = min(p1[1], p2[1])
    maxy = max(p1[1], p2[1])

    is_valid = (minx <= r[0] <= maxx) and (miny <= r[1] <= maxy)

    if is_valid:
        return r
    else:
        return None
Helotism answered 3/2, 2017 at 0:47 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.