the equation right hand side is wrong

[Citronnier]

I went through and trialed the code and found an error:

// For the Poisson equation the divergence of the guidance field is necessary.
        cv::Mat vxx, vyy;
        cv::Mat kernelx = (cv::Mat_<float>(1, 3) << -0.5, 0, 0.5);
        cv::Mat kernely = (cv::Mat_<float>(3, 1) << -0.5, 0, 0.5);
        cv::filter2D(vx, vxx, CV_32F, kernelx);
        cv::filter2D(vy, vyy, CV_32F, kernely);
        
        cv::Mat f = vxx + vyy;

From clone.cpp line 167-174

Here the divergence of the guidance field (i.e. Laplacian of g) is added to the right hand side of the linear system. However, the addend should be the sum of all (g_p - g_q), which is not the same with the Laplacian.

[cheind]

Hi, I'm not sure I can follow (and I'm not sure I can remember everything about my code :) Do you have references for your suggestion? Are you referring to the addition of the second order derivatives in x and y direction vxx + vyy? AFAIR the divergence is the addition of those derivatives.

[Citronnier]

Hi, I'm referring to the addition of the second order derivative vxx+vyy, and this sum is, of course, the divergence of the guidance field. The problem is it should not be the right hand side addend of the linear system here:

poisson-solver.cpp, line 203-204

// Add f to rhs.
                rhs.row(pid) += Eigen::Map<Eigen::VectorXf>(f.ptr<float>(p.y, p.x), channels);

The correct addend should be the sum of the projection of the guidance field on all vector pq. You could refer to the paper you implemented this algorithm from mentioned in readme. The addend is specified therein as equation (11) in section 3.

I'll try to submit a pull request if it helps :)

[cheind]

Ah, now I think I know what you are refering to. Equation 11 however refers to the numerical estimation of the gradient (see equation 9) in an image through subtracting neighboring pixels in the source image. This is already accounted for in the kernels for the first order and second order derivatives

cv::Mat kernelx = (cv::Mat_<float>(1, 3) << -0.5, 0, 0.5);

It differs from equation 11 in that the gradient for a pixel p in x-direction is computed symmetrically from image g as follows -g((p.x-1), p.y)*0.5 + g((p.x+1), p.y)*0.5. Not sure why I introduced the scaling of 0.5 though.

In case you are concerned about the projection on the oriented edge, right beneath equation 11 it says that the gradient of the source image is already considered the projection on the oriented edge, which is why, I think, I haven't considered any other projection.

It seems in the original paper they suggest to compute the discrete gradient in every direction, i.e for every p take all q pixels in the neighborhood of p. That's probably something I haven't followed very strictly.

[Citronnier]

Even without the 0.5 scaling, it would still be incorrect. The reason is the neighbors of p are on different sides of p, thus the sum of all these projections, according to equation 11, is 4*g(p) - \Sigma_q(g(q)). Note that for the 'northern' and 'eastern' neighbor, this projection is actually the opposite of the gradient at p. While in the source code, it is the value of the variable "f" at point p, which is the value of the Laplacian.

I wonder how your implementation turns out, as it doesn't really work when I tried, which led me to trace the error and found what I believe is the problem. After I made this correction, the code works just fine. :)

[cheind]

Yes you are right! Interesting, never experienced an issue in my tests.

[NengQian]

Actually, the direction of 'northern' and 'eastern' neighbor is still correct by taking that kernel:
cv::Mat kernelx = (cv::Mat_<float>(1, 3) << -0.5, 0, 0.5);
because:
f = vxx + vyy
and,
vxx(p.x, p.y) = 0.5*(-vx(p.(x-1), p.y) + vx(p.(x+1), p.y))
since,
vx(p.(x-1), p.y) = 0.5(-g(p.(x-2), p.y) + g(p.x, p.y))
vx(p.(x+1), p.y) = 0.5(-g(p.x, p.y) + g(p.(x+2), p.y))
Therefore
vxx(p.x, p.y) = 0.5*( -0.5(-g(p.(x-2), p.y) + g(p.x, p.y)) + 0.5( -g(p.x, p.y) + g(p.(x+2), p.y)) )
vxx(p.x, p.y) = 0.25*( g(p.(x-2), p.y) - g(p.x, p.y) - g(p.x, p.y) + g(p.(x+2), p.y) )
vxx(p.x, p.y) = 0.25*( - 2*g(p.x, p.y) + g(p.(x-2), p.y) + g(p.(x+2), p.y) )
With the same procedure,
vyy(p.x, p.y) = 0.25*( - 2*g(p.x, p.y) + g(p.x, p.(y-2)) + g(p.x, p.(y+2)) )
Therefore,
f = 0.25*( - 4*g(p.x, p.y) + g(p.x, p.(y-2)) + g(p.x, p.(y+2)) + g(p.(x-2), p.y) + g(p.(x+2), p.y) )

The direction is still correct, the only issue is the 0.25 scale and, it is not exactly the 4-neighbor of p but the +2 and -2 neighbor. But since the 4-neighbor itself is just an approximation of discrete gradient, the result from this code is quite similar to the paper.
A better way to calculate this gradient and laplacian staff is by setting different kernels for gradient and laplacian computation. For example, the gradient kernel should be:
cv::Mat kernelx = (cv::Mat_<float>(1, 3) << 0 , -1, 1);
and the laplacian kernel should be:
cv::Mat kernelx = (cv::Mat_<float>(1, 3) << -1 , 1, 0);
For the y direction, it is the same.
By doing this,
f = 4*g(p) - sum of 4-neighbor of p

ref: the implementation of seamless_cloning_impl in OpenCV: