Astro Processing.Rmd

---
jupyter:
  jupytext:
    formats: ipynb,Rmd
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.13.7
  kernelspec:
    display_name: venv
    language: python
    name: venv
---

```{python}
import os
import sys
```

```{python}
# %matplotlib inline

from matplotlib import rc
import matplotlib.pyplot as plt 
from matplotlib.pyplot import imshow, quiver
```

```{python}
import scipy
from scipy import ndimage as ndi
import networkx as nx
```

```{python}
import numpy as np
from numpy.random import randn
```

```{python}
from tqdm.auto import tqdm
```

```{python}
import napari
```

```{python}
import astromorpho as astro
```

```{python}
data_place = '/home/levtg/astro-morpho/data/'
```

# Image processing


![image-2.png](attachment:image-2.png)

```{python}

```

# Image Data


Линии

```{python}
img_lines = np.zeros((100,100))
img_lines[2::10] = 1
img_lines[3::10] = 1
img_lines = img_lines[20:-20,20:-20]

img_lines = np.clip(img_lines, 0, 1)
crop = (slice(20,-20), slice(20,-20))
img = img_lines
imshow(img, cmap='gray')
```

```{python}
SNR = 5
img_noisy = SNR*img + randn(*img.shape)
imshow(img_noisy, cmap='gray')
```

Сетка

```{python}
img_lines = np.zeros((100,100))
img_lines[2::20] = 1
img_lines[3::20] = 1
img_lines[4::20] = 1

img_lines[..., 2::20] = 1
img_lines[..., 3::20] = 1
img_lines[..., 4::20] = 1

# img_lines = img_lines[20:-20,20:-20]

img_lines = np.clip(img_lines, 0, 1)
crop = (slice(20,-20), slice(20,-20))
img = img_lines[crop]
imshow(img, cmap='gray')
```

```{python}
SNR = 3
img_noisy = SNR*img + randn(*img.shape)
imshow(img_noisy, cmap='gray')
```

Косая Сетка

```{python}
img_lines = np.zeros((100,100))
img_lines[2::20] = 1
img_lines[3::20] = 1
img_lines[4::20] = 1

img_lines = ndi.rotate(img_lines, 24)

img_lines[..., 2::20] = 1
img_lines[..., 3::20] = 1
img_lines[..., 4::20] = 1

img_lines = img_lines[20:-20,20:-20]

img_lines = np.clip(img_lines, 0, 1)
crop = (slice(20,-20), slice(20,-20))
img = img_lines[crop]
imshow(img, cmap='gray')
```

```{python}
SNR = 3
img_noisy = SNR*img + randn(*img.shape)
imshow(img_noisy, cmap='gray')
```

Повернутый крест

```{python}
img_cross = np.zeros((100,100))
img_cross[50:52, 23:78] = 1
img_cross[23:78,50:52] = 1
img_cross = ndi.rotate(img_cross, 24)
img_cross = img_cross[20:-20,20:-20]

img_cross = np.clip(img_cross, 0, 1)
crop = (slice(20,-20), slice(20,-20))
img = img_cross[crop]
imshow(img, cmap='gray', interpolation='nearest')
```

```{python}
SNR = 5
img_noisy = SNR*img + randn(*img.shape)
imshow(img_noisy, cmap='gray',interpolation='nearest')
```

Астроцит (реальный)

```{python}
import ccdb
```

```{python}
def read_image(filename):
    stack, meta = ccdb.read_pic(filename)
    dims = ccdb.get_axes(meta)
    if len(dims):
        zoom = dims[0][0]/dims[-1][0]
    else:
        zoom = 4
    return ndi.zoom(stack.astype(float),(1, zoom,zoom))
```

```{python}
filename = '/home/levtg/astro-morpho/data/3wk-both1-grn-raw.pic'
```

```{python}
stack, meta = ccdb.read_pic(filename)
```

```{python}
stack.shape
```

```{python}
plt.imshow(stack.max(0), cmap='gray')
```

# Image preprocesssing


## CLAHE

```{python}
import cv2

clahe = cv2.createCLAHE(clipLimit =2.0, tileGridSize=(8,8))
```

### 2D

```{python}
img_clahe = clahe.apply(img_noisy.astype('uint8'))
```

```{python}
imshow(img_clahe, cmap='gray')
```

### 3D

```{python}
IMG_SHAPE = stack.shape
stack_flat = stack.reshape((IMG_SHAPE[0], -1))
img_flat = clahe.apply(stack_flat)
img_clahe = img_flat.reshape(IMG_SHAPE)
```

```{python}
plt.imshow(image_clahe.max(0), cmap='gray')
```

 

