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

  • $\begingroup$ 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. $\endgroup$ – trichoplax May 12 '16 at 17:04
  • $\begingroup$ 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 $\endgroup$ – Tarta May 12 '16 at 20:01
  • 1
    $\begingroup$ 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 $\endgroup$ – trichoplax May 12 '16 at 22:18

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& \ \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);
    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 &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;
  • $\begingroup$ 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. $\endgroup$ – RichieSams May 12 '16 at 23:09
  • 2
    $\begingroup$ 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. $\endgroup$ – Nathan Reed May 13 '16 at 2:01
  • $\begingroup$ Fixed and updated. $\endgroup$ – RichieSams May 13 '16 at 2:36
  • $\begingroup$ Also in MathWorld article: Marsaglia (1972) - point in disk method should trump trig based. $\endgroup$ – MB Reynolds May 17 '16 at 9:18
  • $\begingroup$ I would need to do some profiling, but I'm not entirely convinced that rejection sampling would beat trig. $\endgroup$ – RichieSams May 17 '16 at 11:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.