{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "64cdaa1c", "metadata": {}, "source": [ "# Process geometry" ] }, { "cell_type": "code", "execution_count": 1, "id": "12040cb6-7143-404a-8e0c-5da728d0a8a4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: vtk in i:\\iisc-project\\venv\\lib\\site-packages (9.2.6)\n", "Requirement already satisfied: matplotlib>=2.0.0 in i:\\iisc-project\\venv\\lib\\site-packages (from vtk) (3.8.0)\n", "Requirement already satisfied: contourpy>=1.0.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.0.7)\n", "Requirement already satisfied: cycler>=0.10 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (0.11.0)\n", "Requirement already satisfied: fonttools>=4.22.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (4.38.0)\n", "Requirement already satisfied: kiwisolver>=1.0.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.4.4)\n", "Requirement already satisfied: numpy<2,>=1.21 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (1.24.2)\n", "Requirement already satisfied: packaging>=20.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (23.0)\n", "Requirement already satisfied: pillow>=6.2.0 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (9.4.0)\n", "Requirement already satisfied: pyparsing>=2.3.1 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (3.0.9)\n", "Requirement already satisfied: python-dateutil>=2.7 in i:\\iisc-project\\venv\\lib\\site-packages (from matplotlib>=2.0.0->vtk) (2.8.2)\n", "Requirement already satisfied: six>=1.5 in i:\\iisc-project\\venv\\lib\\site-packages (from python-dateutil>=2.7->matplotlib>=2.0.0->vtk) (1.16.0)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n", "[notice] A new release of pip is available: 23.0 -> 24.0\n", "[notice] To update, run: python.exe -m pip install --upgrade pip\n" ] } ], "source": [ "!pip install vtk" ] }, { "cell_type": "code", "execution_count": 2, "id": "d78f6cc0", "metadata": { "vscode": { "languageId": "plaintext" } }, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'cv2'", "output_type": "error", "traceback": [ "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", "\u001B[1;31mModuleNotFoundError\u001B[0m Traceback (most recent call last)", "Cell \u001B[1;32mIn[2], line 18\u001B[0m\n\u001B[0;32m 16\u001B[0m \u001B[38;5;28;01mimport\u001B[39;00m \u001B[38;5;21;01mmatplotlib\u001B[39;00m\n\u001B[0;32m 17\u001B[0m \u001B[38;5;28;01mimport\u001B[39;00m \u001B[38;5;21;01mmath\u001B[39;00m\n\u001B[1;32m---> 18\u001B[0m \u001B[38;5;28;01mimport\u001B[39;00m \u001B[38;5;21;01mcv2\u001B[39;00m\n\u001B[0;32m 20\u001B[0m \u001B[38;5;28;01mfrom\u001B[39;00m \u001B[38;5;21;01mtyping\u001B[39;00m \u001B[38;5;28;01mimport\u001B[39;00m Dict, Any, List\n\u001B[0;32m 21\u001B[0m \u001B[38;5;28;01mfrom\u001B[39;00m \u001B[38;5;21;01mpathlib\u001B[39;00m \u001B[38;5;28;01mimport\u001B[39;00m Path\n", "\u001B[1;31mModuleNotFoundError\u001B[0m: No module named 'cv2'" ] } ], "source": [ "\"\"\"\n", "- computing quantification results of a network model via VMTK\n", "- created by Deepa\n", "- 25/06/2020\n", "- interpreter condapy36\n", "\n", "The output vtk files network_length and network_radius can be imported in Paraview\n", "to color the branch segments and visualize the length and radii distributions\n", "\n", "References:\n", "----------\n", "https://mail.google.com/mail/u/0/?tab=rm&ogbl#search/vmtk/FMfcgxwHMsVGfTPNDjDdffjtNQnSKbkC\n", "In short: load the vtk file, use tube filter\n", "Properties tab: Color edges by: Lengths/radius\n", "\"\"\"\n", "import os\n", "import vtk\n", "import vtk.numpy_interface.dataset_adapter as dsa\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "from vmtk import pypes\n", "from pprint import pprint\n", "from vmtk import vmtkscripts\n", "from vtk.util.numpy_support import vtk_to_numpy\n", "from collections import OrderedDict\n", "from typing import List\n", "from matpancreas.settings_model import RESULTS_DIR\n", "from matpancreas.utils.io_utils_py import write_output\n", "from matpancreas.utils.graph_utils_py import plot_histogram\n", "\n", "\n", "def write_ploydata_output(surface, file: str) -> None:\n", " \"\"\"\n", " write polyData output: for visualization in paraview\n", " :param surface:\n", " :param file:\n", " :return:\n", " \"\"\"\n", " surfacewriter = vmtkscripts.vmtkSurfaceWriter()\n", " surfacewriter.Surface = surface\n", " surfacewriter.OutputFileName = f\"network_{file}.vtp\"\n", " surfacewriter.Execute()\n", "\n", "\n", "def extract_segment_lengths(networkExtraction) -> List:\n", " \"\"\"\n", " Generates model output containing segment lengths in CellData\n", " :param newtork:\n", " :return:\n", " \"\"\"\n", " network = networkExtraction.Network\n", "\n", " centerlinesGeometry = vmtkscripts.vmtkCenterlineGeometry()\n", " centerlinesGeometry.Centerlines = network\n", " centerlinesGeometry.Execute()\n", " centerline = centerlinesGeometry.Centerlines\n", "\n", " # output\n", " write_ploydata_output(surface=centerline, file=\"length\")\n", "\n", " # table data\n", " centerline_dsa = dsa.WrapDataObject(centerline)\n", " length = []\n", " for edgeId in range(dsa.WrapDataObject(network).GetNumberOfCells()):\n", " length.append(centerline_dsa.CellData['Length'][edgeId])\n", "\n", " # plot_histogram(data=length, label='length', binwidth=5)\n", " return length\n", "\n", "\n", "def extract_segment_radius(networkExtraction) -> List:\n", " \"\"\"\n", " Generates model output containing segment radii in CellData\n", " :param newtork:\n", " :return:\n", " \"\"\"\n", " graph = networkExtraction.GraphLayout\n", "\n", " # output\n", " write_ploydata_output(surface=graph, file=\"radius\")\n", "\n", " # table data\n", " graph_dsa = dsa.WrapDataObject(graph)\n", " radius = []\n", " for edgeId in range(dsa.WrapDataObject(graph).GetNumberOfCells()):\n", " radius.append(graph_dsa.CellData['Radius'][edgeId])\n", " # plot_histogram(data=radius, label='radius', binwidth=1)\n", " return radius\n", "\n", "\n", "def get_coordinates(networkExtraction):\n", " \"\"\"\n", " get the x,y,z coordinates of head and tail nodes in graph\n", " :param vtkarray:\n", " :param graph:\n", " :return: df containing xyz coordinates\n", " \"\"\"\n", " graph = networkExtraction.GraphLayout\n", "\n", " cleaner = vtk.vtkCleanPolyData()\n", " cleaner.SetInputData(graph) # set the data to processed\n", " cleaner.Update() # execute the algorithm\n", " cleanedGraph = cleaner.GetOutput()\n", "\n", " cleanedGraph_dsa = dsa.WrapDataObject(graph)\n", "\n", " start = []\n", " end = []\n", "\n", " for edgeId in range(cleanedGraph.GetNumberOfCells()):\n", " edge = cleanedGraph.GetCell(edgeId) # this will be a vtkLine, which only has two points\n", "\n", " edgeNode0 = edge.GetPointId(0)\n", " edgeNode1 = edge.GetPointId(1)\n", " start.append(edgeNode0)\n", " end.append(edgeNode1)\n", "\n", " g = pd.DataFrame({'tail': start, 'head': end})\n", "\n", " pos = vtk_to_numpy(cleanedGraph_dsa.Points)\n", " nbranch = cleanedGraph.GetNumberOfCells()\n", " pos = pos.reshape(nbranch, 2, 3)\n", "\n", " xyz = OrderedDict()\n", " for i in range(len(end)):\n", " xyz[f'{start[i]}'] = pos[i][0]\n", " xyz[f'{end[i]}'] = pos[i][1]\n", "\n", " pos = pd.DataFrame(xyz)\n", " pos = pos.transpose()\n", " pos.columns = ['xpos', 'ypos', 'zpos']\n", " return g, pos\n", "\n", "\n", "if __name__ == '__main__':\n", "\n", " # previous input: LN1_crop1.stl islet.stl\n", " # cut open\n", " # vmtkCommand = '''vmtksurfaceclipper -ifile isletv3.stl -ofile isletv3_clipped.vtp '''\n", " #\n", " # p = pypes.Pype()\n", " # p.SetArgumentsString(vmtkCommand)\n", " # p.ParseArguments()\n", " # p.Execute()\n", " # exit()\n", " reader = vmtkscripts.vmtkSurfaceReader()\n", " reader.InputFileName = 'tumor_clipped.vtp'\n", " reader.Execute()\n", "\n", " # network extraction\n", " networkExtraction = vmtkscripts.vmtkNetworkExtraction()\n", " networkExtraction.Surface = reader.Surface\n", " networkExtraction.Execute()\n", "\n", " # get coordinates of terminal/free ends and N-furcation points in network\n", " g, pos = get_coordinates(networkExtraction=networkExtraction)\n", " pprint(pos)\n", "\n", " # generate ployData: lengths\n", " g['l'] = extract_segment_lengths(networkExtraction=networkExtraction)\n", "\n", " # generate ployData: radii\n", " g['r'] = extract_segment_radius(networkExtraction=networkExtraction)\n", "\n", " write_output(d={'graph': g, 'pos': pos}, file='data_vmtk')\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.11.1" } }, "nbformat": 4, "nbformat_minor": 5 }