{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":71549,"databundleVersionId":8561470,"sourceType":"competition"},{"sourceId":9426995,"sourceType":"datasetVersion","datasetId":5726454}],"dockerImageVersionId":30762,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Step 1: Load and Preprocess the Data\nFirst, we'll load the CSV file containing the labels and read the images from the directory. We'll process each image to:\n\nCollect up to 10 coordinate pairs (x, y) with their corresponding severity levels.\nPad with (0, 0) and severity 0 if fewer than 10 coordinates.\nNormalize the images.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\n# Load the labels CSV file\nlabels_csv_path = '/kaggle/input/rsna-consolidated-training-data/RSNA_consolidated_train_data.csv'\ndf = pd.read_csv(labels_csv_path)\n\n# Drop rows with missing values if any\ndf = df.dropna()\n\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:05.335084Z","iopub.execute_input":"2024-09-24T15:22:05.335402Z","iopub.status.idle":"2024-09-24T15:22:05.943297Z","shell.execute_reply.started":"2024-09-24T15:22:05.335362Z","shell.execute_reply":"2024-09-24T15:22:05.942347Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Create a Mapping for Severity Levels\nWe'll create a mapping to convert the string severity levels into integers:","metadata":{}},{"cell_type":"code","source":"severity_mapping = {\n    \"Normal/Mild\": 1,\n    \"Moderate\": 2,\n    \"Severe\": 3\n}\n\n\n# Map the severity strings to integers\ndf['severity'] = df['severity'].map(severity_mapping)\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:11.548189Z","iopub.execute_input":"2024-09-24T15:22:11.548577Z","iopub.status.idle":"2024-09-24T15:22:11.565222Z","shell.execute_reply.started":"2024-09-24T15:22:11.548540Z","shell.execute_reply":"2024-09-24T15:22:11.564299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.2 Aggregate Labels for Each Image\nWe'll group the labels by study_id, series_id, and instance_number (which uniquely identify each image).","metadata":{}},{"cell_type":"code","source":"# Group labels by image\ngrouped = df.groupby(['study_id', 'series_id', 'instance_number'])\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:12.177885Z","iopub.execute_input":"2024-09-24T15:22:12.178260Z","iopub.status.idle":"2024-09-24T15:22:12.183274Z","shell.execute_reply.started":"2024-09-24T15:22:12.178225Z","shell.execute_reply":"2024-09-24T15:22:12.182195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"grouped.head()","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:12.498247Z","iopub.execute_input":"2024-09-24T15:22:12.498626Z","iopub.status.idle":"2024-09-24T15:22:12.533608Z","shell.execute_reply.started":"2024-09-24T15:22:12.498587Z","shell.execute_reply":"2024-09-24T15:22:12.532778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.3 Prepare the Data List\nWe'll create a list where each element corresponds to an image and contains:\n\nThe image path.\nThe series description.\nA list of up to 10 (x, y) coordinate pairs.\nA list of up to 10 severity levels.","metadata":{}},{"cell_type":"code","source":"import os\n\nimage_dir = '/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/train_images'\n\ndata_list = []\n\nfor (study_id, series_id, instance_number), group in grouped:\n    image_path = os.path.join(image_dir, str(study_id), str(series_id), f\"{instance_number}.dcm\")\n    series_description = group['series_description'].iloc[0]\n\n    # Collect coordinates and severity levels\n    coordinates = group[['x', 'y']].values.tolist()\n    severities = group['severity'].tolist()\n\n    # Pad coordinates and severities to length 10\n    num_coords = len(coordinates)\n    if num_coords < 10:\n        padding = 10 - num_coords\n        coordinates.extend([[0, 0]] * padding)\n        severities.extend([0] * padding)\n    else:\n        # If more than 10 coordinates, truncate the list\n        coordinates = coordinates[:10]\n        severities = severities[:10]\n\n    data_list.append({\n        'image_path': image_path,\n        'series_description': series_description,\n        'coordinates': coordinates\n    })\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:13.125803Z","iopub.execute_input":"2024-09-24T15:22:13.126147Z","iopub.status.idle":"2024-09-24T15:22:26.683106Z","shell.execute_reply.started":"2024-09-24T15:22:13.126116Z","shell.execute_reply":"2024-09-24T15:22:26.682043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.4 Filter Series Descriptions\nSince you only want to include images with specific series descriptions, we'll filter them.","metadata":{}},{"cell_type":"code","source":"# Allowed series descriptions\nallowed_series_descriptions = [\"Sagittal T1\", \"Axial T2\", \"Sagittal T2/STIR SCS\"]\n\n# Filter data_list based on series_description\ndata_list = [item for item in data_list if item['series_description'] in allowed_series_descriptions]\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:26.684866Z","iopub.execute_input":"2024-09-24T15:22:26.685186Z","iopub.status.idle":"2024-09-24T15:22:26.699687Z","shell.execute_reply.started":"2024-09-24T15:22:26.685144Z","shell.execute_reply":"2024-09-24T15:22:26.698819Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"1.5 Normalize the Images\nWe'll normalize the images by dividing by 255 and then by the mean and standard deviation. We'll compute the mean and std from the dataset.\n\nFirst, compute mean and std (you can precompute this and save it).","metadata":{}},{"cell_type":"code","source":"import random\nimport pydicom\nimport numpy as np\n\ndef compute_mean_std(data_list, series_description_field='series_description', batch_size=100):\n    # Step 1: Group the data by series description\n    grouped_data = {}\n    for item in data_list:\n        series_description = item[series_description_field]\n        if series_description not in grouped_data:\n            grouped_data[series_description] = []\n        grouped_data[series_description].append(item)\n\n    # Step 2: Take 45% of data from each series description group, ensuring 15% of total from each group\n    sampled_data = []\n    for series_desc, group in grouped_data.items():\n        # Shuffle the group to get random sampling\n        random.shuffle(group)\n        # Get 15% of the total data length from each series description group\n        sample_size = int(len(data_list) * 0.15)\n        sampled_data.extend(group[:sample_size])\n\n    # Step 3: Compute mean and std incrementally in batches\n    pixel_sum = 0.0\n    pixel_sum_sq = 0.0\n    pixel_count = 0\n\n    for i in range(0, len(sampled_data), batch_size):\n        batch = sampled_data[i:i+batch_size]\n        for item in batch:\n            image_path = item['image_path']\n            # Load image and extract pixel values\n            dicom_data = pydicom.dcmread(image_path)\n            image = dicom_data.pixel_array.astype(np.float32) / 255.0\n            image_pixels = image.flatten()\n\n            # Update running sums for mean and standard deviation calculation\n            pixel_sum += np.sum(image_pixels)\n            pixel_sum_sq += np.sum(np.square(image_pixels))\n            pixel_count += len(image_pixels)\n\n    # Step 4: Compute mean and std\n    mean = pixel_sum / pixel_count\n    std = np.sqrt((pixel_sum_sq / pixel_count) - (mean ** 2))\n    \n    return mean, std\n\n\n# Example usage:\n# mean, std = compute_mean_std(data_list)\n# print(f\"Mean: {mean}, Std: {std}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:26.701128Z","iopub.execute_input":"2024-09-24T15:22:26.701840Z","iopub.status.idle":"2024-09-24T15:22:26.895663Z","shell.execute_reply.started":"2024-09-24T15:22:26.701794Z","shell.execute_reply":"2024-09-24T15:22:26.894702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# For the purpose of this example, let's assume some values\nmean = 0.9222420911188985  # Replace with actual mean\nstd = 1.720927786896051  # Replace with actual std\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:26.897712Z","iopub.execute_input":"2024-09-24T15:22:26.898052Z","iopub.status.idle":"2024-09-24T15:22:26.902486Z","shell.execute_reply.started":"2024-09-24T15:22:26.898020Z","shell.execute_reply":"2024-09-24T15:22:26.901557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step 2: Split the Data\nWe'll split the data into training and validation sets, ensuring equal representation of all severity levels.\n\n2.1 Create a DataFrame from data_list","metadata":{}},{"cell_type":"code","source":"data_df = pd.DataFrame(data_list)\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:26.903651Z","iopub.execute_input":"2024-09-24T15:22:26.903979Z","iopub.status.idle":"2024-09-24T15:22:26.936038Z","shell.execute_reply.started":"2024-09-24T15:22:26.903949Z","shell.execute_reply":"2024-09-24T15:22:26.935355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"2.2 Assign a Primary Severity Level to Each Image\nFor balancing, we'll assign a primary severity level to each image based on the maximum severity level in that image.","metadata":{}},{"cell_type":"markdown","source":"2.3 Stratified Split\nWe'll perform a stratified split based on the primary_severity to ensure equal representation.","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\ntrain_df, val_df = train_test_split(\n    data_df,\n    test_size=0.4,\n    random_state=42,\n)\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:22:26.937094Z","iopub.execute_input":"2024-09-24T15:22:26.937455Z","iopub.status.idle":"2024-09-24T15:22:27.573549Z","shell.execute_reply.started":"2024-09-24T15:22:26.937413Z","shell.execute_reply":"2024-09-24T15:22:27.572770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now both train_df and val_df have approximately equal representations of all severity levels.","metadata":{}},{"cell_type":"markdown","source":"Step 3: Create an Efficient Data Loader\nWe'll use TensorFlow's tf.data API to create an efficient data loader.\n\n3.1 Define the Dataset Class","metadata":{}},{"cell_type":"code","source":"import os\nimport cv2\nimport numpy as np\nimport tensorflow as tf\n\n# Example mapping for series descriptions to integer values\nseries_description_mapping = {\n    \"Sagittal T1\": 1,\n    \"Axial T2\": 2,\n    \"Sagittal T2/STIR SCS\": 3\n}\n\nclass LumbarSpineDataset(tf.keras.utils.Sequence):\n    def __init__(self, dataframe, batch_size=32, shuffle=True, image_size=(224, 224),\n                 mean=0.5, std=0.5, augment=False):\n        \"\"\"\n        Initializes the dataset object.\n\n        :param dataframe: Pandas dataframe with columns 'image_path', 'series_description', and 'coordinates'\n        :param batch_size: Number of samples per batch\n        :param shuffle: Boolean to determine whether to shuffle the data\n        :param image_size: Target size to resize images to\n        :param mean: Mean value for image normalization\n        :param std: Standard deviation for image normalization\n        :param augment: Boolean to apply augmentations\n        \"\"\"\n        self.dataframe = dataframe\n        self.batch_size = batch_size\n        self.shuffle = shuffle\n        self.image_size = image_size\n        self.mean = mean\n        self.std = std\n        self.augment = augment\n        self.normalizer = Normalise(mean, std)\n        self.rotator = Rotate_Flip()\n        self.indexes = np.arange(len(self.dataframe))\n        self.on_epoch_end()\n\n    def __len__(self):\n        \"\"\" Returns the number of batches per epoch \"\"\"\n        return len(self.dataframe) // self.batch_size\n\n    def on_epoch_end(self):\n        \"\"\" Updates indexes after each epoch \"\"\"\n        if self.shuffle:\n            np.random.shuffle(self.indexes)\n\n    def __getitem__(self, idx):\n        \"\"\" Generates a batch of data \"\"\"\n        batch_indexes = self.indexes[idx * self.batch_size:(idx + 1) * self.batch_size]\n        batch_data = self.dataframe.iloc[batch_indexes]\n\n        images, series_descriptions, coordinates = [], [], []\n\n        for _, row in batch_data.iterrows():\n            image_path = row['image_path']\n            series_description = row['series_description']\n            coords = row['coordinates']\n\n            # Load image and recalculate coordinates\n            image, recalculated_coords1 = self.load_image(image_path, coords)\n\n#             # Apply augmentations if needed\n#             if self.augment:\n#                 image, recalculated_coords2 = self.apply_augmentations(image, recalculated_coords1)\n            \n            # Ensure the image has a channel dimension (grayscale)\n            if len(image.shape) == 2:  # If image has no channel dimension\n                image = np.expand_dims(image, axis=-1)  # Add channel dimension\n\n            # Print image shape to debug\n#             print(f\"Loaded image shape: {image.shape}\")  # Debugging the image shape\n            \n            original_height, original_width = image.shape[:2]\n            recalculated_coords = []\n            \n            for x, y in recalculated_coords1:\n                recalculated_x = x / 224\n                recalculated_y = y / 224\n                recalculated_coords.append([recalculated_x, recalculated_y])\n                # Print image shape to debug\n#                 print(f\"recalculated_x, recalculated_y: {recalculated_x, recalculated_y}\")  # Debugging the image shape\n\n            # Convert the series_description string to an integer using the mapping\n            series_description_int = series_description_mapping[series_description]\n\n            images.append(image)\n            series_descriptions.append(series_description_int)  # Use integer label\n            coordinates.append(recalculated_coords)\n\n        \n        # Convert lists to NumPy arrays and then to Tensors\n        images = tf.convert_to_tensor(np.array(images), dtype=tf.float32)  # Ensure all images have the same shape\n        series_descriptions = tf.convert_to_tensor(np.array(series_descriptions, dtype=np.int32))\n        coordinates = tf.convert_to_tensor(np.array(coordinates), dtype=tf.float32)\n\n        return (images, series_descriptions), coordinates\n\n\n\n    def load_image(self, image_path, coordinates):\n        \"\"\" \n        Load and preprocess image, and recalculate coordinate positions based on the resized image.\n        :param image_path: Path to the DICOM image file.\n        :param coordinates: Original coordinates before resizing.\n        :return: Tuple of (preprocessed image, recalculated coordinates).\n        \"\"\"\n        # Load DICOM image\n        dicom_data = pydicom.dcmread(image_path)\n        image = dicom_data.pixel_array.astype(np.float32)\n\n        # Get original image dimensions\n        original_height, original_width = image.shape[:2]\n\n        # Resize the image\n        image_resized = cv2.resize(image, self.image_size)\n\n        # Normalize the resized image\n        image_normalized = self.normalizer(image_resized)\n\n        # Add a channel dimension to the image (for grayscale, channel=1)\n        image_normalized = np.expand_dims(image_normalized, axis=-1)\n\n        # Recalculate coordinates for the resized image\n        recalculated_coords = []\n        for x, y in coordinates:\n            recalculated_x = (x / original_width) * self.image_size[1]  # Scale x-coordinate\n            recalculated_y = (y / original_height) * self.image_size[0]  # Scale y-coordinate\n            recalculated_coords.append([recalculated_x, recalculated_y])\n\n        # Return preprocessed image and recalculated coordinates\n        return image_normalized, np.array(recalculated_coords, dtype=np.float32)\n\n\n\n    def normalize_coordinates(self, coordinates, image):\n        \"\"\" Normalize coordinates to be within [0, 1] \"\"\"\n        original_height, original_width = image.shape[:2]\n        recalculated_coords = []\n        for x, y in coordinates:\n            recalculated_x = x / original_width\n            recalculated_y = y / original_height\n            recalculated_coords.append([recalculated_x, recalculated_y])\n        return np.array(recalculated_coords, dtype=np.float32)\n\n    def apply_augmentations(self, image, coordinates):\n        \"\"\" Apply augmentations to the image and coordinates \"\"\"\n        image, coordinates = self.rotator(image, coordinates)\n        return image, coordinates\n\n# Augmentation Classes\n\nclass Rotate_Flip:\n    \"\"\" Randomly rotate image by 90 degrees and randomly flip image horizontally and vertically \"\"\"\n    def __call__(self, image, coordinates):\n        \"\"\"\n        Apply the same transformations to both the image and the coordinates.\n        :param image: The input image.\n        :param coordinates: The input coordinates (normalized).\n        :return: Transformed image and coordinates\n        \"\"\"\n        h, w = image.shape[:2]\n\n        # Apply horizontal flip\n        if np.random.randint(2):\n            image = cv2.flip(image, 1)  # Horizontal flip\n            coordinates[:, 0] = 1 - coordinates[:, 0]  # Flip x-coordinate\n\n        # Apply vertical flip\n        if np.random.randint(2):\n            image = cv2.flip(image, 0)  # Vertical flip\n            coordinates[:, 1] = 1 - coordinates[:, 1]  # Flip y-coordinate\n\n        # Random rotation by 0, 90, 180, or 270 degrees using OpenCV\n        angle = np.random.randint(0, 4) * 90\n        if angle > 0:\n            # Get the rotation matrix for the image\n            center = (w // 2, h // 2)\n            rotation_matrix = cv2.getRotationMatrix2D(center, angle, 1.0)\n\n            # Rotate the image\n            image = cv2.warpAffine(image, rotation_matrix, (w, h))\n\n            # Adjust the coordinates according to the rotation\n            coordinates = self.rotate_coordinates(coordinates, angle, h, w)\n\n        return image, coordinates\n\n    def rotate_coordinates(self, coordinates, angle, height, width):\n        \"\"\"\n        Rotate the coordinates based on the angle using OpenCV rotation matrix.\n        :param coordinates: Normalized coordinates (x, y).\n        :param angle: Angle of rotation.\n        :param height: Height of the image.\n        :param width: Width of the image.\n        :return: Rotated coordinates\n        \"\"\"\n        # Convert normalized coordinates to absolute pixel values\n        coords_abs = coordinates * np.array([width, height])\n\n        # Create a rotation matrix based on the angle\n        if angle == 90:\n            coords_rotated = np.array([[y, width - x] for x, y in coords_abs])  # (x, y) -> (y, w-x)\n        elif angle == 180:\n            coords_rotated = np.array([[width - x, height - y] for x, y in coords_abs])  # (x, y) -> (w-x, h-y)\n        elif angle == 270:\n            coords_rotated = np.array([[height - y, x] for x, y in coords_abs])  # (x, y) -> (h-y, x)\n        else:\n            coords_rotated = coords_abs  # No rotation\n\n        # Normalize the rotated coordinates\n        coords_rotated = coords_rotated / np.array([width, height])\n\n        return coords_rotated\n\nclass Normalise:\n    \"\"\" Normalize image pixels \"\"\"\n    def __init__(self, mean, std):\n        self.mean = mean\n        self.std = std\n\n    def __call__(self, image):\n        image = image / 255\n        image = (image - self.mean) / self.std\n        return image\n\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:28:58.471321Z","iopub.execute_input":"2024-09-24T15:28:58.471861Z","iopub.status.idle":"2024-09-24T15:28:58.504256Z","shell.execute_reply.started":"2024-09-24T15:28:58.471820Z","shell.execute_reply":"2024-09-24T15:28:58.503426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"3.2 Create the Datasets","metadata":{}},{"cell_type":"code","source":"# For the purpose of this example, let's assume some values\nmean = 0.9222420911188985  # Replace with actual mean\nstd = 1.720927786896051  # Replace with actual std\n\n# Define batch size, mean, and std for normalization\nbatch_size = 32\n\n# Create training dataset using the custom Sequence-based data loader\ntrain_dataset = LumbarSpineDataset(train_df, batch_size=batch_size, augment=True, mean=mean, std=std)\n\n# Create validation dataset using the custom Sequence-based data loader\nval_dataset = LumbarSpineDataset(val_df, batch_size=batch_size, augment=False, mean=mean, std=std)\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:28:59.485866Z","iopub.execute_input":"2024-09-24T15:28:59.486698Z","iopub.status.idle":"2024-09-24T15:28:59.492770Z","shell.execute_reply.started":"2024-09-24T15:28:59.486660Z","shell.execute_reply":"2024-09-24T15:28:59.491790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example batch usage in the LumbarSpineDataset\n\n# Get a batch of images and coordinates\n(images, series_descriptions), coordinates = train_dataset.__getitem__(0)\n\n# Check the shape of the images\nprint(f\"Batch image shape: {images.shape}\")  # Should be (batch_size, height, width, 1)\nprint(f\"Batch coordinates shape: {coordinates.shape}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:29:00.424622Z","iopub.execute_input":"2024-09-24T15:29:00.425014Z","iopub.status.idle":"2024-09-24T15:29:01.225170Z","shell.execute_reply.started":"2024-09-24T15:29:00.424977Z","shell.execute_reply":"2024-09-24T15:29:01.224215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport cv2\n\ndef denormalize_image(image, mean, std):\n    \"\"\"\n    Denormalizes the image by reversing the normalization steps applied during preprocessing.\n    \"\"\"\n    denormalized_image = (image * std) + mean  # Reverse the normalization\n    denormalized_image = np.clip(image * 255, 0, 255).astype(np.uint8)  # Rescale to 0-255\n    return denormalized_image\n\ndef denormalize_coordinates(normalized_coords, image_size):\n    \"\"\"\n    Denormalizes the coordinates by reversing the normalization applied during preprocessing.\n    \"\"\"\n    denormalized_coords = []\n    for coord in normalized_coords:\n        print(coord)\n        x, y = coord[:2]  # Fix: Take only the first two values (x, y), ignore other values if present.\n        # Skip (0, 0) coordinates (ignore all-zero coordinates)\n        if x == 0 and y == 0:\n            continue\n        denormalized_x = x * 224  # Multiply by image width\n        denormalized_y = y * 224  # Multiply by image height\n        denormalized_coords.append([denormalized_x, denormalized_y])\n    return denormalized_coords\n\ndef plot_image_with_coordinates(ax, image, coordinates, title=\"Image with Coordinates\"):\n    \"\"\"\n    Plots the denormalized image and overlays the denormalized coordinates on it within a specific axis.\n    \n    Args:\n    - ax: The matplotlib axis to plot on.\n    - image: The denormalized image to plot.\n    - coordinates: The denormalized (x, y) coordinates to overlay on the image.\n    - title: The title of the plot.\n    \"\"\"\n    image = cv2.resize(image, (224,224), interpolation=cv2.INTER_LINEAR)\n    ax.imshow(image, cmap='gray', aspect='auto')  # Set aspect='auto' to avoid stretching\n    coordinates = np.array(coordinates)\n    if len(coordinates) > 0:  # Check if there are valid coordinates\n        ax.scatter(coordinates[:, 0], coordinates[:, 1], c='r', marker='x')  # Red 'x' marks for coordinates\n    ax.set_title(title)\n    ax.set_xlabel('X')\n    ax.set_ylabel('Y')\n\ndef visualize_samples_in_grid(dataset, mean, std, image_size=(224, 224), num_samples=6, grid_size=(2, 3)):\n    \"\"\"\n    Visualizes multiple samples from the dataset in a grid along with their denormalized coordinates.\n    \n    Args:\n    - dataset: The TensorFlow dataset from which to pick samples.\n    - mean: The mean used during normalization of the images.\n    - std: The standard deviation used during normalization of the images.\n    - image_size: The fixed size to which the images were resized.\n    - num_samples: Number of samples to visualize.\n    - grid_size: Tuple indicating grid dimensions (rows, cols).\n    \"\"\"\n    dataset_iter = iter(dataset)\n    \n    # Create a grid of subplots\n    fig, axes = plt.subplots(grid_size[0], grid_size[1], figsize=(15, 10))\n    axes = axes.ravel()  # Flatten the grid to easily iterate over\n    \n    for i in range(num_samples):\n        (image_batch, _), coordinates_batch = next(dataset_iter)\n        \n        # Get the first image and coordinates from the batch\n        image = image_batch[0, :, :, 0]  # First image, grayscale (height, width)\n        coordinates = coordinates_batch[0]  # First set of normalized coordinates\n        \n        # Denormalize the image and coordinates for visualization\n        denormalized_image = denormalize_image(image, mean, std)\n        denormalized_coords = denormalize_coordinates(coordinates, image_size)\n        \n        # Plot the image and coordinates in the grid\n        plot_image_with_coordinates(axes[i], denormalized_image, denormalized_coords, title=f\"Sample {i + 1}\")\n    \n    # Adjust layout and show the plot\n    plt.tight_layout()\n    plt.show()\n\n# Visualize 6 samples from the dataset in a 2x3 grid\nvisualize_samples_in_grid(train_dataset, mean, std, image_size=(224, 224), num_samples=6, grid_size=(2, 3))\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:29:01.306372Z","iopub.execute_input":"2024-09-24T15:29:01.306687Z","iopub.status.idle":"2024-09-24T15:29:07.273332Z","shell.execute_reply.started":"2024-09-24T15:29:01.306656Z","shell.execute_reply":"2024-09-24T15:29:07.272454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"3.3 Map Series Description Labels to Integers\nEnsure that allowed_series_descriptions is defined and consistent.","metadata":{}},{"cell_type":"code","source":"allowed_series_descriptions = [\"Sagittal T1\", \"Axial T2\", \"Sagittal T2/STIR SCS\"]\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:29:07.274813Z","iopub.execute_input":"2024-09-24T15:29:07.275100Z","iopub.status.idle":"2024-09-24T15:29:07.279437Z","shell.execute_reply.started":"2024-09-24T15:29:07.275069Z","shell.execute_reply":"2024-09-24T15:29:07.278439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Step 4: Update the Model to Accommodate the New Data\nSince the model now receives images and series descriptions (as integer labels), and outputs coordinates and severity levels, we need to adjust the model accordingly.\n\n4.1 Update the Model Inputs","metadata":{}},{"cell_type":"markdown","source":"A common approach to predict uncertainty in regression tasks is to model the coordinates as probability distributions rather than point estimates. This can be done by predicting both the mean and the variance (standard deviation) of the coordinate predictions. You can then treat the coordinate predictions as samples from a Gaussian distribution, where the predicted mean represents the actual coordinate, and the variance represents the uncertainty.\n\nModifying the model for uncertainty prediction:\nCoordinate head outputs both mean and variance: Instead of directly predicting the coordinates, the model will predict two values per coordinate: the mean and the variance (uncertainty).\nLoss function: You will need to modify the loss function to account for the uncertainty by using a negative log-likelihood loss for the Gaussian distribution.\nReparameterization trick: To ensure the variance is positive, you can output log_variance and then exponentiate it.","metadata":{}},{"cell_type":"code","source":"import tensorflow as tf\n\n# Input for the image (224, 224, 1)\nimage_input = tf.keras.layers.Input(shape=(224, 224, 1), name='image_input')\n\n# ResNet50-based feature extractor (without pre-trained weights and no top layer)\nbase_model = tf.keras.applications.ResNet50(\n    include_top=False,\n    weights=None,  # No pre-trained weights since it's a grayscale image\n    input_shape=(224, 224, 1)\n)\n\n# Extract features from the image\nx = base_model(image_input)\n\n# Flatten the features from ResNet\nx = tf.keras.layers.Flatten()(x)\n\n# Input for the series description (integer-encoded, assume 3 possible descriptions)\nseries_description_input = tf.keras.layers.Input(shape=(), dtype='int32', name='series_description')\n\n# Embedding layer for the series description\nseries_embedding = tf.keras.layers.Embedding(\n    input_dim=3,  # Assuming 3 different series descriptions\n    output_dim=16\n)(series_description_input)\n\n# Concatenate the image features and the embedded series description\nconcat_features = tf.keras.layers.Concatenate()([x, series_embedding])\n\n# Coordinate prediction head (predicts x_mean, y_mean, log_x_variance, log_y_variance)\ncoord_pred = tf.keras.layers.Dense(512, activation='relu')(concat_features)\ncoord_pred = tf.keras.layers.Dense(256, activation='relu')(coord_pred)\ncoord_pred = tf.keras.layers.Dense(10 * 4)(coord_pred)  # Output: 10 (x_mean, y_mean, log_x_variance, log_y_variance)\n\n# Reshape the output to (10, 4)\ncoord_pred = tf.keras.layers.Reshape((10, 4))(coord_pred)\n\n# Split the predictions into mean and log variance\ncoord_mean = coord_pred[..., :2]  # (x_mean, y_mean)\ncoord_log_variance = coord_pred[..., 2:]  # (log_x_variance, log_y_variance)\n\n# Define the model with two inputs and two outputs\nmodel = tf.keras.models.Model(inputs=[image_input, series_description_input], outputs=[coord_mean, coord_log_variance])\n\n# Summary of the model\nmodel.summary()\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:29:14.729268Z","iopub.execute_input":"2024-09-24T15:29:14.729651Z","iopub.status.idle":"2024-09-24T15:29:15.417096Z","shell.execute_reply.started":"2024-09-24T15:29:14.729614Z","shell.execute_reply":"2024-09-24T15:29:15.416135Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"epochs = 10  # Set the number of epochs\noptimizer = tf.keras.optimizers.Adam(learning_rate=0.001)\n\n# Define custom loss functions\ndef gaussian_nll_loss(y_true_coords, y_pred):\n    print(y_pred)\n    coord_mean = y_pred[0]  # Predicted mean (x_mean, y_mean)\n    coord_log_variance = y_pred[1]  # Predicted log variance (log_x_variance, log_y_variance)\n\n    # Ensure that y_true_coords has shape [batch_size, 10, 2]\n    y_true_coords = tf.reshape(y_true_coords, [-1, 10, 2])\n\n    # Clamp log_variance to avoid extreme values\n    coord_log_variance = tf.clip_by_value(coord_log_variance, -10, 10)\n    \n    # Compute the variance\n    coord_variance = tf.exp(coord_log_variance)\n    \n    # Add a small epsilon to prevent division by zero\n    epsilon = 1e-6\n    coord_variance = coord_variance + epsilon\n\n    # Compute the loss with uncertainty\n    loss = tf.reduce_mean(\n        0.5 * tf.reduce_sum(((coord_mean - y_true_coords) ** 2 / coord_variance) + coord_log_variance, axis=-1)\n    )\n    \n    return loss\n\n\n# Severity loss (assuming you use sparse categorical cross-entropy)\nseverity_loss_fn = tf.keras.losses.SparseCategoricalCrossentropy(from_logits=True)\n\nimport tensorflow as tf\n\n# Directory where the model will be saved\ncheckpoint_dir = '/kaggle/working/best_model.keras'\n\n# Create a callback to save the best model based on validation loss\ncheckpoint_callback = tf.keras.callbacks.ModelCheckpoint(\n    filepath=checkpoint_dir,\n    monitor='val_loss',  # Save based on validation loss\n    save_best_only=True,  # Save only the best model\n    verbose=1\n)\n\n# Compile the model\nmodel.compile(\n    optimizer=optimizer,\n    loss=[gaussian_nll_loss, gaussian_nll_loss],  # Add other losses if required\n)\n\n# Train the model with the checkpoint callback\nhistory = model.fit(\n    train_dataset,\n    validation_data=val_dataset,\n    epochs=1000,\n    verbose=1,\n    callbacks=[checkpoint_callback]  # Include the checkpoint callback\n)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T15:29:16.103641Z","iopub.execute_input":"2024-09-24T15:29:16.104530Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import FileLink\n\n# Create a link to download the saved model\nFileLink('/kaggle/working/best_model.h5')\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T08:45:09.616508Z","iopub.status.idle":"2024-09-24T08:45:09.616908Z","shell.execute_reply.started":"2024-09-24T08:45:09.616720Z","shell.execute_reply":"2024-09-24T08:45:09.616741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport numpy as np\nimport tensorflow as tf\nimport cv2\nimport pydicom\n\ndef denormalize_image(image, mean, std):\n    \"\"\"\n    Denormalizes the image by reversing the normalization steps applied during preprocessing.\n    \"\"\"\n    denormalized_image = (image * std) + mean  # Reverse the normalization\n    denormalized_image = np.clip(denormalized_image * 255, 0, 255).astype(np.uint8)  # Rescale to 0-255\n    return denormalized_image\n\ndef denormalize_coordinates(normalized_coords, image_size):\n    \"\"\"\n    Denormalizes the coordinates by reversing the normalization applied during preprocessing.\n    \"\"\"\n    denormalized_coords = []\n    for coord in normalized_coords:\n        x, y = coord[:2]  # Take only x and y coordinates\n        if x == 0 and y == 0:\n            continue  # Ignore (0, 0) coordinates\n        denormalized_x = x * image_size[1]  # Multiply by image width\n        denormalized_y = y * image_size[0]  # Multiply by image height\n        denormalized_coords.append([denormalized_x, denormalized_y])\n    return denormalized_coords\n\ndef plot_image_with_predicted_coordinates(image, coordinates, title=\"Prediction Result\"):\n    \"\"\"\n    Plots the image and overlays the predicted coordinates on it.\n    \"\"\"\n    plt.imshow(image, cmap='gray')\n    \n    coordinates = np.array(coordinates)\n    if len(coordinates) > 0:\n        plt.scatter(coordinates[:, 0], coordinates[:, 1], c='r', marker='x')  # Red 'x' marks for coordinates\n    \n    plt.title(title)\n    plt.xlabel('X')\n    plt.ylabel('Y')\n    plt.show()\n\ndef preprocess_image(image, mean, std, image_size=(224, 224)):\n    \"\"\"\n    Resizes and normalizes the input image for the model.\n    \"\"\"\n    resized_image = cv2.resize(image, image_size)\n    normalized_image = (resized_image / 255.0 - mean) / std\n    preprocessed_image = np.expand_dims(normalized_image, axis=-1)  # Add channel dimension\n    return preprocessed_image\n\ndef preprocess_series_description(series_description, allowed_series_descriptions):\n    \"\"\"\n    Converts the series description into its corresponding integer label.\n    \"\"\"\n    return allowed_series_descriptions.index(series_description)\n\ndef predict_and_display(model, image, series_description, allowed_series_descriptions, mean, std, image_size=(224, 224)):\n    \"\"\"\n    Preprocesses the input image and series description, makes a prediction using the model, and displays the result.\n    \n    Args:\n    - model: The trained TensorFlow model.\n    - image: The original image to make a prediction on.\n    - series_description: The series description for the image.\n    - allowed_series_descriptions: The list of all possible series descriptions (used to convert to an integer).\n    - mean: The mean used during normalization.\n    - std: The standard deviation used during normalization.\n    - image_size: The size to which the image was resized for the model.\n    \"\"\"\n    # Preprocess the image for prediction\n    preprocessed_image = preprocess_image(image, mean, std, image_size)\n    \n    # Preprocess the series description (convert to integer label)\n    preprocessed_series_description = preprocess_series_description(series_description, allowed_series_descriptions)\n    \n    # Expand dimensions to simulate batch size 1\n    preprocessed_image = np.expand_dims(preprocessed_image, axis=0)  # Shape: (1, 224, 224, 1)\n    preprocessed_series_description = np.expand_dims(preprocessed_series_description, axis=0)  # Shape: (1,)\n    \n    # Make predictions (assume the model outputs coordinates as the second output)\n    predicted_coordinates = model.predict([preprocessed_image, preprocessed_series_description])[1][0]  # Get the first prediction\n    \n    # Denormalize the image and coordinates\n    denormalized_image = denormalize_image(image, mean, std)\n    denormalized_coords = denormalize_coordinates(predicted_coordinates, image_size)\n    \n    # Display the image with predicted coordinates\n    plot_image_with_predicted_coordinates(denormalized_image, denormalized_coords, title=\"Predicted Coordinates\")\n\n# Example usage\n# Load the model (replace 'model_path' with the actual path to your saved model)\n# Define the custom objects\ncustom_objects = {\n    'MultiTaskLumbarModelWithUncertainty': MultiTaskLumbarModelWithUncertainty,\n    'gaussian_nll_loss': gaussian_nll_loss  # If you are using a custom loss function\n}\n\n# Load the model with custom objects\nmodel = tf.keras.models.load_model('/kaggle/working/best_model.keras', custom_objects=custom_objects)\n\n# Load a test image (replace 'test_image_path' with the actual path)\ndicom_data = pydicom.dcmread('/kaggle/input/rsna-2024-lumbar-spine-degenerative-classification/test_images/44036939/2828203845/1.dcm')\ntest_image = dicom_data.pixel_array.astype(np.float32)\n\n# Series description of the test image (you should replace this with the actual description)\nseries_description = \"sAxial T2\"  # Replace this with the actual series description\n\n# List of all possible series descriptions (this should be defined according to your training data)\nallowed_series_descriptions = [\"Sagittal T1\", \"Axial T2\", \"Sagittal T2/STIR SCS\"]  # Replace with actual descriptions\n\n# Make prediction and display results\npredict_and_display(model, test_image, series_description, allowed_series_descriptions, mean, std, image_size=(224, 224))\n","metadata":{"execution":{"iopub.status.busy":"2024-09-24T11:10:02.413767Z","iopub.execute_input":"2024-09-24T11:10:02.414386Z","iopub.status.idle":"2024-09-24T11:10:02.634548Z","shell.execute_reply.started":"2024-09-24T11:10:02.414344Z","shell.execute_reply":"2024-09-24T11:10:02.633142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}