```{python}
dims = ccdb.get_axes(meta)#[::-1]
if len(dims):
    zoom = dims[0][0]/dims[-1][0]
else:
    zoom = 4
```

```{python}
img_noisy = ndi.zoom(img_clahe.astype(float),(1, zoom,zoom))
```

## Фильтрация изображения

```{python}
from skimage.filters import threshold_li, threshold_minimum
from skimage.morphology import remove_small_objects
```

```{python}
def filter_image(image, filter_func, **kwargs):
    threshold = filter_func(image)
    img_filt = np.where(image > threshold, image, 0)
    binary_clean = remove_small_objects(image >= threshold, 5, connectivity=3)
    return np.where(binary_clean, img_filt, 0)
```

```{python}
# img_clear = filter_image(img_noisy, threshold_li)
img_clear = filter_image(img_noisy, threshold_minimum)
```

### 2D

```{python}
imshow(img_clear, cmap='gray')
```

```{python}
final_image = img_clear
```

### 3D

```{python}
plt.imshow(img_clear.max(0), cmap='gray')
```

```{python}
final_image = img_clear
```

# Soma segmentation


## Erosion

```{python}
from skimage.morphology import erosion, dilation
```

```{python}
w = napari.view_image(final_image)
```

```{python}
to_erose = final_image[:]
```

```{python}
for i in range(5):
    to_erose = erosion(to_erose)
eroded = remove_small_objects(to_erose > 0, 5, connectivity=3)
soma_mask = np.where(dilation(eroded), True, False)
```

```{python}
w = napari.view_image(final_image)
w.add_image(soma_mask, blending='additive', colormap='cyan')
```

```{python}
def soma_extraction(image):
    for i in range(5):
        image = erosion(image)
    eroded = remove_small_objects(image > 0, 5, connectivity=3)
    return np.where(dilation(eroded), True, False)
```

```{python}
data_path = '/home/levtg/astro-morpho/data/'
```

```{python}
# names = !ls $data_path/*.pic
names = [name for name in names if not 'hm' in name or not 'red' in name]
```

```{python}
len(names)/5
```

```{python}
w = napari.view_image(np.array([[0]]))
for imagename in tqdm(names[::10]):
    stack = read_image(imagename)
    image = filter_image(stack, threshold_li)
    eimage = np.where(image>0, np.exp(image), 0)
    w.add_image(image, colormap='magenta', name=os.path.basename(imagename))
    soma = soma_extraction(eimage)
#     soma = filter_image(e_image, lambda x: np.max(x)*0.8)
    w.add_image(soma, blending='additive', colormap='cyan', name=os.path.basename(imagename) + '_soma')
```

# Tree segmentation


## Вектора Гессе - 2d

```{python}
alpha = np.where(final_image > 0, 1, 0)
# alpha = 1
```

```{python}
fig, axes = plt.subplots(2, 2, figsize=(20, 20))

sigma = 1
sato, Vf = astro.hessian.sato2d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)

ax = axes[0][0]
ax.imshow(final_image, cmap='gray', origin='lower')

weights = astro.enh.percentile_rescale(sato)**0.5 * alpha
i = 0
Vfx = Vf[...,i][...,::-1]
V = Vfx[...,0] # row directions (Y)
U = Vfx[...,1] # col directions (X)
h = ax.quiver(U*weights, V*weights, weights, scale=25, cmap='inferno')


sigma = 2
sato, Vf = astro.hessian.sato2d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)

ax = axes[0][1]
ax.imshow(final_image, cmap='gray', origin='lower')

weights = astro.enh.percentile_rescale(sato)**0.5 * alpha
i = 0
Vfx = Vf[...,i][...,::-1]
V = Vfx[...,0] # row directions (Y)
U = Vfx[...,1] # col directions (X)
h = ax.quiver(U*weights, V*weights, weights, scale=25, cmap='inferno')


sigma = 3
sato, Vf = astro.hessian.sato2d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)

ax = axes[1][0]
ax.imshow(final_image, cmap='gray', origin='lower')

weights = astro.enh.percentile_rescale(sato)**0.5 * alpha
i = 0
Vfx = Vf[...,i][...,::-1]
V = Vfx[...,0] # row directions (Y)
U = Vfx[...,1] # col directions (X)
h = ax.quiver(U*weights, V*weights, weights, scale=25, cmap='inferno')


sigma = 4
sato, Vf = astro.hessian.sato2d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)

ax = axes[1][1]
ax.imshow(final_image, cmap='gray', origin='lower')

weights = astro.enh.percentile_rescale(sato)**0.5 * alpha
i = 0
Vfx = Vf[...,i][...,::-1]
V = Vfx[...,0] # row directions (Y)
U = Vfx[...,1] # col directions (X)
h = ax.quiver(U*weights, V*weights, weights, scale=25, cmap='inferno')

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.tight_layout()
```

