{"cells":[{"metadata":{"_uuid":"d28fa6267504a734ee96e9f9bd7859c1e44d2f65"},"cell_type":"markdown","source":"# Overview\nA simple notebook to read in a few example datasets and package the rest of the datasets together for making training models easier (HDF5 instead of dicom mess). \n\nThe code uses code borrows heavily from the tutorial provided by Booz (https://github.com/booz-allen-hamilton/DSB2) in order to preprocess the data. You can make new kernels by having this dataset as a starting point"},{"metadata":{"trusted":true,"_uuid":"a4cd5e1cbc5edbcc0dfc4c379206ca75542ac947"},"cell_type":"code","source":"%matplotlib inline\nimport cv2\nimport numpy as np\nimport pydicom as dicom\nimport json\nimport os\nimport shutil\nimport sys\nimport random\nfrom matplotlib import image\nimport matplotlib.pyplot as plt\nimport re\nfrom skimage.util.montage import montage2d\nmontage3d = lambda x, **k: montage2d(np.stack([montage2d(y, **k) for y in x],0))\nimport h5py\n\nimport time\nstart= time.time()\nimport numpy as np \nimport pandas as pd \nimport pydicom\n# import dicom\nimport os\nimport scipy.ndimage\nimport matplotlib.pyplot as plt\n\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\nprint(\"done\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan_test(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n    for s in slices:\n        s.SliceThickness = slice_thickness\n    return slices\n\nprint(\"done\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"ALL_DATA_DIR =  os.path.join('..', 'input', 'train', 'train','140','study','2ch_6')\nprint\n# tData = Dataset(base_path,'140')\nprint(os.listdir('/input/second-annual-data-science-bowl/train/train/1'))\n\n# ../input/second-annual-data-science-bowl/train/train/140/study/2ch_6/'IM-0699-0001.dcm'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"# load_scan_test('../input/second-annual-data-science-bowl/test/test/1000/study/2ch_59')\nALL_DATA_DIR =  os.path.join('..', 'input', 'train', 'train','140','study','2ch_6')\nslices = [pydicom.read_file(ALL_DATA_DIR + '/' + s) for s in os.listdir(ALL_DATA_DIR)]\nprint(len(slices))\nslices.sort(key = lambda x: int(x.InstanceNumber))\ntry:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\nexcept:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\nfor s in slices:\n        s.SliceThickness = slice_thickness\nprint(slices)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"# load_scan_test('../input/second-annual-data-science-bowl/test/test/1000/study/2ch_59')\nslices = [pydicom.read_file(ALL_DATA_DIR + '/' + 'IM-0699-0001.dcm')]\nslices.sort(key = lambda x: int(x.InstanceNumber))\n# try:\n#         slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n# except:\n#         slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \ntry:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\nexcept:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \nfor s in slices:\n        s.SliceThickness = slice_thickness\nprint(slices)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load_scan_test('../input/second-annual-data-science-bowl/test/test/1000/study/2ch_59')\n# ALL_DATA_DIR =  os.path.join('..', 'input', 'train', 'train','140','study','2ch_6')\nALL_DATA_DIR = '/kaggle/input/second-annual-data-science-bowl/train/train/140/study/2ch6'\nprint(ALL_DATA_DIR)\npixel_array = [pydicom.read_file(ALL_DATA_DIR + '/' + 'IM-0699-0001.dcm').pixel_array]\nprint(pixel_array)\nh5 = h5py.File(\"output1.h5\")\nh5.create_dataset(\"data\", data=pixel_array)\nh5.close()\n# print()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# os.chdir(os.path.join('..', 'input', 'train', 'train','140','study','2ch_6'))\nprint(os.listdir('/kaggle/input/second-annual-data-science-bowl/train/train/140/study/sax_36'))\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir('/kaggle/working'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nfile_name = '/kaggle/working/out_final.h5'\nhf = h5py.File(file_name)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nprint('Keys:', list(hf.keys()))\nprint('Attrs:', dict(hf.attrs))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"######### testing dicom to h5 conversion starts here","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"%matplotlib inline\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport pydicom\n# import dicom\nimport os\nimport scipy.ndimage\nimport matplotlib.pyplot as plt\n\nfrom skimage import measure, morphology\nfrom mpl_toolkits.mplot3d.art3d import Poly3DCollection\n\n# Some constants \nINPUT_FOLDER = '/kaggle/second-annual-data-science-bowl/train/train/140/study/2ch6'\npatients = os.listdir(INPUT_FOLDER)\npatients.sort()\nprint(\"done\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_scan(path):\n    slices = [pydicom.read_file(path + '/' + s) for s in os.listdir(path)]\n    slices.sort(key = lambda x: int(x.InstanceNumber))\n    \n    try:\n        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])\n    except:\n        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n        \n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# conversion of dicom to h5 ends here","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# load_scan_test('../input/second-annual-data-science-bowl/test/test/1000/study/2ch_59')\nALL_DATA_DIR =  os.path.join('..', 'input', 'train', 'train','140','study','2ch_6')\n\nprint(ALL_DATA_DIR)\npixel_array = [pydicom.read_file(ALL_DATA_DIR + '/' + s).pixel_array  for s in os.listdir(ALL_DATA_DIR)]\n# print(pixel_array)\n# h5 = h5py.File(\"out.h5\")\nhf=h5py.File(h5py.File(\"out.h5\"))\n# h5.create_dataset(\"data\", data=pixel_array)\n# h5.close()\n\nprint('Keys:', list(hf.keys()))\nprint('Attrs:', dict(hf.attrs))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# OUT_DIR =  os.path.join('..','output', 'kaggle', 'working')\nOUT_DIR =  os.path.join('./')\nprint(OUT_DIR)\nfile_name = OUT_DIR + '/' + 'out.h5' \n\n# hf = [h5py.File(file_name) for file_name in os.listdir(OUT_DIR)]\nout=h5py.File(file_name)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"287d98f5de5a3359ca9719c7caf6839240e81f07"},"cell_type":"code","source":"# number of bins to use in histogram for gaussian regression\nNUM_BINS = 100\n# number of standard deviatons past which we will consider a pixel an outlier\nSTD_MULTIPLIER = 2\n# number of points of our interpolated dataset to consider when searching for\n# a threshold value; the function by default is interpolated over 1000 points,\n# so 250 will look at the half of the points that is centered around the known\n# myocardium pixel\nTHRESHOLD_AREA = 250\n# number of pixels on the line within which to search for a connected component\n# in a thresholded image, increase this to look for components further away\nCOMPONENT_INDEX_TOLERANCE = 20\n# number of angles to search when looking for the correct orientation\nANGLE_SLICES = 36\nALL_DATA_DIR =  os.path.join('..', 'input', 'train', 'train')\nX_DIM, Y_DIM = 64, 64\nX_DIM, Y_DIM = 128, 128\nT_DIM = 30","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e2cf6bd7b43dddc8625137d2665d8e39b4240f58"},"cell_type":"code","source":"class Dataset(object):\n    dataset_count = 0\n\n    def __init__(self, directory, subdir):\n        # deal with any intervening directories\n        while True:\n            subdirs = next(os.walk(directory))[1]\n            if len(subdirs) == 1:\n                directory = os.path.join(directory, subdirs[0])\n            else:\n                break\n\n        slices = []\n        for s in subdirs:\n            m = re.match(\"sax_(\\d+)\", s)\n            if m is not None:\n                slices.append(int(m.group(1)))\n\n        slices_map = {}\n        first = True\n        times = []\n        for s in slices:\n            files = next(os.walk(os.path.join(directory, \"sax_%d\" % s)))[2]\n            offset = None\n\n            for f in files:\n                m = re.match(\"IM-(\\d{4,})-(\\d{4})\\.dcm\", f)\n                if m is not None:\n                    if first:\n                        times.append(int(m.group(2)))\n                    if offset is None:\n                        offset = int(m.group(1))\n\n            first = False\n            slices_map[s] = offset\n\n        self.directory = directory\n        self.time = sorted(times)\n        self.slices = sorted(slices)\n        self.slices_map = slices_map\n        self.name = subdir\n\n    def _filename(self, s, t):\n        return os.path.join(self.directory,\"sax_%d\" % s, \"IM-%04d-%04d.dcm\" % (self.slices_map[s], t))\n\n    def _read_dicom_image(self, filename):\n        d = dicom.read_file(filename)\n        img = d.pixel_array\n        return np.array(img)\n\n    def _read_all_dicom_images(self):\n        f1 = self._filename(self.slices[0], self.time[0])\n        d1 = dicom.read_file(f1)\n        (x, y) = d1.PixelSpacing\n        (x, y) = (float(x), float(y))\n        f2 = self._filename(self.slices[1], self.time[0])\n        d2 = dicom.read_file(f2)\n\n        # try a couple of things to measure distance between slices\n        try:\n            dist = np.abs(d2.SliceLocation - d1.SliceLocation)\n        except AttributeError:\n            try:\n                dist = d1.SliceThickness\n            except AttributeError:\n                dist = 8  # better than nothing...\n\n        self.images = np.array([[self._read_dicom_image(self._filename(d, i))\n                                 for i in self.time]\n                                for d in self.slices])\n        self.dist = dist\n        self.area_multiplier = x * y\n\n    def load(self):\n        self._read_all_dicom_images()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d3578526b2475f036e441874aabc9a6ded768a5f"},"cell_type":"markdown","source":"# Load a test patient\nHere we load patient 140 as an example"},{"metadata":{"trusted":true,"_uuid":"253fecfe97612ab143d6b2866beb86d928c973e1"},"cell_type":"code","source":"base_path = os.path.join(ALL_DATA_DIR, '140')\ntData = Dataset(base_path,'140')\ntData.load()\nprint(tData.slices)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"39e8ae1213bbd5e2c93fc4ef14a804f99f076533"},"cell_type":"code","source":"%matplotlib inline\nplt.imshow(tData.images[0,0,:,:], cmap = 'bone')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"ac7a1663d1670e47173f301fb8acdc073e5bcd2d","collapsed":true},"cell_type":"code","source":"fig, ax1 = plt.subplots(1,1, figsize = (20,20))\nax1.imshow(montage3d(tData.images), cmap = 'bone')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2cc286e0d4328fbacebd453d003bf0c7a12260fc"},"cell_type":"code","source":"from scipy.ndimage import zoom\nrezoom = lambda in_data: zoom(in_data.images, [1, \n                                               T_DIM/in_data.images.shape[1], \n                                               X_DIM/in_data.images.shape[2], \n                                               Y_DIM/in_data.images.shape[3]])\nimage_stack = rezoom(tData)\nprint(image_stack.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c2f7ab57adec22d0ebf224a445eca82b26cf423b","collapsed":true},"cell_type":"code","source":"fig, ax1 = plt.subplots(1,1, figsize = (20,20))\nax1.imshow(montage3d(image_stack), cmap = 'bone')","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"1d24816d483fc1820421db15e1b17d18d3fa5d86"},"cell_type":"code","source":"from glob import glob\nbase_path = os.path.join(ALL_DATA_DIR, '*')\nall_series = glob(base_path)\nfrom warnings import warn\ndef read_and_process(in_path):\n    try:\n        cur_data = Dataset(in_path,\n                           os.path.basename(in_path))\n        cur_data.load()\n        if cur_data.time is not None:\n            zoom_time = zoom(cur_data.time, [T_DIM/len(cur_data.time)])\n        else:\n            zoom_time = range(T_DIM)\n        return [in_path, zoom_time, cur_data.area_multiplier, rezoom(cur_data)]\n    except Exception as e:\n        warn('{}'.format(e), RuntimeWarning)\n        return None","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"692f4857f44e170742868f904a196cdcc0152526","collapsed":true},"cell_type":"code","source":"%%time\na,d,b,c = read_and_process(all_series[-100])\nprint(c.shape)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"08f7b3711a598aa53819e39d001205495f2e68c1"},"cell_type":"code","source":"import dask\nimport dask.diagnostics as diag\nfrom bokeh.io import output_notebook\nfrom bokeh.resources import CDN\nfrom dask import bag as dbag\nfrom multiprocessing.pool import ThreadPool","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7d4d2eebdd490ad9c60294fda31237a39fac8d4f"},"cell_type":"markdown","source":"# Final Processing\nHere we randomly select half of the patients for further processing (since @Kaggle only lets us output 1GB or so of data)"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"86a738e284f690a020960adc0cfca81313952c61"},"cell_type":"code","source":"np.random.seed(2018)\npath_bag = dbag.from_sequence(np.random.choice(all_series, 170))\nimage_bag = path_bag.map(read_and_process)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":true,"_uuid":"daa7410bc2d938899271c0da413aa4991589c415","collapsed":true},"cell_type":"code","source":"with diag.ProgressBar(), diag.Profiler() as prof, diag.ResourceProfiler(0.5) as rprof, dask.set_options(pool = ThreadPool(4)):\n    all_img_data = image_bag.compute()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"85978f9dd72db7f8d70fdb6168328b8ebba74f59","collapsed":true},"cell_type":"code","source":"output_notebook(CDN, hide_banner=True)\ndiag.visualize([prof, rprof])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9407edfc4024532b05d77e9cdd63ab350961cbb1","collapsed":true},"cell_type":"code","source":"im_stack = np.concatenate([x[-1] for x in all_img_data if x is not None],0)\nprint(im_stack.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"61b7ee587e8e42f99218c2ee034863ccaa355028","collapsed":true},"cell_type":"code","source":"am_stack = np.concatenate([ [x[2]]*x[-1].shape[0] for x in all_img_data if x is not None],0)\nprint(am_stack.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"d8ebe389daed355e11189c029808feac8cb644f5","collapsed":true},"cell_type":"code","source":"path_stack = np.concatenate([ [os.path.basename(x[0])]*x[-1].shape[0] for x in all_img_data if x is not None],0)\nprint(path_stack.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8c088f14abe256303f3ddf100ddc01751d30561d","collapsed":true},"cell_type":"code","source":"time_stack = np.concatenate([ [x[1]]*x[-1].shape[0] for x in all_img_data if x is not None],0)\nprint(time_stack.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e88b61f7e2b48eb45b9d144b8ab87f907678644a","collapsed":true},"cell_type":"code","source":"import pandas as pd\ntrain_targets = {int(k['Id']): k for k in pd.read_csv('../input/train.csv').T.to_dict().values()}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"290cc03aa1d7d03ec12cba902a3ea3d9cc09b7b0","collapsed":true},"cell_type":"code","source":"import h5py\nwith h5py.File('train_mri_{}_{}.h5'.format(X_DIM, Y_DIM), 'w') as w:\n    w.create_dataset('image', data = im_stack, compression = 9)\n    w.create_dataset('systole', data = [train_targets[int(c_id)]['Systole'] for c_id in path_stack])\n    w.create_dataset('diastole', data = [train_targets[int(c_id)]['Diastole'] for c_id in path_stack])\n    w.create_dataset('id', data = [int(c_id) for c_id in path_stack])\n    w.create_dataset('area_multiplier', data = am_stack)\n    w.create_dataset('time', data = time_stack)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8df91079402d488f674a8434ff99cff57b3ca3ee","collapsed":true},"cell_type":"code","source":"!ls -lh *","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a52b35a4c43cc44f3e78d6ddbf474a09c3467a84"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}