# Introduction

This tutorial will guide you through estimating normal vectors and principal curvatures from 3D point clouds using classic PCA and Jet fitting methods as well as our recent DeepFit method. 


Surface normals and curvatures are a very important properties in shape analysis and are widely used in different fields like computer graphics, computer vision and more. 


### The math
3D point clouds are represented as a matrix of $(x, y, z)$ coordinates of points in 3D space. 
If you are unfamiliar with plane fitting and jets, this section provides a short background. For more in-depth information see the refrences at the bottom. 

The following methods share the following: 
* Input: 3D point cloud  + querry point $q_i$
* Find $q_i$'s k nearest neighbors.  
* Do something fancey :)
* Output: Normal vector at  query point $N_{q_i}$

For estimating the normal at each point we use each point in the point cloud as the query point. 


#### PCA
In this method we estimate the tangent plane to the underlying surface at the query point. This boils down to  solving the eigenvalue and eigenvector decomposition of the covariance matrix created from the points nearest neighbors:

\begin{equation}
C = \frac{1}{k}\sum_{i=1}^{k}(p_i-\hat{p})(p_i-\hat{p})^T
\end{equation}

Here $p_i$ are the neighboring points and $\hat{p}$ is the neighbours centroid. 
The normal vector is the eugenvector associated with the smallest eigenvalue. 

Using this method we cannot directly estimate the principal curvatures (the principal curvatures of a plane are 0).


#### Jet fitting

 An $n$-jet of the height function over a surface is given by:

\begin{equation}
    f(x,y)=J_{\beta,n}(x,y)= \sum_{k=0}^{n}\sum_{j=0}^{k}\beta_{k-j,j}x^{k-j}y^j
\end{equation}

Here $\beta$ is the jet coefficients vector that consists of $N_n=(n+1)(n+2)/2$ terms.

We require that every point satisfy the equation above, yielding the system of linear equations:
\begin{equation}
    M\beta = B
\end{equation}

It is well known that the solution can be expressed in closed-form as: 
\begin{equation}
    \beta = (M^TM)^{-1}M^TB
\end{equation}

Here $M=(1, x_i, y_i, ..., x_i y_i^{n-1}, y_i^n)_{i=1,...,N_p}\in \mathbb{R}^{N_p \times N_n}$ is the Vandermonde matrix and the height function vector $B=(z_1, z_2,...z_{N_p})^T \in \mathbb{R}^N_p$. Both represent the sampled points.


### The code


```python
import sys  
sys.path.insert(0, '../utils')
sys.path.insert(0, '../models')
sys.path.insert(0, '../trained_models')
import DeepFit
import tutorial_utils as tu
import torch
import ipyvolume as ipv
import ipywidgets as widgets
from IPython.display import display
import functools
```

### Generating data
3D point cloud data can be obtained from 3D sensors like LiDAR or RGBD cameras. For this tutorial we will generate a synthetic example. This way we have true ground truth parametric surface as reference for evaluating our results. 


```python
jet_order_data = 3
n_points = 4096
point_cloud_dataset = tu.SyntheticPointCloudDataset(n_points, jet_order_data, points_per_patch=128)
dataloader = torch.utils.data.DataLoader(point_cloud_dataset, batch_size=256, num_workers=8)

points = point_cloud_dataset.points
gt_normals = point_cloud_dataset.gt_normals
```

### Normal estimation


```python
jet_order_fit = 3
for batchind, points in enumerate(dataloader):
    pca_beta, n_est, _ = DeepFit.fit_Wjet(points, torch.ones_like(points[:, 0]), order=1,
                               compute_neighbor_normals=False)
    pca_normals = n_est if batchind==0 else torch.cat([pca_normals, n_est], 0)
    
    jet_beta, n_est, _ = DeepFit.fit_Wjet(points, torch.ones_like(points[:, 0]), order=jet_order_fit,
                               compute_neighbor_normals=False)
    jet_normals = n_est if batchind==0 else torch.cat([jet_normals, n_est], 0)
```

The normals are unoriented so we will now flip them upwards (assuming that the positive z axis points up)


```python
n_sign = torch.sign(torch.sum(pca_normals*
                              torch.tensor([0., 0., 1.]).repeat([n_points, 1]), dim=1)).unsqueeze(-1)
pca_normals = n_sign * pca_normals
pca_normals = pca_normals.detach().cpu().numpy()

n_sign = torch.sign(torch.sum(jet_normals*
                              torch.tensor([0., 0., 1.]).repeat([n_points, 1]), dim=1)).unsqueeze(-1)
jet_normals = n_sign * jet_normals
jet_normals = jet_normals.detach().cpu().numpy()
```

### Visualization