## Graph creation


В качестве веса ребер будем брать косинусное расстояние между ними


x, y, vectors, weights, img_noisy

```{python}
from tqdm.auto import tqdm
```

```{python}
def weight_of(a, b):
    try:
        cos_dist = scipy.spatial.distance.cosine(a, b)
    except:
        cos_dist = 0
    return np.abs(1 - cos_dist) * np.mean([np.linalg.norm(a), np.linalg.norm(b)])
```

### 2d

```{python}
sigma = 1
```

```{python}
sato, Vf = astro.hessian.sato2d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)
```

```{python}
plt.figure(figsize=(11,9))
imshow(final_image, cmap='gray', origin='lower')

weights = astro.enh.percentile_rescale(sato)**0.5

i = 0
Vfx = Vf[...,i][...,::-1]

V = Vfx[...,0] # row directions (Y)
U = Vfx[...,1] # col directions (X)
plt.setp(plt.gca(), xticks=[], yticks=[])
h = quiver(U*weights, V*weights, weights, scale=25, cmap='inferno')
plt.colorbar(h, ax=plt.gca())
plt.tight_layout()
```

```{python}
lengths = astro.enh.percentile_rescale(sato)**0.5
```

```{python}
vectors = np.stack((U*lengths, V*lengths), axis=2)
```

```{python}
nr, nc = (1, U.shape[0]) if U.ndim == 1 else U.shape
x, y = np.meshgrid(np.arange(nc), np.arange(nr))
```

```{python}
k_h, k_w = 3, 3
```

```{python}
G = nx.Graph()
nodes = {}
```

```{python}
for row in tqdm(range(vectors.shape[0])):
    for col in range(vectors.shape[1]):
        cur = (row, col)
        cur_node = (x[cur], y[cur])
        nodes[cur_node] = cur_node
        G.add_node(cur_node, size=final_image[cur])
        for i in range(k_w):
            for j in range(k_h):
                pos = (np.clip(row + i - k_w//2, 0, U.shape[0] - 1), np.clip(col + j - k_h//2, 0, U.shape[1] - 1)) 
                pos_node = (x[pos], y[pos])
                if G.has_edge(cur_node, pos) or cur_node == pos_node:
                    continue
                weight = weight_of(vectors[cur], vectors[pos])
                if weight < 0.1:
                    continue
                nodes[pos_node] = pos_node
                G.add_edge(cur_node, pos_node, weight=1 - weight)
                
```

Считаем, что точки находятся на одной линии если угол между ними не превышает 45 градусов, т.е. значение весов >= 0.7

```{python}
# # %matplotlib notebook
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in G.edges(data=True)])
# width = np.where(width < 0.3, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(G, nodes, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)
labels = nx.draw_networkx_labels(G, nodes, font_size=6)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
# imshow(clean_img, cmap='gray', origin='lower')
```

Построение пути

```{python}
paths = nx.all_shortest_paths(G, (3, 50), (23, 20), weight='weight')
```

```{python}
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in G.edges(data=True)])
width = np.where(width > 0.7, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(G, nodes, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)
for path in paths:
    path_edges = nx.draw_networkx_nodes(G, nodes, nodelist=path, alpha=0.5)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
# imshow(clean_img, cmap='gray', origin='lower')
```

### 3d 

```{python}
w = napari.view_image(final_image)
```

```{python}
sigma = 2
```

```{python}
sato, Vf = astro.hessian.sato3d(final_image, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)
```

```{python}
Vfx = Vf[...,0][...,::-1]
V = Vfx[..., 0]
U = Vfx[..., 1]
C = Vfx[..., 2]
```

```{python}
lengths = astro.enh.percentile_rescale(sato)**0.5
```

```{python}
vectors = np.stack((U*lengths, V*lengths, C*lengths), axis=3)
```

```{python}
vectors.shape
```

```{python}
import hessian_vectors as hv
```

```{python}
hv.add_hessian_vectors(w, *hv.sato2napari_vectors(sato, Vf))
```

```{python}
nr, nc, nz = (1, U.shape[0]) if U.ndim == 1 else U.shape
x, y, z = np.meshgrid(np.arange(nc), np.arange(nr), np.arange(nz))
```

```{python}
k_h, k_w, k_s = 3, 3, 3
```

```{python}
G = nx.Graph()
nodes = {}
```

