{"cells":[{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"# Preprocessing DICOM files to PNG with normalized colors and sizes\n\nReference : https://www.kaggle.com/donkeys/preprocessing-images-to-normalize-colors-and-sizes/notebook\n\n3D scaling and generation of 3D image arrays"},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"import pydicom\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom PIL import Image\nimport glob\nimport re\nimport os\nimport numpy as np\nimport pandas as pd\n\nimport pydicom\nimport matplotlib.pyplot as plt\nfrom tqdm.auto import tqdm\n\nimport scipy\nfrom scipy import ndimage\n\ntqdm.pandas()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"os.listdir('/kaggle/input/')","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/rsna-str-pulmonary-embolism-detection\"\n#DATA_DIR = \"data\"\n# !ls -l data/train | wc -l/","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(f\"{DATA_DIR}/train.csv\")\ndf_train.head()\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"df_test = pd.read_csv(f\"{DATA_DIR}/test.csv\")\ndf_test.head()\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"there are varying number of slices under each patient. each numbered file is a 2d-slice of a 3d image"},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true,"collapsed":true},"cell_type":"code","source":"# !ls /kaggle/input/rsna-str-pulmonary-embolism-detection/train/0003b3d648eb/d2b2960c2bbf\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"# Example Patient\n\nAs an example, let's take a look at one patients data and a few images. All the slices are under the patient ID directory, and can be sorted in numerical order to have the slices in the order they appear in the 3D construct. I viewed them as an animation in Aseprite when playing a bit and cheking pixel values. Would make nice gifs to animate if anyone wants to try."},{"metadata":{"trusted":true},"cell_type":"code","source":"\nTRAIN_DIR = \"/kaggle/input/rsna-str-pulmonary-embolism-detection/train\"\nflist1 = os.listdir(TRAIN_DIR)\ndicom_list = []\n\nfor i in range(len(flist1)):\n    temp_files = len(glob.glob(f'{TRAIN_DIR}/{flist1[i]}/*/*.dcm'))\n    dicom_list.append(temp_files)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"flist1[i]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dicom_list[0:10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(max(dicom_list), min(dicom_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nplt.plot(dicom_list)","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"files = glob.glob(f\"{DATA_DIR}/train/0003b3d648eb/*/*.dcm\")\n# files","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n# files","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"## Patient metadata\n\nThe DICOM file contains both the image data as well as metadata about the system that genrated it, and how the image is formatted. This metadata also contains the key values to normalize the values to same color space, at least close to. Probably some values helpful for rotation etc as well as I explored in my other EDA kernel.\n\nLook at the metadata for an image or two:"},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"import types\n\nds = pydicom.dcmread(f'{DATA_DIR}/train/0003b3d648eb/d2b2960c2bbf/2cb4facfafee.dcm')\nattrs = dir(ds)\nfor attr in attrs:\n    if attr.startswith(\"_\"):\n        continue\n    if attr == \"PixelData\" or attr == \"pixel_array\":\n        continue\n    var_type = type(getattr(ds,attr))\n    if var_type == types.MethodType:\n        continue\n    print(f\"{attr}: {getattr(ds,attr)}\")","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"# Rescaling algorithm\n\nThe following is an algorithm I pieced together from some Stackover flow posts, and some mentions in the DICOM specs.\n\n- [Rescale type attribute](https://dicom.innolitics.com/ciods/ct-image/ct-image/00281054)\n- [Window Center](https://dicom.innolitics.com/ciods/ct-image/voi-lut/00281050)\n- [RescaleIntercept](https://dicom.innolitics.com/ciods/ct-image/ct-image/00281052)\n\nThere seem to be quite a few StackOverflow posts discussing various tricks, most of which don't really work across all the images. I tried by running them and then eyeballing some images for all the results. WHat I have here worked best, and seemed to make sense. If you know better, let me know.\n\nSome of the more helpful StackOverflow posts:\n\nhttps://stackoverflow.com/questions/10193971/rescale-slope-and-rescale-intercept\nhttps://stackoverflow.com/questions/8756096/window-width-and-center-calculation-of-dicom-image/8765366#8765366"},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"import math\nfrom collections.abc import Iterable\nimport pydicom\n\ndef show_img(img_path, colormap = None, extra_brightness=0):\n    ds = pydicom.dcmread(img_path)\n    shape = ds.pixel_array.shape\n    target = 255\n\n    # Convert to float to avoid overflow or underflow losses.\n    image_2d = ds.pixel_array.astype(float)\n    img_data = image_2d\n    print(f\"data min: {img_data.min()}, max: {img_data.max()}\")\n    print(f\"window center: {ds.WindowCenter}, rescale intercept: {ds.RescaleIntercept}\")\n    multival = isinstance(ds.WindowCenter, Iterable)\n    if multival:\n        scale_center = -ds.WindowCenter[0]\n    else:\n        scale_center = -ds.WindowCenter\n    intercept = scale_center+ds.RescaleIntercept+extra_brightness\n    print(f\"final intercept: {intercept}\")\n    image_2d += intercept\n    print(f\"after applying intercept, min: {image_2d.min()}, max: {image_2d.max()}\")\n\n    # Rescaling grey scale between 0-255\n    image_2d_scaled = (np.maximum(image_2d,0) / image_2d.max()) * 255.0\n    print(f\"after scaling to 0-255, min: {image_2d_scaled.min()}, max: {image_2d_scaled.max()}\")\n\n    # Convert to uint\n    image_2d_scaled = np.uint8(image_2d_scaled)\n\n    plt.figure(figsize=(12,8))\n    plt.imshow(image_2d_scaled, cmap=colormap)\n    plt.show()\n\n#show_img(f'{DATA_DIR}/train/ID00011637202177653955184/1.dcm', colormap=plt.cm.bone) <-image 0 below\nshow_img(f'{DATA_DIR}/train/0003b3d648eb/d2b2960c2bbf/2cb4facfafee.dcm', colormap=plt.cm.bone) #image 1\nshow_img(files[0], colormap=plt.cm.bone) #image 2\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"import types\n\nds = pydicom.dcmread(files[0])\nattrs = dir(ds)\nfor attr in attrs:\n    if attr.startswith(\"_\"):\n        continue\n    if attr == \"PixelData\" or attr == \"pixel_array\":\n        #skip printing the long arrays as they will just spam the output too much with hex code\n        continue\n    var_type = type(getattr(ds,attr))\n    if var_type == types.MethodType:\n        continue\n    print(f\"{attr}: {getattr(ds,attr)}\")","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"## Different colormaps and extra brightness adjustments\n\nSomeone asked if the colors could be brighter. Lets try a different colormap:"},{"metadata":{"trusted":true},"cell_type":"code","source":"#show_img(f'{DATA_DIR}/train/ID00011637202177653955184/1.dcm', colormap=None)\nshow_img(f'{DATA_DIR}/train/0003b3d648eb/d2b2960c2bbf/2cb4facfafee.dcm', colormap=None)\nshow_img(files[0], colormap=None)","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"And bone colormap with some extra brightness on top:"},{"metadata":{"trusted":true},"cell_type":"code","source":"#show_img(f'{DATA_DIR}/train/ID00011637202177653955184/1.dcm', colormap=plt.cm.bone, extra_brightness=1000)\nshow_img(f'{DATA_DIR}/train/0003b3d648eb/d2b2960c2bbf/2cb4facfafee.dcm', colormap=plt.cm.bone, extra_brightness=1000)\nshow_img(files[0], colormap=plt.cm.bone, extra_brightness=1000)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"My guess is, a CNN will \"see\" the pixel values fine either way, but it would be interesting to experiment. Easier perhaps to look at some of the images with human eye if brightness is a bit higher."},{"metadata":{},"cell_type":"markdown","source":"# Scaling Z-Axis"},{"metadata":{"trusted":true},"cell_type":"code","source":"def resample_z(image, scan, z, new_spacing=[1, 1, 1]):\n    #This is probably a bit more complex implementation than is needed, since we just want to scale Z to constant\n    #But its what I ended up and it works, so I left it as is..\n    # Determine current pixel spacing\n    spacing = np.array([scan[0][\"SliceThickness\"]] + [scan[0][\"PixelSpacing_0\"], scan[0][\"PixelSpacing_1\"]], dtype=np.float32)\n\n    resize_factor = spacing / new_spacing\n    new_real_shape = image.shape * resize_factor\n    new_shape = np.round(new_real_shape)\n    real_resize_factor = new_shape / image.shape\n    new_spacing = spacing / real_resize_factor\n    #factor is (z,y,x)\n    real_resize_factor = (z/image.shape[0],1,1)\n\n    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor, mode='nearest')\n\n    return image, new_spacing","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def dump_scaled(np_arrays, df_name, patient_id, save_png=True, save3d=False):\n    full_3d = []\n    full_3d_256 = []\n    full_3d_192 = []\n    full_3d_128 = []\n    for idx, img in enumerate(np_arrays):\n        file = f\"img{idx}\"\n        processed_filename = os.path.basename(file)\n        processed_filename = processed_filename.split(\".\")[0]+\".png\"\n        processed_dir = f\"scaled_png/{df_name}/{patient_id}\"\n        os.makedirs(processed_dir, exist_ok=True)\n        processed_path = f\"{processed_dir}/{processed_filename}\"\n\n        im = Image.fromarray(img)\n#        im = detect_border(im, file)\n        if save_png:\n            im.save(processed_path)\n        full_3d.append(np.array(im))\n#         im256 = im.resize((256,256))\n#         im192 = im.resize((192,192))\n#         im128 = im.resize((128,128))\n#         full_3d_256.append(np.array(im256))\n#         full_3d_192.append(np.array(im192))\n#         full_3d_128.append(np.array(im128))\n    full_3d = np.array(full_3d)\n#     full_3d_256 = np.array(full_3d_256)\n#     full_3d_192 = np.array(full_3d_192)\n#     full_3d_128 = np.array(full_3d_128)\n    if save3d:\n        np.save(f\"scaled_png/{df_name}/{patient_id}/full_3d_512\", full_3d)\n        #np.save(f\"scaled_png/{df_name}/{patient_id}/full_3d_256\", full_3d_256)\n        #np.save(f\"scaled_png/{df_name}/{patient_id}/full_3d_192\", full_3d_192)\n        #np.save(f\"scaled_png/{df_name}/{patient_id}/full_3d_128\", full_3d_128)\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"# Processing all files"},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"Now, turning all these into functions to process all the images in a directory:"},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"import os\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\n\n#this function is similar to the one used higher above to illustrate the intercept scaling on color values\n#this is just tailored to be run on all files at once\ndef scale_to_png(ds, file, df_name, patient_id):\n    image_2d = ds.pixel_array.astype(float)\n    multival = isinstance(ds.WindowCenter, Iterable)\n    if multival:\n        scale_center = -ds.WindowCenter[0]\n    else:\n        scale_center = -ds.WindowCenter\n    intercept = scale_center+ds.RescaleIntercept\n    image_2d += intercept\n\n    # Rescaling grey scale between 0-255\n    image_2d_scaled = (np.maximum(image_2d,0) / image_2d.max()) * 255.0\n    img_data = image_2d_scaled\n\n    image_2d_scaled = np.uint8(image_2d_scaled)\n\n    processed_filename = os.path.basename(file)\n    processed_filename = processed_filename.split(\".\")[0]+\".png\"\n    processed_dir = f\"png/{df_name}/{patient_id}\"\n    os.makedirs(processed_dir, exist_ok=True)\n    processed_path = f\"{processed_dir}/{processed_filename}\"\n\n    im = Image.fromarray(image_2d_scaled)\n    im = detect_border(im, file)\n    im.save(processed_path)\n\n    shape = ds.pixel_array.shape\n\n    return im, processed_filename\n\nfrom PIL import Image, ImageChops, ImageOps\n\n#most common filesize seemed to be 512x512, so converting all files to that size\ndef resize_to_512(im: Image, image_name: str):\n    width, height = im.size\n    if width == 512 and height == 512:\n        return im\n    if width != height:\n        if width < height:\n            pad_w = height - width\n            pad_w /= 2\n            pad_h = 0\n        else:\n            pad_w = 0\n            pad_h = width-height\n            pad_h /= 2\n        padding = (pad_w, pad_h, pad_w, pad_h )\n        print(\"WARN: resizing image {image_name}\")\n        #we should not come here but if we do, this should resize to square\n        ImageOps.expand(im, padding, Image.ANTIALIAS)\n\n    im2 = im.resize((512, 512))\n    return im2\n\n#some images in the dataset have a grayish border around the actual image data, \n#with the actual image data being 512x512 size. \n#this removes the grey border and keeps the actual data\ndef detect_border(im, filepath):\n    bg = Image.new(im.mode, im.size, im.getpixel((0,0)))\n    diff = ImageChops.difference(im, bg)\n    diff = ImageChops.add(diff, diff, 1.0, -0)\n    bbox = diff.getbbox()\n    if bbox:\n        im2 = im.crop(bbox)\n        width, height = im2.size\n        if width == height:\n            im = im2\n        else:\n            im = im\n    im = resize_to_512(im, filepath)\n    return im\n\ndef process_dataset(base_dir, df_data, df_name):\n    ds_files = []\n    attribute_names = set()\n\n    #patients = df_data[\"Patient\"].unique()\n    patients = df_data[\"StudyInstanceUID\"].unique()\n    processed_patients = set()\n    num_patients = patients.shape[0]\n    for idx, patient_id in tqdm(enumerate(patients), total=num_patients):\n        print(f\"processing patient: {patient_id}\")\n        if patient_id in processed_patients:\n            continue\n        processed_patients.add(patient_id)\n        #the following is just to skip the first 5 images on Kaggle, as it does not have GDCM installed\n        if len(processed_patients) < 5:\n            continue\n        patient_dir = f'{base_dir}/{patient_id}/'\n        files = glob.glob(f\"{patient_dir}/*/*.dcm\")\n        files.sort(key=lambda f: int(re.sub('\\D', '', f)))\n        ds_size_before = len(ds_files)\n        patient_images = []\n        patient_meta = []\n        for file in files:\n            ds = pydicom.dcmread(file)\n            attrs = dir(ds)\n            metadata = {}\n            metadata[\"PatientId\"] = patient_id\n            metadata[\"file_idx\"] = idx\n            ds_files.append(metadata)\n            im, filename = scale_to_png(ds, file, df_name, patient_id)\n            metadata[\"filename\"] = filename\n            patient_images.append(np.array(im))\n            patient_meta.append(metadata)\n\n            for attr in attrs:\n                if attr.startswith(\"_\"):\n                    continue\n                if attr == \"PixelData\" or attr == \"pixel_array\" or attr == \"fromkeys\" or attr == \"copy\":\n                    continue\n                var_type = type(getattr(ds,attr))\n                if var_type == types.MethodType:\n                    continue\n                if attr not in attribute_names:\n                    print(f\"{attr}: {var_type}\")\n                    attribute_names.add(attr)\n                value = getattr(ds,attr)\n                if type(value) is list or type(value) is pydicom.multival.MultiValue:\n                    for sub_idx, sub_value in enumerate(value):\n                        metadata[f\"{attr}_{sub_idx}\"] = f\"{sub_value}\"\n                else:\n                    metadata[attr] = f\"{value}\"\n            del ds\n            del attrs\n        np_images = np.array(patient_images)\n        resampled_pixels_64, spacing = resample_z(np_images, patient_meta, 64)\n        resampled_pixels_200, spacing = resample_z(np_images, patient_meta, 200)\n        #the PNG files are not really needed outside checking the code works, since the 3D arrays contain it all\n        #but I leave it here so we can see the results\n        dump_scaled(resampled_pixels_64, df_name+\"_64\", patient_id, save_png=True, save3d=True)\n        dump_scaled(resampled_pixels_200, df_name+\"_200\", patient_id, save_png=True, save3d=True)\n        ds_size_after = len(ds_files)\n        ds_diff = ds_size_after - ds_size_before\n        #this is here just to cap the number of files processed in Kaggle\n        if len(processed_patients) > 10:\n            break\n    return ds_files\n\nbase_dir = f'{DATA_DIR}/train/'\nds_train_files = process_dataset(base_dir, df_train, \"train\")\n\nbase_dir = f'{DATA_DIR}/test/'\nds_test_files = process_dataset(base_dir, df_test, \"test\")\n\n#NOTE: below all different attributes and their types found in the dataset will be printed as well","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Did we generate files? "},{"metadata":{"trusted":true},"cell_type":"code","source":"# !ls png/train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#also check the 3D arrays. here we generated ones with 20 and 30 images, so should see separate dirs for each\n# !ls scaled_png/","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# !ls scaled_png/train_30/","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%% md\n"}},"cell_type":"markdown","source":"Above processing functions also collected all the image metadata into two lists. Here we can convert those lists into dataframes and write those dataframes to disk as CSV files. This way no need to go through all the files another time on Kaggle."},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"# df_test_file_meta = pd.DataFrame(ds_test_files)\n# df_test_file_meta.to_csv(\"meta_files_test.csv\")\n# df_test_file_meta.head()\n","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"# df_train_file_meta = pd.DataFrame(ds_train_files)\n# df_train_file_meta.to_csv(\"meta_files_train.csv\")\n# df_train_file_meta.head()","execution_count":null,"outputs":[]},{"metadata":{"pycharm":{"name":"#%%\n"},"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"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}