{"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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom matplotlib import pyplot as plt\nimport os\nimport cv2\nimport pydicom\nimport timeit\nimport seaborn as sns\nplt.style.use('seaborn-talk')\nimport glob\nfrom pathlib import Path\nfrom tqdm import tqdm\nfrom matplotlib import animation, rc\nimport joblib\nimport difflib\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-08-26T15:08:16.582472Z","iopub.execute_input":"2021-08-26T15:08:16.582924Z","iopub.status.idle":"2021-08-26T15:08:18.092986Z","shell.execute_reply.started":"2021-08-26T15:08:16.582832Z","shell.execute_reply":"2021-08-26T15:08:18.091805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def normalize_by_img_max(img):\n#     print(img.min(), img.max(), img.dtype)\n    img_max = max(1, img.max())\n    return (img / img_max * 255).astype('uint8')\n\ndef read_dcm(dcm_name, normalize=True):\n    np_img = pydicom.read_file(dcm_name).pixel_array\n    if normalize:\n        np_img = normalize_by_img_max(np_img)\n    return np_img\n\nread_dcm('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/FLAIR/Image-273.dcm', normalize=False)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.094841Z","iopub.execute_input":"2021-08-26T15:08:18.095282Z","iopub.status.idle":"2021-08-26T15:08:18.125306Z","shell.execute_reply.started":"2021-08-26T15:08:18.095235Z","shell.execute_reply":"2021-08-26T15:08:18.124130Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def compare_two_dicom_headers(file1, file2):\n    datasets = tuple([pydicom.dcmread(filename, force=True)\n                  for filename in (file1, file2)])\n\n    rep = []\n    for dataset in datasets:\n        lines = str(dataset).split(\"\\n\")\n        lines = [line + \"\\n\" for line in lines]  # add the newline to end\n        rep.append(lines)\n\n\n    diff = difflib.Differ()\n    for line in diff.compare(rep[0], rep[1]):\n        if line[0] != \"?\" and line[0] in ['-', '+']:\n            print(line)\n    \nfile_name1 = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/FLAIR/Image-273.dcm'\nfile_name2 = '../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/FLAIR/Image-272.dcm'\n\ncompare_two_dicom_headers(file_name1, file_name2)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.127554Z","iopub.execute_input":"2021-08-26T15:08:18.127982Z","iopub.status.idle":"2021-08-26T15:08:18.163954Z","shell.execute_reply.started":"2021-08-26T15:08:18.127938Z","shell.execute_reply":"2021-08-26T15:08:18.163044Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pydicom.read_file('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00031/T1w/Image-33.dcm')","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.165477Z","iopub.execute_input":"2021-08-26T15:08:18.165997Z","iopub.status.idle":"2021-08-26T15:08:18.184378Z","shell.execute_reply.started":"2021-08-26T15:08:18.165964Z","shell.execute_reply":"2021-08-26T15:08:18.183481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# <h1> Resample volume  1x1x1mm volumes","metadata":{}},{"cell_type":"code","source":"data_path = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00031/T1w\"\ng = glob.glob(data_path + '/*.dcm')\n\n# Print out the first 5 file names to verify we're in the right folder.\nprint (\"Total of %d DICOM images.\\nFirst 5 filenames:\" % len(g))\nprint(g[:5])","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.185724Z","iopub.execute_input":"2021-08-26T15:08:18.186303Z","iopub.status.idle":"2021-08-26T15:08:18.195397Z","shell.execute_reply.started":"2021-08-26T15:08:18.186267Z","shell.execute_reply":"2021-08-26T15:08:18.194539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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    slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)\n    print(slice_thickness)\n    for s in slices:\n        s.SliceThickness = slice_thickness\n        \n    return slices\n\ndef get_pixels_hu(scans):\n    image = np.stack([s.pixel_array for s in scans])\n    # Convert to int16 (from sometimes int16), \n    # should be possible as values should always be low enough (<32k)\n    image = image.astype(np.int16)\n\n    # Set outside-of-scan pixels to 1\n    # The intercept is usually -1024, so air is approximately 0\n    image[image == -2000] = 0\n    \n    # Convert to Hounsfield units (HU)\n    intercept = scans[0].RescaleIntercept\n    slope = scans[0].RescaleSlope\n    \n    if slope != 1:\n        image = slope * image.astype(np.float64)\n        image = image.astype(np.int16)\n        \n    image += np.int16(intercept)\n    \n    return np.array(image, dtype=np.int16)\n\nid=0\npatient = load_scan(data_path)\nimgs = get_pixels_hu(patient)\nprint(imgs.shape)\nnp.save(\"fullimages_%d.npy\" % (id), imgs)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.196505Z","iopub.execute_input":"2021-08-26T15:08:18.196938Z","iopub.status.idle":"2021-08-26T15:08:18.559843Z","shell.execute_reply.started":"2021-08-26T15:08:18.196907Z","shell.execute_reply":"2021-08-26T15:08:18.558651Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_used=\"fullimages_%d.npy\" % id\n\nimgs_to_process = np.load(file_used).astype(np.float64) \n\nplt.hist(imgs_to_process.flatten(), bins=50, color='c')\nplt.xlabel(\"Hounsfield Units (HU)\")\nplt.ylabel(\"Frequency\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:18.561347Z","iopub.execute_input":"2021-08-26T15:08:18.561794Z","iopub.status.idle":"2021-08-26T15:08:19.124105Z","shell.execute_reply.started":"2021-08-26T15:08:18.561744Z","shell.execute_reply":"2021-08-26T15:08:19.123328Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Slice Thickness: \", patient[0].SliceThickness)\nprint(\"Pixel Spacing (row, col):\",  (patient[0].PixelSpacing[0], patient[0].PixelSpacing[1]))","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:19.126227Z","iopub.execute_input":"2021-08-26T15:08:19.126670Z","iopub.status.idle":"2021-08-26T15:08:19.133226Z","shell.execute_reply.started":"2021-08-26T15:08:19.126637Z","shell.execute_reply":"2021-08-26T15:08:19.132247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import scipy\nid = 0\nimgs_to_process = np.load('fullimages_{}.npy'.format(id))\ndef resample(image, scan, new_spacing=[1,1,1]):\n    # Determine current pixel spacing\n    spacing = map(float, ([scan[0].SliceThickness] + list(scan[0].PixelSpacing)))\n    spacing = np.array(list(spacing))\n    resize_factor = spacing / new_spacing\n    print(resize_factor)\n    new_real_shape = image.shape * resize_factor\n    print(new_real_shape)\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    \n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)\n    \n    return image, new_spacing\n\nprint(\"Shape before resampling\\t\", imgs_to_process.shape)\nimgs_after_resamp, spacing = resample(imgs_to_process, patient, [1,1,1])\nprint(\"Shape after resampling\\t\", imgs_after_resamp.shape)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:19.135066Z","iopub.execute_input":"2021-08-26T15:08:19.135384Z","iopub.status.idle":"2021-08-26T15:08:24.091874Z","shell.execute_reply.started":"2021-08-26T15:08:19.135354Z","shell.execute_reply":"2021-08-26T15:08:24.090663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spacing","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:24.093439Z","iopub.execute_input":"2021-08-26T15:08:24.093880Z","iopub.status.idle":"2021-08-26T15:08:24.101980Z","shell.execute_reply.started":"2021-08-26T15:08:24.093835Z","shell.execute_reply":"2021-08-26T15:08:24.100424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_volume_axis(volume):\n    fig = plt.figure(figsize=(20, 20))\n    ax1 = fig.add_subplot(1, 3, 1)\n    ax1.imshow(volume[volume.shape[0]//2], cmap='gray')\n    \n    ax2 = fig.add_subplot(1, 3, 2)\n    ax2.imshow(volume[:, volume.shape[1]//2], cmap='gray')\n    \n    ax3 = fig.add_subplot(1, 3, 3)\n    ax3.imshow(volume[:, :, volume.shape[2]//2], cmap='gray')\n    fig.show()\n\nplot_volume_axis(imgs_after_resamp)","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:24.103115Z","iopub.execute_input":"2021-08-26T15:08:24.103461Z","iopub.status.idle":"2021-08-26T15:08:24.795352Z","shell.execute_reply.started":"2021-08-26T15:08:24.103429Z","shell.execute_reply":"2021-08-26T15:08:24.794348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skimage import measure\nimport plotly\nimport plotly.express as px\nfrom plotly.offline import  iplot\nfrom plotly.tools import FigureFactory as FF\n\n\ndef make_mesh(image, threshold=-300, step_size=1):\n\n    print(\"Transposing surface\")\n    p = image.transpose(2,1,0)\n    \n    print(\"Calculating surface\")\n    verts, faces, norm, val = measure.marching_cubes(p, threshold, step_size=step_size, allow_degenerate=True) \n    return verts, faces\n\ndef plotly_3d(verts, faces):\n    x,y,z = zip(*verts) \n    \n    print(\"Drawing\")\n    # Make the colormap single color since the axes are positional not intensity. \n#    colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']\n    colormap= px.colors.sequential.Brwnyl\n    backgroundcolor = 'slategray'\n    fig = FF.create_trisurf(x=x,\n                        y=y, \n                        z=z, \n                        plot_edges=False,\n                        colormap=colormap,\n                        simplices=faces,\n                        backgroundcolor=backgroundcolor,\n                        title=\"Interactive Visualization\")\n    iplot(fig)\n\n    \nv, f = make_mesh(imgs_after_resamp, threshold=10, step_size=10)\nplotly_3d(v, f)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:34:18.858970Z","iopub.execute_input":"2021-08-26T15:34:18.859384Z","iopub.status.idle":"2021-08-26T15:34:19.121733Z","shell.execute_reply.started":"2021-08-26T15:34:18.859346Z","shell.execute_reply":"2021-08-26T15:34:19.120359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<h1> Volume preprocessing Pipeline","metadata":{}},{"cell_type":"code","source":"def get_image_plane(data):\n    x1, y1, _, x2, y2, _ = [round(j) for j in data.ImageOrientationPatient]\n    cords = [x1, y1, x2, y2]\n\n    if cords == [1, 0, 0, 0]:\n        return 'Coronal'\n    elif cords == [1, 0, 0, 1]:\n        return 'Axial'\n    elif cords == [0, 1, 0, 0]:\n        return 'Sagittal'\n    else:\n        return 'Unknown'\n    \ndcm_img = pydicom.read_file('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00031/T1w/Image-22.dcm')\nplane = get_image_plane(dcm_img)\nplt.imshow(dcm_img.pixel_array, cmap='gray')\nplane","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:31.410667Z","iopub.execute_input":"2021-08-26T15:08:31.411104Z","iopub.status.idle":"2021-08-26T15:08:31.675677Z","shell.execute_reply.started":"2021-08-26T15:08:31.411064Z","shell.execute_reply":"2021-08-26T15:08:31.674601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_volume(study_dir):\n    imgs = []\n    dcm_dir = Path(study_dir)\n    dcm_paths = sorted(dcm_dir.glob(\"*.dcm\"), key=lambda x: int(x.stem.split(\"-\")[-1]))\n    positions = []\n    \n    for dcm_path in dcm_paths:\n        img = pydicom.dcmread(str(dcm_path))\n        imgs.append(img.pixel_array)\n        positions.append(img.ImagePositionPatient)\n        \n    plane = get_image_plane(img)\n    volume = np.stack(imgs)\n    \n    # reorder planes if needed and rotate volume\n    if plane == \"Coronal\":\n        if positions[0][1] < positions[-1][1]:\n            volume = volume[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        volume = volume.transpose((1, 0, 2))\n    elif plane == \"Sagittal\":\n        if positions[0][0] < positions[-1][0]:\n            volume = volume[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        volume = volume.transpose((1, 2, 0))\n        volume = np.rot90(volume, 2, axes=(1, 2))\n    elif plane == \"Axial\":\n        if positions[0][2] > positions[-1][2]:\n            volume = volume[::-1]\n            print(f\"{study_id} {scan_type} {plane} reordered\")\n        volume = np.rot90(volume, 2)\n    else:\n        raise ValueError(f\"Unknown plane {plane}\")\n    return volume, plane","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:08:31.676979Z","iopub.execute_input":"2021-08-26T15:08:31.677300Z","iopub.status.idle":"2021-08-26T15:08:31.689667Z","shell.execute_reply.started":"2021-08-26T15:08:31.677269Z","shell.execute_reply":"2021-08-26T15:08:31.688192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from skimage.transform import rescale, resize, downscale_local_mean\n\n\ndef calc_padding_inds(keep, padding):\n    s_edge_ind = np.argmax(keep)\n    e_edge_ind = np.argmax(keep[::-1])\n    keep_s = max(0, s_edge_ind-padding)\n    keep_e = min(len(keep)-e_edge_ind+padding, len(keep))\n    return keep_s, keep_e\n\ndef crop_volume(volume, padding=40):\n    if volume.sum() == 0:\n        return volume\n    keep = (volume.mean(axis=(0, 1)) > 0)\n    keep_s, keep_e = calc_padding_inds(keep, padding)\n    volume = volume[:, :, keep_s:keep_e]\n    \n    \n    keep = (volume.mean(axis=(0, 2)) > 0)\n    keep_s, keep_e = calc_padding_inds(keep, padding)\n    volume = volume[:, keep_s:keep_e]\n    \n    \n    keep = (volume.mean(axis=(1, 2)) > 0)\n    keep_s, keep_e = calc_padding_inds(keep, padding)\n\n    volume = volume[keep_s:keep_e]\n    return volume\n\ndef resize_volume(volume, sz=128, interpolation=cv2.INTER_LINEAR):\n    output = np.zeros((sz, sz, sz), dtype=np.float32)\n\n    if np.argmax(volume.shape) == 0:\n        for i, s in enumerate(np.linspace(0, volume.shape[0] - 1, sz)):\n            output[i] = cv2.resize(volume[int(s)], (sz, sz), interpolation=interpolation)\n\n    elif np.argmax(volume.shape) == 1:\n        for i, s in enumerate(np.linspace(0, volume.shape[1] - 1, sz)):\n            output[:, i] = cv2.resize(volume[:, int(s)], (sz, sz), interpolation=interpolation)\n\n    elif np.argmax(volume.shape) == 2:\n        for i, s in enumerate(np.linspace(0, volume.shape[2] - 1, sz)):\n            output[:, :, i] = cv2.resize(volume[:, :, int(s)], (sz, sz), interpolation=interpolation)\n\n    return output","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:22:36.102671Z","iopub.execute_input":"2021-08-26T15:22:36.103059Z","iopub.status.idle":"2021-08-26T15:22:36.119807Z","shell.execute_reply.started":"2021-08-26T15:22:36.103026Z","shell.execute_reply":"2021-08-26T15:22:36.118748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install --upgrade https://github.com/VincentStimper/mclahe/archive/numpy.zip\nfrom mclahe import mclahe","metadata":{"execution":{"iopub.status.busy":"2021-08-26T15:19:52.835791Z","iopub.execute_input":"2021-08-26T15:19:52.836198Z","iopub.status.idle":"2021-08-26T15:20:01.988061Z","shell.execute_reply.started":"2021-08-26T15:19:52.836135Z","shell.execute_reply":"2021-08-26T15:20:01.986848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<h2><b><i> Pipeline </i> </b></h2>","metadata":{}},{"cell_type":"code","source":"volume, _ = get_volume(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00045/FLAIR\")\nprint(volume.shape, volume.max())\n\nvolume = crop_volume(volume, padding=5)\nvolume = resize_volume(volume)\n\nprint(volume.shape, volume.max())\n\nvolume = mclahe(volume, n_bins=64)\nprint(volume.shape, volume.max())\n\nplot_volume_axis(volume)\nv, f = make_mesh(volume, threshold=0, step_size=2)\nplotly_3d(v, f)\n","metadata":{"execution":{"iopub.status.busy":"2021-08-26T16:10:06.231487Z","iopub.execute_input":"2021-08-26T16:10:06.231896Z","iopub.status.idle":"2021-08-26T16:10:11.022172Z","shell.execute_reply.started":"2021-08-26T16:10:06.231853Z","shell.execute_reply":"2021-08-26T16:10:11.020984Z"},"trusted":true},"execution_count":null,"outputs":[]}]}