```{python}
def create_graph(image, vectors):
    G = nx.Graph()
    nodes = {}
    shape = vectors.shape
    for row in tqdm(range(shape[0])):
        for col in tqdm(range(shape[1])):
            for stk in tqdm(range(shape[2])):
                cur = (row, col, stk)
                cur_node = (x[cur], y[cur], z[cur])
                nodes[cur_node] = cur_node
                G.add_node(cur_node, size=image[cur])
                for i in range(k_w):
                    for j in range(k_h):
                        for k in range(k_s):
                            pos = (np.clip(row + i - k_w//2, 0, shape[0] - 1), 
                                   np.clip(col + j - k_h//2, 0, shape[1] - 1),
                                   np.clip(stk + k - k_s//2, 0, shape[2] - 1)) 
                            pos_node = (x[pos], y[pos], z[pos])
                            if G.has_edge(cur_node, pos) or cur_node == pos_node:
                                continue
                            weight = weight_of(vectors[cur], vectors[pos])
                            if weight < 10e-3:
                                continue
                            nodes[pos_node] = pos_node
                            G.add_edge(cur_node, pos_node, weight=1 - weight)
    return G, nodes
  
```

Считается ооооочень долго ~2.4 часа

```{python}
# G, nodes = create_graph(final_image, vectors)
```

## ОПТИМИЗАЦИЯ


Вместо того, чтобы отдельно обходить каждый пиксель будем смещать весь массив и считать веса через взаимодействие двух массивов.
Всего получитя 18 массивов для объемных изображений и 8 для плоских.


### 2d


lu, uu, ru

ll,........ rr

ld, dd, rd

```{python}
arr = vectors[:]
```

```{python}
lu = np.roll(arr, 1, axis=(0, 1))
lu[0] = 0
lu[:, 0] = 0
uu = np.roll(arr, 1, axis=0)
uu[0] = 0
ru = np.roll(arr, (1, -1), axis=(0, 1))
ru[0] = 0
ru[:, -1] = 0

ll = np.roll(arr, 1, axis=1)
ll[:, 0] = 0
rr = np.roll(arr, -1, axis=1)
rr[:, -1] = 0

ld = np.roll(arr, (-1, 1), axis=(0, 1))
ld[-1] = 0
ld[:, 0] = 0
dd = np.roll(arr, -1, axis=0)
dd[-1] = 0
rd = np.roll(arr, -1, axis=(0, 1))
rd[:, -1] = 0
rd[-1] = 0

```

```{python}
lu.shape
```

Для рассчетов нам понадобится переделать функцию косинусного расстояния:
$1 - \frac{uv}{||u||_2||v||_2}$


Скалярное произведение двух масиивов с векторами

```{python}
mults = np.einsum('...ij,...ij->...i', arr,uu)
```

```{python}
mults.shape
```

Нормы

```{python}
norms = np.linalg.norm(arr, axis=-1)
```

```{python}
norms
```

```{python}
np.linalg.norm(arr[0][0])
```

```{python}
norms_uu = np.linalg.norm(uu, axis=-1)
```

```{python}
cosines = 1 - mults / ((norms * norms_uu))
```

```{python}
arr[1][0]
```

```{python cell_style="center"}
uu[1][0]
```

```{python cell_style="split"}
cosines[1][0]
```

```{python cell_style="split"}
scipy.spatial.distance.cosine(arr[1][0], uu[1][0])
```

Yey!!

```{python}
norms[1][0], norms_uu[1][0]
```

```{python cell_style="split"}
np.mean([norms[1][0], norms_uu[1][0]])
```

```{python cell_style="split"}
np.mean([norms, norms_uu], axis=0)[1][0]
```

```{python}
def weight_of_arr(a, b):
    dprod = np.einsum('...ij,...ij->...i', a, b)
    norm_a = np.linalg.norm(a, axis=-1)
    norm_b = np.linalg.norm(b, axis=-1)
    cos_dist = np.nan_to_num(dprod / ((norm_a * norm_b)), nan=0)
    cos_dist = 1 - cos_dist
    return np.abs(1 - cos_dist) * np.mean([norm_a, norm_b], axis=0)
```

Для определения весов вычитаем полученное значение из единицы тк нужные нам ребра должны быть меньше по весу

```{python}
weights = 1 - weight_of_arr(lu, arr)
```

Индексы

```{python}
i, j = np.indices(arr.shape[:-1])
```

```{python}
idxs = np.stack((i,j), axis=2)
```

Слайсы для выделения нужных данных


![image.png](attachment:image.png)


![image.png](attachment:image.png)


Всего 8 вариантов слайсов

```{python}
lu_crop = (slice(1, None), slice(1, None))
rd_crop = (slice(None, -1), slice(None, -1))

uu_crop = (slice(1, None))
dd_crop = (slice(None, -1))

ru_crop = (slice(1, None), slice(None, -1))
ld_crop = (slice(None, -1), slice(1, None))

rr_crop = (slice(None), slice(None, -1))
ll_crop = (slice(None), slice(1, None))
```

