{"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":"I've noticed that some of the min-max normalized images are too dark or too bright with low-contrast.  \nIn this notebook, I tried to find better method to normalize pixel_arrays.","metadata":{}},{"cell_type":"code","source":"!conda install gdcm -c conda-forge -y","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-05-25T06:38:57.889625Z","iopub.execute_input":"2021-05-25T06:38:57.890019Z","iopub.status.idle":"2021-05-25T06:40:03.010093Z","shell.execute_reply.started":"2021-05-25T06:38:57.889939Z","shell.execute_reply":"2021-05-25T06:40:03.008778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pathlib import Path\nfrom tqdm import tqdm\nimport pydicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport cv2\nimport numpy as np\nimport matplotlib.pyplot as plt\n\ntrain_dir = Path('/kaggle/input/siim-covid19-detection/train')\nfilenames = sorted(train_dir.glob('**/*.dcm'))","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:43:48.179933Z","iopub.execute_input":"2021-05-25T06:43:48.180626Z","iopub.status.idle":"2021-05-25T06:43:57.288834Z","shell.execute_reply.started":"2021-05-25T06:43:48.180569Z","shell.execute_reply":"2021-05-25T06:43:57.287761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Is lookup table or windowing available at all?","metadata":{}},{"cell_type":"code","source":"def is_lut_available(ds):\n    '''\n    Original from pydicom/util.py (apply_voi_lut)\n    '''\n    valid_voi = False\n    if 'VOILUTSequence' in ds:\n        ds.VOILUTSequence = cast(List[\"Dataset\"], ds.VOILUTSequence)\n        valid_voi = None not in [\n            ds.VOILUTSequence[0].get('LUTDescriptor', None),\n            ds.VOILUTSequence[0].get('LUTData', None)\n        ]\n    valid_windowing = None not in [\n        ds.get('WindowCenter', None),\n        ds.get('WindowWidth', None)\n    ]\n    return valid_voi or valid_windowing\n\nany_lut_available = False\nfor filename in tqdm(filenames):\n    dcm = pydicom.dcmread(filename, stop_before_pixels=True)\n    any_lut_available |= is_lut_available(dcm)\n\nprint('At least one lut is available.' if any_lut_available else 'There is no lut available.')","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:43:57.290268Z","iopub.execute_input":"2021-05-25T06:43:57.290547Z","iopub.status.idle":"2021-05-25T06:44:27.984856Z","shell.execute_reply.started":"2021-05-25T06:43:57.290521Z","shell.execute_reply":"2021-05-25T06:44:27.983751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So I can't rely on lookup tables or pre-defined windows.\n\n## Min-max normalization","metadata":{"execution":{"iopub.status.busy":"2021-05-24T12:01:34.734547Z","iopub.execute_input":"2021-05-24T12:01:34.73496Z","iopub.status.idle":"2021-05-24T12:01:34.74023Z","shell.execute_reply.started":"2021-05-24T12:01:34.734924Z","shell.execute_reply":"2021-05-24T12:01:34.738716Z"}}},{"cell_type":"code","source":"N = 10\nIMAGE_SIZE = 128\n\ndef resize(img):\n    return cv2.resize(img, (IMAGE_SIZE, IMAGE_SIZE))\n\ndef read_xray(path, voi_lut = True, fix_monochrome = True):\n    # Original from: https://www.kaggle.com/raddar/convert-dicom-to-np-array-the-correct-way\n    dicom = pydicom.read_file(path)\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        \n    return data","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:44:27.987285Z","iopub.execute_input":"2021-05-25T06:44:27.987752Z","iopub.status.idle":"2021-05-25T06:44:27.996525Z","shell.execute_reply.started":"2021-05-25T06:44:27.987707Z","shell.execute_reply":"2021-05-25T06:44:27.995289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imgs = np.stack([resize(read_xray(filename)) for filename in tqdm(filenames[:N**2])])\ntile = np.concatenate(np.concatenate(imgs.reshape((N, N, IMAGE_SIZE, IMAGE_SIZE)), axis=1), axis=1)\nplt.figure(figsize=(20, 20))\nplt.imshow(tile, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:44:27.998209Z","iopub.execute_input":"2021-05-25T06:44:27.998547Z","iopub.status.idle":"2021-05-25T06:44:55.217064Z","shell.execute_reply.started":"2021-05-25T06:44:27.998504Z","shell.execute_reply":"2021-05-25T06:44:55.215393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are quite a few too dark/bright images.\n\n## Too dark/bright images\nI picked up some samples of too dark/bright images for the inspection","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\ndark_bright_filenames = ['1e96d5eb4c91/703a9f4c4ffb/7f1924880cf8.dcm', '1e96d5eb4c91/3928035a1e6d/3d12cb6aad8b.dcm', '1a53a4506f10/db8d01dc5e4f/95bf5f4f6153.dcm', '047b450939fd/76eb6433f31d/cc54237ef4db.dcm']\n\nfor filename in dark_bright_filenames:\n    dcm = pydicom.dcmread(train_dir / filename)\n    pixel_array = dcm.pixel_array\n    print(filename)\n    print('min, max =',pixel_array.min(), pixel_array.max())\n    plt.figure(figsize=(10, 5))\n    plt.subplot(1,2,1)\n    plt.imshow(read_xray(train_dir / filename), cmap='gray')\n    plt.subplot(1,2,2)\n    sns.histplot(pixel_array.ravel())\n    plt.show()\n","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:44:55.218494Z","iopub.execute_input":"2021-05-25T06:44:55.218806Z","iopub.status.idle":"2021-05-25T06:45:32.606551Z","shell.execute_reply.started":"2021-05-25T06:44:55.218774Z","shell.execute_reply":"2021-05-25T06:45:32.605428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Apparently, these images contain some outliers in pixel values.\n\n## Robust pixel_array scaling\nThe basic idea here is to exclude outlier pixels in the scaling process.","metadata":{}},{"cell_type":"code","source":"def robust_read_xray(path, voi_lut = True, fix_monochrome = True):\n    # Original from: https://www.kaggle.com/raddar/convert-dicom-to-np-array-the-correct-way\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    q1,q2,q3 = np.quantile(data, [.25,.5,.75])\n    iqr = q3 - q1\n    multiplier = 2\n    mask = ((q2 - multiplier * iqr) < data) & (data < (q2 + multiplier * iqr))\n    p = .001\n    data = data.astype(np.float32) - np.quantile(data[mask], p)\n    data = data / np.quantile(data[mask], 1-p)\n    data = np.clip(data, 0, 1)\n    return data","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:45:32.608308Z","iopub.execute_input":"2021-05-25T06:45:32.608725Z","iopub.status.idle":"2021-05-25T06:45:32.617516Z","shell.execute_reply.started":"2021-05-25T06:45:32.608683Z","shell.execute_reply":"2021-05-25T06:45:32.616450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"imgs = np.stack([resize(robust_read_xray(filename)) for filename in tqdm(filenames[:N**2])])\ntile = np.concatenate(np.concatenate(imgs.reshape((N,N,IMAGE_SIZE, IMAGE_SIZE)), axis=1), axis=1)\nplt.figure(figsize=(20,20))\nplt.imshow(tile, cmap='gray')","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:45:32.619200Z","iopub.execute_input":"2021-05-25T06:45:32.619886Z","iopub.status.idle":"2021-05-25T06:45:59.369256Z","shell.execute_reply.started":"2021-05-25T06:45:32.619841Z","shell.execute_reply":"2021-05-25T06:45:59.368198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Side by side comparison","metadata":{}},{"cell_type":"code","source":"for filename in dark_bright_filenames:\n    print(filename)\n    plt.subplot(1,2,1)\n    plt.imshow(read_xray(train_dir / filename), cmap='gray')\n    plt.title('min-max')\n    plt.subplot(1,2,2)\n    plt.imshow(robust_read_xray(train_dir / filename), cmap='gray')\n    plt.title('robust')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2021-05-25T06:45:59.370559Z","iopub.execute_input":"2021-05-25T06:45:59.370856Z","iopub.status.idle":"2021-05-25T06:46:04.570120Z","shell.execute_reply.started":"2021-05-25T06:45:59.370825Z","shell.execute_reply":"2021-05-25T06:46:04.569103Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There maybe better method or parameters (multiplier and p in read_xray_robust) for the pixel_array normalization.","metadata":{}}]}