You can use the form factor of a disc pointing towards sample position:
qVec3 diff = ePos - rPos; // emitter (disc) - receiving sample position
float dR = dot(rDir, diff); // rDir = receiver normal
if (dR > FP_EPSILON)
{
float dE = -dot(eDir, diff); // eDir = disc normal
if (dE > FP_EPSILON)
{
float sqD = dot(diff, diff);
float formFactor = (dR * dE) / (sqD * (PI * sqD + discArea)) * discArea;
qVec3 light = eLuminosity * eCol;
sampleReceivedLight += formFactor * light;
}
}
If that works for you, you could try to optimize for a sphere.
The code snippet below should help to understand the math.
It projects a sphere to the hemisphere and then to the normal plane of the sample (We want to know the area if the projected shape).
The disc is approximated then by a simple dot product (assuming orthonormal projection and ignoring perspective).
A disc intersecting the normal plane is also approximated only with a dot product -> some error when disc is large and close to sample.
inline void ProjectDisc (
qVec3 &ray, float &hemiSphereProjectedAreaByPI,
const qVec3 &relPos, const qVec3 &discDir, const qMat3 &local) // todo: remove matrix and give localRay instead
{
// Note: Disc direction is assumed to point towards the sample position!
float sqDist = relPos.SqL() + FP_TINY;
ray = local.Unrotate(relPos);
hemiSphereProjectedAreaByPI = -(ray[2] * discDir.Dot(relPos)) /
(sqDist * (PI * sqDist + discDir[3])) * discDir[3]; // discDir[3] = disc area
ray /= sqrt(sqDist);
// unoptimized math:
//
// float unitRad = radius / sqrt(radius*radius + dist*dist); // sin(solidAngle); // radius / Hypotenuse
// float unitArea = unitRad*unitRad*PI;
// planeArea = unitArea * ray[2];
// planeArea *= fabs (discDir.Dot(ray));
// planeAreaByPI = planeArea / PI;
//
// float solidAngle = atan (radius * ooDist);
// RenderCone (qVec3(0,0,0), (qVec3&)ray, solidAngle, 0,0.5,1, 1,0,0);
// RenderCircle (sin(solidAngle), (qVec3&)(ray * cos(solidAngle)), (qVec3&)ray, 0,0.5,1);
}