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?

CodePudding user response：

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

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.