Bad accuracy in IntrEllipse2Ellipse2
darkstarx opened this issue · 13 comments
There is a problem with accuracy which leads to wrong results.
Here is an example. We have one circle and one ellipse. The ellipse intersects the circle in two points from the left and in one point from the right,
Ellipse2<double> ellipses[2];
ellipses[0].center = { 100.0, 100.0 };
ellipses[0].extent = { 100.0, 100.0 };
ellipses[1].center = { 75.0, 100.0 };
ellipses[1].extent = { 125.0, 50.0 };
FIQuery<double, Ellipse2<double>, Ellipse2<double>> query;
auto area = query.AreaOfIntersection(ellipses[0], ellipses[1]);
auto result = area.findResult;
if (result.intersect) {
std::cout << "!!!=== " << (result.intersect ? "Intersect" : "No intersection") << std::endl;
std::cout << "!!!=== Number of points: " << result.numPoints << std::endl;
std::cout << "!!!=== Area: " << area.area << std::endl;
if (result.numPoints < 5) {
for (int i = 0; i < result.numPoints; ++i) {
const auto &p = result.points[i];
std::cout << "!!!=== " << p[0] << ", " << p[1] << std::endl;
}
}
}
But the intersector give 4 points:
!!!=== Intersect
!!!=== Number of points: 4
!!!=== Area: 16598.173372451663
!!!=== 199.99999999999997, 100.00000238418579
!!!=== 199.99999999999997, 99.999997615814209
!!!=== 9.5238095238095184, 142.59177099999599
!!!=== 9.5238095238095184, 57.408229000004013
In another case, when the ellipse is oriented vertically, it gives only 2 points instead of 3 points:
ellipses[0].center = { 100.0, 100.0 };
ellipses[0].extent = { 100.0, 100.0 };
ellipses[1].center = { 100.0, 125.0 };
ellipses[1].extent = { 50.0, 125.0 };
!!!=== Intersect
!!!=== Number of points: 2
!!!=== Area: 16598.173372451656
!!!=== 57.408229000003651, 190.47619047619355
!!!=== 142.5917709999963, 190.47619047619355
I can understand if the problem shows up on fractional numbers, say if we are trying to position the ellipse using trigonometric functions. But these are strong integer numbers and the correct result must obviously consist of 3 points.
Could you suggest any solution to this problem, please?
And one more question.. Are there any unit tests for that? They might help at this problem, I think.
The problem is fundamental in computing. You are expecting results as if real-valued arithmetic is used. However, on a computer you are using floating-point arithmetic. Floating-point arithmetic has rounding errors. When the two ellipses are tangent at a point of intersection, rounding errors can lead to the presence of two intersection points that are very close to each other. This is what happens in your first example. Rounding errors can lead to no intersection points. This is what happens in your second example. The rounding errors make it appear as if the ellipses have been slightly perturbed.
My implementation has a step that converts the center-axes-extents representation to coefficients of a quadratic equation Q(x,y) = 0. This implicit equation defines the ellipse. The conversion has (very small) rounding errors because of using floating-point arithmetic. If you step into the query with a debugger, set a breakpoint on line 454 of IntrEllipse2Ellipse2.h. The member variables mA and mB are those coefficients. Notice that mA = {1000,-200,-200,1,0}, and the rounding errors have not caused the computed result to differ from the theoretical result. However, notice that mB = {8400,-23.999999999999996,-200,0.15999999999999998,0}. The rounding errors have caused the computed result to differ from the theoretical result {8400,-24,-200,0.16,0}.
After the coefficients are computed, I reduce the problem to computing the roots of a quartic polynomial of 1 variable. The logic is quite detailed. I assume you read the PDF Intersection of Ellipses? This document describes the root finding. The implementation also has some rounding errors. In the end, if area computation is what you want, the computed area has some rounding errors, but the fact that the computed intersection points do not agree with the theoretical intersection points does not lead to catastrophic errors in the computed area. See the code below. The theoretical results were computed using Mathematica.
#include <Mathematics/IntrEllipse2Ellipse2.h>
#include <iostream>
using namespace gte;
int main()
{
Ellipse2<double> ellipses[2];
FIQuery<double, Ellipse2<double>, Ellipse2<double>> query{};
FIQuery<double, Ellipse2<double>, Ellipse2<double>>::AreaResult result{};
ellipses[0].center = { 100.0, 100.0 };
ellipses[0].extent = { 100.0, 100.0 };
ellipses[1].center = { 75.0, 100.0 };
ellipses[1].extent = { 125.0, 50.0 };
result = query.AreaOfIntersection(ellipses[0], ellipses[1]);
// Computed solutions are
// (x0,y0) = (199.99999999999997, 100.00000238418579)
// (x1,y1) = (199.99999999999997, 99.999997615814209)
// (x2,y2) = (9.5238095238095184, 142.59177099999599)
// (x3,y3) = (9.5238095238095184, 57.408229000004013)
// area = 16598.173372451663
//
// Theoretical solutions are the right-hand side of the first
// equalities. The computations of these also use floating-point
// arithmetic and the results are the right-hand side of the
// second equalities
// (x0,y0)
// = (200, 100)
// (x2,y2)
// = (200, 2100 + 400 * sqrt(5)) / 21
// = (9.5238095238095237, 142.59177099999602)
// (x3,y3)
// = (200, 2100 - 400 * sqrt(5)) / 21
// = (9.5238095238095237, 57.408229000004006)
// area
// = -(2500/21) (4 Sqrt[5]-105 ArcCot[Sqrt[5]/4]-168 ArcTan[1/(2 Sqrt[5])])
// = 16598.173372451660
ellipses[0].center = { 100.0, 100.0 };
ellipses[0].extent = { 100.0, 100.0 };
ellipses[1].center = { 100.0, 125.0 };
ellipses[1].extent = { 50.0, 125.0 };
result = query.AreaOfIntersection(ellipses[0], ellipses[1]);
// Computed solutions are
// (x0,y0) = (57.408229000003651, 190.47619047619355)
// (x1,y1) = (142.59177099999630, 190.47619047619355)
// area = 16598.173372451656
//
// Theoretical solutions are the right-hand side of the first
// equalities. The computations of these also use floating-point
// arithmetic and the results are the right-hand side of the
// second equalities
// (x0,y0)
// = (2100 - 400 * sqrt(5), 200) / 21
// = (57.408229000004006, 9.5238095238095237)
// (x1,y1)
// = (2100 + 400 * sqrt(5), 200) / 21
// = (142.59177099999602, 9.5238095238095237)
// (x2,y2)
// = (200, 100)
// area
// = -(2500/21) (4 Sqrt[5]-105 ArcCot[Sqrt[5]/4]-168 ArcTan[1/(2 Sqrt[5])])
// = 16598.173372451663
return 0;
}
It is possible to classify exactly the type of intersection point (transverse or tangential), but you must use exact rational arithmetic to compute the coefficients of the two quadratic polynomials that define the ellipses: P0(x,y) = 0, P1(x,y) = 0. Eliminate the y-variable to obtain a quartic polynomial P(x) = 0, also done using exact rational arithmetic. The discriminant of the quartic polynomial is used for classification. See Low-Degree Polynomial Roots.
Regarding my unit tests, many of them involve using Mathematica (including the ellipse-ellipse intersection algorithm), so I doubt I can help you with that.
Ok, got it. Thanks a lot for your detailed response!
I've read this paper and tried change double to Rational and I've got much more precise result with correct decisions about number of intersection points. Unfortunately I had to make some changes to the code where the <cmath> functions (like std::sqrt, std::atan2 etc.) are used, in order to make the code compilable (may be I forgot some defines). But it works in general, thanks for the great work!
You do not need to change the std::sqrt, std::atan2 functions. At the top of your source file add the line
#include <Mathematics/ArbitraryPrecision.h>
The namespace problem has to do with the 2-pass system the compiler uses. Having this included header first avoids the problem.
@davideberly could you help me one more time? I noticed oversizing the Rational number for rotated ellipses:
typedef BSRational<UIntegerAP32> Rational;
Ellipse2<Rational> ellipses[2];
ellipses[0].center = { 170.0, 144.0 };
ellipses[0].extent = { 64.0, 43.0 };
ellipses[0].axis[0] = { 0.97, 0.26 };
ellipses[0].axis[1] = { -0.26, 0.97 };
ellipses[1].center = { 180.0, 94.0 };
ellipses[1].extent = { 72.0, 29.0 };
ellipses[1].axis[0] = { 0.81, -0.59 };
ellipses[1].axis[1] = { 0.59, 0.81 };
FIQuery<Rational, Ellipse2<Rational>, Ellipse2<Rational>> query;
auto area = query.AreaOfIntersection(ellipses[0], ellipses[1]); // HERE the calculation is very slow
void Mul(UInteger const& n0, UInteger const& n1) has to perform many iterations due to the very large number of operand bits.
Is there any way to use fixed precision arithmetic but with rounding the value instead of assertion "N not large enough to store number of bits." like UIntegerFP32 does? Or may be you can suggest any other idea.
Thanks a lot in advance!
Running this on an Intel i9-10900 in Release x64 (without a debugger attached) required 15.236 seconds. The maximum number of words (the maximum N value) is 86384. Using Intel VTune, the main bottleneck as expected is Mul. The call stacks show that most of the time spent in Mul occurs in two functions, RootsPolynomial::SolveDepressed{Cubic,Quartic}. VTune measured 12.791 seconds in Mul. A secondary bottleneck is std::vector<uint32_t>::operator[] called in Mul. VTune measured 2.471 seconds. Sub measured at 0.016 seconds, resizing std::vector<uint32_t> measured at 0.016 seconds, and ShiftLeft measured at 0.016s.
Using Rational, the area computation is not exact because of the root finding, although the accuracy is improved over a floating-point type for the template parameter T.
I need to think about the shortest path to improving the performance for the area computation. The rational root finding allows you to correctly classify the number of intersection points of the 2 ellipses. I want to preserve this, but once I know the number and type (transverse intersection or tangential intersection), it might be possible to switch to 'double' and use an iterative algorithm to approximate the intersection point. I can examine the ensuing computations and apply rounding of BSRational numbers to a user-specified precision. I already have code for this. The rounding is explicit--you decide when to round. This is different from rounding after each add, sub, mul or div (see item (2) in my list below). Specifically, you can execute multiple BSRational operations followed by a rounding operation.
I will try these ideas and then follow-up with another comment in this thread.
=====
Long-term items on my Infinite TODO List are the following with (4) having the largest development time.
(1) Hack in the Karatsuba algorithm to obtain a somewhat faster Mul.
(2) Implement classes similar to BSNumber and BSRational with user-specified precision and have the operations round each time, effectively doing what an IEEE floating-point system would do for a user-specified number of bits.
(3) Implement a class using the GNU Multiprecision Arithmetic Library for integer operations, effectively giving you another option after UIntegerAP32 or UIntegerFP32.
(4) Replace the O(n^2) Mul by the Schönhage-Strassen algorithm. I have read through some computer algebra papers and books about this topic, but they tend to have a lot of mathematics and not much information on implementation issues. (Who would have ever thought I would complain about too much mathematics...)
Reopening the issue so I can more easily track it.
Thank you for the quick response.
It would be totally great if you find any solution to optimize the algorithm, even if it would be just switching between Rational and double during computation.
Now I decided to check if the rotation is zero (or multiple of 90 degree). If so, then I use Rational, otherwise I use double because it's hardly possible to specify precise position of the rotated ellipses so that they have tangential intersection point (unless one of the ellipses is circle with the same center). I mean it's not so obvious how many intersection points should be for rotated ellipses, and such results seem to be suitable.
Of course, this is a temporary solution, but I have no other idea at the moment. I'll be waiting for algo optimisation.
Anyway, thanks for the great job!
Sorry to take so long on this. Contracting and travel always seem to interfere with maintenance of my own code.
I modified the area-computing code for the intersection of 2 ellipses. I now convert the rational-valued coefficients of the quadratic and quartic polynomials to double precision and then re-convert to rational. For your test example, the compute time on my development computer was reduced from 15 seconds to about 0.02 seconds.
Feel free to re-open the issue if you run into other problems with the code.
Sorry to take so long on this. Contracting and travel always seem to interfere with maintenance of my own code.
I modified the area-computing code for the intersection of 2 ellipses. I now convert the rational-valued coefficients of the quadratic and quartic polynomials to double precision and then re-convert to rational. For your test example, the compute time on my development computer was reduced from 15 seconds to about 0.02 seconds.
Feel free to re-open the issue if you run into other problems with the code.
Awesome, thanks a lot!
I'll check my code and give you feedback asap.
Hello. Sorry for the delay, was busy with another project.
Unfortunately, it's not working for me.
ellipses[0].center = { 100.0, 100.0 };
ellipses[0].extent = { 100.0, 100.0 };
ellipses[0].axis[0] = { 1.0, 0.0 };
ellipses[0].axis[1] = { 0.0, 1.0 };
ellipses[1].center = { 100.0, 100.0 };
ellipses[1].extent = { 100.0, 50.0 };
ellipses[1].axis[0] = { 0.985, 0.174 }; // 10 degrees
ellipses[1].axis[1] = { -0.174, 0.985 };Now it returns
0 (1.51922, 82.6352)
1 (1.51923, 82.6352)
2 (198.481, 117.365)
3 (198.481, 117.365)
instead of two points, and my tests fail.
So, I have to stay on the previous version, it's more accurate in most cases. And the problem of performance is solved by choosing the type of number with next conditions:
const bool sameCenter = offset2DGetX(center1) == offset2DGetX(center2)
&& offset2DGetY(center1) == offset2DGetY(center2);
const bool useRational = sameCenter
? (xAxis1X == 0 || xAxis1Y == 0 || yAxis1X == 0 || yAxis1Y == 0)
|| (xAxis2X == 0 || xAxis2Y == 0 || yAxis2X == 0 || yAxis2Y == 0)
: (xAxis1X == 0 || xAxis1Y == 0 || yAxis1X == 0 || yAxis1Y == 0)
&& (xAxis2X == 0 || xAxis2Y == 0 || yAxis2X == 0 || yAxis2Y == 0);
if (useRational) {
applyResult<Rational>();
} else {
applyResult<Real>();
}
In other words, I use Rational only if one of the ellipses is rotated on the angle multiple of 90 degrees, or if both of them are (in the case when they have different centers). In that case there is no performance problem, so I can safely use Rational. In other cases I can use double 'cause the error doesn't seem obvious.
Thanks.
I have rolled back my changes to IntrEllipse2Ellipse2. Until I wrap GMP with GTE classes for rational arithmetic, the multiplication performance will be slow when the number of words is large. I plan to revisit the algorithm in case there is a geometric way to determine the classification of intersections (transverse or tangential) using a smaller number of words.