[Abandoned] Implementation of the Simpson's rule and the length of an Ellipse(-arc)

[CVH]

Andrew, all,

REllipse::getSimpsonLength(a, b) in REllipse.cpp (Current) is basically the implementation of Calculation of Ellipse Arc Length found online.

Both methods use 20 intervals.
As per Simpson's rule the factor q is 1.0 for the first and last interval, 4.0 when even, otherwise a factor of 2.0 is used.
Minor difference: In QCAD the interval designation is zero based, in the linked reference it starts with 1, odd/even interval factors are swapped.

An important note for the Java applet is Scope : First Quadrant, 0 ≤ x ≤ R
Where R is the length of the Semi-Axis lying on the X-axis aka the Major Radius.

Another important fact is that yi is a positive root, the sum of the factored y values is positive thus the df / 3.0 part must be positive.
Meaning that df must be positive and that a2 > a1 so that df = (a2-a1)/20 is larger than zero.

The use of getSimpsonLength(a, b) is limited to REllipse::getLength() for returning an estimate of the length of the ellipse(-arc).
It is thus important that the implementation in getLength() complies with the above.


Studied the implementation in depth knowing that the length of an Ellipse-arc may fail in some situations (FS#2693 >> FS#2565).
It doesn't comply with the first quadrant limit.
Estimates of a half full ellipse circumference are used, then also over no more than 20 intervals.
It fails with for example: An ellipse starting at 45°, passing 0° in CW sense and ends at (-)180°... =2Pi is the issue here.

Not that it matters a lot to know the exact length of an ellipse-arc but when returning NaN drawing tools like D2 may fail

4 years ago a 'faster' algorithm was implemented for the circumference of a full ellipse. ('faster' is debatable, elegant too )
Not referenced but the Padé 3/3 approximation is said to be very practical and more efficient than the earlier used Ramanujan II approximation.
Padé 3/3 doesn't use a root and the overall accuracy is increased by a 1.5 notch on the log scale. (Review of known formulae)
My 2 main reasons for proposing it.

It is not the limit of what can be achieved.
The Cantrell-Ramanujan approximation (C6) has for example a maximum error of 14.5ppm.
Even better approximations exist at the cost of complexity. (e.g. 152ppb Stan's Library)
(Also see notes in wire Centroids but these calculate a nearly exact full circumference by well suited infinite series expansion.)


Knowing the 'almost' exact circumference of a full ellipse one can easily figure out a half or a quarter ellipse circumference.
Because it is symmetrical in regards with the X and/or Y-axis simply divide for example the Padé 3/3 approximation by 2 or 4.

Then we could chop up an arbitrary ellipse-arc and estimate parts or remainders in ... or mirrored to the first quadrant.
Our best option is then to estimate the more flatter part and subtract that from a quarter ellipse circumference if required.
Already limited to call getSimpsonLength(a, b) at most twice in combination with the 'fast' method for a full ellipse.
It will always be an estimated length without implementing elliptic integrals of the first and second kind.


Next would be to implement all this under QCAD in ECMAScript and report back on the success.
I can then compare the summed length of some arbitrary ellipse-arcs with the nearly exact full circumference from Centroids.
e.g. 3x 1/3 full ellipse-arcs should match with a full ellipse and so on.

Regards,
CVH

[CVH]

All,

I can already share my first observation concerning the estimation of half the full ellipse circumference.
In detail the use of getSimpsonLength(0.0, Pi) or getSimpsonLength(Pi, 2Pi) in REllipse::getSimpsonLength(a, b).
(Second main issue)

Rather good for nearly round ellipses (ratio ≈1.0 or Eccentricity near zero).
The more the ellipse is elongated, flatter, the poorer the result.

2 times a quarter ellipse (getSimpsonLength(0.0, Pi/2)) is a far better estimation for a half ellipse.
In the mid range with for example ratio = 0.5 it is over 150000 times more accurate.
But that rapidly declines to 6.5 times more accurate for a 0.05 ratio or when almost flat.

For rather elongated ellipses the reason must be obvious.
We are applying the Simpson's rule on sections with high to very high curvature near the ellipse major points.
For a half ellipse twice near both major points and that with coarser steps than 2 times the quarter ellipse approach.

Then 1/4 or 1/2 of the Padé 3/3 approximation of a full ellipse is always more accurate.


For the record: Although never attempted to estimate for in REllipse::getSimpsonLength(a, b)
An ellipse length from Pi/2 to 3Pi/2 exhibits similar poor results as those found for 0-Pi or Pi-2Pi.
What supports the statement to avoid the extremities and rather estimate the flatter part if convenient.

Regards,
CVH

[CVH]

All,

I abandoned the Simpson's rule for approximating the length of an Ellipse in general.
Although:
- I first fixed the problems with NaN lengths, see post 49156.
- I figured out a method for a near exact approximation of a complete Ellipse circumference by an adaptive amount of intervals.

The fix is already committed and will be beneficial for tools that rely on a defined length.
In these cases it doesn't matter how correct that is.

Adaptive intervals is no solution for arbitrary ellipse arc's.
The biggest problem for me was comparing any of the returned results with well-substantiated methods and proven values.
Give or take some minute floating point inaccuracy.

Most online resources work well for fairly round ellipses and are nothing more than approximations, just like the Simpson's rule.
They all start to suffer from inaccuracies below a ratio of 0.75.
The accuracy can drop below 0.1%, 0.5% or even 1% in edge cases.
The method implemented in QCAD is no exception.

Already mentioned in the above post 49156:

• The real solution is implementing the elliptic integrals.
  Especially those of the second kind are of interest to us.
  E(m) and E(ø|m) where m is the ellipse eccentricity squared.

Then it is good to know that the eccentricity is equal to sqrt(1-b²/a²) or sqrt(1-ratio²).
Because the ratio is a REllipse Public Attribute itself, it must not be calculated.
Both the parameters a and ratio are stored and fixed values for a certain ellipse.
The modulus m is thus equal to 1-ratio².

Another observation is that elliptic integrals their reference base is the first minor point.
Once you have wrapped your mind around that, it all starts to explain itself.

E(m) is the arc length of 1/4 unit ellipse going from a minor to a major point or vice versa.
Then is a*E(m) the arc length of 1/4 ellipse with semi-major radius a.
And 4a*E(m) is the circumference of an ellipse with semi-major radius a.

An arbitrary arc length is a bit more complicated because it is referenced from t=Pi/2.
a*E(m) + a*E(ø-Pi/2|m) is the arc length from t=zero to t=ø with semi-major radius a.
Then is a*E(ø2-Pi/2|m) - a*E(ø1-Pi/2|m) the arc length from t=ø1 to t=ø2.
In this last equation we eliminate a*E(m) - a*E(m) = zero in the subtraction of two zero based arcs.


It sounds not very fair to disguise an unsolvable function as f(x) = E(..).
Simply said, there is no exact solution and that is why they have their own name.
But there are well established methods to approximate E(..) nearly exact in any number system.
For a double in floating point we must settle for 15-16 significant digits with some reserves for the last, the last two.


In QCAD we are also interested in the inverse, the parameter angle ø for a given length.
Or for what ø the function a*E(ø|m) - length becomes zero ... The root of that function.
This is approximated with Newton's method and converges quickly to a final result.
An initial guess its length is corrected by E() divided by its first derivative and repeated until the length fits.
The final argument ø then relates directly to the given length.
Some pre and post conditioning must be preformed to work in the first quadrant.


With all this implemented (JS based) and in working order I can now start the process of verifying things.
- The summation of ellipse arcs, subtracting, when the sweep is 2Pi it must be 4a*E(m) long.
I have several routes for verifying 4a*E(m).
- Going forth and back between parameter angles and length.
- ....
Preliminary results look fine and I don't expect major flaws.

Regards,
CVH