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

IntrSegment2Segment2 precison problem #88

Closed owai1980 closed 6 months ago

owai1980 commented 6 months ago

Hello,

I have a problem while calculating the intersection of 2 "Segment2<double>"

i know it's a precision problem, and maybe the answer is "don't use the double version" ... but i still want to explain:

(they make a "T" shape like this:) segment0 = (100, 1900) , (900,1900) segment1 = (550, 601) , (550.00000000000034,1900)

as you can see, i have a small precision error on the "x1" coordinate of segment1

when i call the Result() function in IntrSegment2Segment2.h (line 215), the result is "intersect = false" (line 244), probably because of the GetCenteredForm() before...

I understand that the precison of the "double" is difficult to manage, but i though there was a "tolerance range" somewhere in GeometricTools (?) maybe i'm wrong...

edit : the "tolerance" parameter is "epsilon", used in Arc2.h ... maybe we should use it in Segment2 too... what do you think?

thank you

Johan

davideberly commented 6 months ago

Adding "epsilon" in the algorithms is a rabbit hole I choose not to go down. Infrequently, there might be a need and I use an epsilon, but not as a general solution to dealing with floating-point error. Is the epsilon relative or absolute? How small (or large) must it be? If the epsilon is absolute, the magnitude of epsilon depends on the input data, so the user must decide what epsilon to choose based on an understanding of the characteristics of the input.

Regarding your example, my choice of GetCenteredForm is unfortunate. I tried to formulate the segments in the manner of an oriented bounding box (center+direction+extent), but yes, the rounding errors can be a problem. I am eliminating the GetCenteredForm in GTL. For code that is not time critical (for example, not intended for real-time applications and not used in a very large number of queries), I prefer the rational number approach. If you use BSRational as the rational type, you get "intersect = true". However, you can use the alternative member function "Exact". For your example, you get "intersect = true".

    using Rational = BSRational<UIntegerAP32>;

    Segment2<double> segment0{}, segment1{};
    segment0.p[0] = { 100.0, 1900.0 };
    segment0.p[1] = { 900.0, 1900.0 };
    segment1.p[0] = { 550.0, 601.0 };
    segment1.p[1] = { 550.00000000000034, 1900.0 };
    FIQuery<double, Segment2<double>, Segment2<double>> query{};
    auto result = query(segment0, segment1);
    // result.intersect = false
    // result.numIntersections = 0
    // result.segment0Parameter = { 0.0, 0.0 }
    // result.segment1Parameter = { 0.0, 0.0 }
    // result.point[0] = { 0.0, 0.0 }
    // result.point[1] = { 0.0, 0.0 }

    result = query.Exact(segment0, segment1);
    // result.intersect = true
    // result.numIntersections = 1
    // result.segment0Parameter = { 0.56250000000000044, 0.56250000000000044 }
    // result.segment1Parameter = { 1.0, 1.0 }
    // result.point[0] = { 550.00000000000034, 1900.0000000000000 }
    // result.point[1] = { 550.00000000000034, 1900.0000000000000 }

    Segment2<Rational> rsegment0{}, rsegment1{};
    rsegment0.p[0] = { 100.0, 1900.0 };
    rsegment0.p[1] = { 900.0, 1900.0 };
    rsegment1.p[0] = { 550.0, 601.0 };
    rsegment1.p[1] = { 550.00000000000034, 1900.0 };
    FIQuery<Rational, Segment2<Rational>, Segment2<Rational>> rquery{};
    auto rresult = rquery(rsegment0, rsegment1);

    FIQuery<double, Segment2<double>, Segment2<double>>::Result convert{};
    convert.intersect = rresult.intersect;
    convert.numIntersections = rresult.numIntersections;
    convert.segment0Parameter[0] = rresult.segment0Parameter[0];
    convert.segment0Parameter[1] = rresult.segment0Parameter[1];
    convert.segment1Parameter[0] = rresult.segment1Parameter[0];
    convert.segment1Parameter[1] = rresult.segment1Parameter[1];
    convert.point[0][0] = rresult.point[0][0];
    convert.point[0][1] = rresult.point[0][1];
    convert.point[1][0] = rresult.point[1][0];
    convert.point[1][1] = rresult.point[1][1];
    // convert.intersect = true
    // convert.numIntersections = 1
    // convert.segment0Parameter = { 50.000000000000341, 50.000000000000341 }
    // convert.segment1Parameter = { 649.50000000000000, 649.50000000000000 }
    // convert.point[0] = { 550.00000000000034, 1900.0000000000000 }
    // convert.point[1] = { 550.00000000000034, 1900.0000000000000 }
owai1980 commented 6 months ago

Once again, thank you for your fast and complete answer.

I didn't know there was a "exact()" function! I will test it tomorrow morning! 👍

