How to generate random points within an octahedron without discarding?

Question

I need random points within an octahedron, uniformly distributed. I am defining an octahedron as the volume where all points satisfy abs(x) + abs(y) + abs(z) <= 1 where abs gives absolute value. IE: each of the six vertices would be on an axis, 1 away from 0,0,0. Maybe you could call it a unit octahedron.

With the definition in mind, I can naively generate a point like so:

val x: Double = nextDouble() // 0-1 range
val y = nextDouble(1.0 -x) // 1-x is upper bound, probably <1
val z = nextDouble(1.0 -(x+y))

The problem is that this leans toward small y values, and smaller z values. Clearly not an even distribution. Also clearly, all these points are in just one of eight quadrants.

I'm avoiding the discard method because this function will be called a lot, and it just seems like I should be able to do better than throwing out the majority of points.

Note that the dual of the octahedron is the cube. Because of this, I have an inkling that there might exist a simple function to translate any point within a cube to be within the octahedron, but that's just an intuition I'm still exploring.

Answer 1

You know how to choose points in a cube with uniform distribution, and a cube can be dissected into eight square pyramids. (Sorry I can't provide graphics.)

I'd start with a cube: abs(x) <= 1; abs(y) <= 1; abs(z) <= 1

Pick a point in it (a column vector, (x, y, z)), then reflect to bring it into the "top and bottom" pyramids:

if abs(x) > abs(z), swap x and z. Equivalently, multiply by

0 0 1
0 1 0
1 0 0

if abs(y) > abs(z), swap y and z. Equivalently, multiply by 

1 0 0
0 0 1
0 1 0

Then invert the two pyramids to make an octahedron:

if z>0
 z = 1-z

if z<0
 z = -1-z

Then rotate and scale:

multiply by

1/2 -1/2  0
1/2  1/2  0
  0    0  1

Answer 2

Here is the idea - sample points from Dirichlet distribution in D+1, select D points such that you're uniform in the simplex

x₀+x₁+x₂ <= 1, x_i >= 0

To make octahedron, randomly select octant to put your point.

Code in Python

import math
import random

def Dirichlet():
    """sample 4d Dirichlet"""
    x0 = -math.log(1.0 - random.random()) # exponential
    x1 = -math.log(1.0 - random.random()) # exponential
    x2 = -math.log(1.0 - random.random()) # exponential
    x3 = -math.log(1.0 - random.random()) # exponential
    s = 1.0/(x0+x1+x2+x3) # scaling

    return (x0*s, x1*s, x2*s, x3*s)

def OctahedronSampling():

    x0, x1, x2, _ = Dirichlet()

    octant = random.randint(0, 7)

    if octant == 0:
        return (x0, x1, x2)
    elif octant == 1:
        return (x0, -x1, x2)
    elif octant == 2:
        return (x0, x1, -x2)
    elif octant == 3:
        return (x0, -x1, -x2)
    elif octant == 4:
        return (-x0, x1, x2)
    elif octant == 5:
        return (-x0, -x1, x2)
    elif octant == 6:
        return (-x0, x1, -x2)
    elif octant == 7:
        return (-x0, -x1, -x2)

    return None

for _ in range(0, 2000):
    x0, x1, x2 = OctahedronSampling()

    print(f"{x0}   {x1}   {x2}")

And here is quick graph with 2K points