davideberly / GeometricTools

A collection of source code for computing in the fields of mathematics, geometry, graphics, image analysis and physics.
Boost Software License 1.0
1.14k stars 214 forks source link

New distance queries #66

Open radevgit opened 1 year ago

radevgit commented 1 year ago

Hi, Do you have any plans to implement distance queries for other geometric objects like: SEG-ARC, ARC-ARC ?

davideberly commented 1 year ago

The distance between 2 segments in any dimension is implemented in DistSegmentSegment.h.

I am currently implementing distance algorithms between line/ray/segment and circle/arc. The line/ray/segment-circle code is finished, including unit testing. I will post this shortly. I need to implement line/ray/segment-arc. This should take me another 1 or 2 days.

owai1980 commented 1 year ago

Waw! This is a very good news! Thank you! 👍

Le lun. 31 juil. 2023 à 15:34, David Eberly @.***> a écrit :

The distance between 2 segments in any dimension is implemented in DistSegmentSegment.h.

I am currently implementing distance algorithms between line/ray/segment and circle/arc. The line/ray/segment-circle code is finished, including unit testing. I will post this shortly. I need to implement line/ray/segment-arc. This should take me another 1 or 2 days.

— Reply to this email directly, view it on GitHub https://github.com/davideberly/GeometricTools/issues/66#issuecomment-1658383734, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG3BAQN3FCBO7KDKTTJLVA3XS6X5LANCNFSM6AAAAAA247QVHU . You are receiving this because you are subscribed to this thread.Message ID: @.***>

radevgit commented 1 year ago

Thanks David, I checked before the DistSegmentSegment.h, and I was surprised how complex is the algorithm. My naive thinking was, that is should be simple than intersection algorithm. Thanks again for the fast response.

davideberly commented 1 year ago

The document Robust Computation of Distance Between Line Segmants has a couple of examples that show the straightforward algorithm, when computed using floating-point arithmetic, can have gross errors because of subtractive cancellation. In particular, see Example 2. If you use the straightforward algorithm computing with rational arithmetic, all is well.

The file I mentioned has a function ComputeRobust(...), where I attempt to avoid the problem with subtractive cancellation. The idea is reasonable in the floating-point era when you do not have a fused-multiply-add instruction. That instruction computes fma(x,y,z) = x*y+z with 1 rounding step. Without such an instruction, w = x*y has 1 rounding step and w+z has 1 rounding step for a total of 2 rounding steps, which generally leads to more rounding error.

With the fma instruction, you can "fix" the subtractive cancellation in expressions x*y-z*w. The implementation is in my Math.h file, functions RobustDOP. The idea is due to Professor William Kahan On the Cost of Floating-Point Computation Without Extra-Precise Arithmetic and succinctly mentioned by Matt Pharr Accurate Differences of Products with Kahan’s Algorithm. At some time I will revisit the segment-segment distance code and modify it to use this paradigm.

radevgit commented 1 year ago

Thanks for the references. I learned a lot. Regarding the fused-multiply-add instruction, it should be available std::fma. Also, regarding such numerical issues, I noted somewhere in the code the use of dot(v0, v1) and sqrt for the length of a vector. There are claims that hypot uses robust algorithm for the same. It looks from some comments, that the hypot also use fma operations.

davideberly commented 1 year ago

The IEEE Standard for Floating-Point Arithmetic 754-2019 includes the requirement that implementations provide fusedMultiplyAdd. For hypotenuse in 2D (squared length of vectors), you can compute xx+yy using the same idea I mentioned. Implementation is RobustSOP (Robust Sum Of Products) in Math.h. The replacement of DOP and SOP will take some time in my code.

mgamito commented 1 year ago

If I may add something to this discussion, I found at work that although std:fma is always supported by the standard, it may not necessarily map to a machine instruction on older architectures and needs to be emulated by software. This turns out to be woefully slow and definitely slower than just writing y = a*x + b. I found this while working with NanoVdb that uses fused-multiply-adds and we had to create a local patched version that removes those functions, so that the older machines on our render farm could stay efficient.

On Tue, Aug 1, 2023 at 4:34 PM David Eberly @.***> wrote:

The IEEE Standard for Floating-Point Arithmetic 754-2019 includes the requirement that implementations provide fusedMultiplyAdd. For hypotenuse in 2D (squared length of vectors), you can compute xx+yy using the same idea I mentioned. Implementation is RobustSOP (Robust Sum Of Products) in Math.h. The replacement of DOP and SOP will take some time in my code.

