I'm trying to preprocess a height map image with a helmholtz-type equation as described in this paper.
The equation is: $$ddx(h') + ddy(h') + y(h'-h) = 0$$
I solved for h and got: $$\dfrac{ddx(h')}{y} + \dfrac{ddy(h')}{y} + h' = h$$
I then turned this into an equation of the form $Ax=b$ where $b$ is the source image and A is a matrix with elements around the diagonal corresponding to the central finite difference approximation of the second derivative of $h'$ divided by $y$ plus an extra $1$ added on the diagonal. I then attempted to solve this sparse system with a ConjugateGradient solver.
Unfortunately, the results don't look anything like a smoothed version of the original image. For very high y, the resulting image is mostly similar to h, but for lower y it quickly degrades into noise and eventually just a black image.
I suspect part of the issue is the way I'm dealing with the boundary conditions?
Here is my source code using Eigen and C++. If anyone can guide me in the correct description I would be very grateful!