And i agree with you about the "rabbit hole"...

Thank you!

Johan

Le mar. 26 mars 2024 à 18:07, David Eberly @.***> a écrit :

Adding "epsilon" in the algorithms is a rabbit hole I choose not to go down. Infrequently, there might be a need and I use an epsilon, but not as a general solution to dealing with floating-point error. Is the epsilon relative or absolute? How small (or large) must it be? If the epsilon is absolute, the magnitude of epsilon depends on the input data, so the user must decide what epsilon to choose based on an understanding of the characteristics of the input.

Regarding your example, my choice of GetCenteredForm is unfortunate. I tried to formulate the segments in the manner of an oriented bounding box (center+direction+extent), but yes, the rounding errors can be a problem. I am eliminating the GetCenteredForm in GTL. For code that is not time critical (for example, not intended for real-time applications and not used in a very large number of queries), I prefer the rational number approach. If you use BSRational as the rational type, you get "intersect = true". However, you can use the alternative member function "Exact". For your example, you get "intersect = true".

using Rational = BSRational<UIntegerAP32>;

Segment2<double> segment0{}, segment1{};
segment0.p[0] = { 100.0, 1900.0 };
segment0.p[1] = { 900.0, 1900.0 };
segment1.p[0] = { 550.0, 601.0 };
segment1.p[1] = { 550.00000000000034, 1900.0 };
FIQuery<double, Segment2<double>, Segment2<double>> query{};
auto result = query(segment0, segment1);
// result.intersect = false
// result.numIntersections = 0
// result.segment0Parameter = { 0.0, 0.0 }
// result.segment1Parameter = { 0.0, 0.0 }
// result.point[0] = { 0.0, 0.0 }
// result.point[1] = { 0.0, 0.0 }

result = query.Exact(segment0, segment1);
// result.intersect = true
// result.numIntersections = 1
// result.segment0Parameter = { 0.56250000000000044, 0.56250000000000044 }
// result.segment1Parameter = { 1.0, 1.0 }
// result.point[0] = { 550.00000000000034, 1900.0000000000000 }
// result.point[1] = { 550.00000000000034, 1900.0000000000000 }

Segment2<Rational> rsegment0{}, rsegment1{};
rsegment0.p[0] = { 100.0, 1900.0 };
rsegment0.p[1] = { 900.0, 1900.0 };
rsegment1.p[0] = { 550.0, 601.0 };
rsegment1.p[1] = { 550.00000000000034, 1900.0 };
FIQuery<Rational, Segment2<Rational>, Segment2<Rational>> rquery{};
auto rresult = rquery(rsegment0, rsegment1);

FIQuery<double, Segment2<double>, Segment2<double>>::Result convert{};
convert.intersect = rresult.intersect;
convert.numIntersections = rresult.numIntersections;
convert.segment0Parameter[0] = rresult.segment0Parameter[0];
convert.segment0Parameter[1] = rresult.segment0Parameter[1];
convert.segment1Parameter[0] = rresult.segment1Parameter[0];
convert.segment1Parameter[1] = rresult.segment1Parameter[1];
convert.point[0][0] = rresult.point[0][0];
convert.point[0][1] = rresult.point[0][1];
convert.point[1][0] = rresult.point[1][0];
convert.point[1][1] = rresult.point[1][1];
// convert.intersect = true
// convert.numIntersections = 1
// convert.segment0Parameter = { 50.000000000000341, 50.000000000000341 }
// convert.segment1Parameter = { 649.50000000000000, 649.50000000000000 }
// convert.point[0] = { 550.00000000000034, 1900.0000000000000 }
// convert.point[1] = { 550.00000000000034, 1900.0000000000000 }

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

owai1980 commented 6 months ago

I did some tests...

The exact() function does work well for yesterday's example... But doesn't work if the values are different.

Anyway, i will not use the segment2segment2 intersection, I will use the line2line2 intersection then test if the point is ok.

I understand that the double precision is always a problem...

Johan

Le mar. 26 mars 2024 à 18:07, David Eberly @.***> a écrit :

Adding "epsilon" in the algorithms is a rabbit hole I choose not to go down. Infrequently, there might be a need and I use an epsilon, but not as a general solution to dealing with floating-point error. Is the epsilon relative or absolute? How small (or large) must it be? If the epsilon is absolute, the magnitude of epsilon depends on the input data, so the user must decide what epsilon to choose based on an understanding of the characteristics of the input.

Regarding your example, my choice of GetCenteredForm is unfortunate. I tried to formulate the segments in the manner of an oriented bounding box (center+direction+extent), but yes, the rounding errors can be a problem. I am eliminating the GetCenteredForm in GTL. For code that is not time critical (for example, not intended for real-time applications and not used in a very large number of queries), I prefer the rational number approach. If you use BSRational as the rational type, you get "intersect = true". However, you can use the alternative member function "Exact". For your example, you get "intersect = true".

using Rational = BSRational<UIntegerAP32>;

Segment2<double> segment0{}, segment1{};
segment0.p[0] = { 100.0, 1900.0 };
segment0.p[1] = { 900.0, 1900.0 };
segment1.p[0] = { 550.0, 601.0 };
segment1.p[1] = { 550.00000000000034, 1900.0 };
FIQuery<double, Segment2<double>, Segment2<double>> query{};
auto result = query(segment0, segment1);
// result.intersect = false
// result.numIntersections = 0
// result.segment0Parameter = { 0.0, 0.0 }
// result.segment1Parameter = { 0.0, 0.0 }
// result.point[0] = { 0.0, 0.0 }
// result.point[1] = { 0.0, 0.0 }

result = query.Exact(segment0, segment1);
// result.intersect = true
// result.numIntersections = 1
// result.segment0Parameter = { 0.56250000000000044, 0.56250000000000044 }
// result.segment1Parameter = { 1.0, 1.0 }
// result.point[0] = { 550.00000000000034, 1900.0000000000000 }
// result.point[1] = { 550.00000000000034, 1900.0000000000000 }

Segment2<Rational> rsegment0{}, rsegment1{};
rsegment0.p[0] = { 100.0, 1900.0 };
rsegment0.p[1] = { 900.0, 1900.0 };
rsegment1.p[0] = { 550.0, 601.0 };
rsegment1.p[1] = { 550.00000000000034, 1900.0 };
FIQuery<Rational, Segment2<Rational>, Segment2<Rational>> rquery{};
auto rresult = rquery(rsegment0, rsegment1);

FIQuery<double, Segment2<double>, Segment2<double>>::Result convert{};
convert.intersect = rresult.intersect;
convert.numIntersections = rresult.numIntersections;
convert.segment0Parameter[0] = rresult.segment0Parameter[0];
convert.segment0Parameter[1] = rresult.segment0Parameter[1];
convert.segment1Parameter[0] = rresult.segment1Parameter[0];
convert.segment1Parameter[1] = rresult.segment1Parameter[1];
convert.point[0][0] = rresult.point[0][0];
convert.point[0][1] = rresult.point[0][1];
convert.point[1][0] = rresult.point[1][0];
convert.point[1][1] = rresult.point[1][1];
// convert.intersect = true
// convert.numIntersections = 1
// convert.segment0Parameter = { 50.000000000000341, 50.000000000000341 }
// convert.segment1Parameter = { 649.50000000000000, 649.50000000000000 }
// convert.point[0] = { 550.00000000000034, 1900.0000000000000 }
// convert.point[1] = { 550.00000000000034, 1900.0000000000000 }

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

davideberly commented 6 months ago

The Exact function produces the theoretically correct values when using a rational type. Is the rational performance an issue for you?

For line2-line2 intersection, what do you expect the results to be when the lines are "nearly" parallel. I can imagine constructing an example for which the point of intersection has NaN or infinities because the theoretically correct value cannot be represented with the precision of the floating-point type you use. But perhaps for your data, this approach will work well.

You might try DistSegmentSegment.h instead to compute the distance between two line segments. When the segments theoretically intersect (distance = 0), rounding errors can lead to distance being a positive number that is nearly zero. You can use an epsilon test in your application code. Note that DistSegmentSegment has an operator() query that produces the theoretically correct squared distance when using Rational. It also has ComputeRobust which uses a conjugate gradient approach to try to be more robust.

owai1980 commented 6 months ago

Thank you...

Maybe you are right and i could convert my segment2 from double to rationals just for the time of the intersection call. I'm afraid to start big scale changes in the whole app...

I will save all your advices...

Thank you.

Johan

Le mer. 27 mars 2024 à 23:11, David Eberly @.***> a écrit :

The Exact function produces the theoretically correct values when using a rational type. Is the rational performance an issue for you?

For line2-line2 intersection, what do you expect the results to be when the lines are "nearly" parallel. I can imagine constructing an example for which the point of intersection has NaN or infinities because the theoretically correct value cannot be represented with the precision of the floating-point type you use. But perhaps for your data, this approach will work well.

You might try DistSegmentSegment.h instead to compute the distance between two line segments. When the segments theoretically intersect (distance = 0), rounding errors can lead to distance being a positive number that is nearly zero. You can use an epsilon test in your application code. Note that DistSegmentSegment has an operator() query that produces the theoretically correct squared distance when using Rational. It also has ComputeRobust which uses a conjugate gradient approach to try to be more robust.

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

davideberly commented 6 months ago

If you need to re-open this, no worries.