Points and Vectors

I occasionaly write scripts for geometry optimization, therefore I have to wrestle with geometry. Graphically very easy problems become Ye Royale Paine In Ye Buttocks to solve analytically. At first I wrestle with piles of points, lines, determinants, trigonometric functions and whatchamacallits. Then replace all of that with a few points and vectors. Here are a few simple functions I either found online or wrote myself to help me with geometric calculations.

Angle Between Vectors

I found this on Stackoverflow. Since it works, I had no reason to fix it.

import numpy as np
import scipy as sp

def vector(x, y, z):
    return np.array([x, y, z])

def norm(v):
    return np.linalg.norm(v)

def unit_vector(vector):
    return vector / norm(vector)

def angle_between(v1, v2):
    """ Kudos: https://stackoverflow.com/questions/2827393/
    Returns the angle in radians between vectors 'v1' and 'v2':
   
    >>> angle_between((1, 0, 0), (0, 1, 0))
    1.5707963267948966
    >>> angle_between((1, 0, 0), (1, 0, 0))
    0.0
    >>> angle_between((1, 0, 0), (-1, 0, 0))
    3.141592653589793
    """

    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)

    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))

Point Rotation

It’s a generic problem but I can still share the code that worked so far.

def rotate(point, angle, axis='x'):
    """ rotates a point around an axis by specified angle """
    if axis == 'y':
        rot_matrix = np.array([
            [np.cos(angle),  0, np.sin(angle)],
            [0,              1, 0            ],
            [-np.sin(angle), 0, np.cos(angle)],
        ])
    elif axis == 'z':
        rot_matrix = np.array([
            [np.cos(angle), -np.sin(angle), 0],
            [np.sin(angle), np.cos(angle),  0],
            [0,             0,              1],
        ])
    else: # == 'x':
        rot_matrix = np.array([
                [1, 0, 0],
                [0, np.cos(angle), -np.sin(angle)],
                [0, np.sin(angle), np.cos(angle)],
            ])

    return point.dot(rot_matrix)

De-rotate a Point Back to XY-Plane

This is not a projection of the point to XY plane but a rotation around x-axis to that plane. There are two options: calculate the angle and rotate the point by that angle or transform cartesian coordinates to cylindrical and set angle to zero. The latter is the more elegant but the first works just as well (depending on your precision requirements).

def de_rotate(p):
    """ rotate the point p around x-axis so that z-axis=0 """
    # measure angle between normal and point
    normal = np.array([0, 1, 0])
    point = np.array([0, p[1], p[2]])
    phi = angle_between(point, normal)
    
    # rotate the oritinal point around x-axis by measured angle
    rotated = rotate(p, phi)
    return rotated, phi

def de_rotate(p):
    # alternative solution:
    # calculate r in cylindrical coordinate system;
    # also return the angle in cylindrical system
    return np.array([p[0], (p[1]**2 + p[2]**2)**0.5, 0]), np.arctan2(p[2], p[1])

Intersection of Two Lines

This is a 2D-only problem. We have two lines, each is given by two points. The objective is to calculate the intersection between these two lines.

This code is just a transcription of dictation of a recipe from Wikipedia.

def xy_line_intersection(p_1, p_2, p_3, p_4):
    """ p_1 and p_2 define the first line, p_3 and p_4 define the second; 
        return a point of intersection between these two lines in x-y plane

        Kudos: https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
    """
    # only take x and y coordinates
    x1 = p_1[0]
    y1 = p_1[1]

    x2 = p_2[0]
    y2 = p_2[1]

    x3 = p_3[0]
    y3 = p_3[1]

    x4 = p_4[0]
    y4 = p_4[1]

    def det(p1, p2, p3, p4):
        return np.linalg.det(np.array([[p1, p2], [p3, p4]]))

    Dx1 = det(x1, y1, x2, y2)
    Dx2 = det(x1,  1, x2,  1)
    Dx3 = det(x3, y3, x4, y4)
    Dx4 = det(x3,  1, x4,  1)
    Dx5 = Dx2
    Dx6 = det(y1,  1, y2,  1)
    Dx7 = Dx4
    Dx8 = det(y3,  1, y4,  1)

    # x-coordinate
    Px = det(Dx1, Dx2, Dx3, Dx4)/det(Dx5, Dx6, Dx7, Dx8)

    # y-coordinate
    Dy1 = Dx1
    Dy2 = Dx6
    Dy3 = Dx3
    Dy4 = Dx8
    Dy5 = Dx2
    Dy6 = Dx6
    Dy7 = Dx7
    Dy8 = Dx8

    Py = det(Dy1, Dy2, Dy3, Dy4)/det(Dy5, Dy6, Dy7, Dy8)
    return vector(Px, Py, 0)

