Optimize this compass directions code

General Help
1

I thought this might be fun for someone to improve on some code I just wrote.

I’m converting bearings to compass directions. It needs to be in the most efficient way possible as this code will be called millions of times as part of an already time consuming process.

The following code appears to be working and I’m pretty happy with it but perhaps someone can suggest an improvement, either in terms of performance or readability.

The case statement is designed to hande most of the calls, each of its clauses represents the integer values that the rounded bearing can safely be assumed to be within.

The else clause of the case statement does an actual accurate check against the floating point boundary values for the cases when the bearing did not fall within the range of safe integer values.

Suggestions?

TCompassDirection = (cdNorth, cdNNE, cdNE, cdENE, cdEast, cdESE, cdSE, cdSSE, cdSouth, cdSSW, cdSW, cdWSW, cdWest, cdWNW, cdNW, cdNNW);

function CompassDirectionOf(const inBearing : double) : TCompassDirection;
const
  DEGREES_PER_DIRECTION = 360 / (Ord(High(TCompassDirection)) + 1);
  ANTI_CLOCKWISE_OFFSET = DEGREES_PER_DIRECTION / 2;
  BOUNDARY_MARGIN = 1;
var
  cd : TCompassDirection;
begin

  case Round(inBearing) of

    0
    ..Round(((Ord(cdNorth) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdNorth;

    Round((Ord(cdNNE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdNNE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdNNE;

    Round((Ord(cdNE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdNE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdNE;

    Round((Ord(cdENE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdENE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdENE;

    Round((Ord(cdEast) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdEast) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdEast;

    Round((Ord(cdESE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdESE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdESE;

    Round((Ord(cdSE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdSE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdSE;

    Round((Ord(cdSSE) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdSSE) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdSSE;

    Round((Ord(cdSouth) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdSouth) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdSouth;

    Round((Ord(cdSSW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdSSW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdSSW;

    Round((Ord(cdSW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdSW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdSW;

    Round((Ord(cdWSW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdWSW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdWSW;

    Round((Ord(cdWest) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdWest) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdWest;

    Round((Ord(cdWNW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdWNW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdWNW;

    Round((Ord(cdNW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdNW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdNW;

    Round((Ord(cdNNW) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN
    ..Round(((Ord(cdNNW) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) - BOUNDARY_MARGIN:
      Result := cdNNW;

    Round((Ord(cdNorth) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) + BOUNDARY_MARGIN + 360 { negative degree values converted to almost 360 degree values }
    ..360:
      Result := cdNorth;

  else
    { bearing is between the integer boundary margins, so do
    an exact check against the floating point boundaries }
    Result := cdNorth; { edited from original post - Thanks Mark Griffiths }
    for cd := cdNorth to cdNNW do begin
      if inBearing <= (((Ord(cd) + 1) * DEGREES_PER_DIRECTION) - ANTI_CLOCKWISE_OFFSET) then begin
        Result := cd;
        Exit;
      end;
    end;
  end;
end;
1 Like
2

A pre calculated lookup table should be a bit quicker.

If you multiply the input value by 2 before rounding, I think that you could also avoid the need for the slower check.

Just checking, what value do you get from this code for a bearing value that is just North of the NNE/N boundary?

3

North - the ANTI_CLOCKWISE_OFFSET drags everything back 11.25 degrees so the actual true direction is in the midpoint of its range.

The lookup table idea is worth trialling, thanks for that suggestion.

Doubling the bearing to remove the slower check in the else clause, interesting, not sure it would work but I’ll investigate.

4

Are you sure? My reading of the code is that it’s only checking between negative values and 11.25 degrees and that values between 348.75 degrees and 360 degrees aren’t being checked, leaving the result undefined.

Is the compiler warning you that the result could be undefined?

I could be misreading your code of course.

I’m assuming that you’ll never be passing in values < 0 or > 360?

5

How does inBearing change? Can it be any direction or is inBearing (n+1) related to inBearing (n) somehow? If it can only change +/- a certain number of degrees could you track the last bearing and limit your search to a subset?

6

0 <= inBearing <= 360

7

You’ve missed the last clause of the case statement just before the else.

The checking for North is split between the first clause and that last clause. If the clauses for a case statement are in order, you get more efficient code generation. Testing for those high values at the top of the statement would mean less efficient code gen for the whole case statement.

8

I thought I’d better double check that assumption. It’s something I heard from Danny Thorpe (RIP) a long time ago so I’ve always taken it as gospel. I just moved that last clause from the bottom to the top so they’re now out of order and the assembly that was generated didn’t change in any meaningful way.

I’m going to assume that there were more optimizations done to the case statement after Danny left, rather than he was mistaken.

9

I’m not meaning the main part of the case statement - I’m talking about the part in the else statement.

You didn’t answer my question as to whether the compiler warns about the return value possibly being undefined, so I just tried it and confirmed that it does give a warning.

One additional consideration with my suggestion to double the value and use a lookup table is that by my understanding, Delphi uses bankers rounding which means that when inBearing is right on the border, some will fall to the left and some to the right. Whether that’s an issue or not for you is up to you to decide. There would be a fix for that if it is important.

10

Nice catch, will fix.

Bankers rounding and other ambiguities are why I went with the else clause slower check for those bearings closer than a degree to the boundary. None of the alternative rounding modes looked like they were going to really improve that.

Those close to boundary bearings and their slower check may not be a big issue, I’ll do some performance testing and see whether the slow check requires urgent optimization or if I just come back and try the doubling method for fun another time.

Thanks for your suggestions and I hope it was an interesting diversion for you.

11

If you want to do classical rounding, you can use:

Trunc(SimpleRoundTo(Value, 0))

It’s not going to be as quick as regular rounding though, so if speed is critical, I’d suggest using a lookup table with regular rounding and for those values where the rounding method matters, have a special value which you can then use to flag the requirement to do a normal rounding and do a lookup in a 2nd table.

I’d generate the lookup tables at runtime rather than having large constant arrays that would be hard to verify.

This would be a good candidate for unit testing I think - use the slow code to generate the correct answers and verify all results between 0 and 360 in 1/4 or even 1/8 degree increments.

No worries. It was an interesting diversion, thanks :slight_smile:

12

SimpleRoundTo() won’t be allowed in the case statement clause criteria though so I’d be using different rounding in those from the input bearing. I could always replace those with literals though so it might be worth a go.

The lookup table implementation will be easy to try so I’ll let you know how it compares.

13

Your case statement has lots of calculations in it. You can minimize calculations like this:

function CompassDirectionOf(const inBearing: double): TCompassDirection;
const
  DEGREES_PER_DIRECTION = 360 div (Ord(High(TCompassDirection)) + 1);
  ANTI_CLOCKWISE_OFFSET = DEGREES_PER_DIRECTION div 2;
begin
  Result := TCompassDirection((Round(inBearing) + ANTI_CLOCKWISE_OFFSET) div DEGREES_PER_DIRECTION);
end;
14

I like your thinking. It doesn’t look like it would handle the edge cases properly but that could be fixed.

15

Nice one @Scott_Sedgwick.

The calculations in my switch statement don’t really matter as they’re all evaluated at compile time but this is a much more succinct and readable solution.

I’m a little worried about the integer division but with the limited test dataset I have at the moment it’s returning the same results as my implementation.

I’ll try it out with a broader set of test data, compare the two methods and report back.

16

Try this:

Result := TCompassDirection(
trunc((bearing+11.25)/22.5) mod 16
)
Regards
Peter Lovelace

17

I had a fairly cuckoo idea, and had a bit of fun following it up. And learnt a bit about doubles too. :slight_smile:

I had a look to see what the binary format of the compass point angles were …
( these are the bytes in reverse order 7 → 0 )

image
image1479×1457 183 KB

They are “exact binary” numbers, so I wanted to try comparing bytes 6 & 5 from the double against these values.
I made 10,000 random doubles from 0 to 365.0, and running through them 1000 times gave a run time around the 200-400 ms mark.

// - - -

I ended up measuring CompassDirectionOf at around 245 ms.
This was about equal to straight up comparing inBearing to an array of doubles ( 22.5, 45.0, … )

CompassDirectionOf2 from @Scott_Sedgwick was simple and fast, around 180 ms.
It’s pretty fascinating how efficient the straight calculation is, and how many things can cost performance.

I eventually got mine down to 140-150 ms, but the code would probably be judged kinda ugly.
An interesting exercise, and it makes me appreciate the compiler.

// - - -

Getting access to the two bytes of the double and swapping their order was a speed bump.
Double record helper double.Bytes was ok.
Array[0..7] of Byte absolute inBearing; was faster.

I was a bit surprised I had to unroll

  for var cd := cdNorth to cdNNW do                                           // 175 ms
      if I < CompPoints[cd] then begin
         Result := cd;
         Exit;
      end;

to

    if I < cp[ 0] then begin Result := TCompassDirection( 0); exit; end       // 150 ms
    else
    if I < cp[ 1] then begin Result := TCompassDirection( 1); exit; end
    else 
    if I < cp[ 2] then begin Result := TCompassDirection( 2); exit; end
    else ...

to gain about 20-30 ms.

My various attempts are here: GitHub - pmcgee69/double-precision-to-enum-perf-test

Inlining sped things up, but I didn’t record any numbers using it.
Starting with Result := cdNorth; sometimes helped and sometimes killed performance.

2 Likes
18

I’m a complete Delphi newbie; I joined this group from a post by Paul McGee on FB.
I’ve always used C++ Builder, rather than Delphi.
Anyway, I’m reading your code and trying to work out what you want as an output? I guess this routine is just a small part of a larger program. Are you just looking for the nearest of the 16 compass points? I’m curious what the application of that would be and in that case, why the need for such accuracy?

1 Like
19

Welcome Colin. C++ Developers are welcome here, though you’re always more likely to get a response if you keep the C++ as simple as possible (you probably know more Delphi than most of us here know C++).

You’re right about the output and purpose of the code. It’s a function that returns one of the values of the enumerated type TCompassDirection declared at the top of the code.

I can’t say too much about the application other than it’s analysing GPS data from trackers worn by sports people.

As for the need for accuracy, well… why not be as accurate as possible?

Eventually someone would spot that the boundary for the North range was inconsistent with the range for NE for example. That may not be a big issue in itself but it will raise questions in their mind about what other calculations may have inconsistencies.

20

Thanks Lachlan.
The code seems really complex to me just to generate a direction, although it is probably nicely written Delphi code.
My approach would have been to not bother with rounding at all and then simply compare the double value with a range either in a switch or an if…else.
For example, a bearing of 027° would fall in the range of >=011.25° and <033.75° and would give NNE.
If the bearing was (say) 033.76° it would fail that test and be caught by the next step >=033.75° and < 056.25 and give NE. That way you don’t have to worry about either side of the N and S values and the changes of E/W hemispheres.
Of course, if your code is used in other parts of the program then this is not relevant, I was just looking at a simple direction.
(I spent 50 years at sea and aviation and GIS systems are my hobby so compasses and bearings have been a big part of my life, hence the interest in your question :slightly_smiling_face:)