```python
color_n_pca = tu.normal2rgb(pca_normals) #convert normal vectors to RGB
color_n_jet = tu.normal2rgb(jet_normals)
color_n_gt = tu.normal2rgb(gt_normals)

# make some buttons to toggle between the colors 
btn_pca = widgets.Button(description='PCA')
btn_pca.style.button_color='lightgray'
btn_jet = widgets.Button(description='Jet')
btn_jet.style.button_color='lightgray'
btn_solid = widgets.Button(description='solid')
btn_solid.style.button_color='lightgray'
btn_gt = widgets.Button(description='GT')
btn_gt.style.button_color='lightgray'

def update_pc_normal_color(b, scatter_h):
    """    
        this function is linked to the buttons and updates the the point cloud color
    """
    if b.description == 'PCA':
        scatter_h.color = color_n_pca
    elif b.description == 'Jet': 
        scatter_h.color = color_n_pca
    elif b.description == 'solid':
        scatter_h.color = 'red'
    elif b.description == 'GT':
        scatter_h.color = color_n_gt


#plot
fig_h = ipv.figure()
scatter_h = ipv.pylab.scatter(point_cloud_dataset.points[:, 0], 
                              point_cloud_dataset.points[:, 1], 
                              point_cloud_dataset.points[:, 2], size=1, marker="sphere", color='red')
ipv.pylab.xyzlim(-1, 1)
ipv.style.use('minimal')
# ipv.show()

btn_pca.on_click(functools.partial(update_pc_normal_color, scatter_h=scatter_h))
btn_jet.on_click(functools.partial(update_pc_normal_color, scatter_h=scatter_h))
btn_gt.on_click(functools.partial(update_pc_normal_color, scatter_h=scatter_h))
btn_solid.on_click(functools.partial(update_pc_normal_color, scatter_h=scatter_h))

display(widgets.VBox((widgets.HBox((btn_pca, btn_jet, btn_gt, btn_solid)), fig_h)))
```


    VBox(children=(HBox(children=(Button(description='PCA', style=ButtonStyle(button_color='lightgray')), Button(d…


### Loading the data


```python
point_cloud_dataset = tu.SinglePointCloudDataset('./Boxy_smooth100k.xyz', points_per_patch=128)
print(point_cloud_dataset.points)

```

### Fit an n-jet (and compute normal vector)
This may take a while, depending on the number of points


```python
jet_order =3
dataloader = torch.utils.data.DataLoader(point_cloud_dataset, batch_size=256, num_workers=8)

for batchind, data in enumerate(dataloader, 0):
    points = data[0]
    data_trans = data[1]
    beta, n_est, neighbors_n_est = DeepFit.fit_Wjet(points, torch.ones_like(points[:, 0]), order=jet_order,
                               compute_neighbor_normals=False)
    n_est = torch.bmm(n_est.unsqueeze(1), data_trans.transpose(2, 1)).squeeze(dim=1) # cancel out pca
    normals = n_est if batchind==0 else torch.cat([normals, n_est], 0)
```

The normals are unoriented so we will now flip them outwards (assuming that the origin is an internal point and the shape is simple)


```python
n_sign = torch.sign(torch.sum(normals*torch.tensor(point_cloud_dataset.points), dim=1)).unsqueeze(-1)
normals = n_sign * normals
normals = normals.detach().cpu().numpy()
```

### Visualize the point cloud
We plot the point cloud and allow for 3 color overlays:
* Solid - (all points have the same color)
* Normals - We map the normal vectors to the RGB cube and use these values for coloring the point cloud.
* Curvatures - We map the principal curvatures to the XXX colormap to RGB. 


```python
color_n = tu.normal2rgb(normals) #convert normal vectors to RGB

# make some buttons to toggle between the colors 
btn_pc = widgets.Button(description='Solid')
btn_pc.style.button_color='lightgray'
btn_n = widgets.Button(description='normals')
btn_n.style.button_color='lightgray'
btn_c = widgets.Button(description='Curvatures TBD')
btn_c.style.button_color='lightgray'

def update_pc_color_to_solid(b, scatter_h):
    """    
        this function is linked to the buttons and updates the the point cloud color
    """
    scatter_h.color = 'red'

    
def update_pc_color_to_normal(b, scatter_h):
    """    
        this function is linked to the buttons and updates the the point cloud color
    """
    scatter_h.color = color_n

#plot
fig_h = ipv.figure()
scatter_h = ipv.pylab.scatter(point_cloud_dataset.points[:, 0], 
                              point_cloud_dataset.points[:, 1], 
                              point_cloud_dataset.points[:, 2], size=1, marker="sphere", color='red')
ipv.pylab.xyzlim(-1, 1)
ipv.style.use('minimal')
# ipv.show()

btn_pc.on_click(functools.partial(update_pc_color_to_solid, scatter_h=scatter_h))
btn_n.on_click(functools.partial(update_pc_color_to_normal, scatter_h=scatter_h))

display(widgets.VBox((widgets.HBox((btn_pc, btn_n)), fig_h)))
```

### References
* [PCL normal estimation using PCA](http://pointclouds.org/documentation/tutorials/normal_estimation.php)
* [CGAL normal estimation using Jets](https://doc.cgal.org/latest/Jet_fitting_3/index.html#Jet_fitting_3Mathematical)
* [DeepFit paper](https://arxiv.org/pdf/2003.10826.pdf)
* [Jet fitting paper](https://graphics.stanford.edu/courses/cs468-03-fall/Papers/cazals_jets.pdf)


```python

```

6bc907dca0e4ac9b8012548b91416307c5e74194
ipynb
Jupyter Notebook
tutorial/.ipynb_checkpoints/DeepFit_tutorial-checkpoint.ipynb
sitzikbs/DeepF
e6106953ff6d47cb126e1d2ca488a6030bdc9f81
MIT
2020-07-21T09:09:03.000Z
2022-03-28T09:29:22.000Z
tutorial/.ipynb_checkpoints/DeepFit_tutorial-checkpoint.ipynb
sitzikbs/DeepF
e6106953ff6d47cb126e1d2ca488a6030bdc9f81
MIT
2021-01-09T21:44:20.000Z
2021-11-19T08:08:09.000Z
tutorial/.ipynb_checkpoints/DeepFit_tutorial-checkpoint.ipynb
sitzikbs/DeepF
e6106953ff6d47cb126e1d2ca488a6030bdc9f81
MIT
2020-08-09T03:35:24.000Z
2022-03-25T07:00:41.000Z
Qwen/Qwen-72B
1. YES
2. YES