{"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":"# Description\nIn new version, all JPEG2000 images will be processed with DALI if GPU is avaliable, because processing is much faster this way.\nIn this notebook will be showed all functionalits that does all the processing: \n* windowing - performs windowing operation so that there is better contrast between soft and hard tissue\n* shuffling - shuffles input data\n* resizing - resize to preffered dimensions\n* saving - save images to output directory\n* converting images to one side unifomrly (left or right).\n\nAlso, this class can be used as a part of the pipeline, to preprocess all input images and loads them to memory.\n\nNote: dataset can be found on link: https://www.kaggle.com/competitions/rsna-breast-cancer-detection/data\n<br>or by running a comamand: kaggle competitions download -c rsna-breast-cancer-detection\n\n","metadata":{}},{"cell_type":"code","source":"# Install required libraries\n!pip install -qU python-gdcm pydicom pylibjpeg\n!pip install -qU dicomsdl\n\n!pip install --extra-index-url https://developer.download.nvidia.com/compute/redist/nightly --upgrade nvidia-dali-nightly-cuda110","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:39:26.922983Z","iopub.execute_input":"2023-04-07T16:39:26.923376Z","iopub.status.idle":"2023-04-07T16:40:32.036768Z","shell.execute_reply.started":"2023-04-07T16:39:26.923298Z","shell.execute_reply":"2023-04-07T16:40:32.035565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Imports","metadata":{}},{"cell_type":"code","source":"import os\nimport time\nimport numpy as np\nimport pandas as pd\n# image manipulation\nimport cv2\nimport PIL\nfrom PIL import Image\n# dicom image \nimport pydicom\nimport dicomsdl\nfrom pydicom.pixel_data_handlers import apply_windowing\n\n# visualisation\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# helpers\nfrom tqdm.notebook import tqdm\nimport time\nimport copy\nimport gc\nfrom enum import Enum\nfrom joblib import Parallel, delayed\nimport glob\nimport shutil\n\n\n# for cnn\nfrom torchvision import transforms\n\nimport torch\nimport nvidia.dali.fn as fn\nimport nvidia.dali.types as types\nfrom nvidia.dali import pipeline_def\nfrom nvidia.dali.types import DALIDataType\nfrom pydicom.filebase import DicomBytesIO\n","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:32.039435Z","iopub.execute_input":"2023-04-07T16:40:32.040523Z","iopub.status.idle":"2023-04-07T16:40:34.586604Z","shell.execute_reply.started":"2023-04-07T16:40:32.040475Z","shell.execute_reply":"2023-04-07T16:40:34.585616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#%%writefile /opt/conda/lib/python3.7/site-packages/nvidia/dali/plugin/pytorch.py\n\n# Copyright (c) 2017-2022, NVIDIA CORPORATION & AFFILIATES. All rights reserved.\n#\n# Licensed under the Apache License, Version 2.0 (the \"License\");\n# you may not use this file except in compliance with the License.\n# You may obtain a copy of the License at\n#\n#     http://www.apache.org/licenses/LICENSE-2.0\n#\n# Unless required by applicable law or agreed to in writing, software\n# distributed under the License is distributed on an \"AS IS\" BASIS,\n# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n# See the License for the specific language governing permissions and\n# limitations under the License.\n\nfrom nvidia.dali.backend import TensorGPU, TensorListGPU\nfrom nvidia.dali.pipeline import Pipeline\nimport nvidia.dali.ops as ops\nfrom nvidia.dali import types\nfrom nvidia.dali.plugin.base_iterator import _DaliBaseIterator\nfrom nvidia.dali.plugin.base_iterator import LastBatchPolicy\nimport torch\nimport torch.utils.dlpack as torch_dlpack\nimport ctypes\nimport numpy as np\n\nto_torch_type = {\n    types.DALIDataType.FLOAT:   torch.float32,\n    types.DALIDataType.FLOAT64: torch.float64,\n    types.DALIDataType.FLOAT16: torch.float16,\n    types.DALIDataType.UINT8:   torch.uint8,\n    types.DALIDataType.INT8:    torch.int8,\n    types.DALIDataType.UINT16:  torch.int16,\n    types.DALIDataType.INT16:   torch.int16,\n    types.DALIDataType.INT32:   torch.int32,\n    types.DALIDataType.INT64:   torch.int64\n}\n\n\ndef feed_ndarray(dali_tensor, arr, cuda_stream=None):\n    \"\"\"\n    Copy contents of DALI tensor to PyTorch's Tensor.\n\n    Parameters\n    ----------\n    `dali_tensor` : nvidia.dali.backend.TensorCPU or nvidia.dali.backend.TensorGPU\n                    Tensor from which to copy\n    `arr` : torch.Tensor\n            Destination of the copy\n    `cuda_stream` : torch.cuda.Stream, cudaStream_t or any value that can be cast to cudaStream_t.\n                    CUDA stream to be used for the copy\n                    (if not provided, an internal user stream will be selected)\n                    In most cases, using pytorch's current stream is expected (for example,\n                    if we are copying to a tensor allocated with torch.zeros(...))\n    \"\"\"\n    dali_type = to_torch_type[dali_tensor.dtype]\n\n    assert dali_type == arr.dtype, (\"The element type of DALI Tensor/TensorList\"\n                                    \" doesn't match the element type of the target PyTorch Tensor: \"\n                                    \"{} vs {}\".format(dali_type, arr.dtype))\n    assert dali_tensor.shape() == list(arr.size()), \\\n        (\"Shapes do not match: DALI tensor has size {0}, but PyTorch Tensor has size {1}\".\n            format(dali_tensor.shape(), list(arr.size())))\n    cuda_stream = types._raw_cuda_stream(cuda_stream)\n\n    # turn raw int to a c void pointer\n    c_type_pointer = ctypes.c_void_p(arr.data_ptr())\n    if isinstance(dali_tensor, (TensorGPU, TensorListGPU)):\n        stream = None if cuda_stream is None else ctypes.c_void_p(cuda_stream)\n        dali_tensor.copy_to_external(c_type_pointer, stream, non_blocking=True)\n    else:\n        dali_tensor.copy_to_external(c_type_pointer)\n    return arr\n\n\nclass DALIGenericIterator(_DaliBaseIterator):\n    \"\"\"\n    General DALI iterator for PyTorch. It can return any number of\n    outputs from the DALI pipeline in the form of PyTorch's Tensors.\n\n    Parameters\n    ----------\n    pipelines : list of nvidia.dali.Pipeline\n                List of pipelines to use\n    output_map : list of str\n                List of strings which maps consecutive outputs\n                of DALI pipelines to user specified name.\n                Outputs will be returned from iterator as dictionary\n                of those names.\n                Each name should be distinct\n    size : int, default = -1\n                Number of samples in the shard for the wrapped pipeline (if there is more than\n                one it is a sum)\n                Providing -1 means that the iterator will work until StopIteration is raised\n                from the inside of iter_setup(). The options `last_batch_policy` and\n                `last_batch_padded` don't work in such case. It works with only one pipeline inside\n                the iterator.\n                Mutually exclusive with `reader_name` argument\n    reader_name : str, default = None\n                Name of the reader which will be queried to the shard size, number of shards and\n                all other properties necessary to count properly the number of relevant and padded\n                samples that iterator needs to deal with. It automatically sets `last_batch_policy`\n                to PARTIAL when the FILL is used, and `last_batch_padded` accordingly to match\n                the reader's configuration\n    auto_reset : string or bool, optional, default = False\n                Whether the iterator resets itself for the next epoch or it requires reset() to be\n                called explicitly.\n\n                It can be one of the following values:\n\n                * ``\"no\"``, ``False`` or ``None`` - at the end of epoch StopIteration is raised\n                  and reset() needs to be called\n                * ``\"yes\"`` or ``\"True\"``- at the end of epoch StopIteration is raised but reset()\n                  is called internally automatically\n\n    dynamic_shape : any, optional,\n                Parameter used only for backward compatibility.\n    fill_last_batch : bool, optional, default = None\n                **Deprecated** Please use ``last_batch_policy`` instead\n\n                Whether to fill the last batch with data up to 'self.batch_size'.\n                The iterator would return the first integer multiple\n                of self._num_gpus * self.batch_size entries which exceeds 'size'.\n                Setting this flag to False will cause the iterator to return\n                exactly 'size' entries.\n    last_batch_policy: optional, default = LastBatchPolicy.FILL\n                What to do with the last batch when there are not enough samples in the epoch\n                to fully fill it. See :meth:`nvidia.dali.plugin.base_iterator.LastBatchPolicy`\n    last_batch_padded : bool, optional, default = False\n                Whether the last batch provided by DALI is padded with the last sample\n                or it just wraps up. In the conjunction with ``last_batch_policy`` it tells\n                if the iterator returning last batch with data only partially filled with\n                data from the current epoch is dropping padding samples or samples from\n                the next epoch. If set to ``False`` next\n                epoch will end sooner as data from it was consumed but dropped. If set to\n                True next epoch would be the same length as the first one. For this to happen,\n                the option `pad_last_batch` in the reader needs to be set to True as well.\n                It is overwritten when `reader_name` argument is provided\n    prepare_first_batch : bool, optional, default = True\n                Whether DALI should buffer the first batch right after the creation of the iterator,\n                so one batch is already prepared when the iterator is prompted for the data\n\n    Example\n    -------\n    With the data set ``[1,2,3,4,5,6,7]`` and the batch size 2:\n\n    last_batch_policy = LastBatchPolicy.PARTIAL, last_batch_padded = True  -> last batch = ``[7]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.PARTIAL, last_batch_padded = False -> last batch = ``[7]``,\n    next iteration will return ``[2, 3]``\n\n    last_batch_policy = LastBatchPolicy.FILL, last_batch_padded = True   -> last batch = ``[7, 7]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.FILL, last_batch_padded = False  -> last batch = ``[7, 1]``,\n    next iteration will return ``[2, 3]``\n\n    last_batch_policy = LastBatchPolicy.DROP, last_batch_padded = True   -> last batch = ``[5, 6]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.DROP, last_batch_padded = False  -> last batch = ``[5, 6]``,\n    next iteration will return ``[2, 3]``\n    \"\"\"\n\n    def __init__(self,\n                 pipelines,\n                 output_map,\n                 size=-1,\n                 reader_name=None,\n                 auto_reset=False,\n                 fill_last_batch=None,\n                 dynamic_shape=False,\n                 last_batch_padded=False,\n                 last_batch_policy=LastBatchPolicy.FILL,\n                 prepare_first_batch=True):\n\n        # check the assert first as _DaliBaseIterator would run the prefetch\n        assert len(set(output_map)) == len(output_map), \"output_map names should be distinct\"\n        self._output_categories = set(output_map)\n        self.output_map = output_map\n\n        _DaliBaseIterator.__init__(self,\n                                   pipelines,\n                                   size,\n                                   reader_name,\n                                   auto_reset,\n                                   fill_last_batch,\n                                   last_batch_padded,\n                                   last_batch_policy,\n                                   prepare_first_batch=prepare_first_batch)\n\n        self._first_batch = None\n        if self._prepare_first_batch:\n            try:\n                self._first_batch = DALIGenericIterator.__next__(self)\n                # call to `next` sets _ever_consumed to True but if we are just calling it from\n                # here we should set if to False again\n                self._ever_consumed = False\n            except StopIteration:\n                assert False, \"It seems that there is no data in the pipeline. This may happen \" \\\n                       \"if `last_batch_policy` is set to PARTIAL and the requested batch size is \" \\\n                       \"greater than the shard size.\"\n\n    def __next__(self):\n        self._ever_consumed = True\n        if self._first_batch is not None:\n            batch = self._first_batch\n            self._first_batch = None\n            return batch\n\n        # Gather outputs\n        outputs = self._get_outputs()\n\n        data_batches = [None for i in range(self._num_gpus)]\n        for i in range(self._num_gpus):\n            dev_id = self._pipes[i].device_id\n            # initialize dict for all output categories\n            category_outputs = dict()\n            # segregate outputs into categories\n            for j, out in enumerate(outputs[i]):\n                category_outputs[self.output_map[j]] = out\n\n            # Change DALI TensorLists into Tensors\n            category_tensors = dict()\n            category_shapes = dict()\n            for category, out in category_outputs.items():\n                category_tensors[category] = out.as_tensor()\n                category_shapes[category] = category_tensors[category].shape()\n\n            category_torch_type = dict()\n            category_device = dict()\n            torch_gpu_device = None\n            torch_cpu_device = torch.device('cpu')\n            # check category and device\n            for category in self._output_categories:\n                category_torch_type[category] = to_torch_type[category_tensors[category].dtype]\n                if type(category_tensors[category]) is TensorGPU:\n                    if not torch_gpu_device:\n                        torch_gpu_device = torch.device('cuda', dev_id)\n                    category_device[category] = torch_gpu_device\n                else:\n                    category_device[category] = torch_cpu_device\n\n            pyt_tensors = dict()\n            for category in self._output_categories:\n                pyt_tensors[category] = torch.empty(category_shapes[category],\n                                                    dtype=category_torch_type[category],\n                                                    device=category_device[category])\n\n            data_batches[i] = pyt_tensors\n\n            # Copy data from DALI Tensors to torch tensors\n            for category, tensor in category_tensors.items():\n                if isinstance(tensor, (TensorGPU, TensorListGPU)):\n                    # Using same cuda_stream used by torch.zeros to set the memory\n                    stream = torch.cuda.current_stream(device=pyt_tensors[category].device)\n                    feed_ndarray(tensor, pyt_tensors[category], cuda_stream=stream)\n                else:\n                    feed_ndarray(tensor, pyt_tensors[category])\n\n        self._schedule_runs()\n\n        self._advance_and_check_drop_last()\n\n        if self._reader_name:\n            if_drop, left = self._remove_padded()\n            if np.any(if_drop):\n                output = []\n                for batch, to_copy in zip(data_batches, left):\n                    batch = batch.copy()\n                    for category in self._output_categories:\n                        batch[category] = batch[category][0:to_copy]\n                    output.append(batch)\n                return output\n\n        else:\n            if self._last_batch_policy == LastBatchPolicy.PARTIAL and (\n                                          self._counter > self._size) and self._size > 0:\n                # First calculate how much data is required to return exactly self._size entries.\n                diff = self._num_gpus * self.batch_size - (self._counter - self._size)\n                # Figure out how many GPUs to grab from.\n                numGPUs_tograb = int(np.ceil(diff / self.batch_size))\n                # Figure out how many results to grab from the last GPU\n                # (as a fractional GPU batch may be required to bring us\n                # right up to self._size).\n                mod_diff = diff % self.batch_size\n                data_fromlastGPU = mod_diff if mod_diff else self.batch_size\n\n                # Grab the relevant data.\n                # 1) Grab everything from the relevant GPUs.\n                # 2) Grab the right data from the last GPU.\n                # 3) Append data together correctly and return.\n                output = data_batches[0:numGPUs_tograb]\n                output[-1] = output[-1].copy()\n                for category in self._output_categories:\n                    output[-1][category] = output[-1][category][0:data_fromlastGPU]\n                return output\n\n        return data_batches\n\n\nclass DALIClassificationIterator(DALIGenericIterator):\n    \"\"\"\n    DALI iterator for classification tasks for PyTorch. It returns 2 outputs\n    (data and label) in the form of PyTorch's Tensor.\n\n    Calling\n\n    .. code-block:: python\n\n       DALIClassificationIterator(pipelines, reader_name)\n\n    is equivalent to calling\n\n    .. code-block:: python\n\n       DALIGenericIterator(pipelines, [\"data\", \"label\"], reader_name)\n\n    Parameters\n    ----------\n    pipelines : list of nvidia.dali.Pipeline\n                List of pipelines to use\n    size : int, default = -1\n                Number of samples in the shard for the wrapped pipeline (if there is more than\n                one it is a sum)\n                Providing -1 means that the iterator will work until StopIteration is raised\n                from the inside of iter_setup(). The options `last_batch_policy` and\n                `last_batch_padded` don't work in such case. It works with only one pipeline inside\n                the iterator.\n                Mutually exclusive with `reader_name` argument\n    reader_name : str, default = None\n                Name of the reader which will be queried to the shard size, number of shards and\n                all other properties necessary to count properly the number of relevant and padded\n                samples that iterator needs to deal with. It automatically sets `last_batch_policy`\n                to PARTIAL when the FILL is used, and `last_batch_padded` accordingly to match\n                the reader's configuration\n    auto_reset : string or bool, optional, default = False\n                Whether the iterator resets itself for the next epoch or it requires reset() to be\n                called explicitly.\n\n                It can be one of the following values:\n\n                * ``\"no\"``, ``False`` or ``None`` - at the end of epoch StopIteration is raised\n                  and reset() needs to be called\n                * ``\"yes\"`` or ``\"True\"``- at the end of epoch StopIteration is raised but reset()\n                  is called internally automatically\n\n    dynamic_shape : any, optional,\n                Parameter used only for backward compatibility.\n    fill_last_batch : bool, optional, default = None\n                **Deprecated** Please use ``last_batch_policy`` instead\n\n                Whether to fill the last batch with data up to 'self.batch_size'.\n                The iterator would return the first integer multiple\n                of self._num_gpus * self.batch_size entries which exceeds 'size'.\n                Setting this flag to False will cause the iterator to return\n                exactly 'size' entries.\n    last_batch_policy: optional, default = LastBatchPolicy.FILL\n                What to do with the last batch when there are not enough samples in the epoch\n                to fully fill it. See :meth:`nvidia.dali.plugin.base_iterator.LastBatchPolicy`\n    last_batch_padded : bool, optional, default = False\n                Whether the last batch provided by DALI is padded with the last sample\n                or it just wraps up. In the conjunction with ``last_batch_policy`` it tells\n                if the iterator returning last batch with data only partially filled with\n                data from the current epoch is dropping padding samples or samples from\n                the next epoch. If set to ``False`` next\n                epoch will end sooner as data from it was consumed but dropped. If set to\n                True next epoch would be the same length as the first one. For this to happen,\n                the option `pad_last_batch` in the reader needs to be set to True as well.\n                It is overwritten when `reader_name` argument is provided\n    prepare_first_batch : bool, optional, default = True\n                Whether DALI should buffer the first batch right after the creation of the iterator,\n                so one batch is already prepared when the iterator is prompted for the data\n\n    Example\n    -------\n    With the data set ``[1,2,3,4,5,6,7]`` and the batch size 2:\n\n    last_batch_policy = LastBatchPolicy.PARTIAL, last_batch_padded = True  -> last batch = ``[7]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.PARTIAL, last_batch_padded = False -> last batch = ``[7]``,\n    next iteration will return ``[2, 3]``\n\n    last_batch_policy = LastBatchPolicy.FILL, last_batch_padded = True   -> last batch = ``[7, 7]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.FILL, last_batch_padded = False  -> last batch = ``[7, 1]``,\n    next iteration will return ``[2, 3]``\n\n    last_batch_policy = LastBatchPolicy.DROP, last_batch_padded = True   -> last batch = ``[5, 6]``,\n    next iteration will return ``[1, 2]``\n\n    last_batch_policy = LastBatchPolicy.DROP, last_batch_padded = False  -> last batch = ``[5, 6]``,\n    next iteration will return ``[2, 3]``\n    \"\"\"\n\n    def __init__(self,\n                 pipelines,\n                 size=-1,\n                 reader_name=None,\n                 auto_reset=False,\n                 fill_last_batch=None,\n                 dynamic_shape=False,\n                 last_batch_padded=False,\n                 last_batch_policy=LastBatchPolicy.FILL,\n                 prepare_first_batch=True):\n        super(DALIClassificationIterator, self).__init__(pipelines, [\"data\", \"label\"],\n                                                         size,\n                                                         reader_name=reader_name,\n                                                         auto_reset=auto_reset,\n                                                         fill_last_batch=fill_last_batch,\n                                                         dynamic_shape=dynamic_shape,\n                                                         last_batch_padded=last_batch_padded,\n                                                         last_batch_policy=last_batch_policy,\n                                                         prepare_first_batch=prepare_first_batch)\n\n\nclass TorchPythonFunction(ops.PythonFunctionBase):\n    schema_name = \"TorchPythonFunction\"\n    ops.register_cpu_op('TorchPythonFunction')\n    ops.register_gpu_op('TorchPythonFunction')\n\n    def _torch_stream_wrapper(self, function, *ins):\n        with torch.cuda.stream(self.stream):\n            out = function(*ins)\n        self.stream.synchronize()\n        return out\n\n    def torch_wrapper(self, batch_processing, function, device, *args):\n        func = function if device == 'cpu' else \\\n               lambda *ins: self._torch_stream_wrapper(function, *ins)\n        if batch_processing:\n            return ops.PythonFunction.function_wrapper_batch(func,\n                                                             self.num_outputs,\n                                                             torch.utils.dlpack.from_dlpack,\n                                                             torch.utils.dlpack.to_dlpack,\n                                                             *args)\n        else:\n            return ops.PythonFunction.function_wrapper_per_sample(func,\n                                                                  self.num_outputs,\n                                                                  torch_dlpack.from_dlpack,\n                                                                  torch_dlpack.to_dlpack,\n                                                                  *args)\n\n    def __call__(self, *inputs, **kwargs):\n        pipeline = Pipeline.current()\n        if pipeline is None:\n            Pipeline._raise_no_current_pipeline(\"TorchPythonFunction\")\n        if self.stream is None:\n            self.stream = torch.cuda.Stream(device=pipeline.device_id)\n        return super(TorchPythonFunction, self).__call__(*inputs, **kwargs)\n\n    def __init__(self, function, num_outputs=1, device='cpu', batch_processing=False, **kwargs):\n        self.stream = None\n        super(TorchPythonFunction, self).__init__(impl_name=\"DLTensorPythonFunctionImpl\",\n                                                  function=lambda *ins:\n                                                  self.torch_wrapper(batch_processing,\n                                                                     function, device,\n                                                                     *ins),\n                                                  num_outputs=num_outputs, device=device,\n                                                  batch_processing=batch_processing, **kwargs)\n\n\nops._wrap_op(TorchPythonFunction, \"fn\", __name__)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-04-07T16:40:34.588242Z","iopub.execute_input":"2023-04-07T16:40:34.588909Z","iopub.status.idle":"2023-04-07T16:40:34.636844Z","shell.execute_reply.started":"2023-04-07T16:40:34.588871Z","shell.execute_reply":"2023-04-07T16:40:34.635805Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Defining helper classes","metadata":{}},{"cell_type":"code","source":"class BoundingBox(object):\n    def __init__(self, top_left, bottom_right):\n        self.top_left = top_left\n        self.bottom_right = bottom_right","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:34.640424Z","iopub.execute_input":"2023-04-07T16:40:34.640836Z","iopub.status.idle":"2023-04-07T16:40:34.646293Z","shell.execute_reply.started":"2023-04-07T16:40:34.6408Z","shell.execute_reply":"2023-04-07T16:40:34.645189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Side(Enum):\n    LEFT = 1\n    RIGHT = 2","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:34.647864Z","iopub.execute_input":"2023-04-07T16:40:34.648331Z","iopub.status.idle":"2023-04-07T16:40:34.654601Z","shell.execute_reply.started":"2023-04-07T16:40:34.648297Z","shell.execute_reply":"2023-04-07T16:40:34.653623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Main class preprocessor\nNote: it has only one public method, and all others are private, but since the purpose of this notebook is to provide demo on the class functionality, those functions will be converted to public.","metadata":{}},{"cell_type":"code","source":"class Preprocessor(object):\n    \"\"\"Class that preprocess input dataset in dicom format for breast cancer recognition. \n    \"\"\"\n    def __init__(self, imgs_path: str, csv_path: str, save_dir: str, device: str):\n        \"\"\"\n        Args:\n            imgs_path (str): Path to images\n            csv_path (str): Path to csv metadata\n            save_dir (str): Path to save directory\n        \"\"\"\n        self.imgs_path = imgs_path\n        self.csv_path = csv_path\n        \n        # make saving directory\n        os.makedirs(save_dir, exist_ok=True)\n        self.save_dir = save_dir\n        self.device = device\n\n        \n        \n    def convert_dicom_to_j2k(self, file, save_folder):\n        \"\"\" Decodes dicom image to JPEG2000 if DICOM image is in JPEG2000 format\n        \"\"\"\n        patient = file.split('/')[-2]\n        image = file.split('/')[-1][:-4]\n        dcmfile = pydicom.dcmread(file)\n\n        if dcmfile.file_meta.TransferSyntaxUID == '1.2.840.10008.1.2.4.90':\n            with open(file, 'rb') as fp:\n                raw = DicomBytesIO(fp.read())\n                ds = pydicom.dcmread(raw)\n            offset = ds.PixelData.find(b\"\\x00\\x00\\x00\\x0C\")  #<---- the jpeg2000 header info we're looking for\n            hackedbitstream = bytearray()\n            hackedbitstream.extend(ds.PixelData[offset:])\n            with open(save_folder + f\"{patient}_{image}.jp2\", \"wb\") as binary_file:\n                binary_file.write(hackedbitstream)\n\n            \n    @pipeline_def\n    def j2k_decode_pipeline(self, j2kfiles):\n        \"\"\"DALI pipeline for decoding images\"\"\"\n        jpegs, _ = fn.readers.file(files=j2kfiles)\n        images = fn.experimental.decoders.image(jpegs, device='mixed', output_type=types.ANY_DATA, dtype=DALIDataType.UINT16)\n        return images\n    \n    \n    \n    def make_batches(self, data : list, batch_size : int):\n        \"\"\"Makes N batche from list of strings that represents file paths\"\"\"\n        retVal = []\n        batch = []        \n        for i, el in enumerate(data):\n\n            batch.append(el)\n\n            if (i + 1) % batch_size == 0:\n                retVal.append(batch) \n                batch = []\n\n        # append final non-full batch\n        if len(batch) != 0:\n            retVal.append(batch)\n            batch = []\n\n\n        return retVal\n    \n    \n    def process_all_j2k(self, paths,  windowing = True, shuffle = False, size = (512, 512), convertToOneSide = False, side = Side.LEFT):\n        \"\"\"Main function for processing all JPEG2000 images.\n        Args:\n            n (int): Number of images to process\n            windowing (bool): Windowing technique to be applied. Defaults to True.\n            shuffle (bool): Weather images will be shuffled. Defaults to False.\n            size (tuple): output image size. Defaults to (512, 512)\n            convertToOneSide (bool): Weather images will be converted to one side, either left or right. Defaults to False.\n            side (Side): Which side (left or right) are images oriented.\n        \"\"\"\n        \n        batches = self.make_batches(paths, 32)\n        J2K_FOLDER = \"/tmp/j2k/\"\n\n        for batch in tqdm(batches):\n            gc.collect()\n\n            os.makedirs(J2K_FOLDER, exist_ok=True)\n\n            _ = Parallel(n_jobs=2)(\n                delayed(self.convert_dicom_to_j2k)(img, save_folder=J2K_FOLDER)\n                for img in batch\n            )\n\n            j2kfiles = glob.glob(J2K_FOLDER + \"*.jp2\")\n\n            if not len(j2kfiles):\n                continue\n\n            pipe = self.j2k_decode_pipeline(j2kfiles, batch_size=1, num_threads=2, device_id=0, debug=True)\n            pipe.build()\n\n            for i, path in enumerate(j2kfiles):\n                patient, image = path.split('/')[-1][:-4].split('_')\n                dicom = pydicom.dcmread(f\"{self.imgs_path}/{patient}/{image}.dcm\")\n\n                out = pipe.run()\n\n                # Dali -> Torch\n                img = out[0][0]\n                img_torch = torch.empty(img.shape(), dtype=torch.int16, device=\"cuda\")\n                feed_ndarray(img, img_torch, cuda_stream=torch.cuda.current_stream(device=0))\n                img = img_torch.float()\n\n                # Scale, invert, apply windowing on GPU ! (this probably won't be done on GPU, because for cropping\n                # we need OpenCV)\n                if windowing:\n                    img = self.windowing(img, dicom)\n                    \n                min_, max_ = img.min(), img.max()\n                img = (img - min_) / (max_ - min_)\n\n                # resizing will be done after cropping\n                #if dims:\n                #    img = F.interpolate(img.view(1, 1, img.size(0), img.size(1)), (dims, dims), mode=\"bilinear\")[0, 0]\n\n                if dicom.PhotometricInterpretation == \"MONOCHROME1\":\n                    img = 1 - img\n\n                # Back to CPU\n                img = (img * 255).cpu().numpy().astype(np.uint8)\n\n\n                self.process_image(img, path, size = size, convertToOneSide = convertToOneSide, side = side)\n\n\n            shutil.rmtree(J2K_FOLDER)\n            gc.collect()\n        \n    def windowing(self, img, dicom):\n        \"\"\" This windowing function is designed for GPU processing.\"\"\"\n        center = dicom.WindowCenter        \n        width = dicom.WindowWidth\n        bits_stored = dicom.BitsStored\n        \n        # some images doesn't have this attribute\n        if hasattr(dicom, 'VOILUTFunction'):\n            function = dicom.VOILUTFunction\n        else:\n            function = 'LINEAR'\n        \n        # MultiValue instead of List\n        if isinstance(center, pydicom.multival.MultiValue):\n            center = torch.Tensor([center[0]]).to(device = 'cuda')\n        else:\n            center = torch.Tensor([center]).to(device = 'cuda')\n        \n        if isinstance(width, pydicom.multival.MultiValue):\n            \n            width = torch.Tensor([width[0]]).to(device = 'cuda')\n        else:\n            width = torch.Tensor([width]).to(device = 'cuda')\n        \n        y_range = torch.Tensor([2**bits_stored - 1]).to(device = 'cuda')\n        \n        if function == 'SIGMOID':\n            img = y_range / (1 + torch.exp(-4 * (img - center) / width))\n        else: # LINEAR\n            center -= 0.5\n            width -= 1\n            below = img <= (center - width / 2)\n            above = img > (center + width / 2)\n            between = torch.logical_and(~below, ~above)\n            img[below] = 0\n            img[above] = y_range\n            img[between] = ((img[between] - center) / width + 0.5) * y_range\n        return img\n    \n    # ------------------------------ MAIN FUNCTION --------------------------------------#\n    def process_images(self, n, windowing = True, shuffle = False, size = (512, 512), convertToOneSide = False, side = Side.LEFT):\n\n        \"\"\" Main function that does all the processing: windowing, shuffling, resizing, saving and converting images to one side unifomrly (left or right)\n        Args:\n            n (int): Number of images to process\n            windowing (bool): Windowing technique to be applied. Defaults to True.\n            shuffle (bool): Weather images will be shuffled. Defaults to False.\n            size (tuple): output image size. Defaults to (512, 512)\n            convertToOneSide (bool): Weather images will be converted to one side, either left or right. Defaults to False.\n            side (Side): Which side (left or right) are images oriented.\n        \"\"\"\n        all_paths = self.get_file_paths(self.imgs_path, self.csv_path, n, shuffle = shuffle)\n        \n        # IF CUDA is the selected device, than do processing with DALI. Process with CPU otherwise.\n        if self.device == torch.device('cuda'):\n            self.process_all_j2k(all_paths, windowing = windowing, size = size, convertToOneSide = convertToOneSide, side = side)\n        \n        # exclude all images that are already processed\n        processed_imgs_paths = glob.glob(f\"{self.save_dir}*.png\")\n        nonprocessed_paths = self.find_nonprocessed_imgs(all_paths, processed_imgs_paths)\n\n        _ = Parallel(n_jobs = 1)(delayed(self.process_image)(None, path = path, windowing = windowing, size = size,\n                                     convertToOneSide = convertToOneSide, side = side) for path in tqdm(nonprocessed_paths))\n    # ----------------------------------------------------------------------------------------#\n    \n    def process_image(self, img, path : str, windowing = True, size = (512, 512), convertToOneSide = False, side = Side.LEFT):\n\n        \"\"\" Process a single image given as numpy array.\n        Args:\n            img (numpy.ndarray): Input image.\n            img_name (string): Image name for saving\n            shuffle (bool): Weather images will be shuffled. Defaults to False.\n            size (tuple): output image size. Defaults to (512, 512)\n            convertToOneSide (bool): Weather images will be converted to one side, either left or right. Defaults to False.\n            side (Side): Which side (left or right) are images oriented.        \n        \"\"\"\n        # if image is None, it means that it should be loaded from path\n        # This is in case that image is loaded with CPU, so that this process of retrieving pixels from DICOM\n        # is also parallelized.\n\n        if img is None:\n            img = self.get_pixels(path, windowing = windowing)\n            # format image name in slightly different manner\n            image_id = path.split('/')[-1].split('.')[0]\n            patient_id = path.split('/')[-2]\n            img_name = f\"{patient_id}_{image_id}.png\"\n        else:\n            # format image name in slightly different manner\n            patient_id, image_id = path.split('/')[-1].split('.')[0].split('_')\n            img_name = f\"{patient_id}_{image_id}.png\"\n    \n\n        \n        # determine side before cropping\n        img_side = self.determine_side(img)\n\n        bbox, _ = self.get_bbox(img, 5)\n        img = self.crop_image(img, bbox)\n            \n        if convertToOneSide:\n            # if image is not on specified side\n            if img_side != side: \n                # Use Flip code 1 to flip horisontally\n                img = cv2.flip(img, 1)\n                \n                \n        img = cv2.resize(img, size)\n        \n        cv2.imwrite(os.path.join(self.save_dir, img_name), img)\n        return 0\n\n         \n    def get_file_paths(self, imgs_path: str, csv_path: str, n: int, shuffle = False):\n        \"\"\"_summary_\n\n        Args:\n            imgs_path (str): path to directory contatining images\n            csv_path (str): path to directory contatining csv data\n            n (int): number of images to load\n            shuffle (bool): Weather images will be shuffled. Defaults to False.\n\n        Returns:\n            List: list of paths to file\n        \"\"\"\n        \n        df = pd.read_csv(csv_path)\n        if shuffle:\n            data = np.random.permutation(df[['image_id', 'patient_id']].values)\n            img_id = data[:, 0]\n            patient_id = data[:, 1]\n        else:\n            data = df[['image_id', 'patient_id']].values\n            img_id = data[:, 0]\n            patient_id = data[:, 1]    \n        file_names = []\n        for i in range(n):\n            file_names.append(f'{imgs_path}/{str(patient_id[i])}/{str(img_id[i])}.dcm')\n\n        return file_names\n    \n    def find_nonprocessed_imgs(self, all_paths, processed_imgs_paths):\n        \"\"\" This function excludes paths of images that have already been processed with DALIi (All JPEG2000 images)\"\"\"\n        img_ids_processed = []\n        for p in processed_imgs_paths:\n            img_ids_processed.append(p.split('/')[-1].split('_')[-1].split('.')[0])\n        \n        n = len(all_paths)\n        for i, el in enumerate(reversed(all_paths)):\n            if el.split('/')[-1].split('.')[0] in img_ids_processed:\n                del all_paths[n - 1 - i]\n\n        return all_paths\n\n    def get_pixels(self, dcm_file : str, windowing = True):\n        \"\"\"  Gets image as numpy array in specified datatype from dicom format\n\n        Args:\n            dcm_file (str): Path to dcm file\n            windowing (bool): Windowing technique to be applied.. Defaults to True.\n            to_dtype (string): Output image format. Defaults to 'uint8'.\n\n        Raises:\n            Exception: For invalid arguments\n\n        Returns:\n            np.array: Image\n        \"\"\"\n\n        scan = pydicom.dcmread(dcm_file)\n\n        data = scan.pixel_array\n\n        # apply windowing to the original image to make contrast between soft and\n        # hard tissue more visible\n        if windowing:\n            data = apply_windowing(data, scan)\n\n        # normalize image to \n        if scan.PhotometricInterpretation == \"MONOCHROME1\":\n            data = np.amax(data) - data\n        else:\n            data = data - np.min(data)\n\n        if np.max(data) != 0:\n            data = data / np.max(data)\n\n        data=(data * 255).astype(dtype = 'uint8')\n\n        return data\n    \n    \n    def get_bbox(self, img, threshold = 5):\n        \"\"\"\n        img: input image\n        threashold: should be 5, experimentally derived\n        \"\"\"\n        if img.dtype == 'uint8':\n            ret, thresh = cv2.threshold(img, threshold, 255, 0)\n        else:\n            img = (255*img).astype(dtype = 'uint8')\n            ret, thresh = cv2.threshold(img, threshold, 255, 0)\n\n        contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n        contour = max(contours, key=cv2.contourArea)\n\n        contour = np.squeeze(contour)\n        bbox = BoundingBox(top_left = (min(contour[:, 0]), min(contour[:, 1])), bottom_right = (max(contour[:, 0]), max(contour[:, 1])))\n        return bbox, img # return bbox and image in uint8 representation\n    \n\n    def determine_side(self, img, threshold = 5):\n        \"\"\"\n        img: input image\n        threshold: for binirizing image, should be 5\n        Side is determined simply by finding more white side of the image.\n        \"\"\"\n\n        if img.dtype == 'uint8':\n            ret, thresh = cv2.threshold(img, threshold, 255, 0)\n        else:\n            img = (255*img).astype(dtype = 'uint8')\n            ret, thresh = cv2.threshold(img, threshold, 255, 0)\n\n        width = img.shape[1]\n        # take all columns up to half image (in width), sumarize them and compare with other half\n        if sum(sum(thresh[:, :width // 2])) > sum(sum(thresh[:, width // 2:])): \n            return Side.LEFT\n        else:\n            return Side.RIGHT\n    \n    \n    def crop_image(self, img, bbox):\n        \"\"\"\n        img: image to crop\n        bbox: bounding box that defines new image\n        \"\"\"\n        return img[bbox.top_left[1]:bbox.bottom_right[1], bbox.top_left[0]:bbox.bottom_right[0]]","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:34.656509Z","iopub.execute_input":"2023-04-07T16:40:34.657212Z","iopub.status.idle":"2023-04-07T16:40:34.704424Z","shell.execute_reply.started":"2023-04-07T16:40:34.657178Z","shell.execute_reply":"2023-04-07T16:40:34.703387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessor demo\nIn this section will be demonstrated some possibilites of this Preprocessor class.","metadata":{}},{"cell_type":"code","source":"device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f'Current device is {device}')","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:34.707193Z","iopub.execute_input":"2023-04-07T16:40:34.707496Z","iopub.status.idle":"2023-04-07T16:40:34.775008Z","shell.execute_reply.started":"2023-04-07T16:40:34.707459Z","shell.execute_reply":"2023-04-07T16:40:34.773794Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"csv_path = '/kaggle/input/rsna-breast-cancer-detection/train.csv'\nimgs_path = '/kaggle/input/rsna-breast-cancer-detection/train_images'\nsave_dir = '/kaggle/working/temp/'\ndf = pd.read_csv(csv_path)\n\npreprocessor = Preprocessor(imgs_path, csv_path, save_dir, device)\npaths = preprocessor.get_file_paths(imgs_path, csv_path, n = 5, shuffle = False)\n\nf, ax = plt.subplots(5, 4, figsize = (25, 25))\nimgs = []\nfor i, path in enumerate(paths):\n    img_windowed = preprocessor.get_pixels(path, windowing = True)\n    plt.subplot(5, 4, 4*i + 1)\n    plt.title('With windowing')\n    plt.imshow(img_windowed, cmap = 'gray')\n    \n    img = preprocessor.get_pixels(path, windowing = False)\n    plt.subplot(5, 4, 4*i + 2)\n    plt.imshow(img, cmap = 'gray')\n    plt.title('Without windowing')\n    \n    \n    bbox, img = preprocessor.get_bbox(img_windowed, 5)\n    # convert to colored and plot bounding box\n    img = np.repeat(img[:, :, np.newaxis], 3, axis=2)\n    # plot bounding box with green\n    boxed_img = cv2.rectangle(img, bbox.top_left, bbox.bottom_right, (0, 255, 0), 10)\n    plt.subplot(5, 4, 4*i + 3)\n    plt.imshow(boxed_img)\n    plt.title('With bounding box')\n    \n    \n    final_img = preprocessor.crop_image(img_windowed, bbox)\n    plt.subplot(5, 4, 4*i + 4)\n    plt.imshow(final_img, cmap = 'gray')\n    plt.title('Final cropped image')\n    \n    side = preprocessor.determine_side(img_windowed)\n    print(side)","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:40:34.776657Z","iopub.execute_input":"2023-04-07T16:40:34.777807Z","iopub.status.idle":"2023-04-07T16:41:25.131511Z","shell.execute_reply.started":"2023-04-07T16:40:34.777768Z","shell.execute_reply":"2023-04-07T16:41:25.130646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Making a Dataset","metadata":{}},{"cell_type":"code","source":"csv_path = '/kaggle/input/rsna-breast-cancer-detection/train.csv'\nimgs_path = '/kaggle/input/rsna-breast-cancer-detection/train_images'\nsave_dir = '/kaggle/working/temp/'\ndf = pd.read_csv(csv_path)\n\npreprocessor = Preprocessor(imgs_path, csv_path, save_dir, device)\n\nstart = time.time()\n\nn = 1000 # this is just for sample, n shoudl be len(df) for full dataset\nimgs = preprocessor.process_images(n = len(df), windowing = True, shuffle = False, size = (512, 512), convertToOneSide = False)\nduration = time.time() - start\nprint(f\"Processing of {n} images lasted {duration}\")\nprint(f\"Processing of 54000 images from training data should last {duration*54/60} minutes, or {duration*54/60/60} hours\")","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:54:59.079589Z","iopub.execute_input":"2023-04-07T16:54:59.080617Z","iopub.status.idle":"2023-04-08T00:29:34.096832Z","shell.execute_reply.started":"2023-04-07T16:54:59.080579Z","shell.execute_reply":"2023-04-08T00:29:34.095082Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Plot processed image examples","metadata":{}},{"cell_type":"code","source":"processed_imgs_paths = glob.glob(f\"{save_dir}*.png\")\n\nf, ax = plt.subplots(5, 2, figsize = (20, 50))\ni = 1\nfor img_path in processed_imgs_paths:\n    img = cv2.imread(img_path)\n    plt.subplot(5, 2, i)\n    plt.imshow(img, cmap = 'gray')\n    i += 1\n    if i > 10:\n        break","metadata":{"execution":{"iopub.status.busy":"2023-04-08T16:03:11.447464Z","iopub.execute_input":"2023-04-08T16:03:11.447887Z","iopub.status.idle":"2023-04-08T16:03:11.52699Z","shell.execute_reply.started":"2023-04-08T16:03:11.447804Z","shell.execute_reply":"2023-04-08T16:03:11.525527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusion\nThis Preprocessor can easily be used when submiting the solution to the competition, because it runs very fastly","metadata":{}},{"cell_type":"code","source":"x = cv2.imread(img_path)\nlen(x)","metadata":{"execution":{"iopub.status.busy":"2023-04-07T16:53:48.055031Z","iopub.execute_input":"2023-04-07T16:53:48.05541Z","iopub.status.idle":"2023-04-07T16:53:48.064809Z","shell.execute_reply.started":"2023-04-07T16:53:48.055379Z","shell.execute_reply":"2023-04-07T16:53:48.06354Z"},"trusted":true},"execution_count":null,"outputs":[]}]}