I am trying to reconstruct a 3D surface given the normals of the unknown surface. Reading through this paper on section 4 they say
[...] denote the surface by $z(x,y)$. The directions of the normals are given by $n(x,y) = (p,q,-1)^T$, with $p = z_x$ and $q = z_y$ which denote the partial derivatives of $z$ with respect to $x$ and $y$ respectively. The recovered scaled surface normal $(n_x, n_y, n_z)$ roughly satisfies $$(n_x,n_y,n_z)=\frac{\rho}{\sqrt{p^2+q^2+1}}(p,q,-1)$$
From there, I understand the fact that $p = \frac{-n_x}{n_z}$ and $q = \frac{-ny}{nz}$ and since $p$ and $q$ are the partial derivatives of $z$ with respect to $x$ and $y$ respectively we can approximate them by
$$p = z(x+1,y)-z(x,y)$$ $$q = z(x,y+1)-z(x,y)$$
Now, what I don't really get is how to solve for $z$. All the PDEs I've seen so far are like $u_x+u_y = 0$ or $\Delta u = u_{xx}+u_{yy} = 0$ or something of that style. But here I don't see how can I integrate, this seems to be a system of equations (I'm not sure about this).
Can someone explain me how can I solve this numerically? should I just build 2 linear systems? since there is a term $z(x,y)$ in both equations can I equate them so to have only one and then solve the linear system? I believe there's something I am not considering/seeing/forgetting and it could be solved easily.
EDIT
I found some useful information in another forum. It is very related to what I am trying to do. Given the gradient field $\vec N$, in my case the normals (correct me if I'm wrong) I want to find a function $S(x,y,z(x,y))$ such that $\vec N = \nabla S$ (pretty much the same problem stated in the forum). I understand the way it could be done analytically, but I still cannot get my head around when going into the discrete formulation which is what I have.
Analytically, first I integrate with respect to $x$ and that result I need to get the derivative in order to plug it into the equation where I need to integrate over $y$. How to connect these ideas to a numerical integration scheme?