{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport glob\nfrom tqdm import tqdm\nfrom random import randint\n\nimport numpy as np\nimport pydicom\n\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport matplotlib.animation as anim\nimport matplotlib.patches as mpatches\nimport matplotlib.gridspec as gridspec\n\n\nimport imageio\nfrom skimage.transform import resize\nfrom skimage.util import montage\n\nfrom IPython.display import Image as show_gif\n\nimport warnings\nwarnings.simplefilter(\"ignore\")","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:50:13.790379Z","iopub.execute_input":"2021-09-13T12:50:13.790819Z","iopub.status.idle":"2021-09-13T12:50:15.422584Z","shell.execute_reply.started":"2021-09-13T12:50:13.790733Z","shell.execute_reply":"2021-09-13T12:50:15.421376Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_id = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000'\npath_x = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/T2w/*.dcm'\nstart = len('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/T2w/Image-')\nend = len('.dcm')\npath_to_slices = sorted(glob.glob(path_x), key= lambda x: int(x[start:-end]))\n","metadata":{"execution":{"iopub.status.busy":"2021-09-13T12:50:15.425169Z","iopub.execute_input":"2021-09-13T12:50:15.425647Z","iopub.status.idle":"2021-09-13T12:50:15.514319Z","shell.execute_reply.started":"2021-09-13T12:50:15.425612Z","shell.execute_reply":"2021-09-13T12:50:15.513294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The easiest way to render a 3d image is to draw its slices one by one in 2d , for easy viewing we can make a GIF from each slice","metadata":{}},{"cell_type":"code","source":"class ImageToGIF:\n    \"\"\"Create GIF without saving image files.\"\"\"\n    def __init__(self,\n                 size=(500, 500), \n                 xy_text=(80, 30),\n                 dpi=100, \n                 cmap='CMRmap'):\n\n        self.fig = plt.figure()\n        self.fig.set_size_inches(size[0] / dpi, size[1] / dpi)\n        self.xy_text = xy_text\n        self.cmap = cmap\n        \n        self.ax = self.fig.add_axes([0, 0, 1, 1])\n        self.ax.set_xticks([])\n        self.ax.set_yticks([])\n        self.images = []\n \n    def add(self, image, label, with_mask=False):\n        plt.set_cmap(self.cmap)\n        plt_img = self.ax.imshow(image, animated=True)\n        plt_text = self.ax.text(*self.xy_text, label, color='red')\n        to_plot = [plt_img, plt_text]\n        self.images.append(to_plot)\n        plt.close()\n \n    def save(self, filename, fps):\n        animation = anim.ArtistAnimation(self.fig, self.images)\n        animation.save(filename, writer='imagemagick', fps=fps)\n        \n\nsample_data_gif = ImageToGIF()\nlabel = sample_id.replace('/', '.').split('.')[-2]\nfilename = f'{label}_3d_2d.gif'\n\nfor i in range(len(path_to_slices)):\n    image = pydicom.read_file(path_to_slices[i]).pixel_array\n    #mask = np.clip(np.rot90(sample_mask[i]), 0, 1)\n    sample_data_gif.add(image, label=f'{label}_{str(i)}')\n\n    \nsample_data_gif.save(filename, fps=15)\nshow_gif(filename, format='png')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-09-13T12:50:15.516678Z","iopub.execute_input":"2021-09-13T12:50:15.516977Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## From pixels to voxels","metadata":{}},{"cell_type":"code","source":"class Image3dToGIF3d:\n    \"\"\"\n    Displaying 3D images in 3d axes.\n    Parameters:\n        img_dim: shape of cube for resizing.\n        figsize: figure size for plotting in inches.\n    Step by step explanation - https://terbium.io/2017/12/matplotlib-3d/\n    \"\"\"\n    def __init__(self, \n                 img_dim: tuple = (55, 55, 55),\n                 figsize: tuple = (15, 10),\n                 binary: bool = False,\n                 normalizing: bool = True,\n                ):\n        \"\"\"Initialization.\"\"\"\n        self.img_dim = img_dim\n        print(img_dim)\n        self.figsize = figsize\n        self.binary = binary\n        self.normalizing = normalizing\n\n    def _explode(self, data: np.ndarray):\n        \"\"\"\n        Takes: array and return an array twice as large in each dimension,\n        with an extra space between each voxel.\n        \"\"\"\n        shape_arr = np.array(data.shape)\n        size = shape_arr[:3] * 2 - 1\n        exploded = np.zeros(np.concatenate([size, shape_arr[3:]]),\n                            dtype=data.dtype)\n        exploded[::2, ::2, ::2] = data\n        return exploded\n\n    def _expand_coordinates(self, indices: np.ndarray):\n        \"\"\" \n        Parameters:\n            indices: coordinats of array with only original values\n            (before explode transformaion)\n        \n        Returns:\n        The arrays of values each dimensions (x y z) as arguments needed \n        for the plt.figure.voxels functionto extend coordinates\n        for the rendering only colored voxels\"\"\"\n        x, y, z = indices\n        x[1::2, :, :] += 1\n        y[:, 1::2, :] += 1\n        z[:, :, 1::2] += 1\n        return x, y, z\n    \n    def _normalize(self, arr: np.ndarray):\n        \"\"\"Normilize image value between 0 and 1.\"\"\"\n        arr_min = np.min(arr)\n        return (arr - arr_min) / (np.max(arr) - arr_min)\n\n    \n    def _scale_by(self, arr: np.ndarray, factor: int = 2):\n        \"\"\"\n        Scale 3d Image to factor (guesstimated transformation).\n        Parameters:\n            arr: 3d image for scalling.\n            factor: factor for scalling.\n        \"\"\"\n        mean = np.mean(arr)\n        return (arr - mean) * factor + mean\n    \n    def get_transformed_data(self, data: np.ndarray):\n        \"\"\"Data transformation: normalization, scaling, resizing.\"\"\"\n        if self.binary:\n            resized_data = resize(data, self.img_dim, preserve_range=True)\n            return np.clip(resized_data.astype(np.uint8), 0, 1).astype(np.float32)\n            \n        norm_data = np.clip(self._normalize(data)-0.1, 0, 1) ** 0.4\n        scaled_data = np.clip(self._scale_by(norm_data) - 0.1, 0, 1)\n        resized_data = resize(scaled_data, self.img_dim, preserve_range=True)\n        \n        return resized_data\n    \n    def plot_cube(self,\n                  cube,\n                  title: str = '', \n                  init_angle: int = 0,\n                  make_gif: bool = False,\n                  path_to_save: str = 'filename.gif'\n                 ):\n        \"\"\"\n        Plot 3d data.\n        -> Take array \n        -> return an array twice as large in each dimension,\n        (with an extra space between each voxel)\n        -> expand coordinates of each voxel for core rendering (to remove gaps)\n        because additional fake voxels have been added\n        -> set each voxel’s transparency equal to its value.       \n        Parameters:\n            cube: 3d data\n            title: title for figure.\n            init_angle: angle for image plot (from 0-360).\n            make_gif: if True create gif from every 5th frames from 3d image plot.\n            path_to_save: path to save GIF file.\n            \"\"\"\n\n \n        if self.normalizing:\n            cube = self._normalize(cube)\n            \n        facecolors = cm.gist_stern(cube)          \n        facecolors[:,:,:,-1] = cube\n        facecolors = self._explode(facecolors)\n\n        filled = facecolors[:,:,:,-1] != 0\n        x, y, z = self._expand_coordinates(np.indices(np.array(filled.shape) + 1))\n\n        with plt.style.context(\"dark_background\"):\n\n            fig = plt.figure(figsize=self.figsize)\n            ax = fig.gca(projection='3d')\n\n            ax.view_init(30, init_angle)\n            ax.set_xlim(right = self.img_dim[0] * 2)\n            ax.set_ylim(top = self.img_dim[1] * 2)\n            ax.set_zlim(top = self.img_dim[2] * 2)\n            ax.set_title(title, fontsize=18, y=1.05)\n\n            ax.voxels(x, y, z, filled, facecolors=facecolors, shade=False)\n\n            if make_gif:\n                images = []\n                for angle in tqdm(range(0, 360, 5)):\n                    ax.view_init(30, angle)\n                    fname = str(angle) + '.png'\n\n                    plt.savefig(fname, dpi=120, format='png', bbox_inches='tight')\n                    images.append(imageio.imread(fname))\n                    #os.remove(fname)\n                imageio.mimsave(path_to_save, images)\n                plt.close()\n\n            else:\n                plt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tensor = np.zeros((512, 512, len(path_to_slices)))\nfor i in range(len(path_to_slices)):\n    image = pydicom.read_file(path_to_slices[i]).pixel_array\n    tensor[:,:,i] = image","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"title = sample_id.replace(\".\", \"/\").split(\"/\")[-1]\nfilename = title+\"_3d.gif\"\n\ndata_to_3dgif = Image3dToGIF3d(img_dim = (120, 120, 78))\ntransformed_data = data_to_3dgif.get_transformed_data(np.moveaxis(np.flipud(tensor), [0, 1, 2], [-1, -2, -3]))\ndata_to_3dgif.plot_cube(\n    transformed_data[:77, :100, :55],\n    title=title,\n    make_gif=True,\n    path_to_save=filename\n)\n\nshow_gif(filename, format='png')","metadata":{"_kg_hide-input":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can also compare images\nwith T1-weighted, T1CE-weighted, T2-weighted and FLAIR-weighted, and notice that their range of values is different","metadata":{}},{"cell_type":"code","source":"NUM = '00012'\n\npath_flair = f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/FLAIR/*.dcm'\nstart_flair = len(f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/FLAIR/Image-')\nend_flair = len('.dcm')\npath_to_slices_flair = sorted(glob.glob(path_flair), key= lambda x: int(x[start_flair:-end_flair]))\n\npath_t1w = f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T1w/*.dcm'\nstart_t1w = len(f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T1w/Image-')\nend_t1w = len('.dcm')\npath_to_slices_t1w = sorted(glob.glob(path_t1w), key= lambda x: int(x[start_t1w:-end_t1w]))\n\npath_t1wce = f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T1wCE/*.dcm'\nstart_t1wce = len(f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T1wCE/Image-')\nend_t1wce = len('.dcm')\npath_to_slices_t1wce = sorted(glob.glob(path_t1wce), key= lambda x: int(x[start_t1wce:-end_t1wce]))\n\npath_t2w = f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T2w/*.dcm'\nstart_t2w = len(f'../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/{NUM}/T2w/Image-')\nend_t2w = len('.dcm')\npath_to_slices_t2w = sorted(glob.glob(path_t2w), key= lambda x: int(x[start_t2w:-end_t2w]))\n\n\ndef slices_dcm_to_tensor(pth: str) -> np.ndarray:\n    tensor = np.zeros((512, 512,len(pth)))\n    for i in range(len(pth)):\n        image = pydicom.read_file(pth[i]).pixel_array\n        tensor[:,:,i] = image\n    return tensor","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flair_arr = slices_dcm_to_tensor(path_to_slices_flair)\nt1w_arr = slices_dcm_to_tensor(path_to_slices_t1w)\nt1wce_arr = slices_dcm_to_tensor(path_to_slices_t1wce)\nt2w_arr = slices_dcm_to_tensor(path_to_slices_t2w)\n\nprint(flair_arr.shape, t1w_arr.shape, t1wce_arr.shape, t2w_arr.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If we look at the color scales for the random slices, we can see how many points reach the high or low values (greater than, say, 4000, or less than 500).","metadata":{}},{"cell_type":"code","source":"fig = plt.figure(figsize=(20, 10))\n\ngs = gridspec.GridSpec(nrows=2, ncols=4, height_ratios=[1, 1.5])\n\n#  Varying density along a streamline\nax0 = fig.add_subplot(gs[0, 0])\nflair = ax0.imshow(flair_arr[:,:,150], cmap='bone')\nax0.set_title(\"FLAIR\", fontsize=18, weight='bold', y=-0.2)\nfig.colorbar(flair)\n\n#  Varying density along a streamline\nax1 = fig.add_subplot(gs[0, 1])\nt1 = ax1.imshow(t1w_arr[:,:,15], cmap='bone')\nax1.set_title(\"T1\", fontsize=18, weight='bold', y=-0.2)\nfig.colorbar(t1)\n\n#  Varying density along a streamline\nax2 = fig.add_subplot(gs[0, 2])\nt1ce = ax2.imshow(t1wce_arr[:,:,160], cmap='bone')\nax2.set_title(\"T1 contrast\", fontsize=18, weight='bold', y=-0.2)\nfig.colorbar(t1ce)\n\n#  Varying density along a streamline\nax3 = fig.add_subplot(gs[0, 3])\nt2 = ax3.imshow(t2w_arr[:,:,175], cmap='bone')\nax3.set_title(\"T2\", fontsize=18, weight='bold', y=-0.2)\nfig.colorbar(t2)\n\nplt.show();","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]}]}