{"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":"# Compare different DICOM conversions\n\nDifferent methods of converting DICOM to images are compared.\n\n### There is a **significant** difference between images that do and do not require GDCM. \n\n### There is additionally a significant difference between lossless images that have the exact same TransferSyntaxUID name.\n\nImages with a TransferSyntaxUID that includes 'JPEG Lossless' may need to be converted differently from 'Implicit'. The function dcm_apply_windows is the most consistent between the three. By playing with the window and level, the bone is better seen. The value from the dicom metadata is overridden by the preset values. \n\nWould love to hear ideas from the community about ways to address this.","metadata":{}},{"cell_type":"code","source":"! pip install pylibjpeg -q\n! pip install python-gdcm -q\n! pip install pylibjpeg-libjpeg -q","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:18:16.278365Z","iopub.execute_input":"2022-09-16T21:18:16.279626Z","iopub.status.idle":"2022-09-16T21:18:53.775516Z","shell.execute_reply.started":"2022-09-16T21:18:16.279525Z","shell.execute_reply":"2022-09-16T21:18:53.774428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport numpy as np\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\n\nfrom PIL import Image\n\nimport matplotlib.image as mpimg\nimport cv2","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:18:53.777761Z","iopub.execute_input":"2022-09-16T21:18:53.778192Z","iopub.status.idle":"2022-09-16T21:18:54.230784Z","shell.execute_reply.started":"2022-09-16T21:18:53.778149Z","shell.execute_reply":"2022-09-16T21:18:54.229728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"root = Path('../input/rsna-2022-cervical-spine-fracture-detection')\ntrain_folder = root/'train_images'\ntrain_df = pd.read_csv(\"../input/rsna-2022-cervical-spine-fracture-detection/train.csv\")\nstudy_ids = train_df.StudyInstanceUID.unique()","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:18:54.23218Z","iopub.execute_input":"2022-09-16T21:18:54.233294Z","iopub.status.idle":"2022-09-16T21:18:54.26344Z","shell.execute_reply.started":"2022-09-16T21:18:54.233252Z","shell.execute_reply":"2022-09-16T21:18:54.262242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_fn(study_id, slice):\n    return train_folder/study_id/f'{slice}.dcm'","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:18:54.265819Z","iopub.execute_input":"2022-09-16T21:18:54.266324Z","iopub.status.idle":"2022-09-16T21:18:54.271789Z","shell.execute_reply.started":"2022-09-16T21:18:54.266286Z","shell.execute_reply":"2022-09-16T21:18:54.270578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fn = get_fn('1.2.826.0.1.3680043.10001',100)\ndcm = pydicom.dcmread(fn)\narr = dcm.pixel_array\nplt.hist(arr)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:19:41.672931Z","iopub.execute_input":"2022-09-16T21:19:41.673495Z","iopub.status.idle":"2022-09-16T21:19:51.94068Z","shell.execute_reply.started":"2022-09-16T21:19:41.673449Z","shell.execute_reply":"2022-09-16T21:19:51.939566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conversion Functions","metadata":{}},{"cell_type":"code","source":"#my function\n\ndef dcm_apply_windows(dcm, width = None, level = None):\n    arr = dcm.pixel_array\n    #slope, intercept\n    slope = 1\n    intercept = 0\n    if \"RescaleIntercept\" in dcm and \"RescaleSlope\" in dcm:\n        intercept = int(dcm.RescaleIntercept)\n        slope = int(dcm.RescaleSlope)\n\n    arr = arr * slope + intercept\n\n    #window\n    if width is None or level is None:\n        width,level = get_window_from_dicom(dcm)\n    arr = np.clip(arr, level - width // 2, level + width // 2)\n\n    #scale\n    #this was fixed\n    arr = (((arr - np.min(arr)) / np.max(arr)) * 256).astype(np.uint8) \n    return arr\n\n\n# def dcm_apply_windows_old(dcm, scale = False):\n#     \"\"\"\n#     Applies Intercept/Slope and Window Center/Width\n#     \"\"\"\n#     arr = dcm.pixel_array\n#     #slope, intercept\n#     slope = 1\n#     intercept = 0\n#     if \"RescaleIntercept\" in dcm and \"RescaleSlope\" in dcm:\n#         intercept = int(dcm.RescaleIntercept)\n#         slope = int(dcm.RescaleSlope)\n        \n#     arr = arr * slope + intercept\n    \n#     #window\n#     width,level = get_window_from_dicom(dcm)\n#     if width is not None and level is not None:\n#         arr = np.clip(arr, level - width // 2, level + width // 2)\n        \n#      #scale\n#     if scale:\n#         arr = (((arr - np.min(arr)) / np.max(arr)) * 256).astype(np.uint8)\n        \n#     return arr\n\ndef get_window_from_dicom(dcm, default = (2000, 500)):\n    \"\"\"\n    Returns window width and window center values or first example if MultiValue\n    Strips comma from value if present (seen in a different dataset)\n    If no window width/level is provided or available, returns default.\n    \"\"\"\n    width, level = default\n\n    if \"WindowWidth\" in dcm:\n        width = dcm.WindowWidth\n        if isinstance(width, pydicom.multival.MultiValue):\n            width = float(width[0])\n        else:\n            width = float(str(width).replace(',', ''))\n\n    if \"WindowCenter\" in dcm:\n        level = dcm.WindowCenter\n        if isinstance(level, pydicom.multival.MultiValue):\n            level = float(level[0])\n        else:\n            level = float(str(level).replace(',', ''))\n            \n    return width, level","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:31:31.143441Z","iopub.execute_input":"2022-09-16T21:31:31.143861Z","iopub.status.idle":"2022-09-16T21:31:31.16049Z","shell.execute_reply.started":"2022-09-16T21:31:31.143828Z","shell.execute_reply":"2022-09-16T21:31:31.159237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_dicom(path, voi_lut = True, fix_monochrome = True):\n    '''ref: https://www.kaggle.com/code/raddar/convert-dicom-to-np-array-the-correct-way\n    '''\n    dicom = pydicom.read_file(path)\n    \n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom.pixel_array, dicom)\n    else:\n        data = dicom.pixel_array\n               \n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n        \n    data = data - np.min(data)\n    data = data / np.max(data)\n    data = (data * 255).astype(np.uint8)\n    return data, dicom","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:28:55.500979Z","iopub.execute_input":"2022-09-16T21:28:55.502069Z","iopub.status.idle":"2022-09-16T21:28:55.510006Z","shell.execute_reply.started":"2022-09-16T21:28:55.502025Z","shell.execute_reply":"2022-09-16T21:28:55.508758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Code from:\n#https://www.kaggle.com/code/thedevastator/tf-rsna-efficient-net-baseline\n#https://www.kaggle.com/code/realneuralnetwork/rsna-efficientnet-infer\ndef load_dicom_from_path(path, size = 512):\n    try:\n        dcm = pydicom.dcmread(path)\n        dcm.PhotometricInterpretation = 'YBR_FULL'\n        data=dcm.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 cv2.cvtColor(data.reshape(size, size), cv2.COLOR_GRAY2RGB)\n    except:     \n        print('Error')\n        return np.zeros((size, size, 3))","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:28:56.034351Z","iopub.execute_input":"2022-09-16T21:28:56.035723Z","iopub.status.idle":"2022-09-16T21:28:56.046146Z","shell.execute_reply.started":"2022-09-16T21:28:56.035667Z","shell.execute_reply":"2022-09-16T21:28:56.044539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def show_conversions(path):\n    plt.figure(figsize=(20, 20))\n\n    r = 5\n    c = 5\n\n    ax = plt.subplot(r,c, 1)\n    ax.set_title('read_dicom, no lut')\n    arr, _ = read_dicom(path, voi_lut = False)\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 2)\n    ax.set_title('read_dicom, + lut')\n    arr, _ = read_dicom(path, voi_lut = True)\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 3)\n    ax.set_title('load_dicom_from_path')\n    im = load_dicom_from_path(path)\n    plt.axis('off')   \n    plt.imshow(im[:,:,0], cmap='bone')\n\n    ax = plt.subplot(r,c, 4)\n    ax.set_title('dcm_apply_windows')\n    dcm = pydicom.dcmread(path)\n    arr = dcm_apply_windows(dcm)\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 5)\n    ax.set_title('dcm_apply_windows + 2000/600')\n    dcm = pydicom.dcmread(path)\n    arr = dcm_apply_windows(dcm,width=2000,level=600)\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 6)\n    ax.set_title('raw')\n    arr = dcm.pixel_array\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 7)\n    ax.set_title('slope & intercept only')\n    arr = dcm.pixel_array\n    slope = 1\n    intercept = 0\n    if \"RescaleIntercept\" in dcm and \"RescaleSlope\" in dcm:\n        intercept = int(dcm.RescaleIntercept)\n        slope = int(dcm.RescaleSlope)\n    arr = arr * slope + intercept\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')\n\n    ax = plt.subplot(r,c, 8)\n    ax.set_title('WW & WC only')\n    arr = dcm.pixel_array\n    width,level = get_window_from_dicom(dcm)\n    arr = np.clip(arr, level - width // 2, level + width // 2)\n    im = Image.fromarray(arr)\n    plt.axis('off')   \n    plt.imshow(im,cmap='bone')","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:50:33.060988Z","iopub.execute_input":"2022-09-16T21:50:33.061411Z","iopub.status.idle":"2022-09-16T21:50:33.077861Z","shell.execute_reply.started":"2022-09-16T21:50:33.061376Z","shell.execute_reply":"2022-09-16T21:50:33.076654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Comparison","metadata":{}},{"cell_type":"markdown","source":"# Image without GDCM requirement","metadata":{}},{"cell_type":"code","source":"filename_1 = get_fn('1.2.826.0.1.3680043.10001',100)\nshow_conversions(filename_1)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:50:34.701526Z","iopub.execute_input":"2022-09-16T21:50:34.701948Z","iopub.status.idle":"2022-09-16T21:50:35.614831Z","shell.execute_reply.started":"2022-09-16T21:50:34.70191Z","shell.execute_reply":"2022-09-16T21:50:35.613626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image with GDCM requirement - JPEG Lossless","metadata":{}},{"cell_type":"code","source":"#requires GDCM\nfilename_2 = get_fn(study_ids[100],100)\nshow_conversions(filename_2)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:50:35.956756Z","iopub.execute_input":"2022-09-16T21:50:35.957916Z","iopub.status.idle":"2022-09-16T21:50:36.912026Z","shell.execute_reply.started":"2022-09-16T21:50:35.957869Z","shell.execute_reply":"2022-09-16T21:50:36.911079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Different JPEG Lossless","metadata":{}},{"cell_type":"code","source":"study_id = '1.2.826.0.1.3680043.23904'\n#requires GDCM\nfilename_3 = get_fn(study_id,180)\nshow_conversions(filename_3)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T21:50:37.492054Z","iopub.execute_input":"2022-09-16T21:50:37.493488Z","iopub.status.idle":"2022-09-16T21:50:38.409469Z","shell.execute_reply.started":"2022-09-16T21:50:37.493395Z","shell.execute_reply":"2022-09-16T21:50:38.408171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DICOM values for the sample images","metadata":{}},{"cell_type":"code","source":"pydicom.dcmread(filename_1, stop_before_pixels=True)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T18:12:11.901364Z","iopub.execute_input":"2022-09-16T18:12:11.902451Z","iopub.status.idle":"2022-09-16T18:12:11.9135Z","shell.execute_reply.started":"2022-09-16T18:12:11.902414Z","shell.execute_reply":"2022-09-16T18:12:11.912025Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pydicom.dcmread(filename_2, stop_before_pixels=True)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T18:12:11.91539Z","iopub.execute_input":"2022-09-16T18:12:11.916196Z","iopub.status.idle":"2022-09-16T18:12:11.928563Z","shell.execute_reply.started":"2022-09-16T18:12:11.916157Z","shell.execute_reply":"2022-09-16T18:12:11.927201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pydicom.dcmread(filename_3, stop_before_pixels=True)","metadata":{"execution":{"iopub.status.busy":"2022-09-16T18:12:11.932064Z","iopub.execute_input":"2022-09-16T18:12:11.932447Z","iopub.status.idle":"2022-09-16T18:12:11.944958Z","shell.execute_reply.started":"2022-09-16T18:12:11.932405Z","shell.execute_reply":"2022-09-16T18:12:11.94338Z"},"trusted":true},"execution_count":null,"outputs":[]}]}