Как генерировать случайные точки внутри октаэдра, не отбрасывая?

Question

Мне нужны случайные точки внутри октаэдра, равномерно распределенные. Я определяю октаэдр как объем, где все точки удовлетворяют abs(x) + abs(y) + abs(z) <= 1, где abs дает абсолютное значение. IE: каждая из шести вершин будет на оси, 1 от 0,0,0. Может быть, вы могли бы назвать это единичным октаэдром.

Имея в виду это определение, я могу наивно сгенерировать такую ​​точку:

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))

Проблема в том, что это склоняется к небольшим значениям yи меньшие значения z. Очевидно, что это не равномерное распределение. Также ясно, что все эти точки находятся только в одном из восьми квадрантов.

Я избегаю метода отбрасывания , потому что эта функция будет вызываться много раз, и, похоже, мне следуетможет сделать лучше, чем выбрасывать большинство точек.

Обратите внимание, что двойник октаэдра является кубом. Из-за этого у меня есть предположение, что может существовать простая функция для перевода любой точки внутри куба в октаэдр, но это всего лишь интуиция, которую я до сих пор изучаю.

Answer 1

Вот идея: выборочные точки из распределения Дирихле в D + 1, выберите D точек так, чтобы вы были едины в симплексе

x ₀ + x ₁ + x ₂ <= 1, x <sub>i > = 0

Чтобы создать октаэдр, случайным образом выберите октант, чтобы поставить точку.

Код на 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}")

А вот быстрый график с 2K точками

Answer 2

Вы знаете, как выбирать точки в кубе с равномерным распределением, и куб можно разбить на восемь квадратных пирамид. (Извините, я не могу предоставить графику.)

Я бы начал с куба: abs(x) <= 1; abs(y) <= 1; abs(z) <= 1

Укажите точку в нем (вектор столбца, (x, y, z)), затем отразитечтобы привести его в «верхнюю и нижнюю» пирамиды:

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

Затем переверните две пирамиды, чтобы создать октаэдр:

if z>0
 z = 1-z

if z<0
 z = -1-z

Затем поверните и масштабируйте:

multiply by

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

Answer 3

Редактировать: Как отмечает @Tonci, сумма 3 равномерных случайных величин не является равномерной случайной величиной. Фактически для a,b,c,d вероятность быть около 0 примерно в 1,5 раза больше вероятности около 1 (или около -1). Поэтому ответ ниже является лишь приблизительным. Я оставляю это здесь, потому что, возможно, это может дать кому-то правильное вдохновение.

Октаэдр определяется 8 разграничивающими плоскостямиУравнения этих плоскостей являются обобщением abs(x) + abs(y) + abs(z) <= 1, где abs каждый раз может быть либо положительным, либо отрицательным. Итак:

x + y + z <= 1; x + y - z <= 1; x - y + z <= 1; x - y - z <= 1;
-x + y + z <= 1; -x + y - z <= 1; -x - y + z <= 1; -x - y - z <= 1

Объединение уравнений противоположной плоскости:

-1 <= x + y + z <= 1; -1 <= x + y - z <= 1; -1 <= x - y + z <= 1; -1 <= x - y - z <= 1

Теперь выберите значения для каждого из этих выражений:

x + y + z = a; x - y + z = b; x - y - z = c; x + y - z = d

Или для каждого случайного x,у, z внутри октаэдра у нас есть 4 значения:

-1 <= a <= 1; -1 <= b <= 1; -1 <= c <= 1; -1 <= d <= 1

Работая в обратном порядке, если у нас есть значения для a, b, c, d, мы можем найти соответствующие x, y, z. Эти уравнения линейно зависимы, случайным образом выбирая 3 из a, b, c, d, мы можем вычислить четвертое, а также x, y, z.

Поэтому выбираем случайные числа a, b, c в трех направлениях. Все в диапазоне от -1 до 1:

• x + y + z = a
• x - y + z = b
• x - y - z = c

Решение для x, y, z: x: a/2 + c/2, y: a/2 - b/2, z: b/2 - c/2.

В коде:

val a: Double = 2*nextDouble()-1  // range -1..1
val b: Double = 2*nextDouble()-1  // range -1..1
val c: Double = 2*nextDouble()-1  // range -1..1
val x: Double = (a+c) / 2
val y: Double = (a-b) / 2
val z: Double = (b-c) / 2

или с меньшим количеством операций:

val a: Double = nextDouble()  // range 0..1
val b: Double = nextDouble()  // range 0..1
val c: Double = nextDouble()  // range 0..1
val x: Double = a+c-1
val y: Double = a-b
val z: Double = b-c