{"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":"# New crop and hist-scaled crop with Dali Tensor\n\nThis is the optimized version of **kena0ki,tivfrvqhs5 and christofhenkel**'s Notebook.\n\n### The changes I made are the followings.\n\n- **add hist-scaled methods for replacing 0-1 normalization**\n- **a new crop method for MLO and adjust the CC image corp**\n\n#### Reference:\n\n[SE-ResNeXt50 GPU optimized](https://www.kaggle.com/code/tivfrvqhs5/se-resnext50-gpu-optimized)\n\n[SE-ResNeXt50 full GPU decoding](https://www.kaggle.com/code/christofhenkel/se-resnext50-full-gpu-decoding)\n\n[[Optimized] RSNA ConvNextV2 Inference Tensorflow](https://www.kaggle.com/code/kena0ki/optimized-rsna-convnextv2-inference-tensorflow)","metadata":{}},{"cell_type":"code","source":"%%capture\n# Source: https://www.kaggle.com/code/remekkinas/fast-dicom-processing-1-6-2x-faster?scriptVersionId=113360473\n!pip install /kaggle/input/rsnamodules/dicomsdl-0.109.1-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl \n\ntry:\n    import pylibjpeg\nexcept:\n   !pip install /kaggle/input/rsna-2022-whl/{pylibjpeg-1.4.0-py3-none-any.whl,python_gdcm-3.0.15-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl}","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:53:54.796046Z","iopub.execute_input":"2023-03-07T03:53:54.796982Z","iopub.status.idle":"2023-03-07T03:54:57.649645Z","shell.execute_reply.started":"2023-03-07T03:53:54.796881Z","shell.execute_reply":"2023-03-07T03:54:57.648353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Install Keras CV Attention Model Pip Package for ConvNextV2 Models\n!pip install --no-deps /kaggle/input/keras-cv-attention-models/keras_cv_attention_models-1.3.9-py3-none-any.whl","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:54:57.652614Z","iopub.execute_input":"2023-03-07T03:54:57.653011Z","iopub.status.idle":"2023-03-07T03:55:19.989138Z","shell.execute_reply.started":"2023-03-07T03:54:57.65297Z","shell.execute_reply":"2023-03-07T03:55:19.988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install -q /kaggle/input/nvidia-dali-nightly-cuda110-1230dev/nvidia_dali_nightly_cuda110-1.23.0.dev20230203-7187866-py3-none-manylinux2014_x86_64.whl","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:55:19.990904Z","iopub.execute_input":"2023-03-07T03:55:19.992558Z","iopub.status.idle":"2023-03-07T03:56:02.032703Z","shell.execute_reply.started":"2023-03-07T03:55:19.992513Z","shell.execute_reply":"2023-03-07T03:56:02.031575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pylibjpeg\nimport pydicom\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport tensorflow as tf\n\nimport torch\nimport torch.nn.functional as F\n\nfrom joblib import Parallel, delayed\nfrom tqdm.notebook import tqdm\nfrom multiprocessing import cpu_count\nfrom keras_cv_attention_models import convnext\n\nimport cv2\nimport glob\nimport importlib\nimport os\nimport joblib\nimport time\nimport dicomsdl\nimport gc\n\n# Tensorflow and CV2 set number of threads to 1 for speedup in parallell function mapping\ntf.config.threading.set_inter_op_parallelism_threads(num_threads=1)\ncv2.setNumThreads(1)\n\n# Pandas DataFrame Display Options\npd.options.display.max_colwidth = 99","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:02.036341Z","iopub.execute_input":"2023-03-07T03:56:02.036758Z","iopub.status.idle":"2023-03-07T03:56:09.740379Z","shell.execute_reply.started":"2023-03-07T03:56:02.036722Z","shell.execute_reply":"2023-03-07T03:56:09.739436Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"IS_INTERACTIVE = os.environ['KAGGLE_KERNEL_RUN_TYPE'] == 'Interactive'\n\nTARGET_HEIGHT = 1344\nTARGET_WIDTH = 768\nN_CHANNELS = 1\nINPUT_SHAPE = (TARGET_HEIGHT, TARGET_WIDTH, N_CHANNELS)\nTARGET_HEIGHT_WIDTH_RATIO = TARGET_HEIGHT / TARGET_WIDTH\nTHRESHOLD_BEST = 0.55\n\nCLAHE = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(32, 32))\n\nCROP_IMAGE = True\nAPPLY_CLAHE = False\nAPPLY_EQ_HIST = False\nIMAGE_FORMAT = 'jpg'","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2023-03-07T03:56:09.741656Z","iopub.execute_input":"2023-03-07T03:56:09.742326Z","iopub.status.idle":"2023-03-07T03:56:09.755005Z","shell.execute_reply.started":"2023-03-07T03:56:09.742271Z","shell.execute_reply":"2023-03-07T03:56:09.753769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# VOI LUT","metadata":{}},{"cell_type":"code","source":"# Source: https://www.kaggle.com/code/bobdegraaf/dicomsdl-voi-lut\ndef voi_lut(image, dicom, util=np):\n    # Additional Checks\n    if 'WindowWidth' not in dicom.getPixelDataInfo() or 'WindowWidth' not in dicom.getPixelDataInfo():\n        return image\n    \n    # Load only the variables we need\n    center = dicom['WindowCenter']\n    width = dicom['WindowWidth']\n    bits_stored = dicom['BitsStored']\n    voi_lut_function = dicom['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        image = y_range / (1 + util.exp(-4 * (image - 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 = image <= (center - width / 2)\n        above = image > (center + width / 2)\n        between = util.logical_and(~below, ~above)\n\n        image[below] = y_min\n        image[above] = y_max\n        if between.any():\n            image[between] = (\n                ((image[between] - center) / width + 0.5) * y_range + y_min\n            )\n\n    return image","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:09.759Z","iopub.execute_input":"2023-03-07T03:56:09.759766Z","iopub.status.idle":"2023-03-07T03:56:09.780044Z","shell.execute_reply.started":"2023-03-07T03:56:09.75972Z","shell.execute_reply":"2023-03-07T03:56:09.778924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Crop Image","metadata":{}},{"cell_type":"code","source":"xoffset_times = 1.10\nn_threshold_y_top = 0.9  # devide this value to keep more bottom pixels\nn_threshold_y_bottom = 0.35  #  to find choose_last \n                             # (the last index of lines for regression of the cut line)\nthreshold_slope = -0.35 # threshold for control the slope the regressed cut line\nmove_factor = 0.15  # to control the moving amount of intercept when slope is smaller than threshold_slope\naccount_factor = 3  # to get choose_1st with the choose_last \nrotate_v = 1.3      # rotate the slope\nrotate_inter= 1\nglobal_max_check = 0.2   # control if we modify the value of choose_1st and choose_last \n# indexmax_begin = 100\n\n\n# Smooth vector used to smoothen sums/stds of axes\ndef smooth(l, util):\n    # kernel size is 1% of vector\n    kernel_size = int(len(l) * 0.01)\n    if hasattr(util, 'convolve'):\n        kernel = util.ones(kernel_size) / kernel_size\n        return util.convolve(l, kernel, mode='same')\n    else:\n        kernel = util.ones(kernel_size,device='cuda') / kernel_size\n        l = l.view(1, 1, l.size(0))\n        kernel = kernel.view(1, 1, kernel.size(0))\n        return F.conv1d(l.type(util.float), kernel, padding='same')[0][0]\n\n# X Crop offset based on first column with sum below 5% of maximum column sums*std\ndef get_x_offset(image, max_col_sum_ratio_threshold=0.05, debug=None, util=np,MLO=True):\n    # Image Dimensions\n    H, W = image.shape\n    # Percentual margin added to offset\n    margin = int(image.shape[1] * 0.00)\n    # Threshold values based on smoothed sum x std to capture varying intensity columns\n    vv = smooth(image.sum(axis=0).squeeze(), util) * smooth(image.std(axis=0).squeeze(), util)\n    # Find maximum sum in first 75% of columns\n    vv_argmax = vv[:int(image.shape[1] * 0.75)].argmax()\n    # Threshold value\n    vv_threshold = vv.max() * max_col_sum_ratio_threshold\n    \n    # Find first column after maximum column below threshold value\n    for offset, v in enumerate(vv):\n        # Start searching from vv_argmax\n        if offset < vv_argmax:\n            continue\n        \n        # Column below threshold value found\n        if v < vv_threshold:\n            offset = min(W, offset + margin)\n            break\n    offset = int(offset*xoffset_times)\n    \n    if isinstance(debug, np.ndarray):\n        if isinstance(image, torch.Tensor):\n            debug[1].imshow(image.cpu().numpy())\n            vv = vv.cpu().numpy()\n            vv_argmax = vv_argmax.cpu().numpy()\n            vv_threshold = vv_threshold.cpu().numpy()\n        else:\n            debug[1].imshow(image)\n        debug[1].set_title('X Offset')\n        vv_scale = H / vv.max() * 0.90\n        # Values\n        debug[1].plot(H - vv * vv_scale , c='red', label='vv')\n        # Threshold\n        debug[1].hlines(H - vv_threshold * vv_scale, 0, W -1, colors='orange', label='threshold')\n        # Max Value\n        debug[1].scatter(vv_argmax, H - vv[vv_argmax] * vv_scale, c='blue', s=100, label='Max', zorder=np.PINF)\n        # First Column Below Threshold\n        debug[1].scatter(offset, H - vv[offset] * vv_scale, c='purple', s=100, label='Offset', zorder=np.PINF)\n        debug[1].set_ylim(H, 0)\n        debug[1].legend()\n        debug[1].axis('off')\n        \n    return offset\n\n\n# Y Crop offset based on first bottom and top rows with sum below 10% of maximum row sum*std\ndef get_y_offsets_top2(image, max_row_sum_ratio_threshold=0.10, debug=None, util=np,MLO=True):\n    # Image Dimensions\n    H, W = image.shape\n    # Margin to add to offsets\n    margin = 0\n    # Threshold values based on smoothed sum x std to capture varying intensity columns\n    vv = smooth(image.sum(axis=1).squeeze(), util) * smooth(image.std(axis=1).squeeze(), util)\n    # Find maximum sum * std row in inter quartile rows\n    vv_argmax = int(image.shape[0] * 0.25) + vv[int(image.shape[0] * 0.25):int(image.shape[0] * 0.75)].argmax()\n    # Threshold value\n    vv_threshold = vv.max() * max_row_sum_ratio_threshold\n    # Default crop offsets\n    offset_bottom = 0\n\n    # Bottom offset, search from argmax to bottom\n    for offset in reversed(range(0, vv_argmax)):\n        v = vv[offset]\n        if v < vv_threshold:\n            offset_bottom = offset\n            break\n    return max(0, offset_bottom - margin)      \n\n# Y Crop offset based on first bottom and top rows with sum below 10% of maximum row sum*std\ndef get_y_offsets(image, max_row_sum_ratio_threshold=0.10, debug=None, util=np,MLO=True):\n    # Image Dimensions\n    H, W = image.shape\n    # Margin to add to offsets\n    margin = 0\n    # Threshold values based on smoothed sum x std to capture varying intensity columns\n    vv = smooth(image.sum(axis=1).squeeze(), util) * smooth(image.std(axis=1).squeeze(), util)\n    # Find maximum sum * std row in inter quartile rows\n    vv_argmax = int(image.shape[0] * 0.25) + vv[int(image.shape[0] * 0.25):int(image.shape[0] * 0.75)].argmax()\n    # Threshold value\n    vv_threshold = vv.max() * max_row_sum_ratio_threshold\n    # Default crop offsets\n    offset_bottom = 0\n    offset_top = H\n\n    # Bottom offset, search from argmax to bottom\n    for offset in reversed(range(0, vv_argmax)):\n        v = vv[offset]\n        if v < vv_threshold:\n            offset_bottom = offset\n            break\n    \n    vv_scale = W / vv.max() * 0.90\n    if isinstance(vv, torch.Tensor):\n        vs = (vv * vv_scale).cpu().numpy()\n    else:\n        vs =  (vv * vv_scale).astype(np.int16)\n            \n    \n    if MLO: \n        global_max = vs.argmax()  # list(vs).index(np.max(vs)) \n        choose_last = int(len(vs)*n_threshold_y_bottom)\n        \n        if global_max< (global_max_check*W):\n            choose_1st = global_max\n        else:                   \n            if choose_last>vv_argmax:\n                choose_last = vv_argmax\n            choose_1st = int(choose_last/account_factor)\n        choose_vs = vs[choose_1st:choose_last]\n        \n        try:\n            z1 = np.polyfit(range(len(choose_vs)),choose_vs,1)\n            # slope, intercept\n#             print(f'z1{z1}')\n            if z1[0]>threshold_slope:\n                z1[0] = threshold_slope\n                z1[1] = np.mean(choose_vs)+move_factor*W\n#                 print(f'modified z1{z1}')\n        except:\n            print('z1 = (0,0)')\n            z1 = (0,0)\n    \n    if isinstance(debug, np.ndarray):\n        if isinstance(image, torch.Tensor):\n            debug[2].imshow(image.cpu().numpy())\n            vv = vv.cpu().numpy()\n            vv_argmax = vv_argmax.cpu().numpy()\n            vv_threshold = vv_threshold.cpu().numpy()\n        else:\n            debug[2].imshow(image)\n        debug[2].set_title('Y Bottom Offset')\n        vv_scale = W / vv.max() * 0.90\n        # Values\n        debug[2].plot(vv * vv_scale, np.arange(H), c='red', label='vv')\n        # Threshold\n        debug[2].vlines(vv_threshold * vv_scale, 0, H -1, colors='orange', label='threshold')\n        # Max Value\n        debug[2].scatter(vv[vv_argmax] * vv_scale, vv_argmax, c='blue', s=100, label='Max', zorder=np.PINF)\n        # First Column Below Threshold\n        debug[2].scatter(vv[offset_bottom] * vv_scale, offset_bottom, c='purple', s=100, label='Offset', zorder=np.PINF)\n        debug[2].set_ylim(H, 0)\n        debug[2].legend()\n        debug[2].axis('off')\n            \n    # Top offset, search from argmax to top\n    for offset in range(vv_argmax, H):\n        v = vv[offset]\n        if v < vv_threshold:\n            offset_top = offset\n            break\n    \n    if MLO:\n        offset_top = int(offset_top/n_threshold_y_top)\n            \n    if isinstance(debug, np.ndarray):\n        if isinstance(image, torch.Tensor):\n            debug[3].imshow(image.cpu().numpy())\n#             vv = vv.cpu().numpy()\n#             vv_argmax = vv_argmax.cpu().numpy()\n#             vv_threshold = vv_threshold.cpu().numpy()\n        else:\n            debug[3].imshow(image)\n        if MLO:\n            debug[3].plot(np.arange(W)*z1[0]*rotate_v+z1[1]*rotate_inter, np.arange(W), c='blue', label='new')\n            \n        debug[3].set_title('Y Top Offset')\n        vv_scale = W / vv.max() * 0.90\n        # Values\n        debug[3].plot(vv * vv_scale, np.arange(H) , c='red', label='vv')\n        if MLO:\n            debug[3].plot(vs[choose_1st:choose_last],\n                          np.arange(H)[choose_1st:choose_last], c='yellow', label='vv check')\n        # Threshold\n        debug[3].vlines(vv_threshold * vv_scale, 0, H -1, colors='orange', label='threshold')\n        # Max Value\n        debug[3].scatter(vv[vv_argmax] * vv_scale, vv_argmax, c='blue', s=100, label='Max', zorder=np.PINF)\n        # First Column Below Threshold\n        debug[3].scatter(vv[offset_top] * vv_scale, offset_top, c='purple', s=100, label='Offset', zorder=np.PINF)\n        debug[2].set_ylim(H, 0)\n        debug[3].legend()\n        debug[3].axis('off')\n    if MLO:\n        return max(0, offset_bottom - margin), min(image.shape[0], offset_top + margin),z1\n    else:\n        return max(0, offset_bottom - margin), min(image.shape[0], offset_top + margin)\n            \n#     return max(0, offset_bottom - margin), min(image.shape[0], offset_top + margin)\n\n# Crop image and pad offsets to target image height/width ratio to preserve information\ndef crop(image, size=None, debug=False, util=np,MLO=True):\n    # Image dimensions\n    H, W = image.shape\n    # Compute x/bottom/top offsets\n#     x_offset = get_x_offset(image, debug=debug, util=util,MLO=True)\n#     offset_bottom, offset_top = get_y_offsets(image[:,:x_offset], debug=debug, util=util)\n    \n    x_offset = get_x_offset(image, debug=debug, util=util,MLO=MLO)\n    if MLO:\n        offset_bottom, offset_top,z1 = get_y_offsets(image[:,:x_offset], debug=debug, util=util,MLO = MLO)\n    else:\n        offset_bottom, offset_top = get_y_offsets(image[:,:x_offset], debug=debug, util=util,MLO=MLO)\n        \n    # Crop Height and Width\n    h_crop = offset_top - offset_bottom\n    w_crop = x_offset\n    \n    # Pad crop offsets to target aspect ratio\n    if size is not None:\n        # Height too large, pad x offset\n        if (h_crop / w_crop) > TARGET_HEIGHT_WIDTH_RATIO:\n            x_offset += int(h_crop / TARGET_HEIGHT_WIDTH_RATIO - w_crop)\n        else:\n            # Height too small, pad bottom/top offsets\n            offset_bottom -= int(0.50 * (w_crop * TARGET_HEIGHT_WIDTH_RATIO - h_crop))\n            offset_bottom_correction = max(0, -offset_bottom)\n            offset_bottom += offset_bottom_correction\n\n            offset_top += int(0.50 * (w_crop * TARGET_HEIGHT_WIDTH_RATIO - h_crop))\n            offset_top += offset_bottom_correction\n        \n        # Crop Image\n    if MLO and not (z1[0]==0 and z1[1]==0):\n        image = image[:offset_top,:x_offset]\n                    \n        pts = np.array([[0,0], [0,int(-z1[1]/z1[0])],[int(z1[1]),0],[0,0],])\n#         print(pts)\n        pts = np.array([pts])\n        mask = np.zeros(image.shape[:2], np.uint8)\n        cv2.polylines(mask, pts, 1, 255)    \n        cv2.fillPoly(mask, pts, 255)    \n        if isinstance(image, torch.Tensor):\n            image = image.cpu().numpy()\n            image = cv2.bitwise_and(image, image, mask=255-mask)\n            image = torch.from_numpy(image)\n        else:\n            image = cv2.bitwise_and(image, image, mask=255-mask)\n\n#         image = up\n        offset_bottom2 = get_y_offsets_top2(image, max_row_sum_ratio_threshold=0.10, debug=None,MLO=True) \n#         print(offset_bottom2)\n        image = image[offset_bottom2:,:]\n    else:\n        image = image[offset_bottom:offset_top:,:x_offset]\n        \n    return image","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:09.781603Z","iopub.execute_input":"2023-03-07T03:56:09.782115Z","iopub.status.idle":"2023-03-07T03:56:09.829478Z","shell.execute_reply.started":"2023-03-07T03:56:09.782081Z","shell.execute_reply":"2023-03-07T03:56:09.828478Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# DALI","metadata":{}},{"cell_type":"code","source":"from nvidia.dali import fn, math, pipeline_def, types\nfrom nvidia.dali.plugin.pytorch import feed_ndarray as feed_ndarray\n\n@pipeline_def\ndef jpeg_pipeline():\n    jpegs = fn.external_source(device=\"cpu\", name=\"jpeg\")\n    images = fn.experimental.decoders.image(jpegs, device='mixed', output_type=types.ANY_DATA, dtype=types.UINT16)\n    images = fn.cast(images, dtype=types.FLOAT)\n    return images\n\n\ndef read_encoded_stream(filename):\n    dcmfile = pydicom.dcmread(filename)   \n    if dcmfile.file_meta.TransferSyntaxUID == '1.2.840.10008.1.2.4.90':\n        offset = dcmfile.PixelData.find(b\"\\x00\\x00\\x00\\x0C\")   #<---- the jpeg2000 header\n    else:\n        offset = dcmfile.PixelData.find(b\"\\xff\\xd8\") #<---- the jpeg lossless header    \n\n    buff = np.array(bytearray(dcmfile.PixelData[offset:]), dtype=np.uint8)\n    return buff\n\ndef decode_jpg(file_path, dicom, pipe):\n    buff = read_encoded_stream(file_path)\n    pipe.feed_input(\"jpeg\", [buff])\n    out = pipe.run()\n    dali_img = out[0][0]\n\n    image = torch.empty(dali_img.shape(), dtype=torch.float, device=\"cuda\")\n    feed_ndarray(dali_img, image, cuda_stream=torch.cuda.current_stream(device=0))\n\n    return image.squeeze()","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:09.830765Z","iopub.execute_input":"2023-03-07T03:56:09.831232Z","iopub.status.idle":"2023-03-07T03:56:10.008237Z","shell.execute_reply.started":"2023-03-07T03:56:09.831194Z","shell.execute_reply":"2023-03-07T03:56:10.007319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Process Image","metadata":{}},{"cell_type":"code","source":"def process(file_path, size=(TARGET_WIDTH, TARGET_HEIGHT), crop_image=CROP_IMAGE, \n            apply_clahe=APPLY_CLAHE, apply_eq_hist=APPLY_EQ_HIST, \n            debug=False, save=True, pipe=None,MLO_type = \"MLO\"):\n    gpu = pipe is not None\n    # Read Dicom File\n    dicom = dicomsdl.open(file_path)\n    if gpu:\n        image = decode_jpg(file_path, dicom, pipe)\n    else:\n        image = dicom.pixelData()\n\n    if gpu:\n        util = torch\n    else:\n        util = np\n    \n    # Save original image for debug purposes\n    if debug:\n        fig, axes = plt.subplots(1, 5, figsize=(20,10))\n        if gpu:\n            image0 = image.cpu().numpy()\n        else:\n            image0 = np.copy(image)\n        axes[0].imshow(image0)\n        axes[0].set_title('Original Image')\n        axes[0].axis('off')\n    else:\n        axes = False\n    \n    # voi_lut\n    try:\n        image = voi_lut(image, dicom, util)\n    except:\n        pass\n    \n    # Some images have 0 values as highest intensity and need to be inverted\n    if dicom.getPixelDataInfo()['PhotometricInterpretation'] == 'MONOCHROME1':\n        image = util.max(image) - image\n\n    # Normalize [0,1] range\n    image = (image - image.min()) / (image.max() - image.min())\n\n    # Convert to uint8 image in range [0, 255]\n    if hasattr(image,'astype'):\n        image = (image * 255).astype(util.uint8)\n    else:\n        image = (image * 255).type(util.uint8).type(util.float)\n    \n    # Flip T0 Left/Right Orientation\n    h0, w0 = image.shape\n    if image[:,int(-w0 * 0.10):].sum() > image[:,:int(w0 * 0.10)].sum():\n        image = util.flip(image, (1,))\n    \n    # Crop Image\n    if crop_image:\n        MLO = False\n        if MLO_type=='MLO':\n            MLO = True\n        image = crop(image, debug=axes, util=util,MLO = MLO)\n        \n    # Resize\n    if size is not None:\n        # Pad black pixels to make square image\n        h, w = image.shape\n        if (h / w) > TARGET_HEIGHT_WIDTH_RATIO:\n            pad = int(h / TARGET_HEIGHT_WIDTH_RATIO - w)\n            if hasattr(util, 'pad'):\n                image = util.pad(image, [[0,0], [0, pad]])\n            else:\n                image = F.pad(image,(0,pad))\n            h, w = image.shape\n        else:\n            pad = int(0.50 * (w * TARGET_HEIGHT_WIDTH_RATIO - h))\n            if hasattr(util, 'pad'):\n                image = util.pad(image, [[pad, pad], [0,0]])\n            else:\n                image = F.pad(image,(pad,pad,0,0))\n            h, w = image.shape\n        # Resize\n        if gpu:\n            # The result of torch's interpolate and cv2 are slightly different.\n            # If you take this change serious,\n            # you may want to convert image to numpy ndarray before resizing.\n            image = F.interpolate(image.view(1, 1, h, w), (size[1],size[0]), mode=\"area\")[0, 0]\n        else:\n            image = cv2.resize(image, size, interpolation=cv2.INTER_AREA)\n\n    if gpu:\n        image = image.type(util.uint8).cpu().numpy()\n        \n    # Apply CLAHE contrast enhancement\n    if apply_clahe:\n        image = CLAHE.apply(image)\n        \n     # Apply Histogram Equalization\n    if apply_eq_hist:\n        image = cv2.equalizeHist(image)\n        \n    # Show Processed Image    \n    if debug:\n        axes[4].imshow(image)\n        axes[4].set_title('Processed Image')\n        axes[4].axis('off')\n        plt.show()\n        \n    # Save Only\n    if save:\n        image_id = file_path.split('/')[-1].split('.')[0]\n        if IMAGE_FORMAT == 'png':\n            cv2.imwrite(f'{image_id}.png', image)\n        else:\n            cv2.imwrite(f'{image_id}.jpg', image, [cv2.IMWRITE_JPEG_QUALITY, 95])","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:10.009556Z","iopub.execute_input":"2023-03-07T03:56:10.009891Z","iopub.status.idle":"2023-03-07T03:56:10.029752Z","shell.execute_reply.started":"2023-03-07T03:56:10.009855Z","shell.execute_reply":"2023-03-07T03:56:10.028514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def process2(file_path, size=(TARGET_WIDTH, TARGET_HEIGHT), crop_image=CROP_IMAGE, \n            apply_clahe=APPLY_CLAHE, apply_eq_hist=APPLY_EQ_HIST, \n            debug=False, save=True, pipe=None,MLO_type = \"MLO\"):\n    gpu = pipe is not None\n    # Read Dicom File\n    dicom = dicomsdl.open(file_path)\n    if gpu:\n        image = decode_jpg(file_path, dicom, pipe)\n    else:\n        image = dicom.pixelData()\n\n    if gpu:\n        util = torch\n    else:\n        util = np\n    \n    # Save original image for debug purposes\n    if debug:\n        fig, axes = plt.subplots(1, 5, figsize=(20,10))\n        if gpu:\n            image0 = image.cpu().numpy()\n        else:\n            image0 = np.copy(image)\n        axes[0].imshow(image0)\n        axes[0].set_title('Original Image')\n        axes[0].axis('off')\n    else:\n        axes = False\n    \n    # voi_lut\n    try:\n        image = voi_lut(image, dicom, util)\n    except:\n        pass\n    \n    # Some images have 0 values as highest intensity and need to be inverted\n    if dicom.getPixelDataInfo()['PhotometricInterpretation'] == 'MONOCHROME1':\n        image = util.max(image) - image\n\n    # Normalize [0,1] range\n    image = (image - image.min()) / (image.max() - image.min())\n\n    # Convert to uint8 image in range [0, 255]\n    if hasattr(image,'astype'):\n        image = (image * 255).astype(util.uint8)\n    else:\n        image = (image * 255).type(util.uint8).type(util.float)\n    \n    # Flip T0 Left/Right Orientation\n    h0, w0 = image.shape\n    if image[:,int(-w0 * 0.10):].sum() > image[:,:int(w0 * 0.10)].sum():\n        image = util.flip(image, (1,))\n    \n    # Crop Image\n    if crop_image:\n        MLO = False\n        if MLO_type=='MLO':\n            MLO = True\n        image = crop(image, debug=axes, util=util,MLO = MLO)\n        \n    # Resize\n    if size is not None:\n        # Pad black pixels to make square image\n        h, w = image.shape\n        if (h / w) > TARGET_HEIGHT_WIDTH_RATIO:\n            pad = int(h / TARGET_HEIGHT_WIDTH_RATIO - w)\n            if hasattr(util, 'pad'):\n                image = util.pad(image, [[0,0], [0, pad]])\n            else:\n                image = F.pad(image,(0,pad))\n            h, w = image.shape\n        else:\n            pad = int(0.50 * (w * TARGET_HEIGHT_WIDTH_RATIO - h))\n            if hasattr(util, 'pad'):\n                image = util.pad(image, [[pad, pad], [0,0]])\n            else:\n                image = F.pad(image,(pad,pad,0,0))\n            h, w = image.shape\n        # Resize\n        if gpu:\n            # The result of torch's interpolate and cv2 are slightly different.\n            # If you take this change serious,\n            # you may want to convert image to numpy ndarray before resizing.\n            image = F.interpolate(image.view(1, 1, h, w), (size[1],size[0]), mode=\"area\")[0, 0]\n        else:\n            image = cv2.resize(image, size, interpolation=cv2.INTER_AREA)\n\n    if gpu:\n        image = image.type(util.uint8).cpu().numpy()\n        \n    # Apply CLAHE contrast enhancement\n    if apply_clahe:\n        image = CLAHE.apply(image)\n        \n     # Apply Histogram Equalization\n    if apply_eq_hist:\n        image = cv2.equalizeHist(image)\n        \n    # Show Processed Image    \n    if debug:\n        axes[4].imshow(image)\n        axes[4].set_title('Processed Image')\n        axes[4].axis('off')\n        plt.show()\n        \n    # Save Only\n    if save:\n        image_id = file_path.split('/')[-1].split('.')[0]\n        if IMAGE_FORMAT == 'png':\n            cv2.imwrite(f'{image_id}.png', image)\n        else:\n            cv2.imwrite(f'{image_id}.jpg', image, [cv2.IMWRITE_JPEG_QUALITY, 95])\n    return image","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:10.031433Z","iopub.execute_input":"2023-03-07T03:56:10.03186Z","iopub.status.idle":"2023-03-07T03:56:10.052118Z","shell.execute_reply.started":"2023-03-07T03:56:10.031824Z","shell.execute_reply":"2023-03-07T03:56:10.050979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Example Preprocessing","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/train.csv')\n    \ndef get_file_path(args):\n    patient_id, image_id = args\n    return f'/kaggle/input/rsna-breast-cancer-detection/train_images/{patient_id}/{image_id}.dcm'\n    \ntrain['file_path'] = train[['patient_id', 'image_id']].apply(get_file_path, axis=1)","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:10.056822Z","iopub.execute_input":"2023-03-07T03:56:10.057107Z","iopub.status.idle":"2023-03-07T03:56:10.588663Z","shell.execute_reply.started":"2023-03-07T03:56:10.057063Z","shell.execute_reply":"2023-03-07T03:56:10.587663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"p = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\np.build()\n","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:10.590193Z","iopub.execute_input":"2023-03-07T03:56:10.59058Z","iopub.status.idle":"2023-03-07T03:56:10.80109Z","shell.execute_reply.started":"2023-03-07T03:56:10.590542Z","shell.execute_reply":"2023-03-07T03:56:10.800112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"path = '/kaggle/input/rsna-breast-cancer-detection/train_images/32292/924198868.dcm'\n\ntry:\n#         process(fp[-1],pipe=p,MLO_type=fp[4])\n    image1 = process2(path, crop_image=True, size=(TARGET_WIDTH, TARGET_HEIGHT), \n        debug=True, save=False,pipe=p,MLO_type='CC')\nexcept Exception as e:\n    print('Fall back to CPU')\n#         process(fp[-1],MLO_type=fp[4])\n    image2 = process2(path, crop_image=True, size=(TARGET_WIDTH, TARGET_HEIGHT), \n        debug=True, save=False,MLO_type='CC')\n    p = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\n    p.build()","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:10.802706Z","iopub.execute_input":"2023-03-07T03:56:10.80308Z","iopub.status.idle":"2023-03-07T03:56:24.574479Z","shell.execute_reply.started":"2023-03-07T03:56:10.803043Z","shell.execute_reply":"2023-03-07T03:56:24.57354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for histscale norm\nfrom fastai.vision.all import tensor\nimport torch\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\n\ndef array_freqhist_bins(image, n_bins=100):\n    \"A numpy based function to split the range of pixel values into groups, such that each group has around the same number of pixels\"\n    imsd = np.sort(image.flatten())\n    t = np.array([0.001])\n    t = np.append(t, np.arange(n_bins)/n_bins+(1/2/n_bins))\n    t = np.append(t, 0.999)\n    t = (len(imsd)*t+0.5).astype(int)\n    return np.unique(imsd[t])\n\ndef freqhist_bins(image, n_bins=100):\n    \"A function to split the range of pixel values into groups, such that each group has around the same number of pixels\"\n    imsd = image.view(-1).sort()[0]\n    t = torch.cat([tensor([0.001]),\n                   torch.arange(n_bins).float()/n_bins+(1/2/n_bins),\n                   tensor([0.999])])\n    t = (len(imsd)*t).long()\n    return imsd[t].unique()","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:24.575558Z","iopub.execute_input":"2023-03-07T03:56:24.57591Z","iopub.status.idle":"2023-03-07T03:56:25.557159Z","shell.execute_reply.started":"2023-03-07T03:56:24.575875Z","shell.execute_reply":"2023-03-07T03:56:25.556093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def hist_scaled_pt(image, brks=None):\n    # Pytorch-only version - switch to this if/when interp_1d can be optimized\n    if brks is None: brks = freqhist_bins(image)\n    brks = brks.to(device)\n    ys = torch.linspace(0., 1., len(brks)).to(device)\n    return image.flatten().interp_1d(brks, ys).reshape(image.shape).clamp(0.,1.)\n\ndef hist_scaled(image, brks=None):\n#     \"Scales a tensor using `freqhist_bins` to values between 0 and 1\"\n    if device.type=='cuda': return hist_scaled_pt(image,brks)\n    if brks is None: brks = freqhist_bins(image)\n    ys = np.linspace(0., 1., len(brks))\n    x = self.numpy().flatten()\n    x = np.interp(x, brks.numpy(), ys)\n    return tensor(x).reshape(self.shape).clamp(0.,1.)","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:25.558707Z","iopub.execute_input":"2023-03-07T03:56:25.559358Z","iopub.status.idle":"2023-03-07T03:56:25.568785Z","shell.execute_reply.started":"2023-03-07T03:56:25.55932Z","shell.execute_reply":"2023-03-07T03:56:25.567562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test","metadata":{}},{"cell_type":"code","source":"file_path = '/kaggle/input/rsna-breast-cancer-detection/train_images/32292/924198868.dcm'\ndicom = dicomsdl.open(file_path)\nimage = decode_jpg(file_path, dicom, p)","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:25.570781Z","iopub.execute_input":"2023-03-07T03:56:25.571672Z","iopub.status.idle":"2023-03-07T03:56:25.698452Z","shell.execute_reply.started":"2023-03-07T03:56:25.571629Z","shell.execute_reply":"2023-03-07T03:56:25.697319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"t_bin = freqhist_bins(image,n_bins=100)\n\nt_bin2 = t_bin.cpu().numpy()\nplt.hist(t_bin2, bins=t_bin2, color='c');\nplt.show()\nplt.plot(t_bin2, np.linspace(0,1,len(t_bin2)));\nplt.show()\na = plt.hist(image.flatten().cpu().numpy(), bins=100)","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:25.700214Z","iopub.execute_input":"2023-03-07T03:56:25.700945Z","iopub.status.idle":"2023-03-07T03:56:26.70492Z","shell.execute_reply.started":"2023-03-07T03:56:25.700904Z","shell.execute_reply":"2023-03-07T03:56:26.703985Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test for single image","metadata":{}},{"cell_type":"code","source":"# if not register the p, uncomment these two lines\n# p = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\n# p.build()\n\npath = '/kaggle/input/rsna-breast-cancer-detection/train_images/32292/924198868.dcm'\n\ndicom = dicomsdl.open(path)\nimage = decode_jpg(path, dicom, p)\n# image.type   # <function Tensor.type>\n\ndata_scaled = hist_scaled(image)\nplt.hist(image.flatten().cpu().numpy());\nplt.show()\nplt.hist(data_scaled.flatten().cpu().numpy());\nplt.show()\n\nplt.figure(figsize=(6, 4),dpi=100)\nplt.imshow(image.cpu().numpy(), cmap=plt.cm.bone)\nplt.title(\"Original\")\nplt.axis('off')\nplt.show()\n\nplt.figure(figsize=(6, 4),dpi=100)\nplt.imshow(data_scaled.cpu().numpy(), cmap=plt.cm.bone)\nplt.title(\"Processed\")\nplt.axis('off')\nplt.show()\n\nprint(image.min(),image.max(),image.mean(),image.std(),'\\n',\ndata_scaled.min(),data_scaled.max(),data_scaled.mean(),data_scaled.std())\n","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:26.706458Z","iopub.execute_input":"2023-03-07T03:56:26.706813Z","iopub.status.idle":"2023-03-07T03:56:29.257663Z","shell.execute_reply.started":"2023-03-07T03:56:26.706778Z","shell.execute_reply":"2023-03-07T03:56:29.256617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# hist-scaled method comparision","metadata":{}},{"cell_type":"code","source":"paths = ['/kaggle/input/rsna-breast-cancer-detection/train_images/25639/12367200.dcm',\n    '/kaggle/input/rsna-breast-cancer-detection/train_images/32292/924198868.dcm',\n         '/kaggle/input/rsna-breast-cancer-detection/train_images/3765/193353300.dcm',\n         '/kaggle/input/rsna-breast-cancer-detection/train_images/10042/102733848.dcm',\n         '/kaggle/input/rsna-breast-cancer-detection/train_images/41660/144622564.dcm',\n         \n         ]\n\n# p = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\n# p.build()\n\n\nfig, axs = plt.subplots(2,len(paths),figsize=(10, 4))\nfig.set_dpi(300)\nimages = []\nfor ii,i in enumerate(paths):\n    dicom = dicomsdl.open(i)\n    image = decode_jpg(i, dicom, p)\n    print(image.mean(),image.std(),image.max(),image.min())\n    data_scaled = hist_scaled(image)\n    print(data_scaled.mean(),data_scaled.std(),\n          data_scaled.max(),data_scaled.min())\n\n    axs[0,ii].imshow(image.cpu().numpy())\n    axs[0,ii].axis('off')\n    axs[1,ii].imshow(data_scaled.cpu().numpy())\n    axs[1,ii].axis('off')\nfig.suptitle('fastai hist_scaled normalization \\n upper: original,below: hist_scaled')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:29.258892Z","iopub.execute_input":"2023-03-07T03:56:29.259244Z","iopub.status.idle":"2023-03-07T03:56:38.717127Z","shell.execute_reply.started":"2023-03-07T03:56:29.259208Z","shell.execute_reply":"2023-03-07T03:56:38.71629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The background is 0, and the first (or two) on the left is a common normal image. \n\nThe other four pairs on the right need to be processed","metadata":{}},{"cell_type":"markdown","source":"## testing of the new crop method ","metadata":{}},{"cell_type":"code","source":"# GPU\n\nN = 8\n\np = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\np.build()\n\nfor fp in tqdm((train.head(N)).values):  # for fp in tqdm(train['file_path'].head(N)):\n    process(fp[-1], crop_image=True, size=(TARGET_WIDTH, TARGET_HEIGHT), \n            debug=True, save=False, pipe=p,MLO_type=fp[4])","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:56:38.718412Z","iopub.execute_input":"2023-03-07T03:56:38.719513Z","iopub.status.idle":"2023-03-07T03:57:03.962161Z","shell.execute_reply.started":"2023-03-07T03:56:38.719473Z","shell.execute_reply":"2023-03-07T03:57:03.96114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# GPU Image Processing","metadata":{}},{"cell_type":"code","source":"test = pd.read_csv('/kaggle/input/rsna-breast-cancer-detection/test.csv')\n\ndef get_file_path(args):\n    patient_id, image_id = args\n    return f'/kaggle/input/rsna-breast-cancer-detection/test_images/{patient_id}/{image_id}.dcm'\n    \ntest['file_path'] = test[['patient_id', 'image_id']].apply(get_file_path, axis=1)\n\ndisplay(test.info())\ndisplay(test.head())","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:57:03.963669Z","iopub.execute_input":"2023-03-07T03:57:03.964111Z","iopub.status.idle":"2023-03-07T03:57:04.007014Z","shell.execute_reply.started":"2023-03-07T03:57:03.964075Z","shell.execute_reply":"2023-03-07T03:57:04.005908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\np = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\np.build()\n\nfor fp in tqdm(test.values):\n    try:\n        process(fp[-1],pipe=p,MLO_type=fp[4])\n    except Exception as e:\n        print('Fall back to CPU')\n        process(fp[-1],MLO_type=fp[4])\n        p = jpeg_pipeline(batch_size=1, num_threads=2, device_id=0, prefetch_queue_depth=1)\n        p.build()","metadata":{"execution":{"iopub.status.busy":"2023-03-07T03:57:04.009411Z","iopub.execute_input":"2023-03-07T03:57:04.010016Z","iopub.status.idle":"2023-03-07T03:57:05.142335Z","shell.execute_reply.started":"2023-03-07T03:57:04.009977Z","shell.execute_reply":"2023-03-07T03:57:05.141328Z"},"trusted":true},"execution_count":null,"outputs":[]}]}