I have an image, which is a rectangular array of pixel values. I have the Laplacian of the image, which is computed as $\Delta I = I_{xx} + I_{yy}$ where $I_{xx},I_{yy}$ refer to the second derivatives of the image in the x and y directions respectively. The image derivatives are computed as: $I_{x}(i,j) = I(i+1,j) - I(i,j)$
$I_{y}(i,j) = I(i,j+1)-I(i,j)$
I am trying to obtain the original image $I$. Is this possible? What information is needed (in addition to the Laplacian of the image)?
You need the boundary conditions. The Dirichlet condition ($0$ value on the boundary) is not a natural thing to impose on an image, so the Neumann condition ($0$ derivative on the boundary) is preferred.
One can solve Discrete Poisson equation as a big linear system (part of which is boundary conditions). The system takes some effort to set up for 2D problems. There is also an approach via discrete Fourier transform: see Poisson equation with Neumann boundary conditions which has Python code.
Two things to check:
The computation of Laplacian, as you described it, looks bad to me. One shouldn't compute the second derivative by applying the same first-derivative formula twice. Applying $I(i+1,j)−I(i,j)$ twice gets you $I(i+2,j)-2I(i+1,j)+I(i,j)$ while the correct expression for $xx$-derivative at $(i,j)$ is $$I(i+1,j)-2I(i,j)+I(i+1,j)$$
Also check how the boundary was treated when the Laplacian was computed. Is $\Delta u$ only given at the interior points of the grid, or was $\Delta u$ also computed at the boundary, using "ghost points" outside of the grid? This affects the setup of the linear system. In the latter case, the boundary condition is already built into the computation of $\Delta u$.