{"cells":[{"metadata":{},"cell_type":"markdown","source":"<H1>Dicom Dataset EDA</H1>\n\nThe images in this contest are supplied in a number of ways including jpegs, tfRecords, and Dicom datasets. But as this contest is posed as a medical imaging problem, one has to wonder if it wouldn't be best to use the medical images themselves.  I have worked with dicom images once before as they were the only data provided for the 2017 Data Science Bowl contest on Kaggle.\n\nThe EDA below is a brief exploration of these images and using them with the pydicom library which you can install easily with \"pip install pydicom\"","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install p_tqdm","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# ----------------------------------------\n# imports\n# ----------------------------------------\nimport os\nimport pydicom\nimport cv2\nfrom PIL import Image\nfrom matplotlib import pyplot as plt\nfrom tqdm.notebook import tqdm\n# uncomment if you want to run the dicom data collection code\nfrom multiprocessing import pool, cpu_count\nfrom multiprocessing.dummy import Pool as ThreadPool\nfrom p_tqdm import p_map #pip install p_tqdm\nimport pandas as pd\nimport numpy as np","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ----------------------------------------\n# constants\n# ----------------------------------------\nroot_path = '/kaggle/input/siim-isic-melanoma-classification'\ndata_path = root_path\noutput_path = '/kaggle/working'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H2>Dicom Metadata Exploration</H2>\nDicom datasets come with a variety of metadata items embedded within them.  Dicom is a standard and each metadata item has a specific meaning an purpose.  Here I'll summarize a few, but you can look up anything you feel you need to know for this challenge at <A HREF=\"https://dicom.innolitics.com/ciods/segmentation/general-image/00080008\">this website.</A>  As you check this out note the right hand column which is of the form (xxxx, xxxx).  Make sure to use these references to make sure you are looking at the right definition as the text descriptions are reused across different image types.\n<br>\n<br>\nThe <B>Image Type</B> shown here indicates that this is an image whose pixel values have been derived in some manner from the pixel value of one or more other images.  Also, this was not created as a direct result of a patient examination. Instead this image is a SECONDARY image; an image created after the initial patient examination.\n<br>\n<br>\nThe <B>Modality</B> is a code indicating the equipment used to acquire the image.  In this case, XC means that the image was acquired by external camera.\n<br>\n<br>\nOne nuance to note is the interaction between <B>Samples Per Pixel</B> and <B>Photometric Interpretation</B>.  Samples Per Pixel is the number of separate planes in this image.  So no real surprise its a 3 plane image (think RGB, YUV, YCbCr).  But note that the Photometric Interpretation is  YBR_FULL_422 which means that this image is a YCbCr representation but Cb and Cr values are sampled horizontally at half the Y rate and as a result there are half as many Cb and Cr values as Y values.  Therefore, the Samples per Pixel describes the nominal number of channels (i.e., 3), and does not reflect that two chrominance samples are shared between four luminance samples.\n<br>\n<br>\nMany of the data elements are self explanatory.  So <B>Rows</B> and <B>Columns</B> are just image dimensions, and of course the image itself is called <B>Pixel Data</B>.\n<br>\n<br>\nI'll leave the rest for you to research yourself.  If you find something useful or helpful to us all, please post it!","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# ------------------------------------------------------\n# get a list of all of the training dicom files\n# ------------------------------------------------------\ndcom_list = [f for f in os.listdir(f'{data_path}/train')]\ndcom_test_list = [f for f in os.listdir(f'{data_path}/test')]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ------------------------------------------------------\n# read in a dicom image and all of its metadata\n# ------------------------------------------------------\ndataset = pydicom.dcmread(f'{data_path}/train/{dcom_list[0]}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# ------------------------------------------------------\n# print out a full listing\n# ------------------------------------------------------\nfor element in dataset:\n    print(element)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H2>Metadata Statistics</H2>\nI know that the competition organizers supplied us with a CSV to go along with the JPEG images. But this file only includes a few of the items that are available to us in the Dicom images.  Lets have a look at a broader set of metadata elements.  \n<br>\nNote that the following code block may take a little time as we have to open every Dicom dataset, extract the data elements we want and build a list.  I've saved out a csv if you want to skip running this.  That said, I did parallelize so it runs as quickly as possible.  Your mileage may vary.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# -------------------------------------------------------------------------------------------\n# helper function to open a dicom image, extract some meta data and return a list of values\n# this works with the thread pooler below to reduce wait time\n# to get both test and train you will have to switch the data_type argument to the desired\n# folder\n# -------------------------------------------------------------------------------------------\ndef collect_dcom_data(dcom_file, data_path=data_path):\n    ds = pydicom.dcmread(f'{data_path}/train/{dcom_file}')\n    return [ds.PatientID, ds.PatientSex, ds.PatientAge, ds.BodyPartExamined, ds.InstitutionName, \n            ds.ImageType, ds.Modality, ds.PhotometricInterpretation, ds.Rows, ds.Columns]\ndef collect_dcom_data_test(dcom_file, data_path=data_path):\n    ds = pydicom.dcmread(f'{data_path}/test/{dcom_file}')\n    return [ds.PatientID, ds.PatientSex, ds.PatientAge, ds.BodyPartExamined, ds.InstitutionName, \n            ds.ImageType, ds.Modality, ds.PhotometricInterpretation, ds.Rows, ds.Columns]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# -------------------------------------------------------------------------------------------\n# Multithreaded extract of meta data\n# I was using a GCP instance which only provides vCPUs (one thread per CPU)\n# your mileage will vary\n#\n# I've saved csv's of the results so you just skip this part if you want\n# -------------------------------------------------------------------------------------------\npool = ThreadPool(cpu_count())  \nmeta_data = p_map(collect_dcom_data, dcom_list)   \ndf = pd.DataFrame(meta_data, columns=['PatientID', 'PatientSex', 'PatientAge', 'BodyPartExamined', 'InstitutionName', \n                                      'ImageType', 'Modality', 'PhotometricInterpretation', \n                                      'Rows', 'Columns'])\ndf.to_csv(output_path + '/' + 'metadata.csv', index=False)\n\n\nmeta_data = p_map(collect_dcom_data_test, dcom_test_list)   \ndf_test = pd.DataFrame(meta_data, columns=['PatientID', 'PatientSex', 'PatientAge', 'BodyPartExamined', 'InstitutionName',\n                                           'ImageType', 'Modality', 'PhotometricInterpretation', \n                                           'Rows', 'Columns'])\ndf_test.to_csv(output_path + '/' + 'metadata_test.csv', index=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv(output_path + '/' + 'metadata.csv')\ndf_test = pd.read_csv(output_path + '/' + 'metadata_test.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Modality</H3>\nAll of the images are \"XC\" meaning they were all taken via an external camera.  This is true of both the train and test data.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# create the first chart\ndf_mods = df.Modality.value_counts()\nax1 = df_mods.plot.pie(ax=ax1)\n\n# create the second chart\ndf_mods = df_test.Modality.value_counts()\nax2 = df_mods.plot.pie(ax=ax2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>ImageType</H3>\nAll of the images contain pixel values have been derived in some manner from the pixel value of one or more other images.  This is also true for both train and test.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# create the first chart\ndf_type = df.ImageType.value_counts()\nax1 = df_type.plot.pie(ax=ax1, labels=None)\n\n# create the second chart\ndf_type = df_test.ImageType.value_counts()\nax2 = df_type.plot.pie(ax=ax2, labels=None)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Image Width</H3>\nThere are 88 combinations of height and width for these images.  The chart below shows the top 10 width and the related image counts with those widths.  The distribution between train and test are similar, but given the resizes that will be necessary, there are some differences to note.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# create the first chart\ndf_cols = df.Columns.value_counts()\ndf_cols[:10].plot.bar(ax=ax1)\n\n# create the second chart\ndf_cols = df_test.Columns.value_counts()\ndf_cols[:10].plot.bar(ax=ax2)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Image Height</H3>\nThe chart below shows the top 10 heights and the related image counts.  Same resize concerns are relevant here.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# create the first chart\ndf_rows = df.Rows.value_counts()\ndf_rows[:10].plot.bar(ax=ax1)\n\n# create the second chart\ndf_rows = df_test.Rows.value_counts()\ndf_rows[:10].plot.bar(ax=ax2)\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Gender Distribution</H3>\nThere are slightly more men than women in the training data and even more imbalance in the test set.  There are also a few records for which gender is not known in the training set. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# charts\nax1 = df.PatientSex.value_counts().plot.pie(ax=ax1)\nax2 = df_test.PatientSex.value_counts().plot.pie(ax=ax2)\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Age Distribution</H3>\nThe age distributions are also reasonably similar although the test set is a little notchy towards the right tail. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# set up the fig\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))\nax1.set_title('Train')\nax2.set_title('Test')\n\n# create the first chart\nages = df.PatientAge.unique()\nages.sort()\nax1 = df.PatientAge.value_counts().reindex(ages.tolist()).plot.bar(ax=ax1)\n\n# create the second chart\nages = df_test.PatientAge.unique()\nages.sort()\nax2 = df_test.PatientAge.value_counts().reindex(ages.tolist()).plot.bar(ax=ax2)\n\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H2>Dicom Images</H2>\nAs noted above, the Dicom images are natively a modified version of the YCrCb color space. Lets have a look at some of the differences between the JPEGs and the Dicom images.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# you can display a dcom image easily with matplotlib\nplt.imshow(dataset.pixel_array)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get a list of all of the training jpeg files\njpeg_list = [f for f in os.listdir(f'{data_path}/jpeg/train')]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# -----------------------------------------------\n# Displays differences for a single image\n# -----------------------------------------------\ndef display_dif_compare(image_name, nrows=4, ncols=3, figsize=15):\n    \n    # clean file name and open image in variety of styles\n    file_name, _ = os.path.splitext(image_name)\n    ds = pydicom.dcmread(f'{data_path}/train/{file_name}.dcm')\n\n    dcom_image = ds.pixel_array\n    jpeg = cv2.imread(f'{data_path}/jpeg/train/{file_name}.jpg', cv2.IMREAD_UNCHANGED)\n    #jpeg_rgb = jpeg\n    jpeg_rgb = cv2.cvtColor(jpeg, cv2.COLOR_BGR2RGB)\n    jpeg_ycc = cv2.cvtColor(jpeg, cv2.COLOR_BGR2YCR_CB)\n    dcom_rgb = cv2.cvtColor(dcom_image, cv2.COLOR_YCR_CB2RGB)\n\n    jpeg_ycc_pil = Image.open(f'{data_path}/jpeg/train/{file_name}.jpg')\n    jpeg_ycc_pil.draft('YCbCr', None)\n    jpeg_ycc_pil.load() \n\n    # plot the figures (apologies for the brute force inelegance of this block)\n    fig, axarr = plt.subplots(nrows=nrows, ncols=ncols, figsize=(figsize, figsize))\n\n    axarr[0, 0].set_title('Raw Dicom')\n    axarr[0, 0].imshow(dcom_image)\n    axarr[0, 1].set_title('JPEG Raw via cv2')\n    axarr[0, 1].imshow(jpeg_rgb)\n    axarr[0, 2].set_title('Dif: ' + str(np.sum(dcom_image - jpeg_rgb)//1000000) + ' x 10^5')\n    axarr[0, 2].imshow(abs(dcom_image - jpeg_rgb))\n    axarr[0, 0].set(ylabel = file_name)\n    \n    axarr[1, 0].set_title('Raw Dicom')\n    axarr[1, 0].imshow(dcom_image)\n    axarr[1, 1].set_title('JPEG Raw->YCrCb via cv2')\n    axarr[1, 1].imshow(jpeg_ycc)\n    axarr[1, 2].set_title('Dif: ' + str(np.sum(dcom_image - jpeg_ycc)//1000000) + ' x 10^5')\n    axarr[1, 2].imshow(abs(dcom_image - jpeg_ycc))\n    axarr[1, 0].set(ylabel = file_name)\n    \n    axarr[2, 0].set_title('Dicom Raw->RGB via cv2')\n    axarr[2, 0].imshow(dcom_rgb)\n    axarr[2, 1].set_title('JPEG Raw via cv2')\n    axarr[2, 1].imshow(jpeg_rgb)\n    axarr[2, 2].set_title('Dif: ' + str(np.sum(dcom_rgb - jpeg_rgb)//1000000) + ' x 10^5')\n    axarr[2, 2].imshow(abs(dcom_rgb - jpeg_rgb))\n    axarr[2, 0].set(ylabel = file_name)\n    \n    axarr[3, 0].set_title('Dicom Raw')\n    axarr[3, 0].imshow(dcom_image)\n    axarr[3, 1].set_title('JPEG Raw-YCrCb via PIL \"draft\"')\n    axarr[3, 1].imshow(jpeg_ycc_pil)\n    axarr[3, 2].set_title('Dif: ' + str(np.sum(dcom_image - jpeg_ycc_pil)//1000000) + ' x 10^5')\n    axarr[3, 2].imshow(abs(dcom_image - jpeg_ycc_pil))\n    axarr[3, 0].set(ylabel = file_name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Image Handling Will Be Important</H3>\nHere we compare a single image between the Dicom pixel array and the jpeg image.  Depending upon how you handle the opening and conversion of each image format, you will get various differences that may effect your model.  Recall that the Raw Dicom image is a modified YCbCr image, while the JPEG when opened by cv2 is opened in BGR format (not RGB).  Note also that if you open the image via cv2 in the normal imread fashion.  The image will be manipulated through a number of steps (see <A HREF=\"https://www.graphicsmill.com/docs/gm5/UnderstandingofJPEGEncodingParameters.htm\">here</A>. ","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"image_name = np.random.choice(jpeg_list)\ndisplay_dif_compare(image_name)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>RGB Differences</H3>\nHere you can see the differences between the two images in each color channel.  Its kind of hard to imagine this will have no impact on model results.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"file_name, _ = os.path.splitext(image_name)\nds = pydicom.dcmread(f'{data_path}/train/{file_name}.dcm')\nimage1 = cv2.cvtColor(ds.pixel_array, cv2.COLOR_YCR_CB2RGB)\nimage2 = cv2.imread(f'{data_path}/jpeg/train/{image_name}', cv2.IMREAD_UNCHANGED)\n\n# tuple to select colors of each channel line\ncolors = (\"r\", \"g\", \"b\")\nchannel_ids = (0, 1, 2)\n\nplt.figure(figsize=(20,5))\n\n# create the histogram plot, with three lines, one for\n# each color for each image\nfor channel_id, c in zip(channel_ids, colors):\n    histogram1, bin_edges = np.histogram(image1[:, :, channel_id], bins=256, range=(0, 256))\n\nfor channel_id, c in zip(channel_ids, colors):\n    histogram2, bin_edges = np.histogram(image2[:, :, channel_id], bins=256, range=(0, 256))\n\n\n# plot the reds for each image\nax1 = plt.subplot(1,3,1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='r', label='Dicom-RGB via cv2')\nplt.plot(bin_edges[0:-1], histogram2, color='r', linestyle='dashed', label='JPEG Raw via cv2')\nplt.legend()\nplt.xlabel(\"Color value\")\nplt.ylabel(\"Pixels\")\n\n# plot the greens for each image\nax2 = plt.subplot(1,3,2, sharey=ax1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='g', label='Dicom-RGB via cv2')\nplt.plot(bin_edges[0:-1], histogram2, color='g', linestyle='dashed', label='JPEG Raw via cv2')\nplt.legend()\nplt.xlabel(\"Color value\")\nplt.ylabel(\"Pixels\")\n\n# plot the greens for each image\nax3 = plt.subplot(1,3,3, sharey=ax1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='b', label='Dicom-RGB via cv2')\nplt.plot(bin_edges[0:-1], histogram2, color='b', linestyle='dashed', label='JPEG Raw via cv2')\nplt.legend()\nplt.xlabel(\"Color value\")\nplt.ylabel(\"Pixels\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Luminance and Chrominance Differences</H3>\nHere you can see the differences between the two images in the Y, Cb, and Cr channels.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"image1 = ds.pixel_array\nimage2 = cv2.cvtColor(image2, cv2.COLOR_BGR2YCR_CB)\n\n# tuple to select colors of each channel line\ncolors = (\"y\", \"Cb\", \"Cr\")\nchannel_ids = (0, 1, 2)\n\nplt.figure(figsize=(20,5))\n\n# create the histogram plot, with three lines, one for\n# each color for each image\nfor channel_id, c in zip(channel_ids, colors):\n    histogram1, bin_edges = np.histogram(image1[:, :, channel_id], bins=256, range=(0, 256))\n\nfor channel_id, c in zip(channel_ids, colors):\n    histogram2, bin_edges = np.histogram(image2[:, :, channel_id], bins=256, range=(0, 256))\n\n\n# plot the reds for each image\nax1 = plt.subplot(1,3,1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='y', label='Dicom Raw')\nplt.plot(bin_edges[0:-1], histogram2, color='y', linestyle='dashed', label='JPEG Raw-YCrCb via PIL \"draft\"')\nplt.legend()\nplt.xlabel(\"Luminance\")\nplt.ylabel(\"Pixels\")\n\n# plot the greens for each image\nax2 = plt.subplot(1,3,2, sharey=ax1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='b', label='Dicom Raw')\nplt.plot(bin_edges[0:-1], histogram2, color='b', linestyle='dashed', label='JPEG Raw-YCrCb via PIL \"draft\"')\nplt.legend()\nplt.xlabel(\"Cb\")\nplt.ylabel(\"Pixels\")\n\n# plot the greens for each image\nax3 = plt.subplot(1,3,3, sharey=ax1)\nplt.xlim([0, 256])\nplt.plot(bin_edges[0:-1], histogram1, color='r', label='Dicom Raw')\nplt.plot(bin_edges[0:-1], histogram2, color='r', linestyle='dashed', label='JPEG Raw-YCrCb via PIL \"draft\"')\nplt.legend()\nplt.xlabel(\"Cr\")\nplt.ylabel(\"Pixels\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"<H3>Conclusion</H3>\nHopefully you find this information on the Dicom format and comparisons to the JPEG format helpful.  If you do, please give me an upvote!  Thanks.","execution_count":null}],"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":4,"nbformat_minor":4}