python - test if two segments are roughly collinear using numpy-CodePudding

I want to test if two segments are roughly collinear 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 answer by Anderson Oliveira.

CodePudding user response：

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.

CodePudding user response：

Cross products are only really defined in 3D. Dot products are defined everywhere. The angle between two vectors is defined as the arccosine 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) < threshold or abs(angle) > np.pi - 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 cthresh = np.cos(threahold) once:

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