This can also be done with vectors. Instead of line a vector between two points is created. Then each is scaled so that the new vectors produce the same point. Here’s some ASCII math for clarification (vectors are in bold):

P1(P1x, P1y), P2(P2x, P2y) - points that define first line
P3(P3x, P3y), P4(P4x, P4y) - points that define the second line
P1 + k1*(P2 - P1) = P3 + k2*(P4 - P3)

Now the above equation  contains vectors so it is in fact two equations, one for x-component and one for y.

P1x + k1*(P2x - P1x) = P3x + k2*(P4x - P3x)
P1y + k1*(P2y - P1y) = P3y + k2*(P4y - P3y)

Above is a linear system of two equations and two unknowns. By solving it we get k1 and k2 and in theory, the intersecting point can be calculated with either the left or the right side of the first equation. In practice (my specific case) the left side works (within 1e-16 of the previous method) but right doesn’t (more like within 1e-5 to 1e-3). I am not sure why but it might be due to positions of the points P3 and P4. Any comment would be most welcome.

def xy_line_intersection(p_1, p_2, p_3, p_4):
    # alternative solution with vectors
    A = np.array([
        [p_2[0] - p_1[0], p_4[0] - p_3[0]],
        [p_2[1] - p_1[1], p_4[1] - p_3[1]],
    ])
    
    b = np.array([p_3[0] - p_1[0], p_3[1] - p_1[1]])
    
    k1k2 = np.linalg.solve(A, b)
    k1 = k1k2[0]
    k2 = k1k2[1]
    
    va = vector(
        p_1[0] + k1*(p_2[0] - p_1[0]),
        p_1[1] + k1*(p_2[1] - p_1[1]),
        0
    )
    
    vb = vector(
        p_3[0] + k2*(p_4[0] - p_3[0]),
        p_3[1] + k2*(p_4[1] - p_3[1]),
        0
    )
    
    # print(P-va, P-vb)
    return va

Extend a Vector to Specified y

In this problem we look for the point where the line given by two points meets given y-value. This is a quick hack and doesn’t work with y=constant lines. But works in all other cases.

def extend_to_y(p_1, p_2, y):
    """ return a point that lies on a line defined by p_1 and p_2 and on y=y """
    fk_3 = lambda k: p_1[1] + k*(p_2 - p_1)[1] - y
    k_3 = sp.optimize.newton(fk_3, 0)

    return p_1 + k_3*(p_2 - p_1)

Check If a Point is on a Line

We have three points, A, B, and C. Points A and B define a line and we would like to check if point C lies on that line within certain tolerances.

This problem was the most Paineful in me Buttocks until I switched to vector notation. The recipe is now quite simple (again, mind the bold vectors):

AB = B - A
AC = C - A
k = |AC| / |AB|
C' = A + k*AB

Now C‘ is a new point that we got to using another route. If it is on the same line it should be the same as point C. Within certain tolerance:

|C - C'| < tol ---> the point is on line AB.

And the code:

def is_between(A, B, C, tol):
    # returns True if point C lies on a line between A and B;

    # take an alternative path to point C:
    AB = B - A
    AC = C - A
    k = np.linalg.norm(AC) / np.linalg.norm(AB)
    C__ = A + k*AB

    # it's the same point if the two are close enough
    return np.linalg.norm(C - C__) < tol

If this is of any use to anyone, I will very pleased and even more surprised. If this is the case, I wish you all the best and good luck with your points and vectors.