Быстрые равномерно распределенные случайные точки на поверхности единичного полушария

Ответ 1

Самый простой способ создания равномерного распределения на единичной сфере (независимо от его размерности) состоит в том, чтобы нарисовать независимые нормальные распределения и нормализовать результирующий вектор.

Действительно, например, в размерности 3 e ^ (- x ^ 2/2) e ^ (- y ^ 2/2) e ^ (- z ^ 2/2) = e ^ (- (x ^ 2 + y ^ 2 + z ^ 2)/2), поэтому совместное распределение инвариантно вращением.

Это быстро, если вы используете быстрый генератор распределения (Ziggurat или Ratio-Of-Uniforms) и быструю нормализацию (google для "быстрого обратного квадратного корня" ). Не требуется трансцендентного вызова функции.

Кроме того, Марсалья не является однородным на полусфере. У экватора будет больше точек, так как точка соответствия на 2D-диске, точка на полусфере не является изометрической. Последнее кажется правильным, хотя (однако я не делал расчет, чтобы убедиться в этом).

Ответ 2

Если вы берете горизонтальный срез единичной сферы с высотой h, ее площадь поверхности равна 2 pi h. (Так Архимед вычислил площадь поверхности сферы.) Таким образом, z-координата равномерно распределена в [0,1]:

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
osRay.c._z = dsfmt_genrand_close_open(&dsfmtt);

xyproj = sqrt(1 - osRay.c._z*osRay.c._z);
osRay.c._x = xyproj*cos(azimuthal); 
osRay.c._y = xyproj*sin(azimuthal);

Кроме того, вы можете сэкономить некоторое время, вычислив cos(azimuthal) и sin(azimuthal) вместе - см. этот вопрос о стеке для обсуждения для обсуждения.

Отредактировано для добавления: ОК, теперь я вижу, что это всего лишь небольшая настройка вашего третьего метода. Но он вырезает шаг.

Ответ 3

Вы пытались избавиться от asin?

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
sin2_zenith = dsfmt_genrand_close_open(&dsfmtt);
sin_zenith = sqrt(sin2_zenith);

// Calculate the cartesian point
osRay.c._x = sin_zenith*cos(azimuthal); 
osRay.c._y = sin_zenith*sin(azimuthal);
osRay.c._z = sqrt(1 - sin2_zenith);

Ответ 4

Это должно быть быстрым, если у вас быстрый RNG:

// RNG::draw() returns a uniformly distributed number between -1 and 1.

void drawSphereSurface(RNG& rng, double& x1, double& x2, double& x3)
{
    while (true) {
        x1 = rng.draw();
        x2 = rng.draw();
        x3 = rng.draw();
        const double radius = sqrt(x1*x1 + x2*x2 + x3*x3);
        if (radius > 0 && radius < 1) {
            x1 /= radius;
            x2 /= radius;
            x3 /= radius;
            return;
        }
    }   
}

Чтобы ускорить его, вы можете переместить вызов sqrt внутри блока if.

Ответ 5

Я думаю, что проблема, которую вы испытываете с неравномерными результатами, состоит в том, что в полярных координатах случайная точка на окружности неравномерно распределена на радиальной оси. Если вы посмотрите на область на [theta, theta+dtheta]x[r,r+dr], для фиксированных theta и dtheta, область будет отличаться от разных значений r. Intuitivly, есть "больше области" дальше от центра. Таким образом, для этого необходимо масштабировать свой случайный радиус. У меня нет доказательства, лежащего вокруг, но масштабирование r=R*sqrt(rand), при этом r является радиусом круга и rand начинает случайное число.

Ответ 6

Второй и третий методы фактически создают равномерно распределенные случайные точки на поверхности сферы со вторым методом (Marsaglia 1972) обеспечивая наивысшие времена работы примерно в два раза быстрее на четырехъядерном процессоре Intel Xeon с тактовой частотой 2,8 ГГц.

Как отмечено Alexandre C, существует дополнительный метод, использующий нормальное распределение, которое расширяется до n-сфер лучше, чем методы, которые я представил.

Эта ссылка даст вам дополнительную информацию о выборе равномерно распределенных случайных точек на поверхности сферы.

Мой первоначальный метод, обозначенный TonyK, не создает равномерно распределенных точек и, скорее, смещает полюсы при создании случайных точек. Это необходимо из-за проблемы, которую я пытаюсь решить, но я просто предположил, что она будет генерировать равномерно случайные точки. Как предложено Pablo, этот метод можно оптимизировать, удалив вызов asin(), чтобы сократить время выполнения примерно на 20%.

Ответ 7

Первая попытка (неправильная)

point=[rand(-1,1),rand(-1,1),rand(-1,1)];
len = length_of_vector(point);

Редакция:

А что?

while(1)
 point=[rand(-1,1),rand(-1,1),rand(-1,1)];
 len = length_of_vector(point);
 if( len > 1 )
     continue;
 point = point / len
     break

Acception здесь приблизительно 0,4. Это означает, что вы откажетесь от 60% решений.