compas-dev / compas

Core packages of the COMPAS framework.
https://compas.dev/compas/
MIT License
309 stars 105 forks source link

Line closest_point potential bug #1226

Closed petrasvestartas closed 6 months ago

petrasvestartas commented 9 months ago

I think the line.closest_point method gives bad results for finding the closest parameter on a line.

Example:

Given a line: Line(Point(x=-6.0, y=0.0, z=0.0), Point(x=7.5, y=0.0, z=0.0)) And a point:Point(x=-5.0, y=0.0, z=0.0) The result is 1.0

But it should be: 0.07407407407407407

After looking at OpenNurbs C++ implementation, this is an equivalent in Python:

        def closest_point(line, point):
            t = 0.0 

            v = line.end - line.start
            vov = v.length**2

            if vov > 0.0:
                if (point - line.start).length**2 <= (point - line.end).length**2:
                    t = v.dot(point - line.start) / vov
                else:
                    t = 1.0 + v.dot(point - line.end) / vov
            else:
                t = 0.0 # degenerate line

            return line.point_at(t), t

this is original one in compas:

        a = self.start
        vector = point - a
        direction = self.direction
        t = vector.dot(direction)
        c = a + direction * t
        if return_parameter:
            return c, t
        return c
tomvanmele commented 9 months ago

i would double check your test because the result seems correct to me...

from compas.geometry import Point, Line
from compas_view2.app import App

line = Line(Point(-6.0, 0.0, 0.0), Point(7.5, 0.0, 0.0))
point = Point(-5, 0, 0)
closest, t = line.closest_point(point, True)

print(t)

viewer = App()
viewer.add(line)
viewer.add(closest)

viewer.run()
Screenshot 2023-12-15 at 20 32 20
tomvanmele commented 9 months ago

t is expressed wrt to the direction vector, which is a unit vector. therefore 1.0 seems to make sense, since Point(5, 0, 0) is Point(-6, 0, 0) + Vector(1, 0, 0) * 1.0...

tomvanmele commented 9 months ago

the difference with the C++ implementation lies mostly in the division by the squared length of the vector of the line, but that is implicitly included in the Python implementation because of the use of the unit direction vector.

petrasvestartas commented 9 months ago

I see, it is not normalized parameter, but is a scale value of line.direction.

I thought that the t parameter can be directly used in point_on_line = line.point_at(t) method. This was the confusion.

Which means that instead, it has to be used like this: point_on_line = line.start+line.direction*t

In this case it is not a bug, but intended behavior. I guess I was too much used to Rhino implementation.

tomvanmele commented 9 months ago

no, you're right. it should be possible to combine the two. it was not clear to me from your question/statements that this was the problem...

tomvanmele commented 9 months ago

the domain in also defined as 0.0 => start and 1.0 => end so we should indeed update the implementation.

if don't know why the second "if" is necessary though. t = v.dot(point - line.start) / vov and t = 1.0 + v.dot(point - line.end) / vov seem equivalent to me...

start = self.start
vector = point - start
direction = self.vector
t = vector.dot(direction) / self.length**2
closest = start + direction * t
if return_parameter:
    return closest, t
return closest
tomvanmele commented 9 months ago

i will add a degeneracy check as part of the tolerance module update

tomvanmele commented 9 months ago

if we don't care about the lookups, can be simplified to

vector = point - self.start
t = vector.dot(self.vector) / self.length**2
closest = self.start + self.vector * t
if return_parameter:
    return closest, t
return closest
petrasvestartas commented 9 months ago

if don't know why the second "if" is necessary though. t = v.dot(point - line.start) / vov and t = 1.0 + v.dot(point - line.end) / vov seem equivalent to me...

I checked different scenarios, I agree, there is no need to do that because the result is always the same. The dot product projects one vector to the other, using the same starting point, so it is the same. Also, there is never a need to know to which end the point is closer.

Thank you for changing to the normalized result: t = vector.dot(self.vector) / self.length**2 I think it is more clear when using normalization.

jf--- commented 6 months ago

@petrasvestartas || @petrasvestartas this seems resolved?

tomvanmele commented 6 months ago

yes indeed