Test if two segments are roughly collinear (on the same line)

Question

I want to test if two segments are roughly collinear (on the same line) using numpy.cross. I have the coordinates in meters of the segments.

import numpy as np

segment_A_x1 = -8020537.5158307655
segment_A_y1 = 5674541.918222183
segment_A_x2 = -8020547.42095263
segment_A_y2 = 5674500.781350276

segment_B_x1 = -8020556.569040865
segment_B_y1 = 5674462.788207927
segment_B_x2 = -8020594.740831952
segment_B_y2 = 5674328.095911447

a = np.array([[segment_A_x1, segment_A_y1], [segment_A_x2, segment_A_y2]])
b = np.array([[segment_B_x1, segment_B_y1], [segment_B_x2, segment_B_y2]])
crossproduct = np.cross(a, b)

>>>array([7.42783487e+08, 1.65354844e+09])

The crossproduct values are pretty high even if I would say those two segments are roughly collinear. Why?

How can I determine if the segments are colinear with the crossproduct result?

Is there a possibility of using a tolerance in meters to tell if the segments are roughly collinear?

Answer 1

The cross-product of two vector a and b can be defined defined as: ||a|| ||b|| sin(θ) where θ is the angle between the two vectors. The thing is the even if θ is small, the product of the two norms can be huge and can result in a big cross-product. This is your case:

na = np.linalg.norm(a, axis=1)  # array([9824940.10284589, 9824924.42969893])
nb = np.linalg.norm(b, axis=1)  # array([9824909.95439353, 9824863.32407282])
nprod = na * nb                 # array([9.65291518e+13, 9.65285397e+13])
sinθ = crossproduct / nprod     # array([7.69491364e-06, 1.71301508e-05])

Indeed, sinθ is very small. This means θ is very small too (in fact, θ is roughly equal to sinθ value since the derivative of sin(θ) with θ = 1 is 1).

How can I determine if the segments are colinear with the crossproduct result?

You can compute the value np.abs(sinθ) base on the above code, set an arbitrary threshold and check if the value is smaller than the chosen threshold. The best value for the threshold depends of your actual application/input (eg. statistical noise, precision of the input, accuracy of the computations, etc.). If you do not know what to use, you can start with 1e-4 for example.

Is there a possibility of using a tolerance in meters to tell if the segments are roughly collinear?

Note that assuming the input vectors are in meters, then the value of the cross-product is in square meters based on the above formula. Thus, a tolerance in meters does not make much sense. More generally, an setting an absolute tolerance is certainly a bad idea. The above solution works using a tolerance relative to the input vectors.

Answer 2

The problem with your approach is that the cross product value depends on the measurement scale.

Maybe the most intuitive measure of collinearity is the angle between the line segments. Let's calculate it:

import math

def slope(line): 
    """Line slope given two points"""
    p1, p2 = line
    return (p2[1] - p1[1]) / (p2[0] - p1[0])

def angle(s1, s2): 
    """Angle between two lines given their slopes"""
    return math.degrees(math.atan((s2 - s1) / (1 + (s2 * s1))))

ang = angle(slope(b), slope(a))
print('Angle in degrees = ', ang)

Angle in degrees = 2.2845

I made use of an answer by Anderson Oliveira.

Answer 3

Collinearity has two aspects: angle between the segments, and proximity. You can test these aspects separately, and short-circuit the check based on either. I propose a proximity check that removes the need for an angle check entirely. The angle check below is shown for legacy reasons, since it can be useful for other reasons.

Angle

Cross products are only really defined in 3D. Dot products are defined everywhere. The angle between two vectors is defined as the arc-cosine of the dot product of their unit vectors. That applies in 2D as it does in 10D.

You can define two vectors from the origin as

a = np.array([segment_A_x2 - segment_A_x1,
              segment_A_y2 - segment_A_y1])
b = np.array([segment_B_x2 - segment_B_x1,
              segment_B_y2 - segment_B_y1])

