AutoCenterlines-Take2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Compute the associated average outward flux\n",
    "\n",
    "To have a good quality gradient of the distance map, perform a light smooth over it. Define convolution kernels in each direction and use them recursively. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib ipympl\n",
    "import matplotlib.pyplot as plt\n",
    "from vmtk import vmtkscripts\n",
    "from vmtk import vtkvmtk\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "from vtk.numpy_interface import dataset_adapter as dsa\n",
    "import random\n",
    "import networkx as nx\n",
    "import vtk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from vtkdatawidgets.vtk_binding import VtkJupyterBridge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "bridge = VtkJupyterBridge()\n",
    "sphere = vtk.vtkSphereSource()\n",
    "bridge.SetInputConnection(sphere.GetOutputPort())\n",
    "bridge.Update()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "imgr = vmtkscripts.vmtkImageReader()\n",
    "imgr.InputFileName = os.path.join(os.getcwd(), 'test-surf-normal-hole.mha')\n",
    "imgr.Execute()\n",
    "\n",
    "mc = vmtkscripts.vmtkMarchingCubes()\n",
    "mc.Image = imgr.Image\n",
    "mc.Level = 1000.0\n",
    "mc.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading STL surface file.\n"
     ]
    }
   ],
   "source": [
    "surfr = vmtkscripts.vmtkSurfaceReader()\n",
    "surfr.InputFileName = os.path.join(os.getcwd(), 'realtest.stl')\n",
    "surfr.Execute()\n",
    "\n",
    "surf = surfr.Surface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "fedges = vtk.vtkFeatureEdges()\n",
    "fedges.BoundaryEdgesOn()\n",
    "fedges.FeatureEdgesOff()\n",
    "fedges.ManifoldEdgesOff()\n",
    "fedges.SetInputData(surf)\n",
    "fedges.Update()\n",
    "\n",
    "ofedges = fedges.GetOutput()\n",
    "\n",
    "numEdges = ofedges.GetNumberOfPoints()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if numEdges != 0:\n",
    "    tempcapper = vmtkscripts.vmtkSurfaceCapper()\n",
    "    tempcapper.Surface = surf\n",
    "    tempcapper.Interactive = 0\n",
    "    tempcapper.Execute()\n",
    "    \n",
    "    networkSurface = tempcapper.Surface\n",
    "else:\n",
    "    networkSurface = surf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Progress: 100%Progress: 100%\n"
     ]
    }
   ],
   "source": [
    "numCells = surf.GetNumberOfCells()\n",
    "cellToDelete = random.randrange(0, numCells-1)\n",
    "\n",
    "networkSurface.BuildLinks()\n",
    "networkSurface.DeleteCell(cellToDelete)\n",
    "networkSurface.RemoveDeletedCells()\n",
    "\n",
    "net = vmtkscripts.vmtkNetworkExtraction()\n",
    "net.Surface = networkSurface\n",
    "net.AdvancementRatio = 1.001\n",
    "net.Execute()\n",
    "\n",
    "network = net.Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wrapping vtkPolyData object\n",
      "converting cell data: \n",
      "Topology\n",
      "converting points\n",
      "converting point data: \n",
      "Radius\n",
      "converting cell connectivity list\n"
     ]
    }
   ],
   "source": [
    "convert = vmtkscripts.vmtkCenterlinesToNumpy()\n",
    "convert.Centerlines = network\n",
    "convert.Execute()\n",
    "ad = convert.ArrayDict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellDataTopology = ad['CellData']['Topology']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(102, 2)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cellDataTopology.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "fullTopology = np.array(cellDataTopology)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(102, 2)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fullTopology.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "topologyToKeep = np.where(fullTopology[:,0] != 0)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "topol = fullTopology[topologyToKeep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(101, 2)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topol.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2],\n",
       "       [ 2, -1],\n",
       "       [ 2,  3],\n",
       "       [ 3,  4],\n",
       "       [ 3,  5],\n",
       "       [ 4, -1],\n",
       "       [ 5,  6],\n",
       "       [ 5, -1],\n",
       "       [ 6,  8],\n",
       "       [ 6, -1],\n",
       "       [ 7,  9],\n",
       "       [ 7, 52],\n",
       "       [ 8, 10],\n",
       "       [ 9, 11],\n",
       "       [ 9, 12],\n",
       "       [10, 14],\n",
       "       [10, 15],\n",
       "       [11, 16],\n",
       "       [11, -1],\n",
       "       [12, 17],\n",
       "       [12, 18],\n",
       "       [13, 19],\n",
       "       [13, 20],\n",
       "       [14, -1],\n",
       "       [14, -1],\n",
       "       [15, 21],\n",
       "       [16, 22],\n",
       "       [16, -1],\n",
       "       [17, 23],\n",
       "       [17, 24],\n",
       "       [18, -1],\n",
       "       [18, -1],\n",
       "       [19, -1],\n",
       "       [19, -1],\n",
       "       [20, -1],\n",
       "       [21, 25],\n",
       "       [21, -1],\n",
       "       [22, -1],\n",
       "       [22, -1],\n",
       "       [23, 26],\n",
       "       [24, 27],\n",
       "       [24, -1],\n",
       "       [25, 29],\n",
       "       [25, 28],\n",
       "       [26, 31],\n",
       "       [26, -1],\n",
       "       [27, 32],\n",
       "       [28, 33],\n",
       "       [29, -1],\n",
       "       [30, 34],\n",
       "       [30, 35],\n",
       "       [31, 36],\n",
       "       [31, 37],\n",
       "       [32, 38],\n",
       "       [32, -1],\n",
       "       [34, 39],\n",
       "       [34, -1],\n",
       "       [35, -1],\n",
       "       [35, 40],\n",
       "       [36, 41],\n",
       "       [36, 42],\n",
       "       [37, -1],\n",
       "       [37, -1],\n",
       "       [38, 44],\n",
       "       [38, -1],\n",
       "       [39, -1],\n",
       "       [40, -1],\n",
       "       [40, -1],\n",
       "       [41, 45],\n",
       "       [42, -1],\n",
       "       [42, 46],\n",
       "       [43, 47],\n",
       "       [43, -1],\n",
       "       [44, 53],\n",
       "       [44, 54],\n",
       "       [45, 48],\n",
       "       [48, -1],\n",
       "       [48, -1],\n",
       "       [49, -1],\n",
       "       [49, -1],\n",
       "       [51, -1],\n",
       "       [51, -1],\n",
       "       [ 1, -1],\n",
       "       [ 4,  7],\n",
       "       [ 8, 13],\n",
       "       [20, -1],\n",
       "       [23, 30],\n",
       "       [28, 55],\n",
       "       [29, 39],\n",
       "       [33, -1],\n",
       "       [36, -1],\n",
       "       [43, 56],\n",
       "       [-1, 57],\n",
       "       [45, 49],\n",
       "       [47, -1],\n",
       "       [47, -1],\n",
       "       [50, 58],\n",
       "       [-1, 59],\n",
       "       [50, 51],\n",
       "       [50, -1],\n",
       "       [-1, 60]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "topol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "nodes = np.arange(1, topol.shape[0]+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,\n",
       "       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,\n",
       "       35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n",
       "       52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68,\n",
       "       69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85,\n",
       "       86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "topolWithNodes = np.copy(topol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "endEdgesToReplace = np.where(topolWithNodes[:, :] == -1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7, 10, 12, 13, 14, 15, 16, 17, 18, 19, 20, 29, 33, 36, 38, 39, 40,\n",
       "       43, 45, 47, 49, 51, 52, 54, 55, 56, 57, 60, 63, 68, 70, 72, 75, 78,\n",
       "       79, 81, 85, 87, 88, 89, 90, 91, 92, 93, 94, 96])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "endEdgesToReplace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx in endEdgesToReplace:\n",
    "    if topolWithNodes[idx, 1] == -1:\n",
    "        topolWithNodes[idx, 1] = nodes[idx]\n",
    "    elif topolWithNodes[idx, 0] == -1:\n",
    "        topolWithNodes[idx, 0] = nodes[idx]\n",
    "    else:\n",
    "        print('bad')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.Graph()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "G.add_edges_from(topolWithNodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "97"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "G.number_of_edges()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rick/miniconda3/envs/vmtkplayground/lib/python3.6/site-packages/matplotlib/__init__.py:942: MatplotlibDeprecationWarning: nbagg.transparent is deprecated and ignored. Use figure.facecolor instead.\n",
      "  mplDeprecation)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6516bbac645b4b7191c8a403ac1f862f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pos = nx.nx_agraph.graphviz_layout(G)\n",
    "nx.draw(G, pos=pos, with_labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for nodeidx in G:\n",
    "    connections = [i for i in nx.all_neighbors(G, nodeidx)]\n",
    "    if len(connections) < 2:\n",
    "        print(f'{nodeidx}: {connections}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for node in G.nodes():\n",
    "    print(f'{node}, {nx.}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[i for i in C]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = nx.nx_agraph.graphviz_layout(C)\n",
    "nx.draw(C, pos=pos, with_labels=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endCellIds = []\n",
    "startCellIds = []\n",
    "for cellId, cellTop in enumerate(cellDataTopology):\n",
    "    if cellTop[1] == -1:\n",
    "        endCellIds.append(cellId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endCellIds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endPointIds = []\n",
    "for cellId in endCellIds:  \n",
    "    endPointIds.append(ad['CellData']['CellPointIds'][cellId][-2])\n",
    "\n",
    "endPointCoordinates = []\n",
    "for pointId in endPointIds:\n",
    "    endPointCoordinates.append(ad['Points'][pointId])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endPointCoordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endPointCoordinates = np.array(endPointCoordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "endPointCoordinates.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sourcePoints = endPointCoordinates[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sourcePoints.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flatSourcePoints = sourcePoints.flatten(order='C').tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sourcePoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targetPoints = endPointCoordinates[0:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targetPoints.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flatTargetPoints = targetPoints.flatten(order='C').tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targetPoints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl = vmtkscripts.vmtkCenterlines()\n",
    "cl.Surface = surf\n",
    "cl.SeedSelectorName = 'pointlist'\n",
    "cl.SourcePoints = flatSourcePoints\n",
    "cl.TargetPoints = flatTargetPoints\n",
    "cl.Execute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numberOfCells = surfr.Surface.GetNumberOfCells()\n",
    "cellToRemove = numberOfCells - 10\n",
    "cellIds = vtk.vtkIdTypeArray()\n",
    "cellIds.SetNumberOfComponents(1)\n",
    "cellIds.InsertNextValue(cellToRemove)\n",
    "\n",
    "selectionNode = vtk.vtkSelectionNode()\n",
    "selectionNode.SetFieldType(selectionNode.CELL)\n",
    "selectionNode.SetContentType(selectionNode.INDICES)\n",
    "selectionNode.SetSelectionList(cellIds)\n",
    "\n",
    "selection = vtk.vtkSelection()\n",
    "selection.AddNode(selectionNode)\n",
    "\n",
    "extractSelection = vtk.vtkExtractSelection()\n",
    "extractSelection.SetInputData(0, surf)\n",
    "extractSelection.SetInputData(1, selection)\n",
    "extractSelection.Update()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "select = vmtkscripts.vmtkSurfaceConnectivitySelector()\n",
    "select.Surface = surfr.Surface\n",
    "select.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sv = vmtkscripts.vmtkSurfaceViewer()\n",
    "sv.Surface = surf\n",
    "sv.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edit = vmtkscripts.vmtkNetworkEditor()\n",
    "edit.Network = net.Network\n",
    "edit.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv = vmtkscripts.vmtkCenterlineViewer()\n",
    "cv.Centerlines = network\n",
    "cv.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rend = vmtkscripts.vmtkRenderer()\n",
    "rend.Execute()\n",
    "\n",
    "view = vmtkscripts.vmtkSurfaceViewer()\n",
    "view.vmtkRenderer = rend.vmtkRenderer\n",
    "view.Opacity = 0.3\n",
    "view.Surface = surf\n",
    "view.Execute()\n",
    "\n",
    "view2 = vmtkscripts.vmtkCenterlineViewer()\n",
    "view2.vmtkRenderer = rend.vmtkRenderer\n",
    "view2.Centerlines = cl.Centerlines\n",
    "view2.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write = vmtkscripts.vmtkSurfaceWriter()\n",
    "write.Surface = net.Network\n",
    "write.OutputFileName = os.path.join(os.getcwd(), 'test-multi-hole.vtp')\n",
    "write.Execute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}