day1-tut02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WMLMIA 2018\n",
    "## Tutorial 2 : DICOM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import scipy.io as sio\n",
    "import numpy as np\n",
    "from PIL import Image\n",
    "from myFcnCart2Pol import myFcnCart2Pol\n",
    "from skimage.segmentation import random_walker\n",
    "import skimage\n",
    "import scipy.signal as signal\n",
    "import scipy.ndimage.morphology as morph\n",
    "from myFcnPol2Cart import myFcnPol2Cart\n",
    "import pickle\n",
    "import os\n",
    "import dicom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imagePath = 'data/day1/tut01/Originaltotal/'\n",
    "gtPath = 'data/day1/tut01/GTtotal/'\n",
    "\n",
    "beta = 2000\n",
    "\n",
    "# Size of image in polar coordinates\n",
    "szPol = 256\n",
    "\n",
    "# Window size for feature extraction\n",
    "Kernelmat = np.array([3,5,7,13,17,21,25,30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading trained Ranfom-Forest model for segmentation\n",
    "pkl_file = open('data/day1/tut01/myfile.pkl','rb')\n",
    "rf = pickle.load(pkl_file,encoding='latin1')\n",
    "pkl_file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N = len(os.listdir(imagePath))\n",
    "for fileIdx in range(1,N+1):\n",
    "    start_time = time.time()      \n",
    "    if fileIdx < 10:    \n",
    "        filename = 'frame_01_000'+str(fileIdx)+'_003.png'\n",
    "    else:\n",
    "        filename = 'frame_01_00'+str(fileIdx)+'_003.png'\n",
    "    print('Image: '+str(fileIdx)+'/'+str(N+1))\n",
    "    OrgImage = np.asarray(Image.open(imagePath+filename).convert('L'));\n",
    "    GtImage = np.asarray(Image.open(gtPath+filename).convert('L'));\n",
    "    (szCart1, szCart2) = OrgImage.shape;\n",
    "    \n",
    "    # Creating volume tesnsor\n",
    "    if fileIdx==1:\n",
    "        Vol = np.empty((szCart1, szCart2, 3, N), dtype=np.uint8)        \n",
    "    \n",
    "    polarOrgImage = myFcnCart2Pol(OrgImage, szPol, szPol);\n",
    "    polarGtImage = myFcnCart2Pol(GtImage, szPol, szPol);\n",
    "    markers = np.zeros(polarOrgImage.shape, dtype=np.uint);   \n",
    "    markers[0,:] = 2;\n",
    "    markers[szPol-1,:] = 1;\n",
    "    labels = random_walker(skimage.img_as_float(polarOrgImage), markers, beta, mode='bf', copy=True, return_full_prob=True);   \n",
    "    for scaleCtr in range(0,len(Kernelmat)):\n",
    "        hScale = np.ones((Kernelmat[scaleCtr], Kernelmat[0]));\n",
    "        hScale = hScale/np.sum(hScale);\n",
    "        ConvMean  = signal.convolve2d(polarOrgImage, hScale, mode='same');\n",
    "        SqMean = signal.convolve2d(polarOrgImage**2, hScale, mode='same');\n",
    "        ConvVer = SqMean - ConvMean**2;\n",
    "        if scaleCtr == 0:\n",
    "            MeanSet = np.reshape(ConvMean, (1,np.product(ConvMean.shape)))[0];\n",
    "            VerSet = np.reshape(ConvVer, (1,np.product(ConvVer.shape)))[0];\n",
    "            ConfSet = np.reshape(labels[1,:,:], (1,np.product(labels[1,:,:].shape)))[0];\n",
    "        else:\n",
    "            MeanSet = np.vstack((MeanSet,np.reshape(ConvMean, (1,np.product(ConvMean.shape)))[0]));\n",
    "            VerSet = np.vstack((VerSet,np.reshape(ConvVer, (1,np.product(ConvVer.shape)))[0]));\n",
    "    TestingFeatureSet = np.zeros((np.product(polarOrgImage.shape),17));\n",
    "    MeanSet = MeanSet.transpose();\n",
    "    VerSet = VerSet.transpose();\n",
    "    TestingFeatureSet[:,0:16] = np.hstack((MeanSet,VerSet));\n",
    "    TestingFeatureSet[:,16] = ConfSet;\n",
    "    Pred = rf.predict(TestingFeatureSet[:,0:17]);\n",
    "    PredImage = np.reshape(Pred,(szPol,szPol));\n",
    "    \n",
    "    \n",
    "    PredBwLum = np.zeros(PredImage.shape);     \n",
    "    PredBwLum[PredImage == 1] = 1;\n",
    "    PredBwMed = np.zeros(PredImage.shape);     \n",
    "    PredBwMed[(PredImage ==1) | (PredImage ==2)] = 1;\n",
    "    \n",
    "    LumBw = skimage.morphology.dilation((skimage.morphology.erosion(PredBwLum, skimage.morphology.disk(2))), skimage.morphology.disk(2))\n",
    "    LumBw = skimage.morphology.label(LumBw); \n",
    "    MedBw = skimage.morphology.dilation((skimage.morphology.erosion(PredBwMed, skimage.morphology.disk(2))), skimage.morphology.disk(2))\n",
    "    MedBw = skimage.morphology.label(MedBw);\n",
    "    LumBw[np.nonzero(LumBw >1)] = 0\n",
    "    MedBw[np.nonzero(MedBw >1)] = 0\n",
    "    LumBw = morph.binary_fill_holes(LumBw)\n",
    "    MedBw = morph.binary_fill_holes(MedBw)\n",
    "    \n",
    "    Lum = myFcnPol2Cart(LumBw, szCart1, szCart2);\n",
    "    Med = myFcnPol2Cart(MedBw, szCart1, szCart2);\n",
    "    k = 15\n",
    "    \n",
    "    Lum = skimage.morphology.erosion((skimage.morphology.dilation(Lum, skimage.morphology.disk(k))), skimage.morphology.disk(k))\n",
    "    Lum = skimage.morphology.label(Lum); \n",
    "    Lum[np.nonzero(Lum >1)] = 0\n",
    "    Lum = morph.binary_fill_holes(Lum)\n",
    "    Med = skimage.morphology.erosion((skimage.morphology.dilation(Med, skimage.morphology.disk(k))), skimage.morphology.disk(k))\n",
    "    Med = morph.binary_fill_holes(Med)\n",
    "    \n",
    "    Media = Med ^ Lum\n",
    "    Org = np.empty((szCart1, szCart2, 3), dtype=np.uint8)\n",
    "    Org[:, :, 0] = Org[:, :, 1] = Org[:, :, 2] = OrgImage\n",
    "    Org[Media == 1,0] = 255; Org[Media == 1,1] = 0; Org[Media == 1,2] = 0;\n",
    "    Vol[:,:,:,fileIdx-1] = Org    \n",
    "    print(\"Time taken for the execution of the code:  %s seconds\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = dicom.read_file('data/day1/tut02/testVolume.dcm')  # loading a dummy .dcm file for meta data information\n",
    "\n",
    "Vol = Vol.transpose(3,0,1,2)\n",
    "A.PixelData = np.uint8(Vol).tostring()\n",
    "A.NumberofFrames = str(np.shape(Vol)[0])\n",
    "A.Rows = np.shape(Vol)[1]\n",
    "A.Columns = np.shape(Vol)[2]\n",
    "A.SamplesperPixel = 3\n",
    "A.PhotometricInterpretation = 'RGB'\n",
    "# Saving volume as dicom file\n",
    "A.save_as('data/day1/tut02/IVUSVolume.dcm') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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
}