FS#2703 - RVector::getAngle() can return (exactly) 2Pi where zero is expected.
Andrew,
In an ideal world the angle of a horizontal vector pointing to the right would be zero.
But the numerical system is not ideal, is it?
Breaking it down in the current reference:
double dp = getDotProduct(*this, RVector(1.0, 0.0));
On itself a of an overkill ... It is simply dp = x
Then handled as 3 cases:
- x / m >= 1.0 the intermediate result is 0.0
- x / m < -1.0 the intermediate result is Pi
- else the intermediate result is acos(x / m), a result limited in quadrant I or II
Not really critical: Case 2 is missing an equality sign.
In a last stage this value is projected to the III and IV quadrant when y is negative.
Assuming x is equal to m and the intermediate result is 0.0
Then y must be zero and the final result is zero.
Even if y differs from zero by a tiny bit, sqrt(x²+y²) can not be equal to x.
Theoretically.
In floating point this would occur when x is large enough and y small enough.
So that x’²+y’² = x’² and x’ / m = 1 holds.
Because of a very small negative value of y the result is then 2Pi - 0.0 = 2Pi.
The flaw: X is normalized to a unit circle, Y not!
Test reveals that y/m < 0.0 will not trigger on nearly negative zero.
The solution would be to normalize the result of for example:
- RVector::getAngle()
- RVector::getAngleTo(other)
- RPolyline::getSegmentAt(i) where i points to a bulging segment
+ many other methods based on getAngle()
Each time something of the sort is used.
But applying RMath::getNormalizedAngle(a) each time to filter out a few rare
cases among hundred of thousands based on getAngle() is very time consuming.
If you can already like in the example case of the getSegmentAt(i) method.
Regards,
CVH