--- a
+++ b/plot_3d_CT.ipynb
@@ -0,0 +1,270 @@
+{
+ "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
+}