## Octahedral Angles

Thursday, September 10th, 2009 at 11:34 am

This is a follow-on from the Diamond Angle article about useful encoding of plane angles without trigonometry.

The same trick can be applied to 3D vectors using a regular octahedron instead of a diamond.

The corners of the regular octahedron are:

```  { (1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1) }
```

Another way of writing it that shows off how rubbish mathematical notation can be sometimes is:

```  { (x, y, z) | x, y, z in {-1, 0, 1} where 2 of the values are 0 }
```

The corners of the cube, on the other hand, can be expressed as:

```  { (x, y, z) | x, y, z in {-1, 1} }
```

Let’s start with the conventional angular encoding of a vector into latitude and longitude:

```def RadianPolarAngle(x, y, z)
r = sqrt(x*x + y*y)
return (atan2(x, y), atan2(r, z))
```

This is obviously decoded by:

```    (x, y, z) = (sin(a)*sin(b), cos(a)*sin(b), cos(b))
```

The ranges for the two encoding angles are 0 < a < 360 and 0 < b < 180, but the disadvantages include the use of slow trigonometry, and the existence of two poles where the encoding density is high.

An uneven density means an unnecessary loss of precision:

Suppose I encode this pair of values into two bytes (16 bits) by rescaling them like so:

```   int(a/360*256) * 256 + int(b/180) *256
```

Then the directions near b=0 and b=180 will be preserved far more accurately than the directions near b=90. This represents an inefficient waste of my small and finite number of encoding values.

Here is a picture of the point distribution stepping every 5 degrees: The three blue lines are the x, y, and z (pointing up) axes.

The slow trig functions problem can easily be solved by substituting DiamondAngle in place of atan2 so that 0 < a < 4 and 0 < b < 2. This gives a picture like so:  The right hand picture plots the vectors without normalization so you can see how they lie on the surface of an octahedron. But looking at the left hand picture, you can see a very slight higher density in line with the x and y axes (due to the DiamongArg encoding), as well as the very high density near the points (0,0,1) and (0,0,-1) axis.

A better solution is to extend the concept of the DiamondAngle to the OctahedronAngle by breaking the 3D direction into octants (rather than quadrants) and encoding the intersection with the triangle spanning the relevant octant.

We can number the octants between 0 and 7 like so:

```    octantnumber(x, y, z) = (x > 0 ? 0 : 1) + (y > 0 ? 0 : 2) + (z < 0 ? 0 : 4)
```

and then, given x, y, z > 0, we can encode the vector’s intersection with the triangle spanning the octant:

```    triangle[ (1,0,0), (0,1,0), (0, 0, 1) ]
```

like so:

```    (a, b) = ( x/(x+y+z), y(x+y+z) )
```

where a, b, > 0 and a + b < 1

The inverse function (subject to normalization) is:

```   (x, y, z) = (a, b, 1 - a - b)
```

So, the first version of our encoding is obvious: We can encode the number of the octant in 3 bits, and the direction within the chosen octant with two positive values whose sum is less than 1.

But we can make this more compact.

There are 8 such right-angled triangles in the square defined by -1 < a, b < 1, so we could use the number of the triangle to number the octant:

```    trianglenumber(a, b) = (x > 0 ? 0 : 1) + (y > 0 ? 2 : -2) + (abs(a) + abs(b) < 1 ? 0 : 4)
```

and substituting a valid range for a and b like so:

```    (a, b) = (abs(a) + abs(b) < 1 ? (abs(a), abs(b)) : (1 – abs(a), 1 – abs(b)) )
```

What we are working towards here is a wrapping of the 2×2 square onto the octahedron, and trying to have the same effect we wrapped the line segment [0, 4] around the diamond. But can we make it continuous?

We can get continuity on the top pyramid of four triangles easily. For z > 0

```   OctahedronArg(x, y, z) = (x / (abs(x) + abs(y) + z), y / (abs(x) + abs(y) + z))
```

And the inverse (subject to normalization) is:

```   (x, y, z) = (a, b, 1 – abs(a) – abs(b))
```

Here is the picture of it. You can see the DiamondArg style slight elevation in density in the planes of the x, y, z axes, but it is nothing like the huge density near the point (0,0,1) axis when we do it the traditional way.

Maybe a slightly better visualization involves imposing a checkerboard pattern by discarding the points where:

```   int((a + 2) * 10 + 0.4) + int((b + 2) * 10 + 0.3) mod 2 == 1
```

