# Multigrid Poisson blending + Jacobi relaxation

## Recommended Posts

irreversible    2860
I'll start off by noting that I'm rather math-unsavvy, which is why I think I need some hand-holding with this. Regardless, I'm trying to do everything from paper (that is to say, I don't want to use anyone else's source code as a basis; nor do I want to use any libraries for this).

At this point I'm having some difficulty interpreting a whitepaper description and translating it into actual code, largely because I'm confused by the terminology and and mathematical notation. [url="http://graphics.cs.cmu.edu/projects/gradient-paint/grad.light.r2226.pdf"]Here's the paper in question[/url]. Point 6 discusses the integration step and provides a very rough outline of the algorithm in pseudo-code. I can understand much of it, but there are a number of questions I can't seem to figure out. I'll outline the algorithm the way I understand it:

1) VCycle(fh) = multigrid call per resolution level
2) f2h <- Rfh: convolve the original image with R and downsample it
3) u2h <- VCycle(f2h): recurse, downsample the downsampled image
4) uh <- Pu2h: upsample the downsampled image and blend it with the original after convolving it with the same filter kernel used above
5) uh <- Relax(uh, fh, xh0)*
6) uh <- Relax(uh, fh, xh1)*

* this here's where the confusion starts: interpreting the Jacobi relaxation and calculating the L in the formula

Here's my code for the Jacobi relaxation bit (have I understood this correctly?):

[code]
//this should work out the constants and factors - it doesn't do blending at this stage

void Relax(INOUT float * data, IN int iWidth, IN int iHeight, IN TReal fRelaxationFactor, IN int iIteration)
{
TReal h = 1.0f / (TReal)iIteration;
TReal h2 = h * h;
TReal mult = 1.0f / (3.0 * h2);
TReal m = (-8 * h2 - 4) * mult - fRelaxationFactor;
TReal e = (h2 + 2) * mult;
TReal c = (h2 - 1) * mult;

float fKernel[] = {
c, e, c,
e, m, e,
c, e, c };

//do simple NxN<->3x3 convolution
Convolve(data, iWidth, iHeight, fKernel);

//TODO: calculate the error and blend the result with the original image
}

//where
....

TReal h = 1.0f / (TReal)iIteration;
TReal h_m1 = pow(h, -1);
TReal h_m2 = pow(h, -2);

TReal fRelaxationFactor0 = -2.1532f + 1.5070f * h_m1 + 0.5882f * h_m2;
TReal fRelaxationFactor1 = 0.1138f + 0.9529f * h_m1 + 1.5605f * h_m2;

Relax(data2, iWidth, iHeight, fRelaxationFactor0, iIteration);
Relax(data2, iWidth, iHeight, fRelaxationFactor1, iIteration);
....
[/code]

Granted, the above code decidedly incorrect as I can't decipher the relaxation function line. Furthermore the capital [i]I [/i]is an unknown to me in that line.

Is my interpretation that I can reduce the Jacobi pass to a simple localized convolution correct? I mean, together with the multigrid part it should work either way; the question is rather whether it is factually correct. A further question I have is regarding blending the diffused downscale with the original, which again ties to how the relaxation function expands (particularly how the error L is handled from one iteration to the next and how h is supposed to be expressed (I've expressed it as 1/resolution, which I'm afraid is incorrect)). Can anyone help me chew through this and work out the particulars?

The result I'm striving towards uses [url="http://artis.imag.fr/Publications/2008/OBWBTS08/diffusion_curves.pdf"]gradient-domain diffusion curves[/url] as a basis so my primary goal, for now, is to get only the diffusion integration part working properly.

## Create an account

Register a new account