<p>Firstly, as @trichoplax correctly pointed out, your randomPoint function calculates a point in a cube, then uses rejection sampling to return <strong><em>all</em></strong> points that are <strong><em>inside</em></strong> a unit sphere. In order to return points on a sphere, you would need to change the greater than to an equals. That said, rejection sampling is very inefficient. </p>

<p>A better way to sample a sphere, is to sample in spherical space, then transform to cartesian. ie:</p>

<pre><code>float phi = PI * randf();
float theta = 2 * PI * randf();    

float x = radius * std::sinf(phi) * std::sinf(theta);
float y = radius * std::cosf(phi);
float z = radius * std::sinf(phi) * std::cosf(theta);
</code></pre>

<p>However, this is not uniform, and will cause samples to clump at the poles. To prevent this, we transform phi:</p>

<pre><code>float phi = std::acosf(2.0f * randf() - 1.0f);
float theta = 2 * PI * randf();    
</code></pre>

<p>We can use some trig identities to make it a bit more efficient to calculate with the computer:
$$\begin{align*}
\sin(\cos^{-1}(x)) \equiv&amp; \ \sqrt{1 - x^2}\\
\cos(\cos^{-1}(x)) \equiv&amp; \ x\\
\\
\phi =&amp; \ \cos^{-1}(2\xi - 1)\\
\sin(\phi) =&amp; \ sin(\cos^{-1}(2 \xi - 1))\\
=&amp; \ \sqrt{1 - (2 \xi - 1)^2}\\
\cos(\phi) =&amp; \ 2 \xi - 1
\end{align*}$$</p>

<p>So the full change is:</p>

<pre><code>float cosPhi = 2.0f * randf() - 1.0f;
float sinPhi = std::sqrt(1.0f - cosPhi * cosPhi);
float theta = 2 * PI * randf();    

float x = radius * sinPhi * std::sinf(theta);
float y = radius * cosPhi;
float z = radius * sinPhi * std::cosf(theta);
</code></pre>

<p><br/></p>

<p>However, in the case of next event estimation, uniform sampling the whole sphere is inefficient, because a ray can only 'see' half of the sphere at a time. So if we generate a point on the 'back' of the sphere, the sphere will occlude the point, and your calculation will be wasted.</p>

<p>Instead, you generate samples in a cone, which covers the great circle of the sphere, as viewed from the starting point of the ray.</p>

<pre><code>// Sample sphere uniformly inside subtended cone
float rand1 = randf();
float rand2 = randf();

// Compute theta and phi values for sample in cone
float distanceSquared = DistanceSquared(rayOrigin, sphereCenter);

// We use geometry to calculate the angle of the cone (aka, the maximum phi can be when we sample)
// It's easier / cheaper to use geometry to calculate sin/cos phi directly, rather than generating phi and using sin/cos
float sinPhiMaxSquared = radius * radius / distanceSquared;
float cosPhiMax = std::sqrt(1.0f - sinPhiMaxSquared);
float cosPhi = (1.0f - rand1) + rand1 * cosPhiMax;
float sinPhi = std::sqrt(1.0f - cosPhi * cosPhi);

// Phi can be anything in 2 PI
float theta = 2 * PI * rand2;

float x = radius * sinPhi * std::sinf(theta);
float y = radius * cosPhi;
float z = radius * sinPhi * std::cosf(theta);
</code></pre>

<p><strong>NOTE: x, y, z are in LOCAL coordinate space.</strong></p>

<p>So we need to transform to world coordinates:</p>

<pre><code>float3 zAxis = normalize(sphereCenter - rayOrigin);
float3 xAxis;
if (std::abs(zAxis.x) &gt; std::abs(zAxis.y))
    xAxis = float3(-zAxis.z, 0.0f, zAxis.x) / std::sqrt(zAxis.x * zAxis.x + zAxis.z * zAxis.z);
else
    xAxis = float3(0.0f, zAxis.z, -zAxis.y) / std::sqrt(zAxis.y * zAxis.y + zAxis.z * zAxis.z);
yAxis = cross(zAxis, xAxis);


float3 samplePoint = x * xAxis + y * yAxis + z * zAxis;
</code></pre>

<p>And the pdf is calculated as follows:</p>

<pre><code>float pdf = 1 / (2 * PI * (1 - cosPhiMax));
</code></pre>

<p>The code here is heavily influenced / copied from <a href="https://github.com/mmp/pbrt-v3/blob/master/src/shapes/sphere.cpp" rel="noreferrer">PBRT v3</a>. They have a series of classes and functions for sampling from shapes.</p>

