I've decided to ditch the the lightmapping library I was using previously, since it ran unbelievably slowly and it had some strange artifacts that I didn't want to deal with. In it's place, I've just gone with integrating embree and doing ray-tracing on the CPU. The basic procedure for what I'm doing so far is this:

direct_buff is a black image

for each texel in the lightmap:
    Cast a ray originating from the texel, in the direction of the directional light source
    if the ray was not occluded:
        influence = max(dot(light_dir, texel_norm), 0)
        direct_buff[texel] = influence * light_color
    endif
endfor

indirect_buff = copy(direct_buff)

for each texel in the lightmap:
    Compute a TBN matrix for this texel
    Create 128 random unit vectors on the hemisphere represented by the TBN matrix
    Cast a ray originating from the texel, in the direction of each vector
    for each ray that hit something:
        Sample albedo value stored at the hit coordinate
        Sample direct lighting value stored at the hit coordinate
        influence = dot(ray_dir, texel_norm) / (distance * distance)
        indirect_buff[texel] += influence * albedo * direct
    endfor
endfor


return indirect_buff

It's not ideal, since I'm using a physically based rendering pipeline but this has no consideration for any of the usual PBR stuff, other than albedo. However, as long as it looks decent I'm fine with it.

So I whipped up a simple scene similar to the Cornell Box, when unlit it looks like this:

Computing just the direct lighting value (with the light source outside the box) and multiplying with albedo, I get this:

Ok, looks pretty good. Then as a test, I output just the indirect influence (not multiplying by albedo, also greatly reducing brightness):

Makes sense, I think. Surfaces that are close to other surfaces get more indirect ray hits and influence than those that aren't. Should it be that bright around the corners though? I'm not sure...

Anyway, here's what it looks like with the final output terms (and multiplying by albedo):

Ehh, yeah that doesn't look quite right. Why is there such an insane amount of light in the corners, even where there's no direct lighting? I'm getting some nice color bleeding from the red wall onto the left cube, but that doesn't distract from the rest of it. Here's another shot, taken from behind the right box:

Yeah, there definitely shouldn't be that much light there. I'll check my math, but if anyone could give me any tips for where I might be going wrong, I'd appreciate it. I can post the code, but it's pretty gross at the moment.

Thanks!

you need to multiply your resulting sum by 1 / NumSamples, and also multiply each sample by 1 / p(x)

Wouldn't that make the expression:

float3 result = sum * (1.0f / (numSamples * 2 * Pi));

Edit: Reading that again, I'm actually not sure. There is definitely too much light when it's not done as sum / Area, though.

Double Edit: Yeah, you're definitely right. The reason I was getting too much light was I was dividing by the number of batch samples (8), rather than the number of total samples (128).

Anyway, thank you! It looks much better now.

I'm still getting slightly too much light in areas where you wouldn't expect (around the corners even in the dark), I'm gonna check my direct buffer sampling code to make sure I didn't make any mistakes there.

Update: I was accidentally swizzling the barycentric coordinates for the ray hit location. It now looks like this:

I looks like cosine falloff is still missing. light from normal surface direction contributes more than light from a tangent direction, the weight is simply dot(normal, ray).
You see the discontinuities at the ground in front bottom corners of the cube, they're caused by swithing between too much / no indirect light coming from the lit cube.

You can either multiply the incoming light by the dot product (but i'm not 100% sure),
or change the distribution of the random ray directions so that there are more rays in normal direction (needs less rays for similar quality).

Edit:
I think to achieve the latter you can do it this way:
Generate random 2D samples in a unit circle with uniform distribution. (can also use a grid and reject points outside the circle for a quick test)
Make the samples 3D and unit vectors by setting z = 1 - (x*x + y*y).
Use those 3D samples as your ray directions.

As a result the cosine falloff should happen automatically in a similiar way the attenuation happens.

be very careful with your Pi's. It's really easy to miss one (or have one too many), which will subtly mess up your results

Agree to that - main source of bugs. Especially true for ray tracing because rays don't have area.
it's easy to get confused from references talking either of an area based approach (e.g. subdividing patches radiosity solver)
or an area less approach (discretization with rays, where we either need to weight or distribute properly).

I like the simple example of projecting a spherical light source to the hemisphere and then to the sample plane (assumes the sphere is above the sample plane so completely visible).
This helps to understand the math, is simple enough to explain diffuse lighting, and also serves as a testcase to proof the correctness when aproximating the same with a number of rays:

float unitRad = sphereRad / sphereDist; // disc radius of sphere projected to unitsphere (usually referred as hemisphere)
//float solidAngle = asin (unitRad); // angle of the cone from sample position enclosing the sphere

float unitArea = unitRad*unitRad * PI; // area of unitsphere projection

float planeArea = unitArea * sampleNormal.Dot(rayToSphere); // projecting from unitsphere down to sample plane (the dot product i've talked about)

float const areaOfUnitCircle = PI; // intersection area of unitsphere and sample plane
float formFactor = planeArea / areaOfUnitCircle; // how much of the samples 'view' is occupied by the sphere

sampleReceivedLight += formFactor * sphereOutGoingLight;