[9173ee]: / examples / view_thickness_calculations / thick_calc_visuals.ipynb

Download this file

290 lines (289 with data), 8.6 kB

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymskt as mskt\n",
    "from pymskt.mesh import BoneMesh, CartilageMesh\n",
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "location_seg = '../../data/right_knee_example.nrrd'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INTERSECTION IS:  2\n"
     ]
    }
   ],
   "source": [
    "fem_cart = CartilageMesh(\n",
    "    path_seg_image=location_seg,  # path to the segmentation iamge being used\n",
    "    label_idx=1,                  # what is the label of this bone\n",
    ")\n",
    "\n",
    "fem_cart.create_mesh()\n",
    "fem_cart.resample_surface(clusters=10000)\n",
    "\n",
    "# initiate the bone mesh object\n",
    "femur = BoneMesh(path_seg_image=location_seg,  # path to the segmentation iamge being used\n",
    "                 label_idx=5,                  # what is the label of this bone\n",
    "                #  list_cartilage_labels=[1])    # a list of labels for cartialge associated w/ this bone\n",
    "                list_cartilage_meshes=[fem_cart]\n",
    ")\n",
    "\n",
    "# Create the bone mesh\n",
    "femur.create_mesh(\n",
    "    smooth_image_var=1.0       # This is the variance of the gaussian filter applied to binary image b4 meshing\n",
    ")\n",
    "\n",
    "# Resample the bone surface to have a specified number of nodes. \n",
    "femur.resample_surface(\n",
    "    clusters=10000             # This is the number of nodes/vertices on the surface. It might vary slightly\n",
    ")\n",
    "\n",
    "# Calcualte cartialge thickness for the cartialge meshes associated with the bone\n",
    "femur.calc_cartilage_thickness(\n",
    "    image_smooth_var_cart=0.5   # variance for the gaussian filter applied to the binary cart meshes b4 smoothing. \n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_points = 100\n",
    "\n",
    "# Get 10 points on bone surface with cartilage thickness\n",
    "thickness = femur.get_scalar('thickness (mm)')\n",
    "cart_points = np.where(thickness > 0.0)\n",
    "# get 10 rand points. \n",
    "cart_points = np.random.choice(cart_points[0], n_points)\n",
    "\n",
    "# get points xyz positions \n",
    "bone_origin = femur.point_coords[cart_points]\n",
    "\n",
    "# Project vector normal to surface for these points. \n",
    "from pymskt.mesh.utils import get_surface_normals\n",
    "normals = get_surface_normals(femur.mesh)\n",
    "normals = normals.GetOutput().GetPointData().GetNormals()\n",
    "normals = np.asarray([normals.GetTuple(idx) for idx in cart_points])\n",
    "# Get intersection points on cartilage for these points.\n",
    "from pymskt.mesh.utils import get_obb_surface, get_intersect, get_arrow\n",
    "import pyvista as pv\n",
    "obb_cartilage = get_obb_surface(femur.list_cartilage_meshes[0].mesh)\n",
    "\n",
    "dict_vectors = {}\n",
    "\n",
    "length = 10\n",
    "length_opposite = 1\n",
    "\n",
    "for idx, point in enumerate(bone_origin):\n",
    "    vector = normals[idx]\n",
    "    start_point_ray = point - length_opposite * vector\n",
    "    end_point_ray = point + length * vector\n",
    "    points_intersect, cell_ids_intersect = get_intersect(obb_cartilage, start_point_ray, end_point_ray)\n",
    "\n",
    "    arrow = get_arrow(\n",
    "        vector,\n",
    "        start_point_ray,\n",
    "        scale=length + length_opposite,\n",
    "        tip_length=0.25,\n",
    "        tip_radius=0.1,\n",
    "        tip_resolution=20, \n",
    "        shaft_radius=0.01,\n",
    "        shaft_resolution=20\n",
    "    )\n",
    "\n",
    "    spheres = [pv.Sphere(center=point, radius=0.5) for point in points_intersect]\n",
    "\n",
    "    dict_vectors[idx] = {\n",
    "        'start': start_point_ray, \n",
    "        'end': end_point_ray, \n",
    "        'points_intersect': points_intersect,\n",
    "        'arrow': arrow,\n",
    "        'spheres': spheres\n",
    "    }\n",
    "\n",
    "\n",
    "# draw line for the vector between these two points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e1527c4403b84b4cb2ccbd730db0c247",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from itkwidgets import Viewer\n",
    "\n",
    "# Note that we have to get the mesh from the femur object `femur.mesh` to input it into itkwidgets. \n",
    "# In this notebook, you should be able to swap between scalars (thickness & labels)\n",
    "# thickness is the cartilage thickness in mm, labels is the cartilage label of the associated thickness value. \n",
    "Viewer(geometries=[femur.mesh, femur.list_cartilage_meshes[0].mesh, dict_vectors[0]['arrow']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4806d95894804b7bad37c63dfc78e29f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# View with just vectors\n",
    "list_geometries = [femur.mesh, femur.list_cartilage_meshes[0].mesh]\n",
    "list_arrows = [dict_vectors[idx]['arrow'] for idx in dict_vectors]\n",
    "list_geometries.extend(list_arrows)\n",
    "Viewer(geometries=list_geometries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fb56bc7726ab4593acad287437061148",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_vectors = 10\n",
    "\n",
    "# View with just vectors\n",
    "list_geometries = [femur.mesh, femur.list_cartilage_meshes[0].mesh]\n",
    "list_arrows = [dict_vectors[idx]['arrow'] for idx in range(n_vectors)]\n",
    "list_spheres = [dict_vectors[idx]['spheres'][i] for idx in range(n_vectors) for i in range(2)]\n",
    "list_geometries.extend(list_arrows)\n",
    "list_geometries.extend(list_spheres)\n",
    "\n",
    "Viewer(geometries=list_geometries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b1779dcbfd4c4f4191fe4d35daed58d2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Viewer(geometries=[{'vtkClass': 'vtkPolyData', 'points': {'vtkClass': 'vtkPoints', 'name': '_points', 'numberO…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Viewer(geometries=[femur.mesh, femur.list_cartilage_meshes[0].mesh])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "imaging",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}