<p><br/></p>

<p>Finally, the pdf. In monte-carlo integration, we need to combine the pdf's of each integration we do. In path tracing, we can integrate over many many things. For example, the general rendering equation integrates the incoming light over the hemisphere, depth of field can be treated as an integration over a focal distance, etc.</p>

<p>For next event estimation, you explicitly split the rendering equation into two integrands, direct lighting, and indirect lighting.</p>

<p>Standard rendering equation:
$$ L_{\text{o}}(p, \omega_{\text{o}}) =  L_{e}(p, \omega_{\text{o}}) \ + \  \int_{\Omega} f(p, \omega_{\text{o}}, \omega_{\text{i}}) L_{\text{i}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}} $$</p>

<p>Next Event Estimation:
$$ L_{\text{o}}(p, \omega_{\text{o}}) =  L_{e}(p, \omega_{\text{o}}) \ + \ \int_{\Omega} f(p, \omega_{\text{o}}, \omega_{\text{i}}) L_{\text{i, direct}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}}\ \ + \ \int_{\Omega} f(p, \omega_{\text{o}}, \omega_{\text{i}}) L_{\text{i, indirect}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | d\omega_{\text{i}} $$</p>

<p>In simple naive forward path tracing, everything is is treated as indirect light. In next event estimation, we directly calculate the direct lighting and add it to the indirect lighting. And if we hit a light, we ignore the contribution, since we're calculating the direct lighting.</p>

<p>Since we have two integrations, each will have its own pdf. Aka:
$$L_{\text{o}}(p, \omega_{\text{o}}) =  L_{e}(p, \omega_{\text{o}}) \ \ + \ \sum_{k=0}^{\infty } \frac{f(p, \omega_{\text{o}}, \omega_{\text{i}}) L_{\text{i, direct}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | }{pdf_{direct}}\ \ + \ \sum_{k=0}^{\infty } \frac{f(p, \omega_{\text{o}}, \omega_{\text{i}}) L_{\text{i, indirect}}(p, \omega_{\text{i}}) \left | \cos \theta_{\text{i}} \right | }{pdf_{indirect}}$$</p>

<p>If you want to see how this is implemented, you can check out my implementation <a href="https://github.com/RichieSams/lantern/blob/master/source/renderer/renderer.cpp#L92" rel="noreferrer">here</a>. Note: my light sampling is a bit more complicated, since it does multiple importance sampling. But it can be as simple as:</p>

<pre><code>float3 Renderer::EstimateDirect(Light *light, UniformSampler *sampler, float3a &amp;surfacePos, float3a &amp;surfaceNormal, float3a &amp;wo, Material *material) const {
    float pdf;
    float3 wi;

    // Sample a point on the light and get the pdf
    float3 Li = light-&gt;SampleLi(sampler, m_scene, surfacePos, surfaceNormal, &amp;wi, &amp;pdf);

    // Calculate the brdf value
    float f = material-&gt;Eval(wi, wo, surfaceNormal);
    return f * Li * / pdf;
}
</code></pre>

CC BY-SA 3.0
2016-05-12T23:09:45.017Z
I believe you are correct. I will check my PBRT book when I get home and correct the code. IIRC, I need to add a square root.
RichieSams
CC BY-SA 3.0
2016-05-13T02:01:09.303Z
Yes, @trichoplax is right, sampling both angles uniformly will lead to points bunching at the poles. The simplest way to fix it is sample theta and Y uniformly, then set the XZ radius = $\sqrt{r^2 - y^2}$. Which, BTW, is exactly what you're doing in the second code snippet when you sample `cosPhi` (which is proportional to Y) uniformly. See also [this MathWorld article](http://mathworld.wolfram.com/SpherePointPicking.html).
Nathan Reed
CC BY-SA 3.0
2016-05-17T09:18:09.863Z
Also in MathWorld article: Marsaglia (1972) - point in disk method should trump trig based.
MB Reynolds
CC BY-SA 3.0
2016-05-17T11:56:21.650Z
I would need to do some profiling, but I'm not entirely convinced that rejection sampling would beat trig.
RichieSams
CC BY-SA 3.0
2016-05-18T14:38:39.343Z
If I may ask you, if you feel like of course, I have another question that might definitely solve my project.. Thanks in any case for your precious help so far:
http://computergraphics.stackexchange.com/questions/2461/resulting-probabilty-density-in-path-tracer-for-paths-using-next-event-estimatio
Tarta
CC BY-SA 3.0
2016-05-12T21:31:21.500Z
2016-05-13T02:36:03.227Z
2016-05-13T02:36:03.227Z
RichieSams
<p>I am trying to implement my own path tracer but before arriving to the question I want to give you a short overview:<br>
In the implementation of the rendering equation I use some particular technique in order to sample surfaces. For example: when one of my rays hits a diffuse surface, the next ray bouncing from that surface will be calculated using a <em>Cosine-weighted Random Direction</em>. This implies that in my rendering equation I have to take into account the PDF that this specific random distribution implies and specifically divide for this PDF in my equation. So far everything is okey.</p>

<p>Now, my question. I want to implement a particular technique called "next event estimation" which simply samples the lights. In order to do so, I want to pick a <strong>random point</strong> on my light, which is spherical, by using the following code (C++):</p>

<pre><code>Vec3&lt;float&gt; randomPoint() const
{
    float x;
    float y;
    float z;

        // random vector in unit sphere
    do
    {
        x = ((float)rand() / (RAND_MAX)) * 2 - 1;           // random coordinates between -1 and 1
        y = ((float)rand() / (RAND_MAX)) * 2 - 1;
        z = ((float)rand() / (RAND_MAX)) * 2 - 1;
    } while (pow(x, 2) + pow(y, 2) + pow(z, 2) &gt; 1);        // simple rejection sampling

    return centerOfSphere + Vec3&lt;float&gt;(x, y, z) * radius;
}
</code></pre>

<p>As far as I understand (please correct me if wrong), this implements a <strong>uniformly random sampling</strong>. Also, as far as I have read, this kind of sampling has a PDF = 1 / (b - a) . My question: do I have to use this PDF as well like I do for my <em>Cosine-weighted Random Distribution</em> or not, since it is a uniform distribution? if yes, what's the range (b - a) that PDF talks about? Thanks in advance!</p>

CC BY-SA 3.0
2016-05-12T17:04:48.740Z
The code appears to return a point uniformly selected from the *interior volume* of the unit sphere. This will give different results from the surface of the sphere. Is this what you intend? I ask because you mention picking a point "on" the light, which sounds like a point on the surface.
trichoplax is on Codidact now
CC BY-SA 3.0
2016-05-12T20:01:05.323Z
But.. _centerOfSphere + Vec3<float>(x, y, z) * radius;_ should give me a point on the surface of the sphere! right? cause I move from the center towards the surface of _radius_
Tarta
CC BY-SA 3.0
2016-05-12T22:18:24.223Z
If the points are already on the surface of the sphere, then multiplying by radius will give the surface of a sphere of that radius. If the points are already in the interior volume of the sphere, then multiplying by radius will put them in the interior volume of a sphere of that radius. Multiplying by, for example, radius 2, will give a sphere twice as large. However, whether the points are on the surface or in the volume of this larger sphere depends on whether they were on the surface or in the volume of the original unit sphere
trichoplax is on Codidact now
<rendering><c++><pathtracing><distribution>
Role of PDF of Uniform Random Sampling in a path tracer
Tarta
# Role of PDF of Uniform Random Sampling in a path tracer

(Asked by: username_0 on 2016-05-12)

I am trying to implement my own path tracer but before arriving to the question I want to give you a short overview:  
In the implementation of the rendering equation I use some particular technique in order to sample surfaces. For example: when one of my rays hits a diffuse surface, the next ray bouncing from that surface will be calculated using a *Cosine-weighted Random Direction*. This implies that in my rendering equation I have to take into account the PDF that this specific random distribution implies and specifically divide for this PDF in my equation. So far everything is okey.

Now, my question. I want to implement a particular technique called "next event estimation" which simply samples the lights. In order to do so, I want to pick a **random point** on my light, which is spherical, by using the following code (C++):

```
Vec3<float> randomPoint() const
{
    float x;
    float y;
    float z;

        // random vector in unit sphere
    do
    {
        x = ((float)rand() / (RAND_MAX)) * 2 - 1;           // random coordinates between -1 and 1
        y = ((float)rand() / (RAND_MAX)) * 2 - 1;
        z = ((float)rand() / (RAND_MAX)) * 2 - 1;
    } while (pow(x, 2) + pow(y, 2) + pow(z, 2) > 1);        // simple rejection sampling

    return centerOfSphere + Vec3<float>(x, y, z) * radius;
}

```

As far as I understand (please correct me if wrong), this implements a **uniformly random sampling**. Also, as far as I have read, this kind of sampling has a PDF = 1 / (b - a) . My question: do I have to use this PDF as well like I do for my *Cosine-weighted Random Distribution* or not, since it is a uniform distribution? if yes, what's the range (b - a) that PDF talks about? Thanks in advance!

**Tags:** rendering, c++, pathtracing, distribution

> **Comments:**
> **username_1**: The code appears to return a point uniformly selected from the *interior volume* of the unit sphere. This will give different results from the surface of the sphere. Is this what you intend? I ask because you mention picking a point "on" the light, which sounds like a point on the surface.  
> **username_0**: But.. _centerOfSphere + Vec3<float>(x, y, z) * radius;_ should give me a point on the surface of the sphere! right? cause I move from the center towards the surface of _radius_  
> **username_1**: If the points are already on the surface of the sphere, then multiplying by radius will give the surface of a sphere of that radius. If the points are already in the interior volume of the sphere, then multiplying by radius will put them in the interior volume of a sphere of that radius. Multiplying by, for example, radius 2, will give a sphere twice as large. However, whether the points are on the surface or in the volume of this larger sphere depends on whether they were on the surface or in the volume of the original unit sphere  

---

## 1 answers

### Answer by username_2 on 2016-05-12. Score: 7 | Accepted answer

Firstly, as @username_1 correctly pointed out, your randomPoint function calculates a point in a cube, then uses rejection sampling to return ***all*** points that are ***inside*** a unit sphere. In order to return points on a sphere, you would need to change the greater than to an equals. That said, rejection sampling is very inefficient.

A better way to sample a sphere, is to sample in spherical space, then transform to cartesian. ie:

```
float phi = PI * randf();
float theta = 2 * PI * randf();    

float x = radius * std::sinf(phi) * std::sinf(theta);
float y = radius * std::cosf(phi);
float z = radius * std::sinf(phi) * std::cosf(theta);

```

However, this is not uniform, and will cause samples to clump at the poles. To prevent this, we transform phi:

```
float phi = std::acosf(2.0f * randf() - 1.0f);
float theta = 2 * PI * randf();    

```

We can use some trig identities to make it a bit more efficient to calculate with the computer:
$$\begin{align\*}
\sin(\cos^{-1}(x)) \equiv& \ \sqrt{1 - x^2}\\
\cos(\cos^{-1}(x)) \equiv& \ x\\
\\
\phi =& \ \cos^{-1}(2\xi - 1)\\
\sin(\phi) =& \ sin(\cos^{-1}(2 \xi - 1))\\
=& \ \sqrt{1 - (2 \xi - 1)^2}\\
\cos(\phi) =& \ 2 \xi - 1
\end{align\*}$$

So the full change is:

```
float cosPhi = 2.0f * randf() - 1.0f;
float sinPhi = std::sqrt(1.0f - cosPhi * cosPhi);
float theta = 2 * PI * randf();    

float x = radius * sinPhi * std::sinf(theta);
float y = radius * cosPhi;
float z = radius * sinPhi * std::cosf(theta);

```

However, in the case of next event estimation, uniform sampling the whole sphere is inefficient, because a ray can only 'see' half of the sphere at a time. So if we generate a point on the 'back' of the sphere, the sphere will occlude the point, and your calculation will be wasted.

Instead, you generate samples in a cone, which covers the great circle of the sphere, as viewed from the starting point of the ray.

```
// Sample sphere uniformly inside subtended cone
float rand1 = randf();
float rand2 = randf();

// Compute theta and phi values for sample in cone
float distanceSquared = DistanceSquared(rayOrigin, sphereCenter);

// We use geometry to calculate the angle of the cone (aka, the maximum phi can be when we sample)
// It's easier / cheaper to use geometry to calculate sin/cos phi directly, rather than generating phi and using sin/cos
float sinPhiMaxSquared = radius * radius / distanceSquared;
float cosPhiMax = std::sqrt(1.0f - sinPhiMaxSquared);
float cosPhi = (1.0f - rand1) + rand1 * cosPhiMax;
float sinPhi = std::sqrt(1.0f - cosPhi * cosPhi);

// Phi can be anything in 2 PI
float theta = 2 * PI * rand2;

float x = radius * sinPhi * std::sinf(theta);
float y = radius * cosPhi;
float z = radius * sinPhi * std::cosf(theta);

```

**NOTE: x, y, z are in LOCAL coordinate space.**

So we need to transform to world coordinates:

```
float3 zAxis = normalize(sphereCenter - rayOrigin);
float3 xAxis;
if (std::abs(zAxis.x) > std::abs(zAxis.y))
    xAxis = float3(-zAxis.z, 0.0f, zAxis.x) / std::sqrt(zAxis.x * zAxis.x + zAxis.z * zAxis.z);
else
    xAxis = float3(0.0f, zAxis.z, -zAxis.y) / std::sqrt(zAxis.y * zAxis.y + zAxis.z * zAxis.z);
yAxis = cross(zAxis, xAxis);

float3 samplePoint = x * xAxis + y * yAxis + z * zAxis;

```

And the pdf is calculated as follows:

```
float pdf = 1 / (2 * PI * (1 - cosPhiMax));

```

The code here is heavily influenced / copied from [PBRT v3](https://github.com/mmp/pbrt-v3/blob/master/src/shapes/sphere.cpp). They have a series of classes and functions for sampling from shapes.

Finally, the pdf. In monte-carlo integration, we need to combine the pdf's of each integration we do. In path tracing, we can integrate over many many things. For example, the general rendering equation integrates the incoming light over the hemisphere, depth of field can be treated as an integration over a focal distance, etc.

For next event estimation, you explicitly split the rendering equation into two integrands, direct lighting, and indirect lighting.

Standard rendering equation:
$$ L\_{\text{o}}(p, \omega\_{\text{o}}) = L\_{e}(p, \omega\_{\text{o}}) \ + \ \int\_{\Omega} f(p, \omega\_{\text{o}}, \omega\_{\text{i}}) L\_{\text{i}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | d\omega\_{\text{i}} $$

Next Event Estimation:
$$ L\_{\text{o}}(p, \omega\_{\text{o}}) = L\_{e}(p, \omega\_{\text{o}}) \ + \ \int\_{\Omega} f(p, \omega\_{\text{o}}, \omega\_{\text{i}}) L\_{\text{i, direct}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | d\omega\_{\text{i}}\ \ + \ \int\_{\Omega} f(p, \omega\_{\text{o}}, \omega\_{\text{i}}) L\_{\text{i, indirect}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | d\omega\_{\text{i}} $$

In simple naive forward path tracing, everything is is treated as indirect light. In next event estimation, we directly calculate the direct lighting and add it to the indirect lighting. And if we hit a light, we ignore the contribution, since we're calculating the direct lighting.

Since we have two integrations, each will have its own pdf. Aka:
$$L\_{\text{o}}(p, \omega\_{\text{o}}) = L\_{e}(p, \omega\_{\text{o}}) \ \ + \ \sum\_{k=0}^{\infty } \frac{f(p, \omega\_{\text{o}}, \omega\_{\text{i}}) L\_{\text{i, direct}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | }{pdf\_{direct}}\ \ + \ \sum\_{k=0}^{\infty } \frac{f(p, \omega\_{\text{o}}, \omega\_{\text{i}}) L\_{\text{i, indirect}}(p, \omega\_{\text{i}}) \left | \cos \theta\_{\text{i}} \right | }{pdf\_{indirect}}$$

If you want to see how this is implemented, you can check out my implementation [here](https://github.com/username_2/lantern/blob/master/source/renderer/renderer.cpp#L92). Note: my light sampling is a bit more complicated, since it does multiple importance sampling. But it can be as simple as:

```
float3 Renderer::EstimateDirect(Light *light, UniformSampler *sampler, float3a &surfacePos, float3a &surfaceNormal, float3a &wo, Material *material) const {
    float pdf;
    float3 wi;

    // Sample a point on the light and get the pdf
    float3 Li = light->SampleLi(sampler, m_scene, surfacePos, surfaceNormal, &wi, &pdf);

    // Calculate the brdf value
    float f = material->Eval(wi, wo, surfaceNormal);
    return f * Li * / pdf;
}

```

> **Comments:**
> **username_2**: I believe you are correct. I will check my PBRT book when I get home and correct the code. IIRC, I need to add a square root.  
> **username_3**: Yes, @username_1 is right, sampling both angles uniformly will lead to points bunching at the poles. The simplest way to fix it is sample theta and Y uniformly, then set the XZ radius = $\sqrt{r^2 - y^2}$. Which, BTW, is exactly what you're doing in the second code snippet when you sample `cosPhi` (which is proportional to Y) uniformly. See also [this MathWorld article](http://mathworld.wolfram.com/SpherePointPicking.html).  
> **username_4**: Also in MathWorld article: Marsaglia (1972) - point in disk method should trump trig based.  
> **username_2**: I would need to do some profiling, but I'm not entirely convinced that rejection sampling would beat trig.  
> **username_0**: If I may ask you, if you feel like of course, I have another question that might definitely solve my project.. Thanks in any case for your precious help so far: http://computergraphics.stackexchange.com/questions/2461/resulting-probabilty-density-in-path-tracer-for-paths-using-next-event-estimatio  

---

