{"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":"code","source":"import os\nos.environ['CUDA_MODULE_LOADING']='LAZY'\n\ntry: \n    import dicomsdl\n    \nexcept:\n    !pip install /kaggle/input/rsna-2022-whl/pylibjpeg-1.4.0-py3-none-any.whl\n    !pip install /kaggle/input/rsna-2022-whl/python_gdcm-3.0.15-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl\n    !pip install /kaggle/input/dicomsdl-offline-installer/dicomsdl-0.109.1-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl\n    !cp /kaggle/input/easy-load-the-image-with-nvjpeg2000/nvjpeg2k.so ./\n    \nimport torch\nprint(torch.__version__)\n\nimport sys\nsys.path.append('/kaggle/input/rsna-breast-mammography-00')\n\nprint('install ok')","metadata":{"execution":{"iopub.status.busy":"2023-02-09T19:00:40.002011Z","iopub.execute_input":"2023-02-09T19:00:40.002491Z","iopub.status.idle":"2023-02-09T19:02:13.449162Z","shell.execute_reply.started":"2023-02-09T19:00:40.002412Z","shell.execute_reply":"2023-02-09T19:02:13.447272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here i opensource my dicom converter which is one of the fastest available after combining the hard works of various kagglers. ","metadata":{}},{"cell_type":"code","source":"#this is dicom_reader_v3a.py\n\nimport pandas as pd\nimport numpy as np\nimport cv2\n\nfrom timeit import default_timer as timer\nfrom joblib import Parallel, delayed\nfrom glob import glob\n##from tqdm import tqdm\nfrom tqdm.notebook import tqdm\n\nimport os\nimport sys\n\nimport pydicom\n\nimport dicomsdl\nimport nvjpeg2k\n\n###########################################################################################\n#from pydicom.pixel_data_handlers.util import apply_voi_lut\n\nfrom typing import (\n    Dict, Optional, Union, List, Tuple, TYPE_CHECKING, cast, Iterable,\n    ByteString\n)\nfrom pydicom.valuerep import VR\n\ndef apply_voi_lut(\n    arr: \"np.ndarray\",\n    ds: \"Dataset\",\n    index: int = 0,\n    prefer_lut: bool = True\n) -> \"np.ndarray\":\n    \"\"\"Apply a VOI lookup table or windowing operation to `arr`.\n\n    .. versionadded:: 1.4\n\n    .. versionchanged:: 2.1\n\n        Added the `prefer_lut` keyword parameter\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the VOI LUT or windowing operation\n        to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,3010) *VOI LUT Sequence* is present then returns an array\n        of ``np.uint8`` or ``np.uint16``, depending on the 3rd value of\n        (0028,3002) *LUT Descriptor*. If (0028,1050) *Window Center* and\n        (0028,1051) *Window Width* are present then returns an array of\n        ``np.float64``. If neither are present then `arr` will be returned\n        unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n    prefer_lut : bool\n        When the VOI LUT Module contains both *Window Width*/*Window Center*\n        and *VOI LUT Sequence*, if ``True`` (default) then apply the VOI LUT,\n        otherwise apply the windowing operation.\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied VOI LUT or windowing operation.\n\n    Notes\n    -----\n    When the dataset requires a modality LUT or rescale operation as part of\n    the Modality LUT module then that must be applied before any windowing\n    operation.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_voi`\n    :func:`~pydicom.pixel_data_handlers.util.apply_windowing`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    valid_voi = False\n    if ds.get('VOILUTSequence'):\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\n    if valid_voi and valid_windowing:\n        if prefer_lut:\n            return apply_voi(arr, ds, index)\n\n        return apply_windowing(arr, ds, index)\n\n    if valid_voi:\n        return apply_voi(arr, ds, index)\n\n    if valid_windowing:\n        return apply_windowing(arr, ds, index)\n\n    return arr\n\n\ndef apply_voi(\n    arr: \"np.ndarray\", ds: \"Dataset\", index: int = 0\n) -> \"np.ndarray\":\n    \"\"\"Apply a VOI lookup table to `arr`.\n\n    .. versionadded:: 2.1\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the VOI LUT to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,3010) *VOI LUT Sequence* is present then returns an array\n        of ``np.uint8`` or ``np.uint16``, depending on the 3rd value of\n        (0028,3002) *LUT Descriptor*, otherwise `arr` will be returned\n        unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied VOI LUT.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_windowing`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    if not ds.get('VOILUTSequence'):\n        return arr\n\n    if not np.issubdtype(arr.dtype, np.integer):\n        print(#warnings.warn\n            \"Applying a VOI LUT on a float input array may give \"\n            \"incorrect results\"\n        )\n\n    # VOI LUT Sequence contains one or more items\n    item = cast(List[\"Dataset\"], ds.VOILUTSequence)[index]\n    lut_descriptor = cast(List[int], item.LUTDescriptor)\n    nr_entries = lut_descriptor[0] or 2**16\n    first_map = lut_descriptor[1]\n\n    # PS3.3 C.8.11.3.1.5: may be 8, 10-16\n    nominal_depth = lut_descriptor[2]\n    if nominal_depth in list(range(10, 17)):\n        dtype = 'uint16'\n    elif nominal_depth == 8:\n        dtype = 'uint8'\n    else:\n        raise NotImplementedError(\n            f\"'{nominal_depth}' bits per LUT entry is not supported\"\n        )\n\n    # Ambiguous VR, US or OW\n    unc_data: Iterable[int]\n    if item['LUTData'].VR == VR.OW:\n        endianness = '<' if ds.is_little_endian else '>'\n        unpack_fmt = f'{endianness}{nr_entries}H'\n        unc_data = unpack(unpack_fmt, cast(bytes, item.LUTData))\n    else:\n        unc_data = cast(List[int], item.LUTData)\n\n    lut_data: \"np.ndarray\" = np.asarray(unc_data, dtype=dtype)\n\n    # IVs < `first_map` get set to first LUT entry (i.e. index 0)\n    clipped_iv = np.zeros(arr.shape, dtype=dtype)\n    # IVs >= `first_map` are mapped by the VOI LUT\n    # `first_map` may be negative, positive or 0\n    mapped_pixels = arr >= first_map\n    clipped_iv[mapped_pixels] = arr[mapped_pixels] - first_map\n    # IVs > number of entries get set to last entry\n    np.clip(clipped_iv, 0, nr_entries - 1, out=clipped_iv)\n\n    return cast(\"np.ndarray\", lut_data[clipped_iv])\n\n\ndef apply_windowing(\n    arr: \"np.ndarray\", ds: \"Dataset\", index: int = 0\n) -> \"np.ndarray\":\n    \"\"\"Apply a windowing operation to `arr`.\n\n    .. versionadded:: 2.1\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the windowing operation to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,1050) *Window Center* and (0028,1051) *Window Width* are\n        present then returns an array of ``np.float64``, otherwise `arr` will\n        be returned unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied windowing operation.\n\n    Notes\n    -----\n    When the dataset requires a modality LUT or rescale operation as part of\n    the Modality LUT module then that must be applied before any windowing\n    operation.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_voi`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    if \"WindowWidth\" not in ds and \"WindowCenter\" not in ds:\n        return arr\n\n    if ds.PhotometricInterpretation not in ['MONOCHROME1', 'MONOCHROME2']:\n        raise ValueError(\n            \"When performing a windowing operation only 'MONOCHROME1' and \"\n            \"'MONOCHROME2' are allowed for (0028,0004) Photometric \"\n            \"Interpretation\"\n        )\n\n    # May be LINEAR (default), LINEAR_EXACT, SIGMOID or not present, VM 1\n    voi_func = cast(str, getattr(ds, 'VOILUTFunction', 'LINEAR')).upper()\n    # VR DS, VM 1-n\n    elem = ds['WindowCenter']\n    center = (\n        cast(List[float], elem.value)[index] if elem.VM > 1 else elem.value\n    )\n    center = cast(float, center)\n    elem = ds['WindowWidth']\n    width = cast(List[float], elem.value)[index] if elem.VM > 1 else elem.value\n    width = cast(float, width)\n\n    # The output range depends on whether or not a modality LUT or rescale\n    #   operation has been applied\n    ds.BitsStored = cast(int, ds.BitsStored)\n    y_min: float\n    y_max: float\n    if ds.get('ModalityLUTSequence'):\n        # Unsigned - see PS3.3 C.11.1.1.1\n        y_min = 0\n        item = cast(List[\"Dataset\"], ds.ModalityLUTSequence)[0]\n        bit_depth = cast(List[int], item.LUTDescriptor)[2]\n        y_max = 2**bit_depth - 1\n    elif ds.PixelRepresentation == 0:\n        # Unsigned\n        y_min = 0\n        y_max = 2**ds.BitsStored - 1\n    else:\n        # Signed\n        y_min = -2**(ds.BitsStored - 1)\n        y_max = 2**(ds.BitsStored - 1) - 1\n\n    slope = ds.get('RescaleSlope', None)\n    intercept = ds.get('RescaleIntercept', None)\n    if slope is not None and intercept is not None:\n        ds.RescaleSlope = cast(float, ds.RescaleSlope)\n        ds.RescaleIntercept = cast(float, ds.RescaleIntercept)\n        # Otherwise its the actual data range\n        y_min = y_min * ds.RescaleSlope + ds.RescaleIntercept\n        y_max = y_max * ds.RescaleSlope + ds.RescaleIntercept\n\n    y_range = y_max - y_min\n    arr = arr.astype('float32')\n    #arr = arr.astype('float64')\n\n    if voi_func in ['LINEAR', 'LINEAR_EXACT']:\n        # PS3.3 C.11.2.1.2.1 and C.11.2.1.3.2\n        if voi_func == 'LINEAR':\n            if width < 1:\n                raise ValueError(\n                    \"The (0028,1051) Window Width must be greater than or \"\n                    \"equal to 1 for a 'LINEAR' windowing operation\"\n                )\n            center -= 0.5\n            width -= 1\n        elif width <= 0:\n            raise ValueError(\n                \"The (0028,1051) Window Width must be greater than 0 \"\n                \"for a 'LINEAR_EXACT' windowing operation\"\n            )\n\n        below = arr <= (center - width / 2)\n        above = arr > (center + width / 2)\n        between = np.logical_and(~below, ~above)\n\n        arr[below] = y_min\n        arr[above] = y_max\n        if between.any():\n            arr[between] = (\n                ((arr[between] - center) / width + 0.5) * y_range + y_min\n            )\n    elif voi_func == 'SIGMOID':\n        # PS3.3 C.11.2.1.3.1\n        if width <= 0:\n            raise ValueError(\n                \"The (0028,1051) Window Width must be greater than 0 \"\n                \"for a 'SIGMOID' windowing operation\"\n            )\n\n        arr = y_range / (1 + np.exp(-4 * (arr - center) / width)) + y_min\n    else:\n        raise ValueError(\n            f\"Unsupported (0028,1056) VOI LUT Function value '{voi_func}'\"\n        )\n\n    return arr\n###########################################################################################\n\ndef time_to_str(t, mode='min'):\n    if mode=='min':\n        t  = int(t)/60\n        hr = t//60\n        min = t%60\n        return '%2d hr %02d min'%(hr,min)\n\n    elif mode=='sec':\n        t   = int(t)\n        min = t//60\n        sec = t%60\n        return '%2d min %02d sec'%(min,sec)\n\n    else:\n        raise NotImplementedError\n\n\n###########################################################################################\ndef read_image(df, image_dir):\n    image = []\n    for t,d in df.iterrows():\n        image_file = f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png'\n        m = cv2.imread(image_file,cv2.IMREAD_ANYDEPTH)\n        image.append(m)\n    return image\n\ndef make_transfer_syntax_uid(df, dcm_dir):\n    machine_id_to_transfer = {}\n    machine_id = df.machine_id.unique()\n    for i in machine_id:\n        d = df[df.machine_id == i].iloc[0]\n        f = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dicom = pydicom.dcmread(f)\n        machine_id_to_transfer[i] = dicom.file_meta.TransferSyntaxUID\n    return machine_id_to_transfer\n\ndef normalised_to_8bit(image, photometric_interpretation):\n    xmin = image.min()\n    xmax = image.max()\n\n    norm = np.empty_like(image, dtype=np.uint8)\n    dicomsdl.util.convert_to_uint8(image, norm, xmin, xmax)\n    if photometric_interpretation == 'MONOCHROME1':\n        norm = 255 - norm\n    return norm\n\ndef resize_image_to_height(image, image_height):\n    h, w = image.shape[:2]\n    s = image_height/h\n    if image_height!=h:\n        image = cv2.resize(image, dsize=None, fx=s, fy=s, interpolation=cv2.INTER_LINEAR)\n    return image\n\n#----------------------------------------------------------------\n# dicomsdl reader\ndef dicomsdl_to_numpy_image(ds, index=0):\n    # https://stackoverflow.com/questions/44659924/returning-numpy-arrays-via-pybind11\n    info = ds.getPixelDataInfo()\n    if info['SamplesPerPixel'] != 1:\n        raise RuntimeError('SamplesPerPixel != 1')\n\n    shape = [info['Rows'], info['Cols']]\n    dtype = info['dtype']\n    outarr = np.empty(shape, dtype=dtype)\n    ds.copyFrameData(index, outarr)\n    return outarr\n\ndef dicomsdl_parallel_process(d, dcm_dir, image_dir, image_height, is_voi_lut):\n    dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n    ds = dicomsdl.open(dcm_file)\n    image = dicomsdl_to_numpy_image(ds)\n    image = resize_image_to_height(image, image_height)\n\n    if is_voi_lut:\n        dc = pydicom.dcmread(dcm_file)\n        image = apply_voi_lut(image, dc)\n        image = image.astype(np.float32)\n    image = normalised_to_8bit(image, ds.PhotometricInterpretation)  # +1\n\n    # save as png\n    os.makedirs(f'{image_dir}/{d.machine_id}/{d.patient_id}', exist_ok=True)\n    cv2.imwrite(f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png', image)\n\ndef process_non_j2k(df, dcm_dir, image_dir, image_height, n_jobs, is_voi_lut=True):\n    #https://stackoverflow.com/questions/56659294/does-joblib-parallel-keep-the-original-order-of-data-passed\n    #Parallel(n_jobs=2, backend='multiprocessing')(\n    Parallel(n_jobs=n_jobs)(\n        delayed(dicomsdl_parallel_process)(d, dcm_dir, image_dir, image_height, is_voi_lut)\n        for t,d in tqdm(df.iterrows())\n    )\n\n#----------------------------------------------------------------\n# nvjpeg2k reader\n\n'''\nTransferSyntaxUID\n1.2.840.10008.1.2.4.70 = JPEG Lossless, Nonhierarchical, First- Order Prediction (Processes 14)\n1.2.840.10008.1.2.4.90 = JPEG 2000 Image Compression (Lossless Only)\n'''\nj2k_decoder = nvjpeg2k.Decoder()\n\ndef process_j2k(df, dcm_dir, image_dir, image_height, is_voi_lut=True):\n    for t, d in tqdm(df.iterrows()):\n        dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dc = pydicom.dcmread(dcm_file)\n        offset = dc.PixelData.find(b'\\x00\\x00\\x00\\x0C')\n        jpeg_stream = bytearray(dc.PixelData[offset:])\n        image = j2k_decoder.decode(jpeg_stream)\n        image = resize_image_to_height(image, image_height)\n        #     jpeg_stream = bytearray(ds.PixelData[offset:])\n        #     jpeg_stream = np.array(bytearray(ds.PixelData[offset:]),np.uint8)\n\n        if is_voi_lut:\n            image = apply_voi_lut(image, dc)\n            image = image.astype(np.float32)\n        image = normalised_to_8bit(image, dc.PhotometricInterpretation)\n\n        # save as png\n        os.makedirs(f'{image_dir}/{d.machine_id}/{d.patient_id}', exist_ok=True)\n        cv2.imwrite(f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png', image)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-02-09T19:02:13.451395Z","iopub.execute_input":"2023-02-09T19:02:13.452266Z","iopub.status.idle":"2023-02-09T19:02:16.46449Z","shell.execute_reply.started":"2023-02-09T19:02:13.452226Z","shell.execute_reply":"2023-02-09T19:02:16.463412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#this is dicom_reader_v3a.py\n\nimport pandas as pd\nimport numpy as np\nimport cv2\n\nfrom timeit import default_timer as timer\nfrom joblib import Parallel, delayed\nfrom glob import glob\n##from tqdm import tqdm\nfrom tqdm.notebook import tqdm\n\nimport os\nimport sys\n\nimport pydicom\n\nimport dicomsdl\nimport nvjpeg2k\n\n###########################################################################################\n#from pydicom.pixel_data_handlers.util import apply_voi_lut\n\nfrom typing import (\n    Dict, Optional, Union, List, Tuple, TYPE_CHECKING, cast, Iterable,\n    ByteString\n)\nfrom pydicom.valuerep import VR\n\ndef apply_voi_lut(\n    arr: \"np.ndarray\",\n    ds: \"Dataset\",\n    index: int = 0,\n    prefer_lut: bool = True\n) -> \"np.ndarray\":\n    \"\"\"Apply a VOI lookup table or windowing operation to `arr`.\n\n    .. versionadded:: 1.4\n\n    .. versionchanged:: 2.1\n\n        Added the `prefer_lut` keyword parameter\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the VOI LUT or windowing operation\n        to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,3010) *VOI LUT Sequence* is present then returns an array\n        of ``np.uint8`` or ``np.uint16``, depending on the 3rd value of\n        (0028,3002) *LUT Descriptor*. If (0028,1050) *Window Center* and\n        (0028,1051) *Window Width* are present then returns an array of\n        ``np.float64``. If neither are present then `arr` will be returned\n        unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n    prefer_lut : bool\n        When the VOI LUT Module contains both *Window Width*/*Window Center*\n        and *VOI LUT Sequence*, if ``True`` (default) then apply the VOI LUT,\n        otherwise apply the windowing operation.\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied VOI LUT or windowing operation.\n\n    Notes\n    -----\n    When the dataset requires a modality LUT or rescale operation as part of\n    the Modality LUT module then that must be applied before any windowing\n    operation.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_voi`\n    :func:`~pydicom.pixel_data_handlers.util.apply_windowing`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    valid_voi = False\n    if ds.get('VOILUTSequence'):\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\n    if valid_voi and valid_windowing:\n        if prefer_lut:\n            return apply_voi(arr, ds, index)\n\n        return apply_windowing(arr, ds, index)\n\n    if valid_voi:\n        return apply_voi(arr, ds, index)\n\n    if valid_windowing:\n        return apply_windowing(arr, ds, index)\n\n    return arr\n\n\ndef apply_voi(\n    arr: \"np.ndarray\", ds: \"Dataset\", index: int = 0\n) -> \"np.ndarray\":\n    \"\"\"Apply a VOI lookup table to `arr`.\n\n    .. versionadded:: 2.1\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the VOI LUT to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,3010) *VOI LUT Sequence* is present then returns an array\n        of ``np.uint8`` or ``np.uint16``, depending on the 3rd value of\n        (0028,3002) *LUT Descriptor*, otherwise `arr` will be returned\n        unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied VOI LUT.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_windowing`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    if not ds.get('VOILUTSequence'):\n        return arr\n\n    if not np.issubdtype(arr.dtype, np.integer):\n        print(#warnings.warn\n            \"Applying a VOI LUT on a float input array may give \"\n            \"incorrect results\"\n        )\n\n    # VOI LUT Sequence contains one or more items\n    item = cast(List[\"Dataset\"], ds.VOILUTSequence)[index]\n    lut_descriptor = cast(List[int], item.LUTDescriptor)\n    nr_entries = lut_descriptor[0] or 2**16\n    first_map = lut_descriptor[1]\n\n    # PS3.3 C.8.11.3.1.5: may be 8, 10-16\n    nominal_depth = lut_descriptor[2]\n    if nominal_depth in list(range(10, 17)):\n        dtype = 'uint16'\n    elif nominal_depth == 8:\n        dtype = 'uint8'\n    else:\n        raise NotImplementedError(\n            f\"'{nominal_depth}' bits per LUT entry is not supported\"\n        )\n\n    # Ambiguous VR, US or OW\n    unc_data: Iterable[int]\n    if item['LUTData'].VR == VR.OW:\n        endianness = '<' if ds.is_little_endian else '>'\n        unpack_fmt = f'{endianness}{nr_entries}H'\n        unc_data = unpack(unpack_fmt, cast(bytes, item.LUTData))\n    else:\n        unc_data = cast(List[int], item.LUTData)\n\n    lut_data: \"np.ndarray\" = np.asarray(unc_data, dtype=dtype)\n\n    # IVs < `first_map` get set to first LUT entry (i.e. index 0)\n    clipped_iv = np.zeros(arr.shape, dtype=dtype)\n    # IVs >= `first_map` are mapped by the VOI LUT\n    # `first_map` may be negative, positive or 0\n    mapped_pixels = arr >= first_map\n    clipped_iv[mapped_pixels] = arr[mapped_pixels] - first_map\n    # IVs > number of entries get set to last entry\n    np.clip(clipped_iv, 0, nr_entries - 1, out=clipped_iv)\n\n    return cast(\"np.ndarray\", lut_data[clipped_iv])\n\n\ndef apply_windowing(\n    arr: \"np.ndarray\", ds: \"Dataset\", index: int = 0\n) -> \"np.ndarray\":\n    \"\"\"Apply a windowing operation to `arr`.\n\n    .. versionadded:: 2.1\n\n    Parameters\n    ----------\n    arr : numpy.ndarray\n        The :class:`~numpy.ndarray` to apply the windowing operation to.\n    ds : dataset.Dataset\n        A dataset containing a :dcm:`VOI LUT Module<part03/sect_C.11.2.html>`.\n        If (0028,1050) *Window Center* and (0028,1051) *Window Width* are\n        present then returns an array of ``np.float64``, otherwise `arr` will\n        be returned unchanged.\n    index : int, optional\n        When the VOI LUT Module contains multiple alternative views, this is\n        the index of the view to return (default ``0``).\n\n    Returns\n    -------\n    numpy.ndarray\n        An array with applied windowing operation.\n\n    Notes\n    -----\n    When the dataset requires a modality LUT or rescale operation as part of\n    the Modality LUT module then that must be applied before any windowing\n    operation.\n\n    See Also\n    --------\n    :func:`~pydicom.pixel_data_handlers.util.apply_modality_lut`\n    :func:`~pydicom.pixel_data_handlers.util.apply_voi`\n\n    References\n    ----------\n    * DICOM Standard, Part 3, :dcm:`Annex C.11.2\n      <part03/sect_C.11.html#sect_C.11.2>`\n    * DICOM Standard, Part 3, :dcm:`Annex C.8.11.3.1.5\n      <part03/sect_C.8.11.3.html#sect_C.8.11.3.1.5>`\n    * DICOM Standard, Part 4, :dcm:`Annex N.2.1.1\n      <part04/sect_N.2.html#sect_N.2.1.1>`\n    \"\"\"\n    if \"WindowWidth\" not in ds and \"WindowCenter\" not in ds:\n        return arr\n\n    if ds.PhotometricInterpretation not in ['MONOCHROME1', 'MONOCHROME2']:\n        raise ValueError(\n            \"When performing a windowing operation only 'MONOCHROME1' and \"\n            \"'MONOCHROME2' are allowed for (0028,0004) Photometric \"\n            \"Interpretation\"\n        )\n\n    # May be LINEAR (default), LINEAR_EXACT, SIGMOID or not present, VM 1\n    voi_func = cast(str, getattr(ds, 'VOILUTFunction', 'LINEAR')).upper()\n    # VR DS, VM 1-n\n    elem = ds['WindowCenter']\n    center = (\n        cast(List[float], elem.value)[index] if elem.VM > 1 else elem.value\n    )\n    center = cast(float, center)\n    elem = ds['WindowWidth']\n    width = cast(List[float], elem.value)[index] if elem.VM > 1 else elem.value\n    width = cast(float, width)\n\n    # The output range depends on whether or not a modality LUT or rescale\n    #   operation has been applied\n    ds.BitsStored = cast(int, ds.BitsStored)\n    y_min: float\n    y_max: float\n    if ds.get('ModalityLUTSequence'):\n        # Unsigned - see PS3.3 C.11.1.1.1\n        y_min = 0\n        item = cast(List[\"Dataset\"], ds.ModalityLUTSequence)[0]\n        bit_depth = cast(List[int], item.LUTDescriptor)[2]\n        y_max = 2**bit_depth - 1\n    elif ds.PixelRepresentation == 0:\n        # Unsigned\n        y_min = 0\n        y_max = 2**ds.BitsStored - 1\n    else:\n        # Signed\n        y_min = -2**(ds.BitsStored - 1)\n        y_max = 2**(ds.BitsStored - 1) - 1\n\n    slope = ds.get('RescaleSlope', None)\n    intercept = ds.get('RescaleIntercept', None)\n    if slope is not None and intercept is not None:\n        ds.RescaleSlope = cast(float, ds.RescaleSlope)\n        ds.RescaleIntercept = cast(float, ds.RescaleIntercept)\n        # Otherwise its the actual data range\n        y_min = y_min * ds.RescaleSlope + ds.RescaleIntercept\n        y_max = y_max * ds.RescaleSlope + ds.RescaleIntercept\n\n    y_range = y_max - y_min\n    arr = arr.astype('float32')\n    #arr = arr.astype('float64')\n\n    if voi_func in ['LINEAR', 'LINEAR_EXACT']:\n        # PS3.3 C.11.2.1.2.1 and C.11.2.1.3.2\n        if voi_func == 'LINEAR':\n            if width < 1:\n                raise ValueError(\n                    \"The (0028,1051) Window Width must be greater than or \"\n                    \"equal to 1 for a 'LINEAR' windowing operation\"\n                )\n            center -= 0.5\n            width -= 1\n        elif width <= 0:\n            raise ValueError(\n                \"The (0028,1051) Window Width must be greater than 0 \"\n                \"for a 'LINEAR_EXACT' windowing operation\"\n            )\n\n        below = arr <= (center - width / 2)\n        above = arr > (center + width / 2)\n        between = np.logical_and(~below, ~above)\n\n        arr[below] = y_min\n        arr[above] = y_max\n        if between.any():\n            arr[between] = (\n                ((arr[between] - center) / width + 0.5) * y_range + y_min\n            )\n    elif voi_func == 'SIGMOID':\n        # PS3.3 C.11.2.1.3.1\n        if width <= 0:\n            raise ValueError(\n                \"The (0028,1051) Window Width must be greater than 0 \"\n                \"for a 'SIGMOID' windowing operation\"\n            )\n\n        arr = y_range / (1 + np.exp(-4 * (arr - center) / width)) + y_min\n    else:\n        raise ValueError(\n            f\"Unsupported (0028,1056) VOI LUT Function value '{voi_func}'\"\n        )\n\n    return arr\n###########################################################################################\n\ndef time_to_str(t, mode='min'):\n    if mode=='min':\n        t  = int(t)/60\n        hr = t//60\n        min = t%60\n        return '%2d hr %02d min'%(hr,min)\n\n    elif mode=='sec':\n        t   = int(t)\n    \n        min = t//60\n        sec = t%60\n        return '%2d min %02d sec'%(min,sec)\n\n    else:\n        raise NotImplementedError\n    \ndef _binarize(img, threshold):\n    return (img > threshold).astype(np.uint8)\n\n# Get contour points of the breast\ndef _extract_contour(bin_img):\n    contours, _ = cv2.findContours(\n        bin_img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n    contour = max(contours, key=cv2.contourArea)\n    return contour\n\n\n# Set to background pixels of the image to zero\ndef _erase_background(img, contour):\n    mask = np.zeros(img.shape, np.uint8)\n    cv2.drawContours(mask, [contour], -1, 255, cv2.FILLED)\n    output = cv2.bitwise_and(img, mask)\n    return output\n    \n # Crop the useless background of the image\ndef _crop(img):\n    bin_img = _binarize(img, threshold=5)\n    contour = _extract_contour(bin_img)\n    img = _erase_background(img, contour)\n    x1, x2 = np.min(contour[:, :, 0]), np.max(contour[:, :, 0])\n    y1, y2 = np.min(contour[:, :, 1]), np.max(contour[:, :, 1])\n    x1, x2 = int(0.99 * x1), int(1.01 * x2)\n    y1, y2 = int(0.99 * y1), int(1.01 * y2)\n    return img[y1:y2, x1:x2]   \n\ndef img2roi(img):\n    img_ori = img.copy()\n    # Binarize the image\n#     img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n    # img = img.astype(np.uint8)\n    bin_img = cv2.threshold(img, 20, 255, cv2.THRESH_BINARY)[1]\n    # print('aaa')\n    # Make contours around the binarized image, keep only the largest contour\n    contours, _ = cv2.findContours(bin_img, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n    contour = max(contours, key=cv2.contourArea)\n\n    # Find ROI from largest contour\n    ys = contour.squeeze()[:, 0]\n    xs = contour.squeeze()[:, 1]\n    roi =  img_ori[np.min(xs):np.max(xs), np.min(ys):np.max(ys)]\n    return roi\n\n###########################################################################################\ndef read_image(df, image_dir):\n    image = []\n    for t,d in df.iterrows():\n        image_file = f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png'\n        m = cv2.imread(image_file,cv2.IMREAD_ANYDEPTH)\n        image.append(m)\n    return image\n\ndef make_transfer_syntax_uid(df, dcm_dir):\n    machine_id_to_transfer = {}\n    machine_id = df.machine_id.unique()\n    for i in machine_id:\n        d = df[df.machine_id == i].iloc[0]\n        f = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dicom = pydicom.dcmread(f)\n        machine_id_to_transfer[i] = dicom.file_meta.TransferSyntaxUID\n    return machine_id_to_transfer\n\ndef normalised_to_8bit(image, photometric_interpretation):\n    xmin = image.min()\n    xmax = image.max()\n\n    norm = np.empty_like(image, dtype=np.uint8)\n    dicomsdl.util.convert_to_uint8(image, norm, xmin, xmax)\n    if photometric_interpretation == 'MONOCHROME1':\n        norm = 255 - norm\n    return norm\n\ndef resize_image_to_height(image, image_height):\n    \n    h, w = image.shape[:2]\n    s = image_height/h\n    if image_height!=h:\n        image = cv2.resize(image, dsize=None, fx=s, fy=s, interpolation=cv2.INTER_LINEAR)\n    return image\n\n#----------------------------------------------------------------\n# dicomsdl reader\ndef dicomsdl_to_numpy_image(ds, index=0):\n    # https://stackoverflow.com/questions/44659924/returning-numpy-arrays-via-pybind11\n    info = ds.getPixelDataInfo()\n    if info['SamplesPerPixel'] != 1:\n        raise RuntimeError('SamplesPerPixel != 1')\n\n    shape = [info['Rows'], info['Cols']]\n    dtype = info['dtype']\n    outarr = np.empty(shape, dtype=dtype)\n    ds.copyFrameData(index, outarr)\n    return outarr\n\ndef dicomsdl_parallel_process(d, dcm_dir, image_dir, image_height, is_voi_lut):\n    dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n    ds = dicomsdl.open(dcm_file)\n    image = dicomsdl_to_numpy_image(ds)\n    image = resize_image_to_height(image, image_height)\n\n    if is_voi_lut:\n        dc = pydicom.dcmread(dcm_file)\n        image = apply_voi_lut(image, dc)\n        image = image.astype(np.float32)\n    image = normalised_to_8bit(image, ds.PhotometricInterpretation)  # +1\n    image = _crop(image)\n    # save as png\n    os.makedirs(f'{image_dir}/{d.machine_id}/{d.patient_id}', exist_ok=True)\n    cv2.imwrite(f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png', image)\n\ndef process_non_j2k(df, dcm_dir, image_dir, image_height, n_jobs, is_voi_lut=True):\n    #https://stackoverflow.com/questions/56659294/does-joblib-parallel-keep-the-original-order-of-data-passed\n    #Parallel(n_jobs=2, backend='multiprocessing')(\n    Parallel(n_jobs=n_jobs)(\n        delayed(dicomsdl_parallel_process)(d, dcm_dir, image_dir, image_height, is_voi_lut)\n        for t,d in tqdm(df.iterrows())\n    )\n\n#----------------------------------------------------------------\n# nvjpeg2k reader\n\n'''\nTransferSyntaxUID\n1.2.840.10008.1.2.4.70 = JPEG Lossless, Nonhierarchical, First- Order Prediction (Processes 14)\n1.2.840.10008.1.2.4.90 = JPEG 2000 Image Compression (Lossless Only)\n'''\nj2k_decoder = nvjpeg2k.Decoder()\n\ndef process_j2k(df, dcm_dir, image_dir, image_height, is_voi_lut=True):\n    for t, d in tqdm(df.iterrows()):\n        dcm_file = f'{dcm_dir}/{d.patient_id}/{d.image_id}.dcm'\n        dc = pydicom.dcmread(dcm_file)\n        \n        offset = dc.PixelData.find(b'\\x00\\x00\\x00\\x0C')\n        jpeg_stream = bytearray(dc.PixelData[offset:])\n        image = j2k_decoder.decode(jpeg_stream)\n        image = resize_image_to_height(image, image_height)\n        #     jpeg_stream = bytearray(ds.PixelData[offset:])\n        #     jpeg_stream = np.array(bytearray(ds.PixelData[offset:]),np.uint8)\n\n        if is_voi_lut:\n            image = apply_voi_lut(image, dc)\n            image = image.astype(np.float32)\n        image = normalised_to_8bit(image, dc.PhotometricInterpretation)\n        image = _crop(image)\n        # save as png\n        os.makedirs(f'{image_dir}/{d.machine_id}/{d.patient_id}', exist_ok=True)\n        cv2.imwrite(f'{image_dir}/{d.machine_id}/{d.patient_id}/{d.image_id}.png', image)\n\n","metadata":{"execution":{"iopub.status.busy":"2023-02-09T19:02:16.466163Z","iopub.execute_input":"2023-02-09T19:02:16.466535Z","iopub.status.idle":"2023-02-09T19:02:16.521222Z","shell.execute_reply.started":"2023-02-09T19:02:16.466497Z","shell.execute_reply":"2023-02-09T19:02:16.519965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport cv2\nfrom timeit import default_timer as timer\nfrom tqdm.notebook import tqdm\nfrom joblib import Parallel, delayed\nfrom glob import glob\nfrom sklearn import metrics\nimport gc\n \n\nimport matplotlib\n#matplotlib.use('TkAgg')  \nimport matplotlib.pyplot as plt\n\nimport torch\nfrom torch.utils.data.dataset import Dataset\nfrom torch.utils.data import DataLoader\nfrom torch.utils.data.sampler import *\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.cuda.amp as amp\nprint( 'torch.cuda.device_count() = %d'%torch.cuda.device_count())\nprint( 'torch.cuda.get_device_properties() = %s' % str(torch.cuda.get_device_properties(0))[21:])\n\n\nimport nvjpeg2k\nprint('import ok!')","metadata":{"execution":{"iopub.status.busy":"2023-02-09T19:02:16.523718Z","iopub.execute_input":"2023-02-09T19:02:16.524109Z","iopub.status.idle":"2023-02-09T19:02:16.997286Z","shell.execute_reply.started":"2023-02-09T19:02:16.524072Z","shell.execute_reply":"2023-02-09T19:02:16.996298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"convert_height = 1536\n# image_height = 1536\n# image_width  = 960\n\n# csv_file = '/kaggle/input/rsna-breast-cancer-detection/test.csv'\n# dcm_dir  = '/kaggle/input/rsna-breast-cancer-detection/test_images'\ncsv_file = '/kaggle/input/rsna-breast-cancer-detection/train.csv'\ndcm_dir = '/kaggle/input/rsna-breast-cancer-detection/train_images'\n\ntest_df = pd.read_csv(csv_file)\nmachine_id_to_transfer = make_transfer_syntax_uid(test_df, dcm_dir)\ntest_df.loc[:, 'i'] = np.arange(len(test_df))\ntest_df.loc[:, 'TransferSyntaxUID'] = test_df.machine_id.map(machine_id_to_transfer)\nprint('test_df', test_df.shape)\nprint(test_df)\nprint('')","metadata":{"execution":{"iopub.status.busy":"2023-02-09T19:02:16.99998Z","iopub.execute_input":"2023-02-09T19:02:17.000354Z","iopub.status.idle":"2023-02-09T19:02:17.71689Z","shell.execute_reply.started":"2023-02-09T19:02:17.000318Z","shell.execute_reply":"2023-02-09T19:02:17.715806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# cache all png images\npng_dir = '/kaggle/png' \n\nj2k_df = test_df[test_df.TransferSyntaxUID == '1.2.840.10008.1.2.4.90'].reset_index(drop=True)\nnon_j2k_df = test_df[test_df.TransferSyntaxUID != '1.2.840.10008.1.2.4.90'].reset_index(drop=True)\n\n# print(f'process_j2k(): {len(j2k_df)}')\n# start_timer = timer()\n# process_j2k(j2k_df, dcm_dir, png_dir, convert_height)\n# #process_non_j2k(j2k_df, dcm_dir, png_dir, convert_height, n_jobs=2)\n# print(time_to_str(timer() - start_timer, 'sec'))\n\nprint(f'process_non_j2k(): {len(non_j2k_df)}')\nstart_timer = timer()\nprocess_non_j2k(non_j2k_df, dcm_dir, png_dir, convert_height, n_jobs=2)  \nprint(time_to_str(timer() - start_timer, 'sec'))\n","metadata":{"execution":{"iopub.status.busy":"2023-02-09T20:38:33.773497Z","iopub.execute_input":"2023-02-09T20:38:33.77473Z","iopub.status.idle":"2023-02-10T01:04:31.16859Z","shell.execute_reply.started":"2023-02-09T20:38:33.774602Z","shell.execute_reply":"2023-02-10T01:04:31.166252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import shutil\nshutil.make_archive('/kaggle/working/part2', 'zip', '/kaggle/png')","metadata":{"execution":{"iopub.status.busy":"2023-02-09T20:24:41.850713Z","iopub.execute_input":"2023-02-09T20:24:41.851116Z","iopub.status.idle":"2023-02-09T20:34:35.960538Z","shell.execute_reply.started":"2023-02-09T20:24:41.851082Z","shell.execute_reply":"2023-02-09T20:34:35.959232Z"},"trusted":true},"execution_count":null,"outputs":[]}]}