```{python}
shape = idxs.shape
shape
```

```{python}
# idx_lu = [(i, j) for i, j in idxs[rd_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))]
# idx_rd = [(i, j) for i, j in idxs[lu_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))]

# idx_uu = [(i, j) for i, j in idxs[dd_crop].reshape(((shape[0]-1)*(shape[1]), 2))]
# idx_dd = [(i, j) for i, j in idxs[uu_crop].reshape(((shape[0]-1)*(shape[1]), 2))]

# idx_ru = [(i, j) for i, j in idxs[ld_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))]
# idx_ld = [(i, j) for i, j in idxs[ru_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))]

# idx_ll = [(i, j) for i, j in idxs[rr_crop].reshape(((shape[0]-1)*(shape[1]), 2))]
# idx_rr = [(i, j) for i, j in idxs[ll_crop].reshape(((shape[0]-1)*(shape[1]), 2))]

idx_lu = idxs[lu_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))
idx_rd = idxs[rd_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))

idx_uu = idxs[uu_crop].reshape(((shape[0]-1)*(shape[1]), 2))
idx_dd = idxs[dd_crop].reshape(((shape[0]-1)*(shape[1]), 2))

idx_ru = idxs[ru_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))
idx_ld = idxs[ld_crop].reshape(((shape[0]-1)*(shape[1]-1), 2))

idx_ll = idxs[ll_crop].reshape(((shape[0]-1)*(shape[1]), 2))
idx_rr = idxs[rr_crop].reshape(((shape[0]-1)*(shape[1]), 2))
```

Множество ребер определяют фон и не несут необходимой информации. Чтобы избавиться от них воспользуемся методом Ли

```{python}
weights = weight_of_arr(arr[lu_crop], arr[rd_crop])
weight = weights[lu_crop]
```

```{python}
plt.imshow(weights)
plt.colorbar()
```

```{python}
img = weights

li = threshold_li(img)
img_filt = np.where(img > li, img, 0)
binary_clean = remove_small_objects(img >= li, 5, connectivity=3)
weights_clear = np.where(binary_clean, img_filt, 0)
```

```{python}
li
```

```{python}
plt.imshow(weights_clear)
plt.colorbar()
```

Метод отлично выделил фон и сохранил необходимую информацию об изображении


Можно пользоваться только половиной массивов так как значения весов симметричны

```{python}
Gh = nx.Graph()
nodes = {}
```

```{python}
rd_crop = (slice(1, None), slice(1, None))
lu_crop = (slice(None, -1), slice(None, -1))

dd_crop = (slice(1, None))
uu_crop = (slice(None, -1))

ru_crop = (slice(1, None), slice(None, -1))
ld_crop = (slice(None, -1), slice(1, None))

ll_crop = (slice(None), slice(None, -1))
rr_crop = (slice(None), slice(1, None))
```

```{python}
def calc_edges(arr1, arr2, index1, index2):
    weights = weight_of_arr(arr1, arr2)
    weight = weights.ravel()
    li = threshold_li(weight)
    idx1 = [tuple(i) for i in index1.reshape((-1, index1.shape[-1]))[weight>li]]
    idx2 = [tuple(i) for i in index2.reshape((-1, index2.shape[-1]))[weight>li]]
    return zip(idx1, idx2, 1 - weight[weight>li])
```

```{python}
Gh = nx.Graph()
nodes2 = {}
```

```{python}
# %time Gh.add_weighted_edges_from(calc_edges(arr[lu_crop], arr[rd_crop], idxs[lu_crop], idxs[rd_crop]))
# %time Gh.add_weighted_edges_from(calc_edges(arr[uu_crop], arr[dd_crop], idxs[uu_crop], idxs[dd_crop]))
# %time Gh.add_weighted_edges_from(calc_edges(arr[ru_crop], arr[ld_crop], idxs[ru_crop], idxs[ld_crop]))
# %time Gh.add_weighted_edges_from(calc_edges(arr[ll_crop], arr[rr_crop], idxs[ll_crop], idxs[rr_crop]))
```

```{python}
nodes2 = {n:n for n in Gh.nodes()}
```

```{python cell_style="split"}
# Половина массива
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in Gh.edges(data=True)])
# width = np.where(width > 0.7, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(Gh, nodes2, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)
# labels = nx.draw_networkx_labels(G, nodes, font_size=6)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
```

```{python cell_style="split"}
# Обход по картинке
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in G.edges(data=True)])
# width = np.where(width > 0.7, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(G, nodes, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
```

### 3d

```{python}
arr = vectors[:]
```

```{python}
arr.shape
```

В ходе расчетов было выяснено, что нам понадобятся: слайсы


