2D Discrete laplacian (del2) in C++

Question

I am trying to figure out how to port the del2() function in matlab to C++.

I have a couple of masks that I am working with that are ones and zeros, so I wrote code liket his:

for(size_t i = 1 ; i < nmax-1 ; i++)

{
    for(size_t j = 1 ; j < nmax-1 ; j++)

    {
        transmask[i*nmax+j] = .25*(posmask[(i+1)*nmax + j]+posmask[(i-1)*nmax+j]+posmask[i*nmax+(j+1)]+posmask[i*nmax+(j-1)]);

    }
}

to compute the interior points of the laplacians. I think according to some info in "doc del2" in matlab, the border conditions just use the available info to compute, right? SO i guess I just need to write cases for the border conditions at i,j = 0 and nmax

However, i would think these values from the code I have posted here would be correct for the interior points as is, but it seems like the del2 results are different!

I dug through the del2 source, and I guess I am not enough of a matlab wizard to figure out what is going on with some of the code for the interior computation

Answer 1

You can see the code of del2 by edit del2 or type del2. Note that del2 does cubic interpolation on the boundaries.

Answer 2

The problem is that the line you have there:

transmask[i*nmax+j] = .25*(posmask[(i+1)*nmax + j]+posmask[(i-1)*nmax+j]+posmask[i*nmax+(j+1)]+posmask[i*nmax+(j-1)]);  

isn't the discrete Laplacian at all.

What you have is (I(i+1,j) + I(i-1,j) + I(i,j+1) + I(i,j-1) ) / 4

I dont' know what this mask is, but the discrete Laplacian (assuming the spacing between each pixel in each dimension is 1) is:

(-4 * I(i,j) + I(i+1,j) + I(i-1,j) + I(i,j+1) + I(i,j-1) )

So basically, you missed a term, and you don't need to divide by 4. I suggest going back and rederiving the discrete Laplacian from its definition, which is the second x derivative of the image plus the second y derivative of the image.

Edit: I see where you got the /4 from, as Matlab uses this definition for some reason (even though this isn't standard mathematically).

Answer 3

I think that with the Matlab compiler you can convert the m code into C code. Have you tried that?

I found this link where another methot to convert to C is explained.

http://www.kluid.com/mlib/viewtopic.php?t=337

Good luck.