{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":6128529,"sourceType":"datasetVersion","datasetId":3418372}],"dockerImageVersionId":30028,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Introduction\n\nThis is an exploratory notebook on the HuBMAP challenge introducing the tool rasterio to read and play with tiff\n\n![img](https://i.imgur.com/F4naKqh.png)\n\n## Part-1 : Crops and masks with rasterio\n\n### a. Load an image with rasterio\n\nLoading an image, checking the coordinate system and the linear transformation\n\n*rasterio*\n\n### b. Visualise Anatomical Structures\n\nWe use the annotated json to visualise the differente structues inside an image\n\n*json, rasterio, rasterio.mask*\n\n### c. Visualise glomerulus\n\nWe use the annotated json to generate masks and visualise the glomeruli\n\n*json, rasterio, rasterio.mask*\n\n### d. Generate overlay mask of Glomerulus on Anatomical Structure with rasterio\n\nIn this last part, we create a full downsized image with masks of glomeruli overlayed\n\n*json, rasterio, rasterio.mask*\n\n## Part-2 : EDA\n\n### a. Usefull surface, Glomeruli surface\n\nSome global stats, id per id:\n\n1. We want to check the total surface in pixel of the images\n2. We want to check the total space covered by glomeruli\n3. Rather than ploting the ratio of surface_glomeruli/total_surface, we check only the \"usefull\" surface, in which there is actually glomeruli\n\n*cv2, rasterio, rasterio.mask, json*\n\n### a. Height, width and surface of the glomeruli in an image\n\nIt is important to know how much space take a full glomerulus in order to calibrate correctly the crops\n\n### b. Size and oriention of glomeruli\n\nAs the glomeruli can have an orientation, we would like as well to know the global dimension of the glomeruli. \n\n### c. Biggest and smallest glomeruli in the set\n\nVisualisation of the extrema. Helpfull to identify possible mistakes in labels","metadata":{}},{"cell_type":"code","source":"# ================================\n# Cell 1: Import Necessary Libraries\n# ================================\nimport os\nimport numpy as np\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nimport rasterio\nfrom glob import glob\nfrom sklearn.model_selection import train_test_split\nimport matplotlib.pyplot as plt\nfrom tqdm import tqdm\nfrom skimage.restoration import denoise_bilateral\nfrom skimage import exposure\nimport cv2\nimport time\n\n# Ensure GPU is used if available\ntorch.backends.cudnn.benchmark = True\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(f\"Using device: {device}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:37.610012Z","iopub.execute_input":"2024-10-01T17:23:37.610383Z","iopub.status.idle":"2024-10-01T17:23:43.426434Z","shell.execute_reply.started":"2024-10-01T17:23:37.610349Z","shell.execute_reply":"2024-10-01T17:23:43.425598Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# =========================================\n# Cell 2: Define SAR and Optical Preprocessing Functions\n# =========================================\ndef preprocess_sar_image(image, target_shape):\n    \"\"\"Preprocess SAR bands and ensure they match the target optical shape.\"\"\"\n    processed_channels = []\n    for i in range(image.shape[-1]):  # Iterate over each channel\n        channel = image[:, :, i]\n        filtered_channel = denoise_bilateral(channel, sigma_color=0.05, sigma_spatial=15, multichannel=False)\n        processed_channels.append(filtered_channel)\n    \n    # Stack the processed channels back into a multi-band image\n    image = np.stack(processed_channels, axis=-1)\n\n    # Resampling to Match Optical Resolution\n    image_resampled = cv2.resize(image, (target_shape[1], target_shape[0]), interpolation=cv2.INTER_LINEAR)\n    return image_resampled\n\n# Optical Preprocessing Function\ndef preprocess_optical_image(image):\n    \"\"\"Apply atmospheric correction using Histogram Equalization.\"\"\"\n    corrected_image = exposure.equalize_adapthist(image, clip_limit=0.03)\n    return corrected_image\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.428551Z","iopub.execute_input":"2024-10-01T17:23:43.42884Z","iopub.status.idle":"2024-10-01T17:23:43.437971Z","shell.execute_reply.started":"2024-10-01T17:23:43.42881Z","shell.execute_reply":"2024-10-01T17:23:43.437119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ====================================\n# Cell 3: Define FloodDataset Class for Image Handling\n# ====================================\nclass FloodDataset(Dataset):\n    def __init__(self, image_paths, mask_paths, transform=None):\n        self.image_paths = image_paths\n        self.mask_paths = mask_paths\n        self.transform = transform\n\n    def __len__(self):\n        return len(self.image_paths)\n\n    def __getitem__(self, idx):\n        img_path = self.image_paths[idx]\n        mask_path = self.mask_paths[idx]\n\n        # Load and preprocess the image\n        with rasterio.open(img_path) as src:\n            image = src.read().transpose(1, 2, 0)  # CHW to HWC\n\n            # Separate Optical and SAR Bands\n            sar_bands = image[..., 3:5]  # Last two bands are SAR-like\n            optical_bands = image[..., :3]  # First three bands are Optical\n\n            # Pass the shape of the optical bands to the SAR preprocessing function\n            sar_bands = preprocess_sar_image(sar_bands, target_shape=optical_bands.shape[:2])\n            optical_bands = preprocess_optical_image(optical_bands)\n\n            # Combine the processed bands (dimensions should now match)\n            image = np.concatenate([optical_bands, sar_bands], axis=-1)\n\n        # Load the mask\n        with rasterio.open(mask_path) as src:\n            mask = src.read(1)\n\n        # Normalize image and mask\n        image = image.astype(np.float32) / 255.0\n        mask = mask.astype(np.float32) / 255.0\n\n        # Convert to PyTorch tensors\n        image = torch.from_numpy(image.transpose(2, 0, 1))  # HWC to CHW\n        mask = torch.from_numpy(mask)\n\n        return image, mask\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.439182Z","iopub.execute_input":"2024-10-01T17:23:43.439464Z","iopub.status.idle":"2024-10-01T17:23:43.452364Z","shell.execute_reply.started":"2024-10-01T17:23:43.439437Z","shell.execute_reply":"2024-10-01T17:23:43.451696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_data_loaders(data_dir, batch_size=32):\n    image_paths = glob(os.path.join(data_dir, 'images', '*.tif'))\n    mask_paths = [path.replace('images', 'annotations') for path in image_paths]\n\n    # Print paths to debug\n    print(f\"Found {len(image_paths)} image files.\")\n    print(f\"Found {len(mask_paths)} mask files.\")\n\n    train_img, test_img, train_mask, test_mask = train_test_split(image_paths, mask_paths, test_size=0.2, random_state=42)\n    train_img, val_img, train_mask, val_mask = train_test_split(train_img, train_mask, test_size=0.2, random_state=42)\n\n    train_dataset = FloodDataset(train_img, train_mask)\n    val_dataset = FloodDataset(val_img, val_mask)\n    test_dataset = FloodDataset(test_img, test_mask)\n\n    train_loader = DataLoader(train_dataset, batch_size=batch_size, shuffle=True, num_workers=2, pin_memory=True)\n    val_loader = DataLoader(val_dataset, batch_size=batch_size, shuffle=False, num_workers=2, pin_memory=True)\n    test_loader = DataLoader(test_dataset, batch_size=batch_size, shuffle=False, num_workers=2, pin_memory=True)\n\n    return train_loader, val_loader, test_loader\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.45372Z","iopub.execute_input":"2024-10-01T17:23:43.453987Z","iopub.status.idle":"2024-10-01T17:23:43.468133Z","shell.execute_reply.started":"2024-10-01T17:23:43.453962Z","shell.execute_reply":"2024-10-01T17:23:43.467484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ===============================\n# Cell 5: Define the U-Net Model\n# ===============================\nclass UNet(nn.Module):\n    def __init__(self, n_channels, n_classes):\n        super(UNet, self).__init__()\n        self.n_channels = n_channels\n        self.n_classes = n_classes\n\n        def double_conv(in_channels, out_channels):\n            return nn.Sequential(\n                nn.Conv2d(in_channels, out_channels, 3, padding=1),\n                nn.ReLU(inplace=True),\n                nn.Conv2d(out_channels, out_channels, 3, padding=1),\n                nn.ReLU(inplace=True)\n            )\n\n        self.dc1 = double_conv(n_channels, 64)\n        self.dc2 = double_conv(64, 128)\n        self.dc3 = double_conv(128, 256)\n        self.dc4 = double_conv(256, 512)\n        self.dc5 = double_conv(512, 1024)\n\n        self.up1 = nn.ConvTranspose2d(1024, 512, 2, stride=2)\n        self.dc6 = double_conv(1024, 512)\n        self.up2 = nn.ConvTranspose2d(512, 256, 2, stride=2)\n        self.dc7 = double_conv(512, 256)\n        self.up3 = nn.ConvTranspose2d(256, 128, 2, stride=2)\n        self.dc8 = double_conv(256, 128)\n        self.up4 = nn.ConvTranspose2d(128, 64, 2, stride=2)\n        self.dc9 = double_conv(128, 64)\n        self.final = nn.Conv2d(64, n_classes, 1)\n\n    def forward(self, x):\n        x1 = self.dc1(x)\n        x2 = self.dc2(nn.MaxPool2d(2)(x1))\n        x3 = self.dc3(nn.MaxPool2d(2)(x2))\n        x4 = self.dc4(nn.MaxPool2d(2)(x3))\n        x5 = self.dc5(nn.MaxPool2d(2)(x4))\n\n        x = self.up1(x5)\n        x = self.dc6(torch.cat([x4, x], dim=1))\n        x = self.up2(x)\n        x = self.dc7(torch.cat([x3, x], dim=1))\n        x = self.up3(x)\n        x = self.dc8(torch.cat([x2, x], dim=1))\n        x = self.up4(x)\n        x = self.dc9(torch.cat([x1, x], dim=1))\n        return torch.sigmoid(self.final(x))\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.470621Z","iopub.execute_input":"2024-10-01T17:23:43.470887Z","iopub.status.idle":"2024-10-01T17:23:43.491877Z","shell.execute_reply.started":"2024-10-01T17:23:43.470861Z","shell.execute_reply":"2024-10-01T17:23:43.491237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ===============================\n# Cell 6: Define Loss Function\n# ===============================\ndef dice_loss(pred, target, smooth=1.0):\n    pred = pred.contiguous()\n    target = target.contiguous()\n    intersection = (pred * target).sum(dim=2).sum(dim=2)\n    loss = (1 - ((2. * intersection + smooth) / (pred.sum(dim=2).sum(dim=2) + target.sum(dim=2).sum(dim=2) + smooth)))\n    return loss.mean()\n\ndef calc_loss(pred, target, bce_weight=0.5):\n    bce = nn.BCELoss()(pred, target)\n    dice = dice_loss(pred, target)\n    loss = bce * bce_weight + dice * (1 - bce_weight)\n    return loss\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.493037Z","iopub.execute_input":"2024-10-01T17:23:43.493337Z","iopub.status.idle":"2024-10-01T17:23:43.508784Z","shell.execute_reply.started":"2024-10-01T17:23:43.49331Z","shell.execute_reply":"2024-10-01T17:23:43.508066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# =========================================\n# Cell 6: Training Function and Main\n# =========================================\ndef train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=50):\n    model.to(device)\n    train_losses = []\n    val_losses = []\n\n    for epoch in range(num_epochs):\n        model.train()\n        train_loss = 0.0\n\n        for images, masks in tqdm(train_loader, desc=f\"Epoch {epoch+1}/{num_epochs}\"):\n            images, masks = images.to(device), masks.to(device)\n            optimizer.zero_grad()\n            outputs = model(images)\n            loss = criterion(outputs, masks.unsqueeze(1))\n            loss.backward()\n            optimizer.step()\n            train_loss += loss.item() * images.size(0)\n\n        train_loss = train_loss / len(train_loader.dataset)\n        train_losses.append(train_loss)\n\n        print(f\"Epoch {epoch+1}/{num_epochs}, Train Loss: {train_loss:.4f}\")\n\n    return model, train_losses\n\n\ndef main():\n    data_dir = '/kaggle/input/gf-floodnet'\n    batch_size = 16\n    num_epochs = 50\n    learning_rate = 1e-4\n\n    train_loader, val_loader, _ = get_data_loaders(data_dir, batch_size=batch_size)\n    model = UNet(n_channels=5, n_classes=1)\n    criterion = nn.BCELoss()\n    optimizer = optim.Adam(model.parameters(), lr=learning_rate)\n\n    # Start training\n    model, train_losses = train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=num_epochs)\n    print(\"Training completed successfully.\")\n\n\n# Start the training process\nmain()\n","metadata":{"execution":{"iopub.status.busy":"2024-10-01T17:23:43.510273Z","iopub.execute_input":"2024-10-01T17:23:43.51057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# # =========================================\n# # Cell 7: Define the Main Function to Train the Model\n# # =========================================\n# def main():\n#     data_dir = '/kaggle/input/gf-floodnet'  # Set the path to your Kaggle input data directory\n#     checkpoint_dir = 'checkpoints'\n#     batch_size = 16\n#     num_epochs = 50\n#     learning_rate = 1e-4\n\n#     train_loader, val_loader, test_loader = get_data_loaders(data_dir, batch_size=batch_size)\n#     model = UNet(n_channels=5, n_classes=1)\n\n#     criterion = calc_loss\n#     optimizer = optim.Adam(model.parameters(), lr=learning_rate)\n\n#     trained_model, train_losses, val_losses = train_model(model, train_loader, val_loader, criterion, optimizer, num_epochs=num_epochs)\n#     print(\"Training completed successfully.\")\n\n# # Call main to start training\n# main()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Ensure GPU is used if available\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\nprint(f\"Using device: {device}\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}