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