# Objective: Segmentation of tangrams from an image and a surface


```python
from torch_pdegraph.operators import GradMinusNorm
import torch
import numpy as np
import matplotlib.image as mpimg
from torch_pdegraph.utilities import *
from tqdm import tqdm
import matplotlib.pyplot as plt
import faiss
```


```python
# load the tangram image
img = mpimg.imread("../data/tangram.png")
dispImg(img)
```


```python
# InitialFronts/initial_distances are seeds on the image (on each tangram  and backgroun)
# One can use any image annotation tool to create hese  seeds I have used Pandore by Greyc
f = plt.figure(figsize=(12,12))
Front = []  
for i in range(1,10,1):
    seed = mpimg.imread(f"../data/dist_{i}.png")
    f.add_subplot(3,3,i)
    plt.imshow(seed, cmap="gray")
    initial_dist = getDist(seed)
    Front.append(initial_dist)
```

# Graph construction

- In order to simply show the effectiveness of PDEs on graph,  I am only creating a simple K-NN based graphs. This may or maynot be the best graph for a given problem at hand.

- One can create graph using whatsoever apt approach or one can even use third-party network datasets and run a PDE on that graph.  PDEs are extensible to any given graph/network at hand as long as that graph has edges and weights( edge_index and edge_attr).

Although torch_cluster comes with a knn-graph method. I found it to be limited and  slow when the node-features have high dimensions. We shall be using facebook's faiss library which is blazingly fast for a KNN-graph construction.


```python
p = toPmat(img.shape).astype(np.float32) # position coords vector
texture = toFmat(img).flatten() # gray-scale values

# Create a KNN graph using the (x,y) pixel-coord of the image and assign weights using differnce
# between the texture at the pixel-coords
res = faiss.StandardGpuResources()
index = faiss.IndexFlatL2(p.shape[1])
gpu_index_flat = faiss.index_cpu_to_gpu(res,0,index)
gpu_index_flat.add(p)
k = 10
_, I = gpu_index_flat.search(p,k+1)


#Graph 
edge_index = np.vstack((I[:,1:].flatten(), np.repeat(I[:,0].flatten(),10)))
edge_attr = np.exp(-((texture[I]-texture[I][:,0,None])**2)/1e-04)[:,1:].flatten()
edge_index = torch.tensor(edge_index, dtype=torch.long).to('cuda:0')
edge_attr = torch.tensor(edge_attr, dtype=torch.float32).to('cuda:0')
edge_attr = edge_attr.view(-1,1)
graph = Graph(edge_index, edge_attr)
```

# Run a manually defined PDE:

\begin{equation}
\mathbf{x}^{n+1}_{i}  = \mathbf{x}^{n}_{i} + \Delta t \ (\mathbf{1} - \|\nabla^{-}_{w}\ \mathbf{x}^{n}_{i}\|_{p})\\
\end{equation}

- $\mathbf{x}_{i}$ is the node feature/signal at the $i^{th}$ node
- $\nabla^{-}_{w}$ is the negative directional gradient on weighted graphs

**Example:** 
```python
from torch_pdegraph.operators import GradMinusNorm
ope = GradMinusNorm.OPE(graph, **hp) # create an operator instance

for i in range(itr): # run the manual pde inside a loop
     fr = fr + dt * (1.0 - ope(fr)) 
```

