{"cells":[{"metadata":{},"cell_type":"markdown","source":"# DICOM metadata EDA\n\nWork in progress to analyze the metadata attributes of CT scans in DICOM files. Comments welcome.\n\nThe data was extracted in another notebook by loading every single CT scan and dumping each metadata attribute. I saved it to a CSV file to speed up (extracting takes 5-10 mins). [This kernel](https://www.kaggle.com/anarthal/dicom-metadata-extracting-attributes-to-dataframe) shows how I did it."},{"metadata":{},"cell_type":"markdown","source":"# 1. Imports\n\nGCDM is needed to load certain DICOM files. Not strictly needed here, but just in case."},{"metadata":{"_kg_hide-output":true,"trusted":true},"cell_type":"code","source":"!conda install -c conda-forge gdcm -y","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"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)\nimport seaborn as sns\nfrom matplotlib import pyplot as plt\nimport matplotlib\nimport pydicom\nimport os\nfrom os import path\nimport tqdm\nimport IPython\nimport re\nimport PIL\nimport scipy.misc\nfrom pathlib import Path\n\nsns.set()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 2. Loading the data\n\nLoad the data. Some fields are arrays but are stored as strings containing Python literal arrays."},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"dftrain = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/train.csv')\ndftest = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-progression/test.csv')\ndfmeta = pd.read_csv('/kaggle/input/osic-pulmonary-fibrosis-metadata/meta.csv', converters={\n    'ImageType': eval,\n    'ImagePositionPatient': lambda x: tuple(float(elm) for elm in eval(x)) if x != '' else np.nan\n})\ndfmeta['ImageType'] = dfmeta.ImageType.map(lambda x: tuple(x) if type(x) is list else tuple())\ndfmeta['ImagePositionPatientX'] = dfmeta['ImagePositionPatient'].map(lambda x: x[0] if type(x) is tuple else np.nan)\ndfmeta['ImagePositionPatientY'] = dfmeta['ImagePositionPatient'].map(lambda x: x[1] if type(x) is tuple else np.nan)\ndfmeta['ImagePositionPatientZ'] = dfmeta['ImagePositionPatient'].map(lambda x: x[2] if type(x) is tuple else np.nan)\n\ndcms = []\nfor root, dirs, fnames in os.walk('/kaggle/input/osic-pulmonary-fibrosis-progression/train'):\n    dcms += list(os.path.join(root, f) for f in fnames if f.endswith('.dcm'))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"TBC: explain this"},{"metadata":{"trusted":true},"cell_type":"code","source":"dfpat = dfmeta.groupby('PatientID').first()\ndfpat['SliceCount'] = dfmeta.groupby('PatientID').size()\ndfpat['Span'] = dfmeta.groupby('PatientID')['ImagePositionPatientZ'].apply(lambda x: x.max() - x.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfpat.loc['ID00165637202237320314458']","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 3. DICOM attributes\n\nList of analyzed DICOM metadata (see cells below for graphs and tables supporting my arguments):\n\n- BitsAllocated: https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280100 Number of bits per pixel allocated in the image. Always 16 (DROP).\n- BitsStored: https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280101 Number of bits per pixel stored in the image. Ranging from 12 to 16 - see what does this mean.\n- BodyPartExamined: https://dicom.innolitics.com/ciods/raw-data/general-series/00180015 Obvious meaning, always 'Chest' (DROP)\n- [Rows](https://dicom.innolitics.com/ciods/mr-image/image-pixel/00280010), [Columns](https://dicom.innolitics.com/ciods/mr-image/image-pixel/00280011): exactly what you would expect, equivalent to obj.pixel_array.shape.\n- [ConvolutionKernel](https://dicom.innolitics.com/ciods/ct-image/ct-image/00181210): ??? widely range of values available\n- [DeidentificationMethod](https://dicom.innolitics.com/ciods/ct-image/patient/00120063): how the CT scan was made anonymous. All values are identical (DROP)\n- [GantryDetectorTilt](http://dicomlookup.com/lookup.asp?sw=Tnumber&q=(0018,1120)): the angle of the gantry when the CT scan was made. It's informational and its value is constant (DROP)\n- [HighBit](https://dicom.innolitics.com/ciods/us-image/image-pixel/00280102): the specification says it should be 1 less than BitsStored. Train test agrees with that, so no extra info here (DROP).\n- [InstanceNumber](https://dicom.innolitics.com/ciods/ct-image/general-image/00200013): an identification for the current file. In our case, this matches with the file name (e.g. for a file called 42.dcm, InstanceNumber == 42). Not useful (DROP).\n- [KVP](https://dicom.innolitics.com/ciods/digital-x-ray-image/x-ray-generation/00180060): [Peak kilo-voltage](https://en.wikipedia.org/wiki/Peak_kilovoltage) of the X-ray generator. A measure of the power of the generated X-rays. More power tends to generate better images (TBC: study [this article](https://radiopaedia.org/articles/kilovoltage-peak)). Ranging from 100 to 140, mode 120.\n- [Manufacturer](https://dicom.innolitics.com/ciods/ct-image/general-equipment/00080070): the company that produced the equipment (e.g. SIEMENS, Philips...). There are several different values, although I don't see how these could be important for our purpose (DROP).\n- [ManufacturerModelName](https://dicom.innolitics.com/ciods/ct-image/general-equipment/00081090): the equipment model's name. Same consideration as above (DROP).\n- [Modality](https://www.dicomlibrary.com/dicom/modality/): describes what does the DICOM file contain. In our case, it's always CT, for Computed Tomography scan (DROP).\n- [PatientID](https://dicom.innolitics.com/ciods/ct-image/patient/00100020): self explanatory. In the training set, checked that it's always the same as the directory name the image is in (DROP).\n- [PatientPosition](https://dicom.innolitics.com/ciods/ct-image/general-series/00185100): position of the patient when the CT scan was taken. Seems important as the CT scan seems oriented differently for images with different values (see GIFs in cells below to see what I mean). Standard says it's just informative (TBC: check Patient Orientation Code Sequence)\n    - HFS (Head-First Supine): head goes first, and the patient is in supine potition (meaning that her head is looking upwards, to the ceiling).\n    - FFS (Feet-First Supine): feet go first, and the patient is in supine potition (meaning that her head is looking upwards, to the ceiling).\n    - HFP (Head-First Prone): head goes first, and the patient is in prone potition (meaning that her head is looking downwards, to the floor).\n    - FFP (Feet-First Prone): feet go first, and the patient is in prone potition (meaning that her head is looking downwards, to the floor).\n- PatientSex: seems redundant, we already have sex in the tabular metadata (DROP).\n- [PhotometricInterpretation](https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280004): tells us how to interpret the pixel values in the image. For all CT scans here we have PhotometricInterpretation == 'MONOCHROME2', which means that it's a grayscale image with lower pixel corresponding to darker colors than higher ones (DROP).\n- [PixelRepresentation](https://dicom.innolitics.com/ciods/ct-image/image-pixel/00280103): whether the image pixel values are encoded as signed or unsigned integers. I assume this is already taken into account by pydicom when reading the image (grepping the source codes yields a zound of matches). Our dataset has both signed and unsigned represented images. Too low level (DROP).\n- [StudyInstanceID](https://dicom.innolitics.com/ciods/ct-image/general-study/0020000d): identifier for the study. According to [this stack overflow answer](https://stackoverflow.com/questions/1434918/dicom-whats-the-point-of-sopinstanceuid-tag), DICOM defines a hierarchy patient -> study -> series -> instance. In our case, there is a 1-to-1 correspondence between Patient IDs and StudyInstanceID, so this gives us no extra info (DROP).\n- [SOPInstanceUID](https://dicom.innolitics.com/ciods/ct-image/sop-common/00080018): a unique identifier for each CT scan. (In DICOM, SOP means Service Object Pair. A CT scan is one of the possible type of SOPs the standard defines). Not useful (DROP).\n- [ImageType](https://dicom.innolitics.com/ciods/ct-image/general-image/00080008): information about what does the image contain. It is actually a list of 2 or more values. Patient ID00421637202311550012437 appears to be corrupted, as it has the string literal '1'. The values in here mean the following:\n    - ImageType[0] can be ORIGINAL or DERIVED. Derived images are generated from other ones, while originals are not (from what I've read, derived images entail some post-processing while primary images do not). Most of the images are ORIGINAL, but we also have a few DERIVED ones.\n    - ImageType[1] can be PRIMARY or SECONDARY. Primary images are created directly by the examination of the patient, while secondary ones are created afterwards (from what I've read, a 3D reconstruction would be secondary while a plain scan would be primary). Most of the images are PRIMARY, while a few of them are SECONDARY. It seems like the separation is a little bit blurry, so I'm not sure how much info does this give us.\n    - [ImageType[2]](https://dicom.innolitics.com/ciods/ct-image/ct-image/00080008) is CT scan-specific and may be AXIAL or LOCALIZER. I can only see AXIAL images here. However, there are other two values in our dataset: REFORMATTED and OTHER (no idea what is this about).\n    - Further values are implementation specific, very little info on them.\n    \n    \nThe following is a list of attributes I've seen in several notebooks preprocessing images, like ([this](https://www.kaggle.com/anarthal/osic-autoencoder-training/edit) and (this)[https://www.kaggle.com/gzuidhof/full-preprocessing-tutorial]):\n\n- [ImagePositionPatient](https://dicom.innolitics.com/ciods/ct-image/image-plane/00200032): an array of 3 elements indicating the coordinates (x, y, z) of the upper-left corner of the image. Most of the images have this attribute, except for 3 patients.\n    - Z axis is increasing towards the patients head (slices with greater ImagePositionPatient[2] correspond to the upper part of the lungs). This value can be used to order the slices to reconstruct the 3D volume the CT scan is representing.\n    - X and Y axis origins differ from patient to patient. In general, the X and Y origins are the same for every slice of a scan, except for two patients. TBC: study if these attributes are somehow useful.\n- [SliceLocation](https://dicom.innolitics.com/ciods/ct-image/image-plane/00201041): a single number specifying the location of the image plane with respect to an implementation-defined plane. Contains similar information to ImagePositionPatient[2], but it's not always consistent (sometimes ImagePositionPatient[2] indicates that the increasing z-axis direction is one, while SliceLocation indicates the opposite). All train images that have a value for ImagePositionPatient[2] have also a value for SliceLocation. TBC: investigate further and see if there is any relation with ImageOrientationPatient.\n- [RescaleType](https://dicom.innolitics.com/ciods/ct-image/ct-image/00281054): DICOM images are stored scaled. This attribute is supposed to indicate how to convert between the values in the pixel array (the \"stored values\") and actual meaningful values. In most cases, this attribute is either NaN or HU (meaning Hounsfield Units). Some others are US or UNSPECIFIED (synonyms). It's reasonable to assume everything is in HU.\n- [RescaleIntercept](https://dicom.innolitics.com/ciods/ct-image/ct-image/00281052) and [RescaleSlope](https://dicom.innolitics.com/ciods/ct-image/ct-image/00281053): to convert from the pixel array to HU, we have to apply a linear transformation like:  `HU = RescaleIntercept + RescaleSlope * sv`, where `sv` is the stored value. Note that these values are constant among the different slices within a single CT scan. (Actually, in the training set the only meaningful term is the intercept - the slope is always 1).\n- [SliceThickness](https://dicom.innolitics.com/ciods/ct-image/image-plane/00180050): how thick is the slice along the Z axis that this image represents? \n\n\nTBC:\n\n- ImageOrientationPatient\n- DistanceSourceToDetector\n- DistanceSourceToPatient\n- FocalSpots\n- FrameOfReferenceUID\n- GeneratorPower\n- LargestImagePixelValue\n- Modality\n- PatientOrientation\n- PixelPaddingValue\n- PixelRepresentation\n- PixelSpacing\n- PositionReferenceIndicator\n- RevolutionTime\n- RotationDirection\n- SamplesPerPixel\n- SeriesInstanceUID\n- SingleCollimationWidth\n- SliceThickness\n- SmallestImagePixelValue\n- SpacingBetweenSlices\n- SpatialResolution\n- SpecificCharacterSet\n- SpiralPitchFactor\n- TableFeedPerRotation\n- TableHeight\n- TableSpeed\n- TotalCollimationWidth\n- WindowCenter\n- WindowCenterWidthExplanation\n- WindowWidth\n- XRayTubeCurrent"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.SliceCount * dfmeta.","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.1. NAN count"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(10, 20))\nnans = dfmeta.isna().sum().sort_index()\nsns.barplot(y=nans.index, x=nans, orient='h')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.2. BitsStored"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(dfmeta.BitsStored)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.3. Image sizes"},{"metadata":{"trusted":true},"cell_type":"code","source":"sizes = dfmeta.apply(lambda x: f'{x.Rows}x{x.Columns}', axis=1)\nplt.figure(figsize=(15, 8))\nsns.countplot(sizes)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.4. KVP"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(8, 8))\nsns.countplot(dfmeta.KVP)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.5. Manufacturer"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(15, 8))\nsns.countplot(dfmeta.Manufacturer)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.ManufacturerModelName.value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.6. PatientPosition\n\nAs explained above, this tells us how the patient was positioned in the scan. Let's visualize two positions to understand it. TBC: may be worth considering rotating the images?"},{"metadata":{"trusted":true},"cell_type":"code","source":"tmp = dfmeta.groupby('PatientPosition')['PatientID'].nunique()\nsns.barplot(x=tmp.index, y=tmp.values)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's use GIFs to see it. Note: if a patient has too many slices, the GIF may freeze your browser! Note how the CT scan changes depending on PatientPosition. The GIF creation is inspired in [this notebook](https://www.kaggle.com/andradaolteanu/pulmonary-fibrosis-competition-eda-dicom-prep)."},{"metadata":{"trusted":true},"cell_type":"code","source":"def display_gif(img, duration=10000, low=-2048., high=2048.):\n    def to_pil(slc):\n        lower = low\n        upper = high\n        pixels = (255. * (slc - lower) / (upper - lower)).astype('uint8')\n        pixels = np.clip(pixels, 0, 255)\n        return PIL.Image.fromarray(pixels, 'L')\n\n    pil_imgs = [to_pil(slc) for slc in img]\n    print(f'{len(pil_imgs)} frames for this patient')\n    pil_imgs[0].save('tmp.gif', format='GIF',\n                   append_images=pil_imgs[1:],\n                   save_all=True,\n                   duration=duration//len(pil_imgs), loop=0)\n    IPython.display.display(display(IPython.display.Image('tmp.gif')))\n    \ndef read_patient_img(patient_id):\n    patient_dir = Path('/kaggle/input/osic-pulmonary-fibrosis-progression/train') / patient_id\n    slices = [pydicom.read_file(p) for p in patient_dir.glob('*.dcm')]\n    try:\n        slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n    except AttributeError:\n        pass\n    image = np.stack([s.pixel_array.astype(float) * s.RescaleSlope + s.RescaleIntercept  for s in slices])\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"selected_patients = dfmeta.groupby('PatientPosition').apply(lambda gp: gp.PatientID.value_counts().index[-1])\ndescr = {\n    'FFP': 'Feet-first prone (head down)',\n    'FFS': 'Feet-first supine (head up)',\n    'HFP': 'Head-first prone (head down)',\n    'HFS': 'Head-first supine (head up)',\n}\nfor pos in ['FFP', 'FFS']:\n    print(f'{descr[pos]}')\n    display_gif(read_patient_img(selected_patients[pos]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.7. Pixel representation"},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.PixelRepresentation.value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.8. ImageType"},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta['ImageType0'] = dfmeta.ImageType.map(lambda x: x[0] if len(x) >= 1 else np.nan)\ndfmeta['ImageType1'] = dfmeta.ImageType.map(lambda x: x[1] if len(x) >= 2 else np.nan)\ndfmeta['ImageType2'] = dfmeta.ImageType.map(lambda x: x[2:])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(dfmeta.ImageType0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.countplot(dfmeta.ImageType1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.groupby('PatientID').apply(lambda gp: gp.ImageType0.iloc[0]).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.groupby('PatientID').apply(lambda gp: gp.ImageType1.iloc[0]).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.ImageType.map(lambda x: x[:2]).value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.ImageType2.value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.9. ImagePositionPatient"},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta['ImagePositionPatientX'] = dfmeta.ImagePositionPatient.map(lambda x: x[0] if type(x) is tuple else np.nan)\ndfmeta['ImagePositionPatientY'] = dfmeta.ImagePositionPatient.map(lambda x: x[1] if type(x) is tuple else np.nan)\ndfmeta['ImagePositionPatientZ'] = dfmeta.ImagePositionPatient.map(lambda x: x[2] if type(x) is tuple else np.nan)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta[~dfmeta.ImagePositionPatient.isna()].groupby('PatientID').ImagePositionPatientY.nunique().value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.10. Rescaling (converting to Hounsfield Units)"},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta['RescaleType'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.groupby('PatientID')['RescaleIntercept'].nunique().value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta.groupby('PatientID')['RescaleSlope'].nunique().value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta['RescaleIntercept'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta['RescaleSlope'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 3.11. Voxel size"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.jointplot(data=dfpat, x='SliceThickness', y='SliceCount')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfpat[dfpat['SliceThickness'] > 7]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"patid = 'ID00229637202260254240583'\nprint(f'Patient {patid}, with SliceCount={dfpat.loc[patid, \"SliceCount\"]}, SliceThickness={dfpat.loc[patid, \"SliceThickness\"]}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dfmeta[dfmeta.PatientID.eq(patid)].ImagePositionPatient","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display_gif(read_patient_img('ID00229637202260254240583'))","execution_count":null,"outputs":[]}],"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}