{"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":{"collapsed":true,"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\nimport pandas as pd\nmontage3d = lambda x, **k: montage2d(np.stack([montage2d(y, **k) for y in x],0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"287d98f5de5a3359ca9719c7caf6839240e81f07","collapsed":true},"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')\ntarget_df = pd.read_csv('../input/train.csv')\ntrain_targets = {int(k['Id']): k for k in target_df.T.to_dict().values()}\nZ_DIM, T_DIM, X_DIM, Y_DIM = 10, 20, 128, 128\ntarget_df[['Systole', 'Diastole']].plot.hist(alpha = 0.5)\ntarget_df.sample(5)","execution_count":null,"outputs":[]},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"e2cf6bd7b43dddc8625137d2665d8e39b4240f58"},"cell_type":"code","source":"from scipy.ndimage import zoom\nrezoom = lambda in_data: zoom(in_data.images, [Z_DIM/in_data.images.shape[0], T_DIM/in_data.images.shape[1], \n                                               X_DIM/in_data.images.shape[2], Y_DIM/in_data.images.shape[3]])\nclass Dataset(object):\n    dataset_count = 0\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","collapsed":true},"cell_type":"code","source":"base_path = os.path.join(ALL_DATA_DIR, '140')\ntdata = Dataset(base_path,'140')\ntdata.load()\nprint(tdata.images.shape)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"39e8ae1213bbd5e2c93fc4ef14a804f99f076533","collapsed":true},"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(montage2d(tdata.images[0]), cmap = 'bone')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2cc286e0d4328fbacebd453d003bf0c7a12260fc","collapsed":true},"cell_type":"code","source":"image_stack = rezoom(tdata)","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        z_time = zoom(cur_data.time, 1.0*T_DIM/len(cur_data.time))\n        return [in_path, z_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(d.shape, 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":"df264089b6221240fb940653432dc22a22a02128"},"cell_type":"markdown","source":"# Big Loading\nHere we process many cases using dask and then export them. The limit for the export is the size of the final output"},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"86a738e284f690a020960adc0cfca81313952c61"},"cell_type":"code","source":"path_bag = dbag.from_sequence(np.random.choice(all_series, 250))\nimage_bag = path_bag.map(read_and_process)","execution_count":null,"outputs":[]},{"metadata":{"scrolled":false,"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.stack([x[-1] for x in all_img_data if x is not None],0).swapaxes(1, 2) # t, z, x, y\nprint(im_stack.shape, im_stack.dtype)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"61b7ee587e8e42f99218c2ee034863ccaa355028","collapsed":true},"cell_type":"code","source":"am_stack = np.stack([x[2] 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.stack([os.path.basename(x[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.stack([x[1] 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":"markdown","source":"# Export Volume Cubes as HDF5\nHere we export the volume cubes as hdf5 datasets"},{"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 = 8)\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":"823054c6e1b213cb0409895fa274c2d612111090"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}