Though not exactly, this equation is approximates time depenendent **eikonal equation on weighted graphs which can be used to calculate generalized distances on the graphs**. Check the [Thong](https://hal.archives-ouvertes.fr/file/index/docid/404070/filename/VinhThongTA_09_ssvm.pdf)



```python
# Params
dt = 0.05
itr = 5000
hp = dict(p_val = 1)

# Evolved Front
new_front = []

# Instantiate the operator
ope = GradMinusNorm.OPE(graph, **hp)

# Run explict PDE
for fr in tqdm(Front):
    fr = toFmat(fr)
    fr = torch.tensor(fr, dtype=torch.float32).to('cuda:0')
    for i in range(itr):
        fr = fr + dt * (1.0 - ope(fr))
    newfr = fr.to('cpu')
    new_front.append(newfr.numpy())
```

    100%|██████████| 9/9 [00:11<00:00,  1.25s/it]



```python
# visualize the segmented image
dist_img = []
for fr in new_front:
    tmp_img = toImg(fr, img.shape[0:2])
    dist_img.append(tmp_img)
m = np.argmin(dist_img, axis=0)

import matplotlib.pyplot as plt
f = plt.figure(figsize=(8,8))
f.add_subplot(1,2, 1)
plt.imshow(m) # Segmented image
f.add_subplot(1,2, 2)
plt.imshow(img, cmap="gray") # Orignal image
plt.show(block=True)
```

# Reproduce the result for a pointcloud


```python
# Same thing reproduction for a point clould surface
p = np.load("../data/tangram/position_tangram.npz")
t = np.load("../data/tangram/smooth_texture_tangram.npz")
p = p["position"].astype(np.float32)
signal = t["texture"]


# Create a knn graph
res = faiss.StandardGpuResources()
index = faiss.IndexFlatL2(p.shape[1])
gpu_index_flat = faiss.index_cpu_to_gpu(res,0,index)
gpu_index_flat.add(p)
k = 25
_, I = gpu_index_flat.search(p,k+1)
edge_index = np.vstack((I[:,1:].flatten(), np.repeat(I[:,0].flatten(),k)))
edge_attr = np.exp(-np.sum(((signal[I]-signal[I][:,0,None])**2),axis=2)/(1.3e-02)**2)[:,1:].flatten()
edge_index = torch.tensor(edge_index, dtype=torch.long).to('cuda:0')
edge_attr = torch.tensor(edge_attr, dtype=torch.float32).to('cuda:0')
edge_attr = edge_attr.view(-1,1)
graph = Graph(edge_index, edge_attr)


# load the intial_distance/seed pcd
Front = []  # tiff should be preffered than png
for i in range(1,10,1):
    seed = np.load(f"../data/tangram/{i}_texture.npz")
    seed = seed["texture"]
    initial_dist = getDistPcd(seed)
    initial_dist = np.reshape(initial_dist, (len(initial_dist),1))
    Front.append(initial_dist)

# solve the pde
dt = 0.05
itr = 5000
new_front = []
hp = dict(p_val = 1)
ope = GradMinusNorm.OPE(graph, **hp)
for fr in tqdm(Front):
    fr = torch.tensor(fr, dtype=torch.float32).to('cuda:0')
    for i in range(itr):
        fr = fr + dt * (1.0 - ope(fr))
    newfr = fr.to('cpu')
    new_front.append(newfr.numpy())

```

    100%|██████████| 9/9 [00:20<00:00,  2.33s/it]



```python
#visualzie the segmented pcd using matplotlib colormap
m = np.argmin(new_front, axis=0)
m_t = np.uint8(255*(m/np.max(m)))

import matplotlib.cm as cmx
import matplotlib.colors as colors
jet = plt.get_cmap('jet')
cNorm = colors.Normalize(vmin=0, vmax=np.max(m_t))
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
scalarMap.get_clim()
colorVal = scalarMap.to_rgba(m_t)
val_color = np.squeeze(colorVal, axis=1)
val_color = val_color[:,0:3]
```


```python
displayJSur(**dict(position=p, texture=val_color))#Or use displaySur for no jupyter visualization
```


    JVisualizer with 1 geometries



```python
displayJSur(**dict(position=p,texture=signal))#Or use displaySur for no jupyter visualization
```


    JVisualizer with 1 geometries



```python
# Unlike an image, a pcd display cannot be saved in ipython notebook, hence the notebook must be running to 
# to view the pcd. Expect to see the tangrams segmented from the surface like shown in the following image.
from matplotlib import image as mpimg
from matplotlib import pyplot as plt
#res = mpimg.imread("../gallery/tangram_surface.png")
res = mpimg.imread("../gallery/tangram_surf3.png")
f = plt.figure(figsize=(12,12))
plt.imshow(res)
```


```python

```

e5dd5597f6064077777c33e9705a2e41e34968e1
ipynb
Jupyter Notebook
applications/3_segmentation.ipynb
aGIToz/Pytorch_pdegraph
fade6817e437b606c43221a5ca13bdaeec563fff
MIT
2020-08-24T09:04:48.000Z
2022-03-19T03:46:07.000Z
applications/3_segmentation.ipynb
aGIToz/Pytorch_pdegraph
fade6817e437b606c43221a5ca13bdaeec563fff
MIT
applications/3_segmentation.ipynb
aGIToz/Pytorch_pdegraph
fade6817e437b606c43221a5ca13bdaeec563fff
MIT
2020-08-27T15:53:02.000Z
2021-09-13T08:36:27.000Z
Qwen/Qwen-72B
1. YES
2. YES