Hello,

Hopefully some of you can help me out! I have to use the Full Multigrid Algorithm from chapter 20.6.

For my master's thesis I have to reimplement a paper about gradient domain high dynamic range compression. In this paper the gradient field of a luminance image is manipulated by attenuating the magnitudes of large gradients. Then a low dynamic range image is obtained by solving a Poisson equation on the modified gradient field. The paper uses the Full Multigrid Algorithm from NR3 for this.

The Poisson equation:
V²I = div G
with
V = triangle upside down
I = image
G = gradient field

The Laplacian of the equation is approximated as follows:
V²I(x, y) = I(x + 1, y) + I(x - 1, y) + I(x, y + 1) + I(x, y - 1) - 4 * I(x, y)

The gradient is approximated using forward difference:
VH(x, y) = (H(x + 1, y) - H(x, y), H(x, y + 1) - H(x, y))

div G is approximated using backward difference:
div G = Gx(x, y) - Gx(x - 1, y) + Gy(x, y) - Gy(x, y - 1)

from the paper:
"the difference scheme yields a large system of linear equations - one for each pixel in the image, but the corresponding matrix has only five nonzero elements in each row, since each pixel is coupled only with its four neighbours."

Questions:

I have not the slightest idea how to put these equations into the input matrix for the mglin function :( Which format is necessary for the algorithm to perform correctly?

The comment in the code listing mglin.h says "On input u[0..n-1][0..n-1] contains the right-hand side p,..."

Is p = div G ?

What does the output look like? Does each row of u contains the new luminance value for the corresponding pixel?

Thanks in advance for any help you can give me.

The paper:
http://www.cs.huji.ac.il/%7Edanix/hdr/hdrc.pdf

PS: If you help me, I'll send you a post card from my home town :)

I guess what I am looking for is the appropriate Numerical Recipes Example Program xmglin.c as in V2.08 of the Numerical Recipes source code.

Or can I assume that it hasn't change much in nr3?

I guess what I am looking for is the appropriate Numerical Recipes Example Program xmglin.c as in V2.08 of the Numerical Recipes source code.

Or can I assume that it hasn't change much in nr3?

First of all, nr3 is C++, not C. A number of things are much easier to use for C++ programmers, but might represent a "learning opportunity" for those with no C++ experience.

Secondly, the old 1-base array stuff in NR2 (and older) has mercifully been laid to rest. Maybe this is to the chagrin of the legions of Fortran programmers out there (who only occasionally dabble in C), but it's a great relief to denizens of the world of C/C++. Most of the C++ examples of NR2 had already made that transition.

You happened to pick something that is still a function in nr3, and not a C++ object (so using mglin doesn't require any substantial changes in the driving program; only a couple of minor tweaks).

If you have the nr3 source CD or download, copy xmglin.cpp from the legacy/nr2/CPP_211/examples directory into your nr3 project directory and make the following changes. (Leave everything else as-is.)


/*
* You need these headers and none of the old stuff.
* Change to whatever path you need.
*/
#include "../code/nr3.h"
#include "../code/mglin.h"

// Driver for routine mglin

int main(void)
{
const int NSTEP=4,JMAX=33;
int i,j,midl=JMAX/2;

/* change Mat_DP to MatDoub */
/* Mat_DP f(JMAX,JMAX),u(0.0,JMAX,JMAX); */

MatDoub f(JMAX, JMAX);
MatDoub u(JMAX, JMAX, 0.0); // initialize to zeros

.
.
.
/* The following line was NR::mglin(u, 2) */
Mglin(u, 2);
.
.
.


Output:

MGLIN solution:
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0000 0.0000
0.0000 -0.0001 -0.0001 -0.0002 -0.0002 -0.0002 -0.0001 -0.0001 0.0000
0.0000 -0.0001 -0.0002 -0.0003 -0.0005 -0.0003 -0.0002 -0.0001 0.0000
0.0000 -0.0001 -0.0002 -0.0005 -0.0014 -0.0005 -0.0002 -0.0001 0.0000
0.0000 -0.0001 -0.0002 -0.0003 -0.0005 -0.0003 -0.0002 -0.0001 0.0000
0.0000 -0.0001 -0.0001 -0.0002 -0.0002 -0.0002 -0.0001 -0.0001 0.0000
0.0000 -0.0000 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

Test that solution satisfies difference equations:
0.0001 -0.0001 -0.0004 0.0008 -0.0004 -0.0001 0.0001
-0.0001 -0.0002 -0.0008 0.0028 -0.0008 -0.0002 -0.0001
-0.0004 -0.0008 -0.0027 0.0048 -0.0027 -0.0008 -0.0004
0.0008 0.0028 0.0048 2.0171 0.0048 0.0028 0.0008
-0.0004 -0.0008 -0.0027 0.0048 -0.0027 -0.0008 -0.0004
-0.0001 -0.0002 -0.0008 0.0028 -0.0008 -0.0002 -0.0001
0.0001 -0.0001 -0.0004 0.0008 -0.0004 -0.0001 0.0001


Dave

Hello Dave,

thanks for your answer. I didn't even cross my mind to check for one-based loops.

I applied the changes you suggested and compared it to the output you posted. At least now I know that the algorithm is working.


Thanks again,
Christian