{"cells":[{"metadata":{},"cell_type":"markdown","source":"When you handle medical images, most images have format called DICOM(.dcm).\n\nDICOM is global standard of medical imaging and has some unique properties compared to usual imaging formats, such as .png, .jpeg.\n\nIn this notebook, we will see what factors should we concern, and how to handle them.\n\nFirst, let's load some packages.","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"!pip install natsort\n\nimport os\nimport numpy as np\nimport SimpleITK as sitk\nimport matplotlib.pyplot as plt\nfrom natsort import natsorted","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I know pydicom is more famous than SimpleITK, yet I prefer SimpleITK.\n\nYou can use whatever package you want.\n\nNatsort is natural sort package in python.","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"path_train = os.path.join('../input/osic-pulmonary-fibrosis-progression', 'train')\nptns_train = [os.path.join(path_train, _) for _ in os.listdir(path_train)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ptns_train[:10]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Number of Slices\n\nFirst, let's see how much slices each patient have. This can be counted by reading number of files in one folder (=patient).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"tot_len = []\n\nfor ptn in ptns_train:\n    dcmlist = os.listdir(ptn)\n    tot_len.append(len(dcmlist))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.hist(tot_len, bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Oops! There is high discrepancy between patients!!!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Slice Thickness\n\nWhen you handle DICOM, especially CT or MRI data, data are stored in bulk of one shooting.\n\nTherefore, it can be considered as 3-dimensional data.\n\nHowever, computers can only handle discrete data, thus there is a concept called \"Slice Thickness\" in CT images.\n\nUsually, CT machine takes image of the patient by rotating spirally. Therefore, original CT data can be re-constructed with any slice thickness bigger than 0.6mm (Usually. This can depend on CT machine).\n\nTo see how thick train images are, we can check slice thickness metadata. \n\nThis can be done with .GetSpacing() method or .GetMetaData(\"0018|0050\") method on SimpleITK Image class.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_thickness = []\n\nfor ptn in ptns_train:\n    exdcm = [os.path.join(ptn, _) for _ in os.listdir(ptn)][0]\n    spacing = sitk.ReadImage(exdcm).GetSpacing()\n    slice_thickness.append(spacing[2])\nplt.hist(slice_thickness, bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Every patient has slice thickness with 1mm! How happy we are!!\n\nHowever, for the cross-check, let's see how far between slices are.\n\n# Slice Distance, Pixel Spacing\n\nBy subtracting absolute location of each DICOM files, we can acquire distance between slices.\n\nThis can be done with .GetMetaData(\"0020|0032\") method on SimpleITK Image class.\n\nAnd also, we can know physical length of each pixel by .GetSpacing() method.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"slice_interval = []\nspacing_list = []\n\nimport pydicom\n\nfor ptn in ptns_train:\n    try:\n        exdcm1 = natsorted([os.path.join(ptn, _) for _ in os.listdir(ptn)])[0]\n        exdcm2 = natsorted([os.path.join(ptn, _) for _ in os.listdir(ptn)])[1]\n        location1 = sitk.ReadImage(exdcm1).GetMetaData('0020|0032').split('\\\\')[2]\n        location2 = sitk.ReadImage(exdcm2).GetMetaData('0020|0032').split('\\\\')[2]\n        spacing = sitk.ReadImage(exdcm1).GetSpacing()[0]\n        spacing_list.append(spacing)\n        interval = np.abs(float(location2) - float(location1))\n        slice_interval.append(interval)\n    except:\n        print(ptn)\n    slice_thickness.append(interval)\n\nprint(len(slice_interval))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title(\"Physical size of each pixel\")\nplt.hist(spacing_list, bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title(\"Distance between slices\")\nplt.hist(slice_interval, bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Well, pixel size is very different from patient to patient, as well as distance between slices.\n\nWhy distance has discrepancy though slice thickness is 1mm?\n\nThis can happen when protocols of hospital is \"Reconstruct CT image with 1mm. Then save 1 slice of consecutive 10 slices\" for various reasons.\n\nFor example, slict distance with 20mm is that only one slice is selected of consecutive 20 slices.\n\nHappiness does not go so far...","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ptn1_path = os.path.join(path_train, 'ID00078637202199415319443')\nptn1_dcm = [os.path.join(path_train, 'ID00078637202199415319443', _) for _ in natsorted(os.listdir(ptn1_path))][0]\nnpy1 = sitk.GetArrayFromImage(sitk.ReadImage(ptn1_dcm)).squeeze()\nplt.title('Patient 1')\nplt.imshow(npy1, 'gray')\nplt.show()\n\nptn2_path = os.path.join(path_train, 'ID00128637202219474716089')\nptn2_dcm = [os.path.join(path_train, 'ID00128637202219474716089', _) for _ in natsorted(os.listdir(ptn2_path))][0]\nnpy2 = sitk.GetArrayFromImage(sitk.ReadImage(ptn2_dcm)).squeeze()\nplt.title(\"Patient 2\")\nplt.imshow(npy2, 'gray')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In the above images, there is a significant difference - Patient 1 is surrounded by a circle, yet Patient 2 is not.\n\nWhy is this important?\n\nWhen CT gantry rotates, it has round shape and gets image with circle shape.\n\nHowever, what we really see is square image, therefore what we really see is not exactly what CT machine takes.\n\nSome CT manufacturers consider area outside of the gantry as uncredible area thus masking it with values less than -1024(Hounsfield unit of air), like -3072.\n\nOther CT manufacturers consider area outside of the gantry as noise area and just fills that area with air Hounsfield unit.\n\nTherefore, if we see the min value of npy1, npy2 and draws histogram,","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"print(npy1.min(), npy2.min())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.title(\"Histogram of Patient 1\")\nplt.hist(npy1.flatten(), bins=20)\nplt.show()\n\nplt.title(\"Histogram of Patient 2\")\nplt.hist(npy2.flatten(), bins=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can see, the circled image (Patient 1) has three peaks on histogram - the left peak is outside of circle, middle peak is air area, right peak is body area\n\nHowever, the non-circled image (Patient 2) has two peaks on histogram - the left peak is air area, right peak is body area.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"less_than_3000 = []\nbigger_than_3000 = []\n\nfor ptn in ptns_train:\n    exdcm = natsorted([os.path.join(ptn, _) for _ in os.listdir(ptn)])[0]\n    if sitk.GetArrayFromImage(sitk.ReadImage(exdcm)).min()<-3000:\n        less_than_3000.append(ptn)\n    else:\n        bigger_than_3000.append(ptn)\n    \nprint(\"Number of DICOMS that have HU less than -3000:\", len(less_than_3000))\nprint(\"Else:\", len(bigger_than_3000))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we count CT images that has circle area, there seems to be 44 cases for this case","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In conclusion, DICOM is not that easy as png, jpeg image formats. There are lots to consider when preprocessing.\n\nThis is just start of DICOM handling!\n\nI hope this notebook would be helpful for begineers who start handling DICOM.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}