[9173ee]: / examples / create meshes / Create Muscle Mesh Jan.28.2022.ipynb

Download this file

290 lines (289 with data), 7.4 kB

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f11c51b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymskt as mskt\n",
    "import SimpleITK as sitk\n",
    "import numpy as np\n",
    "import os\n",
    "\n",
    "from pymskt.mesh import Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "20bbd34d",
   "metadata": {},
   "outputs": [],
   "source": [
    "location_seg = '../../data/muscles/segmentation.nii.gz'\n",
    "location_save = os.path.expanduser('~/Downloads')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c43ac4d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Use of `point_arrays` is deprecated. Use `point_data` instead.\n",
      "Use of `cell_arrays` is deprecated. Use `cell_data` instead.\n"
     ]
    }
   ],
   "source": [
    "# Loading the image in first so we can see what the labels are. \n",
    "# This can be skipped and a path provided directly if we know\n",
    "# the labels of interest. \n",
    "seg_image = sitk.ReadImage(location_seg)\n",
    "array = sitk.GetArrayFromImage(seg_image)\n",
    "seg_labels = np.unique(array)\n",
    "\n",
    "mesh1 = Mesh(\n",
    "    seg_image=seg_image,\n",
    "    label_idx=seg_labels[1]\n",
    "    \n",
    ")\n",
    "\n",
    "# Below is the command to create the mesh\n",
    "# all of the defaults inputs parameters are being listed \n",
    "# so that the available options are shown. \n",
    "mesh1.create_mesh(\n",
    "    smooth_image=True,            # I suggest leaving this on to help create smooth surfaces.   \n",
    "    smooth_image_var=0.3125/2,    # this is the variance of a gaussian filter - can probably be bigger for muscles.\n",
    "    marching_cubes_threshold=0.5, # this probably shouldnt be changed. \n",
    "    label_idx=None,               # Can specify this here instead of above if you want. \n",
    "    min_n_pixels=None,            # if there is a minimum number of pixels you want for it to create a mesh. \n",
    ")\n",
    "\n",
    "mesh1.resample_surface(clusters=10000) # Resample surface to be a specified # of vertices (w/ in a small margin of error)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4c3a458",
   "metadata": {},
   "source": [
    "## Save the mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7bd7c341",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh1.save_mesh(os.path.join(location_save, 'mesh1.vtk')) # Save .vtk version of the mesh\n",
    "mesh1.save_mesh(os.path.join(location_save, 'mesh1.stl')) # Save .stl version of the mesh. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "394b11b0",
   "metadata": {},
   "source": [
    "# If you want to iterate over all of the labels and save meshes: "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4c39deb",
   "metadata": {},
   "source": [
    "### First create a dictionary to store all of the meshes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4368b846",
   "metadata": {},
   "outputs": [],
   "source": [
    "dict_meshes = {}\n",
    "for label_idx in seg_labels[1:]:\n",
    "    dict_meshes[label_idx] = mskt.mesh.Mesh(\n",
    "        seg_image=seg_image,\n",
    "        label_idx=label_idx\n",
    "\n",
    "    )\n",
    "    \n",
    "    dict_meshes[label_idx].create_mesh(smooth_image_var=3.0) # I ADDED MORE SMOOTHING FOR THESE\n",
    "    dict_meshes[label_idx].resample_surface(clusters=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc29d8d",
   "metadata": {},
   "source": [
    "### Save each mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "f5d6738e",
   "metadata": {},
   "outputs": [],
   "source": [
    "for key, mesh in dict_meshes.items():\n",
    "    mesh.save_mesh(os.path.join(location_save, f'muscle_mesh{key}.vtk')) # Save vtk version\n",
    "    mesh.save_mesh(os.path.join(location_save, f'muscle_mesh{key}.stl')) # Save stl version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0a66862",
   "metadata": {},
   "source": [
    "# Below is the viewer that should work to see these in jupyter notebooks\n",
    "- This wasnt working on my computer for some reason. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b9f350ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itkwidgets import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "505700e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0, 14, 16, 18, 20, 22, 24, 26], dtype=uint16)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seg_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "52169ea2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5f7605d7f4274c78bbcbcb884feb90f8",
       "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 expects geometries to be in a list: \n",
    "geometries = [\n",
    "    dict_meshes[14].mesh  # the Mesh object has the actual mesh inside of it at Mesh.mesh\n",
    "]\n",
    "view(geometries=geometries)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c6b9c5a",
   "metadata": {},
   "source": [
    "### Show multiple meshes at once: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ddf2bed6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b9275aa069ce432b9ab363874895c7d4",
       "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 expects geometries to be in a list: \n",
    "geometries = [\n",
    "    dict_meshes[14].mesh,  # the Mesh object has the actual mesh inside of it at Mesh.mesh\n",
    "    dict_meshes[16].mesh,\n",
    "    dict_meshes[18].mesh,\n",
    "    dict_meshes[20].mesh,\n",
    "    dict_meshes[22].mesh,\n",
    "    dict_meshes[24].mesh,\n",
    "    dict_meshes[26].mesh\n",
    "]\n",
    "view(geometries=geometries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "397e00d4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}