Introduction

I have recently started my thesis on 3D surface reconstruction. One way to do so is to define a surface as the zero level set of a function $f : \mathbb{R}^3 \to \mathbb{R}$, and then find some way to build this function out of a set of points $\mathcal{P} \subset \mathbb{R}^3$. NCH is a method that allows you to define such a function starting from a point cloud along with their normals, which tell you which way is the inside of the surface.

The NCH algorithm is ashtonishingly simple: it took me an afternoon to implement and see it working with no previous experience in 3D whatsoever (using the PCL, of course!). However, it has two big shortcomings: finding $f$ is $\mathcal{O}(n^2)$, and evaluating $f$ at any point is $\mathcal{O}(n)$; this makes it a really expensive algorithm, considering that most surfaces are well into the many thousands of points.

If you read the paper, it will be clear that it is so simple there is no obvious way to reduce the amount of computation it has to do while also preserving the good theoretical properties it has. The objective of my thesis is quite undefined for now, but it will be about simplification: finding ways to cut down the number of points that you have to feed to the algorithm, in order to make it feasible to run it on larger point clouds.

Whatever I end up doing, I will certainly be using this blog as an archive of proofs and whatever interesting things I learn throughout, with the hope that by the day I have to sit down and write my thesis I will at least have some stuff already written down, or at least a clearer picture.

What follows is the first proof I made, which fills a gap the paper does not, namely, proving that the algorithm is insensitive to rotation and translation of the points it is given.

Definitions

Let $\{p_1, …, p_N\}$ be a finite set of points with associated unit length orientation vectors $\{n_1, …, n_N\}$. Let

$$f(x) = \max_{1 \leq i \leq N} f_i(x)$$

Be the NCH signed distance function (the one whose 0 level set is the surface we want to find), with

$$f_i(x) = n_i^t (x - p_i) - \rho_i || x - p_i ||^2$$

The corresponding basis functions for each point, with the parameter $\rho_i \geq 0$ defined as:

$$\rho_i = \max_{j \in J_i} \frac{n_i^t (p_j - p_i)}{||p_j - p_i||^2}$$

Where $J_i = \{j : n_i^t (p_j - p_i) > 0\}$

Rotation Invariance

Let $p_i' = Q p_i$ and $n_i' = Q n_i$ with $Q$ an orthogonal matrix, and $f', f_i', \rho_i'$ defined as above. Then, if we look at the $\rho_i'$ estimation:

\begin{align} \rho_i' & = \max_{j \in J_i'} \frac{n_i'^t (p_j' - p_i')}{|| p_j' - p_i'||^2} \\
& = \max_{j \in J_i'} \frac{n_i ^t Q^t Q (p_j - p_i)}{|| Q (p_j - p_i) ||^2} \\
& = \max_{j \in J_i'} \frac{n_i ^t (p_j - p_i)}{|| p_j - p_i ||^2} \end{align}

Furthermore, if we look at the set $J_i'$:

\begin{align} J_i' & = \{j : n_i'^t (p_j' - p_i') > 0\} \\
& = \{j : n_i^t Q^t Q (p_j - p_i) > 0\} \\
& = \{j : n_i^t (p_j - p_i) > 0\} \\
& = J_i \end{align}

Which implies that $\rho_i' = \rho_i$. Finally, $f_i'(Qx)$:

\begin{align} f_i'(Qx) & = n_i'^t (Qx - p_i') - \rho_i' || Qx - p_i' ||^2 \\
& = n_i^t Q^t Q (x - p_i) - \rho_i || Q (x - p_i) ||^2 \\
& = n_i^t (x - p_i) - \rho_i || x - p_i ||^2 \\
& = f_i(x) \end{align}

Thus, when we evaluate on a rotated point $Qx$, we get the guarantee that:

$$f'(Qx) = f(x)$$

Translation Invariance

The proof strategy here is all too similar to the previous case, the difference being that now our points are changed to $p_i' = p_i + c$, and $n_i' = n_i$. In this case, $p_j' - p_i' = p_j - p_i$ because the constant cancels out; this trivially implies that $J_i' = J_i$ and $\rho_i' = \rho_i$. Finally, the new guarantee that we are going to get is that:

$$f_i'(x + c) = f_i(x)$$

Derived out of $c$ canceling out in the sum, and finally:

$$f'(x + c) = f(x)$$

Why does this work?

Imagine $f(x) = 0$, that is, $x$ is in the zero level set of $f$ or, in practical terms, $x$ is a point that belongs to the border of the reconstructed surface. In the rotation case, this implies that

$$f'(Qx) = 0$$

Thus, the rotated point $Qx$ will be a part of the surface. In the case of a translation, we have that

$$f'(x + c) = 0$$

Which implies that the moved point $x + c$ is part of the surface. This of course works both ways, because of the equalities we established above.