{"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":"markdown","source":"# DAT 490 - Exploratory Data Analysis\n*by Max Murphy*","metadata":{}},{"cell_type":"markdown","source":"## Import All External Packages\nThe categories of library and the packages we are using are as follows:","metadata":{}},{"cell_type":"code","source":"# Python System Libraries\nimport os\nimport random\nimport re\nimport glob\n\n# Libraries to Import Medical Data\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\n\n# Data Analysis Libraries\nimport numpy as np\nimport pandas as pd\n\n# Data Visualization Libraries\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:12:59.566766Z","iopub.execute_input":"2021-10-08T06:12:59.567093Z","iopub.status.idle":"2021-10-08T06:12:59.573189Z","shell.execute_reply.started":"2021-10-08T06:12:59.567050Z","shell.execute_reply":"2021-10-08T06:12:59.572352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Breaking Down the Key Components of the Given Data\n\nThere are three main components to the RSNA-MICCIA Dataset:\n1. A list of labels with patient IDs that indicate whether or not the patient has the MGMT genetic marker\n2. A file structure containing visual data delineated by patient ID\n3. Lastly, the actual visual information","metadata":{}},{"cell_type":"markdown","source":"### Component 1: MGMT Values","metadata":{}},{"cell_type":"code","source":"train_label_df = pd.read_csv(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train_labels.csv\")\ntrain_label_df.head()","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:12:59.574930Z","iopub.execute_input":"2021-10-08T06:12:59.575780Z","iopub.status.idle":"2021-10-08T06:12:59.603427Z","shell.execute_reply.started":"2021-10-08T06:12:59.575736Z","shell.execute_reply":"2021-10-08T06:12:59.602222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To clarify some of the training data, the \"BraTS21ID\" is the patient identifier, and a value of 1 in the \"MGMT_value\" means the patient actually has the MGMT marker and 0 is where the patient does not.","metadata":{}},{"cell_type":"markdown","source":"#### Key Point 1: Looking at the distribution of MGMT in our patient training set","metadata":{}},{"cell_type":"code","source":"label_replacement = {1: 'Yes', 0: 'No'}\nplt.figure(figsize=(5,7))\nax = train_label_df[\"MGMT_value\"].replace(label_replacement).value_counts().plot(kind=\"bar\", rot=0, xlabel=\"Whether or not the patient has the MGMT genotype\", ylabel=\"Patient Count\")\nax.set_title(\"Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n Patient Counts based on Genotype\")\nrects = ax.patches\n# Make some labels.\nlabels = [f\"{i}\" for i in train_label_df[\"MGMT_value\"].value_counts()]\n\nfor rect, label in zip(rects, labels):\n    height = rect.get_height()\n    ax.text(\n        rect.get_x() + rect.get_width() / 2, height + 5, label, ha=\"center\", va=\"bottom\"\n    )","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:12:59.605230Z","iopub.execute_input":"2021-10-08T06:12:59.605596Z","iopub.status.idle":"2021-10-08T06:12:59.849271Z","shell.execute_reply.started":"2021-10-08T06:12:59.605560Z","shell.execute_reply":"2021-10-08T06:12:59.848324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_label_df[\"MGMT_value\"].replace(label_replacement).value_counts()","metadata":{"_kg_hide-input":false,"execution":{"iopub.status.busy":"2021-10-08T06:12:59.850456Z","iopub.execute_input":"2021-10-08T06:12:59.850882Z","iopub.status.idle":"2021-10-08T06:12:59.859289Z","shell.execute_reply.started":"2021-10-08T06:12:59.850844Z","shell.execute_reply":"2021-10-08T06:12:59.858488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we can see there is a small skew towards a patient having the marker. This may just be due to some randomness when the training sample was created for this dataset. This distribution may change in both sets of our training data.","metadata":{}},{"cell_type":"markdown","source":"### Component 2: Patient File Structure","metadata":{}},{"cell_type":"code","source":"n = 5\nrootdir = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train\"\nfor file in os.listdir(rootdir)[:n]:\n    d = os.path.join(rootdir, file)\n    if os.path.isdir(d):\n        print(d)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:12:59.861080Z","iopub.execute_input":"2021-10-08T06:12:59.861414Z","iopub.status.idle":"2021-10-08T06:12:59.875022Z","shell.execute_reply.started":"2021-10-08T06:12:59.861360Z","shell.execute_reply":"2021-10-08T06:12:59.874214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This is the file structure of the train dataset, as we can see that \"BraTS21ID\" is the title of each folder in the training folder. These folders contain the 4 different brain scans of the patient with that ID.","metadata":{}},{"cell_type":"markdown","source":"#### Key Aspect 1 : How many patients are included in the training dataset?","metadata":{}},{"cell_type":"markdown","source":"Next let's figure out how many patients are included in our training data set. This can be found by counting the number of directories in the \"train\" folder of the dataset.","metadata":{}},{"cell_type":"code","source":"path = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train\"\nsum(os.path.isdir(os.path.join(path, i)) for i in os.listdir(path))","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:12:59.876415Z","iopub.execute_input":"2021-10-08T06:12:59.876857Z","iopub.status.idle":"2021-10-08T06:13:00.093633Z","shell.execute_reply.started":"2021-10-08T06:12:59.876826Z","shell.execute_reply":"2021-10-08T06:13:00.092744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Key Aspect 2: Folder Structure for Each Patient\n\nAnother thing we should note is the internal structure of each of the patient folders for our future reference.","metadata":{}},{"cell_type":"code","source":"rootdir = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000\"\nfor file in os.listdir(rootdir):\n    d = os.path.join(rootdir, file)\n    if os.path.isdir(d):\n        print(d)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:00.094840Z","iopub.execute_input":"2021-10-08T06:13:00.095049Z","iopub.status.idle":"2021-10-08T06:13:00.102250Z","shell.execute_reply.started":"2021-10-08T06:13:00.095025Z","shell.execute_reply":"2021-10-08T06:13:00.101473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we can see that in a given folder there will be 4 sub folders covering the 4 different types of MRI scans that we talked about earlier in the literature review, and will go into more detail in the later sections of this document.","metadata":{}},{"cell_type":"markdown","source":"### Component 3: Patient MRI Scans","metadata":{}},{"cell_type":"markdown","source":"Within the patient MRI scans there are important differences that I outlined in the literature review. I would like to briefly go through some of the key points of information about these scans.","metadata":{}},{"cell_type":"markdown","source":"#### Key Point 1 : Varying MRI Types\n\nWithin each patient MRI folder, we have 4 subfolders. These folders are named according to the identifier associated with which MRI scan type the radiologist used to obtain the images within the folder. Let's take a look at the visual differences of these scans first.","metadata":{}},{"cell_type":"code","source":"def load_dicom(path):\n    dicom = pydicom.read_file(path)\n    data = dicom.pixel_array\n    data = data - np.min(data)\n    if np.max(data) != 0:\n        data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return data\n\ndef visualize_sample(\n    brats21id, \n    slice_i,\n    mgmt_value,\n    types=(\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\")\n):\n    plt.figure(figsize=(16, 5))\n    patient_path = os.path.join(\n        \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/\", \n        str(brats21id).zfill(5),\n    )\n    for i, t in enumerate(types, 1):\n        t_paths = sorted(\n            glob.glob(os.path.join(patient_path, t, \"*\")), \n            key=lambda x: int(x[:-4].split(\"-\")[-1]),\n        )\n        data = load_dicom(t_paths[int(len(t_paths) * slice_i)])\n        plt.subplot(1, 4, i)\n        plt.imshow(data, cmap=\"gray\")\n        plt.title(f\"{t}\", fontsize=16)\n        plt.axis(\"off\")\n\n    plt.suptitle(f\"Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n MGMT_value: {mgmt_value}\", fontsize=16)\n    plt.show()\n\nfor i in random.sample(range(train_label_df.shape[0]), 1):\n    _brats21id = train_label_df.iloc[i][\"BraTS21ID\"]\n    _mgmt_value = train_label_df.iloc[i][\"MGMT_value\"]\n    visualize_sample(brats21id=_brats21id, mgmt_value=_mgmt_value, slice_i=0.5)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:14:00.164733Z","iopub.execute_input":"2021-10-08T06:14:00.165029Z","iopub.status.idle":"2021-10-08T06:14:00.629723Z","shell.execute_reply.started":"2021-10-08T06:14:00.165000Z","shell.execute_reply":"2021-10-08T06:14:00.628773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Code Attribution: https://www.kaggle.com/ihelon/brain-tumor-eda-with-animations-and-modeling","metadata":{}},{"cell_type":"markdown","source":"Here we can see the 4 MRI types clearly. There are 3 T1 MRI scans and a single T2 MRI scan. You'll notice there are small differences in the level of contrast and brightness on different areas of the photo. These differences are because of the different imaging techniques used to take the physical photos.","metadata":{}},{"cell_type":"markdown","source":"**Important Aspect: Differing Image Counts**\n\nSomething that we should also check for is whether or not we have a standard number of images in each MRI type. We can do this by deriving some additional features in our data frame and using those features to answer our question. We need to derive some information about the folder names, the file paths, and then the number of images in each of the patient's subfolders to figure out if there are differing counts of MRI images across our patient base.","metadata":{}},{"cell_type":"code","source":"# Shorten the data frame name down because we are going to use it a lot.\ntrain = train_label_df\n\n# Add an image folder feature to the dataframe\ntrain[\"image_folder_name\"] = ['{0:05d}'.format(s) for s in train[\"BraTS21ID\"]]\n\n# Define where the training data is located\ntrain_file_path = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train\"\n\n# Create a concatenated path in the dataframe for the image location of the data.\ntrain[\"image_folder_path\"] = [os.path.join(train_file_path, x) for x in train['image_folder_name']]\n\n# Scan Types\nmri_types = [\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\"]\n\n# Find the number of files in each of the scan type folder\nfor mri_type in mri_types:\n    train[mri_type + \"_count\"] = [len(os.listdir(os.path.join(train['image_folder_path'].iloc[s], mri_type))) for s in range(len(train))]\n    \n# Print the current training dataframe\ntrain.head(5)","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:00.634100Z","iopub.execute_input":"2021-10-08T06:13:00.634451Z","iopub.status.idle":"2021-10-08T06:13:02.001385Z","shell.execute_reply.started":"2021-10-08T06:13:00.634409Z","shell.execute_reply":"2021-10-08T06:13:02.000566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see my suspicions are true, there are differing numbers of images between MRI types *and* there are differing numbers of images across patients in the same MRI type. Let's dive a bit deeper into how widely these vary, because we may need to adjust our methodology to account for these differences.","metadata":{}},{"cell_type":"code","source":"rename_dict = {\"FLAIR_count\": \"FLAIR\", \"T1w_count\": \"T1w\", \"T1wCE_count\": \"T1wCE\", \"T2w_count\": \"T2w\"}\nplt.figure(figsize=(5,7))\nax = train[['FLAIR_count', 'T1w_count', 'T1wCE_count', \"T2w_count\"]].rename(columns=rename_dict).plot(kind=\"box\")\nax.set_xlabel(\"MRI Imaging Type\")\nax.set_ylabel(\"Image Count per Patient\")\nax.set_title(\"Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n Box Plot of Amount of Images of Varying Type in Patient MRI data\")","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:02.002443Z","iopub.execute_input":"2021-10-08T06:13:02.002952Z","iopub.status.idle":"2021-10-08T06:13:02.288287Z","shell.execute_reply.started":"2021-10-08T06:13:02.002925Z","shell.execute_reply":"2021-10-08T06:13:02.287245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Key Point 2 : MRI Imaging Planes\n\nAnother point that we must take a look at is the differing planes with which the images are taken. There are three planes as described in the literature review: Axial (top to bottom), Coronal (front to back), and Sagittal (left to right). These types are all included in arbitrary patterns among each of the different MRI scan types for the patients, and may serve as an important reason why the counts for each of the scan types are so radically different. Luckily the patient image orientations are included in the DiCOM data and we will extract that from the file below.\n\n![MRI Axial Type Diagram](https://my-ms.org/images/mri_planes_gnu.jpg)\n\n> Image Retrieved from : [https://my-ms.org/images/mri_planes_gnu.jpg](https://my-ms.org/images/mri_planes_gnu.jpg)","metadata":{}},{"cell_type":"code","source":"def findImageOrientation(image):\n    # Pull the image orientation data from the DiCOM file\n    (x1,y1,_,x2,y2,_) = [round(value) for value in image.ImageOrientationPatient]\n    if (x1,x2,y1,y2) == (1,0,0,0):\n        return \"Coronal\"\n    elif (x1,x2,y1,y2) == (1,0,0,1):\n        return \"Axial\"\n    elif (x1,x2,y1,y2) == (0,0,1,0):\n        return \"Sagittal\"\n    else:\n        return \"Unknown\"\n\n# Scan Types\nmri_types = [\"FLAIR\", \"T1w\", \"T1wCE\", \"T2w\"]\n\nfor mri_type in mri_types:\n    train[mri_type + \"_axis\"] = [findImageOrientation(pydicom.read_file(os.path.join(os.path.join(train['image_folder_path'].iloc[s], mri_type), os.listdir(os.path.join(train['image_folder_path'].iloc[s], mri_type))[0]))) for s in range(len(train))]\n\nprint(train[['FLAIR_axis', 'T1w_axis', 'T1wCE_axis', 'T2w_axis']].head())","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:02.290468Z","iopub.execute_input":"2021-10-08T06:13:02.290685Z","iopub.status.idle":"2021-10-08T06:13:08.172429Z","shell.execute_reply.started":"2021-10-08T06:13:02.290659Z","shell.execute_reply":"2021-10-08T06:13:08.170814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Code Inspiration: https://www.kaggle.com/davidbroberts/determining-mr-image-planes","metadata":{}},{"cell_type":"markdown","source":"Now that we have the imaging axis extracted from the DiCOM, let's see the distribution of axes across the scan types. We are going to create bar charts to show this.","metadata":{}},{"cell_type":"code","source":"plt.figure()\nrename = {\"FLAIR_axis\": \"FLAIR\", \"T1w_axis\": \"T1w\", \"T1wCE_axis\": \"T1wCE\", \"T2w_axis\": \"T2w\"}\nmri_types_df = train[[mri_type + '_axis' for mri_type in mri_types]].apply(pd.Series.value_counts).rename(columns=rename)\nmri_types_df.plot(kind=\"bar\", rot=0, xlabel=\"MRI Scanning Axis\", ylabel=\"Image Count\", title=\"Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n Image Count Grouped by MRI Scan Axis and Colored by MRI Type\")","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:08.173515Z","iopub.execute_input":"2021-10-08T06:13:08.173735Z","iopub.status.idle":"2021-10-08T06:13:08.483224Z","shell.execute_reply.started":"2021-10-08T06:13:08.173710Z","shell.execute_reply":"2021-10-08T06:13:08.482339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As you can see some of the scan types vary wildly between axes while there are a few exceptions. The most notable is T1w, which was relatively consistent with the axial axis. The two imaging types with the most variation across axis were T2w and FLAIR. What's notable about these pairs is that they have the same distribution of axes. This may be because they are correlated across patients. ","metadata":{}},{"cell_type":"markdown","source":"#### Key Point 3 : Usable Images\n\nAnother point that I would like to note about the data is the amount of unusable or non-important images in each of the patient-imaging-type sets. The sets typically start and end with large amounts of black space and climax with images that are denoted with a \"brilliant zone.\" These brilliant zones show the tumor really well and will likely be the most useful images to look at because they show the subject of our analysis with the greatest fidelity of detail.\n\nSo how do we actually accomplish this kind of analysis? We can use a useful tool from photography called image histograms. A typical image histogram looks like this and shows each of the color channels of the image and the distribution of pixel luminosity across all pixels in the image. \n\n![Example of Image Histogram with 3 color channels](https://www.dummies.com/wp-content/uploads/332882.image1.jpg)\n\nWe can do the same using our existing data analysis tools, and a nice thing that I'm sure you noticed in the image set is that there is a single color channel making the image monochromatic and cutting the amount of dimensions in our histogram. First let's see the raw pixel data and work towards a meaningful analysis of the image.","metadata":{}},{"cell_type":"code","source":"example_file_path = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00000/FLAIR/Image-200.dcm\"\nprint(\"Raw Pixel Array\")\nprint(pydicom.read_file(example_file_path).pixel_array)\nprint(\"Max Pixel Value: \" + str(pydicom.read_file(example_file_path).pixel_array.max()))","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:08.484641Z","iopub.execute_input":"2021-10-08T06:13:08.484863Z","iopub.status.idle":"2021-10-08T06:13:08.499425Z","shell.execute_reply.started":"2021-10-08T06:13:08.484838Z","shell.execute_reply":"2021-10-08T06:13:08.498630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we can see the pixel values in the raw pixel array. Each sub-array signifies a row of pixels in an image, the rows are all of a standard length. The higher the pixel value the more luminosity the pixel has. Now let's unravel these arrays and plot the imaging values on a histogram.","metadata":{}},{"cell_type":"code","source":"def plot_image_histogram(image):\n    ## Get the pixels all on to a single array\n    pixels = np.array(image.pixel_array).ravel()\n    \n    ## Create a two sided figure\n    fig, (ax_image, ax_hist) = plt.subplots(1, 2, figsize = (20,4), gridspec_kw={'width_ratios': [1, 4]})\n    fig.suptitle(f'Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n scan # ({image.InstanceNumber})')\n    \n    ## Begin working on the histogram portion of the graph\n    ### Filter the zero pixels (greatly messes with the scale of the graph)\n    non_zero_pixels = np.nonzero(pixels)\n    ### Run basic statistics to normalize the pixel data (makes data presentable)\n    mean = np.mean(non_zero_pixels)\n    std = np.std(non_zero_pixels)\n    ### Normalize the pixels\n    norm_non_zero_pixels = (non_zero_pixels - mean)/std\n    if len(non_zero_pixels[0]) != 0:\n        ax_hist.hist(norm_non_zero_pixels.flatten(), bins=200)\n        ax_hist.set_xlim(-5,5)\n        ax_hist.set_ylim(0, 350)\n    \n    ax_image.imshow(image.pixel_array, cmap = plt.cm.gray)\n    ax_image.grid(False)\n    ax_image.axis('off')\n    \n    plt.show()\n\nsample_files_list = os.listdir(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\")[::40]\nsample_files_list.sort(key=lambda f: int(re.sub('\\D', '', f)))\nfor file in range(len(sample_files_list)):\n    sample_files_list[file] = os.path.join(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\", sample_files_list[file])\n    \nfor file in sample_files_list:\n    plot_image_histogram(pydicom.read_file(file))","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:08.500651Z","iopub.execute_input":"2021-10-08T06:13:08.500857Z","iopub.status.idle":"2021-10-08T06:13:12.943126Z","shell.execute_reply.started":"2021-10-08T06:13:08.500833Z","shell.execute_reply":"2021-10-08T06:13:12.942146Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we can see some important information about the data. \n* First, there contains whole images with no visual information whatsoever.\n* Second, when there is an image with little visual information we see the variance of the visual data is higher.\n* Finally, of the images with plentiful visual data, images that showed the tumor well skewed slightly right.\n\nAs you can see we are not done. We need more visual data features to make this data actionable. The question now is which data shows the tumor the best?\n\nWe can figure this out by acknowledging something that was done to the visual data to make the histogram. We normalized the visual data so it all fit on the histogram with a similar x-axis scale. When normalizing this data we are able to see common differences among the pictures. One difference that would be helpful is the number of lumiensent pixels above a common normal. So for example in scan 156, we can see a large mass above the common normal of the image around it. That is the tumor! So now we can mathematically find the images that best capture the tumor.","metadata":{}},{"cell_type":"markdown","source":"First let's establish a level. I'll arbitrarily put this at 2x the standard deviation of our normalized luminosity and place it on our histogram","metadata":{}},{"cell_type":"code","source":"file = \"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR/Image-156.dcm\"\ndef find_normal_level(mean, std):\n    return mean + 1.75 * std\ndef plot_image_histogram_level(image):\n    ## Get the pixels all on to a single array\n    pixels = np.array(image.pixel_array).ravel()\n    \n    ## Create a two sided figure\n    fig, (ax_image, ax_hist) = plt.subplots(1, 2, figsize = (20,4), gridspec_kw={'width_ratios': [1, 4]})\n    fig.suptitle(f'Training Data from RSNA-MICCAI Brain Tumor Radiogenomic Classification Dataset \\n scan # ({image.InstanceNumber})')\n    \n    ## Begin working on the histogram portion of the graph\n    ### Filter the zero pixels (greatly messes with the scale of the graph)\n    non_zero_pixels = np.nonzero(pixels)\n    ### Run basic statistics to normalize the pixel data (makes data presentable)\n    mean = np.mean(non_zero_pixels)\n    std = np.std(non_zero_pixels)\n    ### Normalize the pixels\n    norm_non_zero_pixels = (non_zero_pixels - mean)/std\n    mean = np.mean(norm_non_zero_pixels)\n    std = np.std(norm_non_zero_pixels)\n    ## Graph histogram\n    if len(non_zero_pixels[0]) != 0:\n        ax_hist.hist(norm_non_zero_pixels.flatten(), bins=200)\n        ax_hist.set_xlim(-5,5)\n        ax_hist.set_ylim(0, 350)\n        ax_limits = ax_hist.get_ylim()\n        ax_hist.vlines(find_normal_level(mean, std), ymin=ax_limits[0], ymax=ax_limits[1], colors='b', linestyles='dashed')\n        \n    ax_image.imshow(image.pixel_array, cmap = plt.cm.gray)\n    ax_image.grid(False)\n    ax_image.axis('off')\n    \n    plt.show()\n\nsample_files_list = os.listdir(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\")[::40]\nsample_files_list.sort(key=lambda f: int(re.sub('\\D', '', f)))\nfor file in range(len(sample_files_list)):\n    sample_files_list[file] = os.path.join(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\", sample_files_list[file])\n    \nfor file in sample_files_list:\n    plot_image_histogram_level(pydicom.read_file(file))","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:12.945361Z","iopub.execute_input":"2021-10-08T06:13:12.945740Z","iopub.status.idle":"2021-10-08T06:13:17.386706Z","shell.execute_reply.started":"2021-10-08T06:13:12.945710Z","shell.execute_reply":"2021-10-08T06:13:17.385812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Focus your attention to the y-axis, this shows the amount of data that exceeds the normal in a better light. Because the y-axis is so much larger in the vast majority of these images we can see which image is showing a tumor based on how many pixels have a luminosity above the common mean luminosity. We can see in these graphs that this is the case.","metadata":{}},{"cell_type":"markdown","source":"Now let's try to **maximize** this measure and see if we get images that appropriately match our goal. ","metadata":{}},{"cell_type":"code","source":"def calc_brilliance(image):\n    ## Get the pixels all on to a single array\n    pixels = np.array(image.pixel_array).ravel()\n    non_zero_pixels = np.array(image.pixel_array).ravel()[np.nonzero(np.array(image.pixel_array).ravel())]\n    print(image.InstanceNumber)\n    print(non_zero_pixels)\n    if len(non_zero_pixels) == 0:\n        return 0\n    elif len(non_zero_pixels) < (0.1 * len(pixels)):\n        return 0\n    ### Run basic statistics to normalize the pixel data (makes data presentable)\n    mean = np.mean(non_zero_pixels)\n    std = np.std(non_zero_pixels)\n    ### Normalize the pixels\n    norm_non_zero_pixels = (non_zero_pixels - mean)/std\n#     mean = np.mean(norm_non_zero_pixels)\n#     std = np.std(norm_non_zero_pixels)\n    return np.count_nonzero(image.pixel_array > find_normal_level(mean,std))\n\ndef top_brilliant_image(images):\n    brilliancelist = [calc_brilliance(image) for image in images]\n    pd.set_option(\"display.max_rows\", None)\n    print(pd.DataFrame(brilliancelist))\n    print(np.argsort(brilliancelist)[::-1][:10])\n    top_image = np.argsort(brilliancelist)[::-1][:10]\n    return top_image\n\nsample_files = []\n\nfor root, dirs, files in os.walk(os.path.abspath(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\")):\n    for file in files:\n        sample_files.append(os.path.join(root, file))\n\ntop_brilliant_ids = top_brilliant_image([pydicom.read_file(img) for img in sample_files])\n\nfor top_brilliant_id in top_brilliant_ids:\n    path_to_top_brilliant_image = os.path.join(\"../input/rsna-miccai-brain-tumor-radiogenomic-classification/train/00049/FLAIR\", \"Image-\" + str(top_brilliant_id) + \".dcm\")\n    plot_image_histogram_level(pydicom.read_file(path_to_top_brilliant_image))","metadata":{"execution":{"iopub.status.busy":"2021-10-08T06:13:17.388435Z","iopub.execute_input":"2021-10-08T06:13:17.388945Z","iopub.status.idle":"2021-10-08T06:13:24.967043Z","shell.execute_reply.started":"2021-10-08T06:13:17.388897Z","shell.execute_reply":"2021-10-08T06:13:24.966026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Clearly this needs more refining, but the principal more or less holds. In the top 2 of images with the highest cumulitave normalized luminosity above 2 standard deviations of said luminosity, we were able to find images with tumors over 50% of the time. When considering that there is 296 images in this folder alone. This is a great start for cleaning the data for our future model. A future refined version of this heuristic model may help cull the amount of images we have to model. We can simply take out images of less or no importance, which will help make the model more efficient and effective.","metadata":{}}]}