Слайсы

```{python}
iii_crop = (slice(1, None))

rdi_crop = (slice(1, None), slice(1, None), slice(1, None))
lui_crop = (slice(1, None), slice(None, -1), slice(None, -1))

ddi_crop = (slice(1, None), slice(1, None))
uui_crop = (slice(1, None), slice(None, -1))

rui_crop = (slice(1, None), slice(1, None), slice(None, -1))
ldi_crop = (slice(1, None), slice(None, -1), slice(1, None))

lli_crop = (slice(1, None), slice(None), slice(None, -1))
rri_crop = (slice(1, None), slice(None), slice(1, None))


rdd_crop = (slice(None), slice(1, None), slice(1, None))
luu_crop = (slice(None), slice(None, -1), slice(None, -1))

ddd_crop = (slice(None), slice(1, None))
uuu_crop = (slice(None), slice(None, -1))

ruu_crop = (slice(None), slice(1, None), slice(None, -1))
ldd_crop = (slice(None), slice(None, -1), slice(1, None))

lll_crop = (slice(None), slice(None), slice(None, -1))
rrr_crop = (slice(None), slice(None), slice(1, None))

rdo_crop = (slice(None, -1), slice(1, None), slice(1, None))
luo_crop = (slice(None, -1), slice(None, -1), slice(None, -1))

ddo_crop = (slice(None, -1), slice(1, None))
uuo_crop = (slice(None, -1), slice(None, -1))

ruo_crop = (slice(None, -1), slice(1, None), slice(None, -1))
ldo_crop = (slice(None, -1), slice(None, -1), slice(1, None))

llo_crop = (slice(None, -1), slice(None), slice(None, -1))
rro_crop = (slice(None, -1), slice(None), slice(1, None))

ooo_crop = (slice(None, -1))
```

```{python}
# iii, ooo
# rdi, luo
# lui, rdo 
# ddi, uuo 
# uui, ddo 
# rui, ldo 
# ldi, ruo 
# lli, rro 
# rri, llo 
# rdd, luu 
# ddd, uuu 
# ruu, ldd 
# lll, rrr 

crops = [((slice(1, None)), (slice(None, -1))),
         ((slice(1, None), slice(1, None), slice(1, None)), (slice(None, -1), slice(None, -1), slice(None, -1))),
         ((slice(1, None), slice(None, -1), slice(None, -1)), (slice(None, -1), slice(1, None), slice(1, None))),
         ((slice(1, None), slice(1, None)), (slice(None, -1), slice(None, -1))),
         ((slice(1, None), slice(None, -1)), (slice(None, -1), slice(1, None))),
         ((slice(1, None), slice(1, None), slice(None, -1)), (slice(None, -1), slice(None, -1), slice(1, None))),
         ((slice(1, None), slice(None, -1), slice(1, None)), (slice(None, -1), slice(1, None), slice(None, -1))),
         ((slice(1, None), slice(None), slice(None, -1)), (slice(None, -1), slice(None), slice(1, None))),
         ((slice(1, None), slice(None), slice(1, None)), (slice(None, -1), slice(None), slice(None, -1))),
         ((slice(None), slice(1, None), slice(1, None)), (slice(None), slice(None, -1), slice(None, -1))),
         ((slice(None), slice(1, None)), (slice(None), slice(None, -1))),
         ((slice(None), slice(1, None), slice(None, -1)), (slice(None), slice(None, -1), slice(1, None))),
         ((slice(None), slice(None), slice(None, -1)), (slice(None), slice(None), slice(1, None)))]
```

Индексы

```{python}
i, j, k = np.indices(arr.shape[:-1])
```

```{python}
idx = np.stack((i,j,k), axis=3)
```

```{python}
def calc_edges(arr1, arr2, index1, index2):
    weights = weight_of_arr(arr1, arr2)
    weight = weights.ravel()
    li = threshold_li(weight)
#     li = 0.7
    idx1 = [tuple(i) for i in index1.reshape((-1, index1.shape[-1]))[weight>li]]
    idx2 = [tuple(i) for i in index2.reshape((-1, index2.shape[-1]))[weight>li]]
    return zip(idx1, idx2, 1 - weight[weight>li])
```

Расчет весов и создание графа

```{python}
G3 = nx.Graph()
nodes3 = {}
```

```{python}
for crop, acrop in tqdm(crops):
    G3.add_weighted_edges_from(calc_edges(arr[crop], arr[acrop], idx[crop], idx[acrop]))
```

```{python}
nodes3 = {n:n for n in G3.nodes()}
```

```{python cell_style="split"}
len(nodes3) # li=0.027
```

```{python cell_style="split"}
len(nodes3) # li=0.7
```

