Numerical test of integrability with applications in computer vision In my computer vision course, we are working on extracting a 3D surface from a chain of 2D images taken under several conditions. This procedure is known as Photometric stereo. Prior to extracting the surface $f(x,y)$ numerically, we have to perform what's called a test of integrability :
$$
\frac{\partial^{2}f}{\partial x\partial y}=\frac{\partial^{2}f}{\partial y\partial x}
\tag{1}$$
This sanity test can be tested numerically by discretizing the partial derivative so that it accepts small values (and thus we have to accept that to recover $f$, we recover it starting from small values). Numerically this translates to :
$$
\frac{\partial\displaystyle\left(\frac{\alpha(x,y)}{\gamma(x,y)}\right)}{\partial y}-\frac{\partial\displaystyle\left(\frac{\beta(x,y)}{\gamma(x,y)}\right)}{\partial x}\approx0
\tag{2}$$
where :
$$
\frac{\partial f}{\partial x}=\frac{\alpha(x,y)}{\gamma(x,y)}\qquad\text{and}\qquad\frac{\partial f}{\partial y}=\frac{\beta(x,y)}{\gamma(x,y)}
\tag{3}$$
essentially where $(\alpha(x,y),\beta(x,y),\gamma(x,y))$ are the coordinates of the measured value of the unit-normal vector at some point on the surface $f$.

In short, how can we numerically check if $(2)$ is valid? what algorithm or procedure should we do to discretize the partial differentiation?

I am not familiar with advanced concepts such as discretizing partial differential equations. However, I am familiar with numerical integration techniques.
Source of the discussion: Forsyth, Ponce, Computer Vision: A Modern Approach (Second Edition), 2012
 A: Briefly, the suggestion to consider numerical differentiation is vulnerable against errors in the measurements and this test of integrability is likely to produce many false negatives. Instead we should turn to numerical integration.
The situation is as follows. We are given a $C^1$ vector field $F : \Omega \rightarrow \mathbb{R}^2$ where $\Omega \subseteq \mathbb{R}^2$. We wish to assert that $F$ is integrable, i.e., that there exists a function $\phi : \Omega \rightarrow \mathbb{R}^2$ such that $$F(x,y) = \left(\frac{\partial \phi}{\partial x}(x,y), \frac{\partial \phi}{\partial y}(x,y)\right).$$ Now, if $\Omega$ is an open simply connected subset of $\mathbb{R}^2$, the $F = (f,g)$ is integrable if and only $$\forall (x,y) \in \Omega \: : \: \frac{\partial f}{\partial y}(x,y) = \frac{\partial g}{\partial x}(x,y) \label{a} \tag{1}.$$ What can we do when $f$ and $g$ cannot be subjected to direct mathematical analysis, but can only be measured? Firstly, it is clear that we can only examine finitely many points $(x,y) \in \Omega$. Certainly, derivatives can be approximated using finite differences. In particular, we have the standard space central approximation of the first order derivative, i.e.,
$$ \frac{\partial f}{\partial y}(x,y) = \frac{f(x,y+h)-f(x,y-h)}{2h} + O(h^2), \quad h \rightarrow 0, \quad h \not = 0$$ and $$ \frac{\partial g}{\partial x}(x,y) = \frac{f(x+h,y)-f(x-h,y)}{2h} + O(h^2), \quad h \rightarrow 0, \quad h \not = 0.$$
Here the underlying assumptions are that $f$ and $g$ are of class $C^3$. When $f$ and $g$ can be evaluated at any point, then it is possible to reliably estimate the error and determine if equation $\eqref{a}$ is satisfied to the limit of machine precision. The underlying technique is now as Richardson extrapolation. Deciding when the error estimates are reliable can be done by observing Richardson's fraction. A small sample of this technique is shown in this answer to a question on numerical integration and this answer to a question on sequences.
However, while this is very useful, the finite differences are very vulnerable to measurement errors. Suppose that we cannot measure $f$ exactly, but our instruments return $$\hat{f}(x,y) = f(x,y)(1 + \delta_{(x,y)})$$ for what is hopefully always a small value of $\delta_{(x,y)}$. Then the computed approximation will satisfy
$$ \frac{\hat{f}(x,y+h) - \hat{f}(x,y-h)}{2h} =  \frac{f(x,y+h) - f(x,y-h)}{2h} + \frac{f(x,y+h)\delta_1 - f(x,y-h)\delta_2}{2h}$$
This is disastrous! The measurement error is divided with $h$ and sooner rather than later the computed approximation of the derivative will be utterly useless.
What can we do instead? We want to be sure that $F$ is integrable. We will never certain, but we can look for evidence that supports this conjecture. We compute $$ \oint f(x,y)dx+g(x,y)dy$$ along a selection of closed paths. If $F$ is integrable, then result should be zero. We can use numerical integration to evaluate the integrals. There will still be two contributions to the total error, namely the discretization error and the measurement error. The discretization error will decay as the step size is reduced. However, the measurement error is not amplified as before. If the absolute value of computed value of the integral is less than our error bound, then we have found nothing to suggest that $F$ is not integrable. I suspect that this is as far we can go.
