[contravariant]

I don't think you really need to carve individual pixels. If you've got the logic worked out to ensure the light falling at (x,y) ends up at (u,v) then the part you're missing (assuming we're dealing with sunlight for now) is a function f(x,y) = (u,v) that transforms a constant density into the density corresponding to the image you desire.

Using the change of variables formula this basically means that we want h(f^(-1)(x,y)) |det Df| to be constant. Which is quite easy in 1D (it's just the inverse cumulative density function), but significantly trickier in 2D.

In 2D the problems seems to be underdetermined. One solution would be to first solve the horizontal problem for each row and then solve the vertical part for the total densities of each row. Or the other way around. There might be a way to make this optimal in some sense, but I'm not quite sure what to optimize for.