{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":22990,"databundleVersionId":2048213,"sourceType":"competition"},{"sourceId":1782290,"sourceType":"datasetVersion","datasetId":1059618},{"sourceId":2010866,"sourceType":"datasetVersion","datasetId":979056},{"sourceId":2072010,"sourceType":"datasetVersion","datasetId":1241983},{"sourceId":2180975,"sourceType":"datasetVersion","datasetId":1309339},{"sourceId":2210499,"sourceType":"datasetVersion","datasetId":1327520},{"sourceId":2219451,"sourceType":"datasetVersion","datasetId":1332847}],"dockerImageVersionId":30096,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# HuBMAP - Hacking the Kidney\nLet’s understand the competition first.\nThis competition is about developing a segmentation problem to identify regions with glomeruli in human kidney tissue images. \n* train : 15 images\n* test : 5 images\n* Each of them has ~50k pixel size and is saved as a high-resolution tiff image.\n* Credits:
\n[https://www.kaggle.com/yassinealouini/all-the-segmentation-metrics](http://)
\n[https://www.kaggle.com/iafoss/256x256-images](http://)
\n[https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter](http://)
\n[https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub](http://)
\n[https://www.kaggle.com/harshsharma511/one-stop-understanding-eda-efficientunet](http://)
\n[https://www.kaggle.com/vineeth1999/hubmap-eda-pytorch-efficientunet-offline-training](http://)
\n","metadata":{}},{"cell_type":"code","source":"import cv2\nimport datetime\nimport gc\nimport glob\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport skimage.morphology\nimport sys\nimport tifffile\n\n\nplot_full_image = True\n\n# Number of glomeruli to display for each image\nnum_glom_display = 5\n\n# Number of glomberuli to save as tiff files.\nnum_glom_save = 5\n\nglob_scale = 0.25\nbase_path = '../input/hubmap-kidney-segmentation'\n#Directory Contents\nprint(\"Directory Contents\")\nprint('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' )\nprint('\\n'.join(os.listdir(base_path)))\nprint('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' )\ntrain_files = sorted(glob.glob(os.path.join(base_path, 'train/*.tiff')))\nprint(f'Number of training images: {len(train_files)}')\nprint('\\n'.join(train_files))\nprint('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' )\n\n#Test Images\ntest_files = sorted(glob.glob(os.path.join(base_path, 'test/*.tiff')))\nprint(f'Number of test images: {len(test_files)}')\nprint('\\n'.join(test_files))\n\ndf_train = pd.read_csv(os.path.join(base_path, 'train.csv'))\ndf_info = pd.read_csv(os.path.join(base_path,'HuBMAP-20-dataset_information.csv'))\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:35:25.189929Z","iopub.execute_input":"2024-01-15T13:35:25.190287Z","iopub.status.idle":"2024-01-15T13:35:27.134073Z","shell.execute_reply.started":"2024-01-15T13:35:25.190207Z","shell.execute_reply":"2024-01-15T13:35:27.133286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing the trainging images and thier masks ","metadata":{}},{"cell_type":"markdown","source":"The masks for the training images is in the csv file (train.csv) .It contains RLE encoded (run-length encoding) of the mask for each training image.\nlet's talk about RLE: RLE is lossless compression algorithim that helps to compress when there is a lot of repetation . In this competation the mask is compressed such that the segment is encoded as pairs of values that contain a start position and a run length for example 0 5 10 3 implies pixels 0,1,2,3,4 and 10,11,12 .The pixels are numbered from top to bottom, then left to right.\n","metadata":{}},{"cell_type":"code","source":"#credit https://www.kaggle.com/harshsharma511/one-stop-understanding-eda-efficientunet\ndef rle_to_image(rle_mask, image_shape):\n \"\"\"\n Converts an rle string to an image represented as a numpy array.\n Reference: https://www.kaggle.com/paulorzp/rle-functions-run-lenght-encode-decode\n\n :param rle_mask: string with rle mask.\n :param image_shape: (width, height) of array to return\n :return: Image as a numpy array. 1 = mask, 0 = background.\n \"\"\"\n\n # Processing\n s = rle_mask.split()\n starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n starts -= 1\n ends = starts + lengths\n image = np.zeros(image_shape[0] * image_shape[1], dtype=np.uint8)\n for lo, hi in zip(starts, ends):\n image[lo:hi] = 1\n\n return image.reshape(image_shape).T\n\n\ndef overlay_image_mask(image, mask, mask_color=(0,255,0), alpha=1.0):\n im_f= image.astype(np.float32)\n# if mask.ndim == 2:\n# mask = np.expand_dims(mask,-1) \n mask_col = np.expand_dims(np.array(mask_color)/255.0, axis=(0,1))\n return (im_f + alpha * mask * (np.mean(0.8 * im_f + 0.2 * 255, axis=2, keepdims=True) * mask_col - im_f)).astype(np.uint8)\n\n\ndef overlay_image_mask_original(image, mask, mask_color=(0,255,0), alpha=1.0):\n return np.concatenate((image, overlay_image_mask(image, mask)), axis=1)\n\ndef get_image_id(image_file):\n return os.path.splitext(os.path.split(image_file)[1])[0]\n\n\ndef read_image(image_file, scale=1.0):\n image = tifffile.imread(image_file).squeeze()\n if image.shape[0] == 3:\n image = np.transpose(image, (1,2,0))\n \n orig_shape = image.shape\n if scale != 1.0:\n image = cv2.resize(image, (0,0), fx=scale, fy=scale)\n return image, orig_shape\n\n\ndef read_mask(image_file, image_shape, scale=1.0):\n image_id = get_image_id(image_file)\n train_info = df_train.loc[df_train['id'] == image_id]\n rle = train_info['encoding'].values[0] if len(train_info) > 0 else None\n if rle is not None:\n mask = rle_to_image(rle, (image_shape[1], image_shape[0]))\n if scale != 1.0:\n mask = cv2.resize(mask, (0,0), fx=scale, fy=scale)\n return np.expand_dims(mask,-1)\n else:\n return None \n\n \ndef read_image_mask(image_file, scale=1.0):\n image, image_shape = read_image(image_file, scale)\n mask = read_mask(image_file, image_shape, scale)\n return image, mask\n\n\ndef get_tile(image, mask, x, y, tile_size, scale=1.0):\n x = round(x * scale)\n y = round(y * scale)\n size = int(round(tile_size / 2 * scale))\n image_s = image[y-size:y+size, x-size:x+size, :] \n mask_s = mask[y-size:y+size, x-size:x+size, :]\n return image_s, mask_s\n\n\ndef get_particles(mask, scale=1.0):\n num, labels, stats, centroids = cv2.connectedComponentsWithStats(mask)\n df_particles = pd.DataFrame(dict(zip(['x','y','left','top','width','height','area'],\n [(centroids[1:,0]) / scale,\n (centroids[1:,1]) / scale,\n (stats[1:,cv2.CC_STAT_LEFT]) / scale,\n (stats[1:,cv2.CC_STAT_TOP]) / scale,\n (stats[1:,cv2.CC_STAT_WIDTH]) / scale,\n (stats[1:,cv2.CC_STAT_HEIGHT]) / scale,\n (stats[1:,cv2.CC_STAT_AREA]) / (scale * scale)])))\n df_particles.sort_values(['x','y'], inplace=True, ignore_index=True)\n df_particles['no'] = range(len(df_particles))\n return df_particles\n\n\ndef analyze_image(image_file):\n image_id = get_image_id(image_file)\n image, image_shape = read_image(image_file, glob_scale)\n mask = read_mask(image_file, image_shape, glob_scale)\n\n mask_full = read_mask(image_file, image_shape, scale=1.0)\n df_glom = get_particles(mask_full, scale=1.0)\n df_glom['id'] = image_id\n del mask_full\n gc.collect()\n \n info = df_info[df_info['image_file'] == f'{image_id}.tiff']\n print(f'Image ID: {image_id:}')\n print(f'Image Size: {info[\"width_pixels\"].values[0]} x {info[\"height_pixels\"].values[0]}')\n print(f'Patient No: {info[\"patient_number\"].values[0]}')\n print(f'Sex: {info[\"sex\"].values[0]}')\n print(f'Age: {info[\"age\"].values[0]}')\n print(f'Race: {info[\"race\"].values[0]}')\n print(f'Height: {info[\"height_centimeters\"].values[0]} cm')\n print(f'Weight: {info[\"weight_kilograms\"].values[0]} kg')\n print(f'BMI: {info[\"bmi_kg/m^2\"].values[0]} kg/m^2')\n print(f'Laterality: {info[\"laterality\"].values[0]}')\n print(f'Percent Cortex: {info[\"percent_cortex\"].values[0]} %')\n print(f'Percent Medulla: {info[\"percent_medulla\"].values[0]} %')\n \n # Plot full image\n if plot_full_image:\n scale = 0.1\n image_small = cv2.resize(image, (0,0), fx=scale, fy=scale)\n mask_small = cv2.resize(mask, (0,0), fx=scale, fy=scale)\n mask_small = np.expand_dims(mask_small,-1) \n \n plt.figure(figsize=(16, 16))\n plt.imshow(overlay_image_mask(image_small, mask_small))\n plt.axis('off')\n\n # Plot glomeruli images\n fig_cols = 5\n fig_rows = int(math.ceil(num_glom_display/fig_cols))\n plt.figure(figsize=(4 * fig_cols, 4 * fig_rows))\n if num_glom_save > 0 and not os.path.exists(image_id):\n os.mkdir(image_id)\n for i in range(min(max(num_glom_display, num_glom_save), len(df_glom))):\n image_s, mask_s = get_tile(image,mask, df_glom['x'][i], df_glom['y'][i], 1000, scale=glob_scale)\n ovl = overlay_image_mask(image_s, mask_s)\n if i < num_glom_display:\n plt.subplot(fig_rows, fig_cols, i+1)\n plt.imshow(ovl)\n plt.axis('off')\n if i < num_glom_save:\n cv2.imwrite(f'{image_id}_{i:03}.png', cv2.cvtColor(ovl, cv2.COLOR_RGB2BGR)) \n \n del image, mask\n gc.collect()\n return df_glom\n\n\ndef plot_glom(df, image_id, glom_no):\n image, mask = read_image_mask(os.path.join(base_path, f'train/{image_id}.tiff'), scale=glob_scale)\n glom = df.loc[(df['id'] == image_id) & (df['no'] == glom_no)]\n im, ma = get_tile(image, mask, glom['x'].iloc[0], glom['y'].iloc[0], 1000, scale=glob_scale)\n del image, mask\n gc.collect()\n plt.figure(figsize=(16,8))\n plt.imshow(overlay_image_mask_original(im, ma))\n plt.title(f'Image: {image_id}, Glomeruli No: {glom_no}, Area: {glom[\"area\"].iloc[0]}')","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:35:27.135608Z","iopub.execute_input":"2024-01-15T13:35:27.135985Z","iopub.status.idle":"2024-01-15T13:35:27.174693Z","shell.execute_reply.started":"2024-01-15T13:35:27.135942Z","shell.execute_reply":"2024-01-15T13:35:27.173763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Execute code to define df_glom (output of analyze_image)\ndf_glom = analyze_image(\"/kaggle/input/hubmap-kidney-segmentation/train/1e2425f28.tiff\")\n\n# Now you can call visualize_tissue_depth using df_glom\nvisualize_tissue_depth(df_glom)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:35:29.685274Z","iopub.execute_input":"2024-01-15T13:35:29.685618Z","iopub.status.idle":"2024-01-15T13:36:21.492247Z","shell.execute_reply.started":"2024-01-15T13:35:29.685587Z","shell.execute_reply":"2024-01-15T13:36:21.490861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\n\ndef artistic_mask_representation(mask):\n # Assuming mask is a binary mask\n artistic_mask = np.zeros_like(mask, dtype=np.uint8)\n\n # Your implementation for creating an artistic representation goes here\n # You might consider using image processing techniques, filters, or other methods\n\n return artistic_mask\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:36:30.317216Z","iopub.execute_input":"2024-01-15T13:36:30.317575Z","iopub.status.idle":"2024-01-15T13:36:30.322136Z","shell.execute_reply.started":"2024-01-15T13:36:30.317541Z","shell.execute_reply":"2024-01-15T13:36:30.321048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_novelty(image, mask):\n fig, axes = plt.subplots(1, 5, figsize=(25, 5))\n\n # Original Image\n axes[0].imshow(image)\n axes[0].set_title('Original Image')\n axes[0].axis('off')\n\n # Masked Image\n masked_image = overlay_image_mask(image, mask)\n axes[1].imshow(masked_image)\n axes[1].set_title('Masked Image')\n axes[1].axis('off')\n\n # Artistic Mask Representation\n artistic_mask = artistic_mask_representation(mask)\n axes[2].imshow(artistic_mask, cmap='gray')\n axes[2].set_title('Artistic Mask Representation')\n axes[2].axis('off')\n\n # Unique Mask Visualization\n unique_mask_values = np.unique(mask)\n unique_mask_colors = plt.cm.get_cmap('tab10', len(unique_mask_values))\n colored_mask = color.label2rgb(mask.squeeze(), colors=unique_mask_colors.colors, bg_label=0)\n axes[3].imshow(colored_mask)\n axes[3].set_title('Unique Mask Visualization')\n axes[3].axis('off')\n\n # Outer Layer of the Image (Assuming a standard RGB image)\n axes[4].imshow(image)\n axes[4].set_title('Outer Layer (Full Image)')\n axes[4].axis('off')\n\n plt.show()\n\n# Assuming image and mask are loaded from your data\nimage, mask = read_image_mask(\"/kaggle/input/hubmap-kidney-segmentation/train/b9a3865fc.tiff\", scale=glob_scale)\n\n# Visualize novel visualizations\nvisualize_novelty(image, mask)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:37:11.521378Z","iopub.execute_input":"2024-01-15T13:37:11.521738Z","iopub.status.idle":"2024-01-15T13:38:40.014519Z","shell.execute_reply.started":"2024-01-15T13:37:11.521704Z","shell.execute_reply":"2024-01-15T13:38:40.013548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_novelty(image, mask):\n fig, axes = plt.subplots(1, 5, figsize=(25, 5))\n\n # Original Image\n axes[0].imshow(image)\n axes[0].set_title('Original Image')\n axes[0].axis('off')\n\n # Masked Image\n masked_image = overlay_image_mask(image, mask)\n axes[1].imshow(masked_image)\n axes[1].set_title('Masked Image')\n axes[1].axis('off')\n\n # Artistic Mask Representation\n artistic_mask = artistic_mask_representation(mask)\n axes[2].imshow(artistic_mask, cmap='gray')\n axes[2].set_title('Artistic Mask Representation')\n axes[2].axis('off')\n\n # Unique Mask Visualization\n unique_mask_values = np.unique(mask)\n unique_mask_colors = plt.cm.get_cmap('tab10', len(unique_mask_values))\n colored_mask = color.label2rgb(mask.squeeze(), colors=unique_mask_colors.colors, bg_label=0)\n axes[3].imshow(colored_mask)\n axes[3].set_title('Unique Mask Visualization')\n axes[3].axis('off')\n\n # Depth Map\n depth_map = np.mean(image, axis=-1) # You can customize this based on your image structure\n axes[4].imshow(depth_map, cmap='viridis') # Using viridis colormap for depth\n axes[4].set_title('Depth Map')\n axes[4].axis('off')\n\n plt.show()\n\n# Assuming image and mask are loaded from your data\nimage, mask = read_image_mask(\"/kaggle/input/hubmap-kidney-segmentation/train/b9a3865fc.tiff\", scale=glob_scale)\n\n# Visualize novel visualizations\nvisualize_novelty(image, mask)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:43:19.598215Z","iopub.execute_input":"2024-01-15T13:43:19.59856Z","iopub.status.idle":"2024-01-15T13:44:43.041296Z","shell.execute_reply.started":"2024-01-15T13:43:19.598532Z","shell.execute_reply":"2024-01-15T13:44:43.040387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\n\ndef generate_synthetic_glomeruli(image_size=(512, 512), num_glomeruli=10, max_radius=30):\n image = np.zeros(image_size, dtype=np.uint8)\n\n for _ in range(num_glomeruli):\n center = np.random.randint(max_radius, image_size[0] - max_radius, size=2)\n radius = np.random.randint(5, max_radius)\n y, x = np.ogrid[:image_size[0], :image_size[1]]\n mask = (x - center[1]) ** 2 + (y - center[0]) ** 2 <= radius ** 2\n image[mask] = 1\n\n return image\n\ndef overlay_image_mask(image, mask, mask_color=(0, 255, 0), alpha=0.5):\n overlay = np.zeros_like(image, dtype=np.uint8)\n overlay[mask] = 255 # Set the intensity to 255 where the mask is true\n return cv2.addWeighted(image, 1 - alpha, overlay, alpha, 0)\n\n# Generate synthetic data\nsynthetic_image = generate_synthetic_glomeruli()\nsynthetic_mask = synthetic_image.astype(bool)\n\n# Overlay mask on the original image\noverlay_result = overlay_image_mask(synthetic_image, synthetic_mask)\n\n# Plot the results\nplt.figure(figsize=(12, 4))\n\nplt.subplot(1, 3, 1)\nplt.imshow(synthetic_image, cmap='gray')\nplt.title('Original Image')\nplt.axis('off')\n\nplt.subplot(1, 3, 2)\nplt.imshow(synthetic_mask, cmap='gray')\nplt.title('Synthetic Mask')\nplt.axis('off')\n\nplt.subplot(1, 3, 3)\nplt.imshow(overlay_result, cmap='gray') # Display in grayscale since it's a binary mask\nplt.title('Overlay Result')\nplt.axis('off')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:46:25.703716Z","iopub.execute_input":"2024-01-15T13:46:25.704101Z","iopub.status.idle":"2024-01-15T13:46:25.99674Z","shell.execute_reply.started":"2024-01-15T13:46:25.704071Z","shell.execute_reply":"2024-01-15T13:46:25.995875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport cv2\n\ndef generate_synthetic_glomeruli(image_size=(512, 512), num_glomeruli=10, max_radius=30):\n image = np.zeros(image_size, dtype=np.uint8)\n\n for _ in range(num_glomeruli):\n center = np.random.randint(max_radius, image_size[0] - max_radius, size=2)\n radius = np.random.randint(5, max_radius)\n y, x = np.ogrid[:image_size[0], :image_size[1]]\n mask = (x - center[1]) ** 2 + (y - center[0]) ** 2 <= radius ** 2\n image[mask] = 1\n\n return image\n\ndef overlay_image_mask(image, mask, mask_color=(0, 255, 0), alpha=0.5):\n overlay = np.zeros_like(image)\n \n if len(image.shape) == 2: # Grayscale image\n overlay[mask] = mask_color[0]\n else: # RGB image\n for c in range(image.shape[2]):\n overlay[mask, c] = mask_color[c]\n\n return cv2.addWeighted(image, 1 - alpha, overlay, alpha, 0)\n\n\n\n# Generate synthetic data\nsynthetic_image = generate_synthetic_glomeruli()\nsynthetic_mask = synthetic_image.astype(bool)\n\n# Overlay mask on the original image\noverlay_result = overlay_image_mask(synthetic_image, synthetic_mask)\n\n# Plot the results\nplt.figure(figsize=(12, 4))\n\nplt.subplot(1, 3, 1)\nplt.imshow(synthetic_image, cmap='gray')\nplt.title('Original Image')\nplt.axis('off')\n\nplt.subplot(1, 3, 2)\nplt.imshow(synthetic_mask, cmap='gray')\nplt.title('Synthetic Mask')\nplt.axis('off')\n\nplt.subplot(1, 3, 3)\nplt.imshow(overlay_result, cmap='gray') # Display in grayscale since it's a binary mask\nplt.title('Overlay Result')\nplt.axis('off')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:46:51.039145Z","iopub.execute_input":"2024-01-15T13:46:51.039528Z","iopub.status.idle":"2024-01-15T13:46:51.313071Z","shell.execute_reply.started":"2024-01-15T13:46:51.039497Z","shell.execute_reply":"2024-01-15T13:46:51.312177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = pd.DataFrame()\ndf_glom = df_glom.append(analyze_image(train_files[0]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:47:06.75908Z","iopub.execute_input":"2024-01-15T13:47:06.759423Z","iopub.status.idle":"2024-01-15T13:47:54.987415Z","shell.execute_reply.started":"2024-01-15T13:47:06.759391Z","shell.execute_reply":"2024-01-15T13:47:54.986575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = pd.DataFrame()\ndf_glom = df_glom.append(analyze_image(train_files[0]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:51:23.608334Z","iopub.execute_input":"2024-01-15T13:51:23.608687Z","iopub.status.idle":"2024-01-15T13:51:54.213512Z","shell.execute_reply.started":"2024-01-15T13:51:23.608658Z","shell.execute_reply":"2024-01-15T13:51:54.212589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[1]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:49:30.126281Z","iopub.execute_input":"2024-01-15T13:49:30.126628Z","iopub.status.idle":"2024-01-15T13:50:46.712948Z","shell.execute_reply.started":"2024-01-15T13:49:30.126592Z","shell.execute_reply":"2024-01-15T13:50:46.711914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[1]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.94539Z","iopub.status.idle":"2024-01-13T10:53:13.945871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(5):\n plot_glom(df_glom, df_glom['id'].iloc[i], df_glom['no'].iloc[i])","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:51:54.21512Z","iopub.execute_input":"2024-01-15T13:51:54.215509Z","iopub.status.idle":"2024-01-15T13:52:54.725713Z","shell.execute_reply.started":"2024-01-15T13:51:54.215471Z","shell.execute_reply":"2024-01-15T13:52:54.724835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[2]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:52:54.727496Z","iopub.execute_input":"2024-01-15T13:52:54.727845Z","iopub.status.idle":"2024-01-15T13:53:29.664447Z","shell.execute_reply.started":"2024-01-15T13:52:54.72781Z","shell.execute_reply":"2024-01-15T13:53:29.663561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[3]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-15T13:53:29.665715Z","iopub.execute_input":"2024-01-15T13:53:29.666017Z","iopub.status.idle":"2024-01-15T13:54:57.856915Z","shell.execute_reply.started":"2024-01-15T13:53:29.665987Z","shell.execute_reply":"2024-01-15T13:54:57.856032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[4]), ignore_index=True)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.951614Z","iopub.status.idle":"2024-01-13T10:53:13.952131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[5]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.953073Z","iopub.status.idle":"2024-01-13T10:53:13.953543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[6]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.954523Z","iopub.status.idle":"2024-01-13T10:53:13.955051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[7]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.956076Z","iopub.status.idle":"2024-01-13T10:53:13.956549Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[8]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.957666Z","iopub.status.idle":"2024-01-13T10:53:13.958204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[9]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.959395Z","iopub.status.idle":"2024-01-13T10:53:13.959875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[10]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.96089Z","iopub.status.idle":"2024-01-13T10:53:13.961377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[11]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.962347Z","iopub.status.idle":"2024-01-13T10:53:13.962801Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[12]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.963702Z","iopub.status.idle":"2024-01-13T10:53:13.964171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[13]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.965051Z","iopub.status.idle":"2024-01-13T10:53:13.965491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[14]), ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.966401Z","iopub.status.idle":"2024-01-13T10:53:13.966841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Glomerulis","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.96811Z","iopub.status.idle":"2024-01-13T10:53:13.968732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_glom.to_csv('glomeruli.csv')\ndisplay(df_glom)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.969822Z","iopub.status.idle":"2024-01-13T10:53:13.970316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\n\ndef generate_synthetic_glomeruli(image_size=(512, 512), num_glomeruli=10, max_radius=30):\n image = np.zeros(image_size, dtype=np.uint8)\n\n for _ in range(num_glomeruli):\n center = np.random.randint(max_radius, image_size[0] - max_radius, size=2)\n radius = np.random.randint(5, max_radius)\n y, x = np.ogrid[:image_size[0], :image_size[1]]\n mask = (x - center[1]) ** 2 + (y - center[0]) ** 2 <= radius ** 2\n image[mask] = 1\n\n return image\n\ndef overlay_image_mask(image, mask, mask_color=(0, 255, 0), alpha=0.5):\n overlay = np.zeros_like(image, dtype=np.uint8)\n overlay[mask] = 255 # Set the intensity to 255 where the mask is true\n return cv2.addWeighted(image, 1 - alpha, overlay, alpha, 0)\n\n# Generate synthetic data\nsynthetic_image = generate_synthetic_glomeruli()\nsynthetic_mask = synthetic_image.astype(bool)\n\n# Overlay mask on the original image\noverlay_result = overlay_image_mask(synthetic_image, synthetic_mask)\n\n# Plot the results\nplt.figure(figsize=(12, 4))\n\nplt.subplot(1, 3, 1)\nplt.imshow(synthetic_image, cmap='gray')\nplt.title('Original Image')\nplt.axis('off')\n\nplt.subplot(1, 3, 2)\nplt.imshow(synthetic_mask, cmap='gray')\nplt.title('Synthetic Mask')\nplt.axis('off')\n\nplt.subplot(1, 3, 3)\nplt.imshow(overlay_result, cmap='gray') # Display in grayscale since it's a binary mask\nplt.title('Overlay Result')\nplt.axis('off')\n\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.971237Z","iopub.status.idle":"2024-01-13T10:53:13.971673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this competation i used segmentation-models-pytorch :\n* [https://github.com/qubvel/segmentation_models.pytorch](http://) . \n* https://smp.readthedocs.io/en/latest/install.html \n\nIt contains: \n* 9 models architectures for binary and multi class segmentation such as Unet,FPN,Unet++ ,...\n* 104 available encoders and thier pre-trained weights for faster and better convergence\n* Loss functions such as dice ,BCE,jaccard,focal loss ,... \n \nThe offline instalation is required in this competation because you can not access internet when you submit the kerenl.
\nrefer to this notebook for the offline instalation of segmentation-models-pytorch :
https://www.kaggle.com/raghaw/offline-install-segmentation-model-pytorch\n","metadata":{}},{"cell_type":"code","source":"#!pip install segmentation-models-pytorch\n!pip install ../input/segmentation-models-pytorch-0-1-3/pretrainedmodels-0.7.4/pretrainedmodels-0.7.4\n!pip install ../input/segmentation-models-pytorch-0-1-3/efficientnet_pytorch-0.6.3/efficientnet_pytorch-0.6.3\n!pip install ../input/segmentation-models-pytorch-0-1-3/timm-0.3.2-py3-none-any.whl\n!pip install ../input/segmentation-models-pytorch-0-1-3/segmentation_models.pytorch.0.1.3/segmentation_models.pytorch.0.1.3","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.97253Z","iopub.status.idle":"2024-01-13T10:53:13.972988Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport torch\nfrom torch import nn\nimport torchvision\nimport cv2\nfrom torchvision import transforms\nfrom torch.utils.data import Dataset,DataLoader\nfrom torch.utils.data.sampler import SequentialSampler, RandomSampler\nfrom torch.optim import Adam\nfrom torch.optim.lr_scheduler import ReduceLROnPlateau\nfrom scipy.ndimage.interpolation import zoom\nimport albumentations as A\nfrom torch.nn import functional as F\nfrom sklearn.model_selection import GroupKFold\nfrom tqdm.notebook import tqdm\n\nimport time\nimport random\nfrom albumentations.pytorch import ToTensorV2\n#from segmentation_models_pytorch.unet import Unet","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.973868Z","iopub.status.idle":"2024-01-13T10:53:13.974334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# For REPRODUCIBILITY https://pytorch.org/docs/stable/notes/randomness.html\nimport os\nimport numpy as np\ndef seed_everything(seed):\n \n random.seed(seed)\n os.environ['PYTHONHASHSEED'] = str(seed)\n np.random.seed(seed)\n torch.manual_seed(seed)\n torch.cuda.manual_seed_all(seed)\n \nseed_everything(42)\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.975351Z","iopub.status.idle":"2024-01-13T10:53:13.97583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class CONFIG():\n \n # data paths\n data_path = '../input/hubmap-256x256/' \n image_data_path='../input/hubmap-256x256/train/'\n mask_data_path='../input/hubmap-256x256/masks/'\n test_data='../input/hubmap-kidney-segmentation/test/' \n model_file = '../input/humap-fpn-dice-focal-15-epochs/FOLD-0-model.pth'\n\n reduce = 4 #reduce the original images by 4 times\n tile_size = 256\n padding = (100,100) # Border overlap for prediction\n\n # pytorch model (https://github.com/qubvel/segmentation_models.pytorch)\n encoder_name = 'timm-efficientnet-b5'#\"efficientnet-b4\"#efficientnet-b3\n encoder_weights = None\n in_channels = 3\n classes = 1\n \n # dataloader \n train_batch_size = 16\n val_batch_size = 4\n test_batch_size=64\n no_folds=5\n # prediction threshold\n threshold = 0.43\n #\n epochs=15\ncfg = CONFIG()","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.97678Z","iopub.status.idle":"2024-01-13T10:53:13.977296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Trainig Dataset \nThe images in the train folder have large size to train deep learning model reasonably on them we have to split this image into tiles of specific size for example 256X256 .
\nFor example if the image is 5000x5000 , we split this images into tiles of size 256x256, we divide 5000/256=19.53 then we ceil it to 20 so we have 20x20 grid with tile size 256x256 . we pad any image to be of size 256x256.
\nThe dataset i used for training is 256x256 dataset which is the ouput of spliting the image into grid of tiles of size 256x256
\nRefer to this notebook : https://www.kaggle.com/iafoss/256x256-images
\n\n","metadata":{}},{"cell_type":"code","source":"# pytorch custom dataset if your are new to pytorch dataset \n#refer to this link https://pytorch.org/tutorials/beginner/data_loading_tutorial.html\nclass HuBMAPDataset(Dataset):\n def __init__(self, ids, phase):\n self.ids = ids\n if phase=='train':\n self.transform = get_train_transform()\n else:\n self.transform = get_val_transform()\n \n def __getitem__(self, idx):\n name = self.ids[idx]\n #img = cv2.imread(f\"../input/hubmap-256x256/train/{name}\").astype(\"float32\")\n img = cv2.imread(cfg.image_data_path+name).astype(\"float32\")\n\n img /= 255.\n #mask = cv2.imread(f\"../input/hubmap-256x256/masks/{name}\")[:,:,0:1]\n mask = cv2.imread(cfg.mask_data_path+name)[:,:,0:1]\n\n transformed = self.transform(image=img, mask=mask)\n img = transformed['image']\n mask = transformed['mask']\n img = img.transpose(2,0,1).astype('float32')\n mask = mask.transpose(2,0,1).astype('float32')\n return img, mask\n\n def __len__(self):\n return len(self.ids)\n \n# the albumentations libarary of data augmentation\n# albumentations : https://github.com/albumentations-team/albumentations \ndef get_train_transform():\n return A.Compose([\n A.HorizontalFlip(),\n A.VerticalFlip(),\n A.OneOf([\n A.RandomContrast(),\n A.RandomGamma(),\n A.RandomBrightness(),\n ], p=0.3),\n A.OneOf([\n A.ElasticTransform(alpha=120, sigma=120 * 0.05, alpha_affine=120 * 0.03),\n A.GridDistortion(),\n A.OpticalDistortion(distort_limit=2, shift_limit=0.5),\n ], p=0.3),\n A.ShiftScaleRotate(p=0.2),\n A.OneOf([\n A.HueSaturationValue(10,15,10),\n A.RandomBrightnessContrast(), \n ], p=0.3),\n \n A.Resize(256,256,always_apply=True),\n ],p=1.)\n\ndef get_val_transform():\n return A.Compose([\n A.Resize(256,256,always_apply=True),\n ],p=1.)\n ","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.97841Z","iopub.status.idle":"2024-01-13T10:53:13.978894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The Segmenation Model\nThere are a lot of available options for the model to use in this competation that are available in Segmentation Models library:\n* 📦 Segmentation Models\n* Unet\n* Unet++\n* MAnet\n* Linknet\n* FPN\n* PSPNet\n* PAN\n* DeepLabV3\n* DeepLabV3+ \n\nYour chossing of specific model to experiment with usually depends on your familiarity with problem at hand (for example if your are new to segmentation you will start with Unet) and the available resources and time you have for training the model.
\n\nThere are a lot of available options for the model you can use as encoder for segmentation model in Segmentation Models library:
\n🏔 Available Encoders\n* ResNet\n* ResNeXt\n* ResNeSt\n* Res2Ne(X)t\n* RegNet(x/y)\n* SE-Net\n* SK-ResNe(X)t\n* DenseNet\n* Inception\n* EfficientNet\n* MobileNet\n* DPN\n* VGG\n\n","metadata":{}},{"cell_type":"markdown","source":"I experimented with Unet and FPN as segmentation model \n* Unet Paper : https://arxiv.org/pdf/1505.04597.pdf\n* FPN paper : https://arxiv.org/pdf/1612.03144.pdf\n\nFPN and Unet are so similiar but their is subtle difference :\nThe Feature Pyramid Network (FPN) looks a lot like the [U-net] (https://vitalab.github.io/deep-learning/2017/02/27/unet.html 221). The main difference is that there is multiple prediction layers: one for each upsampling layer. Like the U-Net, the FPN has laterals connection between the bottom-up pyramid (left) and the top-down pyramid (right). But, where U-net only copy the features and append them, FPN apply a 1x1 convolution layer before adding them. This allows the bottom-up pyramid called “backbone” to be pretty much whatever you want.\nrefer to this link : https://forums.fast.ai/t/differences-between-u-net-and-feature-pyramid-networks-fpn/33166/3\n\nI experimented with ResNet and EfficientNet as encoder :\n* ResNet Paper : https://arxiv.org/pdf/1512.03385\n* EfficientNet paper : https://arxiv.org/pdf/1905.11946\n\n\n\n\n\n\n","metadata":{}},{"cell_type":"code","source":"import segmentation_models_pytorch as smp\nclass HuBMAP(nn.Module):\n def __init__(self):\n super(HuBMAP, self).__init__()\n #self.cnn_model = smp.Unet(cfg.encoder_name, classes=cfg.classes, encoder_weights=cfg.encoder_weights)\n self.cnn_model = smp.FPN(encoder_name=cfg.encoder_name,classes=cfg.classes, encoder_weights=cfg.encoder_weights)\n \n \n def forward(self, imgs):\n img_segs = self.cnn_model(imgs)\n return img_segs\n ","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.979893Z","iopub.status.idle":"2024-01-13T10:53:13.98046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loss Function \nrefer to thoes links:\n* https://ilmonteux.github.io/2019/05/10/segmentation-metrics.html\n* https://www.kaggle.com/yassinealouini/all-the-segmentation-metrics\n* https://www.kaggle.com/c/hubmap-kidney-segmentation/discussion/208413\n\nLoss functions that are used in semantic segmentation can generaly be divided into two families:\n* Intersection-Over-Union (IOU) Jaccard and Dice family\n* Cross Entropy family
\nas described in the image below \n\n\n* (IOU) Intersection-Over-Union,Jaccard index and Dice family\n\n\n\n*IoU calculation visualized*\n\n$ Dice(A, B) = \\frac{2 \\|A \\cap B\\|}{\\|A\\|+\\|B\\|}, \\qquad Jaccard(A, B) = \\frac{\\|A \\cap B\\|}{\\|A \\cup B\\|} $\n\n\n\n* Cross Entropy family \n * Binary cross entropy : $ -{(y\\log(p) + (1 - y)\\log(1 - p))} $
\n * Cross entropy : $ -{(log(p_{t}))} $ , where $ p_{t} $ = the probability of ground truth class\n * Focal loss : $ -{\\alpha(1-p_{t})^{\\gamma}(log(p_{t}))} $\n *  \n \n* We can achieve the best result by combining loss functions from thoes families .\n I experimented with \n * (dice,cross entropy) \n * (jaccard,cross entropy) \n * (dice,focal loss) \n * (jaccard,focal loss) \n\n* Loss functions available in Segmentation library
\n refer to this link : https://smp.readthedocs.io/en/latest/losses.html
\n 📉 Losses\n * Constants\n *JaccardLoss\n *DiceLoss\n *FocalLoss\n *LovaszLoss\n *SoftBCEWithLogitsLoss\n *SoftCrossEntropyLoss","metadata":{}},{"cell_type":"code","source":"from segmentation_models_pytorch.losses import LovaszLoss\nfrom segmentation_models_pytorch.losses import JaccardLoss\nfrom segmentation_models_pytorch.losses import DiceLoss\nfrom segmentation_models_pytorch.losses import FocalLoss\nfrom segmentation_models_pytorch.losses import SoftBCEWithLogitsLoss\n\nlovasz=LovaszLoss(mode='binary')\njaccard=JaccardLoss(mode='binary')\ndice=DiceLoss(mode='binary')\nfocal_loss=FocalLoss(mode='binary')\nBCE=SoftBCEWithLogitsLoss()\n\ndef symmetric_lovasz(outputs, targets):\n return 0.5*(lovasz(outputs, targets) + lovasz(-outputs, 1.0 - targets))\n\n\nclass DiceLoss_(nn.Module):\n def __init__(self, weight=None, size_average=True):\n super(DiceLoss_, self).__init__()\n\n def forward(self, inputs, targets, smooth=1):\n \n #comment out if your model contains a sigmoid or equivalent activation layer\n inputs = F.sigmoid(inputs) \n \n #flatten label and prediction tensors\n inputs = inputs.view(-1)\n targets = targets.view(-1)\n \n intersection = (inputs * targets).sum() \n #dice = (2.*intersection + smooth)/(inputs.sum() + targets.sum() + smooth) \n dice = (2.*intersection + smooth)/(inputs.sum() + targets.sum() + smooth) \n \n return 1 - dice\n \n \n \nclass DiceBCELoss(nn.Module):\n # Formula Given above.\n def __init__(self, weight=None, size_average=True):\n super(DiceBCELoss, self).__init__()\n\n def forward(self, inputs, targets, smooth=1):\n \n #comment out if your model contains a sigmoid or equivalent activation layer\n inputs = F.sigmoid(inputs) \n \n #flatten label and prediction tensors\n inputs = inputs.view(-1)\n targets = targets.view(-1)\n \n intersection = (inputs * targets).sum() \n dice_loss = 1 - (2.*intersection + smooth)/(inputs.sum() + targets.sum() + smooth) \n BCE = F.binary_cross_entropy(inputs, targets, reduction='mean')\n Dice_BCE = BCE + dice_loss\n \n return Dice_BCE\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.981565Z","iopub.status.idle":"2024-01-13T10:53:13.982111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef HuBMAPLoss(images, targets, model,loss_fun_0,loss_fun_1,device):\n \n model.to(device)\n \n images = images.to(device)\n targets = targets.to(device)\n \n outputs = model(images)\n loss_0 = loss_fun_0(outputs, targets)\n \n if loss_fun_1 != None:\n loss_1=loss_fun_1(outputs, targets)\n loss=loss_0 + loss_1\n else:\n loss=loss_0\n \n \n return loss, outputs\n\ndef train_one_epoch(epoch, model, device, optimizer, scheduler,loss_fun_0,loss_fun_1,trainloader):\n model.train()\n t = time.time()\n total_loss = 0\n tk0=tqdm(trainloader)\n for step, (images, targets) in enumerate(tk0):\n loss, outputs = HuBMAPLoss(images, targets, model,loss_fun_0,loss_fun_1,device)\n optimizer.zero_grad()\n loss.backward()\n #if ((step+1)%4==0 or (step+1)==len(trainloader)):\n optimizer.step()\n scheduler.step()\n #optimizer.zero_grad()\n loss = loss.detach().item()\n total_loss += loss\n tk0.set_postfix(loss=(total_loss / ((step+1) * trainloader.batch_size)))\n\n #if ((step+1)%10==0 or (step+1)==len(trainloader)):\n epoch_loss=total_loss/(len(trainloader)*trainloader.batch_size) \n elapsed_time=time.time() - t\n return epoch_loss,elapsed_time\n\n \n \ndef valid_one_epoch(epoch, model, device, optimizer, scheduler,loss_fun_0,loss_fun_1, validloader):\n model.eval()\n t = time.time()\n total_loss = 0\n tk0=tqdm(validloader)\n running_jac=0\n running_dice=0\n for step, (images, targets) in enumerate(tk0):\n #with torch.no_grad():\n\n loss, outputs = HuBMAPLoss(images, targets, model,loss_fun_0,loss_fun_1,device)\n loss = loss.detach().item()\n total_loss += loss\n tk0.set_postfix(loss=(total_loss / ((step+1) * validloader.batch_size)))\n jac_metric=jaccard(outputs.to(device),targets.to(device))\n dice_metric=dice(outputs.to(device),targets.to(device))\n jac_metric=1-jac_metric.item()\n dice_metric=1-dice_metric.item()\n running_jac+=jac_metric\n running_dice+=dice_metric\n #if ((step+1)%4==0 or (step+1)==len(validloader)):\n # scheduler.step(total_loss/len(validloader))\n #if ((step+1)%10==0 or (step+1)==len(validloader)):\n val_loss= total_loss/(len(validloader)*validloader.batch_size) \n Jac_metric= running_jac/(len(validloader)*validloader.batch_size) \n Dice_metric= running_dice/(len(validloader)*validloader.batch_size) \n \n elapsed_time=time.time() - t\n return val_loss,Jac_metric,Dice_metric,elapsed_time\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.983179Z","iopub.status.idle":"2024-01-13T10:53:13.9837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef prepare_train_valid_dataloader(df, fold):\n train_ids = df.loc[~df.Folds.isin(fold), \"Image_Paths\"].values\n val_ids = df.loc[df.Folds.isin(fold), \"Image_Paths\"].values\n train_ds = HuBMAPDataset(train_ids, \"train\")\n val_ds = HuBMAPDataset(val_ids, \"val\")\n train_loader = DataLoader(train_ds, batch_size=cfg.train_batch_size, pin_memory=True, shuffle=True, num_workers=4)\n val_loader = DataLoader(val_ds, batch_size=cfg.val_batch_size, pin_memory=True, shuffle=False, num_workers=4)\n return train_loader, val_loader\n ","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.984734Z","iopub.status.idle":"2024-01-13T10:53:13.985237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef load_checkpoint(path,model,optimizer,scheduler,device):\n checkpoint = torch.load(path,map_location=device)\n model.load_state_dict(checkpoint['model_state_dict'])\n if optimizer!=None and scheduler !=None:\n optimizer.load_state_dict(checkpoint['optimizer_state_dict'])\n scheduler.load_state_dict(checkpoint['scheduler_state_dict']) \n epoch = checkpoint['epoch']\n val_loss = checkpoint['val_loss']\n jacard_metric=checkpoint['jacard_metric']\n dice_metric=checkpoint['dice_metric']\n \n print(f'epoch={epoch},val_loss={val_loss} ,jacard_metric={jacard_metric} ,dice_metric={dice_metric}')\n\ndef load_models(MODELS,device):\n models = []\n for path in MODELS:\n \n \n model = HuBMAP()\n load_checkpoint(path,model,None,None,device)\n \n model.eval()\n model.to(device)\n models.append(model)\n\n print(len(models))\n return models\n\ndef save_checkpoint(epoch,val_loss,jac_metric,dice_metric,model,optimizer,scheduler,path):\n checkpoint={\n 'epoch': epoch,\n 'val_loss':val_loss,\n 'jacard_metric':jac_metric,\n 'dice_metric':dice_metric,\n 'model_state_dict': model.state_dict(),\n 'optimizer_state_dict': optimizer.state_dict(),\n 'scheduler_state_dict':scheduler.state_dict()\n }\n torch.save(checkpoint,path)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.986223Z","iopub.status.idle":"2024-01-13T10:53:13.98671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\ndirectory_list = os.listdir(cfg.image_data_path)\ndir_df = pd.DataFrame(directory_list, columns=['Image_Paths'])\nar=dir_df[dir_df.columns[0]].values\ngroups=[id_[:-9] for id_ in ar]\ndir_df\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.98783Z","iopub.status.idle":"2024-01-13T10:53:13.988353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GroupKFold.html\nFOLDS = cfg.no_folds\ngkf = GroupKFold(FOLDS)\ndir_df['Folds'] = 0\nfor fold, (tr_idx, val_idx) in enumerate(gkf.split(dir_df, groups=groups)):\n dir_df.loc[val_idx, 'Folds'] = fold\ndir_df","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.989394Z","iopub.status.idle":"2024-01-13T10:53:13.989861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{}},{"cell_type":"code","source":"'''\ntrain_from_scratch=False\nload_model_path=cfg.model_file #'../input/humap-fpn-dice-bce-10-epochs/FOLD-0-model.pth'\n\n\nfor fold, (tr_idx, val_idx) in enumerate(gkf.split(dir_df, groups=groups)):\n if fold>0:\n break\n \n trainloader, validloader = prepare_train_valid_dataloader(dir_df, [fold])\n\n loss_fun_0=dice #jaccard#dice#symmetric_lovasz\n loss_fun_1=focal_loss \n \n model = HuBMAP().to(device) \n optimizer = Adam(model.parameters(), lr=3e-3)\n scheduler = torch.optim.lr_scheduler.StepLR(optimizer, gamma=0.1, step_size=5000) \n \n if train_from_scratch==False and load_model_path !=None:\n load_checkpoint(load_model_path,model,optimizer,scheduler,device)\n \n num_epochs = 5 #cfg.epochs\n \n best_mertic=0\n \n for epoch in range(num_epochs):\n train_loss,train_elapsed_time=train_one_epoch(epoch, model, device, optimizer, \n scheduler,loss_fun_0,loss_fun_1,trainloader)\n with torch.no_grad():\n val_loss,Jac_metric,Dice_metric,val_elapsed_time=valid_one_epoch(epoch, model, device, optimizer,\n scheduler,loss_fun_0,loss_fun_1,\n validloader)\n \n print(\n \n f'epoch {epoch} \\n' \n f'train loss: {(train_loss):.4f}, ' + \\\n f'train time: {(train_elapsed_time):.4f} \\n' \n f'val loss: {(val_loss):.4f}, ' + \\\n f'Jacard metric: {(Jac_metric):.4f}, ' + \\\n f'Dice metric: {(Dice_metric):.4f}, ' + \\\n\n f'val time: {(val_elapsed_time):.4f}'\n )\n \n metric=(Jac_metric+Dice_metric)/2\n if metric > best_mertic :\n path=f'FOLD-{fold}-model.pth'\n save_checkpoint(epoch,val_loss,Jac_metric,Dice_metric,model,optimizer,scheduler,path) \n best_mertic=metric\n print(f'Saving Model Dice={Dice_metric} , jaccard={Jac_metric} , metric={metric}' ) \n'''","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.990905Z","iopub.status.idle":"2024-01-13T10:53:13.991406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploring the trained model on the validation dataset","metadata":{}},{"cell_type":"markdown","source":"## Predicted Mask VS Actual Mask","metadata":{}},{"cell_type":"code","source":"\n#path=f'FOLD-{0}-model.pth'\npath='../input/humap-fpn-dice-focal-15-epochs/FOLD-0-model.pth'\nmodel = HuBMAP().to(device)\nload_checkpoint(path,model,None,None,device)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.992398Z","iopub.status.idle":"2024-01-13T10:53:13.992888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trainloader, validloader = prepare_train_valid_dataloader(dir_df, [0])\ndl = iter(validloader)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.99399Z","iopub.status.idle":"2024-01-13T10:53:13.994462Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\nimport matplotlib.pyplot as plt\n\njaccard_=JaccardLoss(mode='binary',log_loss=False)\ndice_=DiceLoss(mode='binary',log_loss=False)\ndice__=DiceLoss_()\nimgs,masks = next(dl)\n\nwith torch.no_grad():\n\n model.to(device).eval()\n images = imgs.to(device)\n outputs = model(images)\n f, axs = plt.subplots(4,3,figsize=(20,20))\n for i,(img,mask,out) in enumerate(zip(imgs,masks,outputs)):\n img = ((img.permute(1,2,0))*255.0).numpy().astype(np.uint8)\n \n out = (torch.sigmoid(out)>.43).type(torch.uint8)\n\n print(1-symmetric_lovasz(mask.float().cpu(),mask.cpu()).item(),\n 1-dice_(mask.float().cpu(),mask.cpu()).item(),\n 1-jaccard_(mask.float().cpu(),mask.cpu()).item(),\n dice__(mask.float().cpu(),mask.cpu()).item()\n )\n \n axs[i][0].imshow(img,vmin=0,vmax=255)\n axs[i][1].imshow(mask.squeeze().cpu().numpy())\n axs[i][2].imshow(out.type(torch.uint8).squeeze().cpu().numpy())\n\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.995751Z","iopub.status.idle":"2024-01-13T10:53:13.996238Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n#del trainloader,validloader,imgs,masks,outputs\nimport gc\n\ngc.collect()\ntorch.cuda.empty_cache()\n'''","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.997301Z","iopub.status.idle":"2024-01-13T10:53:13.99779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inference","metadata":{}},{"cell_type":"code","source":"\npath='../input/humap-fpn-dice-focal-15-epochs/FOLD-0-model.pth'\nmodel = HuBMAP().to(device).eval()\nload_checkpoint(path,model,None,None,device)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:13.998725Z","iopub.status.idle":"2024-01-13T10:53:13.999226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sz = cfg.tile_size#256 #the size of tiles\nreduce = cfg.reduce #4 #reduce the original images by 4 times\nTH = cfg.threshold #0.43 #threshold for positive predictions\nDATA = cfg.test_data #'../input/hubmap-kidney-segmentation/test/'\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.000222Z","iopub.status.idle":"2024-01-13T10:53:14.000699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import rasterio\nfrom rasterio.windows import Window\n#functions to convert encoding to mask and mask to encoding\ndef enc2mask(encs, shape):\n img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n for m,enc in enumerate(encs):\n if isinstance(enc,np.float) and np.isnan(enc): continue\n s = enc.split()\n for i in range(len(s)//2):\n start = int(s[2*i]) - 1\n length = int(s[2*i+1])\n img[start:start+length] = 1 + m\n return img.reshape(shape).T\n\ndef mask2enc(mask, n=1):\n pixels = mask.T.flatten()\n encs = []\n for i in range(1,n+1):\n p = (pixels == i).astype(np.int8)\n if p.sum() == 0: encs.append(np.nan)\n else:\n p = np.concatenate([[0], p, [0]])\n runs = np.where(p[1:] != p[:-1])[0] + 1\n runs[1::2] -= runs[::2]\n encs.append(' '.join(str(x) for x in runs))\n return encs\n\n#https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\n#with transposed mask\ndef rle_encode_less_memory(img):\n #the image should be transposed\n pixels = img.T.flatten()\n \n # This simplified method requires first and last pixel to be zero\n pixels[0] = 0\n pixels[-1] = 0\n runs = np.where(pixels[1:] != pixels[:-1])[0] + 2\n runs[1::2] -= runs[::2]\n \n return ' '.join(str(x) for x in runs)\n# https://www.kaggle.com/iafoss/256x256-images\nmean = np.array([0.65459856,0.48386562,0.69428385])\nstd = np.array([0.15167958,0.23584107,0.13146145])\n\ns_th = 40 #saturation blancking threshold\np_th = 1000*(sz//256)**2 #threshold for the minimum number of pixels\nidentity = rasterio.Affine(1, 0, 0, 0, 1, 0)\n\ndef img2tensor(img,dtype:np.dtype=np.float32):\n if img.ndim==2 : img = np.expand_dims(img,2)\n img = np.transpose(img,(2,0,1))\n return torch.from_numpy(img.astype(dtype, copy=False))\n\n\nclass HuBMAPDataset_test(Dataset):\n def __init__(self, idx, sz=sz, reduce=reduce):\n self.data = rasterio.open(os.path.join(DATA,idx+'.tiff'), transform = identity,\n num_threads='all_cpus')\n # some images have issues with their format \n # and must be saved correctly before reading with rasterio\n if self.data.count != 3:\n subdatasets = self.data.subdatasets\n self.layers = []\n if len(subdatasets) > 0:\n for i, subdataset in enumerate(subdatasets, 0):\n self.layers.append(rasterio.open(subdataset))\n self.shape = self.data.shape\n self.reduce = reduce\n self.sz = reduce*sz\n self.pad0 = (self.sz - self.shape[0]%self.sz)%self.sz\n self.pad1 = (self.sz - self.shape[1]%self.sz)%self.sz\n self.n0max = (self.shape[0] + self.pad0)//self.sz\n self.n1max = (self.shape[1] + self.pad1)//self.sz\n \n def __len__(self):\n return self.n0max*self.n1max\n \n def __getitem__(self, idx):\n # the code below may be a little bit difficult to understand,\n # but the thing it does is mapping the original image to\n # tiles created with adding padding, as done in\n # https://www.kaggle.com/iafoss/256x256-images ,\n # and then the tiles are loaded with rasterio\n # n0,n1 - are the x and y index of the tile (idx = n0*self.n1max + n1)\n n0,n1 = idx//self.n1max, idx%self.n1max\n # x0,y0 - are the coordinates of the lower left corner of the tile in the image\n # negative numbers correspond to padding (which must not be loaded)\n x0,y0 = -self.pad0//2 + n0*self.sz, -self.pad1//2 + n1*self.sz\n # make sure that the region to read is within the image\n p00,p01 = max(0,x0), min(x0+self.sz,self.shape[0])\n p10,p11 = max(0,y0), min(y0+self.sz,self.shape[1])\n img = np.zeros((self.sz,self.sz,3),np.uint8)\n # mapping the loade region to the tile\n if self.data.count == 3:\n img[(p00-x0):(p01-x0),(p10-y0):(p11-y0)] = np.moveaxis(self.data.read([1,2,3],\n window=Window.from_slices((p00,p01),(p10,p11))), 0, -1)\n else:\n for i,layer in enumerate(self.layers):\n img[(p00-x0):(p01-x0),(p10-y0):(p11-y0),i] =\\\n layer.read(1,window=Window.from_slices((p00,p01),(p10,p11)))\n \n if self.reduce != 1:\n img = cv2.resize(img,(self.sz//reduce,self.sz//reduce),\n interpolation = cv2.INTER_AREA)\n #check for empty imges\n hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n h,s,v = cv2.split(hsv)\n if (s>s_th).sum() <= p_th or img.sum() <= p_th:\n #images with -1 will be skipped\n return img2tensor((img/255.0 - mean)/std), -1\n else: return img2tensor((img/255.0 - mean)/std), idx\n#iterator like wrapper that returns predicted masks\nclass Model_pred:\n def __init__(self, models, dl, tta:bool=True, half:bool=False):\n self.models = models\n self.dl = dl\n self.tta = tta\n self.half = half\n \n def __iter__(self):\n count=0\n with torch.no_grad():\n for x,y in iter(tqdm(self.dl)):\n if ((y>=0).sum() > 0): #exclude empty images\n x = x[y>=0].to(device)\n y = y[y>=0]\n if self.half: x = x.half()\n py = None\n for model in self.models:\n p = model(x)\n p = torch.sigmoid(p).detach()\n if py is None: py = p\n else: py += p\n if self.tta:\n #x,y,xy flips as TTA\n flips = [[-1],[-2],[-2,-1]]\n for f in flips:\n xf = torch.flip(x,f)\n for model in self.models:\n p = model(xf)\n p = torch.flip(p,f)\n py += torch.sigmoid(p).detach()\n py /= (1+len(flips)) \n py /= len(self.models)\n\n py = F.upsample(py, scale_factor=reduce, mode=\"bilinear\")\n py = py.permute(0,2,3,1).float().cpu()\n \n batch_size = len(py)\n for i in range(batch_size):\n yield py[i],y[i]\n count += 1\n \n def __len__(self):\n return len(self.dl.dataset)\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.001849Z","iopub.status.idle":"2024-01-13T10:53:14.002377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\ndf_sample = pd.read_csv('../input/hubmap-kidney-segmentation/sample_submission.csv')\nmodels=[model]","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.003447Z","iopub.status.idle":"2024-01-13T10:53:14.003933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing the predicted mask of one image from the public testset","metadata":{}},{"cell_type":"code","source":"bs=cfg.test_batch_size#64\nnames,preds = [],[]\nidxs=['2ec3f1bb9' ,'3589adb90','d488c759a','aa05346ff','57512b7f1']\nidx = idxs[3]\nprint(idx)\nds = HuBMAPDataset_test(idx)\n#rasterio cannot be used with multiple workers\ndl = DataLoader(ds,bs,num_workers=0,shuffle=False,pin_memory=True)\nmp = Model_pred(models,dl)\n#generate masks\nmask = torch.zeros(len(ds),ds.sz,ds.sz,dtype=torch.int8)\nfor p,i in iter(mp): mask[i.item()] = p.squeeze(-1) > TH\n\n#reshape tiled masks into a single mask and crop padding\nmask = mask.view(ds.n0max,ds.n1max,ds.sz,ds.sz).\\\n permute(0,2,1,3).reshape(ds.n0max*ds.sz,ds.n1max*ds.sz)\nmask = mask[ds.pad0//2:-(ds.pad0-ds.pad0//2) if ds.pad0 > 0 else ds.n0max*ds.sz,\n ds.pad1//2:-(ds.pad1-ds.pad1//2) if ds.pad1 > 0 else ds.n1max*ds.sz]\n\n#convert to rle\n#https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\nrle = rle_encode_less_memory(mask.numpy())\nnames.append(idx)\npreds.append(rle)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.004911Z","iopub.status.idle":"2024-01-13T10:53:14.005419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mask_np=mask.numpy()\nimage_id = get_image_id(test_files[3])\nimage, image_shape = read_image(test_files[3], glob_scale)\nmask_np=np.expand_dims(mask_np,-1)\nmask_np = cv2.resize(mask_np.astype('float32'), (0,0), fx=glob_scale, fy=glob_scale)\n#mask = np.expand_dims(mask,-1)\n\nscale = 0.1\nimage_small = cv2.resize(image, (0,0), fx=scale, fy=scale)\nmask_small = cv2.resize(mask_np, (0,0), fx=scale, fy=scale)\nmask_small = np.expand_dims(mask_small,-1) \n\nplt.figure(figsize=(16, 16))\nplt.imshow(overlay_image_mask(image_small, mask_small))\nplt.axis('off')\n#del mask, ds, dl\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.006396Z","iopub.status.idle":"2024-01-13T10:53:14.006862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Making the submission file","metadata":{}},{"cell_type":"code","source":"\nnames,preds = [],[]\nbs=cfg.test_batch_size#64\nfor idx,row in tqdm(df_sample.iterrows(),total=len(df_sample)):\n \n idx = row['id']\n print(idx)\n ds = HuBMAPDataset_test(idx)\n #rasterio cannot be used with multiple workers\n dl = DataLoader(ds,bs,num_workers=0,shuffle=False,pin_memory=True)\n mp = Model_pred(models,dl)\n #generate masks\n mask = torch.zeros(len(ds),ds.sz,ds.sz,dtype=torch.int8)\n for p,i in iter(mp): mask[i.item()] = p.squeeze(-1) > TH\n \n #reshape tiled masks into a single mask and crop padding\n mask = mask.view(ds.n0max,ds.n1max,ds.sz,ds.sz).\\\n permute(0,2,1,3).reshape(ds.n0max*ds.sz,ds.n1max*ds.sz)\n mask = mask[ds.pad0//2:-(ds.pad0-ds.pad0//2) if ds.pad0 > 0 else ds.n0max*ds.sz,\n ds.pad1//2:-(ds.pad1-ds.pad1//2) if ds.pad1 > 0 else ds.n1max*ds.sz]\n \n #convert to rle\n #https://www.kaggle.com/bguberfain/memory-aware-rle-encoding\n rle = rle_encode_less_memory(mask.numpy())\n names.append(idx)\n preds.append(rle)\n del mask, ds, dl\n #gc.collect()\n","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.008077Z","iopub.status.idle":"2024-01-13T10:53:14.008574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'id':names,'predicted':preds})\ndf.to_csv('submission.csv',index=False)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.009525Z","iopub.status.idle":"2024-01-13T10:53:14.01002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#from IPython.display import FileLink\n#FileLink(r'FOLD-0-model.pth')","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.010972Z","iopub.status.idle":"2024-01-13T10:53:14.01144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#os.stat('submission.csv').st_size/(1024*1024)","metadata":{"execution":{"iopub.status.busy":"2024-01-13T10:53:14.012397Z","iopub.status.idle":"2024-01-13T10:53:14.012883Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}