Surface normals integration

Question

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?

Answer 1

Your equation basically amounts to solving $$ \begin{bmatrix} p\\ q\end{bmatrix} = \nabla z$$ for $z$.

As already noticed in amlrg's answer, this system is overdetermined and hence, one can go for a least squares solution, i.e. the solution of $$ \min_z \| \begin{bmatrix} p\\ q\end{bmatrix} - \nabla z\|^2. $$ The solution of this is obtained by solving the corresponding normal equations. As the transpose of the gradient operator is the negative divergence you obtain the equation $$ -\textrm{div}\nabla z = -\textrm{div}\begin{bmatrix} p\\ q\end{bmatrix}$$ which is nothing else than $$ -\Delta z = -\partial_x p - \partial_y q $$ and this can be solved, in principle, like any other Poisson equation (e.g. by CG with appropriate preconditioning).

Answer 2

If I am interpreting the question correctly, I think you would basically integrate this like an ODE using the euler method or perhaps something more advanced (backward euler, RK, whatever). Assuming that you know the p and q functions as givens, and are given some initial seed point z(x,y), you could "timestep" your two equations to figure out two new z-points on your surface, z(x+1,y)=z(x,y)+p and z(x,y+1)=z(x,y)+q. By repeating the process from those two new points, you can steadily fill out an entire lattice of z-points for all (x,y).

Answer 3

Supposing you know the boundary conditions, you can of course make a linear system and solve it (be it iteratively or by inverting the matrix). This would seem the most straightforward way to solve the other equations that you mentioned, too: $\Delta u = 0$, where up to a constant $\Delta \approx \begin{pmatrix}-1 & 2 & -1 & \cdots \\ 0 & -1 & 2 & -1 & \cdots\end{pmatrix}$ and so on for the 1D case (if this is not obvious, see this Wikipedia article and the related pages), adding two extra lines for the BCs. The 2D case is not much more complicated: just write the points $u_{x, y}$ as a onedimensional vector $u_i$, where e.g. $i = Nx + y$, and the difference equations for each of those points collected to a matrix. I suggest writing this problem up with a source term $\Delta u = \rho$, and making sure that for some constant boundary conditions your solution indeed is a parabola. This is a good, quick exercise for figuring out how the BCs come into play.

It was suggested by rchilton1980 that you time step, and this indeed works if you have all the initial conditions. From a quick reading of the paper, however, it seems that this might not be the case, but I must confess that I am not a 100% on this.

So for the problem you might write $i = Nx + y$, so that $z_i = z_{Nx+y} = z(x, y)$ (assuming that both directions $x$ and $y$ have $N$ points from $0$ to $N-1$, so $i$ runs from 0 to $N^2-1$).

Now $p = z(x + 1, y) - z(x, y)$ could be written in the form $p_i = z_{N(x+1)+y} - z_{Nx+y} = z_{i+N} - z_{i}$. This is valid for all $i<N^2-1-N$, so it can be written in the matrix form $p = Az$, where $A = \begin{pmatrix}-1 & 0 & \cdots & 0 & 1 & 0 & \cdots \\ 0 & -1 & 0 & \cdots & 0 & 1 & 0 & \cdots \\ \vdots \end{pmatrix}$, where the 1 and -1 are separated by a distance $N$.

Now using the equation for $y$, $q_i = z_{Nx+(y+1)} - z_{Nx+y} = z_{i+1} - z_{i}$, which is valid for $y \neq N-1$, you can get more rows to add to $A$.

So here you have "too many" lines in the matrix, and therefore the authors of the paper use least squares (or the pseudoinverse of $A$, if you like) to solve for $z$.