{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":36363,"databundleVersionId":4050810,"sourceType":"competition"},{"sourceId":4264054,"sourceType":"datasetVersion","datasetId":2406209}],"dockerImageVersionId":31260,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\nENHANCED STEP 1: ADVANCED MULTI-WINDOW DATASET CREATION\n========================================================\nProduction-quality preprocessing with state-of-the-art techniques\n\nNEW FEATURES:\n✅ Multi-window preprocessing (bone + soft tissue + wide)\n✅ CLAHE (Contrast Limited Adaptive Histogram Equalization)\n✅ Advanced quality control with detailed metrics\n✅ Robust error handling with automatic fallback\n✅ Cross-validation split preparation\n✅ Comprehensive data quality reports\n✅ Memory-efficient processing\n✅ Progress checkpointing & resume capability\n\nExpected output: High-quality 3-channel volumes ready for training\nRun time: ~15-25 minutes for 300 patients\n\"\"\"\n\nimport os\nimport gc\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nimport pydicom\nfrom glob import glob\nfrom scipy.ndimage import zoom, gaussian_filter\nfrom skimage import exposure\nfrom tqdm import tqdm\n\nprint(\"=\"*80)\nprint(\"🚀 ENHANCED STEP 1: ADVANCED DATASET CREATION v3.0\")\nprint(\"=\"*80)\n\n# ============================================================================\n# ENHANCED CONFIGURATION\n# ============================================================================\n\nCONFIG = {\n # Dataset selection\n 'num_fracture_patients': 150,\n 'num_normal_patients': 150,\n 'random_seed': 42,\n \n # Paths\n 'train_csv_path': '/kaggle/input/rsna-2022-cervical-spine-fracture-detection/train.csv',\n 'train_images_root': '/kaggle/input/rsna-2022-cervical-spine-fracture-detection/train_images',\n 'output_dir': '/kaggle/working/enhanced_dataset_v3',\n \n # Multi-window preprocessing (KEY IMPROVEMENT!)\n 'use_multi_window': True, # 3-channel output\n 'windows': {\n 'bone': {'center': 400, 'width': 1800}, # Bone structure\n 'soft_tissue': {'center': 40, 'width': 400}, # Soft tissue\n 'wide': {'center': 400, 'width': 4000} # Overall context\n },\n \n # CLAHE enhancement (KEY IMPROVEMENT!)\n 'use_clahe': True,\n 'clahe_clip_limit': 0.01, # Prevents over-enhancement\n \n # Resolution\n 'target_shape': (64, 224, 224), # (D, H, W)\n 'target_spacing': (2.0, 1.25, 1.25), # mm\n \n # Quality control\n 'min_slices': 15,\n 'max_slices': 400,\n 'min_hu': -2000,\n 'max_hu': 4000,\n 'verify_integrity': True,\n \n # Cross-validation preparation (KEY IMPROVEMENT!)\n 'create_cv_splits': True,\n 'n_folds': 5,\n \n # Advanced options\n 'enable_checkpointing': True,\n 'checkpoint_interval': 10,\n 'save_quality_report': True,\n}\n\nprint(f\"\\n📋 Configuration:\")\nprint(f\" • Total patients: {CONFIG['num_fracture_patients'] + CONFIG['num_normal_patients']}\")\nprint(f\" • Multi-window: {CONFIG['use_multi_window']} (3-channel output)\")\nprint(f\" • CLAHE enhancement: {CONFIG['use_clahe']}\")\nprint(f\" • Target shape: {CONFIG['target_shape']}\")\nprint(f\" • CV folds: {CONFIG['n_folds']}\")\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\n# ============================================================================\n# ADVANCED PREPROCESSING FUNCTIONS\n# ============================================================================\n\ndef load_dicom_series_robust(patient_folder):\n \"\"\"\n Robust DICOM loading with comprehensive validation\n \"\"\"\n dicom_files = sorted(glob(os.path.join(patient_folder, \"*.dcm\")))\n \n if len(dicom_files) == 0:\n raise ValueError(\"No DICOM files found\")\n \n # Load all slices\n slices = []\n failed = 0\n \n for dcm_file in dicom_files:\n try:\n ds = pydicom.dcmread(dcm_file)\n \n # Basic validation\n if not hasattr(ds, 'ImagePositionPatient'):\n failed += 1\n continue\n if not hasattr(ds, 'PixelSpacing'):\n failed += 1\n continue\n if ds.pixel_array.size == 0:\n failed += 1\n continue\n \n slices.append(ds)\n except Exception:\n failed += 1\n continue\n \n if len(slices) == 0:\n raise ValueError(f\"No valid slices (failed: {failed})\")\n \n # Sort by position\n slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n \n # Stack volume\n volume = np.stack([s.pixel_array for s in slices])\n metadata = slices[0]\n \n # Calculate spacing\n try:\n slice_thickness = float(metadata.SliceThickness)\n except:\n if len(slices) > 1:\n slice_thickness = abs(\n float(slices[1].ImagePositionPatient[2]) - \n float(slices[0].ImagePositionPatient[2])\n )\n else:\n slice_thickness = 1.0\n \n pixel_spacing = [float(x) for x in metadata.PixelSpacing]\n current_spacing = (slice_thickness, pixel_spacing[0], pixel_spacing[1])\n \n return volume, current_spacing, metadata\n\n\ndef apply_hu_conversion(volume, metadata):\n \"\"\"Convert raw pixel values to Hounsfield Units\"\"\"\n try:\n intercept = float(metadata.RescaleIntercept)\n slope = float(metadata.RescaleSlope)\n except:\n intercept = 0.0\n slope = 1.0\n \n volume_hu = volume.astype(np.float32) * slope + intercept\n return volume_hu\n\n\ndef apply_windowing(volume_hu, center, width):\n \"\"\"\n Apply CT windowing to emphasize specific tissue types\n \"\"\"\n lower = center - width // 2\n upper = center + width // 2\n \n volume_windowed = np.clip(volume_hu, lower, upper)\n volume_normalized = (volume_windowed - lower) / (upper - lower)\n \n return volume_normalized.astype(np.float32)\n\n\ndef apply_clahe_3d(volume, clip_limit=0.01):\n \"\"\"\n Apply CLAHE (Contrast Limited Adaptive Histogram Equalization) in 3D\n KEY IMPROVEMENT: Enhances local contrast while preserving overall structure\n \"\"\"\n enhanced_volume = np.zeros_like(volume)\n \n # Apply CLAHE slice-by-slice\n for i in range(volume.shape[0]):\n slice_2d = volume[i]\n \n # CLAHE expects uint8, so we need to convert\n slice_uint8 = (slice_2d * 255).astype(np.uint8)\n \n # Apply adaptive histogram equalization\n slice_enhanced = exposure.equalize_adapthist(\n slice_uint8,\n clip_limit=clip_limit,\n nbins=256\n )\n \n enhanced_volume[i] = slice_enhanced\n \n return enhanced_volume.astype(np.float32)\n\n\ndef create_multi_window_volume(volume_hu, config):\n \"\"\"\n Create multi-channel volume with different windowing\n KEY IMPROVEMENT: Captures complementary information\n \n Returns: (D, H, W, 3) volume with bone/soft_tissue/wide channels\n \"\"\"\n channels = []\n \n for window_name, window_params in config['windows'].items():\n # Apply windowing\n windowed = apply_windowing(\n volume_hu,\n center=window_params['center'],\n width=window_params['width']\n )\n \n # Apply CLAHE if enabled\n if config['use_clahe']:\n windowed = apply_clahe_3d(windowed, config['clahe_clip_limit'])\n \n channels.append(windowed)\n \n # Stack channels: (D, H, W, 3)\n multi_channel = np.stack(channels, axis=-1)\n \n return multi_channel\n\n\ndef resample_volume(volume, current_spacing, target_spacing):\n \"\"\"\n Resample volume to target spacing\n Handles multi-channel volumes\n \"\"\"\n if volume.ndim == 4: # Multi-channel\n # Calculate resize factor\n resize_factor = np.array(current_spacing) / np.array(target_spacing)\n resize_factor = np.append(resize_factor, 1.0) # Don't resize channel dim\n \n # Resample\n resampled = zoom(volume, resize_factor, order=1)\n else: # Single channel\n resize_factor = np.array(current_spacing) / np.array(target_spacing)\n resampled = zoom(volume, resize_factor, order=1)\n \n return resampled.astype(np.float32)\n\n\ndef crop_or_pad_to_shape(volume, target_shape):\n \"\"\"\n Intelligently crop or pad volume to target shape\n Handles multi-channel volumes\n \"\"\"\n is_multichannel = volume.ndim == 4\n \n if is_multichannel:\n current_shape = volume.shape[:3]\n n_channels = volume.shape[3]\n output = np.zeros((*target_shape, n_channels), dtype=volume.dtype)\n else:\n current_shape = volume.shape\n output = np.zeros(target_shape, dtype=volume.dtype)\n \n target_shape_arr = np.array(target_shape)\n current_shape_arr = np.array(current_shape)\n \n # Calculate slice indices\n slices_vol = []\n slices_out = []\n \n for i in range(3):\n if current_shape_arr[i] >= target_shape_arr[i]:\n # Crop - focus on cervical region for depth\n if i == 0: # Depth\n # Cervical spine typically in upper 20-40% of full spine scan\n start = int(current_shape_arr[i] * 0.15)\n start = max(0, min(start, current_shape_arr[i] - target_shape_arr[i]))\n else: # Height/Width - center crop\n start = (current_shape_arr[i] - target_shape_arr[i]) // 2\n \n slices_vol.append(slice(start, start + target_shape_arr[i]))\n slices_out.append(slice(0, target_shape_arr[i]))\n else:\n # Pad - center the content\n start = (target_shape_arr[i] - current_shape_arr[i]) // 2\n slices_vol.append(slice(0, current_shape_arr[i]))\n slices_out.append(slice(start, start + current_shape_arr[i]))\n \n # Copy data\n if is_multichannel:\n output[slices_out[0], slices_out[1], slices_out[2], :] = \\\n volume[slices_vol[0], slices_vol[1], slices_vol[2], :]\n else:\n output[slices_out[0], slices_out[1], slices_out[2]] = \\\n volume[slices_vol[0], slices_vol[1], slices_vol[2]]\n \n return output\n\n\ndef compute_quality_metrics(volume_hu, volume_final):\n \"\"\"\n Compute comprehensive quality metrics\n \"\"\"\n metrics = {\n # HU statistics\n 'hu_min': float(volume_hu.min()),\n 'hu_max': float(volume_hu.max()),\n 'hu_mean': float(volume_hu.mean()),\n 'hu_std': float(volume_hu.std()),\n \n # Shape info\n 'original_shape': str(volume_hu.shape),\n 'final_shape': str(volume_final.shape[:3] if volume_final.ndim == 4 else volume_final.shape),\n 'num_slices': int(volume_hu.shape[0]),\n \n # Intensity distribution (normalized)\n 'intensity_mean': float(volume_final.mean()),\n 'intensity_std': float(volume_final.std()),\n \n # Quality indicators\n 'empty_slices': int(np.sum(volume_hu.sum(axis=(1,2)) == 0)),\n 'near_empty_slices': int(np.sum(volume_hu.sum(axis=(1,2)) < 100)),\n }\n \n # Multi-channel specific metrics\n if volume_final.ndim == 4:\n metrics['n_channels'] = volume_final.shape[3]\n for ch in range(volume_final.shape[3]):\n channel_names = ['bone', 'soft_tissue', 'wide']\n metrics[f'{channel_names[ch]}_mean'] = float(volume_final[..., ch].mean())\n metrics[f'{channel_names[ch]}_std'] = float(volume_final[..., ch].std())\n \n return metrics\n\n\ndef preprocess_patient_advanced(patient_folder, config):\n \"\"\"\n MAIN PREPROCESSING PIPELINE with all enhancements\n \"\"\"\n # Load DICOM\n volume, current_spacing, metadata = load_dicom_series_robust(patient_folder)\n \n # Convert to HU\n volume_hu = apply_hu_conversion(volume, metadata)\n \n # Quality check\n if volume.shape[0] < config['min_slices']:\n raise ValueError(f\"Too few slices: {volume.shape[0]}\")\n if volume.shape[0] > config['max_slices']:\n raise ValueError(f\"Too many slices: {volume.shape[0]}\")\n \n # Create multi-window or single-window volume\n if config['use_multi_window']:\n volume_processed = create_multi_window_volume(volume_hu, config)\n else:\n # Single window (bone) with optional CLAHE\n volume_processed = apply_windowing(volume_hu, 400, 1800)\n if config['use_clahe']:\n volume_processed = apply_clahe_3d(volume_processed, config['clahe_clip_limit'])\n \n # Resample to target spacing\n volume_resampled = resample_volume(\n volume_processed,\n current_spacing,\n config['target_spacing']\n )\n \n # Crop/pad to target shape\n volume_final = crop_or_pad_to_shape(\n volume_resampled,\n config['target_shape']\n )\n \n # Compute quality metrics\n metrics = compute_quality_metrics(volume_hu, volume_final)\n \n return volume_final, metrics\n\n\n# ============================================================================\n# CHECKPOINT MANAGER\n# ============================================================================\n\nclass CheckpointManager:\n \"\"\"Manages processing checkpoints for resume capability\"\"\"\n \n def __init__(self, output_dir, enable=True):\n self.enable = enable\n self.checkpoint_file = os.path.join(output_dir, 'checkpoint.json')\n self.processed_patients = set()\n self.load()\n \n def load(self):\n if self.enable and os.path.exists(self.checkpoint_file):\n try:\n with open(self.checkpoint_file, 'r') as f:\n data = json.load(f)\n self.processed_patients = set(data.get('processed', []))\n if len(self.processed_patients) > 0:\n print(f\" 📂 Resuming: {len(self.processed_patients)} already processed\")\n except:\n self.processed_patients = set()\n \n def save(self, patient_id):\n if self.enable:\n self.processed_patients.add(str(patient_id))\n try:\n with open(self.checkpoint_file, 'w') as f:\n json.dump({'processed': list(self.processed_patients)}, f)\n except:\n pass\n \n def is_processed(self, patient_id):\n return str(patient_id) in self.processed_patients\n\n\n# ============================================================================\n# CROSS-VALIDATION SPLIT CREATOR\n# ============================================================================\n\ndef create_cv_splits(metadata_df, n_folds=5, random_seed=42):\n \"\"\"\n Create stratified K-fold cross-validation splits\n KEY IMPROVEMENT: Better evaluation than single train/val split\n \"\"\"\n from sklearn.model_selection import StratifiedKFold\n \n skf = StratifiedKFold(n_splits=n_folds, shuffle=True, random_state=random_seed)\n \n splits = []\n for fold_idx, (train_idx, val_idx) in enumerate(skf.split(\n metadata_df.index, metadata_df['has_fracture']\n )):\n train_patients = metadata_df.iloc[train_idx]['patient_id'].tolist()\n val_patients = metadata_df.iloc[val_idx]['patient_id'].tolist()\n \n splits.append({\n 'fold': fold_idx,\n 'train': train_patients,\n 'val': val_patients,\n 'n_train': len(train_patients),\n 'n_val': len(val_patients),\n 'train_fracture': int(metadata_df.iloc[train_idx]['has_fracture'].sum()),\n 'val_fracture': int(metadata_df.iloc[val_idx]['has_fracture'].sum()),\n })\n \n return splits\n\n\n# ============================================================================\n# MAIN PROCESSING FUNCTION\n# ============================================================================\n\ndef create_enhanced_dataset():\n \"\"\"Main processing pipeline\"\"\"\n \n print(f\"\\n{'='*80}\")\n print(\"🔧 INITIALIZING ENHANCED PIPELINE\")\n print(f\"{'='*80}\")\n \n # Create directories\n volumes_dir = os.path.join(CONFIG['output_dir'], 'volumes')\n os.makedirs(volumes_dir, exist_ok=True)\n \n # Initialize checkpoint manager\n checkpoint_mgr = CheckpointManager(CONFIG['output_dir'], CONFIG['enable_checkpointing'])\n \n # ========================================================================\n # PATIENT SELECTION\n # ========================================================================\n \n print(f\"\\n📊 STEP 1: Patient Selection\")\n print(f\"{'='*80}\")\n \n try:\n train_df = pd.read_csv(CONFIG['train_csv_path'])\n print(f\" ✓ Loaded train.csv: {len(train_df)} rows\")\n except Exception as e:\n print(f\" ✗ ERROR: {e}\")\n return None\n \n # Group by patient\n patient_fracture_status = train_df.groupby('StudyInstanceUID')['patient_overall'].first()\n fracture_patients = patient_fracture_status[patient_fracture_status == 1].index.tolist()\n normal_patients = patient_fracture_status[patient_fracture_status == 0].index.tolist()\n \n print(f\" ✓ Available:\")\n print(f\" - Fracture: {len(fracture_patients)}\")\n print(f\" - Normal: {len(normal_patients)}\")\n \n # Random selection\n np.random.seed(CONFIG['random_seed'])\n selected_fracture = np.random.choice(\n fracture_patients,\n size=min(CONFIG['num_fracture_patients'], len(fracture_patients)),\n replace=False\n )\n selected_normal = np.random.choice(\n normal_patients,\n size=min(CONFIG['num_normal_patients'], len(normal_patients)),\n replace=False\n )\n \n selected_patients = list(selected_fracture) + list(selected_normal)\n \n # Filter already processed\n remaining_patients = [p for p in selected_patients \n if not checkpoint_mgr.is_processed(p)]\n \n print(f\"\\n ✓ Selection:\")\n print(f\" - Total: {len(selected_patients)}\")\n print(f\" - Already processed: {len(selected_patients) - len(remaining_patients)}\")\n print(f\" - Remaining: {len(remaining_patients)}\")\n \n if len(remaining_patients) == 0:\n print(f\"\\n ⚠️ All patients processed!\")\n metadata_path = os.path.join(CONFIG['output_dir'], 'metadata.csv')\n if os.path.exists(metadata_path):\n return pd.read_csv(metadata_path)\n return None\n \n # ========================================================================\n # PREPROCESSING\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(f\"⚙️ STEP 2: Advanced Preprocessing\")\n print(f\"{'='*80}\\n\")\n \n metadata_list = []\n successful = 0\n failed = 0\n failed_patients = []\n \n for patient_id in tqdm(remaining_patients, desc=\"Processing\"):\n patient_folder = os.path.join(CONFIG['train_images_root'], str(patient_id))\n \n if not os.path.exists(patient_folder):\n failed += 1\n failed_patients.append({\n 'patient_id': str(patient_id),\n 'reason': 'folder_not_found'\n })\n continue\n \n try:\n # Preprocess with all enhancements\n volume, metrics = preprocess_patient_advanced(patient_folder, CONFIG)\n \n # Save volume\n save_path = os.path.join(volumes_dir, f\"{patient_id}.npy\")\n np.save(save_path, volume)\n \n # Verify if enabled\n if CONFIG['verify_integrity']:\n loaded = np.load(save_path)\n if not np.array_equal(volume, loaded):\n raise ValueError(\"Integrity check failed\")\n \n # Get labels\n patient_labels = train_df[train_df['StudyInstanceUID'] == patient_id].iloc[0]\n \n # Compile metadata\n metadata_entry = {\n 'patient_id': str(patient_id),\n 'has_fracture': int(patient_labels['patient_overall']),\n 'c1': int(patient_labels['C1']),\n 'c2': int(patient_labels['C2']),\n 'c3': int(patient_labels['C3']),\n 'c4': int(patient_labels['C4']),\n 'c5': int(patient_labels['C5']),\n 'c6': int(patient_labels['C6']),\n 'c7': int(patient_labels['C7']),\n 'file_size_mb': os.path.getsize(save_path) / (1024**2),\n **metrics\n }\n \n metadata_list.append(metadata_entry)\n successful += 1\n \n # Checkpoint\n if successful % CONFIG['checkpoint_interval'] == 0:\n checkpoint_mgr.save(patient_id)\n \n # Memory cleanup\n del volume\n if successful % 10 == 0:\n gc.collect()\n \n except Exception as e:\n failed += 1\n failed_patients.append({\n 'patient_id': str(patient_id),\n 'reason': str(e)[:200]\n })\n continue\n \n # ========================================================================\n # SAVE METADATA & CREATE CV SPLITS\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(f\"💾 STEP 3: Saving Metadata & Creating CV Splits\")\n print(f\"{'='*80}\")\n \n if len(metadata_list) == 0:\n print(f\"\\n ✗ No successful preprocessing!\")\n return None\n \n # Create metadata DataFrame\n metadata_df = pd.DataFrame(metadata_list)\n metadata_path = os.path.join(CONFIG['output_dir'], 'metadata.csv')\n \n # Merge with existing if resuming\n if os.path.exists(metadata_path):\n try:\n existing_df = pd.read_csv(metadata_path)\n metadata_df = pd.concat([existing_df, metadata_df], ignore_index=True)\n metadata_df = metadata_df.drop_duplicates(subset=['patient_id'], keep='last')\n except:\n pass\n \n metadata_df.to_csv(metadata_path, index=False)\n print(f\"\\n ✓ Saved metadata: {metadata_path}\")\n \n # Create cross-validation splits\n if CONFIG['create_cv_splits']:\n cv_splits = create_cv_splits(metadata_df, CONFIG['n_folds'], CONFIG['random_seed'])\n cv_path = os.path.join(CONFIG['output_dir'], 'cv_splits.json')\n \n with open(cv_path, 'w') as f:\n json.dump(cv_splits, f, indent=2)\n \n print(f\" ✓ Created {CONFIG['n_folds']}-fold CV splits: {cv_path}\")\n \n # Print split summary\n print(f\"\\n CV Split Summary:\")\n for split in cv_splits:\n print(f\" Fold {split['fold']}: \"\n f\"Train={split['n_train']} ({split['train_fracture']} frac), \"\n f\"Val={split['n_val']} ({split['val_fracture']} frac)\")\n \n # Save config\n config_path = os.path.join(CONFIG['output_dir'], 'config.json')\n with open(config_path, 'w') as f:\n json.dump(CONFIG, f, indent=2)\n \n # Save failed patients\n if failed_patients:\n failed_path = os.path.join(CONFIG['output_dir'], 'failed_patients.json')\n with open(failed_path, 'w') as f:\n json.dump(failed_patients, f, indent=2)\n \n # ========================================================================\n # QUALITY REPORT\n # ========================================================================\n \n if CONFIG['save_quality_report']:\n print(f\"\\n{'='*80}\")\n print(\"📊 STEP 4: Quality Report\")\n print(f\"{'='*80}\")\n \n report = {\n 'total_patients': len(metadata_df),\n 'fracture_patients': int(metadata_df['has_fracture'].sum()),\n 'normal_patients': int((1 - metadata_df['has_fracture']).sum()),\n 'success_rate': f\"{successful/(successful+failed)*100:.1f}%\",\n 'avg_slices': float(metadata_df['num_slices'].mean()),\n 'avg_file_size_mb': float(metadata_df['file_size_mb'].mean()),\n 'total_size_mb': float(metadata_df['file_size_mb'].sum()),\n 'hu_range': f\"[{metadata_df['hu_min'].min():.0f}, {metadata_df['hu_max'].max():.0f}]\",\n 'preprocessing': {\n 'multi_window': CONFIG['use_multi_window'],\n 'clahe': CONFIG['use_clahe'],\n 'target_shape': CONFIG['target_shape'],\n }\n }\n \n if CONFIG['use_multi_window']:\n report['channel_statistics'] = {\n 'bone_mean': float(metadata_df['bone_mean'].mean()),\n 'soft_tissue_mean': float(metadata_df['soft_tissue_mean'].mean()),\n 'wide_mean': float(metadata_df['wide_mean'].mean()),\n }\n \n report_path = os.path.join(CONFIG['output_dir'], 'quality_report.json')\n with open(report_path, 'w') as f:\n json.dump(report, f, indent=2)\n \n print(f\"\\n Quality Metrics:\")\n print(f\" • Success rate: {report['success_rate']}\")\n print(f\" • Avg slices: {report['avg_slices']:.1f}\")\n print(f\" • Avg size: {report['avg_file_size_mb']:.2f} MB\")\n print(f\" • Total size: {report['total_size_mb']:.2f} MB\")\n print(f\" • HU range: {report['hu_range']}\")\n \n if CONFIG['use_multi_window']:\n print(f\"\\n Channel Statistics:\")\n print(f\" • Bone: {report['channel_statistics']['bone_mean']:.4f}\")\n print(f\" • Soft tissue: {report['channel_statistics']['soft_tissue_mean']:.4f}\")\n print(f\" • Wide: {report['channel_statistics']['wide_mean']:.4f}\")\n \n # ========================================================================\n # SUMMARY\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(\"📈 FINAL SUMMARY\")\n print(f\"{'='*80}\")\n \n print(f\"\\n✅ Results:\")\n print(f\" • Successful: {successful}\")\n print(f\" • Failed: {failed}\")\n print(f\" • Dataset size: {len(metadata_df)}\")\n print(f\" • Fracture/Normal: {metadata_df['has_fracture'].sum()}/{(1-metadata_df['has_fracture']).sum()}\")\n \n print(f\"\\n📁 Output:\")\n print(f\" • Location: {CONFIG['output_dir']}/\")\n print(f\" • Metadata: metadata.csv\")\n print(f\" • CV splits: cv_splits.json\")\n print(f\" • Quality report: quality_report.json\")\n print(f\" • Volumes: volumes/*.npy\")\n \n print(f\"\\n🎯 Enhancements Applied:\")\n print(f\" ✅ Multi-window preprocessing (3-channel)\")\n print(f\" ✅ CLAHE contrast enhancement\")\n print(f\" ✅ Robust error handling\")\n print(f\" ✅ Cross-validation splits ready\")\n print(f\" ✅ Comprehensive quality metrics\")\n \n return metadata_df\n\n\n# ============================================================================\n# EXECUTION\n# ============================================================================\n\nif __name__ == \"__main__\":\n print(\"\\n\" + \"=\"*80)\n print(\"🚀 STARTING ENHANCED DATASET CREATION\")\n print(\"=\"*80)\n \n try:\n metadata_df = create_enhanced_dataset()\n \n if metadata_df is not None and len(metadata_df) > 0:\n print(\"\\n\" + \"=\"*80)\n print(\"✅ SUCCESS! Enhanced dataset created\")\n print(\"=\"*80)\n \n print(f\"\\n📊 Sample metadata:\")\n print(metadata_df[['patient_id', 'has_fracture', 'num_slices', 'file_size_mb']].head())\n \n print(f\"\\n🎯 Next Step: Enhanced Training with Multi-Channel Input\")\n print(f\" Run ENHANCED STEP 2 for training\")\n \n else:\n print(\"\\n\" + \"=\"*80)\n print(\"⚠️ PROCESSING INCOMPLETE\")\n print(\"=\"*80)\n print(\"\\nPlease check:\")\n print(\" 1. Input paths are correct\")\n print(\" 2. Sufficient disk space available\")\n print(\" 3. DICOM files are valid\")\n \n except KeyboardInterrupt:\n print(f\"\\n{'='*80}\")\n print(\"⚠️ INTERRUPTED BY USER\")\n print(f\"{'='*80}\")\n print(f\"\\nPartial progress saved. Re-run to resume from checkpoint.\")\n \n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ FATAL ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()\n \n print(f\"\\n💡 Troubleshooting:\")\n print(f\" 1. Check if dataset path is correct\")\n print(f\" 2. Ensure sufficient RAM (need ~8GB)\")\n print(f\" 3. Verify DICOM files are accessible\")\n print(f\" 4. Check disk space (~500MB needed)\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-01-22T08:25:43.775883Z","iopub.execute_input":"2026-01-22T08:25:43.776911Z","iopub.status.idle":"2026-01-22T09:46:01.783669Z","shell.execute_reply.started":"2026-01-22T08:25:43.776867Z","shell.execute_reply":"2026-01-22T09:46:01.782951Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nSTEP 2: COMPREHENSIVE DATA EXPLORATION & QUALITY ANALYSIS\n==========================================================\nVerify your enhanced 3-channel preprocessing and understand your data\n\nFeatures:\n✅ Multi-channel volume visualization\n✅ Statistical analysis & distribution plots\n✅ Quality verification (checks preprocessing integrity)\n✅ Class balance analysis\n✅ Channel comparison (bone/soft_tissue/wide)\n✅ Slice quality heatmaps\n✅ Interactive exploration\n✅ Simple baseline model (single-fold verification)\n\nRun this AFTER Step 1 to verify everything works!\n\"\"\"\n\nimport os\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom matplotlib.gridspec import GridSpec\nfrom tqdm import tqdm\n\nprint(\"=\"*80)\nprint(\"🔍 STEP 2: DATA EXPLORATION & QUALITY ANALYSIS\")\nprint(\"=\"*80)\n\n# ============================================================================\n# CONFIGURATION\n# ============================================================================\n\nCONFIG = {\n 'dataset_dir': '/kaggle/working/enhanced_dataset_v3',\n 'output_dir': '/kaggle/working/step2_exploration',\n \n # Visualization settings\n 'n_samples_to_visualize': 5, # Number of patients to visualize in detail\n 'slice_samples': 9, # Number of slices to show per patient\n \n # Quality checks\n 'run_integrity_checks': True,\n 'check_all_files': True, # Check all .npy files (slow for large datasets)\n \n # Simple baseline\n 'run_baseline': True, # Train simple model on 1 fold\n 'baseline_epochs': 5,\n 'baseline_batch_size': 4,\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\nprint(f\"\\n📋 Configuration:\")\nprint(f\" • Dataset: {CONFIG['dataset_dir']}\")\nprint(f\" • Output: {CONFIG['output_dir']}\")\nprint(f\" • Samples to visualize: {CONFIG['n_samples_to_visualize']}\")\nprint(f\" • Run baseline: {CONFIG['run_baseline']}\")\n\n# ============================================================================\n# 1. LOAD METADATA & BASIC STATISTICS\n# ============================================================================\n\ndef load_and_analyze_metadata():\n \"\"\"Load metadata and perform basic analysis\"\"\"\n \n print(f\"\\n{'='*80}\")\n print(\"📊 PART 1: METADATA ANALYSIS\")\n print(f\"{'='*80}\")\n \n metadata_path = os.path.join(CONFIG['dataset_dir'], 'metadata.csv')\n \n if not os.path.exists(metadata_path):\n print(f\"\\n❌ ERROR: {metadata_path} not found!\")\n print(f\" Please run Step 1 first!\")\n return None\n \n metadata_df = pd.read_csv(metadata_path)\n print(f\"\\n ✓ Loaded metadata: {len(metadata_df)} patients\")\n \n # Basic statistics\n print(f\"\\n 📈 Dataset Statistics:\")\n print(f\" • Total patients: {len(metadata_df)}\")\n print(f\" • Fracture cases: {metadata_df['has_fracture'].sum()} ({metadata_df['has_fracture'].mean()*100:.1f}%)\")\n print(f\" • Normal cases: {(1-metadata_df['has_fracture']).sum()} ({(1-metadata_df['has_fracture'].mean())*100:.1f}%)\")\n \n # Per-vertebra fracture statistics\n vertebrae = ['c1', 'c2', 'c3', 'c4', 'c5', 'c6', 'c7']\n print(f\"\\n 🦴 Per-Vertebra Fracture Distribution:\")\n for vert in vertebrae:\n if vert in metadata_df.columns:\n count = metadata_df[vert].sum()\n pct = count / len(metadata_df) * 100\n print(f\" • {vert.upper()}: {count} ({pct:.1f}%)\")\n \n # Quality metrics\n print(f\"\\n 📏 Volume Statistics:\")\n print(f\" • Avg slices: {metadata_df['num_slices'].mean():.1f} ± {metadata_df['num_slices'].std():.1f}\")\n print(f\" • Slice range: [{metadata_df['num_slices'].min()}, {metadata_df['num_slices'].max()}]\")\n print(f\" • Avg file size: {metadata_df['file_size_mb'].mean():.2f} MB\")\n print(f\" • Total size: {metadata_df['file_size_mb'].sum():.2f} MB\")\n \n # HU statistics\n print(f\"\\n 🔬 Hounsfield Unit (HU) Statistics:\")\n print(f\" • HU min range: [{metadata_df['hu_min'].min():.0f}, {metadata_df['hu_min'].max():.0f}]\")\n print(f\" • HU max range: [{metadata_df['hu_max'].min():.0f}, {metadata_df['hu_max'].max():.0f}]\")\n print(f\" • HU mean: {metadata_df['hu_mean'].mean():.1f} ± {metadata_df['hu_mean'].std():.1f}\")\n \n # Multi-channel statistics (if available)\n if 'bone_mean' in metadata_df.columns:\n print(f\"\\n 🎨 Multi-Channel Statistics:\")\n print(f\" • Bone channel mean: {metadata_df['bone_mean'].mean():.4f} ± {metadata_df['bone_mean'].std():.4f}\")\n print(f\" • Soft tissue mean: {metadata_df['soft_tissue_mean'].mean():.4f} ± {metadata_df['soft_tissue_mean'].std():.4f}\")\n print(f\" • Wide channel mean: {metadata_df['wide_mean'].mean():.4f} ± {metadata_df['wide_mean'].std():.4f}\")\n \n return metadata_df\n\n\n# ============================================================================\n# 2. VISUALIZATION FUNCTIONS\n# ============================================================================\n\ndef visualize_multi_channel_volume(volume, patient_id, label, save_path):\n \"\"\"\n Visualize all 3 channels of a single volume\n \n Args:\n volume: (D, H, W, 3) numpy array\n patient_id: patient identifier\n label: fracture label (0 or 1)\n save_path: where to save the figure\n \"\"\"\n depth = volume.shape[0]\n slice_indices = np.linspace(0, depth-1, CONFIG['slice_samples'], dtype=int)\n \n fig = plt.figure(figsize=(20, 10))\n gs = GridSpec(3, CONFIG['slice_samples'], figure=fig, hspace=0.3, wspace=0.1)\n \n channel_names = ['Bone Window', 'Soft Tissue Window', 'Wide Window']\n \n for ch_idx, ch_name in enumerate(channel_names):\n for col_idx, slice_idx in enumerate(slice_indices):\n ax = fig.add_subplot(gs[ch_idx, col_idx])\n ax.imshow(volume[slice_idx, :, :, ch_idx], cmap='gray', vmin=0, vmax=1)\n \n if col_idx == 0:\n ax.set_ylabel(ch_name, fontsize=11, fontweight='bold')\n \n if ch_idx == 0:\n ax.set_title(f'Slice {slice_idx}', fontsize=9)\n \n ax.axis('off')\n \n label_text = \"FRACTURE\" if label == 1 else \"NORMAL\"\n label_color = \"red\" if label == 1 else \"green\"\n \n fig.suptitle(f'Patient: {patient_id} | Label: {label_text}', \n fontsize=16, fontweight='bold', color=label_color, y=0.98)\n \n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n return save_path\n\n\ndef visualize_channel_comparison(volume, patient_id, save_path):\n \"\"\"\n Compare the same slice across all 3 channels\n \"\"\"\n mid_slice = volume.shape[0] // 2\n \n fig, axes = plt.subplots(1, 3, figsize=(15, 5))\n \n channel_names = ['Bone Window', 'Soft Tissue Window', 'Wide Window']\n \n for idx, (ax, ch_name) in enumerate(zip(axes, channel_names)):\n ax.imshow(volume[mid_slice, :, :, idx], cmap='gray', vmin=0, vmax=1)\n ax.set_title(f'{ch_name}\\nMean: {volume[:, :, :, idx].mean():.4f}', \n fontsize=12, fontweight='bold')\n ax.axis('off')\n \n plt.suptitle(f'Channel Comparison - Patient {patient_id} (Mid Slice)', \n fontsize=14, fontweight='bold')\n plt.tight_layout()\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n\n\ndef plot_statistical_distributions(metadata_df, save_path):\n \"\"\"\n Create comprehensive distribution plots\n \"\"\"\n fig = plt.figure(figsize=(20, 12))\n gs = GridSpec(3, 4, figure=fig, hspace=0.3, wspace=0.3)\n \n # 1. Class distribution\n ax1 = fig.add_subplot(gs[0, 0])\n class_counts = metadata_df['has_fracture'].value_counts()\n colors = ['green', 'red']\n bars = ax1.bar(['Normal', 'Fracture'], class_counts.values, color=colors, edgecolor='black', linewidth=2)\n ax1.set_ylabel('Count', fontsize=11, fontweight='bold')\n ax1.set_title('Class Distribution', fontsize=12, fontweight='bold')\n for bar in bars:\n height = bar.get_height()\n ax1.text(bar.get_x() + bar.get_width()/2., height,\n f'{int(height)}', ha='center', va='bottom', fontweight='bold')\n ax1.grid(True, alpha=0.3, axis='y')\n \n # 2. Per-vertebra distribution\n ax2 = fig.add_subplot(gs[0, 1])\n vertebrae = ['c1', 'c2', 'c3', 'c4', 'c5', 'c6', 'c7']\n vert_counts = [metadata_df[v].sum() if v in metadata_df.columns else 0 for v in vertebrae]\n ax2.bar([v.upper() for v in vertebrae], vert_counts, color='steelblue', edgecolor='black')\n ax2.set_ylabel('Fracture Count', fontsize=11, fontweight='bold')\n ax2.set_title('Per-Vertebra Fractures', fontsize=12, fontweight='bold')\n ax2.grid(True, alpha=0.3, axis='y')\n \n # 3. Number of slices distribution\n ax3 = fig.add_subplot(gs[0, 2])\n ax3.hist(metadata_df['num_slices'], bins=30, color='purple', alpha=0.7, edgecolor='black')\n ax3.axvline(metadata_df['num_slices'].mean(), color='red', linestyle='--', \n linewidth=2, label=f\"Mean: {metadata_df['num_slices'].mean():.1f}\")\n ax3.set_xlabel('Number of Slices', fontsize=11)\n ax3.set_ylabel('Frequency', fontsize=11)\n ax3.set_title('Slice Count Distribution', fontsize=12, fontweight='bold')\n ax3.legend()\n ax3.grid(True, alpha=0.3)\n \n # 4. File size distribution\n ax4 = fig.add_subplot(gs[0, 3])\n ax4.hist(metadata_df['file_size_mb'], bins=30, color='orange', alpha=0.7, edgecolor='black')\n ax4.axvline(metadata_df['file_size_mb'].mean(), color='red', linestyle='--', \n linewidth=2, label=f\"Mean: {metadata_df['file_size_mb'].mean():.2f} MB\")\n ax4.set_xlabel('File Size (MB)', fontsize=11)\n ax4.set_ylabel('Frequency', fontsize=11)\n ax4.set_title('File Size Distribution', fontsize=12, fontweight='bold')\n ax4.legend()\n ax4.grid(True, alpha=0.3)\n \n # 5. HU statistics\n ax5 = fig.add_subplot(gs[1, 0])\n ax5.hist(metadata_df['hu_mean'], bins=30, color='brown', alpha=0.7, edgecolor='black')\n ax5.set_xlabel('Mean HU', fontsize=11)\n ax5.set_ylabel('Frequency', fontsize=11)\n ax5.set_title('HU Mean Distribution', fontsize=12, fontweight='bold')\n ax5.grid(True, alpha=0.3)\n \n ax6 = fig.add_subplot(gs[1, 1])\n ax6.hist(metadata_df['hu_std'], bins=30, color='teal', alpha=0.7, edgecolor='black')\n ax6.set_xlabel('HU Std Dev', fontsize=11)\n ax6.set_ylabel('Frequency', fontsize=11)\n ax6.set_title('HU Std Distribution', fontsize=12, fontweight='bold')\n ax6.grid(True, alpha=0.3)\n \n # 7-9. Multi-channel statistics (if available)\n if 'bone_mean' in metadata_df.columns:\n ax7 = fig.add_subplot(gs[1, 2])\n metadata_df.boxplot(column=['bone_mean', 'soft_tissue_mean', 'wide_mean'], ax=ax7)\n ax7.set_ylabel('Normalized Intensity', fontsize=11)\n ax7.set_title('Channel Intensity Comparison', fontsize=12, fontweight='bold')\n ax7.grid(True, alpha=0.3, axis='y')\n \n # Channel correlation\n ax8 = fig.add_subplot(gs[1, 3])\n channel_cols = ['bone_mean', 'soft_tissue_mean', 'wide_mean']\n corr = metadata_df[channel_cols].corr()\n sns.heatmap(corr, annot=True, fmt='.3f', cmap='coolwarm', center=0,\n square=True, ax=ax8, cbar_kws={'label': 'Correlation'})\n ax8.set_title('Channel Correlation', fontsize=12, fontweight='bold')\n \n # 10. Slice quality metrics\n ax9 = fig.add_subplot(gs[2, 0])\n if 'empty_slices' in metadata_df.columns:\n ax9.hist(metadata_df['empty_slices'], bins=20, color='gray', alpha=0.7, edgecolor='black')\n ax9.set_xlabel('Empty Slices', fontsize=11)\n ax9.set_ylabel('Frequency', fontsize=11)\n ax9.set_title('Empty Slice Distribution', fontsize=12, fontweight='bold')\n ax9.grid(True, alpha=0.3)\n \n # 11. Fracture vs Normal comparison\n ax10 = fig.add_subplot(gs[2, 1:3])\n fracture_slices = metadata_df[metadata_df['has_fracture']==1]['num_slices']\n normal_slices = metadata_df[metadata_df['has_fracture']==0]['num_slices']\n \n ax10.hist([normal_slices, fracture_slices], bins=20, label=['Normal', 'Fracture'],\n color=['green', 'red'], alpha=0.6, edgecolor='black')\n ax10.set_xlabel('Number of Slices', fontsize=11)\n ax10.set_ylabel('Frequency', fontsize=11)\n ax10.set_title('Slice Count: Fracture vs Normal', fontsize=12, fontweight='bold')\n ax10.legend()\n ax10.grid(True, alpha=0.3)\n \n # 12. Summary text\n ax11 = fig.add_subplot(gs[2, 3])\n ax11.axis('off')\n \n summary_text = f\"\"\"\n DATASET SUMMARY\n {'='*30}\n \n Total Patients: {len(metadata_df)}\n Fracture: {metadata_df['has_fracture'].sum()}\n Normal: {(1-metadata_df['has_fracture']).sum()}\n \n Avg Slices: {metadata_df['num_slices'].mean():.1f}\n Avg Size: {metadata_df['file_size_mb'].mean():.2f} MB\n Total Size: {metadata_df['file_size_mb'].sum():.1f} MB\n \n HU Range: [{metadata_df['hu_min'].min():.0f}, \n {metadata_df['hu_max'].max():.0f}]\n \"\"\"\n \n ax11.text(0.1, 0.5, summary_text, fontsize=10, family='monospace',\n verticalalignment='center',\n bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))\n \n plt.suptitle('Comprehensive Statistical Analysis', fontsize=16, fontweight='bold', y=0.98)\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved: {save_path}\")\n\n\n# ============================================================================\n# 3. QUALITY VERIFICATION\n# ============================================================================\n\ndef verify_data_integrity(metadata_df):\n \"\"\"\n Verify that all volumes are correctly saved and loadable\n \"\"\"\n print(f\"\\n{'='*80}\")\n print(\"🔍 PART 2: DATA INTEGRITY VERIFICATION\")\n print(f\"{'='*80}\")\n \n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n issues = []\n verified = 0\n \n n_to_check = len(metadata_df) if CONFIG['check_all_files'] else min(20, len(metadata_df))\n \n print(f\"\\n Checking {n_to_check} volumes...\")\n \n for idx in tqdm(range(n_to_check), desc=\"Verifying\"):\n row = metadata_df.iloc[idx]\n patient_id = row['patient_id']\n volume_path = os.path.join(volumes_dir, f\"{patient_id}.npy\")\n \n try:\n # Check file exists\n if not os.path.exists(volume_path):\n issues.append(f\"Missing: {patient_id}\")\n continue\n \n # Load volume\n volume = np.load(volume_path)\n \n # Check shape\n expected_shape = (64, 224, 224, 3) # (D, H, W, C)\n if volume.shape != expected_shape:\n issues.append(f\"{patient_id}: Wrong shape {volume.shape}, expected {expected_shape}\")\n continue\n \n # Check data type\n if volume.dtype != np.float32:\n issues.append(f\"{patient_id}: Wrong dtype {volume.dtype}, expected float32\")\n continue\n \n # Check value range\n if volume.min() < 0 or volume.max() > 1:\n issues.append(f\"{patient_id}: Values out of range [{volume.min():.3f}, {volume.max():.3f}]\")\n continue\n \n # Check for NaN/Inf\n if np.isnan(volume).any() or np.isinf(volume).any():\n issues.append(f\"{patient_id}: Contains NaN or Inf\")\n continue\n \n verified += 1\n \n except Exception as e:\n issues.append(f\"{patient_id}: Load error - {str(e)[:50]}\")\n \n # Report\n print(f\"\\n 📊 Verification Results:\")\n print(f\" • Verified: {verified}/{n_to_check}\")\n print(f\" • Issues: {len(issues)}\")\n \n if len(issues) > 0:\n print(f\"\\n ⚠️ Issues found:\")\n for issue in issues[:10]: # Show first 10\n print(f\" - {issue}\")\n \n if len(issues) > 10:\n print(f\" ... and {len(issues)-10} more\")\n \n # Save issues to file\n issues_path = os.path.join(CONFIG['output_dir'], 'integrity_issues.txt')\n with open(issues_path, 'w') as f:\n f.write('\\n'.join(issues))\n print(f\"\\n ✓ Full issue list saved: {issues_path}\")\n else:\n print(f\"\\n ✅ All volumes passed integrity checks!\")\n \n return verified, issues\n\n\n# ============================================================================\n# 4. INTERACTIVE EXPLORATION\n# ============================================================================\n\ndef explore_sample_volumes(metadata_df):\n \"\"\"\n Visualize sample volumes in detail\n \"\"\"\n print(f\"\\n{'='*80}\")\n print(\"🎨 PART 3: SAMPLE VOLUME VISUALIZATION\")\n print(f\"{'='*80}\")\n \n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n # Select diverse samples\n fracture_samples = metadata_df[metadata_df['has_fracture']==1].sample(\n n=min(CONFIG['n_samples_to_visualize']//2, metadata_df['has_fracture'].sum()),\n random_state=42\n )\n normal_samples = metadata_df[metadata_df['has_fracture']==0].sample(\n n=min(CONFIG['n_samples_to_visualize']//2, (1-metadata_df['has_fracture']).sum()),\n random_state=42\n )\n \n samples = pd.concat([fracture_samples, normal_samples])\n \n print(f\"\\n Visualizing {len(samples)} sample volumes...\")\n \n for idx, row in tqdm(samples.iterrows(), total=len(samples), desc=\"Creating visualizations\"):\n patient_id = row['patient_id']\n label = row['has_fracture']\n \n volume_path = os.path.join(volumes_dir, f\"{patient_id}.npy\")\n \n if not os.path.exists(volume_path):\n continue\n \n volume = np.load(volume_path)\n \n # Multi-channel visualization\n save_path = os.path.join(CONFIG['output_dir'], f'volume_{patient_id}_all_channels.png')\n visualize_multi_channel_volume(volume, patient_id, label, save_path)\n \n # Channel comparison\n save_path = os.path.join(CONFIG['output_dir'], f'volume_{patient_id}_channel_comparison.png')\n visualize_channel_comparison(volume, patient_id, save_path)\n \n print(f\"\\n ✓ Visualizations saved to: {CONFIG['output_dir']}/\")\n\n\n# ============================================================================\n# 5. SIMPLE BASELINE MODEL\n# ============================================================================\n\ndef run_simple_baseline(metadata_df):\n \"\"\"\n Train a simple 3D CNN on first fold to verify pipeline works\n \"\"\"\n print(f\"\\n{'='*80}\")\n print(\"🧪 PART 4: SIMPLE BASELINE MODEL (1-FOLD VERIFICATION)\")\n print(f\"{'='*80}\")\n \n try:\n import torch\n import torch.nn as nn\n import torch.optim as optim\n from torch.utils.data import Dataset, DataLoader\n from sklearn.metrics import roc_auc_score, accuracy_score\n except ImportError:\n print(\"\\n ⚠️ PyTorch not available, skipping baseline\")\n return None\n \n # Load CV splits\n cv_path = os.path.join(CONFIG['dataset_dir'], 'cv_splits.json')\n if not os.path.exists(cv_path):\n print(f\"\\n ⚠️ CV splits not found, skipping baseline\")\n return None\n \n with open(cv_path, 'r') as f:\n cv_splits = json.load(f)\n \n # Use first fold\n fold_0 = cv_splits[0]\n train_patients = fold_0['train']\n val_patients = fold_0['val']\n \n print(f\"\\n Using Fold 0:\")\n print(f\" Train: {len(train_patients)} patients\")\n print(f\" Val: {len(val_patients)} patients\")\n \n # Simple dataset\n class SimpleDataset(Dataset):\n def __init__(self, patient_ids, metadata_df, volumes_dir):\n self.patient_ids = patient_ids\n self.metadata_df = metadata_df.set_index('patient_id')\n self.volumes_dir = volumes_dir\n \n def __len__(self):\n return len(self.patient_ids)\n \n def __getitem__(self, idx):\n patient_id = self.patient_ids[idx]\n \n # Load volume\n volume = np.load(os.path.join(self.volumes_dir, f\"{patient_id}.npy\"))\n \n # Transpose to (C, D, H, W) for PyTorch\n volume = np.transpose(volume, (3, 0, 1, 2)) # (D,H,W,3) -> (3,D,H,W)\n \n label = self.metadata_df.loc[patient_id, 'has_fracture']\n \n return torch.from_numpy(volume).float(), torch.tensor(label).float()\n \n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n train_dataset = SimpleDataset(train_patients, metadata_df, volumes_dir)\n val_dataset = SimpleDataset(val_patients, metadata_df, volumes_dir)\n \n train_loader = DataLoader(train_dataset, batch_size=CONFIG['baseline_batch_size'], \n shuffle=True, num_workers=2)\n val_loader = DataLoader(val_dataset, batch_size=CONFIG['baseline_batch_size'], \n shuffle=False, num_workers=2)\n \n print(f\"\\n Batches: Train={len(train_loader)}, Val={len(val_loader)}\")\n \n # Simple 3D CNN\n class SimpleCNN3D(nn.Module):\n def __init__(self):\n super().__init__()\n \n self.conv1 = nn.Conv3d(3, 16, kernel_size=3, stride=2, padding=1) # 3 channels!\n self.bn1 = nn.BatchNorm3d(16)\n self.conv2 = nn.Conv3d(16, 32, kernel_size=3, stride=2, padding=1)\n self.bn2 = nn.BatchNorm3d(32)\n self.conv3 = nn.Conv3d(32, 64, kernel_size=3, stride=2, padding=1)\n self.bn3 = nn.BatchNorm3d(64)\n \n self.pool = nn.AdaptiveAvgPool3d(1)\n self.fc = nn.Linear(64, 1)\n \n def forward(self, x):\n x = torch.relu(self.bn1(self.conv1(x)))\n x = torch.relu(self.bn2(self.conv2(x)))\n x = torch.relu(self.bn3(self.conv3(x)))\n x = self.pool(x)\n x = x.view(x.size(0), -1)\n x = self.fc(x)\n return x.squeeze(-1)\n \n device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n print(f\" Device: {device}\")\n \n model = SimpleCNN3D().to(device)\n criterion = nn.BCEWithLogitsLoss()\n optimizer = optim.Adam(model.parameters(), lr=1e-3)\n \n print(f\"\\n Training {CONFIG['baseline_epochs']} epochs...\")\n \n history = {'train_loss': [], 'val_loss': [], 'val_auc': [], 'val_acc': []}\n \n for epoch in range(CONFIG['baseline_epochs']):\n # Train\n model.train()\n train_loss = 0\n for volumes, labels in train_loader:\n volumes, labels = volumes.to(device), labels.to(device)\n \n optimizer.zero_grad()\n outputs = model(volumes)\n loss = criterion(outputs, labels)\n loss.backward()\n optimizer.step()\n \n train_loss += loss.item()\n \n train_loss /= len(train_loader)\n \n # Validate\n model.eval()\n val_loss = 0\n all_preds = []\n all_labels = []\n \n with torch.no_grad():\n for volumes, labels in val_loader:\n volumes, labels = volumes.to(device), labels.to(device)\n \n outputs = model(volumes)\n loss = criterion(outputs, labels)\n val_loss += loss.item()\n \n preds = torch.sigmoid(outputs).cpu().numpy()\n all_preds.extend(preds)\n all_labels.extend(labels.cpu().numpy())\n \n val_loss /= len(val_loader)\n \n all_preds = np.array(all_preds)\n all_labels = np.array(all_labels)\n \n val_auc = roc_auc_score(all_labels, all_preds) if len(np.unique(all_labels)) > 1 else 0.5\n val_acc = accuracy_score(all_labels, (all_preds > 0.5).astype(int))\n \n history['train_loss'].append(train_loss)\n history['val_loss'].append(val_loss)\n history['val_auc'].append(val_auc)\n history['val_acc'].append(val_acc)\n \n print(f\" Epoch {epoch+1}/{CONFIG['baseline_epochs']} - \"\n f\"Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}, \"\n f\"Val AUC: {val_auc:.4f}, Val Acc: {val_acc:.4f}\")\n \n # Plot baseline results\n fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n \n epochs = range(1, CONFIG['baseline_epochs']+1)\n \n axes[0].plot(epochs, history['train_loss'], 'b-', label='Train', linewidth=2)\n axes[0].plot(epochs, history['val_loss'], 'r-', label='Val', linewidth=2)\n axes[0].set_xlabel('Epoch', fontsize=11)\n axes[0].set_ylabel('Loss', fontsize=11)\n axes[0].set_title('Baseline Loss', fontsize=12, fontweight='bold')\n axes[0].legend()\n axes[0].grid(True, alpha=0.3)\n \n axes[1].plot(epochs, history['val_auc'], 'g-', label='AUC', linewidth=2, marker='o')\n axes[1].plot(epochs, history['val_acc'], 'orange', label='Accuracy', linewidth=2, marker='s')\n axes[1].axhline(y=0.5, color='gray', linestyle='--', alpha=0.5)\n axes[1].set_xlabel('Epoch', fontsize=11)\n axes[1].set_ylabel('Score', fontsize=11)\n axes[1].set_title('Baseline Metrics', fontsize=12, fontweight='bold')\n axes[1].legend()\n axes[1].grid(True, alpha=0.3)\n \n plt.tight_layout()\n baseline_path = os.path.join(CONFIG['output_dir'], 'baseline_results.png')\n plt.savefig(baseline_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\"\\n ✓ Baseline training complete!\")\n print(f\" 📊 Final metrics:\")\n print(f\" Val AUC: {history['val_auc'][-1]:.4f}\")\n print(f\" Val Accuracy: {history['val_acc'][-1]:.4f}\")\n \n if history['val_auc'][-1] > 0.55:\n print(f\"\\n ✅ Pipeline verified! AUC > 0.55 shows model is learning\")\n else:\n print(f\"\\n ⚠️ Low AUC - may need more epochs or data\")\n \n return history\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n \"\"\"Run complete exploration pipeline\"\"\"\n \n print(f\"\\n🚀 Starting comprehensive data exploration...\")\n \n # 1. Load and analyze metadata\n metadata_df = load_and_analyze_metadata()\n \n if metadata_df is None:\n return\n \n # 2. Create statistical visualizations\n print(f\"\\n{'='*80}\")\n print(\"📊 Creating statistical visualizations...\")\n print(f\"{'='*80}\")\n \n stats_path = os.path.join(CONFIG['output_dir'], 'statistical_analysis.png')\n plot_statistical_distributions(metadata_df, stats_path)\n \n # 3. Verify data integrity\n if CONFIG['run_integrity_checks']:\n verified, issues = verify_data_integrity(metadata_df)\n \n # 4. Explore sample volumes\n explore_sample_volumes(metadata_df)\n \n # 5. Run baseline model\n if CONFIG['run_baseline']:\n baseline_history = run_simple_baseline(metadata_df)\n \n # Final summary\n print(f\"\\n{'='*80}\")\n print(\"✅ EXPLORATION COMPLETE!\")\n print(f\"{'='*80}\")\n \n print(f\"\\n📁 Generated files in {CONFIG['output_dir']}:\")\n print(f\" • statistical_analysis.png - Comprehensive statistics\")\n print(f\" • volume_*_all_channels.png - Multi-channel visualizations\")\n print(f\" • volume_*_channel_comparison.png - Channel comparisons\")\n \n if CONFIG['run_integrity_checks']:\n print(f\" • integrity_issues.txt - Data quality issues (if any)\")\n \n if CONFIG['run_baseline']:\n print(f\" • baseline_results.png - Simple model verification\")\n \n print(f\"\\n🎯 Key Findings:\")\n print(f\" ✓ Dataset loaded: {len(metadata_df)} patients\")\n print(f\" ✓ Multi-channel preprocessing verified\")\n print(f\" ✓ Class balance: {metadata_df['has_fracture'].mean()*100:.1f}% fracture\")\n \n if CONFIG['run_integrity_checks']:\n if len(issues) == 0:\n print(f\" ✓ All volumes passed integrity checks\")\n else:\n print(f\" ⚠️ {len(issues)} volumes have issues - check integrity_issues.txt\")\n \n if CONFIG['run_baseline']:\n if baseline_history and baseline_history['val_auc'][-1] > 0.55:\n print(f\" ✓ Baseline model works (AUC={baseline_history['val_auc'][-1]:.4f})\")\n else:\n print(f\" ⚠️ Baseline AUC is low - may need tuning\")\n \n print(f\"\\n💡 Next Steps:\")\n print(f\" 1. Review visualizations to understand your data\")\n print(f\" 2. Check quality report and fix any issues\")\n print(f\" 3. If baseline works, proceed to full training (Step 3)\")\n print(f\" 4. Use CV splits for robust evaluation\")\n \n print(f\"\\n{'='*80}\")\n print(\"🎉 Ready for full-scale training!\")\n print(f\"{'='*80}\")\n\n\nif __name__ == \"__main__\":\n try:\n main()\n except KeyboardInterrupt:\n print(f\"\\n\\n⚠️ Interrupted by user\")\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T09:46:01.785173Z","iopub.execute_input":"2026-01-22T09:46:01.785539Z","iopub.status.idle":"2026-01-22T09:47:13.054147Z","shell.execute_reply.started":"2026-01-22T09:46:01.785517Z","shell.execute_reply":"2026-01-22T09:47:13.053264Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nENHANCED STEP 2: ADVANCED MULTI-CHANNEL AUGMENTATION v3.0\n==========================================================\nProduction-grade augmentation matching Step 1's preprocessing\n\nKEY IMPROVEMENTS:\n✅ Multi-channel support (3-channel volumes)\n✅ Elastic deformations (medical imaging essential!)\n✅ MixUp augmentation\n✅ Channel-specific strategies\n✅ Test-Time Augmentation (TTA)\n✅ Advanced intensity transforms\n\nExpected: +10-15% accuracy boost\n\"\"\"\n\nimport os\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom matplotlib.gridspec import GridSpec\nfrom scipy.ndimage import rotate, zoom, gaussian_filter, shift, map_coordinates\nfrom tqdm import tqdm\n\nprint(\"=\"*80)\nprint(\"🚀 ENHANCED MULTI-CHANNEL AUGMENTATION v3.0\")\nprint(\"=\"*80)\n\nCONFIG = {\n 'dataset_dir': '/kaggle/working/enhanced_dataset_v3',\n 'output_dir': '/kaggle/working/augmentation_analysis',\n 'use_elastic': True,\n 'use_mixup': True,\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\n# ============================================================================\n# ENHANCED AUGMENTATION CLASS\n# ============================================================================\n\nclass EnhancedAugmentation3D:\n \"\"\"Medical imaging optimized augmentation for multi-channel volumes\"\"\"\n \n def __init__(self):\n self.train_config = {\n # Spatial\n 'flip_prob': 0.5,\n 'rotation': (-15, 15),\n 'rotation_prob': 0.4,\n 'zoom': (0.85, 1.15),\n 'zoom_prob': 0.3,\n 'shift': (-15, 15),\n 'shift_prob': 0.3,\n 'elastic_alpha': (0, 30),\n 'elastic_prob': 0.3,\n \n # Intensity\n 'brightness': (0.85, 1.15),\n 'brightness_prob': 0.4,\n 'contrast': (0.85, 1.15),\n 'contrast_prob': 0.4,\n 'gamma': (0.8, 1.2),\n 'gamma_prob': 0.3,\n 'noise_std': 0.02,\n 'noise_prob': 0.25,\n 'blur': (0.5, 2.0),\n 'blur_prob': 0.15,\n }\n \n self.tta_config = { # Lighter for test-time\n 'flip_prob': 0.5,\n 'rotation': (-5, 5),\n 'rotation_prob': 0.3,\n 'zoom': (0.95, 1.05),\n 'zoom_prob': 0.2,\n }\n \n def horizontal_flip(self, vol):\n \"\"\"Flip along width\"\"\"\n return np.flip(vol, axis=2).copy()\n \n def random_rotation(self, vol, angle=None):\n \"\"\"Rotate in axial plane\"\"\"\n if angle is None:\n angle = np.random.uniform(*self.train_config['rotation'])\n \n # Handle multi-channel (D,H,W,3)\n if vol.ndim == 4:\n rotated = np.zeros_like(vol)\n for ch in range(vol.shape[3]):\n rotated[..., ch] = rotate(vol[..., ch], angle, \n axes=(1,2), reshape=False, order=1)\n return rotated.astype(np.float32)\n else:\n return rotate(vol, angle, axes=(1,2), reshape=False, order=1).astype(np.float32)\n \n def random_zoom(self, vol, scale=None):\n \"\"\"Zoom spatial dimensions\"\"\"\n if scale is None:\n scale = np.random.uniform(*self.train_config['zoom'])\n \n orig_shape = vol.shape\n \n if vol.ndim == 4: # Multi-channel\n factors = [scale, scale, scale, 1.0] # Don't zoom channels!\n zoomed = zoom(vol, factors, order=1)\n \n if scale > 1.0: # Crop\n start = [(zoomed.shape[i] - orig_shape[i])//2 for i in range(3)]\n result = zoomed[start[0]:start[0]+orig_shape[0],\n start[1]:start[1]+orig_shape[1],\n start[2]:start[2]+orig_shape[2], :]\n else: # Pad\n result = np.zeros(orig_shape, dtype=vol.dtype)\n start = [(orig_shape[i] - zoomed.shape[i])//2 for i in range(3)]\n result[start[0]:start[0]+zoomed.shape[0],\n start[1]:start[1]+zoomed.shape[1],\n start[2]:start[2]+zoomed.shape[2], :] = zoomed\n else:\n zoomed = zoom(vol, scale, order=1)\n if scale > 1.0:\n start = [(zoomed.shape[i] - orig_shape[i])//2 for i in range(3)]\n result = zoomed[start[0]:start[0]+orig_shape[0],\n start[1]:start[1]+orig_shape[1],\n start[2]:start[2]+orig_shape[2]]\n else:\n result = np.zeros(orig_shape, dtype=vol.dtype)\n start = [(orig_shape[i] - zoomed.shape[i])//2 for i in range(3)]\n result[start[0]:start[0]+zoomed.shape[0],\n start[1]:start[1]+zoomed.shape[1],\n start[2]:start[2]+zoomed.shape[2]] = zoomed\n \n return result.astype(np.float32)\n \n def elastic_deformation(self, vol, alpha=None):\n \"\"\"\n CRITICAL AUGMENTATION for medical imaging!\n Simulates realistic anatomical variations\n \"\"\"\n if alpha is None:\n alpha = np.random.uniform(*self.train_config['elastic_alpha'])\n if alpha == 0:\n return vol\n \n sigma = 3\n shape = vol.shape[:3]\n \n # Random displacement fields\n dx = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n dy = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n dz = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n \n # Coordinate grid\n d, h, w = shape\n d_coords, h_coords, w_coords = np.meshgrid(\n np.arange(d), np.arange(h), np.arange(w), indexing='ij'\n )\n \n indices = (\n np.reshape(d_coords + dx, (-1, 1)),\n np.reshape(h_coords + dy, (-1, 1)),\n np.reshape(w_coords + dz, (-1, 1))\n )\n \n if vol.ndim == 4: # Multi-channel\n deformed = np.zeros_like(vol)\n for ch in range(vol.shape[3]):\n deformed_ch = map_coordinates(vol[..., ch], indices, order=1, mode='nearest')\n deformed[..., ch] = deformed_ch.reshape(shape)\n else:\n deformed = map_coordinates(vol, indices, order=1, mode='nearest')\n deformed = deformed.reshape(shape)\n \n return deformed.astype(np.float32)\n \n def random_brightness(self, vol, factor=None):\n \"\"\"Adjust brightness\"\"\"\n if factor is None:\n factor = np.random.uniform(*self.train_config['brightness'])\n return np.clip(vol * factor, 0, 1).astype(np.float32)\n \n def random_contrast(self, vol, factor=None):\n \"\"\"Adjust contrast\"\"\"\n if factor is None:\n factor = np.random.uniform(*self.train_config['contrast'])\n mean = vol.mean()\n adjusted = (vol - mean) * factor + mean\n return np.clip(adjusted, 0, 1).astype(np.float32)\n \n def random_gamma(self, vol, gamma=None):\n \"\"\"Gamma correction (non-linear intensity)\"\"\"\n if gamma is None:\n gamma = np.random.uniform(*self.train_config['gamma'])\n return np.clip(np.power(vol, gamma), 0, 1).astype(np.float32)\n \n def gaussian_noise(self, vol):\n \"\"\"Add noise\"\"\"\n noise = np.random.normal(0, self.train_config['noise_std'], vol.shape)\n return np.clip(vol + noise, 0, 1).astype(np.float32)\n \n def gaussian_blur(self, vol, sigma=None):\n \"\"\"Blur\"\"\"\n if sigma is None:\n sigma = np.random.uniform(*self.train_config['blur'])\n \n if vol.ndim == 4:\n blurred = np.zeros_like(vol)\n for ch in range(vol.shape[3]):\n blurred[..., ch] = gaussian_filter(vol[..., ch], sigma=sigma)\n return blurred\n else:\n return gaussian_filter(vol, sigma=sigma).astype(np.float32)\n \n def mixup(self, vol1, vol2, label1, label2, alpha=0.2):\n \"\"\"\n MixUp: Mix two samples for better generalization\n \"\"\"\n lam = np.random.beta(alpha, alpha) if alpha > 0 else 1.0\n mixed_vol = lam * vol1 + (1 - lam) * vol2\n mixed_label = lam * label1 + (1 - lam) * label2\n return mixed_vol.astype(np.float32), float(mixed_label)\n \n def apply_train_augmentations(self, vol, seed=None):\n \"\"\"Full training augmentation pipeline\"\"\"\n if seed is not None:\n np.random.seed(seed)\n \n aug = vol.copy()\n cfg = self.train_config\n \n # Spatial\n if np.random.random() < cfg['flip_prob']:\n aug = self.horizontal_flip(aug)\n if np.random.random() < cfg['rotation_prob']:\n aug = self.random_rotation(aug)\n if np.random.random() < cfg['zoom_prob']:\n aug = self.random_zoom(aug)\n if np.random.random() < cfg['shift_prob']:\n shifts = [np.random.uniform(*cfg['shift']) for _ in range(3)]\n if aug.ndim == 4:\n shifts.append(0) # Don't shift channels\n aug = shift(aug, shifts, order=1, mode='nearest').astype(np.float32)\n if CONFIG['use_elastic'] and np.random.random() < cfg['elastic_prob']:\n aug = self.elastic_deformation(aug)\n \n # Intensity\n if np.random.random() < cfg['brightness_prob']:\n aug = self.random_brightness(aug)\n if np.random.random() < cfg['contrast_prob']:\n aug = self.random_contrast(aug)\n if np.random.random() < cfg['gamma_prob']:\n aug = self.random_gamma(aug)\n if np.random.random() < cfg['noise_prob']:\n aug = self.gaussian_noise(aug)\n if np.random.random() < cfg['blur_prob']:\n aug = self.gaussian_blur(aug)\n \n return aug\n\n\n# ============================================================================\n# VISUALIZATION\n# ============================================================================\n\ndef visualize_comprehensive_demo(volume, label, augmenter, save_dir):\n \"\"\"Create comprehensive augmentation demonstrations\"\"\"\n \n slice_idx = volume.shape[0] // 2\n \n # ========================================================================\n # DEMO 1: Multi-channel augmentation\n # ========================================================================\n \n fig = plt.figure(figsize=(20, 12))\n gs = GridSpec(4, 4, figure=fig, hspace=0.25, wspace=0.15)\n \n aug_types = [\n ('Original', None),\n ('Spatial Only', 100),\n ('Intensity Only', 200),\n ('Full Pipeline', 300)\n ]\n \n for row_idx, (title, seed) in enumerate(aug_types):\n for ch_idx, ch_name in enumerate(['Bone', 'Soft Tissue', 'Wide']):\n ax = fig.add_subplot(gs[row_idx, ch_idx])\n \n if row_idx == 0:\n img = volume[slice_idx, :, :, ch_idx]\n else:\n # Modify config temporarily for demo\n if seed == 100: # Spatial only\n backup = augmenter.train_config.copy()\n augmenter.train_config['brightness_prob'] = 0\n augmenter.train_config['contrast_prob'] = 0\n augmenter.train_config['gamma_prob'] = 0\n augmenter.train_config['noise_prob'] = 0\n augmenter.train_config['blur_prob'] = 0\n aug = augmenter.apply_train_augmentations(volume, seed=seed)\n augmenter.train_config = backup\n elif seed == 200: # Intensity only\n backup = augmenter.train_config.copy()\n augmenter.train_config['flip_prob'] = 0\n augmenter.train_config['rotation_prob'] = 0\n augmenter.train_config['zoom_prob'] = 0\n augmenter.train_config['shift_prob'] = 0\n augmenter.train_config['elastic_prob'] = 0\n aug = augmenter.apply_train_augmentations(volume, seed=seed)\n augmenter.train_config = backup\n else:\n aug = augmenter.apply_train_augmentations(volume, seed=seed)\n \n img = aug[slice_idx, :, :, ch_idx]\n \n ax.imshow(img, cmap='gray', vmin=0, vmax=1)\n if ch_idx == 0:\n ax.set_ylabel(title, fontsize=12, fontweight='bold')\n if row_idx == 0:\n ax.set_title(ch_name, fontsize=13, fontweight='bold')\n ax.axis('off')\n \n # RGB Composite column\n for row_idx, (_, seed) in enumerate(aug_types):\n ax = fig.add_subplot(gs[row_idx, 3])\n \n if row_idx == 0:\n vol = volume\n else:\n vol = augmenter.apply_train_augmentations(volume, seed=seed)\n \n rgb = np.stack([\n vol[slice_idx, :, :, 0],\n vol[slice_idx, :, :, 1],\n vol[slice_idx, :, :, 2]\n ], axis=-1)\n \n ax.imshow(rgb)\n if row_idx == 0:\n ax.set_title('RGB Composite', fontsize=13, fontweight='bold')\n ax.axis('off')\n \n label_text = 'FRACTURE' if label == 1 else 'NORMAL'\n plt.suptitle(f'Multi-Channel Augmentation Analysis | Label: {label_text}',\n fontsize=16, fontweight='bold')\n \n save_path = os.path.join(save_dir, 'demo1_multichannel.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n print(f\" ✓ Saved: demo1_multichannel.png\")\n \n # ========================================================================\n # DEMO 2: Elastic deformation showcase\n # ========================================================================\n \n fig, axes = plt.subplots(2, 4, figsize=(16, 8))\n \n alphas = [0, 10, 20, 30]\n for idx, alpha in enumerate(alphas):\n # Top row: deformed\n if alpha == 0:\n deformed = volume\n else:\n deformed = augmenter.elastic_deformation(volume, alpha=alpha)\n \n axes[0, idx].imshow(deformed[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n axes[0, idx].set_title(f'α={alpha}', fontsize=12, fontweight='bold')\n axes[0, idx].axis('off')\n \n # Bottom row: overlay with original (red=original, green=deformed)\n overlay = np.stack([\n volume[slice_idx, :, :, 0],\n deformed[slice_idx, :, :, 0],\n np.zeros_like(volume[slice_idx, :, :, 0])\n ], axis=-1)\n axes[1, idx].imshow(overlay)\n axes[1, idx].set_title('Overlay', fontsize=10)\n axes[1, idx].axis('off')\n \n plt.suptitle('Elastic Deformation - Critical for Medical Imaging',\n fontsize=14, fontweight='bold')\n plt.tight_layout()\n save_path = os.path.join(save_dir, 'demo2_elastic.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n print(f\" ✓ Saved: demo2_elastic.png\")\n \n # ========================================================================\n # DEMO 3: Multiple random samples\n # ========================================================================\n \n fig, axes = plt.subplots(3, 6, figsize=(18, 9))\n \n for row in range(3):\n # Original\n axes[row, 0].imshow(volume[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n if row == 0:\n axes[row, 0].set_title('Original', fontweight='bold')\n axes[row, 0].axis('off')\n \n # 5 augmented versions\n for col in range(1, 6):\n aug = augmenter.apply_train_augmentations(volume, seed=row*10+col)\n axes[row, col].imshow(aug[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n if row == 0:\n axes[row, col].set_title(f'Aug {col}', fontweight='bold')\n axes[row, col].axis('off')\n \n plt.suptitle('Augmentation Diversity - 15 Random Samples',\n fontsize=14, fontweight='bold')\n plt.tight_layout()\n save_path = os.path.join(save_dir, 'demo3_diversity.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n print(f\" ✓ Saved: demo3_diversity.png\")\n \n # ========================================================================\n # DEMO 4: Statistical analysis\n # ========================================================================\n \n print(f\" 📊 Generating statistical analysis (100 samples)...\")\n \n orig_stats = {\n 'mean': volume.mean(),\n 'std': volume.std(),\n 'ch0_mean': volume[..., 0].mean(),\n 'ch1_mean': volume[..., 1].mean(),\n 'ch2_mean': volume[..., 2].mean(),\n }\n \n aug_means = []\n aug_stds = []\n aug_ch0_means = []\n aug_ch1_means = []\n aug_ch2_means = []\n \n for i in range(100):\n aug = augmenter.apply_train_augmentations(volume, seed=1000+i)\n aug_means.append(aug.mean())\n aug_stds.append(aug.std())\n aug_ch0_means.append(aug[..., 0].mean())\n aug_ch1_means.append(aug[..., 1].mean())\n aug_ch2_means.append(aug[..., 2].mean())\n \n fig, axes = plt.subplots(2, 3, figsize=(15, 8))\n \n # Overall stats\n axes[0, 0].hist(aug_means, bins=30, alpha=0.7, color='blue', edgecolor='black')\n axes[0, 0].axvline(orig_stats['mean'], color='red', linestyle='--', linewidth=2)\n axes[0, 0].set_title('Overall Mean Distribution', fontweight='bold')\n axes[0, 0].set_xlabel('Mean Intensity')\n axes[0, 0].grid(True, alpha=0.3)\n \n axes[0, 1].hist(aug_stds, bins=30, alpha=0.7, color='green', edgecolor='black')\n axes[0, 1].axvline(orig_stats['std'], color='red', linestyle='--', linewidth=2)\n axes[0, 1].set_title('Overall Std Distribution', fontweight='bold')\n axes[0, 1].set_xlabel('Std Deviation')\n axes[0, 1].grid(True, alpha=0.3)\n \n # Channel-wise\n axes[0, 2].hist(aug_ch0_means, bins=30, alpha=0.7, color='brown', edgecolor='black')\n axes[0, 2].axvline(orig_stats['ch0_mean'], color='red', linestyle='--', linewidth=2)\n axes[0, 2].set_title('Bone Channel Mean', fontweight='bold')\n axes[0, 2].grid(True, alpha=0.3)\n \n axes[1, 0].hist(aug_ch1_means, bins=30, alpha=0.7, color='purple', edgecolor='black')\n axes[1, 0].axvline(orig_stats['ch1_mean'], color='red', linestyle='--', linewidth=2)\n axes[1, 0].set_title('Soft Tissue Channel Mean', fontweight='bold')\n axes[1, 0].grid(True, alpha=0.3)\n \n axes[1, 1].hist(aug_ch2_means, bins=30, alpha=0.7, color='orange', edgecolor='black')\n axes[1, 1].axvline(orig_stats['ch2_mean'], color='red', linestyle='--', linewidth=2)\n axes[1, 1].set_title('Wide Channel Mean', fontweight='bold')\n axes[1, 1].grid(True, alpha=0.3)\n \n # Summary text\n axes[1, 2].axis('off')\n summary = f\"\"\"\n AUGMENTATION STATISTICS\n {'='*30}\n \n Original Volume:\n • Overall: {orig_stats['mean']:.4f}±{orig_stats['std']:.4f}\n • Ch0 (Bone): {orig_stats['ch0_mean']:.4f}\n • Ch1 (Soft): {orig_stats['ch1_mean']:.4f}\n • Ch2 (Wide): {orig_stats['ch2_mean']:.4f}\n \n Augmented (n=100):\n • Mean: {np.mean(aug_means):.4f}±{np.std(aug_means):.4f}\n • Std: {np.mean(aug_stds):.4f}±{np.std(aug_stds):.4f}\n \n ✅ Intensity preserved\n ✅ Diversity achieved\n \"\"\"\n \n axes[1, 2].text(0.1, 0.5, summary, fontsize=10, family='monospace',\n verticalalignment='center',\n bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.3))\n \n plt.suptitle('Statistical Stability Analysis', fontsize=14, fontweight='bold')\n plt.tight_layout()\n save_path = os.path.join(save_dir, 'demo4_statistics.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n print(f\" ✓ Saved: demo4_statistics.png\")\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n print(f\"\\n{'='*80}\")\n print(\"🔍 LOADING DATA\")\n print(f\"{'='*80}\")\n \n metadata_path = os.path.join(CONFIG['dataset_dir'], 'metadata.csv')\n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n if not os.path.exists(metadata_path):\n print(f\"\\n❌ ERROR: {metadata_path} not found!\")\n print(f\" Run Step 1 (preprocessing) first!\")\n return\n \n metadata_df = pd.read_csv(metadata_path)\n print(f\"\\n ✓ Loaded metadata: {len(metadata_df)} patients\")\n \n # Load first fracture and normal case\n fracture_patient = metadata_df[metadata_df['has_fracture']==1].iloc[0]['patient_id']\n normal_patient = metadata_df[metadata_df['has_fracture']==0].iloc[0]['patient_id']\n \n fracture_vol = np.load(os.path.join(volumes_dir, f\"{fracture_patient}.npy\"))\n normal_vol = np.load(os.path.join(volumes_dir, f\"{normal_patient}.npy\"))\n \n print(f\"\\n ✓ Loaded volumes:\")\n print(f\" Fracture: {fracture_patient}, shape={fracture_vol.shape}\")\n print(f\" Normal: {normal_patient}, shape={normal_vol.shape}\")\n \n # ========================================================================\n # RUN DEMONSTRATIONS\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(\"🎨 CREATING AUGMENTATION DEMONSTRATIONS\")\n print(f\"{'='*80}\\n\")\n \n augmenter = EnhancedAugmentation3D()\n \n # Fracture case\n print(\" Processing fracture case...\")\n fracture_dir = os.path.join(CONFIG['output_dir'], 'fracture_case')\n os.makedirs(fracture_dir, exist_ok=True)\n visualize_comprehensive_demo(fracture_vol, 1, augmenter, fracture_dir)\n \n # Normal case\n print(\"\\n Processing normal case...\")\n normal_dir = os.path.join(CONFIG['output_dir'], 'normal_case')\n os.makedirs(normal_dir, exist_ok=True)\n visualize_comprehensive_demo(normal_vol, 0, augmenter, normal_dir)\n \n # ========================================================================\n # MIXUP DEMONSTRATION (if enabled)\n # ========================================================================\n \n if CONFIG['use_mixup']:\n print(f\"\\n Creating MixUp demonstration...\")\n \n mixed_vol, mixed_label = augmenter.mixup(fracture_vol, normal_vol, 1, 0, alpha=0.2)\n \n fig, axes = plt.subplots(1, 4, figsize=(16, 4))\n slice_idx = fracture_vol.shape[0] // 2\n \n axes[0].imshow(fracture_vol[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n axes[0].set_title('Fracture (Label=1)', fontweight='bold', color='red')\n axes[0].axis('off')\n \n axes[1].imshow(normal_vol[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n axes[1].set_title('Normal (Label=0)', fontweight='bold', color='green')\n axes[1].axis('off')\n \n axes[2].imshow(mixed_vol[slice_idx, :, :, 0], cmap='gray', vmin=0, vmax=1)\n axes[2].set_title(f'MixUp (Label={mixed_label:.2f})', fontweight='bold', color='blue')\n axes[2].axis('off')\n \n # Difference map\n diff = np.abs(mixed_vol[slice_idx, :, :, 0] - \n (fracture_vol[slice_idx, :, :, 0] + normal_vol[slice_idx, :, :, 0])/2)\n axes[3].imshow(diff, cmap='hot')\n axes[3].set_title('Difference Map', fontweight='bold')\n axes[3].axis('off')\n \n plt.suptitle('MixUp Augmentation - Better Generalization', \n fontsize=14, fontweight='bold')\n plt.tight_layout()\n save_path = os.path.join(CONFIG['output_dir'], 'demo5_mixup.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n print(f\" ✓ Saved: demo5_mixup.png\")\n \n # ========================================================================\n # FINAL SUMMARY\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(\"✅ AUGMENTATION ANALYSIS COMPLETE!\")\n print(f\"{'='*80}\")\n \n print(f\"\\n📁 Generated files in {CONFIG['output_dir']}:\")\n print(f\" 📂 fracture_case/\")\n print(f\" • demo1_multichannel.png - All 3 channels + augmentation types\")\n print(f\" • demo2_elastic.png - Elastic deformation showcase\")\n print(f\" • demo3_diversity.png - 15 random augmentation samples\")\n print(f\" • demo4_statistics.png - Statistical analysis\")\n print(f\" 📂 normal_case/\")\n print(f\" • (same 4 demos for normal case)\")\n if CONFIG['use_mixup']:\n print(f\" • demo5_mixup.png - MixUp demonstration\")\n \n print(f\"\\n🎯 Key Features Demonstrated:\")\n print(f\" ✅ Multi-channel support (bone/soft tissue/wide)\")\n print(f\" ✅ Elastic deformation (medical imaging essential)\")\n print(f\" ✅ Spatial augmentations (flip, rotate, zoom, shift)\")\n print(f\" ✅ Intensity augmentations (brightness, contrast, gamma)\")\n print(f\" ✅ Noise & blur\")\n if CONFIG['use_mixup']:\n print(f\" ✅ MixUp for better generalization\")\n \n print(f\"\\n📊 Expected Impact:\")\n print(f\" • +10-15% accuracy improvement\")\n print(f\" • Better generalization\")\n print(f\" • Reduced overfitting\")\n print(f\" • More robust predictions\")\n \n print(f\"\\n💡 Next Steps:\")\n print(f\" 1. Review visualizations to understand augmentation effects\")\n print(f\" 2. Integrate EnhancedAugmentation3D into training pipeline\")\n print(f\" 3. Use apply_train_augmentations() in DataLoader\")\n print(f\" 4. Enable MixUp for final 10-20 epochs\")\n print(f\" 5. Use TTA (Test-Time Augmentation) for inference\")\n \n print(f\"\\n🔧 Usage Example:\")\n print(f\"\"\"\n # In your training DataLoader:\n augmenter = EnhancedAugmentation3D()\n \n def __getitem__(self, idx):\n volume = np.load(f\"patient_{{idx}}.npy\")\n label = self.labels[idx]\n \n # Apply augmentation\n if self.training:\n volume = augmenter.apply_train_augmentations(volume)\n \n # Optional: MixUp (15% probability)\n if np.random.random() < 0.15:\n idx2 = np.random.randint(len(self))\n volume2 = np.load(f\"patient_{{idx2}}.npy\")\n label2 = self.labels[idx2]\n volume, label = augmenter.mixup(volume, volume2, label, label2)\n \n return torch.from_numpy(volume), torch.tensor(label)\n \"\"\")\n \n print(f\"\\n{'='*80}\")\n print(\"🎉 READY FOR TRAINING!\")\n print(f\"{'='*80}\")\n\n\nif __name__ == \"__main__\":\n try:\n main()\n except KeyboardInterrupt:\n print(f\"\\n\\n⚠️ Interrupted by user\")\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()\n \n print(f\"\\n💡 Troubleshooting:\")\n print(f\" 1. Make sure Step 1 preprocessing completed successfully\")\n print(f\" 2. Check dataset path: {CONFIG['dataset_dir']}\")\n print(f\" 3. Verify volumes directory exists with .npy files\")\n print(f\" 4. Ensure sufficient RAM (~8GB recommended)\")\n print(f\" 5. Check metadata.csv is valid\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T09:47:13.05616Z","iopub.execute_input":"2026-01-22T09:47:13.056729Z","iopub.status.idle":"2026-01-22T09:55:27.456978Z","shell.execute_reply.started":"2026-01-22T09:47:13.056671Z","shell.execute_reply":"2026-01-22T09:55:27.456048Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nSTEP 3: PRODUCTION TRAINING PIPELINE v3.0\n==========================================\nComplete training system with state-of-the-art techniques\n\nNEW FEATURES:\n✅ 3D ResNet architecture for multi-channel volumes\n✅ 5-fold cross-validation\n✅ Mixed precision training (faster + less memory)\n✅ Cosine annealing with warm restarts\n✅ Gradient accumulation for larger effective batch size\n✅ Early stopping & model checkpointing\n✅ Advanced metrics (AUC, F1, per-vertebra accuracy)\n✅ TensorBoard logging\n✅ Ensemble predictions\n✅ Test-Time Augmentation (TTA)\n\nExpected: 75-85% AUC on validation\nTraining time: ~2-4 hours for 50 epochs (5 folds)\n\"\"\"\n\nimport os\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nfrom torch.cuda.amp import autocast, GradScaler\n\nfrom sklearn.metrics import roc_auc_score, accuracy_score, f1_score, confusion_matrix\nimport seaborn as sns\n\nprint(\"=\"*80)\nprint(\"🚀 STEP 3: PRODUCTION TRAINING PIPELINE v3.0\")\nprint(\"=\"*80)\n\n# ============================================================================\n# CONFIGURATION\n# ============================================================================\n\nCONFIG = {\n # Paths\n 'dataset_dir': '/kaggle/working/enhanced_dataset_v3',\n 'output_dir': '/kaggle/working/training_output',\n \n # Training\n 'n_folds': 5,\n 'epochs': 50,\n 'batch_size': 4, # Small due to 3D volumes\n 'accumulation_steps': 4, # Effective batch size = 16\n 'learning_rate': 1e-3,\n 'weight_decay': 1e-4,\n \n # Architecture\n 'in_channels': 3, # Multi-channel input!\n 'base_filters': 16,\n 'dropout': 0.3,\n \n # Augmentation\n 'use_augmentation': True,\n 'use_mixup': True,\n 'mixup_prob': 0.15,\n \n # Advanced training\n 'use_mixed_precision': True,\n 'gradient_clip': 1.0,\n 'early_stopping_patience': 10,\n \n # Evaluation\n 'use_tta': True,\n 'tta_iterations': 5,\n \n # System\n 'num_workers': 2,\n 'random_seed': 42,\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\nprint(f\"\\n📋 Configuration:\")\nprint(f\" • Folds: {CONFIG['n_folds']}\")\nprint(f\" • Epochs: {CONFIG['epochs']}\")\nprint(f\" • Effective batch size: {CONFIG['batch_size'] * CONFIG['accumulation_steps']}\")\nprint(f\" • Learning rate: {CONFIG['learning_rate']}\")\nprint(f\" • Mixed precision: {CONFIG['use_mixed_precision']}\")\nprint(f\" • TTA: {CONFIG['use_tta']}\")\n\n# Set seeds\ntorch.manual_seed(CONFIG['random_seed'])\nnp.random.seed(CONFIG['random_seed'])\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\" • Device: {device}\")\n\n# ============================================================================\n# AUGMENTATION (from Step 2)\n# ============================================================================\n\nfrom scipy.ndimage import rotate, zoom, gaussian_filter, shift, map_coordinates\n\nclass Augmentation3D:\n \"\"\"Simplified augmentation for training\"\"\"\n \n def __init__(self):\n self.config = {\n 'flip_prob': 0.5,\n 'rotation': (-15, 15),\n 'rotation_prob': 0.4,\n 'zoom': (0.85, 1.15),\n 'zoom_prob': 0.3,\n 'elastic_alpha': (0, 30),\n 'elastic_prob': 0.3,\n 'brightness': (0.85, 1.15),\n 'brightness_prob': 0.4,\n 'contrast': (0.85, 1.15),\n 'contrast_prob': 0.4,\n }\n \n def horizontal_flip(self, vol):\n return np.flip(vol, axis=2).copy()\n \n def random_rotation(self, vol, angle=None):\n if angle is None:\n angle = np.random.uniform(*self.config['rotation'])\n \n if vol.ndim == 4:\n rotated = np.zeros_like(vol)\n for ch in range(vol.shape[3]):\n rotated[..., ch] = rotate(vol[..., ch], angle, axes=(1,2), reshape=False, order=1)\n return rotated.astype(np.float32)\n return rotate(vol, angle, axes=(1,2), reshape=False, order=1).astype(np.float32)\n \n def random_zoom(self, vol, scale=None):\n if scale is None:\n scale = np.random.uniform(*self.config['zoom'])\n \n orig_shape = vol.shape\n \n if vol.ndim == 4:\n factors = [scale, scale, scale, 1.0]\n zoomed = zoom(vol, factors, order=1)\n \n if scale > 1.0:\n start = [(zoomed.shape[i] - orig_shape[i])//2 for i in range(3)]\n result = zoomed[start[0]:start[0]+orig_shape[0],\n start[1]:start[1]+orig_shape[1],\n start[2]:start[2]+orig_shape[2], :]\n else:\n result = np.zeros(orig_shape, dtype=vol.dtype)\n start = [(orig_shape[i] - zoomed.shape[i])//2 for i in range(3)]\n result[start[0]:start[0]+zoomed.shape[0],\n start[1]:start[1]+zoomed.shape[1],\n start[2]:start[2]+zoomed.shape[2], :] = zoomed\n return result.astype(np.float32)\n \n zoomed = zoom(vol, scale, order=1)\n if scale > 1.0:\n start = [(zoomed.shape[i] - orig_shape[i])//2 for i in range(3)]\n return zoomed[start[0]:start[0]+orig_shape[0],\n start[1]:start[1]+orig_shape[1],\n start[2]:start[2]+orig_shape[2]].astype(np.float32)\n else:\n result = np.zeros(orig_shape, dtype=vol.dtype)\n start = [(orig_shape[i] - zoomed.shape[i])//2 for i in range(3)]\n result[start[0]:start[0]+zoomed.shape[0],\n start[1]:start[1]+zoomed.shape[1],\n start[2]:start[2]+zoomed.shape[2]] = zoomed\n return result.astype(np.float32)\n \n def elastic_deformation(self, vol, alpha=None):\n if alpha is None:\n alpha = np.random.uniform(*self.config['elastic_alpha'])\n if alpha == 0:\n return vol\n \n shape = vol.shape[:3]\n sigma = 3\n \n dx = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n dy = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n dz = gaussian_filter((np.random.rand(*shape)*2-1), sigma) * alpha\n \n d, h, w = shape\n d_coords, h_coords, w_coords = np.meshgrid(\n np.arange(d), np.arange(h), np.arange(w), indexing='ij'\n )\n \n indices = (\n np.reshape(d_coords + dx, (-1, 1)),\n np.reshape(h_coords + dy, (-1, 1)),\n np.reshape(w_coords + dz, (-1, 1))\n )\n \n if vol.ndim == 4:\n deformed = np.zeros_like(vol)\n for ch in range(vol.shape[3]):\n deformed_ch = map_coordinates(vol[..., ch], indices, order=1, mode='nearest')\n deformed[..., ch] = deformed_ch.reshape(shape)\n else:\n deformed = map_coordinates(vol, indices, order=1, mode='nearest')\n deformed = deformed.reshape(shape)\n \n return deformed.astype(np.float32)\n \n def apply(self, vol):\n \"\"\"Apply random augmentations\"\"\"\n aug = vol.copy()\n \n if np.random.random() < self.config['flip_prob']:\n aug = self.horizontal_flip(aug)\n if np.random.random() < self.config['rotation_prob']:\n aug = self.random_rotation(aug)\n if np.random.random() < self.config['zoom_prob']:\n aug = self.random_zoom(aug)\n if np.random.random() < self.config['elastic_prob']:\n aug = self.elastic_deformation(aug)\n if np.random.random() < self.config['brightness_prob']:\n factor = np.random.uniform(*self.config['brightness'])\n aug = np.clip(aug * factor, 0, 1).astype(np.float32)\n if np.random.random() < self.config['contrast_prob']:\n factor = np.random.uniform(*self.config['contrast'])\n mean = aug.mean()\n aug = np.clip((aug - mean) * factor + mean, 0, 1).astype(np.float32)\n \n return aug\n\n\n# ============================================================================\n# DATASET\n# ============================================================================\n\nclass SpineDataset(Dataset):\n \"\"\"Dataset for cervical spine fracture detection\"\"\"\n \n def __init__(self, patient_ids, metadata_df, volumes_dir, training=True, augmenter=None):\n self.patient_ids = patient_ids\n self.metadata_df = metadata_df.set_index('patient_id')\n self.volumes_dir = volumes_dir\n self.training = training\n self.augmenter = augmenter\n \n def __len__(self):\n return len(self.patient_ids)\n \n def __getitem__(self, idx):\n patient_id = self.patient_ids[idx]\n \n # Load volume (D, H, W, 3)\n volume = np.load(os.path.join(self.volumes_dir, f\"{patient_id}.npy\"))\n \n # Get label\n label = self.metadata_df.loc[patient_id, 'has_fracture']\n \n # Augmentation\n if self.training and self.augmenter is not None:\n volume = self.augmenter.apply(volume)\n \n # Convert to PyTorch: (D,H,W,3) -> (3,D,H,W)\n volume = torch.from_numpy(volume).permute(3, 0, 1, 2).float()\n label = torch.tensor(label, dtype=torch.float32)\n \n return volume, label\n\n\n# ============================================================================\n# 3D RESNET ARCHITECTURE\n# ============================================================================\n\nclass ResidualBlock3D(nn.Module):\n \"\"\"3D Residual block\"\"\"\n \n def __init__(self, in_channels, out_channels, stride=1):\n super().__init__()\n \n self.conv1 = nn.Conv3d(in_channels, out_channels, kernel_size=3, \n stride=stride, padding=1, bias=False)\n self.bn1 = nn.BatchNorm3d(out_channels)\n self.relu = nn.ReLU(inplace=True)\n \n self.conv2 = nn.Conv3d(out_channels, out_channels, kernel_size=3,\n stride=1, padding=1, bias=False)\n self.bn2 = nn.BatchNorm3d(out_channels)\n \n # Shortcut connection\n if stride != 1 or in_channels != out_channels:\n self.shortcut = nn.Sequential(\n nn.Conv3d(in_channels, out_channels, kernel_size=1, \n stride=stride, bias=False),\n nn.BatchNorm3d(out_channels)\n )\n else:\n self.shortcut = nn.Identity()\n \n def forward(self, x):\n residual = x\n \n out = self.conv1(x)\n out = self.bn1(out)\n out = self.relu(out)\n \n out = self.conv2(out)\n out = self.bn2(out)\n \n out += self.shortcut(residual)\n out = self.relu(out)\n \n return out\n\n\nclass ResNet3D(nn.Module):\n \"\"\"3D ResNet for multi-channel medical imaging\"\"\"\n \n def __init__(self, in_channels=3, base_filters=16, dropout=0.3):\n super().__init__()\n \n # Initial convolution\n self.conv1 = nn.Conv3d(in_channels, base_filters, kernel_size=7, \n stride=2, padding=3, bias=False)\n self.bn1 = nn.BatchNorm3d(base_filters)\n self.relu = nn.ReLU(inplace=True)\n self.maxpool = nn.MaxPool3d(kernel_size=3, stride=2, padding=1)\n \n # Residual blocks\n self.layer1 = self._make_layer(base_filters, base_filters*2, blocks=2, stride=1)\n self.layer2 = self._make_layer(base_filters*2, base_filters*4, blocks=2, stride=2)\n self.layer3 = self._make_layer(base_filters*4, base_filters*8, blocks=2, stride=2)\n \n # Global pooling and classifier\n self.avgpool = nn.AdaptiveAvgPool3d(1)\n self.dropout = nn.Dropout(dropout)\n self.fc = nn.Linear(base_filters*8, 1)\n \n def _make_layer(self, in_channels, out_channels, blocks, stride):\n layers = []\n layers.append(ResidualBlock3D(in_channels, out_channels, stride))\n \n for _ in range(1, blocks):\n layers.append(ResidualBlock3D(out_channels, out_channels, stride=1))\n \n return nn.Sequential(*layers)\n \n def forward(self, x):\n # Initial conv\n x = self.conv1(x)\n x = self.bn1(x)\n x = self.relu(x)\n x = self.maxpool(x)\n \n # Residual blocks\n x = self.layer1(x)\n x = self.layer2(x)\n x = self.layer3(x)\n \n # Classification\n x = self.avgpool(x)\n x = x.view(x.size(0), -1)\n x = self.dropout(x)\n x = self.fc(x)\n \n return x.squeeze(-1)\n\n\n# ============================================================================\n# TRAINING FUNCTIONS\n# ============================================================================\n\nclass FocalLoss(nn.Module):\n \"\"\"Focal loss for class imbalance\"\"\"\n \n def __init__(self, alpha=0.25, gamma=2.0):\n super().__init__()\n self.alpha = alpha\n self.gamma = gamma\n \n def forward(self, inputs, targets):\n bce_loss = nn.functional.binary_cross_entropy_with_logits(\n inputs, targets, reduction='none'\n )\n \n pt = torch.exp(-bce_loss)\n focal_loss = self.alpha * (1 - pt) ** self.gamma * bce_loss\n \n return focal_loss.mean()\n\n\ndef train_epoch(model, loader, criterion, optimizer, scaler, device, accumulation_steps):\n \"\"\"Train for one epoch\"\"\"\n model.train()\n running_loss = 0.0\n all_preds = []\n all_labels = []\n \n optimizer.zero_grad()\n \n for batch_idx, (volumes, labels) in enumerate(loader):\n volumes, labels = volumes.to(device), labels.to(device)\n \n # Mixed precision training\n with autocast(enabled=CONFIG['use_mixed_precision']):\n outputs = model(volumes)\n loss = criterion(outputs, labels)\n loss = loss / accumulation_steps\n \n # Backward pass\n scaler.scale(loss).backward()\n \n # Gradient accumulation\n if (batch_idx + 1) % accumulation_steps == 0:\n scaler.unscale_(optimizer)\n torch.nn.utils.clip_grad_norm_(model.parameters(), CONFIG['gradient_clip'])\n scaler.step(optimizer)\n scaler.update()\n optimizer.zero_grad()\n \n running_loss += loss.item() * accumulation_steps\n \n # Store predictions\n with torch.no_grad():\n preds = torch.sigmoid(outputs).cpu().numpy()\n all_preds.extend(preds)\n all_labels.extend(labels.cpu().numpy())\n \n epoch_loss = running_loss / len(loader)\n epoch_auc = roc_auc_score(all_labels, all_preds) if len(np.unique(all_labels)) > 1 else 0.5\n \n return epoch_loss, epoch_auc\n\n\ndef validate_epoch(model, loader, criterion, device):\n \"\"\"Validate for one epoch\"\"\"\n model.eval()\n running_loss = 0.0\n all_preds = []\n all_labels = []\n \n with torch.no_grad():\n for volumes, labels in loader:\n volumes, labels = volumes.to(device), labels.to(device)\n \n with autocast(enabled=CONFIG['use_mixed_precision']):\n outputs = model(volumes)\n loss = criterion(outputs, labels)\n \n running_loss += loss.item()\n \n preds = torch.sigmoid(outputs).cpu().numpy()\n all_preds.extend(preds)\n all_labels.extend(labels.cpu().numpy())\n \n epoch_loss = running_loss / len(loader)\n \n all_preds = np.array(all_preds)\n all_labels = np.array(all_labels)\n \n metrics = {\n 'loss': epoch_loss,\n 'auc': roc_auc_score(all_labels, all_preds) if len(np.unique(all_labels)) > 1 else 0.5,\n 'accuracy': accuracy_score(all_labels, (all_preds > 0.5).astype(int)),\n 'f1': f1_score(all_labels, (all_preds > 0.5).astype(int), zero_division=0)\n }\n \n return metrics\n\n\n# ============================================================================\n# MAIN TRAINING LOOP\n# ============================================================================\n\ndef train_fold(fold_idx, train_patients, val_patients, metadata_df, volumes_dir):\n \"\"\"Train a single fold\"\"\"\n \n print(f\"\\n{'='*80}\")\n print(f\"📊 TRAINING FOLD {fold_idx + 1}/{CONFIG['n_folds']}\")\n print(f\"{'='*80}\")\n print(f\" Train: {len(train_patients)} patients\")\n print(f\" Val: {len(val_patients)} patients\")\n \n # Create datasets\n augmenter = Augmentation3D() if CONFIG['use_augmentation'] else None\n \n train_dataset = SpineDataset(train_patients, metadata_df, volumes_dir, \n training=True, augmenter=augmenter)\n val_dataset = SpineDataset(val_patients, metadata_df, volumes_dir, \n training=False, augmenter=None)\n \n train_loader = DataLoader(train_dataset, batch_size=CONFIG['batch_size'],\n shuffle=True, num_workers=CONFIG['num_workers'])\n val_loader = DataLoader(val_dataset, batch_size=CONFIG['batch_size'],\n shuffle=False, num_workers=CONFIG['num_workers'])\n \n # Model\n model = ResNet3D(\n in_channels=CONFIG['in_channels'],\n base_filters=CONFIG['base_filters'],\n dropout=CONFIG['dropout']\n ).to(device)\n \n print(f\"\\n Model parameters: {sum(p.numel() for p in model.parameters()):,}\")\n \n # Loss and optimizer\n criterion = FocalLoss(alpha=0.25, gamma=2.0)\n optimizer = optim.AdamW(model.parameters(), lr=CONFIG['learning_rate'],\n weight_decay=CONFIG['weight_decay'])\n \n # Learning rate scheduler\n scheduler = optim.lr_scheduler.CosineAnnealingWarmRestarts(\n optimizer, T_0=10, T_mult=2, eta_min=1e-6\n )\n \n scaler = GradScaler(enabled=CONFIG['use_mixed_precision'])\n \n # Training history\n history = {\n 'train_loss': [],\n 'train_auc': [],\n 'val_loss': [],\n 'val_auc': [],\n 'val_accuracy': [],\n 'val_f1': [],\n }\n \n best_auc = 0.0\n patience_counter = 0\n \n # Training loop\n for epoch in range(CONFIG['epochs']):\n print(f\"\\n Epoch {epoch+1}/{CONFIG['epochs']}\")\n \n # Train\n train_loss, train_auc = train_epoch(\n model, train_loader, criterion, optimizer, scaler, device,\n CONFIG['accumulation_steps']\n )\n \n # Validate\n val_metrics = validate_epoch(model, val_loader, criterion, device)\n \n # Scheduler step\n scheduler.step()\n \n # Save history\n history['train_loss'].append(train_loss)\n history['train_auc'].append(train_auc)\n history['val_loss'].append(val_metrics['loss'])\n history['val_auc'].append(val_metrics['auc'])\n history['val_accuracy'].append(val_metrics['accuracy'])\n history['val_f1'].append(val_metrics['f1'])\n \n print(f\" Train - Loss: {train_loss:.4f}, AUC: {train_auc:.4f}\")\n print(f\" Val - Loss: {val_metrics['loss']:.4f}, AUC: {val_metrics['auc']:.4f}, \"\n f\"Acc: {val_metrics['accuracy']:.4f}, F1: {val_metrics['f1']:.4f}\")\n \n # Save best model\n if val_metrics['auc'] > best_auc:\n best_auc = val_metrics['auc']\n patience_counter = 0\n \n checkpoint_path = os.path.join(CONFIG['output_dir'], f'fold{fold_idx}_best.pth')\n torch.save({\n 'epoch': epoch,\n 'model_state_dict': model.state_dict(),\n 'optimizer_state_dict': optimizer.state_dict(),\n 'auc': best_auc,\n }, checkpoint_path)\n print(f\" 💾 Saved best model (AUC: {best_auc:.4f})\")\n else:\n patience_counter += 1\n \n # Early stopping\n if patience_counter >= CONFIG['early_stopping_patience']:\n print(f\"\\n ⚠️ Early stopping triggered (patience={CONFIG['early_stopping_patience']})\")\n break\n \n print(f\"\\n ✅ Fold {fold_idx + 1} complete! Best AUC: {best_auc:.4f}\")\n \n return history, best_auc, model\n\n\n# ============================================================================\n# VISUALIZATION\n# ============================================================================\n\ndef plot_training_history(histories, save_path):\n \"\"\"Plot training curves for all folds\"\"\"\n \n fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n \n for fold_idx, history in enumerate(histories):\n epochs = range(1, len(history['train_loss']) + 1)\n \n # Loss\n axes[0, 0].plot(epochs, history['train_loss'], label=f'Fold {fold_idx+1} Train', alpha=0.7)\n axes[0, 0].plot(epochs, history['val_loss'], label=f'Fold {fold_idx+1} Val', alpha=0.7, linestyle='--')\n \n # AUC\n axes[0, 1].plot(epochs, history['train_auc'], label=f'Fold {fold_idx+1} Train', alpha=0.7)\n axes[0, 1].plot(epochs, history['val_auc'], label=f'Fold {fold_idx+1} Val', alpha=0.7, linestyle='--')\n \n # Accuracy\n axes[1, 0].plot(epochs, history['val_accuracy'], label=f'Fold {fold_idx+1}', alpha=0.7)\n \n # F1\n axes[1, 1].plot(epochs, history['val_f1'], label=f'Fold {fold_idx+1}', alpha=0.7)\n \n axes[0, 0].set_title('Loss', fontweight='bold', fontsize=12)\n axes[0, 0].set_xlabel('Epoch')\n axes[0, 0].set_ylabel('Loss')\n axes[0, 0].legend()\n axes[0, 0].grid(True, alpha=0.3)\n \n axes[0, 1].set_title('AUC', fontweight='bold', fontsize=12)\n axes[0, 1].set_xlabel('Epoch')\n axes[0, 1].set_ylabel('AUC')\n axes[0, 1].legend()\n axes[0, 1].grid(True, alpha=0.3)\n \n axes[1, 0].set_title('Validation Accuracy', fontweight='bold', fontsize=12)\n axes[1, 0].set_xlabel('Epoch')\n axes[1, 0].set_ylabel('Accuracy')\n axes[1, 0].legend()\n axes[1, 0].grid(True, alpha=0.3)\n \n axes[1, 1].set_title('Validation F1 Score', fontweight='bold', fontsize=12)\n axes[1, 1].set_xlabel('Epoch')\n axes[1, 1].set_ylabel('F1 Score')\n axes[1, 1].legend()\n axes[1, 1].grid(True, alpha=0.3)\n \n plt.suptitle(f'{CONFIG[\"n_folds\"]}-Fold Cross-Validation Training', \n fontsize=14, fontweight='bold')\n plt.tight_layout()\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved training curves: {save_path}\")\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n print(f\"\\n{'='*80}\")\n print(\"🔍 LOADING DATA\")\n print(f\"{'='*80}\")\n \n # Load metadata and CV splits\n metadata_path = os.path.join(CONFIG['dataset_dir'], 'metadata.csv')\n cv_path = os.path.join(CONFIG['dataset_dir'], 'cv_splits.json')\n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n if not os.path.exists(metadata_path):\n print(f\"\\n❌ ERROR: {metadata_path} not found!\")\n return\n \n if not os.path.exists(cv_path):\n print(f\"\\n❌ ERROR: {cv_path} not found! Run Step 1 first.\")\n return\n \n metadata_df = pd.read_csv(metadata_path)\n print(f\"\\n ✓ Loaded metadata: {len(metadata_df)} patients\")\n \n with open(cv_path, 'r') as f:\n cv_splits = json.load(f)\n print(f\" ✓ Loaded {len(cv_splits)} CV splits\")\n \n # Train all folds\n all_histories = []\n fold_aucs = []\n \n for fold_idx, split in enumerate(cv_splits):\n history, best_auc, model = train_fold(\n fold_idx,\n split['train'],\n split['val'],\n metadata_df,\n volumes_dir\n )\n \n all_histories.append(history)\n fold_aucs.append(best_auc)\n \n # ========================================================================\n # FINAL SUMMARY\n # ========================================================================\n \n print(f\"\\n{'='*80}\")\n print(\"✅ TRAINING COMPLETE!\")\n print(f\"{'='*80}\")\n \n print(f\"\\n📊 Cross-Validation Results:\")\n for fold_idx, auc in enumerate(fold_aucs):\n print(f\" Fold {fold_idx+1}: AUC = {auc:.4f}\")\n \n mean_auc = np.mean(fold_aucs)\n std_auc = np.std(fold_aucs)\n print(f\"\\n Mean AUC: {mean_auc:.4f} ± {std_auc:.4f}\")\n \n # Plot training curves\n curves_path = os.path.join(CONFIG['output_dir'], 'training_curves.png')\n plot_training_history(all_histories, curves_path)\n \n # Save results\n results = {\n 'fold_aucs': fold_aucs,\n 'mean_auc': float(mean_auc),\n 'std_auc': float(std_auc),\n 'config': CONFIG\n }\n \n results_path = os.path.join(CONFIG['output_dir'], 'results.json')\n with open(results_path, 'w') as f:\n json.dump(results, f, indent=2)\n \n print(f\"\\n📁 Output:\")\n print(f\" • Model checkpoints: {CONFIG['output_dir']}/fold*_best.pth\")\n print(f\" • Training curves: {curves_path}\")\n print(f\" • Results: {results_path}\")\n \n print(f\"\\n💡 Next Steps:\")\n print(f\" 1. Review training curves for overfitting\")\n print(f\" 2. If AUC < 0.70, try:\")\n print(f\" - Increase epochs to 100\")\n print(f\" - Use more augmentation\")\n print(f\" - Reduce learning rate to 5e-4\")\n print(f\" - Increase base_filters to 32\")\n print(f\" 3. If AUC > 0.75:\")\n print(f\" - Proceed to Step 4 (ensemble predictions)\")\n print(f\" - Enable Test-Time Augmentation\")\n print(f\" 4. For production:\")\n print(f\" - Train on all data (no CV split)\")\n print(f\" - Use best hyperparameters from CV\")\n \n print(f\"\\n🔧 Model Usage Example:\")\n print(f\"\"\"\n # Load best model from a fold\n checkpoint = torch.load('training_output/fold0_best.pth')\n model = ResNet3D(in_channels=3, base_filters=16)\n model.load_state_dict(checkpoint['model_state_dict'])\n model.eval()\n \n # Inference on new patient\n volume = np.load('new_patient.npy') # Shape: (64, 224, 224, 3)\n volume_tensor = torch.from_numpy(volume).permute(3,0,1,2).unsqueeze(0)\n \n with torch.no_grad():\n logit = model(volume_tensor)\n probability = torch.sigmoid(logit).item()\n \n print(f\"Fracture probability: {{probability:.2%}}\")\n \"\"\")\n \n print(f\"\\n{'='*80}\")\n print(\"🎉 TRAINING PIPELINE READY!\")\n print(f\"{'='*80}\")\n\n\nif __name__ == \"__main__\":\n try:\n main()\n except KeyboardInterrupt:\n print(f\"\\n\\n⚠️ Training interrupted by user\")\n print(f\" Partial progress saved in checkpoints\")\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()\n \n print(f\"\\n💡 Troubleshooting:\")\n print(f\" 1. Check CUDA availability: torch.cuda.is_available()\")\n print(f\" 2. Reduce batch_size if OOM error (try batch_size=2)\")\n print(f\" 3. Disable mixed precision: use_mixed_precision=False\")\n print(f\" 4. Check dataset paths are correct\")\n print(f\" 5. Verify all .npy volumes are valid\")\n print(f\" 6. For CPU training: expect ~10x slower\")\n print(f\"\\n Common fixes:\")\n print(f\" • OOM → Reduce batch_size to 2 or 1\")\n print(f\" • Slow training → Enable mixed precision\")\n print(f\" • NaN loss → Reduce learning rate to 5e-4\")\n print(f\" • Poor AUC → Increase epochs or augmentation\") ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T09:55:27.459032Z","iopub.execute_input":"2026-01-22T09:55:27.459358Z","iopub.status.idle":"2026-01-22T13:12:44.309038Z","shell.execute_reply.started":"2026-01-22T09:55:27.459336Z","shell.execute_reply":"2026-01-22T13:12:44.30816Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install pylibjpeg pylibjpeg-libjpeg","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T13:25:41.351665Z","iopub.execute_input":"2026-01-22T13:25:41.352481Z","iopub.status.idle":"2026-01-22T13:25:48.60569Z","shell.execute_reply.started":"2026-01-22T13:25:41.352452Z","shell.execute_reply":"2026-01-22T13:25:48.604878Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nSTEP 4: ENSEMBLE INFERENCE & SUBMISSION v4.0 (ENHANCED)\n========================================================\nIMPROVEMENTS:\n✓ Fixed DICOM decompression issues\n✓ Disabled problematic calibration\n✓ Added weighted ensemble by validation AUC\n✓ Better error handling\n✓ Optimized for speed\n\"\"\"\n\nimport os\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torch.nn as nn\nfrom torch.cuda.amp import autocast\n\nprint(\"=\"*80)\nprint(\"🚀 STEP 4: ENSEMBLE INFERENCE & SUBMISSION v4.0 (ENHANCED)\")\nprint(\"=\"*80)\n\n# ============================================================================\n# CONFIGURATION - OPTIMIZED SETTINGS\n# ============================================================================\n\nCONFIG = {\n # Paths\n 'models_dir': '/kaggle/working/training_output',\n 'test_dir': '/kaggle/input/rsna-2022-cervical-spine-fracture-detection/test_images',\n 'output_dir': '/kaggle/working/inference_output',\n \n # Inference\n 'n_folds': 5,\n 'use_tta': True,\n 'tta_iterations': 4, # Reduced from 5 for speed\n 'batch_size': 4,\n \n # Model architecture (must match training!)\n 'in_channels': 3,\n 'base_filters': 16,\n 'dropout': 0.3,\n \n # Calibration - DISABLED (was causing centered predictions)\n 'calibrate_predictions': False, # CHANGED: Was True\n 'temperature': 1.0, # CHANGED: Was 1.5\n \n # NEW: Weighted ensemble\n 'use_weighted_ensemble': True, # Weight by validation AUC\n 'min_auc_threshold': 0.60, # Exclude folds below this\n \n # System\n 'use_mixed_precision': True,\n 'num_workers': 2,\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"\\n📋 Configuration:\")\nprint(f\" • Device: {device}\")\nprint(f\" • Ensemble size: {CONFIG['n_folds']} folds\")\nprint(f\" • TTA: {CONFIG['use_tta']} ({CONFIG['tta_iterations']} iterations)\")\nprint(f\" • Calibration: {CONFIG['calibrate_predictions']} (DISABLED)\")\nprint(f\" • Weighted ensemble: {CONFIG['use_weighted_ensemble']}\")\n\n# ============================================================================\n# MODEL ARCHITECTURE (same as Step 3)\n# ============================================================================\n\nclass ResidualBlock3D(nn.Module):\n def __init__(self, in_channels, out_channels, stride=1):\n super().__init__()\n self.conv1 = nn.Conv3d(in_channels, out_channels, kernel_size=3, \n stride=stride, padding=1, bias=False)\n self.bn1 = nn.BatchNorm3d(out_channels)\n self.relu = nn.ReLU(inplace=True)\n self.conv2 = nn.Conv3d(out_channels, out_channels, kernel_size=3,\n stride=1, padding=1, bias=False)\n self.bn2 = nn.BatchNorm3d(out_channels)\n \n if stride != 1 or in_channels != out_channels:\n self.shortcut = nn.Sequential(\n nn.Conv3d(in_channels, out_channels, kernel_size=1, \n stride=stride, bias=False),\n nn.BatchNorm3d(out_channels)\n )\n else:\n self.shortcut = nn.Identity()\n \n def forward(self, x):\n residual = x\n out = self.conv1(x)\n out = self.bn1(out)\n out = self.relu(out)\n out = self.conv2(out)\n out = self.bn2(out)\n out += self.shortcut(residual)\n out = self.relu(out)\n return out\n\n\nclass ResNet3D(nn.Module):\n def __init__(self, in_channels=3, base_filters=16, dropout=0.3):\n super().__init__()\n self.conv1 = nn.Conv3d(in_channels, base_filters, kernel_size=7, \n stride=2, padding=3, bias=False)\n self.bn1 = nn.BatchNorm3d(base_filters)\n self.relu = nn.ReLU(inplace=True)\n self.maxpool = nn.MaxPool3d(kernel_size=3, stride=2, padding=1)\n \n self.layer1 = self._make_layer(base_filters, base_filters*2, blocks=2, stride=1)\n self.layer2 = self._make_layer(base_filters*2, base_filters*4, blocks=2, stride=2)\n self.layer3 = self._make_layer(base_filters*4, base_filters*8, blocks=2, stride=2)\n \n self.avgpool = nn.AdaptiveAvgPool3d(1)\n self.dropout = nn.Dropout(dropout)\n self.fc = nn.Linear(base_filters*8, 1)\n \n def _make_layer(self, in_channels, out_channels, blocks, stride):\n layers = []\n layers.append(ResidualBlock3D(in_channels, out_channels, stride))\n for _ in range(1, blocks):\n layers.append(ResidualBlock3D(out_channels, out_channels, stride=1))\n return nn.Sequential(*layers)\n \n def forward(self, x):\n x = self.conv1(x)\n x = self.bn1(x)\n x = self.relu(x)\n x = self.maxpool(x)\n x = self.layer1(x)\n x = self.layer2(x)\n x = self.layer3(x)\n x = self.avgpool(x)\n x = x.view(x.size(0), -1)\n x = self.dropout(x)\n x = self.fc(x)\n return x.squeeze(-1)\n\n\n# ============================================================================\n# TTA AUGMENTATIONS - OPTIMIZED\n# ============================================================================\n\nclass TTATransform:\n \"\"\"Simple Test-Time Augmentation - Optimized\"\"\"\n \n def __init__(self):\n self.transforms = [\n self.identity,\n self.horizontal_flip,\n self.rotate_5,\n self.brightness_up,\n ]\n \n def identity(self, x):\n return x\n \n def horizontal_flip(self, x):\n return np.flip(x, axis=2).copy()\n \n def rotate_5(self, x):\n from scipy.ndimage import rotate\n if x.ndim == 4:\n rotated = np.zeros_like(x)\n for ch in range(x.shape[3]):\n rotated[..., ch] = rotate(x[..., ch], 5, axes=(1,2), reshape=False, order=1)\n return rotated.astype(np.float32)\n return rotate(x, 5, axes=(1,2), reshape=False, order=1).astype(np.float32)\n \n def brightness_up(self, x):\n return np.clip(x * 1.05, 0, 1).astype(np.float32)\n \n def get_transform(self, idx):\n return self.transforms[idx % len(self.transforms)]\n\n\n# ============================================================================\n# ENHANCED ENSEMBLE PREDICTOR\n# ============================================================================\n\nclass EnsemblePredictor:\n \"\"\"Enhanced 5-fold ensemble with weighted predictions\"\"\"\n \n def __init__(self, models_dir, n_folds=5, use_tta=True, tta_iterations=4):\n self.models_dir = models_dir\n self.n_folds = n_folds\n self.use_tta = use_tta\n self.tta_iterations = tta_iterations\n self.tta_transform = TTATransform()\n \n # Load all fold models\n self.models = []\n self.fold_weights = []\n self.fold_aucs = []\n \n print(f\"\\n🔍 Loading {n_folds} fold models...\")\n \n for fold_idx in range(n_folds):\n checkpoint_path = os.path.join(models_dir, f'fold{fold_idx}_best.pth')\n \n if not os.path.exists(checkpoint_path):\n print(f\" ⚠️ Warning: {checkpoint_path} not found, skipping\")\n continue\n \n model = ResNet3D(\n in_channels=CONFIG['in_channels'],\n base_filters=CONFIG['base_filters'],\n dropout=CONFIG['dropout']\n ).to(device)\n \n checkpoint = torch.load(checkpoint_path, map_location=device, weights_only=False)\n model.load_state_dict(checkpoint['model_state_dict'])\n model.eval()\n \n auc = checkpoint.get('auc', 0.5)\n \n # Filter by minimum AUC threshold\n if auc < CONFIG['min_auc_threshold']:\n print(f\" ⚠️ Skipping fold {fold_idx} (AUC: {auc:.4f} < {CONFIG['min_auc_threshold']})\")\n continue\n \n self.models.append(model)\n self.fold_aucs.append(auc)\n \n # Weight by AUC if enabled\n if CONFIG['use_weighted_ensemble']:\n weight = auc\n else:\n weight = 1.0\n \n self.fold_weights.append(weight)\n print(f\" ✓ Loaded fold {fold_idx} (AUC: {auc:.4f}, Weight: {weight:.4f})\")\n \n # Normalize weights\n if len(self.fold_weights) > 0:\n total_weight = sum(self.fold_weights)\n self.fold_weights = [w / total_weight for w in self.fold_weights]\n \n print(f\"\\n Total models loaded: {len(self.models)}\")\n if CONFIG['use_weighted_ensemble']:\n print(f\" Weighted ensemble enabled (weights: {[f'{w:.3f}' for w in self.fold_weights]})\")\n \n def predict_single(self, volume):\n \"\"\"\n Predict with weighted ensemble + TTA\n \n Args:\n volume: numpy array (D, H, W, 3)\n \n Returns:\n probability: float [0, 1]\n uncertainty: float (std of predictions)\n \"\"\"\n all_predictions = []\n \n # Ensemble over folds\n for model_idx, model in enumerate(self.models):\n fold_predictions = []\n \n # TTA\n if self.use_tta:\n for tta_idx in range(self.tta_iterations):\n # Apply TTA transform\n augmented = self.tta_transform.get_transform(tta_idx)(volume)\n \n # Convert to tensor\n volume_tensor = torch.from_numpy(augmented).permute(3,0,1,2).unsqueeze(0).float().to(device)\n \n # Predict\n with torch.no_grad():\n with autocast(enabled=CONFIG['use_mixed_precision']):\n logit = model(volume_tensor)\n prob = torch.sigmoid(logit).cpu().item()\n \n fold_predictions.append(prob)\n \n # Average TTA predictions for this fold\n fold_pred = np.mean(fold_predictions)\n else:\n # No TTA - single prediction\n volume_tensor = torch.from_numpy(volume).permute(3,0,1,2).unsqueeze(0).float().to(device)\n \n with torch.no_grad():\n with autocast(enabled=CONFIG['use_mixed_precision']):\n logit = model(volume_tensor)\n fold_pred = torch.sigmoid(logit).cpu().item()\n \n all_predictions.append(fold_pred)\n \n # Weighted ensemble\n if CONFIG['use_weighted_ensemble'] and len(self.fold_weights) > 0:\n mean_prob = sum(p * w for p, w in zip(all_predictions, self.fold_weights))\n else:\n mean_prob = np.mean(all_predictions)\n \n uncertainty = np.std(all_predictions)\n \n # Temperature scaling calibration (now disabled by default)\n if CONFIG['calibrate_predictions']:\n logit = np.log(mean_prob / (1 - mean_prob + 1e-8))\n calibrated_logit = logit / CONFIG['temperature']\n mean_prob = 1 / (1 + np.exp(-calibrated_logit))\n \n return mean_prob, uncertainty\n\n\n# ============================================================================\n# ENHANCED PREPROCESSING (DICOM FIX)\n# ============================================================================\n\ndef load_and_preprocess_test_patient(patient_folder):\n \"\"\"\n Load and preprocess test patient with BETTER error handling\n \"\"\"\n import pydicom\n from glob import glob\n from scipy.ndimage import zoom\n \n dicom_files = sorted(glob(os.path.join(patient_folder, \"*.dcm\")))\n \n if len(dicom_files) == 0:\n raise ValueError(\"No DICOM files found\")\n \n # Load slices with better error handling\n slices = []\n for dcm_file in dicom_files:\n try:\n # FIX: Force decompression for JPEG Lossless\n ds = pydicom.dcmread(dcm_file, force=True)\n \n # Try to decompress if compressed\n if hasattr(ds, 'decompress'):\n try:\n ds.decompress()\n except:\n pass\n \n if hasattr(ds, 'ImagePositionPatient') and hasattr(ds, 'PixelSpacing'):\n slices.append(ds)\n except Exception as e:\n # Skip problematic slices\n continue\n \n if len(slices) == 0:\n raise ValueError(\"No valid slices after decompression attempts\")\n \n slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n \n # Stack volume\n volume = np.stack([s.pixel_array for s in slices])\n metadata = slices[0]\n \n # Convert to HU\n try:\n intercept = float(metadata.RescaleIntercept)\n slope = float(metadata.RescaleSlope)\n except:\n intercept = 0.0\n slope = 1.0\n \n volume_hu = volume.astype(np.float32) * slope + intercept\n \n # Create multi-window (bone, soft tissue, wide)\n def apply_windowing(vol_hu, center, width):\n lower = center - width // 2\n upper = center + width // 2\n windowed = np.clip(vol_hu, lower, upper)\n normalized = (windowed - lower) / (upper - lower)\n return normalized.astype(np.float32)\n \n bone = apply_windowing(volume_hu, 400, 1800)\n soft = apply_windowing(volume_hu, 40, 400)\n wide = apply_windowing(volume_hu, 400, 4000)\n \n # Stack channels\n multi_channel = np.stack([bone, soft, wide], axis=-1)\n \n # Resample to target shape (64, 224, 224, 3)\n target_shape = (64, 224, 224)\n current_shape = multi_channel.shape[:3]\n \n resize_factor = np.array(target_shape) / np.array(current_shape)\n resize_factor = np.append(resize_factor, 1.0)\n \n resampled = zoom(multi_channel, resize_factor, order=1).astype(np.float32)\n \n # Crop/pad if needed\n final = np.zeros((*target_shape, 3), dtype=np.float32)\n \n min_d = min(resampled.shape[0], target_shape[0])\n min_h = min(resampled.shape[1], target_shape[1])\n min_w = min(resampled.shape[2], target_shape[2])\n \n final[:min_d, :min_h, :min_w, :] = resampled[:min_d, :min_h, :min_w, :]\n \n return final\n\n\n# ============================================================================\n# SUBMISSION GENERATION\n# ============================================================================\n\ndef generate_submission(predictor, test_patients, save_path):\n \"\"\"Generate submission file with better error tracking\"\"\"\n \n print(f\"\\n{'='*80}\")\n print(\"📋 GENERATING PREDICTIONS\")\n print(f\"{'='*80}\")\n \n results = []\n failed_count = 0\n success_count = 0\n \n for patient_id in tqdm(test_patients, desc=\"Processing test patients\"):\n patient_folder = os.path.join(CONFIG['test_dir'], patient_id)\n \n try:\n # Preprocess\n volume = load_and_preprocess_test_patient(patient_folder)\n \n # Predict\n probability, uncertainty = predictor.predict_single(volume)\n \n success_count += 1\n \n # For competition: need per-vertebra predictions\n for vertebra in ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7']:\n results.append({\n 'StudyInstanceUID': patient_id,\n 'prediction_id': f'{patient_id}_{vertebra}',\n 'fractured': probability\n })\n \n results.append({\n 'StudyInstanceUID': patient_id,\n 'prediction_id': f'{patient_id}_patient_overall',\n 'fractured': probability\n })\n \n except Exception as e:\n failed_count += 1\n print(f\"\\n ⚠️ Failed {patient_id}: {str(e)[:80]}\")\n \n # Add default predictions\n for vertebra in ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7']:\n results.append({\n 'StudyInstanceUID': patient_id,\n 'prediction_id': f'{patient_id}_{vertebra}',\n 'fractured': 0.5\n })\n results.append({\n 'StudyInstanceUID': patient_id,\n 'prediction_id': f'{patient_id}_patient_overall',\n 'fractured': 0.5\n })\n \n # Create submission dataframe\n submission_df = pd.DataFrame(results)\n submission_df.to_csv(save_path, index=False)\n \n print(f\"\\n ✓ Saved submission: {save_path}\")\n print(f\" ✓ Total predictions: {len(submission_df)}\")\n print(f\" ✓ Success: {success_count}/{len(test_patients)} ({100*success_count/len(test_patients):.1f}%)\")\n print(f\" ⚠️ Failed: {failed_count}/{len(test_patients)} ({100*failed_count/len(test_patients):.1f}%)\")\n \n return submission_df\n\n\n# ============================================================================\n# ENHANCED VISUALIZATION\n# ============================================================================\n\ndef visualize_predictions(submission_df, save_path):\n \"\"\"Enhanced visualization with more insights\"\"\"\n \n fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n \n overall_mask = submission_df['prediction_id'].str.contains('patient_overall')\n overall_preds = submission_df[overall_mask]['fractured']\n \n # Overall distribution\n axes[0, 0].hist(overall_preds, bins=50, alpha=0.7, color='blue', edgecolor='black')\n axes[0, 0].axvline(overall_preds.mean(), color='red', linestyle='--', linewidth=2, label=f'Mean: {overall_preds.mean():.3f}')\n axes[0, 0].set_title('Patient Overall Prediction Distribution', fontweight='bold')\n axes[0, 0].set_xlabel('Fracture Probability')\n axes[0, 0].set_ylabel('Count')\n axes[0, 0].legend()\n axes[0, 0].grid(True, alpha=0.3)\n \n # Per-vertebra distribution\n vertebra_mask = ~overall_mask\n axes[0, 1].hist(submission_df[vertebra_mask]['fractured'], bins=50, alpha=0.7, \n color='green', edgecolor='black')\n axes[0, 1].set_title('Per-Vertebra Prediction Distribution', fontweight='bold')\n axes[0, 1].set_xlabel('Fracture Probability')\n axes[0, 1].set_ylabel('Count')\n axes[0, 1].grid(True, alpha=0.3)\n \n # Vertebra comparison\n vertebra_means = []\n vertebra_labels = []\n for vert in ['C1', 'C2', 'C3', 'C4', 'C5', 'C6', 'C7']:\n vert_data = submission_df[submission_df['prediction_id'].str.contains(f'_{vert}')]\n if len(vert_data) > 0:\n vertebra_means.append(vert_data['fractured'].mean())\n vertebra_labels.append(vert)\n \n axes[1, 0].bar(vertebra_labels, vertebra_means, color='orange', edgecolor='black')\n axes[1, 0].set_title('Mean Fracture Probability by Vertebra', fontweight='bold')\n axes[1, 0].set_xlabel('Vertebra')\n axes[1, 0].set_ylabel('Mean Probability')\n axes[1, 0].grid(True, alpha=0.3, axis='y')\n \n # Summary statistics\n axes[1, 1].axis('off')\n \n summary_text = f\"\"\"\n PREDICTION SUMMARY (v4.0)\n {'='*40}\n \n Total Patients: {len(overall_preds)}\n \n Overall Fracture Predictions:\n • Mean: {overall_preds.mean():.3f}\n • Median: {overall_preds.median():.3f}\n • Std: {overall_preds.std():.3f}\n • Min: {overall_preds.min():.3f}\n • Max: {overall_preds.max():.3f}\n \n Risk Distribution:\n • High risk (>0.7): {(overall_preds > 0.7).sum()}\n • Medium (0.3-0.7): {((overall_preds >= 0.3) & (overall_preds <= 0.7)).sum()}\n • Low risk (<0.3): {(overall_preds < 0.3).sum()}\n \n Improvements Applied:\n ✓ Calibration disabled\n ✓ Weighted ensemble: {CONFIG['use_weighted_ensemble']}\n ✓ DICOM handling enhanced\n ✓ TTA optimized\n \"\"\"\n \n axes[1, 1].text(0.05, 0.5, summary_text, fontsize=9, family='monospace',\n verticalalignment='center',\n bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.3))\n \n plt.suptitle('Enhanced Ensemble Prediction Analysis v4.0', fontsize=14, fontweight='bold')\n plt.tight_layout()\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved visualization: {save_path}\")\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n print(f\"\\n{'='*80}\")\n print(\"🚀 STARTING ENHANCED ENSEMBLE INFERENCE\")\n print(f\"{'='*80}\")\n \n # Initialize predictor\n predictor = EnsemblePredictor(\n models_dir=CONFIG['models_dir'],\n n_folds=CONFIG['n_folds'],\n use_tta=CONFIG['use_tta'],\n tta_iterations=CONFIG['tta_iterations']\n )\n \n if len(predictor.models) == 0:\n print(f\"\\n❌ ERROR: No models found in {CONFIG['models_dir']}\")\n print(f\" Run Step 3 (training) first!\")\n return\n \n # Get test patients\n if os.path.exists(CONFIG['test_dir']):\n test_patients = [d for d in os.listdir(CONFIG['test_dir']) \n if os.path.isdir(os.path.join(CONFIG['test_dir'], d))]\n print(f\"\\n ✓ Found {len(test_patients)} test patients\")\n else:\n print(f\"\\n ⚠️ Test directory not found: {CONFIG['test_dir']}\")\n return\n \n # Generate submission\n submission_path = os.path.join(CONFIG['output_dir'], 'submission_v4.csv')\n submission_df = generate_submission(predictor, test_patients, submission_path)\n \n # Visualize\n viz_path = os.path.join(CONFIG['output_dir'], 'predictions_analysis_v4.png')\n visualize_predictions(submission_df, viz_path)\n \n # Final summary\n print(f\"\\n{'='*80}\")\n print(\"✅ ENHANCED INFERENCE COMPLETE!\")\n print(f\"{'='*80}\")\n \n print(f\"\\n📊 Results:\")\n print(f\" • Mean probability: {submission_df['fractured'].mean():.3f}\")\n print(f\" • Std: {submission_df['fractured'].std():.3f}\")\n print(f\" • Range: [{submission_df['fractured'].min():.3f}, {submission_df['fractured'].max():.3f}]\")\n \n print(f\"\\n✨ Enhancements Applied:\")\n print(f\" ✓ Calibration disabled (was causing centered predictions)\")\n print(f\" ✓ Weighted ensemble by validation AUC\")\n print(f\" ✓ Improved DICOM decompression\")\n print(f\" ✓ Optimized TTA (4 iterations)\")\n \n print(f\"\\n💡 Next: Submit submission_v4.csv to Kaggle!\")\n\n\nif __name__ == \"__main__\":\n try:\n main()\n except KeyboardInterrupt:\n print(f\"\\n\\n⚠️ Interrupted by user\")\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T13:29:42.23865Z","iopub.execute_input":"2026-01-22T13:29:42.239012Z","iopub.status.idle":"2026-01-22T13:29:58.52173Z","shell.execute_reply.started":"2026-01-22T13:29:42.238984Z","shell.execute_reply":"2026-01-22T13:29:58.521099Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nPIPELINE DIAGNOSTIC & OPTIMIZATION TOOL v1.0\n=============================================\nRun this AFTER completing Steps 1-5 to:\n✓ Verify all components are working\n✓ Identify bottlenecks and issues\n✓ Get optimization recommendations\n✓ Generate comprehensive report\n\nThis will tell you exactly what to do next!\n\"\"\"\n\nimport os\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom matplotlib.gridspec import GridSpec\nfrom glob import glob\n\nprint(\"=\"*80)\nprint(\"🔍 COMPREHENSIVE PIPELINE DIAGNOSTIC v1.0\")\nprint(\"=\"*80)\n\n# ============================================================================\n# CONFIGURATION\n# ============================================================================\n\nCONFIG = {\n 'dataset_dir': '/kaggle/working/enhanced_dataset_v3',\n 'training_dir': '/kaggle/working/training_output',\n 'inference_dir': '/kaggle/working/inference_output',\n 'output_dir': '/kaggle/working/diagnostic_report',\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\n\n# ============================================================================\n# DIAGNOSTIC CHECKS\n# ============================================================================\n\nclass PipelineDiagnostic:\n \"\"\"Comprehensive pipeline health check\"\"\"\n \n def __init__(self):\n self.issues = []\n self.warnings = []\n self.successes = []\n self.recommendations = []\n \n def check_preprocessing(self):\n \"\"\"Verify Step 1: Preprocessing\"\"\"\n print(f\"\\n{'='*80}\")\n print(\"📊 STEP 1: PREPROCESSING CHECK\")\n print(f\"{'='*80}\")\n \n metadata_path = os.path.join(CONFIG['dataset_dir'], 'metadata.csv')\n volumes_dir = os.path.join(CONFIG['dataset_dir'], 'volumes')\n \n # Check metadata exists\n if not os.path.exists(metadata_path):\n self.issues.append(\"❌ metadata.csv not found - Run Step 1 first!\")\n print(\" ❌ CRITICAL: metadata.csv not found\")\n return False\n \n # Load and analyze metadata\n metadata_df = pd.read_csv(metadata_path)\n print(f\"\\n ✓ Found metadata: {len(metadata_df)} patients\")\n \n # Check class balance\n fracture_ratio = metadata_df['has_fracture'].mean()\n print(f\" • Fracture ratio: {fracture_ratio:.2%}\")\n \n if fracture_ratio < 0.3 or fracture_ratio > 0.7:\n self.warnings.append(f\"⚠️ Class imbalance: {fracture_ratio:.2%} fractures\")\n self.recommendations.append(\"Use focal loss or class weights in training\")\n else:\n self.successes.append(\"✓ Good class balance\")\n \n # Check volume files\n volume_files = glob(os.path.join(volumes_dir, \"*.npy\"))\n print(f\" • Volume files: {len(volume_files)}\")\n \n if len(volume_files) != len(metadata_df):\n self.issues.append(f\"❌ Mismatch: {len(metadata_df)} metadata vs {len(volume_files)} volumes\")\n else:\n self.successes.append(\"✓ All volumes present\")\n \n # Verify multi-channel\n if len(volume_files) > 0:\n sample_vol = np.load(volume_files[0])\n print(f\" • Sample shape: {sample_vol.shape}\")\n \n if sample_vol.shape == (64, 224, 224, 3):\n self.successes.append(\"✓ Correct 3-channel multi-window preprocessing\")\n else:\n self.issues.append(f\"❌ Wrong shape: {sample_vol.shape}, expected (64,224,224,3)\")\n \n # Check value range\n print(f\" • Value range: [{sample_vol.min():.3f}, {sample_vol.max():.3f}]\")\n if sample_vol.min() < 0 or sample_vol.max() > 1:\n self.warnings.append(\"⚠️ Values outside [0,1] range\")\n \n # Check quality metrics\n if 'bone_mean' in metadata_df.columns:\n self.successes.append(\"✓ Multi-channel statistics computed\")\n else:\n self.warnings.append(\"⚠️ Missing multi-channel statistics\")\n \n return True\n \n def check_training(self):\n \"\"\"Verify Step 3: Training\"\"\"\n print(f\"\\n{'='*80}\")\n print(\"🧠 STEP 3: TRAINING CHECK\")\n print(f\"{'='*80}\")\n \n results_path = os.path.join(CONFIG['training_dir'], 'results.json')\n \n if not os.path.exists(results_path):\n self.issues.append(\"❌ results.json not found - Training incomplete or not run\")\n print(\" ❌ CRITICAL: Training results not found\")\n return False\n \n # Load results\n with open(results_path, 'r') as f:\n results = json.load(f)\n \n mean_auc = results['mean_auc']\n std_auc = results['std_auc']\n fold_aucs = results['fold_aucs']\n \n print(f\"\\n 📊 Cross-Validation Results:\")\n print(f\" Mean AUC: {mean_auc:.4f} ± {std_auc:.4f}\")\n \n for i, auc in enumerate(fold_aucs):\n status = \"✓\" if auc > 0.65 else \"⚠️\"\n print(f\" Fold {i+1}: {auc:.4f} {status}\")\n \n # Evaluate performance\n if mean_auc > 0.80:\n self.successes.append(f\"✓ EXCELLENT: Mean AUC = {mean_auc:.4f}\")\n print(f\"\\n 🎉 EXCELLENT PERFORMANCE!\")\n elif mean_auc > 0.75:\n self.successes.append(f\"✓ GOOD: Mean AUC = {mean_auc:.4f}\")\n print(f\"\\n ✓ Good performance, ready for submission\")\n elif mean_auc > 0.65:\n self.warnings.append(f\"⚠️ MODERATE: Mean AUC = {mean_auc:.4f}\")\n print(f\"\\n ⚠️ Moderate performance - optimization needed\")\n self.recommendations.extend([\n \"Increase epochs to 100\",\n \"Try larger base_filters (32 or 64)\",\n \"Use stronger augmentation\",\n \"Check for overfitting in training curves\"\n ])\n else:\n self.issues.append(f\"❌ POOR: Mean AUC = {mean_auc:.4f}\")\n print(f\"\\n ❌ Poor performance - serious issues\")\n self.recommendations.extend([\n \"CRITICAL: Verify preprocessing is correct\",\n \"Check data quality - review failed_patients.json\",\n \"Disable augmentation temporarily to debug\",\n \"Verify labels match actual fractures\",\n \"Try simpler model first (reduce complexity)\"\n ])\n \n # Check variance\n if std_auc > 0.1:\n self.warnings.append(f\"⚠️ High variance: ±{std_auc:.4f}\")\n self.recommendations.append(\"Reduce variance: use more data or regularization\")\n else:\n self.successes.append(f\"✓ Low variance: ±{std_auc:.4f}\")\n \n # Check model files\n model_files = glob(os.path.join(CONFIG['training_dir'], \"fold*_best.pth\"))\n print(f\"\\n • Model checkpoints: {len(model_files)}/5\")\n \n if len(model_files) < 5:\n self.warnings.append(f\"⚠️ Missing model checkpoints: {5-len(model_files)}\")\n \n return True\n \n def check_inference(self):\n \"\"\"Verify Step 5: Inference\"\"\"\n print(f\"\\n{'='*80}\")\n print(\"🎯 STEP 5: INFERENCE CHECK\")\n print(f\"{'='*80}\")\n \n submission_path = os.path.join(CONFIG['inference_dir'], 'submission_v4.csv')\n \n if not os.path.exists(submission_path):\n self.warnings.append(\"⚠️ submission_v4.csv not found - Run Step 5\")\n print(\" ⚠️ Submission file not found\")\n return False\n \n # Load submission\n submission_df = pd.read_csv(submission_path)\n print(f\"\\n ✓ Found submission: {len(submission_df)} predictions\")\n \n # Analyze predictions\n overall_mask = submission_df['prediction_id'].str.contains('patient_overall')\n overall_preds = submission_df[overall_mask]['fractured']\n \n mean_pred = overall_preds.mean()\n std_pred = overall_preds.std()\n \n print(f\"\\n 📊 Prediction Statistics:\")\n print(f\" Mean: {mean_pred:.3f}\")\n print(f\" Std: {std_pred:.3f}\")\n print(f\" Range: [{overall_preds.min():.3f}, {overall_preds.max():.3f}]\")\n \n # Check for common issues\n if abs(mean_pred - 0.5) < 0.05:\n self.issues.append(\"❌ CRITICAL: Predictions centered at 0.5 (model not confident)\")\n self.recommendations.extend([\n \"URGENT: Model is not learning properly\",\n \"Check if calibration is disabled (should be False)\",\n \"Verify weighted ensemble is enabled\",\n \"Review training - AUC might be too low\"\n ])\n else:\n self.successes.append(f\"✓ Predictions diverse (mean={mean_pred:.3f})\")\n \n if std_pred < 0.1:\n self.warnings.append(f\"⚠️ Low diversity: std={std_pred:.3f}\")\n self.recommendations.append(\"Predictions too similar - check ensemble diversity\")\n elif std_pred > 0.3:\n self.warnings.append(f\"⚠️ High diversity: std={std_pred:.3f}\")\n self.recommendations.append(\"Predictions very spread - check for outliers\")\n else:\n self.successes.append(f\"✓ Good diversity: std={std_pred:.3f}\")\n \n # Distribution check\n high_risk = (overall_preds > 0.7).sum()\n low_risk = (overall_preds < 0.3).sum()\n medium = len(overall_preds) - high_risk - low_risk\n \n print(f\"\\n 📈 Risk Distribution:\")\n print(f\" High (>0.7): {high_risk} ({100*high_risk/len(overall_preds):.1f}%)\")\n print(f\" Medium: {medium} ({100*medium/len(overall_preds):.1f}%)\")\n print(f\" Low (<0.3): {low_risk} ({100*low_risk/len(overall_preds):.1f}%)\")\n \n return True\n \n def generate_report(self):\n \"\"\"Generate comprehensive diagnostic report\"\"\"\n print(f\"\\n{'='*80}\")\n print(\"📋 GENERATING COMPREHENSIVE REPORT\")\n print(f\"{'='*80}\")\n \n # Create detailed report\n fig = plt.figure(figsize=(16, 20))\n gs = GridSpec(5, 2, figure=fig, hspace=0.4, wspace=0.3)\n \n # Title\n title_ax = fig.add_subplot(gs[0, :])\n title_ax.axis('off')\n title_text = \"\"\"\n PIPELINE DIAGNOSTIC REPORT\n ═══════════════════════════════════════════════════════════\n Generated to identify issues and provide optimization guidance\n \"\"\"\n title_ax.text(0.5, 0.5, title_text, fontsize=14, fontweight='bold',\n ha='center', va='center', family='monospace',\n bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))\n \n # Successes\n success_ax = fig.add_subplot(gs[1, 0])\n success_ax.axis('off')\n success_text = \"✅ SUCCESSES\\n\" + \"=\"*40 + \"\\n\\n\"\n if len(self.successes) > 0:\n success_text += \"\\n\".join(self.successes)\n else:\n success_text += \"No major successes detected\"\n \n success_ax.text(0.05, 0.95, success_text, fontsize=10, family='monospace',\n va='top', bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.3))\n \n # Issues\n issues_ax = fig.add_subplot(gs[1, 1])\n issues_ax.axis('off')\n issues_text = \"❌ CRITICAL ISSUES\\n\" + \"=\"*40 + \"\\n\\n\"\n if len(self.issues) > 0:\n issues_text += \"\\n\".join(self.issues)\n else:\n issues_text += \"No critical issues found!\"\n \n issues_ax.text(0.05, 0.95, issues_text, fontsize=10, family='monospace',\n va='top', bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.3))\n \n # Warnings\n warnings_ax = fig.add_subplot(gs[2, 0])\n warnings_ax.axis('off')\n warnings_text = \"⚠️ WARNINGS\\n\" + \"=\"*40 + \"\\n\\n\"\n if len(self.warnings) > 0:\n warnings_text += \"\\n\".join(self.warnings)\n else:\n warnings_text += \"No warnings\"\n \n warnings_ax.text(0.05, 0.95, warnings_text, fontsize=10, family='monospace',\n va='top', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.3))\n \n # Recommendations\n rec_ax = fig.add_subplot(gs[2, 1])\n rec_ax.axis('off')\n rec_text = \"💡 RECOMMENDATIONS\\n\" + \"=\"*40 + \"\\n\\n\"\n if len(self.recommendations) > 0:\n rec_text += \"\\n\".join([f\"{i+1}. {r}\" for i, r in enumerate(self.recommendations)])\n else:\n rec_text += \"Pipeline looks good!\\nReady for submission.\"\n \n rec_ax.text(0.05, 0.95, rec_text, fontsize=9, family='monospace',\n va='top', bbox=dict(boxstyle='round', facecolor='lightcyan', alpha=0.3))\n \n # Overall status\n status_ax = fig.add_subplot(gs[3, :])\n status_ax.axis('off')\n \n if len(self.issues) == 0 and len(self.warnings) <= 2:\n status_color = 'lightgreen'\n status_emoji = \"🎉\"\n status_msg = \"EXCELLENT - Pipeline Ready for Production!\"\n next_step = \"→ Submit to Kaggle competition\"\n elif len(self.issues) == 0:\n status_color = 'lightyellow'\n status_emoji = \"✓\"\n status_msg = \"GOOD - Minor optimizations recommended\"\n next_step = \"→ Address warnings, then submit\"\n else:\n status_color = 'lightcoral'\n status_emoji = \"⚠️\"\n status_msg = \"NEEDS ATTENTION - Critical issues found\"\n next_step = \"→ Fix critical issues before submission\"\n \n status_text = f\"\"\"\n {status_emoji} OVERALL STATUS: {status_msg}\n {'='*60}\n \n Issues: {len(self.issues)}\n Warnings: {len(self.warnings)}\n Successes: {len(self.successes)}\n \n NEXT STEP: {next_step}\n \"\"\"\n \n status_ax.text(0.5, 0.5, status_text, fontsize=12, fontweight='bold',\n ha='center', va='center', family='monospace',\n bbox=dict(boxstyle='round', facecolor=status_color, alpha=0.5))\n \n # Optimization suggestions\n opt_ax = fig.add_subplot(gs[4, :])\n opt_ax.axis('off')\n \n opt_text = \"\"\"\n 🔧 OPTIMIZATION PRIORITY MATRIX\n ═══════════════════════════════════════════════════════════\n \n IF AUC > 0.80: ✓ Submit now, then try:\n • Ensemble with different architectures\n • Hyperparameter tuning\n • External data augmentation\n \n IF AUC 0.75-0.80: ✓ Submit, but consider:\n • More epochs (100 instead of 50)\n • Larger model (base_filters=32)\n • Learning rate tuning\n \n IF AUC 0.65-0.75: ⚠️ Optimize before submitting:\n • Check training curves for overfitting\n • Increase augmentation strength\n • Try focal loss with different α/γ\n • Verify data quality\n \n IF AUC < 0.65: ❌ Critical debugging needed:\n • Verify preprocessing visually\n • Check label alignment\n • Try simpler baseline first\n • Review failed_patients.json\n • Disable augmentation to isolate issue\n \n COMMON FIXES:\n • Predictions at 0.5 → Disable calibration, check training\n • Low variance → More data or stronger augmentation\n • High variance → More regularization (dropout, weight decay)\n • OOM errors → Reduce batch_size to 2 or 1\n \"\"\"\n \n opt_ax.text(0.05, 0.95, opt_text, fontsize=9, family='monospace',\n va='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))\n \n plt.suptitle('Cervical Spine Fracture Detection - Pipeline Diagnostic Report',\n fontsize=16, fontweight='bold')\n \n report_path = os.path.join(CONFIG['output_dir'], 'diagnostic_report.png')\n plt.savefig(report_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved: {report_path}\")\n \n # Save text report\n text_report_path = os.path.join(CONFIG['output_dir'], 'diagnostic_report.txt')\n with open(text_report_path, 'w') as f:\n f.write(\"PIPELINE DIAGNOSTIC REPORT\\n\")\n f.write(\"=\"*80 + \"\\n\\n\")\n \n f.write(\"SUCCESSES:\\n\")\n f.write(\"-\"*80 + \"\\n\")\n for s in self.successes:\n f.write(f\"{s}\\n\")\n \n f.write(\"\\n\\nISSUES:\\n\")\n f.write(\"-\"*80 + \"\\n\")\n for i in self.issues:\n f.write(f\"{i}\\n\")\n \n f.write(\"\\n\\nWARNINGS:\\n\")\n f.write(\"-\"*80 + \"\\n\")\n for w in self.warnings:\n f.write(f\"{w}\\n\")\n \n f.write(\"\\n\\nRECOMMENDATIONS:\\n\")\n f.write(\"-\"*80 + \"\\n\")\n for idx, r in enumerate(self.recommendations, 1):\n f.write(f\"{idx}. {r}\\n\")\n \n print(f\" ✓ Saved: {text_report_path}\")\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n \"\"\"Run comprehensive diagnostic\"\"\"\n \n diagnostic = PipelineDiagnostic()\n \n # Run all checks\n preprocessing_ok = diagnostic.check_preprocessing()\n if preprocessing_ok:\n training_ok = diagnostic.check_training()\n if training_ok:\n diagnostic.check_inference()\n \n # Generate report\n diagnostic.generate_report()\n \n # Print summary\n print(f\"\\n{'='*80}\")\n print(\"📊 DIAGNOSTIC SUMMARY\")\n print(f\"{'='*80}\")\n \n print(f\"\\n✅ Successes: {len(diagnostic.successes)}\")\n for s in diagnostic.successes[:5]: # Show first 5\n print(f\" {s}\")\n if len(diagnostic.successes) > 5:\n print(f\" ... and {len(diagnostic.successes)-5} more\")\n \n print(f\"\\n❌ Critical Issues: {len(diagnostic.issues)}\")\n for i in diagnostic.issues:\n print(f\" {i}\")\n \n print(f\"\\n⚠️ Warnings: {len(diagnostic.warnings)}\")\n for w in diagnostic.warnings[:5]:\n print(f\" {w}\")\n if len(diagnostic.warnings) > 5:\n print(f\" ... and {len(diagnostic.warnings)-5} more\")\n \n print(f\"\\n💡 Top Recommendations:\")\n for idx, r in enumerate(diagnostic.recommendations[:5], 1):\n print(f\" {idx}. {r}\")\n \n print(f\"\\n{'='*80}\")\n print(\"✅ DIAGNOSTIC COMPLETE!\")\n print(f\"{'='*80}\")\n print(f\"\\nReview the full report at:\")\n print(f\" {CONFIG['output_dir']}/diagnostic_report.png\")\n print(f\" {CONFIG['output_dir']}/diagnostic_report.txt\")\n\n\nif __name__ == \"__main__\":\n try:\n main()\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ DIAGNOSTIC ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T13:35:31.898771Z","iopub.execute_input":"2026-01-22T13:35:31.899472Z","iopub.status.idle":"2026-01-22T13:35:32.643188Z","shell.execute_reply.started":"2026-01-22T13:35:31.899438Z","shell.execute_reply":"2026-01-22T13:35:32.642555Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install -q gdcm pylibjpeg pylibjpeg-libjpeg\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T13:56:28.465582Z","iopub.execute_input":"2026-01-22T13:56:28.466458Z","iopub.status.idle":"2026-01-22T13:56:33.412903Z","shell.execute_reply.started":"2026-01-22T13:56:28.466418Z","shell.execute_reply":"2026-01-22T13:56:33.41207Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport zipfile\n\nfolders_to_zip = [\n \"augmentation_analysis\",\n \"diagnostic_report\",\n \"inference_output\",\n \"step2_exploration\",\n \"training_output\",\n]\n\nbase_path = \"/kaggle/working\"\nzip_path = os.path.join(base_path, \"reference_output.zip\")\n\nwith zipfile.ZipFile(zip_path, \"w\", zipfile.ZIP_DEFLATED) as zipf:\n for folder in folders_to_zip:\n folder_path = os.path.join(base_path, folder)\n if not os.path.exists(folder_path):\n print(f\"⚠️ Skipping (not found): {folder}\")\n continue\n \n for root, dirs, files in os.walk(folder_path):\n for file in files:\n file_path = os.path.join(root, file)\n arcname = os.path.relpath(file_path, base_path)\n zipf.write(file_path, arcname)\n\nprint(\"✅ ZIP created successfully:\")\nprint(zip_path)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T15:28:01.194977Z","iopub.execute_input":"2026-01-22T15:28:01.195625Z","iopub.status.idle":"2026-01-22T15:28:07.619788Z","shell.execute_reply.started":"2026-01-22T15:28:01.195599Z","shell.execute_reply":"2026-01-22T15:28:07.619064Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!zip -r FinalOutput.zip /kaggle/working\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-01-22T14:17:45.695325Z","iopub.execute_input":"2026-01-22T14:17:45.69619Z","iopub.status.idle":"2026-01-22T14:20:05.601603Z","shell.execute_reply.started":"2026-01-22T14:17:45.696155Z","shell.execute_reply":"2026-01-22T14:20:05.600864Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"\nSTEP 5: GRAD-CAM EXPLAINABLE AI (XAI) v1.0\n==========================================\nAdd interpretability to your cervical spine fracture detection model\n\nFEATURES:\n✅ 3D GRAD-CAM for volumetric medical imaging\n✅ Multi-slice visualization\n✅ Overlay heatmaps on original CT scans\n✅ Per-vertebra attention maps\n✅ Batch processing for test set\n✅ Interactive HTML reports\n✅ Quantitative saliency metrics\n\nUse this to:\n- Understand what the model looks at\n- Validate clinical relevance\n- Debug model behavior\n- Build trust with radiologists\n\"\"\"\n\nimport os\nimport json\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nfrom matplotlib import cm\nimport cv2\n\nimport torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nfrom torch.cuda.amp import autocast\n\nprint(\"=\"*80)\nprint(\"🔍 STEP 5: GRAD-CAM EXPLAINABLE AI v1.0\")\nprint(\"=\"*80)\n\n# ============================================================================\n# CONFIGURATION\n# ============================================================================\n\nCONFIG = {\n # Paths\n 'models_dir': '/kaggle/working/training_output',\n 'test_dir': '/kaggle/input/rsna-2022-cervical-spine-fracture-detection/test_images',\n 'output_dir': '/kaggle/working/gradcam_output',\n \n # Model settings\n 'fold_to_use': 0, # Which fold's model to use (0-4)\n 'in_channels': 3,\n 'base_filters': 16,\n 'dropout': 0.3,\n \n # GRAD-CAM settings\n 'target_layer': 'layer3', # Which layer to visualize\n 'num_slices_to_show': 16, # How many slices to visualize per patient\n 'alpha': 0.4, # Overlay transparency\n 'colormap': 'jet', # 'jet', 'hot', 'viridis'\n \n # Processing\n 'batch_size': 1, # Keep 1 for GRAD-CAM\n 'num_patients_to_visualize': 20, # Limit for demo\n 'generate_html_report': True,\n \n # System\n 'use_mixed_precision': False, # Disable for GRAD-CAM\n 'device': 'cuda' if torch.cuda.is_available() else 'cpu',\n}\n\nos.makedirs(CONFIG['output_dir'], exist_ok=True)\nos.makedirs(os.path.join(CONFIG['output_dir'], 'heatmaps'), exist_ok=True)\n\ndevice = torch.device(CONFIG['device'])\nprint(f\"\\n📋 Configuration:\")\nprint(f\" • Device: {device}\")\nprint(f\" • Target layer: {CONFIG['target_layer']}\")\nprint(f\" • Colormap: {CONFIG['colormap']}\")\nprint(f\" • Patients to visualize: {CONFIG['num_patients_to_visualize']}\")\n\n# ============================================================================\n# MODEL ARCHITECTURE (same as before)\n# ============================================================================\n\nclass ResidualBlock3D(nn.Module):\n def __init__(self, in_channels, out_channels, stride=1):\n super().__init__()\n self.conv1 = nn.Conv3d(in_channels, out_channels, kernel_size=3, \n stride=stride, padding=1, bias=False)\n self.bn1 = nn.BatchNorm3d(out_channels)\n self.relu = nn.ReLU(inplace=True)\n self.conv2 = nn.Conv3d(out_channels, out_channels, kernel_size=3,\n stride=1, padding=1, bias=False)\n self.bn2 = nn.BatchNorm3d(out_channels)\n \n if stride != 1 or in_channels != out_channels:\n self.shortcut = nn.Sequential(\n nn.Conv3d(in_channels, out_channels, kernel_size=1, \n stride=stride, bias=False),\n nn.BatchNorm3d(out_channels)\n )\n else:\n self.shortcut = nn.Identity()\n \n def forward(self, x):\n residual = x\n out = self.conv1(x)\n out = self.bn1(out)\n out = self.relu(out)\n out = self.conv2(out)\n out = self.bn2(out)\n out += self.shortcut(residual)\n out = self.relu(out)\n return out\n\n\nclass ResNet3D(nn.Module):\n def __init__(self, in_channels=3, base_filters=16, dropout=0.3):\n super().__init__()\n self.conv1 = nn.Conv3d(in_channels, base_filters, kernel_size=7, \n stride=2, padding=3, bias=False)\n self.bn1 = nn.BatchNorm3d(base_filters)\n self.relu = nn.ReLU(inplace=True)\n self.maxpool = nn.MaxPool3d(kernel_size=3, stride=2, padding=1)\n \n self.layer1 = self._make_layer(base_filters, base_filters*2, blocks=2, stride=1)\n self.layer2 = self._make_layer(base_filters*2, base_filters*4, blocks=2, stride=2)\n self.layer3 = self._make_layer(base_filters*4, base_filters*8, blocks=2, stride=2)\n \n self.avgpool = nn.AdaptiveAvgPool3d(1)\n self.dropout = nn.Dropout(dropout)\n self.fc = nn.Linear(base_filters*8, 1)\n \n def _make_layer(self, in_channels, out_channels, blocks, stride):\n layers = []\n layers.append(ResidualBlock3D(in_channels, out_channels, stride))\n for _ in range(1, blocks):\n layers.append(ResidualBlock3D(out_channels, out_channels, stride=1))\n return nn.Sequential(*layers)\n \n def forward(self, x):\n x = self.conv1(x)\n x = self.bn1(x)\n x = self.relu(x)\n x = self.maxpool(x)\n x = self.layer1(x)\n x = self.layer2(x)\n x = self.layer3(x)\n x = self.avgpool(x)\n x = x.view(x.size(0), -1)\n x = self.dropout(x)\n x = self.fc(x)\n return x.squeeze(-1)\n\n\n# ============================================================================\n# 3D GRAD-CAM IMPLEMENTATION\n# ============================================================================\n\nclass GradCAM3D:\n \"\"\"\n GRAD-CAM for 3D medical imaging\n \n Paper: Grad-CAM: Visual Explanations from Deep Networks \n via Gradient-based Localization (Selvaraju et al., 2017)\n \n Adapted for 3D volumetric data\n \"\"\"\n \n def __init__(self, model, target_layer_name):\n \"\"\"\n Args:\n model: Your ResNet3D model\n target_layer_name: Name of layer to visualize (e.g., 'layer3')\n \"\"\"\n self.model = model\n self.model.eval()\n \n # Find target layer\n self.target_layer = None\n for name, module in model.named_modules():\n if name == target_layer_name:\n self.target_layer = module\n break\n \n if self.target_layer is None:\n raise ValueError(f\"Layer {target_layer_name} not found in model\")\n \n # Storage for activations and gradients\n self.activations = None\n self.gradients = None\n \n # Register hooks\n self.target_layer.register_forward_hook(self._save_activation)\n self.target_layer.register_full_backward_hook(self._save_gradient)\n \n def _save_activation(self, module, input, output):\n \"\"\"Hook to save forward activations\"\"\"\n self.activations = output.detach()\n \n def _save_gradient(self, module, grad_input, grad_output):\n \"\"\"Hook to save backward gradients\"\"\"\n self.gradients = grad_output[0].detach()\n \n def generate_cam(self, input_tensor, target_class=None):\n \"\"\"\n Generate GRAD-CAM heatmap\n \n Args:\n input_tensor: Input volume (1, C, D, H, W)\n target_class: Not used for binary classification\n \n Returns:\n cam: 3D heatmap (D, H, W) in range [0, 1]\n \"\"\"\n # Forward pass\n self.model.zero_grad()\n output = self.model(input_tensor)\n \n # Backward pass\n output.backward()\n \n # Get activations and gradients\n activations = self.activations # (1, C, D', H', W')\n gradients = self.gradients # (1, C, D', H', W')\n \n # Global average pooling of gradients (channel importance weights)\n weights = torch.mean(gradients, dim=(2, 3, 4), keepdim=True) # (1, C, 1, 1, 1)\n \n # Weighted combination of activation maps\n cam = torch.sum(weights * activations, dim=1, keepdim=True) # (1, 1, D', H', W')\n \n # Apply ReLU (only positive contributions)\n cam = F.relu(cam)\n \n # Remove batch and channel dimensions\n cam = cam.squeeze() # (D', H', W')\n \n # Normalize to [0, 1]\n cam = cam - cam.min()\n if cam.max() > 0:\n cam = cam / cam.max()\n \n return cam.cpu().numpy()\n \n def generate_cam_upsampled(self, input_tensor, target_size):\n \"\"\"\n Generate GRAD-CAM and upsample to original input size\n \n Args:\n input_tensor: Input volume (1, C, D, H, W)\n target_size: Tuple (D, H, W) to resize to\n \n Returns:\n cam: Upsampled 3D heatmap (D, H, W)\n \"\"\"\n cam = self.generate_cam(input_tensor)\n \n # Upsample to original size using trilinear interpolation\n cam_tensor = torch.from_numpy(cam).unsqueeze(0).unsqueeze(0) # (1, 1, D', H', W')\n cam_upsampled = F.interpolate(\n cam_tensor,\n size=target_size,\n mode='trilinear',\n align_corners=False\n )\n \n return cam_upsampled.squeeze().numpy()\n\n\n# ============================================================================\n# VISUALIZATION UTILITIES\n# ============================================================================\n\ndef apply_colormap(heatmap, colormap='jet'):\n \"\"\"\n Apply colormap to heatmap\n \n Args:\n heatmap: 2D array [0, 1]\n colormap: matplotlib colormap name\n \n Returns:\n colored: RGB image (H, W, 3)\n \"\"\"\n cmap = cm.get_cmap(colormap)\n colored = cmap(heatmap)[:, :, :3] # Remove alpha channel\n return (colored * 255).astype(np.uint8)\n\n\ndef overlay_heatmap(image, heatmap, alpha=0.4, colormap='jet'):\n \"\"\"\n Overlay heatmap on grayscale image\n \n Args:\n image: 2D grayscale [0, 1]\n heatmap: 2D heatmap [0, 1]\n alpha: Overlay transparency\n colormap: Colormap name\n \n Returns:\n overlay: RGB overlay (H, W, 3)\n \"\"\"\n # Convert grayscale to RGB\n image_rgb = (image * 255).astype(np.uint8)\n image_rgb = cv2.cvtColor(image_rgb, cv2.COLOR_GRAY2RGB)\n \n # Apply colormap to heatmap\n heatmap_colored = apply_colormap(heatmap, colormap)\n \n # Blend\n overlay = cv2.addWeighted(image_rgb, 1-alpha, heatmap_colored, alpha, 0)\n \n return overlay\n\n\ndef visualize_3d_gradcam(volume, cam_3d, prediction, patient_id, save_dir, \n num_slices=16, alpha=0.4, colormap='jet'):\n \"\"\"\n Create multi-slice visualization with GRAD-CAM overlays\n \n Args:\n volume: 3D volume (D, H, W, C) - use channel 0 (bone window)\n cam_3d: 3D GRAD-CAM heatmap (D, H, W)\n prediction: Model prediction probability\n patient_id: Patient identifier\n save_dir: Where to save visualization\n num_slices: Number of slices to show\n alpha: Overlay transparency\n colormap: Heatmap colormap\n \"\"\"\n D, H, W = volume.shape[:3]\n \n # Select evenly spaced slices\n slice_indices = np.linspace(0, D-1, num_slices, dtype=int)\n \n # Create grid layout\n cols = 4\n rows = (num_slices + cols - 1) // cols\n \n fig, axes = plt.subplots(rows, cols, figsize=(16, rows*4))\n axes = axes.flatten() if num_slices > 1 else [axes]\n \n for idx, slice_idx in enumerate(slice_indices):\n # Get slice from volume (use bone window - channel 0)\n image_slice = volume[slice_idx, :, :, 0]\n \n # Get corresponding heatmap slice\n heatmap_slice = cam_3d[slice_idx, :, :]\n \n # Create overlay\n overlay = overlay_heatmap(image_slice, heatmap_slice, alpha, colormap)\n \n # Plot\n axes[idx].imshow(overlay)\n axes[idx].set_title(f'Slice {slice_idx}/{D-1}', fontsize=10, fontweight='bold')\n axes[idx].axis('off')\n \n # Add intensity indicator\n max_intensity = heatmap_slice.max()\n if max_intensity > 0.7:\n marker = \"🔴\" # High attention\n elif max_intensity > 0.4:\n marker = \"🟡\" # Medium attention\n else:\n marker = \"⚪\" # Low attention\n \n axes[idx].text(5, 20, marker, fontsize=20)\n \n # Hide unused subplots\n for idx in range(num_slices, len(axes)):\n axes[idx].axis('off')\n \n # Add title with prediction\n risk_level = \"HIGH\" if prediction > 0.7 else \"MEDIUM\" if prediction > 0.3 else \"LOW\"\n color = 'red' if prediction > 0.7 else 'orange' if prediction > 0.3 else 'green'\n \n fig.suptitle(\n f'GRAD-CAM Analysis: {patient_id}\\n'\n f'Fracture Probability: {prediction:.1%} ({risk_level} RISK)',\n fontsize=14, fontweight='bold', color=color\n )\n \n plt.tight_layout()\n \n save_path = os.path.join(save_dir, f'{patient_id}_gradcam.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n return save_path\n\n\ndef create_max_intensity_projection(cam_3d, axis=0):\n \"\"\"\n Create Maximum Intensity Projection (MIP) of 3D heatmap\n \n Args:\n cam_3d: 3D heatmap (D, H, W)\n axis: Projection axis (0=axial, 1=coronal, 2=sagittal)\n \n Returns:\n mip: 2D projection\n \"\"\"\n return np.max(cam_3d, axis=axis)\n\n\n# ============================================================================\n# PREPROCESSING (from Step 4)\n# ============================================================================\n\ndef load_and_preprocess_test_patient(patient_folder):\n \"\"\"Load and preprocess test patient\"\"\"\n import pydicom\n from glob import glob\n from scipy.ndimage import zoom\n \n dicom_files = sorted(glob(os.path.join(patient_folder, \"*.dcm\")))\n \n if len(dicom_files) == 0:\n raise ValueError(\"No DICOM files found\")\n \n slices = []\n for dcm_file in dicom_files:\n try:\n ds = pydicom.dcmread(dcm_file, force=True)\n if hasattr(ds, 'decompress'):\n try:\n ds.decompress()\n except:\n pass\n if hasattr(ds, 'ImagePositionPatient') and hasattr(ds, 'PixelSpacing'):\n slices.append(ds)\n except:\n continue\n \n if len(slices) == 0:\n raise ValueError(\"No valid slices\")\n \n slices.sort(key=lambda x: float(x.ImagePositionPatient[2]))\n volume = np.stack([s.pixel_array for s in slices])\n metadata = slices[0]\n \n try:\n intercept = float(metadata.RescaleIntercept)\n slope = float(metadata.RescaleSlope)\n except:\n intercept = 0.0\n slope = 1.0\n \n volume_hu = volume.astype(np.float32) * slope + intercept\n \n def apply_windowing(vol_hu, center, width):\n lower = center - width // 2\n upper = center + width // 2\n windowed = np.clip(vol_hu, lower, upper)\n normalized = (windowed - lower) / (upper - lower)\n return normalized.astype(np.float32)\n \n bone = apply_windowing(volume_hu, 400, 1800)\n soft = apply_windowing(volume_hu, 40, 400)\n wide = apply_windowing(volume_hu, 400, 4000)\n \n multi_channel = np.stack([bone, soft, wide], axis=-1)\n \n target_shape = (64, 224, 224)\n current_shape = multi_channel.shape[:3]\n resize_factor = np.array(target_shape) / np.array(current_shape)\n resize_factor = np.append(resize_factor, 1.0)\n \n resampled = zoom(multi_channel, resize_factor, order=1).astype(np.float32)\n \n final = np.zeros((*target_shape, 3), dtype=np.float32)\n min_d = min(resampled.shape[0], target_shape[0])\n min_h = min(resampled.shape[1], target_shape[1])\n min_w = min(resampled.shape[2], target_shape[2])\n final[:min_d, :min_h, :min_w, :] = resampled[:min_d, :min_h, :min_w, :]\n \n return final\n\n\n# ============================================================================\n# MAIN PROCESSING\n# ============================================================================\n\ndef process_patients_with_gradcam(model, gradcam, test_patients, num_patients=20):\n \"\"\"\n Process test patients and generate GRAD-CAM visualizations\n \"\"\"\n \n print(f\"\\n{'='*80}\")\n print(\"🎨 GENERATING GRAD-CAM VISUALIZATIONS\")\n print(f\"{'='*80}\")\n \n results = []\n heatmap_dir = os.path.join(CONFIG['output_dir'], 'heatmaps')\n \n patients_to_process = test_patients[:num_patients]\n \n for patient_id in tqdm(patients_to_process, desc=\"Processing patients\"):\n patient_folder = os.path.join(CONFIG['test_dir'], patient_id)\n \n try:\n # Load and preprocess\n volume = load_and_preprocess_test_patient(patient_folder)\n \n # Convert to tensor\n volume_tensor = torch.from_numpy(volume).permute(3, 0, 1, 2).unsqueeze(0).float().to(device)\n \n # Get prediction\n with torch.no_grad():\n logit = model(volume_tensor)\n prediction = torch.sigmoid(logit).cpu().item()\n \n # Generate GRAD-CAM\n volume_tensor.requires_grad = True\n cam_3d = gradcam.generate_cam_upsampled(\n volume_tensor,\n target_size=(64, 224, 224)\n )\n \n # Visualize\n save_path = visualize_3d_gradcam(\n volume=volume,\n cam_3d=cam_3d,\n prediction=prediction,\n patient_id=patient_id,\n save_dir=heatmap_dir,\n num_slices=CONFIG['num_slices_to_show'],\n alpha=CONFIG['alpha'],\n colormap=CONFIG['colormap']\n )\n \n # Compute saliency metrics\n max_saliency = float(cam_3d.max())\n mean_saliency = float(cam_3d.mean())\n saliency_volume = float((cam_3d > 0.5).sum() / cam_3d.size) # Fraction of highly salient voxels\n \n results.append({\n 'patient_id': patient_id,\n 'prediction': prediction,\n 'max_saliency': max_saliency,\n 'mean_saliency': mean_saliency,\n 'saliency_volume': saliency_volume,\n 'visualization_path': save_path\n })\n \n except Exception as e:\n print(f\"\\n ⚠️ Failed {patient_id}: {str(e)[:80]}\")\n continue\n \n # Save results\n results_df = pd.DataFrame(results)\n results_path = os.path.join(CONFIG['output_dir'], 'gradcam_results.csv')\n results_df.to_csv(results_path, index=False)\n \n print(f\"\\n ✓ Processed {len(results)}/{len(patients_to_process)} patients\")\n print(f\" ✓ Saved results: {results_path}\")\n \n return results_df\n\n\n# ============================================================================\n# HTML REPORT GENERATION\n# ============================================================================\n\ndef generate_html_report(results_df, output_path):\n \"\"\"Generate interactive HTML report\"\"\"\n \n html_template = \"\"\"\n \n \n \n GRAD-CAM Analysis Report\n \n \n \n

🔍 GRAD-CAM Explainable AI Report

\n \n
\n

Summary Statistics

\n
\n
\n
{total_patients}
\n
Patients Analyzed
\n
\n
\n
{mean_pred:.1%}
\n
Mean Fracture Probability
\n
\n
\n
{high_risk_count}
\n
High Risk Patients (>70%)
\n
\n
\n
\n \n

Patient-Level Analysis

\n {patient_sections}\n \n \n \"\"\"\n \n patient_template = \"\"\"\n
\n

{patient_id}

\n
\n
\n
{prediction:.1%}
\n
Fracture Probability
\n
\n
\n
{max_saliency:.3f}
\n
Max Saliency
\n
\n
\n
{saliency_volume:.1%}
\n
Salient Volume
\n
\n
\n \"GRAD-CAM\n
\n \"\"\"\n \n # Generate patient sections\n patient_sections = []\n for _, row in results_df.iterrows():\n risk_class = 'high-risk' if row['prediction'] > 0.7 else 'medium-risk' if row['prediction'] > 0.3 else 'low-risk'\n \n section = patient_template.format(\n patient_id=row['patient_id'],\n prediction=row['prediction'],\n max_saliency=row['max_saliency'],\n saliency_volume=row['saliency_volume'],\n risk_class=risk_class,\n viz_path=os.path.relpath(row['visualization_path'], CONFIG['output_dir'])\n )\n patient_sections.append(section)\n \n # Generate final HTML\n html = html_template.format(\n total_patients=len(results_df),\n mean_pred=results_df['prediction'].mean(),\n high_risk_count=(results_df['prediction'] > 0.7).sum(),\n patient_sections='\\n'.join(patient_sections)\n )\n \n with open(output_path, 'w') as f:\n f.write(html)\n \n print(f\" ✓ Saved HTML report: {output_path}\")\n\n\n# ============================================================================\n# MAIN EXECUTION\n# ============================================================================\n\ndef main():\n print(f\"\\n{'='*80}\")\n print(\"🚀 INITIALIZING GRAD-CAM\")\n print(f\"{'='*80}\")\n \n # Load model\n model_path = os.path.join(CONFIG['models_dir'], f\"fold{CONFIG['fold_to_use']}_best.pth\")\n \n if not os.path.exists(model_path):\n print(f\"\\n❌ ERROR: Model not found: {model_path}\")\n print(f\" Run Step 3 (training) first!\")\n return\n \n print(f\"\\n ✓ Loading model from fold {CONFIG['fold_to_use']}\")\n \n model = ResNet3D(\n in_channels=CONFIG['in_channels'],\n base_filters=CONFIG['base_filters'],\n dropout=CONFIG['dropout']\n ).to(device)\n \n checkpoint = torch.load(model_path, map_location=device, weights_only=False)\n model.load_state_dict(checkpoint['model_state_dict'])\n model.eval()\n \n print(f\" ✓ Model loaded (AUC: {checkpoint.get('auc', 0):.4f})\")\n \n # Initialize GRAD-CAM\n gradcam = GradCAM3D(model, CONFIG['target_layer'])\n print(f\" ✓ GRAD-CAM initialized on layer: {CONFIG['target_layer']}\")\n \n # Get test patients\n if not os.path.exists(CONFIG['test_dir']):\n print(f\"\\n❌ ERROR: Test directory not found: {CONFIG['test_dir']}\")\n return\n \n test_patients = [d for d in os.listdir(CONFIG['test_dir']) \n if os.path.isdir(os.path.join(CONFIG['test_dir'], d))]\n print(f\" ✓ Found {len(test_patients)} test patients\")\n \n # Process patients\n results_df = process_patients_with_gradcam(\n model=model,\n gradcam=gradcam,\n test_patients=test_patients,\n num_patients=CONFIG['num_patients_to_visualize']\n )\n \n # Generate HTML report\n if CONFIG['generate_html_report']:\n report_path = os.path.join(CONFIG['output_dir'], 'gradcam_report.html')\n generate_html_report(results_df, report_path)\n \n # Final summary\n print(f\"\\n{'='*80}\")\n print(\"✅ GRAD-CAM ANALYSIS COMPLETE!\")\n print(f\"{'='*80}\")\n \n print(f\"\\n📊 Results Summary:\")\n print(f\" • Total visualizations: {len(results_df)}\")\n print(f\" • Mean prediction: {results_df['prediction'].mean():.3f}\")\n print(f\" • Mean max saliency: {results_df['max_saliency'].mean():.3f}\")\n print(f\" • High attention patients: {(results_df['max_saliency'] > 0.7).sum()}\")\n \n print(f\"\\n📁 Output Files:\")\n print(f\" • Heatmaps: {CONFIG['output_dir']}/heatmaps/\")\n print(f\" • Results CSV: {CONFIG['output_dir']}/gradcam_results.csv\")\n print(f\" • HTML Report: {CONFIG['output_dir']}/gradcam_report.html\")\n \n print(f\"\\n💡 Insights:\")\n if len(results_df) > 0:\n high_attention = results_df[results_df['max_saliency'] > 0.7]\n if len(high_attention) > 0:\n print(f\" • {len(high_attention)} patients show high attention (max saliency > 0.7)\")\n print(f\" • Average prediction for high-attention cases: {high_attention['prediction'].mean():.1%}\")\n \n # Correlation analysis\n if len(results_df) > 5:\n corr = results_df[['prediction', 'max_saliency', 'saliency_volume']].corr()\n pred_sal_corr = corr.loc['prediction', 'max_saliency']\n print(f\" • Correlation (prediction vs max_saliency): {pred_sal_corr:.3f}\")\n if pred_sal_corr > 0.5:\n print(f\" → Strong positive correlation - model is confident when attention is focused\")\n elif pred_sal_corr < 0:\n print(f\" → Negative correlation - investigate further!\")\n \n print(f\"\\n🔬 Clinical Validation:\")\n print(f\" 1. Review HTML report with radiologist\")\n print(f\" 2. Check if attention focuses on:\")\n print(f\" ✓ Vertebral bodies (good)\")\n print(f\" ✓ Known fracture locations (excellent)\")\n print(f\" ✗ Artifacts or irrelevant regions (bad)\")\n print(f\" 3. Compare high-risk vs low-risk attention patterns\")\n print(f\" 4. Use insights to improve model or data quality\")\n \n print(f\"\\n⚙️ Advanced Options:\")\n print(f\" • Try different layers: 'layer1', 'layer2', 'layer3'\")\n print(f\" • Change colormap: 'hot', 'viridis', 'plasma'\")\n print(f\" • Increase num_slices_to_show for more detail\")\n print(f\" • Process all test patients (set num_patients_to_visualize=-1)\")\n \n print(f\"\\n🎯 Next Steps:\")\n print(f\" 1. Validate clinical relevance with domain experts\")\n print(f\" 2. Create attention pattern comparison plots\")\n print(f\" 3. Identify failure cases (low attention on fractures)\")\n print(f\" 4. Consider attention-guided data augmentation\")\n print(f\" 5. Document findings for publication/presentation\")\n \n print(f\"\\n{'='*80}\")\n\n\n# ============================================================================\n# ADDITIONAL ANALYSIS FUNCTIONS\n# ============================================================================\n\ndef compare_attention_patterns(results_df, save_dir):\n \"\"\"\n Compare attention patterns between high-risk and low-risk predictions\n \"\"\"\n print(f\"\\n{'='*80}\")\n print(\"📊 ATTENTION PATTERN ANALYSIS\")\n print(f\"{'='*80}\")\n \n if len(results_df) < 5:\n print(\" ⚠️ Need at least 5 patients for comparison\")\n return\n \n # Split by risk level\n high_risk = results_df[results_df['prediction'] > 0.7]\n low_risk = results_df[results_df['prediction'] < 0.3]\n \n fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n \n # Max saliency distribution\n if len(high_risk) > 0 and len(low_risk) > 0:\n axes[0, 0].hist(high_risk['max_saliency'], bins=20, alpha=0.7, \n label='High Risk', color='red', edgecolor='black')\n axes[0, 0].hist(low_risk['max_saliency'], bins=20, alpha=0.7, \n label='Low Risk', color='green', edgecolor='black')\n axes[0, 0].set_title('Max Saliency Distribution', fontweight='bold')\n axes[0, 0].set_xlabel('Max Saliency')\n axes[0, 0].set_ylabel('Count')\n axes[0, 0].legend()\n axes[0, 0].grid(True, alpha=0.3)\n \n # Mean saliency distribution\n if len(high_risk) > 0 and len(low_risk) > 0:\n axes[0, 1].hist(high_risk['mean_saliency'], bins=20, alpha=0.7, \n label='High Risk', color='red', edgecolor='black')\n axes[0, 1].hist(low_risk['mean_saliency'], bins=20, alpha=0.7, \n label='Low Risk', color='green', edgecolor='black')\n axes[0, 1].set_title('Mean Saliency Distribution', fontweight='bold')\n axes[0, 1].set_xlabel('Mean Saliency')\n axes[0, 1].set_ylabel('Count')\n axes[0, 1].legend()\n axes[0, 1].grid(True, alpha=0.3)\n \n # Saliency volume distribution\n if len(high_risk) > 0 and len(low_risk) > 0:\n axes[0, 2].hist(high_risk['saliency_volume'], bins=20, alpha=0.7, \n label='High Risk', color='red', edgecolor='black')\n axes[0, 2].hist(low_risk['saliency_volume'], bins=20, alpha=0.7, \n label='Low Risk', color='green', edgecolor='black')\n axes[0, 2].set_title('Salient Volume Distribution', fontweight='bold')\n axes[0, 2].set_xlabel('Salient Volume Fraction')\n axes[0, 2].set_ylabel('Count')\n axes[0, 2].legend()\n axes[0, 2].grid(True, alpha=0.3)\n \n # Scatter plots\n axes[1, 0].scatter(results_df['prediction'], results_df['max_saliency'], \n alpha=0.6, c=results_df['prediction'], cmap='RdYlGn_r')\n axes[1, 0].set_title('Prediction vs Max Saliency', fontweight='bold')\n axes[1, 0].set_xlabel('Fracture Probability')\n axes[1, 0].set_ylabel('Max Saliency')\n axes[1, 0].grid(True, alpha=0.3)\n \n axes[1, 1].scatter(results_df['prediction'], results_df['mean_saliency'], \n alpha=0.6, c=results_df['prediction'], cmap='RdYlGn_r')\n axes[1, 1].set_title('Prediction vs Mean Saliency', fontweight='bold')\n axes[1, 1].set_xlabel('Fracture Probability')\n axes[1, 1].set_ylabel('Mean Saliency')\n axes[1, 1].grid(True, alpha=0.3)\n \n axes[1, 2].scatter(results_df['prediction'], results_df['saliency_volume'], \n alpha=0.6, c=results_df['prediction'], cmap='RdYlGn_r')\n axes[1, 2].set_title('Prediction vs Salient Volume', fontweight='bold')\n axes[1, 2].set_xlabel('Fracture Probability')\n axes[1, 2].set_ylabel('Salient Volume Fraction')\n axes[1, 2].grid(True, alpha=0.3)\n \n plt.suptitle('Attention Pattern Comparison Analysis', fontsize=14, fontweight='bold')\n plt.tight_layout()\n \n save_path = os.path.join(save_dir, 'attention_pattern_analysis.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved analysis: {save_path}\")\n \n # Statistical summary\n print(f\"\\n Statistical Summary:\")\n if len(high_risk) > 0:\n print(f\" High Risk (n={len(high_risk)}):\")\n print(f\" • Max saliency: {high_risk['max_saliency'].mean():.3f} ± {high_risk['max_saliency'].std():.3f}\")\n print(f\" • Mean saliency: {high_risk['mean_saliency'].mean():.3f} ± {high_risk['mean_saliency'].std():.3f}\")\n print(f\" • Salient volume: {high_risk['saliency_volume'].mean():.1%} ± {high_risk['saliency_volume'].std():.1%}\")\n \n if len(low_risk) > 0:\n print(f\" Low Risk (n={len(low_risk)}):\")\n print(f\" • Max saliency: {low_risk['max_saliency'].mean():.3f} ± {low_risk['max_saliency'].std():.3f}\")\n print(f\" • Mean saliency: {low_risk['mean_saliency'].mean():.3f} ± {low_risk['mean_saliency'].std():.3f}\")\n print(f\" • Salient volume: {low_risk['saliency_volume'].mean():.1%} ± {low_risk['saliency_volume'].std():.1%}\")\n\n\ndef create_summary_statistics_plot(results_df, save_dir):\n \"\"\"Create comprehensive summary statistics visualization\"\"\"\n \n fig = plt.figure(figsize=(16, 10))\n gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)\n \n # 1. Prediction distribution\n ax1 = fig.add_subplot(gs[0, :2])\n ax1.hist(results_df['prediction'], bins=30, alpha=0.7, color='steelblue', edgecolor='black')\n ax1.axvline(results_df['prediction'].mean(), color='red', linestyle='--', \n linewidth=2, label=f\"Mean: {results_df['prediction'].mean():.3f}\")\n ax1.axvline(results_df['prediction'].median(), color='orange', linestyle='--', \n linewidth=2, label=f\"Median: {results_df['prediction'].median():.3f}\")\n ax1.set_title('Prediction Distribution', fontweight='bold', fontsize=12)\n ax1.set_xlabel('Fracture Probability')\n ax1.set_ylabel('Frequency')\n ax1.legend()\n ax1.grid(True, alpha=0.3)\n \n # 2. Risk category pie chart\n ax2 = fig.add_subplot(gs[0, 2])\n risk_counts = [\n (results_df['prediction'] < 0.3).sum(),\n ((results_df['prediction'] >= 0.3) & (results_df['prediction'] <= 0.7)).sum(),\n (results_df['prediction'] > 0.7).sum()\n ]\n colors = ['#27ae60', '#f39c12', '#e74c3c']\n labels = [f'Low\\n({risk_counts[0]})', f'Medium\\n({risk_counts[1]})', f'High\\n({risk_counts[2]})']\n ax2.pie(risk_counts, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90)\n ax2.set_title('Risk Distribution', fontweight='bold', fontsize=12)\n \n # 3-5. Saliency metrics\n ax3 = fig.add_subplot(gs[1, 0])\n ax3.hist(results_df['max_saliency'], bins=25, alpha=0.7, color='coral', edgecolor='black')\n ax3.set_title('Max Saliency', fontweight='bold')\n ax3.set_xlabel('Value')\n ax3.set_ylabel('Frequency')\n ax3.grid(True, alpha=0.3)\n \n ax4 = fig.add_subplot(gs[1, 1])\n ax4.hist(results_df['mean_saliency'], bins=25, alpha=0.7, color='lightgreen', edgecolor='black')\n ax4.set_title('Mean Saliency', fontweight='bold')\n ax4.set_xlabel('Value')\n ax4.set_ylabel('Frequency')\n ax4.grid(True, alpha=0.3)\n \n ax5 = fig.add_subplot(gs[1, 2])\n ax5.hist(results_df['saliency_volume'], bins=25, alpha=0.7, color='plum', edgecolor='black')\n ax5.set_title('Salient Volume Fraction', fontweight='bold')\n ax5.set_xlabel('Value')\n ax5.set_ylabel('Frequency')\n ax5.grid(True, alpha=0.3)\n \n # 6. Correlation heatmap\n ax6 = fig.add_subplot(gs[2, :2])\n corr_matrix = results_df[['prediction', 'max_saliency', 'mean_saliency', 'saliency_volume']].corr()\n im = ax6.imshow(corr_matrix, cmap='coolwarm', aspect='auto', vmin=-1, vmax=1)\n ax6.set_xticks(range(len(corr_matrix.columns)))\n ax6.set_yticks(range(len(corr_matrix.columns)))\n ax6.set_xticklabels(['Prediction', 'Max Sal.', 'Mean Sal.', 'Sal. Vol.'], rotation=45, ha='right')\n ax6.set_yticklabels(['Prediction', 'Max Sal.', 'Mean Sal.', 'Sal. Vol.'])\n \n # Annotate correlation values\n for i in range(len(corr_matrix)):\n for j in range(len(corr_matrix)):\n text = ax6.text(j, i, f'{corr_matrix.iloc[i, j]:.2f}',\n ha=\"center\", va=\"center\", color=\"black\", fontsize=10, fontweight='bold')\n \n ax6.set_title('Feature Correlation Matrix', fontweight='bold', fontsize=12)\n plt.colorbar(im, ax=ax6, label='Correlation')\n \n # 7. Summary statistics table\n ax7 = fig.add_subplot(gs[2, 2])\n ax7.axis('off')\n \n stats_text = f\"\"\"\n SUMMARY STATISTICS\n {'='*30}\n \n Total Patients: {len(results_df)}\n \n Predictions:\n • Mean: {results_df['prediction'].mean():.3f}\n • Std: {results_df['prediction'].std():.3f}\n • Min: {results_df['prediction'].min():.3f}\n • Max: {results_df['prediction'].max():.3f}\n \n Max Saliency:\n • Mean: {results_df['max_saliency'].mean():.3f}\n • Std: {results_df['max_saliency'].std():.3f}\n \n Mean Saliency:\n • Mean: {results_df['mean_saliency'].mean():.3f}\n • Std: {results_df['mean_saliency'].std():.3f}\n \n Salient Volume:\n • Mean: {results_df['saliency_volume'].mean():.1%}\n • Std: {results_df['saliency_volume'].std():.1%}\n \"\"\"\n \n ax7.text(0.1, 0.5, stats_text, fontsize=9, family='monospace',\n verticalalignment='center',\n bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.3))\n \n plt.suptitle('GRAD-CAM Comprehensive Summary Statistics', fontsize=14, fontweight='bold')\n \n save_path = os.path.join(save_dir, 'summary_statistics.png')\n plt.savefig(save_path, dpi=150, bbox_inches='tight')\n plt.close()\n \n print(f\" ✓ Saved summary statistics: {save_path}\")\n\n\nif __name__ == \"__main__\":\n try:\n # Run main analysis\n main()\n \n # Load results for additional analysis\n results_path = os.path.join(CONFIG['output_dir'], 'gradcam_results.csv')\n if os.path.exists(results_path):\n results_df = pd.read_csv(results_path)\n \n # Generate additional visualizations\n print(f\"\\n{'='*80}\")\n print(\"📈 GENERATING ADDITIONAL ANALYSES\")\n print(f\"{'='*80}\")\n \n compare_attention_patterns(results_df, CONFIG['output_dir'])\n create_summary_statistics_plot(results_df, CONFIG['output_dir'])\n \n print(f\"\\n✅ All analyses complete!\")\n \n except KeyboardInterrupt:\n print(f\"\\n\\n⚠️ Analysis interrupted by user\")\n except Exception as e:\n print(f\"\\n{'='*80}\")\n print(\"❌ ERROR\")\n print(f\"{'='*80}\")\n print(f\"\\nError: {str(e)}\")\n import traceback\n traceback.print_exc()\n \n print(f\"\\n💡 Troubleshooting:\")\n print(f\" 1. Ensure Step 3 (training) completed successfully\")\n print(f\" 2. Check that model checkpoint exists in {CONFIG['models_dir']}\")\n print(f\" 3. Verify test data is accessible at {CONFIG['test_dir']}\")\n print(f\" 4. Try with fewer patients (num_patients_to_visualize=5)\")\n print(f\" 5. If memory error, reduce num_slices_to_show\")\n print(f\" 6. For gradient issues, ensure model is in eval mode\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}