a /= np.linalg.norm(a)
b /= np.linalg.norm(b)

angle = np.arccos(a.dot(b))

angle will run between -np.pi and np.pi. Angles close to pi indicate anti-parallel lines, while angles close to zero indicate parallel lines. You can implement the test something like

if abs(angle) < angular_threshold or abs(angle) > np.pi - angular_threshold:
   ...

If you find yourself doing this computation often, you can save some cycles by skipping the arccos entirely. As the angle goes to zero or pi, the dot product goes to +1 or -1 respectively. That means you only need to pre-compute dot_threshold = np.cos(angular_threshold) once:

if abs(a.dot(b)) >= dot_threshold:

Proximity

To test proximity, you need to define how distance is measured. For exactly parallel lines, this is unambiguous: the two lines must be within distance_threshold of each-other.

For not-quite-parallel lines you can do something similar: no point on one segment may be further than distance_threshold of the line of the other segment. With this definition, you can absorb the need for a separate angle computation directly into the calculation.

See the figure below for reference:

You will have to check all four endpoints against the other lines. If you choose to, you can compute the angle from the difference of the distances, and possibly apply a scaling factor to the threshold based on the length of the line segment.

You can compute the distances using dot products as well:

def dist(p, p1, p2):
    s = p2 - p1
    q = p1 + (p - p1).dot(s) / s.dot(s) * s
    return np.linalg.norm(p - q)

a1 = np.array([segment_A_x1, segment_A_y1])
a2 = np.array([segment_A_x2, segment_A_y2])
b1 = np.array([segment_B_x1, segment_B_y1])
b2 = np.array([segment_B_x2, segment_B_y2])

dists = np.array([[dist(a1, b1, b2), dist(b1, a1, a2)],
                  [dist(a2, b1, b2), dist(b2, a1, a2)]])

A more robust version of dist is available in a utility library I made, called haggis. You can use haggis.math.segment_distance as follows:

dists = segment_distance([[a1, b1], [a2, b2]], # From point
                         [b1, a1],             # Segment start
                         [b2, a2],             # Segment end
                         axis=-1,              # Axis containing vectors
                         segment=False)        # Distance to entire line

The inputs broadcast together, and axis applies to the broadcasted shape, so you don't need to repeat the endpoints twice.

The simplest version of the test is to just constrain the distance directly:

if (dists < distance_threshold).all():
    ...

You could account for the angle by scaling by the length of the segment. A segment that is 100x longer can deviate by 100x more from the line of the other segment and still be considered collinear-ish. In this case, you define distance_threshold as the farthest that a unit vector can be from the line of the other vector. The number must be less than one to be meaningful:

scale = np.linalg.norm([a2 - a1, b2 - b1], axis=-1)
if (dists < distance_threshold * scale).all():
    ...

This version presupposes the shape we imposed on dists in both versions of the computation, since scale has as many elements as dists has columns.

More complicated schemes that also account for the distance between the line segments are also possible. Such definitions of proximity are left as an exercise for the reader.

Answer 4

If you would like to see whether two vectors are in almost the same direction in an efficient way, you could also try the following that uses numpy is_close for proximity check, which allows you to further configure what kinds of tolerances to use and how tight it is.

import numpy as np

# define your segement start and end coordinates here
# segment1_start = ...
# segemnt2_start = ...
# segemnt1_end = ...
# segment1_end = ...

# we use the convention of end-start
# modify this to suit your conventions
a = segment1_end - segment1_start
b = segment2_end - segment2_start


def is_parallel_to_axes(a, b):
    vector1 = np.asarray(a)
    vector2 = np.asarray(b)
    vector1_norm = np.linalg.norm(vector1)
    vector2_norm = np.linalg.norm(vector2)
    dot_prod = np.dot(vector1, vector2)


    # Check if the vector is parallel to any of the axes
    is_parallel = np.isclose(dot_prod, vector1_norm * vector2_norm)
    
    return is_parallel