g.gui.iclass: replace removed dist_point_to_segment()

[nilason]

This replaces dist_point_to_segment() from matplotlib, which was disabled for matplotlib 3.1.0+.

Addresses: #461.

To test the algorithm [updated]:

#!/usr/bin/env python3

import numpy as np

# from https://matplotlib.org/3.0.0/_modules/matplotlib/mlab.html#dist 
def dist(x, y):
    """
    Return the distance between two points.
    """
    d = x-y
    return np.sqrt(np.dot(d, d))

# from https://matplotlib.org/3.0.0/_modules/matplotlib/mlab.html#dist_point_to_segment
def dist_point_to_segment(p, s0, s1):
    """
    Get the distance of a point to a segment.

      *p*, *s0*, *s1* are *xy* sequences

    This algorithm from
    http://geomalgorithms.com/a02-_lines.html
    """
    p = np.asarray(p, float)
    s0 = np.asarray(s0, float)
    s1 = np.asarray(s1, float)
    v = s1 - s0
    w = p - s0

    c1 = np.dot(w, v)
    if c1 <= 0:
        return dist(p, s0)

    c2 = np.dot(v, v)
    if c2 <= c1:
        return dist(p, s1)

    b = c1 / c2
    pb = s0 + b * v
    return dist(p, pb)

# -------------------------------------------------------------------

def dist_jts(a, b):
    a_x, a_y = a
    b_x, b_y = b

    dx = a_x - b_x
    dy = a_y - b_y

    return np.sqrt(dx * dx + dy * dy)

# function ported from jts https://github.com/locationtech/jts/blob/e5e52b64b5be06ad73932a49a7303a04d89d923a/modules/core/src/main/java/org/locationtech/jts/algorithm/Distance.java#L147
def dist_point_to_segment_jts(p, s0, s1):
    p_x, p_y = p
    s0_x, s0_y = s0
    s1_x, s1_y = s1

    if s0_x == s1_x and s0_y == s1_y:
        return dist_jts(p, s1)

    len2 = (s1_x - s0_x) * (s1_x - s0_x) + (s1_y - s0_y) * (s1_y - s0_y)
    
    r = ((p_x - s0_x) * (s1_x - s0_x) + (p_y - s0_y) * (s1_y - s0_y)) / len2

    if r <= 0.0:
        return dist_jts(p, s0)

    if r >= 1.0:
        return dist_jts(p, s1)

    s = ((s0_y - p_y) * (s1_x - s0_x) - (s0_x - p_x) * (s1_y - s0_y)) / len2

    return np.fabs(s) * np.sqrt(len2)

# -------------------------------------------------------------------

# new function based on https://stackoverflow.com/questions/39840030/distance-between-point-and-a-line-from-two-points/39840218
def dist_point_to_segment_new1(p, s0, s1):
    p = np.asarray(p, float)
    s0 = np.asarray(s0, float)
    s1 = np.asarray(s1, float)

    if (s0 == s1).all():
        d = p - s0
        return np.sqrt(np.dot(d, d))

    d = np.cross(s1 - s0, p - s0) / np.linalg.norm(s1 - s0)

    return np.fabs(d)

# -------------------------------------------------------------------

# new function based on http://www.fundza.com/vectors/point2line/index.html
def dist_point_to_segment_new2(p, s0, s1):
    p = np.asarray(p, float)
    s0 = np.asarray(s0, float)
    s1 = np.asarray(s1, float)

    if (s0 == s1).all():
        d = p - s0
        return np.sqrt(np.dot(d, d))

    line_vec = s1 - s0
    p_vec = p - s0
    line_len = np.sqrt(np.dot(line_vec, line_vec))
    line_unitvec = (line_vec[0]/line_len, line_vec[1]/line_len)
    p_vec_scaled = (p_vec[0] * 1.0/line_len, p_vec[1] * 1.0/line_len)
    t = np.dot(line_unitvec, p_vec_scaled)
    if t < 0.0:
        t = 0.0
    elif t > 1.0:
        t = 1.0

    nearest = (line_vec[0] * t, line_vec[1] * t)
    d = p_vec - nearest
    dist = np.sqrt(np.dot(d, d))

    return dist

# ===================================================================

def main():
    p = -5, 40
    s0 = 0, 0
    s1 = 0, 10

    d = dist_point_to_segment(p, s0, s1)
    print("org  d: %s" %d)

    d2 = dist_point_to_segment_jts(p, s0, s1)
    print("jts  d: %s" % d2)

    d3 = dist_point_to_segment_new1(p, s0, s1)
    print("new  d: %s" %d3)

    d4 = dist_point_to_segment_new2(p, s0, s1)
    print("new2 d: %s" % d4)

if __name__ == "__main__":
    main()

[petrasovaa]

Thank you! Could you please fix the whitespace (add space between operators), the flake8 checks are not enforcing it now, but probably in the future.

[nilason]

updated, hope I got it right :)

[petrasovaa]

Thanks for the test, just be careful, at the end you are printing the same value in both cases.

[petrasovaa]

It can return negative numbers, so you need np.abs(...), that's mentioned also in the SO answer.

[nilason]

Thanks for the test, just be careful, at the end you are printing the same value in both cases.

Whoa, what a blunder! Thanks for noting!

[nilason]

Now, having the test file corrected I did some testing of some other algorithms (and updated the code in first comment). It turned out that the SO solution wasn't as reliable as it was short and neat.

• dist_point_to_segment the original from matplotlib

• dist_point_to_segment_jts port from JTS (also used in GEOS)

• dist_point_to_segment_new1 the, somewhat extended, SO solution
  with present implementation, this is not an acceptable algorithm

• dist_point_to_segment_new2 vectorbased solution by Malcolm Kesson

Any thoughts?

[petrasovaa]

The short one is just solving a different problem (point from line, vs point from segment). I suggest using the algorithm from matplotlib, since it was there before, unless you found some issues with it. I am not entirely sure how the licensing works, but there is pseudo code source linked to it too, so I wouldn't consider it a problem.

[nilason]

I have no strong opinion on which method, matplotlib is good for me. I guess a reference of origin would suffice.

I couldn't find in code base any other use of this function... shall I put it in gui/wxpython/iscatt/plots.py or is there a more general-use file that might be a better place?

[petrasovaa]

Not really, I would suggest creating a new file iscatt/utils.py.

[nilason]

Updated as you suggested, please take a look.

[petrasovaa]

Thank you for getting this done!