271 lines (270 with data), 12.4 kB

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Borrowed some parts from Kaggle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "from os.path import dirname, join\n",
    "from os import listdir\n",
    "from pprint import pprint\n",
    "\n",
    "import pydicom\n",
    "\n",
    "from pydicom.filereader import read_dicomdir\n",
    "\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "from skimage import measure, morphology\n",
    "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n",
    "\n",
    "from skimage.segmentation import clear_border\n",
    "from scipy import ndimage as ndi\n",
    "from skimage.filters import roberts, sobel\n",
    "\n",
    "import cv2\n",
    "import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rescale_patient_images(images_zyx, org_spacing_xyz, target_voxel_mm, is_mask_image=True, verbose=False):\n",
    "    if verbose:\n",
    "        print(\"Spacing: \", org_spacing_xyz)\n",
    "        print(\"Shape: \", images_zyx.shape)\n",
    "\n",
    "    # print \"Resizing dim z\"\n",
    "    resize_x = 1.0\n",
    "    resize_y = float(org_spacing_xyz[2]) / float(target_voxel_mm)\n",
    "    interpolation = cv2.INTER_NEAREST if is_mask_image else cv2.INTER_LINEAR\n",
    "    res = cv2.resize(images_zyx, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)  # opencv assumes y, x, channels umpy array, so y = z pfff\n",
    "    # print \"Shape is now : \", res.shape\n",
    "\n",
    "    res = res.swapaxes(0, 2)\n",
    "    res = res.swapaxes(0, 1)\n",
    "    # print \"Shape: \", res.shape\n",
    "    resize_x = float(org_spacing_xyz[0]) / float(target_voxel_mm)\n",
    "    resize_y = float(org_spacing_xyz[1]) / float(target_voxel_mm)\n",
    "\n",
    "    # cv2 can handle max 512 channels..\n",
    "    if res.shape[2] > 512:\n",
    "        res = res.swapaxes(0, 2)\n",
    "        res1 = res[:256]\n",
    "        res2 = res[256:]\n",
    "        res1 = res1.swapaxes(0, 2)\n",
    "        res2 = res2.swapaxes(0, 2)\n",
    "        res1 = cv2.resize(res1, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)\n",
    "        res2 = cv2.resize(res2, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)\n",
    "        res1 = res1.swapaxes(0, 2)\n",
    "        res2 = res2.swapaxes(0, 2)\n",
    "        res = numpy.vstack([res1, res2])\n",
    "        res = res.swapaxes(0, 2)\n",
    "    else:\n",
    "        res = cv2.resize(res, dsize=None, fx=resize_x, fy=resize_y, interpolation=interpolation)\n",
    "\n",
    "    # channels = cv2.split(res)\n",
    "    # resized_channels = []\n",
    "    # for channel in  channels:\n",
    "    #     channel = cv2.resize(channel, dsize=None, fx=resize_x, fy=resize_y)\n",
    "    #     resized_channels.append(channel)\n",
    "    # res = cv2.merge(resized_channels)\n",
    "    # print \"Shape after resize: \", res.shape\n",
    "    res = res.swapaxes(0, 2)\n",
    "    res = res.swapaxes(2, 1)\n",
    "    if verbose:\n",
    "        print(\"Shape after: \", res.shape)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "filelist = glob.glob(\"*.png\")\n",
    "\n",
    "imlist = [plt.imread(file) for file in filelist]\n",
    "    \n",
    "R_slices = np.stack([im[:,:,0] for im in imlist])\n",
    "G_slices = np.stack([im[:,:,1] for im in imlist])\n",
    "B_slices = np.stack([im[:,:,2] for im in imlist])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(166, 1102, 1102)\n"
     ]
    }
   ],
   "source": [
    "R_slices = np.array(R_slices, dtype=np.int16)\n",
    "print(R_slices.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape before resampling\t (166, 1102, 1102)\n"
     ]
    }
   ],
   "source": [
    "def resample(image, new_spacing=[1,1,1]):\n",
    "    # Determine current pixel spacing\n",
    "    pixel_spacing = ['2.5','0.310546875', '0.310546875']\n",
    "    spacing = map(float, (pixel_spacing))\n",
    "\n",
    "    spacing = np.array(list(spacing))\n",
    "\n",
    "    resize_factor = spacing / new_spacing\n",
    "    new_real_shape = image.shape * resize_factor\n",
    "    new_shape = np.round(new_real_shape)\n",
    "    real_resize_factor = new_shape / image.shape\n",
    "    new_spacing = spacing / real_resize_factor\n",
    "    \n",
    "    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)\n",
    "    \n",
    "    return image, new_spacing\n",
    "\n",
    "print(\"Shape before resampling\\t\", R_slices.shape)\n",
    "imgs_after_resamp, spacing = resample(R_slices, [1,1,1])\n",
    "print(\"Shape after resampling\\t\", imgs_after_resamp.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_mesh(image, threshold=-300, step_size=1):\n",
    "\n",
    "    p = image.transpose(2,1,0)\n",
    "    \n",
    "    verts, faces, norm, val = measure.marching_cubes(p, threshold, step_size=step_size, allow_degenerate=True) \n",
    "    return verts, faces\n",
    "\n",
    "\n",
    "def plt_3d(verts, faces):\n",
    "\n",
    "    x,y,z = zip(*verts) \n",
    "    fig = plt.figure(figsize=(10, 10))\n",
    "    ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "    # Fancy indexing: `verts[faces]` to generate a collection of triangles\n",
    "    mesh = Poly3DCollection(verts[faces], linewidths=0.05, alpha=1)\n",
    "    face_color = [1, 1, 0.9]\n",
    "    mesh.set_facecolor(face_color)\n",
    "    ax.add_collection3d(mesh)\n",
    "\n",
    "    ax.set_xlim(0, max(x))\n",
    "    ax.set_ylim(0, max(y))\n",
    "    ax.set_zlim(0, max(z))\n",
    "    ax.set_axis_bgcolor((0.7, 0.7, 0.7))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Surface level must be within volume data range.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-5f5a4cb3c07e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m#R_slices = rescale_patient_images(R_slices, pixel_spacing, TARGET_VOXEL_MM)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mverts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfaces\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmake_mesh\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mR_slices\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0mplt_3d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mverts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfaces\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-21-7817fcf65383>\u001b[0m in \u001b[0;36mmake_mesh\u001b[0;34m(image, threshold, step_size)\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0mp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranspose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m     \u001b[0mverts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfaces\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmeasure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmarching_cubes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mthreshold\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstep_size\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mstep_size\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mallow_degenerate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mverts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfaces\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/skimage/measure/_marching_cubes_lewiner.py\u001b[0m in \u001b[0;36mmarching_cubes\u001b[0;34m(volume, level, spacing, gradient_direction, step_size, allow_degenerate, use_classic)\u001b[0m\n\u001b[1;32m    174\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    175\u001b[0m     return marching_cubes_lewiner(volume, level, spacing, gradient_direction,\n\u001b[0;32m--> 176\u001b[0;31m                                   step_size, allow_degenerate, use_classic)\n\u001b[0m\u001b[1;32m    177\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    178\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/skimage/measure/_marching_cubes_lewiner.py\u001b[0m in \u001b[0;36mmarching_cubes_lewiner\u001b[0;34m(volume, level, spacing, gradient_direction, step_size, allow_degenerate, use_classic)\u001b[0m\n\u001b[1;32m    197\u001b[0m         \u001b[0mlevel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlevel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    198\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mlevel\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mvolume\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mlevel\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mvolume\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 199\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Surface level must be within volume data range.\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    200\u001b[0m     \u001b[0;31m# spacing\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    201\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspacing\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: Surface level must be within volume data range."
     ]
    }
   ],
   "source": [
    "pixel_spacing = ['0.310546875', '0.310546875', '3.0']\n",
    "TARGET_VOXEL_MM = 1.00\n",
    "\n",
    "\n",
    "#R_slices = rescale_patient_images(R_slices, pixel_spacing, TARGET_VOXEL_MM)\n",
    "verts, faces = make_mesh(R_slices)\n",
    "plt_3d(verts, faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}