{"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":"# Why this notebook?\n\nAfter reading this post about [Pydicom vs DicomSDL (Fast dicom export and processing (1.6-2x faster) 💪)](https://www.kaggle.com/competitions/rsna-breast-cancer-detection/discussion/371033), I checked whether DicomSDL had the option to apply a VOI lookup table or windowing operation like Pydicom.<br>\nIt turns out, they don't. They only have this function in their source code 😅\n\n    def apply_window_center():\n        print(\"Hello\")\n        \n\n## Apply_voi_lut\n\nNext, I dug into the source code of Pydicom to see what they do exactly inside the apply_voi_lut method.<br>\nThe code checks several things, of which many don't apply for the images/datasets in the dicom images in this competition.<br>\nFor example, **none** of the datasets in this competition have a valid **VOILUTSequence**.<br>\nThis means the **apply_voi_lut function** will always redirect to the **apply_windowing**.<br>\nTherefore, the next step was to check this function inside Pydicom's code.\n\n## Apply_windowing\n\nHere, the most important thing the code checks is which VOILUTFunction is used in the dataset.<br>\nIn this competition, 22% of the datasets (12139 of 54706) have the Sigmoid function, the rest uses the LINEAR function.<br>\nAs the code is very different for each of these functions, I recreated/copied the necessary code for both functions.<br>\nPlease note that I removed some code for speed, as many checks in the pydicom code is not not applicable in this competition.\nFor example:\n- Pixelrepresentation and ModalityLUTSequence is always NaN\n- Window width/center is always set correctly\n- RescaleSlope & RescaleIntercept is always the same\n\n\n## Results/conclusion\n\nAs you can see in the results below, using the VOILUTFunction is quite important.<br>\nAround 22% of the images uses SIGMOID instead of LINEAR, and the images/colors for Sigmoid are very different, a lot lighter/more yellow.<br>\nBut once you apply the windowing function, they seem completely similar to images with the Linear VOILUTFunction.<br>\n**Therefore**, if you use DicomSDL instead of Pydicom, you should definitely create your own apply_windowing function like in Pydicom!\n","metadata":{}},{"cell_type":"markdown","source":"### ⬇ Libraries","metadata":{}},{"cell_type":"code","source":"# Installs\n!pip install -qU \"python-gdcm\" pydicom pylibjpeg \"opencv-python-headless\"\n!pip install /kaggle/input/dicomsdl/dicomsdl-0.109.1-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl\n\n# Imports\nimport os\nimport time\nimport pydicom\nimport dicomsdl as pysdl\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom tqdm import tqdm\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-01-15T16:11:21.464361Z","iopub.execute_input":"2023-01-15T16:11:21.465134Z","iopub.status.idle":"2023-01-15T16:11:51.377542Z","shell.execute_reply.started":"2023-01-15T16:11:21.465038Z","shell.execute_reply":"2023-01-15T16:11:51.376245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Loading images with Pydicom & DicomSDL","metadata":{}},{"cell_type":"code","source":"def load_image_pydicom(img_path, voi_lut=False):\n    dataset = pydicom.dcmread(img_path)\n    img = dataset.pixel_array\n    if voi_lut:\n        img = apply_voi_lut(img, dataset)\n    if dataset.PhotometricInterpretation == \"MONOCHROME1\":\n        img = np.amax(img) - img\n    return img\n\n\ndef load_image_dicomsdl(img_path, voi_lut=False):\n    dataset = pysdl.open(img_path)\n    img = dataset.pixelData()\n    \n    if voi_lut:\n        # Load only the variables we need\n        center = dataset[\"WindowCenter\"]\n        width = dataset[\"WindowWidth\"]\n        bits_stored = dataset[\"BitsStored\"]\n        voi_lut_function = dataset[\"VOILUTFunction\"]\n\n        # For sigmoid it's a list, otherwise a single value\n        if isinstance(center, list):\n            center = center[0]\n        if isinstance(width, list):\n            width = width[0]\n\n        # Set y_min, max & range\n        y_min = 0\n        y_max = float(2**bits_stored - 1)\n        y_range = y_max\n\n        # Function with default LINEAR (so for Nan, it will use linear)\n        if voi_lut_function == \"SIGMOID\":\n            img = y_range / (1 + np.exp(-4 * (img - center) / width)) + y_min\n        else:\n            # Checks width for < 1 (in our case not necessary, always >= 750)\n            center -= 0.5\n            width -= 1\n\n            below = img <= (center - width / 2)\n            above = img > (center + width / 2)\n            between = np.logical_and(~below, ~above)\n\n            img[below] = y_min\n            img[above] = y_max\n            if between.any():\n                img[between] = (\n                    ((img[between] - center) / width + 0.5) * y_range + y_min\n                )\n    \n    if dataset[\"PhotometricInterpretation\"] == \"MONOCHROME1\":\n        img = np.amax(img) - img\n\n    return img","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-15T16:11:51.380457Z","iopub.execute_input":"2023-01-15T16:11:51.380886Z","iopub.status.idle":"2023-01-15T16:11:51.393108Z","shell.execute_reply.started":"2023-01-15T16:11:51.380849Z","shell.execute_reply":"2023-01-15T16:11:51.392105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Resulting images with VOILUTFunction -> Linear\n\nAs you can see below, images with VOILUTFunction Linear doesn't change at all.","metadata":{}},{"cell_type":"code","source":"img_path_linear = \"../input/rsna-breast-cancer-detection/train_images/10102/1181635673.dcm\"\nfig, axs = plt.subplots(nrows=1, ncols=4,figsize=(30, 16))\n\nimg = load_image_pydicom(img_path_linear)\nnumpydata = np.asarray(img)\nprint(numpydata.shape)\nplt.subplot(141)\nplt.title(\"Pydicom - Linear without VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_dicomsdl(img_path_linear)\nnumpydata = np.asarray(img)\nprint(numpydata.shape)\nplt.subplot(142)\nplt.title(\"DicomSDL - Linear without VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_pydicom(img_path_linear, True)\nnumpydata = np.asarray(img)\nprint(numpydata.shape)\nplt.subplot(143)\nplt.title(\"Pydicom - Linear with VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_dicomsdl(img_path_linear, True)\nnumpydata = np.asarray(img)\nvalue=numpydata.shape\nprint(numpydata.shape)\nplt.subplot(144)\nplt.title(\"DicomSDL - Linear with VOI-LUT\")\nplt.imshow(img)\n\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-15T16:57:58.765567Z","iopub.execute_input":"2023-01-15T16:57:58.766012Z","iopub.status.idle":"2023-01-15T16:58:06.418091Z","shell.execute_reply.started":"2023-01-15T16:57:58.765979Z","shell.execute_reply":"2023-01-15T16:58:06.417299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"import pywt\nimport pywt.data","metadata":{}},{"cell_type":"code","source":"import pywt\nimport pywt.data","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:58:06.419692Z","iopub.execute_input":"2023-01-15T16:58:06.420527Z","iopub.status.idle":"2023-01-15T16:58:06.424819Z","shell.execute_reply.started":"2023-01-15T16:58:06.420493Z","shell.execute_reply":"2023-01-15T16:58:06.423796Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Wavelet transform of image, and plot approximation and details\n# dicom = pydicom.dcmread(train_images[15])\n# img = dicom.pixel_array\n# print('max and min for image before normalisation',img.max(),img.min())\n# img = (img - img.min()) / (img.max() - img.min())\n# print('max and min for image',img.max(),img.min())\n# if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n#     img = 1 - img\ntitles = ['Approximation', ' Horizontal detail',\n          'Vertical detail', 'Diagonal detail']\n# coeffs2 = pywt.dwt2(img, 'bior1.3')\nimg_array=img\nimg_array /= 255;\ncoeffs2 = pywt.dwt2(img_array, 'haar')\nLL, (LH, HL, HH) = coeffs2\nHH = pywt.threshold(HH, value=-0.4, mode='soft')\nHL = pywt.threshold(HL, value=-1.5, mode='hard')\nLH = pywt.threshold(LH, value=-1.5, mode='hard')\n\nfig = plt.figure(figsize=(15, 5))\nfor i, a in enumerate([LL, LH, HL, HH]):\n    ax = fig.add_subplot(1, 4, i + 1)\n    ax.imshow(a, interpolation=\"nearest\", cmap=plt.cm.magma)\n    ax.set_title(titles[i], fontsize=10)\n    ax.set_xticks([])\n    ax.set_yticks([])\n\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:58:06.426174Z","iopub.execute_input":"2023-01-15T16:58:06.427265Z","iopub.status.idle":"2023-01-15T16:58:07.219015Z","shell.execute_reply.started":"2023-01-15T16:58:06.427232Z","shell.execute_reply":"2023-01-15T16:58:07.217946Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def describe_details(arr):# measures of dispersion\n    min = np.amin(arr)\n    max = np.amax(arr)\n    range = np.ptp(arr)\n    variance = np.var(arr)\n    sd = np.std(arr)\n    mean = np.mean(arr)\n    median = np.median(arr)\n\n    #print(\"Array =\", arr)\n    print(\"Measures of Dispersion\")\n    print(\"Minimum =\", min)\n    print(\"Maximum =\", max)\n    print(\"Range =\", range)\n    print(\"Variance =\", variance)\n    print(\"Standard Deviation =\", sd)\n    print(\"Mean =\", mean)\n    print(\"Median =\", median)","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:52:05.026012Z","iopub.execute_input":"2023-01-15T16:52:05.026971Z","iopub.status.idle":"2023-01-15T16:52:05.034571Z","shell.execute_reply.started":"2023-01-15T16:52:05.026927Z","shell.execute_reply":"2023-01-15T16:52:05.033272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"describe_details(HL)","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:52:05.036022Z","iopub.execute_input":"2023-01-15T16:52:05.036334Z","iopub.status.idle":"2023-01-15T16:52:05.097155Z","shell.execute_reply.started":"2023-01-15T16:52:05.036308Z","shell.execute_reply":"2023-01-15T16:52:05.096003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HH","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:52:05.098428Z","iopub.execute_input":"2023-01-15T16:52:05.098776Z","iopub.status.idle":"2023-01-15T16:52:05.107067Z","shell.execute_reply.started":"2023-01-15T16:52:05.098744Z","shell.execute_reply":"2023-01-15T16:52:05.105959Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Resulting images with VOILUTFunction -> Sigmoid\n\nFor VOILUTFunction Sigmoid, however, it does have an effect. You can clearly see that the default images for Sigmoid without VOILUT are a lot more light/yellow, and once VOILUT is applied the colors seem more similar to images with VOILUTFunction Linear.","metadata":{}},{"cell_type":"code","source":"img_path_sigmoid = \"../input/rsna-breast-cancer-detection/train_images/10006/1459541791.dcm\"\nfig, axs = plt.subplots(nrows=1, ncols=4,figsize=(30, 16))\n\nimg = load_image_pydicom(img_path_sigmoid)\n\nplt.subplot(141)\nplt.title(\"Pydicom - Sigmoid without VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_dicomsdl(img_path_sigmoid)\nplt.subplot(142)\nplt.title(\"DicomSDL - Sigmoid without VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_pydicom(img_path_sigmoid, True)\nplt.subplot(143)\nplt.title(\"Pydicom - Sigmoid with VOI-LUT\")\nplt.imshow(img)\n\nimg = load_image_dicomsdl(img_path_sigmoid, True)\nplt.subplot(144)\nplt.title(\"DicomSDL - Sigmoid with VOI-LUT\")\nplt.imshow(img)\n\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-01-15T16:52:05.10829Z","iopub.execute_input":"2023-01-15T16:52:05.108604Z","iopub.status.idle":"2023-01-15T16:52:21.988772Z","shell.execute_reply.started":"2023-01-15T16:52:05.108577Z","shell.execute_reply":"2023-01-15T16:52:21.987652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Wavelet transform of image, and plot approximation and details\n# dicom = pydicom.dcmread(train_images[15])\n# img = dicom.pixel_array\n# print('max and min for image before normalisation',img.max(),img.min())\n# img = (img - img.min()) / (img.max() - img.min())\n# print('max and min for image',img.max(),img.min())\n# if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n#     img = 1 - img\ntitles = ['Approximation', ' Horizontal detail',\n          'Vertical detail', 'Diagonal detail']\n# coeffs2 = pywt.dwt2(img, 'bior1.3')\nimg_array=img\nimg_array /= 255;\ncoeffs2 = pywt.dwt2(img_array, 'haar')\nLL, (LH, HL, HH) = coeffs2\nHH = pywt.threshold(HH, value=-0.4, mode='soft')\nHL = pywt.threshold(HL, value=-0.4, mode='hard')\nLH = pywt.threshold(LH, value=-1.5, mode='hard')\n\nfig = plt.figure(figsize=(15, 5))\nfor i, a in enumerate([LL, LH, HL, HH]):\n    ax = fig.add_subplot(1, 4, i + 1)\n    ax.imshow(a, interpolation=\"nearest\", cmap=plt.cm.magma)\n    ax.set_title(titles[i], fontsize=10)\n    ax.set_xticks([])\n    ax.set_yticks([])\n\nfig.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-15T16:53:36.683348Z","iopub.execute_input":"2023-01-15T16:53:36.68377Z","iopub.status.idle":"2023-01-15T16:53:38.23396Z","shell.execute_reply.started":"2023-01-15T16:53:36.683735Z","shell.execute_reply":"2023-01-15T16:53:38.232815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}