Role of PDF of Uniform Random Sampling in a path tracer 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&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;
}


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!
Firstly, as @trichoplax 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&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*}$$

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


And the pdf is calculated as follows:

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


The code here is heavily influenced / copied from PBRT v3. 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. 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 &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;
}