to make it out of phase with the axes and plot how it crosses them Extending this down to the lower pyramid is less straight forward. The outer triangles of the quadrants (where abs(a) + abs(b) > 1) have to equate their open edges, preferably without discontinuities.

The x, y > 0 double-octant spans the upper and lower pyramids like so (see middle diagram above):

```   (x, y, z) =  (x + y < 1 ? (a / (x + y + z), b / (x + y + z)) : ((1 – y) / (x + y – z), ((1 – x) / (x + y – z)))
```

The x and y coordinates seem to reverse, because at the x + y = 1 line (x, y) = (1 – y, 1 – x).

The inverse function, for a, b > 0 is:

```   (x, y, z) = (a + b < 1 ? (a, b, 1 – a – b) : (1 – b, 1 – a, 1 – a – b) )
```

It’s possible to continue for the final 3 octants and and have something that is as good as it’s going to get:

```def InverseOctahedronAngle(a, b)
h = abs(a) + abs(b)
if h < 1
return (a, b, 1 – h)

if a > 0 and b > 0
(x, y) = (1 – b, 1 – a)
if a > 0 and b < 0
(x, y) = (1 + b, -1 + a)
if a < 0 and b > 0
(x, y) = (-1 + b, 1 + a)
if a < 0 and b < 0
(x, y) = (-1 – b, -1 – a)

return (x, y, 1 – h)
```

For the sake of completion, the encoding function is:

```def OctahedronAngle(x, y, z)
s = abs(x) + abs(y) + abs(z)
if z > 0
return (x / s, y / s)

if x > 0 and y > 0
return (1 – y / s, 1 – x / s)
if x > 0 and y < 0
return (1 + y / s, -1 + x / s)
if x < 0 and y > 0
return (-1 + y / s, 1 + x / s)
if x < 0 and y < 0
return (-1 - y / s, -1 - x / s)
```

I’ve been using greater than signs (>) rather than greater than or equal signs, because I don’t know the code for them, and because these are the boundaries where it gets interesting.

Let’s see how everything stacks up with regards to continuity across the z = 0, abs(a) + abs(b) = 1 boundary in the InverseOctahedronAngle(a, b) function.

```  a > 0, b > 0  ==> a + b = 1 ==>   (1 – b, 1 – a) = (a, b)
a > 0, b < 0  ==> a - b = 1 ==>   (1 + b, -1 + a) = (a, b)
a < 0, b > 0  ==> -a + b = 1 ==>   (-1 + b, 1 + a) = (a, b)
a < 0, b < 0  ==> -a - b = 1 ==>   (-1 - b, -1 - a) = (a, b)
```

That accounts for all the internal boundaries of the -1 < a, b < 1 square. But the 8 half edges of the square equate to the 4 lower edges of the octahedron, and the four corner points equate to the bottom vertex of the octahedron at (0, 0, -1).

Let’s start with those corners where a, b in { -1, 1 }

```  a > 0, b > 0  ==> (a, b) = (1, 1) ==>   (1 – b, 1 – a) = (0, 0)
a > 0, b < 0  ==> (a, b) = (1, -1) ==>   (1 + b, -1 + a) = (0, 0)
a < 0, b > 0  ==> (a, b) = (-1, 1) ==>   (-1 + b, 1 + a) = (0, 0)
a < 0, b < 0  ==> (a, b) = (-1, -1) ==>   (-1 - b, -1 - a) = (0, 0)
```

Now consider b = 1. For t > 0, I would hope for (t, 1) to equate to (-t, 1):

```  a > 0, b > 0  ==> (a, b) = (t, 1) ==>   (1 – b, 1 – a) = (0, 1 - t)
a < 0, b > 0  ==> (a, b) = (-t, 1) ==>   (-1 + b, 1 + a) = (0, 1 - t)
```

So it works!

The same applies to the halves of the other three sides of the square, proving that we have a good mapping and are not going to have to worry too much about precisely knowing which octant a point is within (as implied by all those missing greater-than-or-equal to signs in the above snippets).

It also suggests we could do a pretty good distance function in this space that would satisfy the metric space rule, though this would need to be checked out.

This function was developed for the scallop bisector algorithm in order to pack the directions into one 32 bit word.

A long time ago, when I was working for NC Graphics, I solved the same problem when designing a depth buffer structure by encoding (x, y, z) into three bytes by mapping each -1 < x < +1 into an integer less than 256. The fourth byte was was the colour.

The following year OpenGL became available and we started writing Machining Strategist.