{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":45867,"databundleVersionId":6688004,"sourceType":"competition"}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Cropping WSIs images\n**Whole Slide Images (WSIs)** are high-resolution images used primarily in the field of pathology. They represent an entire pathology slide and are used for digital pathology reviews, research, and more. Given their high resolution, WSIs can be extremely large, sometimes even exceeding gigabytes in size.\n\n\"Duplicate regions\" in WSIs refer to areas on the slide where the same tissue or sample appears more than once. There are several reasons duplicate regions might exist:\n\n1. **Tissue Folding**: During the tissue preparation and embedding process, the tissue might fold onto itself. When the slide is scanned, this fold can appear as duplicate or even triplicate regions of the same tissue.\n\n2. **Multiple Sections**: Sometimes, a single tissue block is cut into multiple sections and placed on the same slide for comparative analysis. This can create areas that look nearly identical.\n\n3. **Slide Preparation Artifacts**: Mistakes or issues during the slide preparation can lead to unintentional replication of regions.\n\n4. **Scanning Artifacts**: There can be errors during the digital scanning of the slides, leading to repeated regions in the digital image.\n\n5. **Intentional Duplications**: In some cases, particularly in teaching slides or reference materials, duplicate regions might be intentionally created for comparison or demonstration purposes.\n\nIt's essential to be aware of and account for these duplicate regions, especially in computational pathology. When analyzing WSIs for quantitative features or training machine learning models, duplicate regions can introduce bias or skew results.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import os\nimport shutil\n\nimport numpy as np\nimport pandas as pd\nimport torch\nimport matplotlib.pyplot as plt\n\nfrom datasets import load_dataset\n","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:11.480794Z","iopub.execute_input":"2023-10-26T14:29:11.481319Z","iopub.status.idle":"2023-10-26T14:29:17.050304Z","shell.execute_reply.started":"2023-10-26T14:29:11.481275Z","shell.execute_reply":"2023-10-26T14:29:17.048933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(\"/kaggle/input/UBC-OCEAN/train.csv\")\ntest_df = pd.read_csv(\"/kaggle/input/UBC-OCEAN/test.csv\")\n\nBASE_DIR = [\"/kaggle/input/UBC-OCEAN/train_thumbnails/\", \"/kaggle/input/UBC-OCEAN/test_thumbnails/\"]","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:17.052988Z","iopub.execute_input":"2023-10-26T14:29:17.054029Z","iopub.status.idle":"2023-10-26T14:29:17.088644Z","shell.execute_reply.started":"2023-10-26T14:29:17.053973Z","shell.execute_reply":"2023-10-26T14:29:17.087268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def organize_images_by_label(df: pd.DataFrame, source_dir: str) -> None:\n    image_paths=[]\n    for _, row in df.iterrows():\n        image_id = row[\"image_id\"]\n        label = row[\"label\"]\n        source_path = os.path.join(source_dir, f\"{image_id}_thumbnail.png\") \n        try:\n            image_paths.append(os.path.join(source_dir, f\"{image_id}_thumbnail.png\"))\n        except FileNotFoundError:\n            image_paths.append(1)\n            continue\n    return image_paths\n\n\nimage_paths = organize_images_by_label(train_df, BASE_DIR[0])\ntrain_df['image_path'] = image_paths\n","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:17.090437Z","iopub.execute_input":"2023-10-26T14:29:17.090913Z","iopub.status.idle":"2023-10-26T14:29:17.158106Z","shell.execute_reply.started":"2023-10-26T14:29:17.090844Z","shell.execute_reply":"2023-10-26T14:29:17.156652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.signal import find_peaks, savgol_filter\nfrom PIL import Image\nimport cv2\nimport numpy as np\n\n# Randomly choose 5 images from each group\nrandom_5_images = train_df[train_df['image_path']!=1].groupby('label').head(5).sort_values('label')","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:17.164444Z","iopub.execute_input":"2023-10-26T14:29:17.16496Z","iopub.status.idle":"2023-10-26T14:29:17.998512Z","shell.execute_reply.started":"2023-10-26T14:29:17.164914Z","shell.execute_reply":"2023-10-26T14:29:17.996974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Key Functions\n### 1. Signal Smoothing\nThe smooth_signal function smooths a given signal using the Savitzky-Golay filter.","metadata":{}},{"cell_type":"code","source":"def smooth_signal(signal, window_length=5, polyorder=3):\n    return savgol_filter(signal, window_length, polyorder)","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:18.000568Z","iopub.execute_input":"2023-10-26T14:29:18.001119Z","iopub.status.idle":"2023-10-26T14:29:18.010556Z","shell.execute_reply.started":"2023-10-26T14:29:18.001052Z","shell.execute_reply":"2023-10-26T14:29:18.007854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. Peak Detection\ndetect_peaks identifies peaks in the smoothed signal based on a prominence threshold.","metadata":{}},{"cell_type":"code","source":"def detect_peaks(signal, prominence=150):  # Adjust prominence as needed\n    peaks, _ = find_peaks(signal, prominence=prominence)\n    return peaks","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:18.013006Z","iopub.execute_input":"2023-10-26T14:29:18.014052Z","iopub.status.idle":"2023-10-26T14:29:18.022192Z","shell.execute_reply.started":"2023-10-26T14:29:18.01399Z","shell.execute_reply":"2023-10-26T14:29:18.020592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Image Cropping\nThe crop_image function processes the image to detect regions of interest (based on contours) and crops these regions out.\n\n","metadata":{}},{"cell_type":"code","source":"def crop_image(input_image):\n    image = cv2.imread(input_image)\n\n    # Convert to grayscale\n    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)\n\n    # Threshold the image to create a binary image\n    _, binary = cv2.threshold(gray, 1, 255, cv2.THRESH_BINARY)\n\n    # Find the contours (regions) in the binary image\n    contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)\n\n    # Extract bounding boxes around each contour and save them as individual images\n    cropped_images=[]\n    for i, contour in enumerate(contours):\n        x, y, w, h = cv2.boundingRect(contour)\n        cropped_image = image[y:y+h, x:x+w]\n        if np.sum(np.array(cropped_image))>100000:\n            cropped_images.append(cropped_image)\n    return cropped_images","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:18.024647Z","iopub.execute_input":"2023-10-26T14:29:18.025255Z","iopub.status.idle":"2023-10-26T14:29:18.03848Z","shell.execute_reply.started":"2023-10-26T14:29:18.025199Z","shell.execute_reply":"2023-10-26T14:29:18.037275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Signal Scaling\nscale_signal scales a signal to a specified range. This can be useful for visualization purposes.","metadata":{}},{"cell_type":"code","source":"def scale_signal(signal, min_range, max_range):\n    # Find the minimum and maximum of the signal\n    min_value = np.min(signal)\n    max_value = np.max(signal)\n    \n    # Scale the signal\n    scaled_signal = min_range + (signal - min_value) / (max_value - min_value) * (max_range - min_range)\n    \n    return scaled_signal","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:18.039905Z","iopub.execute_input":"2023-10-26T14:29:18.041166Z","iopub.status.idle":"2023-10-26T14:29:18.055665Z","shell.execute_reply.started":"2023-10-26T14:29:18.041117Z","shell.execute_reply":"2023-10-26T14:29:18.05433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Workflow\n## A dataset (assumed to be train_df) is filtered to select random images.\nFor each image:\n* The image is opened and the signal is extracted and smoothed.\n* Peaks in the signal are detected.\n* If multiple peaks are found, the image is assumed to have multiple regions of interest. These regions are then cropped out and visualized alongside the original image and the corresponding signal.\n","metadata":{}},{"cell_type":"code","source":"for i, row in random_5_images.iterrows():\n    if os.path.isfile(row['image_path']):\n        image = Image.open(row['image_path'])\n        smoothed_signal = smooth_signal(np.mean(np.array(image)[:,:,0], axis=0))\n        peaks = detect_peaks(smoothed_signal)\n        if len(peaks)>1: #Check if there is more than one image\n            image_path = row['image_path']  \n            images = crop_image(image_path)\n            fig, axarr = plt.subplots(1, 2+len(images), figsize=(3*(2+len(images)),3))\n            fig.suptitle(\"Image ID: {}, Label: {}\".format(row['image_id'], row['label']))\n            for i, ax in enumerate(axarr.ravel()):\n                if i==0:\n                    ax.imshow(np.array(image))\n                    ax.axis('off')\n                    ax.set_title(\"Original Image\")\n                elif i==1:\n                    _range=axarr[0].get_ylim()\n                    ax.plot(scale_signal(smoothed_signal, _range[1],  _range[0]))\n                    ax.axis('off')\n                    ax.set_title(\"1D Signal of original Image\")\n                else:\n                    ax.imshow(np.array(images)[i-2])\n                    ax.axis('off')\n                    ax.set_title(\"Cropped Image {}\".format(i-2))\n            plt.tight_layout()\n            plt.show()\n            \n#         else:\n#             # Use this if you want to visualize non-multiple images\n#             image_path = row['image_path']  \n#             images = crop_image(image_path)\n#             fig, axarr = plt.subplots(1, 2, figsize=(10*(2),10))\n        \n#             for i, ax in enumerate(axarr.ravel()):\n#                 if i==0:\n#                     ax.imshow(np.array(image))\n#                 elif i==1:\n#                     ax.plot(smoothed_signal)\n#             plt.tight_layout()\n#             plt.show()\n#         ax.imshow(image)\n#         ax.set_title(len(peaks))\n#         ax.axis('off')  # To hide axis labels\n\n","metadata":{"execution":{"iopub.status.busy":"2023-10-26T14:29:18.057402Z","iopub.execute_input":"2023-10-26T14:29:18.057857Z","iopub.status.idle":"2023-10-26T14:29:39.600756Z","shell.execute_reply.started":"2023-10-26T14:29:18.057819Z","shell.execute_reply":"2023-10-26T14:29:39.598828Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# End of the Notebook","metadata":{}}]}