2 minutes ago, matt77hias said:

basic and practical knowledge

that's what I need! Thanks! Did you read this book?

2 hours ago, matt77hias said:

Seems to me for offline (i.e. not real-time) rendering.

To make that in real time, we can use the Instant Radiosity (from Keller (in 1997))
And Final gathering - that much faster than photon-based GI
 

2 hours ago, matt77hias said:

Are you sure this uses the radiosity algorithm?

It based on site of Hugo Elias who documented his radiosity light mapper
http://web.archive.org/web/20160311085440/http://freespace.virgin.net/hugo.elias/radiosity/radiosity.htm
 

2 hours ago, Michael Davidog said:

Can you please advice me some good books/resources?

The chapter in OpenGL Super Bible about compute shaders - it's great and covers all basics, trains the kind of thinking you need to have for massive parallelization. To me this was enlightening and the rest then is mostly technical details you get while doing and from experience. So even if you want to use another language, read it. (OpenCL 1.x, Compute shader in OGL/VK / DX11or12 - it's all the same)

I recommend using OpenCL to start. It's by far the easiest way on the API side. VK/DX12 is too complicated to get going, but you can port your OpenCL 1.x kernels to them later. (I use both OpenCL and VK code paths, both using the same shader code but preprocessed to adapt syntax, VK is almost two times faster.)

2 hours ago, Michael Davidog said:

This is the same principle as "patches" but clusters are hierachical and bigger?

Yes, you can divide your surfaces into surfels, disks, quads, polygons... what ever you want. For form factor calculation you typically use only their area and ignore the shape. Error becomes noticeable only if surfaces come very close together, but at this point you subdivide before this becomes a problem. Hierarchies are the main idea to make things faster. Distant objects can be approximated by a single piece of area, as anywhere else where LOD is used.

I'll write a simple simulator in some hours in another post. This should answer most of your questions...

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

};