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
}