```{python}
nx.write_weighted_edgelist(G3, '{}_sigma{}_graph_edgelist'.format(os.path.basename(filename), sigma))
```

```{python}
def draw_nodes(pos, nodelist):
    return np.asarray([pos[n] for n in nodelist])
```

```{python}
# %time path = nx.shortest_path(G3, (39, 126, 96), (39, 136, 167), weight='weight')
```

```{python}
xyz = draw_nodes(nodes3, list(nodes3.keys()))
w.add_points(xyz, edge_color='none', size=0.5, face_color='y')
```

```{python}
xyz = draw_nodes(nodes3, path)
```

```{python}
w.add_shapes(xyz, shape_type='path', edge_color='r', edge_width=0.1)
```

# Скелетонизация

```{python}
#2d
G3 = G
nodes3 = {n:n for n in G3.nodes()}
w = napari.view_image(final_image)
```

```{python}
len(G3.nodes())/36000*0.1*60
```

```{python}
nodes = G3.nodes()
```

```{python}
node_coords = np.array(list(nodes.keys()))
```

```{python}
# center = (39, 126, 96)
# center = (29, 109, 101)
# center = (22, 123, 122)
# center = (50, 108, 98)
#2d
center = (23, 20)
```

```{python}
#connected_to_center = [n for n in tqdm(node_coords) if nx.has_path(G3, center, tuple(n))]
```

```{python}
import random as rd
```

```{python}
node_dists = np.linalg.norm(center-node_coords, axis=-1)
```

```{python}
cutoff = np.percentile(node_dists, 90)
cutoff
```

```{python}
selected_nodes = node_coords[node_dists > cutoff]
```

```{python}
selected_nodes.shape
```

```{python}
# paths = []
# failed = []
# for node in tqdm(nodes):
#     try:
#         path = nx.shortest_path(G3, tuple(center), tuple(node), weight='weight')
#         paths.append(path)
#     except Exception as exc:
#         failed.append(exc)
```

```{python}
paths_dict = nx.single_source_dijkstra_path(G3, center)
paths = list(paths_dict.values())
```

```{python}
np.save('{}_sigma{}_paths_all'.format(os.path.basename(filename), sigma), np.array(paths, dtype=object))
```

```{python}
# paths = []
# with open('3wk-both1-grn-raw.pic_sigma1_paths_all.txt') as fd:
#     for line in fd.readlines():
#         path = []
#         for point in line.split(';'):
#             p = tuple(map(int, point.split(',')))
#             path.append(p)
#         paths.append(path)
```

```{python}
def count_points_paths(paths):
    acc = {}
    for path in paths:
        for n in path:
            if n in acc.keys():
                acc[n] += 1
            else:
                acc[n] = 1
    return acc
            
```

```{python}
points = count_points_paths(paths)
```

```{python}
count_1 = []
for path in paths:
    count = 0
    for n in path:
        if points[n] == 1:
            count += 1
    count_1.append(count)
```

```{python}
plt.hist(count_1)
plt.show()
```

```{python}
w.add_points(draw_nodes(nodes3, points), size=1, edge_color='transparent', face_color='b', name='count')
```

Большая часть точек который встречаются один раз - это концевые точки. От них можно избавиться чтобы не засорять данные

```{python}
to_del = []
for point, count in points.items():
    if count == 1:
        to_del.append(point)
```

```{python}
from copy import copy
points_count = copy(points)
```

```{python}
for point in to_del:
    del points_count[point]
```

```{python}
w.add_points(draw_nodes(nodes3, points_count), size=1, edge_color='transparent', face_color='r', name='count')
```

```{python}
len(points), len(points_count)
```

```{python}
p = np.asarray(list(points_count.keys()))
```

```{python}
vals = np.asarray(list(points_count.values()))
```

```{python}
qstack = np.zeros(vectors.shape[:-1])
for loc, val in points_count.items():
    qstack[loc] = np.log(val)
```

```{python}
w.add_image(qstack)
```

```{python}

```

```{python}
plt.figure(figsize=(16,8))
plt.hist(np.log(vals[vals>1]), bins=500)
#plt.xlim(-10, 100)
plt.axvline(np.mean(np.log(vals)), color='cyan', ls='--', label='th.mean')
plt.axvline(threshold_li(np.log(vals)), color='tomato', ls='--', label='th.Li')
plt.legend()
plt.show()
```

```{python}
np.percentile(vals, (5, 25, 50, 75, 95))
```

```{python}
li_vals = threshold_li(vals)
```

```{python}
points_filt = [k for k, v in points_count.items() if v > np.mean(vals)]
```

```{python}
w.add_points(draw_nodes(nodes3, points_filt), size=1, edge_color='transparent', face_color='m', name='filt mean')
```

