{ "cells": [ { "cell_type": "markdown", "metadata": { "_cell_guid": "19c13585-0c9e-d802-297f-e8402aa52b8a" }, "source": [ "# Crop, Save and View Nodules in 3D\n", "\n", "Final results:\n", "\n", "![enter image description here][1]\n", "\n", "This kernel shows you how to crop the nodule **with given coordinates** in patient's CT scan, how to save it in .npy and .mhd/.raw file, and how to view it in **3D**. About half of the code is learned from other kernels. (Tutorial: U-Net Segmentation Approach to Cancer Diagnosis by [Jonathan Mulholland and Aaron Sander, Booz Allen Hamilton][2], Full Preprocessing Tutorial by [Guido Zuidhof][3], Candidate Generation and LUNA16 preprocessing by [ArnavJain][4]) The most fun part of this code is also learned from the Internet. I wanted a better way to visualize the region of interest my code generated. Then I searched online and found a method to save .mhd/.raw file posted by [Price Jackson][5], with [source code][6]. The part I found useful was built on the MIT licensed work by [Bing Jian and Baba C. Vemuri][7]. I happened to know [Fiji][8] was a good tool to view volume stacks in 3D, with original intensity values. After integrating them together, I find visualize and check the results is not that much of pain. Actually, it comes with some fun and may give you some insights on the journey. This incites me to open a kernel here, though most of the code are from others.\n", "\n", "Running the code below, with given CT scan and given coordinates, will give you a cropped nodule in [19, 19, 19] dimensional numpy array, with spacing [1, 1, 1]mm. I will use LUNA16 as input data in the code. Because it comes with some annotated nodules. So I cannot run it here but I put the psedo-output in Markdown cells. The code has been tested in Python 3.5.\n", "\n", "\n", " [1]: https://i.gyazo.com/24c585876be7e54a2fc20a40fbf3b2e9.gif\n", " [2]: https://www.kaggle.com/c/data-science-bowl-2017#tutorial\n", " [3]: https://www.kaggle.com/gzuidhof/data-science-bowl-2017/full-preprocessing-tutorial\n", " [4]: https://www.kaggle.com/arnavkj95/data-science-bowl-2017/candidate-generation-and-luna16-preprocessing\n", " [5]: https://sites.google.com/site/pjmedphys/tutorials/medical-images-in-python\n", " [6]: https://sites.google.com/site/pjmedphys/scripts\n", " [7]: https://code.google.com/archive/p/diffusion-mri/\n", " [8]: https://imagej.net/Fiji/Downloads" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "1296685b-2b9b-033c-4f90-c1f0d2e684a4" }, "source": [ "### Read annotation data and define preprocessing method\n", "\n", "Please change the input path to fit your environment" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "3247f6aa-5c88-ff88-2798-b33c0371e6e1" }, "outputs": [], "source": [ "import SimpleITK as sitk\n", "import numpy as np\n", "\n", "from glob import glob\n", "import pandas as pd\n", "import scipy.ndimage\n", "\n", "import mhd_utils_3d\n", "\n", "\n", "## Read annotation data and filter those without images\n", "# Learned from Jonathan Mulholland and Aaron Sander, Booz Allen Hamilton\n", "# https://www.kaggle.com/c/data-science-bowl-2017#tutorial\n", "\n", "# Set input path\n", "# Change to fit your environment\n", "luna_path = './LUNA16/'\n", "luna_subset_path = luna_path + 'subset0_samples/'\n", "file_list = glob(luna_subset_path + \"*.mhd\")\n", "\n", "df_node = pd.read_csv(luna_path+'annotations.csv')\n", "\n", "def get_filename(file_list, case):\n", " for f in file_list:\n", " if case in f:\n", " return(f)\n", "\n", "# map file full path to each record \n", "df_node['file'] = df_node['seriesuid'].map(lambda file_name: get_filename(file_list, file_name))\n", "df_node = df_node.dropna()\n", "\n", "## Define resample method to make images isomorphic, default spacing is [1, 1, 1]mm\n", "# Learned from Guido Zuidhof\n", "# https://www.kaggle.com/gzuidhof/data-science-bowl-2017/full-preprocessing-tutorial\n", "def resample(image, old_spacing, new_spacing=[1, 1, 1]):\n", " \n", " resize_factor = old_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 = old_spacing / real_resize_factor\n", "\n", " image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode = 'nearest')\n", " \n", " return image, new_spacing\n" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "f729e436-e441-8f9e-1cbe-5a1595496d20" }, "source": [ "### Define methods to save data\n", "\n", "I defined a method called \"save_nodule\" in this part. And it is called by:\n", "\n", " save_nodule(nodule_crop, name_index)\n", " # nodule_crop is a 3 dimensional numpy array, name_index is the name of the file. I use the index in annotation.csv.\n", "\n", "If your time is limited, just try to call the method with your data. If it works, then you can focus on building the model. It will save three files: **.npy**, **.mhd**, **.raw**." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "d6d29403-802b-182b-a4e1-8aa556837f7d" }, "outputs": [], "source": [ "#!/usr/bin/env python\n", "#coding=utf-8\n", "\n", "#======================================================================\n", "#Program: Diffusion Weighted MRI Reconstruction\n", "#Link: https://code.google.com/archive/p/diffusion-mri\n", "#Module: $RCSfile: mhd_utils.py,v $\n", "#Language: Python\n", "#Author: $Author: bjian $\n", "#Date: $Date: 2008/10/27 05:55:55 $\n", "#Version: \n", "# $Revision: 1.1 by PJackson 2013/06/06 $\n", "# Modification: Adapted to 3D\n", "# Link: https://sites.google.com/site/pjmedphys/tutorials/medical-images-in-python\n", "# \n", "# $Revision: 2 by RodenLuo 2017/03/12 $\n", "# Modication: Adapted to LUNA2016 data set for DSB2017\n", "# Link: \n", "#======================================================================\n", "\n", "import os\n", "import numpy\n", "import array\n", "\n", "def write_meta_header(filename, meta_dict):\n", " header = ''\n", " # do not use tags = meta_dict.keys() because the order of tags matters\n", " tags = ['ObjectType','NDims','BinaryData',\n", " 'BinaryDataByteOrderMSB','CompressedData','CompressedDataSize',\n", " 'TransformMatrix','Offset','CenterOfRotation',\n", " 'AnatomicalOrientation',\n", " 'ElementSpacing',\n", " 'DimSize',\n", " 'ElementType',\n", " 'ElementDataFile',\n", " 'Comment','SeriesDescription','AcquisitionDate','AcquisitionTime','StudyDate','StudyTime']\n", " for tag in tags:\n", " if tag in meta_dict.keys():\n", " header += '%s = %s\\n'%(tag,meta_dict[tag])\n", " f = open(filename,'w')\n", " f.write(header)\n", " f.close()\n", " \n", "def dump_raw_data(filename, data):\n", " \"\"\" Write the data into a raw format file. Big endian is always used. \"\"\"\n", " #Begin 3D fix\n", " data=data.reshape([data.shape[0],data.shape[1]*data.shape[2]])\n", " #End 3D fix\n", " rawfile = open(filename,'wb')\n", " a = array.array('f')\n", " for o in data:\n", " a.fromlist(list(o))\n", " #if is_little_endian():\n", " # a.byteswap()\n", " a.tofile(rawfile)\n", " rawfile.close()\n", " \n", "def write_mhd_file(mhdfile, data, dsize):\n", " assert(mhdfile[-4:]=='.mhd')\n", " meta_dict = {}\n", " meta_dict['ObjectType'] = 'Image'\n", " meta_dict['BinaryData'] = 'True'\n", " meta_dict['BinaryDataByteOrderMSB'] = 'False'\n", " meta_dict['ElementType'] = 'MET_FLOAT'\n", " meta_dict['NDims'] = str(len(dsize))\n", " meta_dict['DimSize'] = ' '.join([str(i) for i in dsize])\n", " meta_dict['ElementDataFile'] = os.path.split(mhdfile)[1].replace('.mhd','.raw')\n", " write_meta_header(mhdfile, meta_dict)\n", "\n", " pwd = os.path.split(mhdfile)[0]\n", " if pwd:\n", " data_file = pwd +'/' + meta_dict['ElementDataFile']\n", " else:\n", " data_file = meta_dict['ElementDataFile']\n", "\n", " dump_raw_data(data_file, data)\n", " \n", "def save_nodule(nodule_crop, name_index):\n", " np.save(str(name_index) + '.npy', nodule_crop)\n", " write_mhd_file(str(name_index) + '.mhd', nodule_crop, nodule_crop.shape[::-1])" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "00d9f771-d821-a33f-48e3-8520c1f4b831" }, "source": [ "### Read CT scan data, process and save " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "9d1e56c7-9814-1c74-4b27-f97ce8c64dda" }, "outputs": [], "source": [ "## Collect patients with nodule and crop the nodule\n", "# In this code snippet, the cropped nodule is a [19, 19, 19] volume with [1, 1, 1]mm spacing.\n", "# Learned from Jonathan Mulholland and Aaron Sander, Booz Allen Hamilton\n", "# https://www.kaggle.com/c/data-science-bowl-2017#tutorial\n", "\n", "# Change the number in the next line to process more\n", "for patient in file_list[:1]:\n", " print(patient)\n", " \n", " # Check whether this patient has nodule or not\n", " if patient not in df_node.file.values:\n", " print('Patient ' + patient + 'Not exist!')\n", " continue\n", " patient_nodules = df_node[df_node.file == patient]\n", " \n", " full_image_info = sitk.ReadImage(patient)\n", " full_scan = sitk.GetArrayFromImage(full_image_info)\n", " \n", " origin = np.array(full_image_info.GetOrigin())[::-1] # get [z, y, x] origin\n", " old_spacing = np.array(full_image_info.GetSpacing())[::-1] # get [z, y, x] spacing\n", " \n", " image, new_spacing = resample(full_scan, old_spacing)\n", " \n", " print('Resample Done')\n", " \n", "\n", " for index, nodule in patient_nodules.iterrows():\n", " nodule_center = np.array([nodule.coordZ, nodule.coordY, nodule.coordX]) \n", " # Attention: Z, Y, X\n", "\n", " v_center = np.rint( (nodule_center - origin) / new_spacing )\n", " v_center = np.array(v_center, dtype=int)\n", "\n", "# print(v_center)\n", " window_size = 9 # This will give you the volume length = 9 + 1 + 9 = 19\n", " # Why the magic number 19, I found that in \"LUNA16/annotations.csv\", \n", " # the 95th percentile of the nodules' diameter is about 19.\n", " # This is kind of a hyperparameter, will affect your final score.\n", " # Change it if you want.\n", " zyx_1 = v_center - window_size # Attention: Z, Y, X\n", " zyx_2 = v_center + window_size + 1\n", "\n", "# print('Crop range: ')\n", "# print(zyx_1)\n", "# print(zyx_2)\n", "\n", " # This will give you a [19, 19, 19] volume\n", " img_crop = image[ zyx_1[0]:zyx_2[0], zyx_1[1]:zyx_2[1], zyx_1[2]:zyx_2[2] ]\n", " \n", " # save the nodule \n", " save_nodule(img_crop, index)\n", " \n", " print('Done for this patient!\\n\\n')\n", "print('Done for all!')" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "b293cb60-905c-f237-39cf-15b13f4cc89f" }, "source": [ "Sample output:\n", "\n", " ./LUNA16/subset0_samples/1.3.6.1.4.1.14519.5.2.1.6279.6001.109002525524522225658609808059.mhd\n", " Resample Done\n", " Done for this patient!\n", "\n", "\n", " Done for all!" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "acedd386-8653-cc68-ab90-eda3fc710f03" }, "source": [ "### Plot in 2D" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "_cell_guid": "a371dea7-6618-90c5-75e6-8953a70dc392" }, "outputs": [], "source": [ "## Plot volume in 2D\n", "\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "def plot_nodule(nodule_crop):\n", " \n", " # Learned from ArnavJain\n", " # https://www.kaggle.com/arnavkj95/data-science-bowl-2017/candidate-generation-and-luna16-preprocessing\n", " f, plots = plt.subplots(int(nodule_crop.shape[0]/4)+1, 4, figsize=(10, 10))\n", " \n", " for z_ in range(nodule_crop.shape[0]): \n", " plots[int(z_/4), z_ % 4].imshow(nodule_crop[z_,:,:])\n", " \n", " # The last subplot has no image because there are only 19 images.\n", " plt.show()\n", " \n", "# Plot one example\n", "img_crop = np.load('25.npy')\n", "plot_nodule(img_crop)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "e2e7aa3e-eac7-50e2-8d29-92a0684a58f7" }, "source": [ "Sample output:\n", "\n", "[![https://gyazo.com/31240b93448349629652fe56cfd3f48f](https://i.gyazo.com/31240b93448349629652fe56cfd3f48f.png)](https://gyazo.com/31240b93448349629652fe56cfd3f48f)" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "bd702ad1-e400-83fa-937d-69149d370bbf" }, "source": [ "### Plot in 3D by Fiji" ] }, { "cell_type": "markdown", "metadata": { "_cell_guid": "80b960ed-c7a0-3e52-3d40-34f8bfa5d2a3" }, "source": [ "In your output path (default is the same as your working directory), there should be some .mhd/.raw files. Follow the steps below to open it in 3D viewer. (I also made a video (18min) and posted in [this discussion thread][1]. You may want to watch it if you like video tutorials. Excuse me that I'm not fluent in English. So the video is recorded kind of slow, you may want to watch it with 1.5 or 2 times speed.)\n", "\n", "### Steps to open nodules in 3D Viewer \n", "1. Download [Fiji][2]\n", "2. Drag one .mhd file to Fiji **status bar**\n", "\n", " ![enter image description here][3]\n", "\n", "3. Click on the new image window. Press \"control +\" to zoom in, \"control -\" to zoom out.\n", "4. In the **Menubar**, Click \"Image > Stacks > Orthogonal Views\" to see it. Scroll to go through all slices. **Notice that when your cursor moves around in the image window, the Fiji Status Bar shows you the XY coordinates and the value, this value is the same as that in Python numpy array**\n", "5. In the **Menubar**, Click \"Plugins > 3D viewer\", in the \"Add ...\" window, change \"Resampling factor\" to 1, click \"OK\", click \"OK\" to convert to 8-bit.\n", "6. Click on the \"ImageJ 3D Viewer\", In the **Menubar**, click \"Edit > Adjust threshold\". Drag threshold bar to around 150.\n", "7. **Drag** the object to rotate it. Hold \"shift\" and **drag** to move it. Scroll to zoom. \n", "8. In the **Menubar**, click \"View > Start animation\" to activate it.\n", "9. Click on the image stack, in the **Menubar**, click \"Image > Adjust > threshold\", check \"Dark background\", move the first threshold bar to around -400, click \"Apply > OK > Yes\". \n", "10. Either use step 5 to open another 3D viewer, or click on the current \"ImageJ 3D Viewer\", in the **Menubar** click \"add > from image\" to add another object. Click on the object and \"shift\" drag to move it. Thresholding on the raw image then creating 3D object gives you different rendering. \n", "\n", "**Sample visualizations:**\n", "\n", "25 in LUNA16/annotation.csv\n", "\n", "![enter image description here][4]\n", "\n", "26 in LUNA16/annotation.csv\n", "\n", "![enter image description here][5]\n", "\n", "\n", "\n", "\n", "**Please feel free to comment if you have questions or suggestions. Please upvote the above mentioned kernels if you find they are helpful. Please upvote this kernel if it helps you on the journey.**\n", "\n", "\n", "\n", "## Updates:\n", "\n", "Changed \"nodule_crop.shape\" to \"nodule_crop.shape[::-1]\" in method \"save_nodule(nodule_crop, name_index)\"\n", "\n", " [1]: https://www.kaggle.com/c/data-science-bowl-2017/discussion/28502\n", " [2]: https://imagej.net/Fiji/Downloads\n", " [3]: https://imagej.net/_images/6/67/Fiji-main-window.jpg\n", " [4]: https://i.gyazo.com/d9eea0182a10af8b8af8f6cf88c3a0e7.gif\n", " [5]: https://i.gyazo.com/3c69952f966704b4055e25bc07495748.gif" ] } ], "metadata": { "_change_revision": 0, "_is_fork": false, "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.0" } }, "nbformat": 4, "nbformat_minor": 0 }