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

Question

Я пытаюсь создать однородные случайные точки на поверхности единичной сферы для программы трассировки лучей Монте-Карло.Когда я говорю равномерно, я имею в виду, что точки равномерно распределены по площади поверхности.Моя текущая методология состоит в том, чтобы вычислять однородные случайные точки на полушарии, указывающие на положительную ось z и базу в плоскости xy.

Случайная точка на полусфере представляет направление излучения теплового излучения для диффузного серого излучателя.

Я получаю правильный результат, когда использую следующие вычисления:

Примечание: dsfmt * is вернет случайное число от 0 до 1.

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
zenith = asin(sqrt(dsfmt_genrand_close_open(&dsfmtt)));

// Calculate the cartesian point
osRay.c._x = sin(zenith)*cos(azimuthal); 
osRay.c._y = sin(zenith)*sin(azimuthal);
osRay.c._z = cos(zenith);

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

Метод отклонения Marsaglia 1972 *

do {
   x1 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   x2 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   S = x1*x1 + x2*x2;
} while(S > 1.0f);


osRay.c._x = 2.0*x1*sqrt(1.0-S);
osRay.c._y = 2.0*x2*sqrt(1.0-S);
osRay.c._z = abs(1.0-2.0*S);

Расчет аналитических декартовых координат

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
u = 2*dsfmt_genrand_close_open(&dsfmtt) -1;
w = sqrt(1-u*u);

osRay.c._x = w*cos(azimuthal);
osRay.c._y = w*sin(azimuthal);
osRay.c._z = abs(u);

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

Кроме того, последние два метода дают одинаковые конечные результаты, однако я уверен, что они неверны, поскольку я сравниваю их с аналитическим решением.Однако распределение не дает правильного результата.

Ошибка в моей реализации или я пропустил основную идею во втором и третьем методах?

Answer 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 диске <-> точка на полусферене изометрический. Последний кажется правильным (однако я не сделал вычисления, чтобы гарантировать это).

Answer 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) вместе - см. этот вопрос о стекопереработке для обсуждения.

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

Answer 3

Это должно быть быстро, если у вас быстрый 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.

Answer 4

Вы пытались избавиться от 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);

Answer 5

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

Answer 6

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

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

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

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

Answer 7

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

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

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