{"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\nimport pandas as pd\nfrom glob import glob\nimport pydicom","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:36:34.352474Z","iopub.execute_input":"2021-08-03T10:36:34.352914Z","iopub.status.idle":"2021-08-03T10:36:34.359138Z","shell.execute_reply.started":"2021-08-03T10:36:34.352875Z","shell.execute_reply":"2021-08-03T10:36:34.35743Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH = '../input/rsna-miccai-brain-tumor-radiogenomic-classification'\naxis_move = {'sagittal': 0, 'coronal': 1, 'axial': 2}","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:36:34.361636Z","iopub.execute_input":"2021-08-03T10:36:34.362683Z","iopub.status.idle":"2021-08-03T10:36:34.371888Z","shell.execute_reply.started":"2021-08-03T10:36:34.362636Z","shell.execute_reply":"2021-08-03T10:36:34.370648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_dicom_files(cohort, case, mpMRI):\n    files_glob = f'{PATH}/{cohort}/{case}/{mpMRI}/*.dcm'\n    sorted_files = sorted(glob(files_glob),key=lambda f: int(f.split('Image-')[1].split('.')[0]))\n    return [pydicom.read_file(f) for f in sorted_files]\n\ndef image_orientation(dicom):\n    rt = 'unkown'\n    # https://www.kaggle.com/davidbroberts/determining-mr-image-planes\n    (x1,y1,_,x2,y2,_) = [round(v) for v in dicom.ImageOrientationPatient]\n    if (x1,y1,x2,y2) == (1,0,0,0):\n        rt = 'coronal'\n    if (x1,y1,x2,y2) == (1,0,0,1):\n        rt = 'axial'\n    if (x1,y1,x2,y2) == (0,1,0,0):\n        rt = 'sagittal'\n    \n    if rt == 'unkown':\n        raise ValueError(f'unkown ImageOrientationPatient: {dicom.ImageOrientationPatient}')\n        \n    return rt\n\ndef stats_images(images):\n    pixels = images.ravel()\n    noncero_pixels = pixels[np.nonzero(pixels)]\n    mean = np.mean(noncero_pixels) # <-- collect \n    std = np.std(noncero_pixels) # <-- collect \n    return (mean,std)\n\ndef top_brilliant_image(images, mean, std):\n    level = mean + 3*std\n    non_cero_pixels = np.count_nonzero(np.reshape(images, (images.shape[0],-1)) > level,  axis=1)\n    top_image = np.argsort(non_cero_pixels)[::-1][0]\n    return top_image\n\ndef top_brilliant_line(image, mean, std, axis):\n    non_cero_pixels = np.count_nonzero(image > mean+3*std,  axis=axis)\n    top_line = np.argsort(non_cero_pixels)[::-1][0]\n    return top_line\n\ndef normalize_image(image, mean, std):\n    return (image-mean)/std\n\ndef calc_center(dicom_file, r, c):\n    orientation = image_orientation(dicom_file)\n    if orientation == 'coronal':\n        center = [dicom_file.ImagePositionPatient[0] + dicom_file.PixelSpacing[0] * c,\n                  dicom_file.ImagePositionPatient[1],\n                  dicom_file.ImagePositionPatient[2] - dicom_file.PixelSpacing[1] * r]\n\n    if orientation == 'sagittal':\n        center = [dicom_file.ImagePositionPatient[0],\n                  dicom_file.ImagePositionPatient[1] + dicom_file.PixelSpacing[0] * c,\n                  dicom_file.ImagePositionPatient[2] - dicom_file.PixelSpacing[1] * r]\n        \n    if orientation == 'axial':\n        center = [dicom_file.ImagePositionPatient[0] + dicom_file.PixelSpacing[0] * c,\n                  dicom_file.ImagePositionPatient[1] + dicom_file.PixelSpacing[0] * r,\n                  dicom_file.ImagePositionPatient[2]]\n\n    return center\n\ndef find_nearest_scan(dicom_files, center):\n    orientation = image_orientation(dicom_files[0])\n    a = np.array([f.ImagePositionPatient for f in dicom_files])\n    scan = np.argsort(np.abs(a - center),axis=0)[0][axis_move[orientation]]\n    return scan","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:36:34.375724Z","iopub.execute_input":"2021-08-03T10:36:34.376678Z","iopub.status.idle":"2021-08-03T10:36:34.398553Z","shell.execute_reply.started":"2021-08-03T10:36:34.376635Z","shell.execute_reply":"2021-08-03T10:36:34.397331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cohort = 'train'\ncase = '00581'\n\n# To Do: review cases 00266, 00457, 00611, 00631\n#   - case 00266 vs 00737\n\nflair_dicom_files = read_dicom_files(cohort, case, 'FLAIR')\nt1w_dicom_files = read_dicom_files(cohort, case, 'T1w')\nt1wce_dicom_files = read_dicom_files(cohort, case, 'T1wCE')\nt2w_dicom_files = read_dicom_files(cohort, case, 'T2w')","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:37:15.143794Z","iopub.execute_input":"2021-08-03T10:37:15.144274Z","iopub.status.idle":"2021-08-03T10:37:16.739076Z","shell.execute_reply.started":"2021-08-03T10:37:15.144226Z","shell.execute_reply":"2021-08-03T10:37:16.73792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flair_orientation = image_orientation(flair_dicom_files[0])\nflair_nscans = len(flair_dicom_files)\nt1w_orientation = image_orientation(t1w_dicom_files[0])\nt1w_nscans = len(t1w_dicom_files)\nt1wce_orientation = image_orientation(t1wce_dicom_files[0])\nt1wce_nscans = len(t1wce_dicom_files)\nt2w_orientation = image_orientation(t2w_dicom_files[0])\nt2w_nscans = len(t2w_dicom_files)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:38:59.451548Z","iopub.execute_input":"2021-08-03T10:38:59.451962Z","iopub.status.idle":"2021-08-03T10:38:59.459373Z","shell.execute_reply.started":"2021-08-03T10:38:59.45193Z","shell.execute_reply":"2021-08-03T10:38:59.458189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"FLAIR: {flair_orientation}, {flair_nscans} scans\")\nprint(f\"T1w: {t1w_orientation}, {t1w_nscans} scans\")\nprint(f\"T1wce: {t1wce_orientation}, {t1wce_nscans} scans\")\nprint(f\"T2w: {t2w_orientation}, {t2w_nscans} scans\")","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:39:15.362907Z","iopub.execute_input":"2021-08-03T10:39:15.363271Z","iopub.status.idle":"2021-08-03T10:39:15.372986Z","shell.execute_reply.started":"2021-08-03T10:39:15.363239Z","shell.execute_reply":"2021-08-03T10:39:15.371177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(flair_dicom_files[0].PatientID,t1w_dicom_files[0].PatientID,t1wce_dicom_files[0].PatientID,t2w_dicom_files[0].PatientID)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:39:52.826183Z","iopub.execute_input":"2021-08-03T10:39:52.826547Z","iopub.status.idle":"2021-08-03T10:39:52.837295Z","shell.execute_reply.started":"2021-08-03T10:39:52.826516Z","shell.execute_reply":"2021-08-03T10:39:52.835925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flair_images = np.array([s.pixel_array for s in flair_dicom_files])\n#-- tiw_images = np.array([s.pixel_array for s in t1w_dicom_files])","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:40:33.383187Z","iopub.execute_input":"2021-08-03T10:40:33.383627Z","iopub.status.idle":"2021-08-03T10:40:33.39921Z","shell.execute_reply.started":"2021-08-03T10:40:33.383563Z","shell.execute_reply":"2021-08-03T10:40:33.397548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\n\nfig = px.imshow(flair_images, animation_frame=0, binary_string=True, labels=dict(animation_frame=\"scan\"), height=800)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:41:19.071376Z","iopub.execute_input":"2021-08-03T10:41:19.071789Z","iopub.status.idle":"2021-08-03T10:41:25.328482Z","shell.execute_reply.started":"2021-08-03T10:41:19.071748Z","shell.execute_reply":"2021-08-03T10:41:25.327016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(min,max)=(flair_images.min(), flair_images.max()) # <-- collect / normalize","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:42:28.843545Z","iopub.execute_input":"2021-08-03T10:42:28.843981Z","iopub.status.idle":"2021-08-03T10:42:28.859765Z","shell.execute_reply.started":"2021-08-03T10:42:28.843949Z","shell.execute_reply":"2021-08-03T10:42:28.858346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport matplotlib.lines as lines\nfrom scipy import stats\nplt.style.use('ggplot')\n\n\na=flair_images.ravel()\na=a[np.nonzero(a)]\nmode = stats.mode(a)\n\n(mean,std) = stats_images(flair_images)\nprint(mode.mode, mean, std)\n\npixels = flair_images.ravel()\nnoncero_pixels = pixels[np.nonzero(pixels)]\n\nplt.figure(figsize = (20,5))\nplt.hist(noncero_pixels, max)\nplt.vlines(mean, 0, mode.count, colors='b')\nplt.vlines(mean+std, 0, mode.count, colors='b', linestyles='dotted')\nplt.vlines(mean+2*std, 0, mode.count, colors='b', linestyles='dashed')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:43:10.08397Z","iopub.execute_input":"2021-08-03T10:43:10.084335Z","iopub.status.idle":"2021-08-03T10:43:11.683641Z","shell.execute_reply.started":"2021-08-03T10:43:10.084302Z","shell.execute_reply":"2021-08-03T10:43:11.682225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top = top_brilliant_image(flair_images, mean, std)\ntop","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:43:43.53148Z","iopub.execute_input":"2021-08-03T10:43:43.531942Z","iopub.status.idle":"2021-08-03T10:43:43.554018Z","shell.execute_reply.started":"2021-08-03T10:43:43.531853Z","shell.execute_reply":"2021-08-03T10:43:43.552691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pydicom.pixel_data_handlers.util import apply_voi_lut\n\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize = (20,10))\nfig.suptitle('normalized image vs VOI LUT image')\n\nimage = flair_images[top]\nimage = normalize_image(image, mean, std)\nax1.imshow(image, cmap = plt.cm.gray)\n\nim = apply_voi_lut(flair_dicom_files[top].pixel_array, flair_dicom_files[top])\nax2.imshow(im, cmap = plt.cm.gray)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:44:25.552505Z","iopub.execute_input":"2021-08-03T10:44:25.552942Z","iopub.status.idle":"2021-08-03T10:44:25.996429Z","shell.execute_reply.started":"2021-08-03T10:44:25.552894Z","shell.execute_reply":"2021-08-03T10:44:25.995151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"(top,flair_dicom_files[top].ImagePositionPatient)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:45:00.023511Z","iopub.execute_input":"2021-08-03T10:45:00.02394Z","iopub.status.idle":"2021-08-03T10:45:00.033595Z","shell.execute_reply.started":"2021-08-03T10:45:00.023909Z","shell.execute_reply":"2021-08-03T10:45:00.032188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rtop = top_brilliant_line(flair_images[top], mean, std, axis=1)\nctop = top_brilliant_line(flair_images[top], mean, std, axis=0)\n\n(rtop,ctop)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:45:34.322171Z","iopub.execute_input":"2021-08-03T10:45:34.322516Z","iopub.status.idle":"2021-08-03T10:45:34.333323Z","shell.execute_reply.started":"2021-08-03T10:45:34.322485Z","shell.execute_reply":"2021-08-03T10:45:34.331754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"flair_dicom_files[top]","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:46:11.998051Z","iopub.execute_input":"2021-08-03T10:46:11.9985Z","iopub.status.idle":"2021-08-03T10:46:12.009163Z","shell.execute_reply.started":"2021-08-03T10:46:11.998469Z","shell.execute_reply":"2021-08-03T10:46:12.007727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"center = calc_center(flair_dicom_files[top], rtop, ctop)\n\ncenter","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:47:16.589648Z","iopub.execute_input":"2021-08-03T10:47:16.590018Z","iopub.status.idle":"2021-08-03T10:47:16.598441Z","shell.execute_reply.started":"2021-08-03T10:47:16.589987Z","shell.execute_reply":"2021-08-03T10:47:16.59667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scan_t1wce = find_nearest_scan(t1wce_dicom_files, center)\n(t1wce_orientation, scan_t1wce, t1wce_dicom_files[scan_t1wce].ImagePositionPatient)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:47:50.869827Z","iopub.execute_input":"2021-08-03T10:47:50.870239Z","iopub.status.idle":"2021-08-03T10:47:50.888552Z","shell.execute_reply.started":"2021-08-03T10:47:50.870191Z","shell.execute_reply":"2021-08-03T10:47:50.886968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scan_t1w = find_nearest_scan(t1w_dicom_files, center)\n(t1w_orientation, scan_t1w, t1w_dicom_files[scan_t1w].ImagePositionPatient)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:48:26.869631Z","iopub.execute_input":"2021-08-03T10:48:26.870033Z","iopub.status.idle":"2021-08-03T10:48:26.88034Z","shell.execute_reply.started":"2021-08-03T10:48:26.87Z","shell.execute_reply":"2021-08-03T10:48:26.879022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scan_t2w = find_nearest_scan(t2w_dicom_files, center)\n(t2w_orientation, scan_t2w, t2w_dicom_files[scan_t2w].ImagePositionPatient)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:48:56.630972Z","iopub.execute_input":"2021-08-03T10:48:56.631362Z","iopub.status.idle":"2021-08-03T10:48:56.643786Z","shell.execute_reply.started":"2021-08-03T10:48:56.631324Z","shell.execute_reply":"2021-08-03T10:48:56.642423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2, 2, figsize = (20,20))\nfig.suptitle('images')\n\nim = apply_voi_lut(flair_dicom_files[top].pixel_array, flair_dicom_files[top])\nax1.imshow(im, cmap = plt.cm.gray)\nax1.set_title(f'FLAIR #scan {top}')\n\nim = apply_voi_lut(t1w_dicom_files[scan_t1w].pixel_array, t1w_dicom_files[scan_t1w])\nax2.imshow(im, cmap = plt.cm.gray)\nax2.set_title(f'T1w #scan {scan_t1w}')\n\nim = apply_voi_lut(t1wce_dicom_files[scan_t1wce].pixel_array, t1wce_dicom_files[scan_t1wce])\nax3.imshow(im, cmap = plt.cm.gray)\nax3.set_title(f'T1wCE #scan {scan_t1wce}')\n\nim = apply_voi_lut(t2w_dicom_files[scan_t2w].pixel_array, t2w_dicom_files[scan_t2w])\nax4.imshow(im, cmap = plt.cm.gray)\nax4.set_title(f'T2w #scan {scan_t2w}')\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:49:37.115685Z","iopub.execute_input":"2021-08-03T10:49:37.11612Z","iopub.status.idle":"2021-08-03T10:49:37.960697Z","shell.execute_reply.started":"2021-08-03T10:49:37.116077Z","shell.execute_reply":"2021-08-03T10:49:37.959467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def process_case_and_plot(cohort, case):\n    flair_dicom_files = read_dicom_files(cohort, case, 'FLAIR')\n    t1w_dicom_files = read_dicom_files(cohort, case, 'T1w')\n    t1wce_dicom_files = read_dicom_files(cohort, case, 'T1wCE')\n    t2w_dicom_files = read_dicom_files(cohort, case, 'T2w')\n    \n    flair_images = np.array([s.pixel_array for s in flair_dicom_files])\n    (mean,std) = stats_images(flair_images)\n    \n    top = top_brilliant_image(flair_images, mean, std)\n    rtop = top_brilliant_line(flair_images[top], mean, std, axis=1)\n    ctop = top_brilliant_line(flair_images[top], mean, std, axis=0)\n\n    center = calc_center(flair_dicom_files[top], rtop, ctop)\n    \n    flair_image = normalize_image(flair_images[top], mean, std)\n    \n    scan_t1w = find_nearest_scan(t1w_dicom_files, center)\n    scan_t1wce = find_nearest_scan(t1wce_dicom_files, center)\n    scan_t2w = find_nearest_scan(t2w_dicom_files, center)\n\n    t1w_images = np.array([s.pixel_array for s in t1w_dicom_files])\n    (mean,std) = stats_images(t1w_images)\n    t1w_image = normalize_image(t1w_images[scan_t1w], mean, std)\n    \n    t1wce_images = np.array([s.pixel_array for s in t1wce_dicom_files])\n    (mean,std) = stats_images(t1wce_images)\n    t1wce_image = normalize_image(t1wce_images[scan_t1wce], mean, std)\n    \n    t2w_images = np.array([s.pixel_array for s in t2w_dicom_files])\n    (mean,std) = stats_images(t2w_images)\n    t2w_image = normalize_image(t2w_images[scan_t2w], mean, std)\n\n    fig, (ax1,ax2,ax3,ax4) = plt.subplots(1, 4, figsize = (20,5))\n    fig.suptitle(f'Case {case}')\n    \n    ax1.imshow(flair_image, cmap = plt.cm.gray)\n    ax1.set_title(f'FLAIR #scan {top}')\n    ax1.grid(False)\n\n    ax2.imshow(t1w_image, cmap = plt.cm.gray)\n    ax2.set_title(f'T1w #scan {scan_t1w}')\n    ax2.grid(False)\n    \n    ax3.imshow(t1wce_image, cmap = plt.cm.gray)\n    ax3.set_title(f'T1wCE #scan {scan_t1wce}')\n    ax3.grid(False)\n\n    ax4.imshow(t2w_image, cmap = plt.cm.gray)\n    ax4.set_title(f'T2w #scan {scan_t2w}')\n    ax4.grid(False)\n\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:50:46.472061Z","iopub.execute_input":"2021-08-03T10:50:46.472541Z","iopub.status.idle":"2021-08-03T10:50:46.487658Z","shell.execute_reply.started":"2021-08-03T10:50:46.472489Z","shell.execute_reply":"2021-08-03T10:50:46.486589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv', converters = {'BraTS21ID': str,})","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:51:18.529232Z","iopub.execute_input":"2021-08-03T10:51:18.529608Z","iopub.status.idle":"2021-08-03T10:51:18.548294Z","shell.execute_reply.started":"2021-08-03T10:51:18.529565Z","shell.execute_reply":"2021-08-03T10:51:18.547227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cohort = 'train'\nfor case in train.sample(10).BraTS21ID:  \n    process_case_and_plot(cohort, case)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:51:50.638105Z","iopub.execute_input":"2021-08-03T10:51:50.638518Z","iopub.status.idle":"2021-08-03T10:52:37.289835Z","shell.execute_reply.started":"2021-08-03T10:51:50.638485Z","shell.execute_reply":"2021-08-03T10:52:37.288878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test = pd.read_csv('../input/rsna-miccai-brain-tumor-radiogenomic-classification/sample_submission.csv', converters = {'BraTS21ID': str,})","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:52:45.229371Z","iopub.execute_input":"2021-08-03T10:52:45.229797Z","iopub.status.idle":"2021-08-03T10:52:45.24162Z","shell.execute_reply.started":"2021-08-03T10:52:45.229752Z","shell.execute_reply":"2021-08-03T10:52:45.240228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cohort = 'test'\nfor case in test.sample(10).BraTS21ID:\n    process_case_and_plot(cohort, case)","metadata":{"execution":{"iopub.status.busy":"2021-08-03T10:53:17.849245Z","iopub.execute_input":"2021-08-03T10:53:17.849606Z","iopub.status.idle":"2021-08-03T10:54:27.167893Z","shell.execute_reply.started":"2021-08-03T10:53:17.849577Z","shell.execute_reply":"2021-08-03T10:54:27.166635Z"},"trusted":true},"execution_count":null,"outputs":[]}]}