Recommended Posts

Just in case anybody has got write access to RosettaCode:

There is a bug in the Pascal implementation:

It checks the average cosine against eps like this:

if c > eps then
else
// Not meaningful
MeanAngle := 0.0;

But it should actually check the absolute value of c like this:

if Abs(c) > eps then
else
// Not meaningful
MeanAngle := 0.0;
Edited by dummzeuch

You are right, but also note that condition"if Abs(c) <= eps" really does not refer to "Not meaningful"

Only simultaneous zero sums like "if Abs(c)  + Abs(s) <= eps﻿" should be considered as extra case

Edited by Boris Novgorodov

Actually the comparison against eps is wrong and should be removed. Call arctan2 unconditionally. Also the divide by cnt is utterly wasteful.

Bogus comparisons against small numbers drive me mad.

• 1

9 hours ago, David Heffernan said:

Actually the comparison against eps is wrong and should be removed. Call arctan2 unconditionally. Also the divide by cnt is utterly wasteful.

Bogus comparisons against small numbers drive me mad.

So, what do you propose the code should look like instead? (Especially when removing the divide by cnt.) I must admit that I have only compared what it does to some other contributions on the site, not really tried to research the algorithm itself.

55 minutes ago, dummzeuch said:

(Especially when removing the divide by cnt.)

ArcTan2 produces the same result when both of its parameters are multiplied with the same factor. Thus dividing both by cnt before calling ArcTan2  has no influence on the result.

In addition, ArcTan2 is pretty well capable of handling some special values like c = 0 (at least in non-ancient Delphi versions). There is no need to catch this case in advance.

So it should look like this:

function MeanAngle(const a:tAngles;cnt:longInt):double;
// calculates mean angle.
// returns 0.0 if direction is not sure.
var
i : LongInt;
s,c,
Sumsin,SumCos : extended;
begin
IF cnt = 0 then
Begin
Result := 0.0;
EXIT;
end;

SumSin:= 0;
SumCos:= 0;
For i := Cnt-1 downto 0 do
Begin
Sumsin := sumSin+s;
SumCos := sumCos+c;
end;
end;

* no division by cnt

* no check for eps

Edited by dummzeuch

Having read a bit about the theory, I think the check for sumcos (or sumcos/count) being significantly different from 0 is valid. The algorithm converts the angles to vectors and adds these vectors, it then calculates the angle of the resulting vector. To do that, it calculates the arctan of sumsin / sumcos which would be undefined for sumcos = 0 (or very near 0).

But since it uses vector addition, also sumsin = 0 and sumcos = 0 would be a problem because that resulting vector would point to the origin and therefore it would not be possible to calculate an angle. I guess that's what Boris meant here:

21 hours ago, Boris Novgorodov said:

Only simultaneous zero sums like "if Abs(c)  + Abs(s) <= eps﻿" should be considered as extra case

Also, this does not work as expected in some cases:

MeanAngle(0, 0, 30) does not return 10, as I would have expected, but 9.896.

There is a nice explantion of the math here:

Edited by dummzeuch

No, it's using arctan2 rather than arctan.

Just read every single other implementation on that page.

As far as both sumsin and sumcos being very close to 0, I would still advise not using a tolerance.  All you can achieve by that is giving a less inaccurate answer than is possible to reach.

If you want to do better, use one of the algorithms used for dealing with roundoff when calculating sums, for instance Kahan summation.

Edited by David Heffernan

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account