— Reply to this email directly, view it on GitHub https://github.com/davideberly/GeometricTools/issues/66#issuecomment-1660569533, or unsubscribe https://github.com/notifications/unsubscribe-auth/AGCZNWNWH5INAS4W7HZ2O4LXTEOXVANCNFSM6AAAAAA247QVHU . You are receiving this because you are subscribed to this thread.Message ID: @.***>

davideberly commented 1 year ago

I appreciate your input, thank you.

Software implementations of floating-point functions have that curse. To follow the standard, the logic can be quite complicated and expensive to compute without hardware support. My focus lately has been on robust computations, much of this using rational arithmetic (which effectively is floating-point implemented in software), so I suppose much of my code is "woefully slow" when using rational support.

I added to my short-term TODO list the ability to "discard" FMA operations using a preprocessor define GTE_DISCARD_FMA. If not defined, std::fma(x,y,z) is used. If defined, x*y+z is computed instead with 2 floating-point operations. Longer term I could add a FMA software implementation using the ideas of BSNumber<UIntegerFP32> for N sufficiently large and then add a preprocessor define GTE_USE_SOFTWARE_FMA. But of course this will be slow.

Given you are using modern NVIDIA technology and a render farm, I would have thought a primary goal is to replace old hardware with new hardware. It seems strange that the GPU support of NanoVDB accesses GPU hardware (that has fma), but the underlying CPU is too old for fma. However, if you have a large number of old computers, I can see why hardware upgrade is not desireable.

radevgit commented 1 year ago

For hypotenuse in 2D (squared length of vectors), you can compute x_x+y_y using the same idea I mentioned. Implementation is RobustSOP (Robust Sum Of Products) in Math.h. The replacement of DOP and SOP will take some time in my code.

Here, I found the actual implementation in glibc, that seems to use some additional tricks: hypot

davideberly commented 1 year ago

I have posted 2D distance queries for line-circle, ray-circle, segment-circle, line-arc, ray-arc, segment-arc. I have not yet posted arc-arc queries (need to think about this one a bit).

davideberly commented 1 year ago

I pressed the wrong button, did not mean to close the issue. Reopening it. I am currently working on the 2D distance queries for circle-circle, circle-arc, and arc-arc.

owai1980 commented 1 year ago

I guess you will need less than 5 minutes for the "circle - circle" distance 🤣

Le mar. 29 août 2023 à 20:01, David Eberly @.***> a écrit :

Reopened #66 https://github.com/davideberly/GeometricTools/issues/66.

— Reply to this email directly, view it on GitHub https://github.com/davideberly/GeometricTools/issues/66#event-10224298856, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG3BAQMTT6GXOAZZVERARQTXXYVA5ANCNFSM6AAAAAA247QVHU . You are receiving this because you commented.Message ID: @.***>

tmatthey commented 4 months ago

I did an implementation of seg-arc and arc-arc. I calculate first the line/circle-circle distance(s), then remove the distance(s) not on of seg or cricle. Since you can assume for circles that the points are on the circle you can check if the arc's angle is equal to the sum angles of endpoints with the point, same for point on seg with the distance. Then calculate all distances between the end points seg/arc-arc and take the minimum. I can provide some tests if needed.

davideberly commented 3 months ago

Sorry for taking so long to respond. Sure, send me tests if you like. I already have seg-arc distance code. Does this not work for you?

tmatthey commented 3 months ago

Line3Circle3and Circle3Circle3works and cannot point on any issues. Well done! I was referring to examples for Segment3Arc3and Arc3Arc3. I can provide you with my examples, input and distance.

davideberly commented 3 months ago

Ah. I lost the context of the dimension (3 in your case) and saw that I had distance code for 2D segment/arc. Technical support questions (github issues and private emails) have been consuming much of my time, especially when I have to revise the PDF documentation. I finally cleared out that queue and can finish up the 3D segment/arc code I had started locally. Yes, any examples you want to send me can be incorporated into the unit tests. Thank you.

tmatthey commented 3 months ago

Good. I'll send you a list of Line/Circle3 - Circle3 with start & end points and the min distance as C++.

davideberly commented 3 months ago

Sorry, I am confused. I thought you were going to provide unit-test examples for distances for Segment3-Arc3 and Arc3-Arc3? I already have unit-tested implementations for Line3-Circle3 and Circle3-Circle3.

tmatthey commented 3 months ago

My tests are seg/arc - arc in 3d. I was looking for arc3, found only arc2. So the examples are seg / arc ones🙂

tmatthey commented 3 months ago

I found 130 tests together for Segement3-Arc3 and Arc3-Arc3 distance queries.