```{python}
len(points), len(points_count), len(points_filt)
```

```{python cell_style="split"}
np.mean(vals)
```

```{python cell_style="split"}
threshold_li(vals)
```

## Path merging by tree construction

```{python}
import anytree as at
```

```{python}
tree = G3.subgraph(points_count)
leaves = {n:n for n in G3.nodes()}
```

```{python cell_style="split"}
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in tree.edges(data=True)])
# width = np.where(width > 0.7, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(tree, leaves, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)
# labels = nx.draw_networkx_labels(G, nodes, font_size=6)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
```

```{python cell_style="split"}
plt.figure(figsize=(30, 30))
width = np.array([edgedata["weight"] for _, _, edgedata in G3.edges(data=True)])
# width = np.where(width > 0.7, np.zeros(width.shape), width)
edges = nx.draw_networkx_edges(G3, nodes3, width=width*5, edge_cmap = plt.cm.plasma, edge_color=width)
# labels = nx.draw_networkx_labels(G, nodes, font_size=6)

plt.setp(plt.gca(), xticks=[], yticks=[])
plt.colorbar(edges, ax=plt.gca())
plt.tight_layout()
```

```{python}
def draw_edges(pos, edgelist):
    edges = np.asarray([[pos[n1], pos[n2]] for n1, n2 in edgelist])
    return edges
```

```{python}
w.add_shapes(draw_edges(leaves, list(tree.edges())), shape_type='path', edge_color='red', edge_width=0.1)
```

```{python}
# nx.algorithms.threshold.find_threshold_graph(G)
```

## Определение центра

```{python}
import itertools as itt

def percentile_rescale(arr, plow=1, phigh=99.5):

    low, high = np.percentile(arr, (plow, phigh))

    if low == high:

        return np.zeros_like(arr)

    else:

        return np.clip((arr-low)/(high-low), 0, 1)

zstack = final_image[:]
# %time X = np.array(list(itt.product(*map(range, zstack.shape))))
weights_s = percentile_rescale(np.ravel(ndi.gaussian_filter(zstack,5))**2,plow=99.5,phigh=99.99)
center = tuple(map(int, np.sum(X*weights_s[:,None],axis=0)/np.sum(weights_s)))
center
```

## Fat Skeletonization


Проводим скелетонизацию только для толстых отростков

```{python}
from skimage.morphology import flood, skeletonize
```

```{python}
fat_branches = filter_image(img_noisy, threshold_minimum)
```

Используем маску чтобы учитывались только ветки

```{python}
flood_mask = flood(fat_branches > 0, center)
```

```{python}
img = np.where(flood_mask, fat_branches/fat_branches.max(), 0)
skeleton = skeletonize(img)
```

```{python}
fs = napari.view_image(img, colormap='magenta')
```

```{python}
fs.add_image(skeleton, colormap='cyan', blending='additive')
```

# Multi-Scale Vessel Segmentation Using Hessian Matrix Enhancement


Посчитаем гессиан для толстых отростков и проверим работает ли наше предположение, что сигмы центральной линии будут отражать толщину отростка


![image.png](attachment:image.png)
https://content.iospress.com/download/bio-medical-materials-and-engineering/bme1149?id=bio-medical-materials-and-engineering%2Fbme1149

```{python}
def vesselness(img, sigma, alpha= 0.5, beta=0.5, cc=20):
    #H = skf.hessian_matrix(img,sigma)
    #Hl = skf.hessian_matrix_eigvals(H)
    Hl = skridges.compute_hessian_eigenvalues(img, sigma)
    sh = img.shape
    v = np.zeros(sh)
    for r in range(sh[0]):
        for c in range(sh[1]):
            for h in range(sh[2]):
                arr = (Hl[0][r,c,h], Hl[1][r,c,h], Hl[2][r,c,h])

                l1,l2,l3 = sorted(arr, key=abs)
                
                if l2 >= 0 or l3 >=0:
                    v[r,c,h] = 0
                else:
                    Ra = np.abs(l2/l3)
                    Rb = np.abs(l1/np.sqrt(np.abs(l2*l3)))
                    S = np.sqrt(l1**2 + l2*2 + l3**2)
                    v[r,c, h] =  (1 - np.exp(-Ra**2/(2*alpha**2)))*np.exp(-Rb**2/(2*beta**2))*(1 - np.exp(-S**2/(2*c**2)))
    return v

def multiscale_vesselness(img, sigma_start, sigma_stop, nlevels=50, reduce_fn=np.max):
    return reduce_fn([vesselness(img, sigma)*sigma**2 for sigma in np.linspace(sigma_start,sigma_stop, nlevels)],0)

```

```{python}

```