output.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "72f79030",
   "metadata": {
    "papermill": {
     "duration": 0.062546,
     "end_time": "2022-07-21T10:39:19.454002",
     "exception": false,
     "start_time": "2022-07-21T10:39:19.391456",
     "status": "completed"
    },
    "tags": [],
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Next-actions-and-TODOs\" data-toc-modified-id=\"Next-actions-and-TODOs-1\"><strong>Next actions and TODOs</strong></a></span></li><li><span><a href=\"#Параметры-для-запуска\" data-toc-modified-id=\"Параметры-для-запуска-2\">Параметры для запуска</a></span></li><li><span><a href=\"#Считывание-изображения\" data-toc-modified-id=\"Считывание-изображения-3\">Считывание изображения</a></span><ul class=\"toc-item\"><li><ul class=\"toc-item\"><li><span><a href=\"#clahe\" data-toc-modified-id=\"clahe-3.0.1\">clahe</a></span></li></ul></li></ul></li><li><span><a href=\"#Фильтрация-изображения\" data-toc-modified-id=\"Фильтрация-изображения-4\">Фильтрация изображения</a></span></li><li><span><a href=\"#Определение-центра\" data-toc-modified-id=\"Определение-центра-5\">Определение центра</a></span></li><li><span><a href=\"#Выделение-сомы\" data-toc-modified-id=\"Выделение-сомы-6\">Выделение сомы</a></span></li><li><span><a href=\"#Матрица-Гессе\" data-toc-modified-id=\"Матрица-Гессе-7\">Матрица Гессе</a></span></li><li><span><a href=\"#Построение-графа\" data-toc-modified-id=\"Построение-графа-8\">Построение графа</a></span><ul class=\"toc-item\"><li><span><a href=\"#Добавление-точек-сомы-в-граф\" data-toc-modified-id=\"Добавление-точек-сомы-в-граф-8.1\">Добавление точек сомы в граф</a></span></li></ul></li><li><span><a href=\"#Расчет-путей\" data-toc-modified-id=\"Расчет-путей-9\">Расчет путей</a></span></li><li><span><a href=\"#Расчет-встречаемости-и-фильтрация-точек-по-встречаемости\" data-toc-modified-id=\"Расчет-встречаемости-и-фильтрация-точек-по-встречаемости-10\">Расчет встречаемости и фильтрация точек по встречаемости</a></span><ul class=\"toc-item\"><li><span><a href=\"#Ниже-ничего-особо-не-менялось\" data-toc-modified-id=\"Ниже-ничего-особо-не-менялось-10.1\">Ниже ничего особо не менялось</a></span></li></ul></li><li><span><a href=\"#Выделение-тонких-веток\" data-toc-modified-id=\"Выделение-тонких-веток-11\">Выделение тонких веток</a></span></li><li><span><a href=\"#Объединение\" data-toc-modified-id=\"Объединение-12\">Объединение</a></span><ul class=\"toc-item\"><li><span><a href=\"#Выделение-крайних-точек\" data-toc-modified-id=\"Выделение-крайних-точек-12.1\">Выделение крайних точек</a></span></li><li><span><a href=\"#Ближайшие-точки-скелета\" data-toc-modified-id=\"Ближайшие-точки-скелета-12.2\">Ближайшие точки скелета</a></span></li><li><span><a href=\"#Добавление-соединений\" data-toc-modified-id=\"Добавление-соединений-12.3\">Добавление соединений</a></span></li></ul></li><li><span><a href=\"#Многомасштабный-гессиан\" data-toc-modified-id=\"Многомасштабный-гессиан-13\">Многомасштабный гессиан</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6cff623",
   "metadata": {
    "papermill": {
     "duration": 0.06062,
     "end_time": "2022-07-21T10:39:19.573963",
     "exception": false,
     "start_time": "2022-07-21T10:39:19.513343",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# **Next actions and TODOs**\n",
    " - [ ] Test performance on other cells\n",
    " - [ ] Test performace of the approach with more sigma steps (log scale is preferred, i.e. `2.0**np.arange(-1,5,0.5)`)\n",
    " - [ ] Think about a way to regularize vector orientations, using orientations of the neighbours, or at different scales\n",
    " - [ ] Find a best way to skeletonize the qstack-based arrays and masks (as one of the approaches)\n",
    " - [ ] Find a way to \"glue\" together paths, that a close-by and have a similar direction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "481fa7e3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:19.712948Z",
     "iopub.status.busy": "2022-07-21T10:39:19.712784Z",
     "iopub.status.idle": "2022-07-21T10:39:19.717659Z",
     "shell.execute_reply": "2022-07-21T10:39:19.717360Z"
    },
    "papermill": {
     "duration": 0.08576,
     "end_time": "2022-07-21T10:39:19.718889",
     "exception": false,
     "start_time": "2022-07-21T10:39:19.633129",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "675110c5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:19.863975Z",
     "iopub.status.busy": "2022-07-21T10:39:19.863840Z",
     "iopub.status.idle": "2022-07-21T10:39:20.425416Z",
     "shell.execute_reply": "2022-07-21T10:39:20.424998Z"
    },
    "papermill": {
     "duration": 0.640048,
     "end_time": "2022-07-21T10:39:20.426590",
     "exception": false,
     "start_time": "2022-07-21T10:39:19.786542",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fdee9b8d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:20.547286Z",
     "iopub.status.busy": "2022-07-21T10:39:20.547151Z",
     "iopub.status.idle": "2022-07-21T10:39:20.565268Z",
     "shell.execute_reply": "2022-07-21T10:39:20.564847Z"
    },
    "papermill": {
     "duration": 0.0794,
     "end_time": "2022-07-21T10:39:20.566382",
     "exception": false,
     "start_time": "2022-07-21T10:39:20.486982",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import cv2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "500b9fa4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:20.687930Z",
     "iopub.status.busy": "2022-07-21T10:39:20.687785Z",
     "iopub.status.idle": "2022-07-21T10:39:20.787587Z",
     "shell.execute_reply": "2022-07-21T10:39:20.787192Z"
    },
    "papermill": {
     "duration": 0.161417,
     "end_time": "2022-07-21T10:39:20.788756",
     "exception": false,
     "start_time": "2022-07-21T10:39:20.627339",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#import cv2\n",
    "import scipy\n",
    "from scipy import ndimage as ndi\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a8b74af0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:20.909381Z",
     "iopub.status.busy": "2022-07-21T10:39:20.909211Z",
     "iopub.status.idle": "2022-07-21T10:39:20.912281Z",
     "shell.execute_reply": "2022-07-21T10:39:20.911920Z"
    },
    "papermill": {
     "duration": 0.065188,
     "end_time": "2022-07-21T10:39:20.913267",
     "exception": false,
     "start_time": "2022-07-21T10:39:20.848079",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import napari"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a6ba9515",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:21.034389Z",
     "iopub.status.busy": "2022-07-21T10:39:21.034216Z",
     "iopub.status.idle": "2022-07-21T10:39:21.063431Z",
     "shell.execute_reply": "2022-07-21T10:39:21.063057Z"
    },
    "papermill": {
     "duration": 0.091215,
     "end_time": "2022-07-21T10:39:21.064576",
     "exception": false,
     "start_time": "2022-07-21T10:39:20.973361",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from tqdm.auto import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "017174b9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-07-21T10:39:21.185244Z",
     "iopub.status.busy": "2022-07-21T10:39:21.185071Z"
    },
    "lines_to_next_cell": 0,
    "papermill": {
     "duration": 3.874642,
     "end_time": "2022-07-21T10:39:24.998934",
     "exception": false,
     "start_time": "2022-07-21T10:39:21.124292",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Can't load imreg package, affine and homography registrations won't work\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/levtg/astro-morpho/venv/lib64/python3.10/site-packages/imfun/external/tifffile.py:7261: UserWarning: No module named 'imfun.external._tifffile'\n",
      "  Functionality might be degraded or be slow.\n",
      "\n",
      "  warnings.warn(\"%s%s\" % (e, warn))\n",
      "/home/levtg/astro-morpho/venv/lib64/python3.10/site-packages/imfun/external/czifile.py:146: UserWarning: ImportError: No module named '_czifile'. Decoding of JXR and JPEG encoded images will not be available. Czifile.pyx can be obtained at http://www.lfd.uci.edu/~gohlke/\n",
      "  warnings.warn(\n",
      "2022-07-21 13:39:21.923105: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/levtg/astro-morpho/venv/lib/python3.10/site-packages/cv2/../../lib64:/home/levtg/v3d_external/v3d_main/common_lib/lib/\n",
      "2022-07-21 13:39:21.923130: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.\n"
     ]
    }
   ],
   "source": [
    "import ccdb\n",
    "import astromorpho as astro"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "165ac6a7",
   "metadata": {
    "lines_to_next_cell": 0,
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from networx2napari import draw_edges, draw_nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29078ab8",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56c3320e",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def weight_of_arr(a, b):\n",
    "    dprod = np.einsum('...ij,...ij->...i', a, b)\n",
    "    norm_a = np.linalg.norm(a, axis=-1)\n",
    "    norm_b = np.linalg.norm(b, axis=-1)\n",
    "    cos_dist = np.nan_to_num(dprod / ((norm_a * norm_b)), nan=0)\n",
    "    cos_dist = 1 - cos_dist\n",
    "    return np.abs(1 - cos_dist) * np.mean([norm_a, norm_b], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b197c93b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def calc_edges(arr1, arr2, index1, index2):\n",
    "    weights = weight_of_arr(arr1, arr2)\n",
    "    weight = weights.ravel()\n",
    "    li = threshold_li(weight)\n",
    "    idx1 = [tuple(i) for i in index1.reshape((-1, index1.shape[-1]))[weight>li]]\n",
    "    idx2 = [tuple(i) for i in index2.reshape((-1, index2.shape[-1]))[weight>li]]\n",
    "    return zip(idx1, idx2, 1 - weight[weight>li])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96dc6b50",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "def count_points_paths(paths):\n",
    "    acc = defaultdict(int)\n",
    "    for path in paths:\n",
    "        for n in path:\n",
    "            acc[n] += 1\n",
    "    return acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67bbdf04",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage.filters import threshold_li, threshold_minimum\n",
    "from skimage.morphology import remove_small_objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73936a8d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# def remove_small_objects(mask, min_size=25):\n",
    "#     labels, nlab = ndi.label(mask)\n",
    "#     objs = ndi.find_objects(labels)\n",
    "#     out_mask = np.zeros_like(mask)\n",
    "#     for k,o in enumerate(objs):\n",
    "#         cond = labels[o]==(k+1)\n",
    "#         if np.sum(cond) >= min_size:\n",
    "#             out_mask[o][cond] = True\n",
    "#     return out_mask\n",
    "\n",
    "def largest_region(mask):\n",
    "    labels, nlab = ndi.label(mask)\n",
    "    if nlab > 0:\n",
    "        objs = ndi.find_objects(labels)\n",
    "        sizes = [np.sum(labels[o]==k+1) for k,o in enumerate(objs)]\n",
    "        k = np.argmax(sizes)\n",
    "        return labels==k+1\n",
    "    else:\n",
    "        return mask\n",
    "        \n",
    "def crop_image(img, mask=None, margin=0, min_obj_size=0):\n",
    "    if mask is None:\n",
    "        mask = img > 0\n",
    "    if min_obj_size > 0:\n",
    "        mask = remove_small_objects(mask, min_obj_size)\n",
    "    if margin > 0:\n",
    "        mask = ndi.binary_dilation(mask, iterations=margin)\n",
    "    objs = ndi.find_objects(mask)\n",
    "    min_bnds = np.min([[sl.start for sl in o] for o in objs],0)\n",
    "    max_bnds = np.max([[sl.stop for sl in o] for o in objs],0)\n",
    "    crop = tuple(slice(mn,mx) for mn,mx in zip(min_bnds, max_bnds))\n",
    "    return img[crop]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0037ed4a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.rc('figure', dpi=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4174fe7f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Параметры для запуска"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7751d8f3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "add080e4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "#filename = '3wk-both1-grn-raw.pic'\n",
    "data_dir = '/home/levtg/astro-morpho/data/'\n",
    "# data_dir = '/home/brazhe/yadisk/data-shared-comfi/3D-astrocyte-images/selected-for-complexity/'\n",
    "filename = '3wk-both1-grn-raw.pic'\n",
    "\n",
    "verbose = True\n",
    "sigma = 2\n",
    "\n",
    "# Set false to start from console\n",
    "HANDY = True\n",
    "\n",
    "# Set true to save output\n",
    "OUT = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50cdeb31",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": [
     "injected-parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "filename = \"/home/levtg/astro-morpho/data/1wk-both1-red-raw.pic\"\n",
    "sigma = 1\n",
    "HANDY = False\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47f5e45d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "filename = Path(data_dir).joinpath(filename)\n",
    "filename"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f53efa3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Считывание изображения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1f5fb8a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# if HANDY:\n",
    "#     filename = '/home/levtg/astro-morpho/data/3wk-ly10-raw.pic'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db0ebb4b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "stack, meta = ccdb.read_pic(filename)\n",
    "dims = ccdb.get_axes(meta)\n",
    "dims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "913030c0",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if len(dims):\n",
    "    zoom = (dims[-1][0]/dims[0][0])\n",
    "else:\n",
    "    zoom = 4\n",
    "    \n",
    "print(zoom)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84937f04",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### clahe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52c79053",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "clahe = cv2.createCLAHE(clipLimit =2.0, tileGridSize=(8,8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "958fb8a3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "stack_shape = stack.shape\n",
    "#stack_roll = np.reshape(stack, (stack_shape[0],-1))\n",
    "#print(\"Rolled shape: \", stack_roll.shape)\n",
    "#%time img = clahe.apply(stack_roll).reshape(stack.shape)\n",
    "\n",
    "img = np.zeros(stack.shape, np.float32)\n",
    "for k,plane in enumerate(stack):\n",
    "    img[k] = clahe.apply(plane)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dbfb010",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    wi = napari.view_image(stack, ndisplay=3, scale=(zoom, 1,1), name='raw', colormap='cyan')\n",
    "    wi.add_image(img, scale=(zoom,1,1), name='CLAHE',colormap='magenta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7de3778",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "threshold_li(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd54ac48",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "max_proj = img.max(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "988ee624",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "domain_mask = ndi.binary_dilation(largest_region(remove_small_objects(max_proj > 0.5*threshold_li(max_proj))), iterations=3)\n",
    "domain_mask = ndi.binary_closing(domain_mask,iterations=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcd9bb81",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.imshow(max_proj, cmap='gray')\n",
    "plt.contour(domain_mask, colors=['r'], levels=[0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45d9256",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "img_cropped = np.array([crop_image(plane,domain_mask, margin=10) for plane in img])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07b03c7c",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "max_proj_1 = img_cropped.max(1)\n",
    "domain_mask_1 = ndi.binary_dilation(largest_region(remove_small_objects(max_proj_1 > 0.5*threshold_li(max_proj_1))), iterations=3)\n",
    "domain_mask_1 = ndi.binary_closing(domain_mask_1,iterations=3)\n",
    "plt.imshow(max_proj_1, cmap='gray')\n",
    "plt.contour(domain_mask_1, colors=['r'], levels=[0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dca9693",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "img_cropped = np.array([crop_image(img_cropped[:,i],domain_mask_1, margin=10) for i in range(img_cropped.shape[1])]).swapaxes(0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0229bf10",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    napari.view_image(img_cropped)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d073d2a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Важный вопрос, как сделать одинаковым масштаб по осям z и xy. Можно downsample XY, можно upsample (by interpolation) Z. Можно комбинировать. В идеале, наверное, XY не трогать, а сделть upsample по Z. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26bd2348",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "downscale = 2\n",
    "%time img_noisy = ndi.zoom(img_cropped.astype(np.float32), (zoom/downscale, 1/downscale, 1/downscale), order=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b1860bc",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.imshow(img_noisy.max(0), cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f89ad98",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "img.shape, img_noisy.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e069682",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Фильтрация изображения"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af844d47",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def filter_image(image, filter_func):\n",
    "    threshold = filter_func(image)\n",
    "    #img_filt = np.where(image > threshold, image, 0)\n",
    "    binary_clean = remove_small_objects(image >= threshold, 5, connectivity=3)\n",
    "    return np.where(binary_clean, image, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1dbc861",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "img_clear = filter_image(img_noisy, threshold_li)\n",
    "if HANDY:\n",
    "    pass\n",
    "#     img_clear = filter_image(img_noisy, threshold_minimum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61cd8785",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "final_image = img_clear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52cab0bf",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "final_image.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0909870c",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# НЕ ЗАКРЫВАТЬ!!!\n",
    "# allW = napari.view_image(final_image, colormap='gray', ndisplay=3, opacity=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78c4ae07",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Определение центра"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "050cfff4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import itertools as itt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc52c385",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def percentile_rescale(arr, plow=1, phigh=99.5):\n",
    "    low, high = np.percentile(arr, (plow, phigh))\n",
    "    if low == high:\n",
    "        return np.zeros_like(arr)\n",
    "    else:\n",
    "        return np.clip((arr-low)/(high-low), 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8acffa2",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def flat_indices(shape):\n",
    "    idx = np.indices(shape)\n",
    "    return np.hstack([np.ravel(x_)[:,None] for x_ in idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a133719",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "X1a = flat_indices(final_image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddf015fc",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%time weights_s = percentile_rescale(np.ravel(ndi.gaussian_filter(final_image,5))**2,plow=99.5,phigh=99.99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ecbb467",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "center = tuple(map(int, np.sum(X1a*weights_s[:,None],axis=0)/np.sum(weights_s)))\n",
    "center"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57b6e33d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Выделение сомы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51ab9b5d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage.morphology import dilation, skeletonize, flood"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e77930d3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from astromorpho import morpho"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa9a9f5",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "**Альтернативный подход к сегментации сомы**\n",
    "1. Работаем со сглаженным стеком\n",
    "2. делаем первичную маску как flood из центра с толерантностью в 10% разницы между максимальным и минимальным значениями в стеке\n",
    "3. Разрастаем (аналог flood) первичную маску в несколько итераций"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0a833c0",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#soma_mask = largest_region(np.where(dilation(eroded), True, False))\n",
    "#soma_mask = largest_region(final_image >= np.percentile(final_image, 99))\n",
    "\n",
    "smooth_stack = ndi.gaussian_filter(final_image, 3)\n",
    "tol = (smooth_stack.max() - smooth_stack[final_image>0].min())/10\n",
    "\n",
    "print('tol:',tol)\n",
    "%time soma_seed_mask = flood(smooth_stack, center, tolerance=tol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "982459fd",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%time soma_mask = morpho.expand_mask(soma_seed_mask, smooth_stack, iterations = 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00a9d2e7",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    w = napari.view_image(final_image, ndisplay=3, opacity=0.5)\n",
    "    w.add_image(soma_seed_mask, blending='additive', colormap='cyan')\n",
    "    w.add_image(soma_mask, blending='additive', colormap='magenta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a32784ab",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%time soma_all  = morpho.mask2points(soma_mask)\n",
    "%time soma_shell = morpho.inside_boundary_pixels(soma_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1adb16ae",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Матрица Гессе"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7096012",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sigmas = 2.0**np.arange(-1, 4, 1)\n",
    "sigmas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8de6d61a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstacks = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6deb8db4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sato_coll = {}\n",
    "Vf_coll = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1ff1321",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sigma in tqdm(sigmas):\n",
    "    #astro.morpho.sato3d is newer and uses tensorflow (if it's installed)\n",
    "    #optimally, the two variants of sato3d should be merged\n",
    "    sato, Vf = astro.morpho.sato3d(final_image, sigma, hessian_variant='gradient_of_smoothed', do_brightness_correction=False, return_vectors=True)\n",
    "    sato_coll[sigma] = (sato*sigma**2)*(final_image > 0)\n",
    "    #sato_coll[sigma] = final_image*sato*sigma**2\n",
    "    #sato_coll[sigma] = sato*(final_image>0)\n",
    "    Vf_coll[sigma] = Vf[...,0][...,::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ed54f68",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sato_coll.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9b89e19",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "lengths_coll = {sigma: astro.enh.percentile_rescale(sato)**0.5 for sigma, sato in sato_coll.items()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c999466a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "vectors_coll = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50820f4b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sigma in Vf_coll:\n",
    "    Vfx = Vf_coll[sigma]\n",
    "    V = Vfx[..., 0]\n",
    "    U = Vfx[..., 1]\n",
    "    C = Vfx[..., 2]\n",
    "    lengths = lengths_coll[sigma]\n",
    "    vectors_coll[sigma] = np.stack((U*lengths, V*lengths, C*lengths), axis=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5906ce81",
   "metadata": {
    "lines_to_next_cell": 0,
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sato_best = np.argmax([sato_coll[sigma] for sigma in sigmas], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c93a3b8b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "015340d3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from ucats import masks as umasks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd108316",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "ksigma = len(sigmas)-1\n",
    "# largest_sigma = sigmas[ksigma]\n",
    "# largest_sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bc077f2",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# largest_sigma_seed = largest_region(umasks.select_overlapping(sato_best == ksigma, soma_mask))\n",
    "sato = sato_coll[sigmas[1]]#*(final_image)\n",
    "#tol = (sato.max()-sato.min())/5\n",
    "#print('tol:', tol)\n",
    "#largest_sigma_seed2 =  largest_region(sato >= 0.5*np.mean(sato[soma_mask]))\n",
    "threshold = threshold_li(sato[sato>0])\n",
    "mask = remove_small_objects(sato>threshold, int(sigma*64))\n",
    "#largest_sigma_seed = flood(sato, center, tolerance=tol)\n",
    "# np.any(largest_sigma_seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f674937",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "masks = {}\n",
    "for sigma in tqdm(sigmas):\n",
    "    sato = sato_coll[sigma]\n",
    "    threshold = threshold_li(sato[sato>0])*sigma**0.5\n",
    "    print(sigma, threshold)\n",
    "    masks[sigma] = remove_small_objects(sato > threshold, min_size=int(sigma*64))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a94277c",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "masks[sigmas[-1]] = umasks.select_overlapping(masks[sigmas[-1]], soma_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aae809d4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for k in range(len(sigmas)-2,-1,-1):\n",
    "    sigma = sigmas[k]\n",
    "    masks[sigma] = umasks.select_overlapping(masks[sigma], ndi.binary_dilation(masks[sigmas[k+1]], iterations=5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d1fcbba",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    w = napari.view_image(final_image, )\n",
    "    for sigma in masks:\n",
    "        w.add_image(masks[sigma], blending='additive', name=f'σ={sigma:02f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66276c18",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "vectors_best = np.zeros(vectors_coll[sigmas[0]].shape)\n",
    "for k in range(len(sigmas)-1,-1,-1):\n",
    "    sigma = sigmas[k]\n",
    "    mask = masks[sigma]\n",
    "    if k < len(sigmas)-1:\n",
    "        mask = mask & (mask ^ masks[sigmas[k+1]])\n",
    "    vectors_best[mask] = vectors_coll[sigma][mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2df4d8b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from hessian_vectors import add_hessian_vectors, sato2napari_vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9561feb4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def vectors2napari(vectors, index=1):\n",
    "    nr, nc, nd = vectors.shape[:-1]\n",
    "    indexgrid = np.meshgrid(np.arange(nc), np.arange(nr), np.arange(nd))\n",
    "    \n",
    "#     length = np.ravel(np.linalg.norm(vectors, axis=-1))\n",
    "\n",
    "    x, y, z = [np.ravel(a)[::index] for a in indexgrid]\n",
    "    x1, y1, z1 = np.ravel(vectors[...,0])[::index], np.ravel(vectors[...,1])[::index], np.ravel(vectors[...,2])[::index]\n",
    "\n",
    "    \n",
    "    \n",
    "    vectors = np.zeros((*x1.shape, 2, 3))\n",
    "    vectors[..., 0, 0] = y\n",
    "    vectors[..., 0, 1] = x\n",
    "    vectors[..., 1, 0] = y1\n",
    "    vectors[..., 1, 1] = x1\n",
    "\n",
    "    vectors[..., 0, 2] = z\n",
    "    vectors[..., 1, 2] = z1\n",
    "    \n",
    "    return vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1259d7d5",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# vb = napari.view_image(final_image)\n",
    "# add_hessian_vectors(vb, vectors2napari(vectors_best, index=50), np.ravel(sato_best)[::50])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83e8e6fe",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Построение графа"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5614b8f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def prep_crops():\n",
    "    \"makes list of crops for edges\"\n",
    "    num2slice = {1: (slice(1,None), slice(None,-1)), \n",
    "                 0: (slice(None), slice(None)), \n",
    "                -1: (slice(None,-1), slice(1,None))}\n",
    "    shifts = list(itt.product(*[(-1,0,1)]*3))\n",
    "    # we only need one half of that\n",
    "    cut = int(np.ceil(len(shifts)/2))\n",
    "    crops_new = [list(zip(*[num2slice[n] for n in tuple])) for tuple in shifts[cut:]]\n",
    "    return crops_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6d87e43",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "crops_new = prep_crops()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8d5ea31",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "graph_coll = {sigma:nx.Graph() for sigma in sigmas}\n",
    "nodes_coll = {sigma:{} for sigma in sigmas}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b969f80",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "graph_coll['best'] = nx.Graph()\n",
    "nodes_coll['best'] = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88d188b1",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "i, j, k = np.indices(final_image.shape)\n",
    "idx = np.stack((i,j,k), axis=3)\n",
    "idx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6ac95a0",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for sigma in sigmas:\n",
    "    vectors = vectors_coll[sigma]#*(largest_sigma_seed[...,None])\n",
    "    vectors = vectors*masks[sigma][...,None]\n",
    "    #if sigma == largest_sigma:\n",
    "    #    vectors = vectors*largest_sigma_seed[...,None]\n",
    "    for crop, acrop in tqdm(crops_new):\n",
    "         graph_coll[sigma].add_weighted_edges_from(calc_edges(vectors[crop], vectors[acrop], idx[crop], idx[acrop]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b92dd4",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "key='best'\n",
    "vectors = vectors_best\n",
    "graph_coll[key] = nx.Graph()\n",
    "for crop, acrop in tqdm(crops_new):\n",
    "         graph_coll[key].add_weighted_edges_from(calc_edges(vectors[crop], vectors[acrop], idx[crop], idx[acrop]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a01544ed",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Добавление точек сомы в граф"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bd3fda8",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_mask_vals(idxs, mask):\n",
    "    idx_mask = mask[idxs[:,0], idxs[:,1], idxs[:,2]]\n",
    "    return idxs[idx_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51c141cd",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_edges(mask, index1, index2, weight):\n",
    "    idx1 = [tuple(i) for i in get_mask_vals(index1.reshape((-1, index1.shape[-1])), mask)]\n",
    "    idx2 = [tuple(i) for i in get_mask_vals(index2.reshape((-1, index2.shape[-1])), mask)]\n",
    "    return zip(idx1, idx2, np.full(len(idx1), weight))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a85469",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "Gsoma = nx.Graph()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18b8bfe9",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for crop, acrop in tqdm(crops_new):\n",
    "    Gsoma.add_weighted_edges_from(get_edges(soma_mask, idx[crop], idx[acrop], 0.7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1efdb175",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# nx.set_edge_attributes(G, 0.7, name='weight')\n",
    "# nx.set_edge_attributes(G, 'soma', name='area')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "512a6367",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "graph_coll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cac0347",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "\n",
    "#G = graph_coll[largest_sigma]\n",
    "#G = graph_coll[current_sigma]\n",
    "\n",
    "for key in graph_coll:\n",
    "    G = graph_coll[key]\n",
    "    for p1, p2, weight in Gsoma.edges(data=True):\n",
    "        try:\n",
    "            old_weight = G.get_edge_data(p1, p2)['weight']\n",
    "        except Exception as exc:\n",
    "            old_weight = 1\n",
    "        G.add_edge(p1, p2, weight=min(weight['weight'], old_weight))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cdceede",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "nodes_coll = {key:{n:n for n in G.nodes()} for key, G in graph_coll.items()} # just a copy of G3 nodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2833c019",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Расчет путей"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1cf6433",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "\n",
    "# fix negative weights\n",
    "\n",
    "for G in graph_coll.values():\n",
    "#     edges = [(u,v,d) for u,v,d in G.edges(data=True) if d['weight']<0]\n",
    "    for p1,p2, weight in tqdm(G.edges(data=True)):\n",
    "        if weight['weight'] < 0:\n",
    "            print(p1,p2,weight)\n",
    "            G.add_edge(p1,p2, weight=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7d88993",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time \n",
    "\n",
    "# fix negative weights\n",
    "\n",
    "for G in tqdm(graph_coll.values()):\n",
    "    edges = {(u,v): 0 for u,v,d in G.edges(data=True) if d['weight']<0}\n",
    "    nx.set_edge_attributes(G, edges, 'weight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a52cf12",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "len(soma_shell)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45a1bb78",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Расчет встречаемости и фильтрация точек по встречаемости"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "040be23d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from copy import copy\n",
    "\n",
    "def make_qstack(G, targets, min_count=1):\n",
    "    paths_dict = nx.multi_source_dijkstra_path(G, targets, )\n",
    "    paths = list(paths_dict.values())\n",
    "    points = count_points_paths(paths)\n",
    "    #to_del = [point for point,count in points.items() if count <=1]\n",
    "    #points_count = copy(points)\n",
    "    #for point in to_del:\n",
    "    #    del points_count[point]\n",
    "    qstack = np.zeros(vectors.shape[:-1])  #Это встречаемость точек в путях\n",
    "    for p, val in points.items():\n",
    "        if val >= min_count:\n",
    "            qstack[p] = np.log(val)\n",
    "    return qstack, paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe05e524",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#qstacks  = {sigma:make_qstack(graph_coll[sigma], soma_shell)[0] for sigma in tqdm(graph_coll)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a0020d6",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstack_masks = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e4f6069",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sigmas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86105395",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "for ksigma in range(len(sigmas)-1,-1,-1):\n",
    "    #alternative: use shell of a full mask as target\n",
    "    sigma = sigmas[ksigma]\n",
    "    print(ksigma, sigma)\n",
    "    G = graph_coll[sigma]\n",
    "    \n",
    "    if ksigma == len(sigmas)-1:\n",
    "        targets = soma_shell\n",
    "    else:\n",
    "        higher_sigma = sigmas[ksigma+1]        \n",
    "        #targets = astro.morpho.boundary_pixels(masks[higher_sigma])\n",
    "        targets = astro.morpho.mask2points(qstack_masks[higher_sigma])\n",
    "        targets = [t for t in targets if G.has_node(t)]\n",
    "        if not len(targets):\n",
    "            print('No points from targets in the graph')\n",
    "            continue\n",
    "    qstack = make_qstack(G, targets)[0]\n",
    "    qstack_mask = qstack > threshold_li(qstack[qstack>0])\n",
    "    if not np.any(qstack_mask):\n",
    "        print('empty qstack mask!')\n",
    "        if ksigma < len(sigmas)-1:\n",
    "            qstack_mask = qstack_masks[higher_sigma]\n",
    "    qstacks[sigma] = qstack\n",
    "    qstack_masks[sigma] = qstack_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38f8bf56",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "isinstance(2.0, (int, float))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4279f8a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Можно предложить как минимум, два варианта объединения масштабов:\n",
    " 1. [ ] \"Best\" -- это где вектора в каждом вокселе взяты из соответствующих масок для разных масштабов, потом все это сведено в один граф, и во всем графе\n",
    "         ищется путь до поверхности сомы. **NOTE:** по идее, маски должны быть \"исключительными\", то есть каждая область может принадлежать только одной сигме.\n",
    " 2. [ ] \"Combined\" -- скелет и пути задаются итеративно от больших масштабов к маленьким, то есть используется свой граф для каждого масштаба и пути ищутся в дополнение к уже найденым. \n",
    "       Кстати, можно сделать лучше (предположительно), если вектора из qstack_mask старшего масштаба добавлять к графу меньшего масштаба и опять искать пути до сомы. Тогда будут дополнительно \n",
    "       \"тренироваться\" пути вдоль больших веток. \n",
    "       Потом можно брать просто сумму qstacks для разных масштабов, маску можно брать как объединение всех масок на разных уровнях или снова как надпороговые пиксели. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86b59732",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstacks['best'] = make_qstack(graph_coll['best'], soma_shell)[0]\n",
    "qstack_masks['best'] = qstacks['best'] > threshold_li(qstacks['best'][qstacks['best']>0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01387456",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstacks['combined'] = sum(qstacks[sigma] for sigma in sigmas)\n",
    "qstack_masks['combined'] = sum(qstack_masks[sigma] for sigma in sigmas) > 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "524b2989",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#sigma = sigmas[ksigma]\n",
    "#sigma = largest_sigma\n",
    "\n",
    "if verbose:\n",
    "    #w.add_image(qstack)\n",
    "    w = napari.view_image(final_image, ndisplay=3, opacity=0.5)\n",
    "    w.add_image(soma_mask, name='soma', interpolation='nearest', blending='additive',colormap='cyan')\n",
    "    #w.add_image(sato_coll[current_sigma], name=f'sato@{current_sigma:0.2f}', colormap='inferno', blending='additive')\n",
    "    #w.add_image(largest_sigma_seed.astype(np.float32), name=f'largest_sigma_seed', interpolation='nearest', colormap='blue', blending='additive')\n",
    "    for sigma in qstacks:\n",
    "        name = sigma if isinstance(sigma, str) else f'qstack@{sigma:0.2f}'\n",
    "        qstack = qstacks[sigma]\n",
    "        data_range = np.percentile(qstack[qstack>0], (1,99))\n",
    "        w.add_image(qstack, name=name, interpolation='nearest', blending='additive', contrast_limits=data_range, visible=False)\n",
    "        #threshold = threshold_li(qstack[qstack>0])\n",
    "        w.add_image(qstack_masks[sigma], name=name+':qmask', blending='additive', colormap='red')\n",
    "        if isinstance(sigma, (int, float)):\n",
    "            w.add_image(masks[sigma], name=f'mask@{sigma:0.2f}', blending='additive', colormap='green',visible=False)\n",
    "    #w.add_image(qstacks[sigma], name=f'qstack@{sigma:0.2f}', interpolation='nearest', colormap='magenta', blending='additive',contrast_limits=(0,12))\n",
    "    #w.add_image(sato_best == len(sigmas)-1, name='kbest', colormap='red',blending='additive')\n",
    "    #w.add_image(qstacks[sigma]*largest_sigma_seed, name=f'x@{sigma:0.2f}', interpolation='nearest', colormap='green', blending='additive', contrast_limits=(0,12))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6958183",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "combined_mask = sum(qstack_masks[sigma] for sigma in qstack_masks) > 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28257971",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "combined_mask = remove_small_objects(combined_mask, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cec38e5",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "combined_skeleton = skeletonize(combined_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "926e54bd",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "w = napari.view_image(final_image, ndisplay=3, opacity=0.5)\n",
    "w.add_image(combined_mask, colormap='green', blending='additive')\n",
    "w.add_image(combined_skeleton, colormap='red', blending='additive')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd2ec66",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#qstacks[largest_sigma]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3335bee",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstacks[sigmas[-2]] = make_qstack(graph_coll[sigmas[-2]], astro.morpho.mask2points(qstacks[sigmas[-1]]>3.5))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab925f38",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Ниже ничего особо не менялось"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff5e5ddd",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if OUT:\n",
    "    np.save(f'{os.path.basename(filename)}-qstacks.npy',qstacks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7523aba",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "qstack_mask = qstack > 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07018947",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#napari.view_image(skeletonize(qstack_mask))\n",
    "# = napari.view_image(final_image, ndisplay=3)\n",
    "#w.add_image(qstack_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7e4d37",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "p = np.asarray(list(points_count.keys()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d7321fb",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "vals = np.asarray(list(points_count.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9cdd0e1",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,8))\n",
    "plt.hist(np.log(vals[vals>1]), bins=500)\n",
    "#plt.xlim(-10, 100)\n",
    "plt.axvline(np.mean(np.log(vals)), color='cyan', ls='--', label='th.mean')\n",
    "plt.axvline(threshold_li(np.log(vals)), color='tomato', ls='--', label='th.Li')\n",
    "plt.axvline(np.log(np.mean(vals)), color='b', ls='--', label='th.log_mean')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da7e2d10",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,8))\n",
    "plt.hist(vals[vals>1], bins=500)\n",
    "#plt.xlim(-10, 100)\n",
    "plt.axvline(np.mean(vals), color='b', ls='--', label='th.mean')\n",
    "plt.axvline(threshold_li(vals), color='tomato', ls='--', label='th.Li')\n",
    "plt.axvline(np.exp(np.mean(np.log(vals))), color='cyan', ls='--', label='th.log_mean')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9807cadd",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "points_mean = [k for k, v in points_count.items() if v > np.mean(vals)]\n",
    "thr_li = threshold_li(vals)\n",
    "points_li = [k for k, v in points_count.items() if v > thr_li]\n",
    "# points_log_mean = [k for k, v in points_count.items() if np.log(v) > np.mean(np.log(vals))]\n",
    "# points_log_li = [k for k, v in points_count.items() if np.log(v) > threshold_li(np.log(vals))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cfebb6a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    w.add_points(draw_nodes(nodes3, points_mean), size=1, edge_color='transparent', face_color='c', name='filt mean')\n",
    "    w.add_points(draw_nodes(nodes3, points_li), size=1, edge_color='transparent', face_color='g', name='filt li')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70a2a60d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Выделение тонких веток"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac8d738a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "Gsub = G3.subgraph(points_mean)\n",
    "nodes_sub = {n:n for n in Gsub.nodes()}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da77b966",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "props = {'weight': 1 - np.array([edgedata[\"weight\"] for _, _, edgedata in Gsub.edges(data=True)])}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8959575b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "w.add_shapes(draw_edges(nodes_sub, list(Gsub.edges())), shape_type='path', edge_color='weight', edge_width=0.1, edge_colormap='inferno', properties=props)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5487d71f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "c = np.array(Gsub.nodes())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e6640b9",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "graph_mask = np.full(final_image.shape, False)\n",
    "graph_mask[c[:,0], c[:,1], c[:,2]] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de03b999",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "th_branch_mask = np.full(final_image.shape, False)\n",
    "th_branch_mask[c[:,0], c[:,1], c[:,2]] = True\n",
    "th_branch_mask[flood_mask] = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "909e0072",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "th_branch_mask = remove_small_objects(th_branch_mask, 5, connectivity=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "195c9b25",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Объединение"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f033ca",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "branch_mask = np.full(final_image.shape, False)\n",
    "branch_mask[th_branch_mask] = True\n",
    "branch_mask[w_branch_mask] = True\n",
    "branch_mask[soma_mask] = True\n",
    "\n",
    "if verbose:\n",
    "    bm = napari.view_image(final_image, colormap='magenta', ndisplay=3)\n",
    "    bm.add_image(branch_mask, blending='additive', colormap='cyan')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68e32251",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "file:///home/levtg/astro-morpho/data/pictures_nb/join%20plan.jpg![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4226cd2a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Выделение крайних точек"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa8f738c",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage.morphology import cube\n",
    "edge = dilation(flood_mask, cube(3)) ^ flood_mask\n",
    "points = edge & th_branch_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97ce3114",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    bm.add_image(points, blending='additive', colormap='red', gamma=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1add0c09",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    bm.add_image(w_branch_mask, blending='additive')\n",
    "    bm.add_image(th_branch_mask, blending='additive')\n",
    "    bm.add_image(soma_mask, blending='additive')\n",
    "    bm.add_image(flood_mask, opacity=0.5, colormap='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79f339de",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Ближайшие точки скелета"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d798dd6d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "points_coords = idx[points]\n",
    "skelet_coords = idx[w_branch_mask]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de244701",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "points2connect = {}\n",
    "for point in points_coords:\n",
    "    point_dists = np.linalg.norm(point - skelet_coords, axis=-1)\n",
    "    if len(points2connect.keys()) > 0:\n",
    "        point_n_dists = np.linalg.norm(point - np.array(list(points2connect.keys())), axis=-1)\n",
    "        if point_n_dists.min() < 1.5:\n",
    "            continue\n",
    "    if point_dists.min() > 1.5:\n",
    "        coord = skelet_coords[point_dists.argmin()]\n",
    "        points2connect[tuple(point)] = tuple(coord)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c84691",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(points_coords.shape, len(points2connect.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b1b5a8b",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## Добавление соединений"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9af22a8",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# from networkx.algorithms.shortest_paths.weighted import multi_source_dijkstra\n",
    "full_mask = branch_mask.copy()\n",
    "for s, e in points2connect.items():\n",
    "#     d, p = multi_source_dijkstra(G3, [tuple(n) for n in skelet_coords], s, 'weight')\n",
    "    path = nx.shortest_path(G3, s, e, 'weight')\n",
    "    path_coords = np.array(path)\n",
    "    full_mask[path_coords[:,0], path_coords[:, 1], path_coords[:,2]] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8998abc",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "bm.add_image(full_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48292645",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# Многомасштабный гессиан"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83cdda28",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from skimage.filters import ridges as skridges"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8911cc7",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ad9747a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "https://content.iospress.com/download/bio-medical-materials-and-engineering/bme1149?id=bio-medical-materials-and-engineering%2Fbme1149\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f951311",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def vesselness(img, sigma, alpha= 0.5, beta=0.5, cc=20):\n",
    "    #H = skf.hessian_matrix(img,sigma)\n",
    "    #Hl = skf.hessian_matrix_eigvals(H)\n",
    "    Hl = skridges.compute_hessian_eigenvalues(img, sigma)\n",
    "    sh = img.shape\n",
    "    v = np.zeros(sh)\n",
    "    for r in range(sh[0]):\n",
    "        for c in range(sh[1]):\n",
    "            for h in range(sh[2]):\n",
    "                arr = (Hl[0][r,c,h], Hl[1][r,c,h], Hl[2][r,c,h])\n",
    "\n",
    "                l1,l2,l3 = sorted(arr, key=abs)\n",
    "                \n",
    "                if l2 >= 0 or l3 >=0:\n",
    "                    v[r,c,h] = 0\n",
    "                else:\n",
    "                    Ra = np.abs(l2/l3)\n",
    "                    Rb = np.abs(l1/np.sqrt(np.abs(l2*l3)))\n",
    "                    S = np.sqrt(l1**2 + l2**2 + l3**2)\n",
    "                    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*cc**2)))\n",
    "    return v\n",
    "\n",
    "\n",
    "def multiscale_vesselness(img, sigma_start, sigma_stop, nlevels=50, reduce_fn=np.max):\n",
    "    return reduce_fn([vesselness(img, sigma)*sigma**2 for sigma in np.linspace(sigma_start,sigma_stop, nlevels)],0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10589b35",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def multiscale_sigma(img, sigma_start, sigma_stop, nlevels=50):\n",
    "    sout = np.zeros(img.shape)\n",
    "    hout = np.zeros(img.shape)\n",
    "    vout = np.zeros((*img.shape, 3, 3))\n",
    "    for sigma in tqdm(np.linspace(sigma_start, sigma_stop, nlevels)):\n",
    "        hcurr, vcurr = astro.hessian.sato3d(img, sigma, hessian_variant='gradient_of_smoothed', return_vectors=True)\n",
    "        hcurr *= sigma**2\n",
    "        mask = hcurr > hout\n",
    "        \n",
    "        hout[mask] = hcurr[mask]\n",
    "        sout[mask] = sigma\n",
    "        vout[mask] = vcurr[mask]\n",
    "    return hout, sout, vout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2498bae",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "h, s, v = multiscale_sigma(final_image, 0.5, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8d32cf3",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hs = napari.view_image(s, ndisplay=3)\n",
    "hs.add_image(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c9d1a6f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from matplotlib.colors import ListedColormap, LinearSegmentedColormap\n",
    "colors = [\"black\", \"\"]\n",
    "cmap1 = LinearSegmentedColormap.from_list(\"mycmap\", colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74f7318d",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hs.add_image(np.where(full_mask^soma_mask, s, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ada3c051",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "hs.add_image(final_image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e365cfe",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from hessian_vectors import add_hessian_vectors, sato2napari_vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5330065f",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "add_hessian_vectors(hs, *sato2napari_vectors(h, v))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "030f665a",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "v_n, h_n = sato2napari_vectors(h, v)\n",
    "v_n.shape\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e8e93f8",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "add_hessian_vectors(hs, v_n.reshape((*full_mask.shape,2,3))[full_mask], h_n.reshape(full_mask.shape)[full_mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38f0beba",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sigma_skel = np.where(full_mask, s, 0)\n",
    "hs.add_image(np.where(sigma_skel < 1, sigma_skel, 0), blending='additive', name='0.5-1')\n",
    "hs.add_image(np.where(np.where(sigma_skel > 1, sigma_skel, 0) < 2, sigma_skel, 0), blending='additive', name='1-2')\n",
    "hs.add_image(np.where(np.where(sigma_skel > 2, sigma_skel, 0) < 3, sigma_skel, 0), blending='additive', name='2-3')\n",
    "hs.add_image(np.where(sigma_skel > 3, sigma_skel, 0), blending='additive', name='3-4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "648ba0be",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "w.add_image(full_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffb92cd0",
   "metadata": {
    "papermill": {
     "duration": null,
     "end_time": null,
     "exception": null,
     "start_time": null,
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "jupytext": {
   "formats": "ipynb,Rmd"
  },
  "kernelspec": {
   "display_name": "venv",
   "language": "python",
   "name": "venv"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 7.489341,
   "end_time": "2022-07-21T10:39:25.710994",
   "environment_variables": {},
   "exception": null,
   "input_path": "Skeletonization.ipynb",
   "output_path": "output.ipynb",
   "parameters": {
    "HANDY": false,
    "filename": "/home/levtg/astro-morpho/data/1wk-both1-red-raw.pic",
    "sigma": 1
   },
   "start_time": "2022-07-21T10:39:18.221653",
   "version": "2.3.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}