Here is the code. It shows a cylinder with grey / sides and a small area of emitting light.
The scene looks boring and ugly but indirect lighting / color bleeding shows up - assuming i've done nothing wrong it's ground truth. Because there can't happen any occlusion, i simply set visibility to 1. You could model a cornell box without the boxes in the middle (so just the walls and the light) to get a nicer image. Adding brute force raytracing to support occlusion should be a matter of minutes.
The main problem here is that form factor calculation is optimized and it's impossible to understand / explain the math. I'll do this tomorrow...
But you already can see how simple it is.
struct Radiosity
{
typedef sVec3 vec3;
inline vec3 cmul (const vec3 &a, const vec3 &b)
{
return vec3 (a[0]*b[0], a[1]*b[1], a[2]*b[2]);
}
struct AreaSample
{
vec3 pos;
vec3 dir;
float area;
vec3 color;
vec3 received;
float emission; // using just color * emission to save memory
};
AreaSample *samples;
int sampleCount;
void InitScene ()
{
// simple cylinder
int nU = 144;
int nV = int( float(nU) / float(PI) );
float scale = 2.0f;
float area = (2 * scale / float(nU) * float(PI)) * (scale / float(nV) * 2);
sampleCount = nU*nV;
samples = new AreaSample[sampleCount];
AreaSample *sample = samples;
for (int v=0; v<nV; v++)
{
float tV = float(v) / float(nV);
for (int u=0; u<nU; u++)
{
float tU = float(u) / float(nU);
float angle = tU * 2.0f*float(PI);
vec3 d (sin(angle), 0, cos(angle));
vec3 p = (vec3(0,tV*2,0) + d) * scale;
sample->pos = p;
sample->dir = -d;
sample->area = area;
sample->color = ( d[0] < 0 ? vec3(0.7f, 0.7f, 0.7f) : vec3(0.0f, 1.0f, 0.0f) );
sample->received = vec3(0,0,0);
sample->emission = ( (d[0] < -0.97f && tV > 0.87f) ? 35.0f : 0 );
sample++;
}
}
}
void SimulateOneBounce ()
{
for (int rI=0; rI<sampleCount; rI++)
{
vec3 rP = samples[rI].pos;
vec3 rD = samples[rI].dir;
vec3 accum (0,0,0);
for (int eI=0; eI<sampleCount; eI++)
{
vec3 diff = samples[eI].pos - rP;
float cosR = rD.Dot(diff);
if (cosR > FP_EPSILON)
{
float cosE = -samples[eI].dir.Dot(diff);
if (cosE > FP_EPSILON)
{
float visibility = 1.0f; // todo: trace a ray from receiver to emitter and set to zero if any hit (or use multiple rays for accuracy)
if (visibility > 0)
{
float area = samples[eI].area;
float d2 = diff.Dot(diff) + FP_TINY;
float formFactor = (cosR * cosE) / (d2 * (float(PI) * d2 + area)) * area;
vec3 reflect = cmul (samples[eI].color, samples[eI].received);
vec3 emit = samples[eI].color * samples[eI].emission;
accum += (reflect + emit) * visibility * formFactor;
}
}
}
}
samples[rI].received = accum;
}
}
void Visualize ()
{
for (int i=0; i<sampleCount; i++)
{
vec3 reflect = cmul (samples[i].color, samples[i].received);
vec3 emit = samples[i].color * samples[i].emission;
vec3 color = reflect + emit;
//float radius = sqrt (samples.area / float(PI));
//Vis::RenderCircle (radius, samples.pos, samples.dir, color[0],color[1],color[2]);
float radius = sqrt(samples[i].area * 0.52f);
Vis::RenderDisk (radius, samples[i].pos, samples[i].dir, color[0],color[1],color[2], 4); // this renders a quad
}
}
};