Calculating the average angle

[dummzeuch 1390]

Just in case anybody has got write access to RosettaCode:

 

There is a bug in the Pascal implementation:

https://rosettacode.org/wiki/Averages/Mean_angle#Pascal

 

It checks the average cosine against eps like this:

 

  if c > eps then
    MeanAngle := RadToDeg(arctan2(s,c))
  else
    // Not meaningful
    MeanAngle := 0.0;

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

  if Abs(c) > eps then
    MeanAngle := RadToDeg(arctan2(s,c))
  else
    // Not meaningful
    MeanAngle := 0.0;

Edited September 3, 2019 by dummzeuch

[Boris Novgorodov 10]

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 September 3, 2019 by Boris Novgorodov

[David Heffernan 2298]

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. 

[dummzeuch 1390]

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.

[Uwe Raabe 1882]

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.

[dummzeuch 1390]

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
    sincos(DegToRad(a[i]),s,c);
    Sumsin := sumSin+s;
    SumCos := sumCos+c;
  end;
  Result := RadToDeg(arctan2(SumSin,SumCos));
end;

* no division by cnt

* no check for eps

Edited September 4, 2019 by dummzeuch

[dummzeuch 1390]

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:

Averaging Angles - math doctors

Edited September 4, 2019 by dummzeuch

[David Heffernan 2298]

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 September 4, 2019 by David Heffernan