What I still don't know is wheter I should take the block based approach or going the scanline based way.
About that, I implemented this one or two years ago with the excellent Nicolas Capens article, (I think you both should give him more credits in your source files).
Any way, back to the topic, I used 8x4 blocks and made some extra optimizations over yours. One that I miss is that you are calculating the plane equation at every corner of the block with multiplications and it can be done incrementally adding the proper deltas to the previous plane equation evaluation.
In my current code there is no multiplication inside the blocks loop all is done using increments.
// Evaluate half-space functions
bool a00 = C1 + DX12 * y0 - DY12 * x0 > 0;
bool a10 = C1 + DX12 * y0 - DY12 * x1 > 0;
bool a01 = C1 + DX12 * y1 - DY12 * x0 > 0;
bool a11 = C1 + DX12 * y1 - DY12 * x1 > 0;
Here you really only need to sub DY12 from one block to the nect on the right, also you dont need to calculate for th four corners as the next block right plane equations are current block plane equations.
Same applies for perspective correct interpolands, and for one block an the one in the same X position on the next line.
This two things speed up a lot all the process.
Also something a it "clever" than;
// Loop through blocks
for(int y = miny; y < maxy; y += BLOCK_SIZE)
for(int x = minx; x < maxx; x += BLOCK_SIZE)
Will also help you.
On the cache think I had my zbuffer stored by tiles instead of scans and that helped a lot to the cache coherency, also using tilled textures helps.
For "pixel shaders" having a function pointer was a problem in the gp32 were I developed that. but I used a template functor and the call was for a block plus a bit field of affected pixels in a block.
Well hope this helps a bit.