{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":91249,"databundleVersionId":11294684,"sourceType":"competition"}],"dockerImageVersionId":31012,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-04-16T16:10:40.866903Z","iopub.execute_input":"2025-04-16T16:10:40.867662Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -*- coding: utf-8 -*-\n\"\"\"\nKaggle Notebook: BYU Flagellar Motor Localization (Revised for Directory Structure & Multi-Motor)\n\nGoal: Detect the presence and (x, y, z) coordinates of bacterial\n      flagellar motors in 3D cryo-ET tomograms stored as slice directories.\n\nInspired Concepts: Analogies from Non-Linear Electrodynamics / String Theory\n    - Tomogram Intensity: Scalar field (like potential)\n    - Gradient (∇I): Vector field (analogous to Electric Field E)\n    - Structure Tensor (∇I ⊗ ∇I averaged): Symmetric Tensor (analogous to Stress Tensor T_µν or Metric Perturbation S_µν)\n    - Hessian (∇²I): Tensor capturing curvature (related to field changes)\n    - Invariants: Rotationally invariant features derived from gradient/structure tensor/Hessian (analogous to Lorentz invariants x, y)\n    - Anisotropic Diffusion/Filtering: Guided propagation of information (analogous to wave propagation influenced by background fields/metrics)\n    - Multi-Metric Idea: Using different analysis scales/features (like g_µν vs G_µν) - e.g., global search vs. local refinement.\n\nMODIFICATION: This version SKIPS training data generation and model training entirely.\n            It focuses prediction ONLY on specified IDs ['tomo_003acc', 'tomo_00e047', 'tomo_00e463']\n            using NO trained model (will predict NO_MOTOR_COORD).\n            The final output is a submission file template for ALL test IDs found.\n\"\"\"\n\n# %% [markdown]\n# # 1. Setup and Imports\n#\n# Load necessary libraries and define constants.\n\n# %%\n# Attempt installation only if needed (e.g., in Kaggle environment)\ntry:\n    import cv2\n    import mrcfile # Though likely unused for loading now\n    print(\"Libraries opencv-python, mrcfile seem available.\")\nexcept ImportError:\n    print(\"Installing required libraries: opencv-python, mrcfile\")\n    # Use %pip instead of !pip in notebooks for better integration\n    %pip install opencv-python mrcfile --quiet\n    import cv2\n    import mrcfile\n    print(\"Libraries installed.\")\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport cv2 # Import OpenCV\n# import mrcfile # Keep if other parts might use it, but not for loading slices\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom scipy import ndimage as ndi\nfrom skimage.feature import structure_tensor, structure_tensor_eigenvalues, hessian_matrix, hessian_matrix_eigvals\nfrom skimage.filters import gaussian, median\nfrom skimage.measure import regionprops, label\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\nfrom sklearn.metrics import fbeta_score, mean_squared_error, make_scorer\nimport gc # Garbage collection\nimport time # For timing\n\n# Optional: Deep Learning Libraries (if used)\n# import torch\n# import torch.nn as nn\n# from torch.utils.data import Dataset, DataLoader\n# import monai # Medical Imaging AI library\n\n# Constants\nCOMPETITION_DIR = '/kaggle/input/byu-locating-bacterial-flagellar-motors-2025'\nTRAIN_DIR = os.path.join(COMPETITION_DIR, 'train')\nTEST_DIR = os.path.join(COMPETITION_DIR, 'test')\nTRAIN_LABELS_PATH = os.path.join(COMPETITION_DIR, 'train_labels.csv')\nSAMPLE_SUB_PATH = os.path.join(COMPETITION_DIR, 'sample_submission.csv')\n\n# Default Voxel spacing (used for test set or if lookup fails)\nDEFAULT_VOXEL_SPACING_A = 10.0\nprint(f\"Using default voxel spacing: {DEFAULT_VOXEL_SPACING_A} Angstroms/pixel (will try to read from labels for training)\")\n\n# Evaluation Threshold (Angstroms)\nDISTANCE_THRESHOLD_A = 1000.0\n# DISTANCE_THRESHOLD_PX will be calculated dynamically based on actual voxel spacing\n\n# F-beta Score Beta value\nF_BETA = 2.0\n\n# For submission: Coordinate value indicating no motor found\nNO_MOTOR_COORD = -1.0 # Ensure float\n\n# Seed for reproducibility\nSEED = 42\nnp.random.seed(SEED)\n\n# %% [markdown]\n# # 2. Load Data and Metadata\n#\n# Load the training labels and find the data directories. Check for consistency.\n\n# %%\nstart_time = time.time()\nprint(\"Loading labels...\")\ntry:\n    train_labels_df = pd.read_csv(TRAIN_LABELS_PATH)\n    print(f\"Training labels shape: {train_labels_df.shape}\")\n    print(train_labels_df.head())\nexcept FileNotFoundError:\n    print(f\"ERROR: Training labels file not found at {TRAIN_LABELS_PATH}\")\n    train_labels_df = pd.DataFrame() # Empty dataframe\n\n# --- Debug: List directories to confirm paths ---\nprint(\"\\n--- Directory Listing ---\")\nprint(f\"Listing {COMPETITION_DIR}:\")\ntry:\n    print(os.listdir(COMPETITION_DIR))\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {COMPETITION_DIR}\")\n\nprint(f\"\\nListing {TRAIN_DIR}:\")\ntry:\n    train_contents = os.listdir(TRAIN_DIR)\n    print(f\"  Found {len(train_contents)} items in train dir. First 10:\")\n    print(train_contents[:10])\n    if len(train_contents) > 10: print(\"  ...\")\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {TRAIN_DIR}\")\n    train_contents = []\n\nprint(f\"\\nListing {TEST_DIR}:\")\ntry:\n    test_contents = os.listdir(TEST_DIR)\n    print(f\"  Found {len(test_contents)} items in test dir. First 10:\")\n    print(test_contents[:10])\n    if len(test_contents) > 10: print(\"  ...\")\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {TEST_DIR}\")\n    test_contents = []\nprint(\"--- End Directory Listing ---\\n\")\n\n\n# --- Find Tomogram Directories ---\nprint(\"Finding tomogram directories...\")\ntrain_dirs = sorted(glob.glob(os.path.join(TRAIN_DIR, 'tomo_*')))\ntest_dirs = sorted(glob.glob(os.path.join(TEST_DIR, 'tomo_*')))\n\n# Extract UNIQUE tomo_ids present in the directories found\ntrain_ids_found = sorted([os.path.basename(d) for d in train_dirs if os.path.isdir(d)]) # Ensure it's a directory\ntest_ids_found = sorted([os.path.basename(d) for d in test_dirs if os.path.isdir(d)]) # Ensure it's a directory\n\nprint(f\"Found {len(train_ids_found)} potential training tomogram directories.\")\nprint(f\"First 5 found: {train_ids_found[:5]}\")\nprint(f\"Found {len(test_ids_found)} potential testing tomogram directories.\")\nprint(f\"First 5 found: {test_ids_found[:5]}\")\n\n# --- Filter Labels based on Found Directories ---\nif not train_labels_df.empty:\n    # Check if 'tomo_id' column exists\n    if 'tomo_id' not in train_labels_df.columns:\n        print(\"ERROR: 'tomo_id' column not found in labels CSV. Cannot proceed.\")\n        train_labels_df_filtered = pd.DataFrame()\n        unique_train_ids_with_data = []\n        voxel_spacing_map = {}\n    else:\n        train_labels_df_filtered = train_labels_df[train_labels_df['tomo_id'].isin(train_ids_found)].copy()\n        print(f\"\\nOriginal label rows: {len(train_labels_df)}\")\n        print(f\"Label rows after filtering by found train directories: {len(train_labels_df_filtered)}\")\n\n        # Get unique tomogram IDs that have both labels AND a found directory\n        unique_train_ids_with_data = sorted(train_labels_df_filtered['tomo_id'].unique())\n        print(f\"Number of unique train tomograms with labels AND data: {len(unique_train_ids_with_data)}\")\n\n        # Create map for Voxel Spacing from labels (use first entry per tomo_id)\n        if 'Voxel spacing' in train_labels_df_filtered.columns and not train_labels_df_filtered.empty:\n            voxel_spacing_map = train_labels_df_filtered.drop_duplicates(subset='tomo_id').set_index('tomo_id')['Voxel spacing'].to_dict()\n            print(f\"Example voxel spacing: {list(voxel_spacing_map.items())[:5]}\")\n        else:\n             print(\"WARNING: 'Voxel spacing' column not found or no valid labels. Cannot create voxel spacing map.\")\n             voxel_spacing_map = {}\n\nelse:\n    print(\"WARNING: Training labels dataframe is empty. Cannot filter or create maps.\")\n    train_labels_df_filtered = pd.DataFrame()\n    unique_train_ids_with_data = []\n    voxel_spacing_map = {}\n\n# Map tomo_id to DIRECTORY path (even if labels are missing)\ntrain_dir_map = {os.path.basename(d): d for d in train_dirs if os.path.isdir(d)}\ntest_dir_map = {os.path.basename(d): d for d in test_dirs if os.path.isdir(d)}\n\n\nprint(f\"Data loading setup took {time.time() - start_time:.2f} seconds.\")\n\nif not unique_train_ids_with_data and not train_labels_df.empty:\n    print(\"\\n\\nCRITICAL WARNING: No intersection between labels and found training directories. Training cannot proceed.\")\nelif not train_ids_found and not test_ids_found:\n     print(\"\\n\\nCRITICAL WARNING: No training or testing directories found. Check data paths.\")\n\n# %% [markdown]\n# # 3. Exploratory Data Analysis (EDA) & Concept Visualization\n#\n# Understand the data distribution, visualize tomograms and labels, and see how our physics-inspired concepts manifest.\n# (Keeping these utility functions defined, even if visualization is partially skipped)\n\n# %%\n# --- Tomogram Loading Function ---\ndef load_tomogram(tomo_id, dir_map):\n    \"\"\"Loads a tomogram by reading and stacking slices from a directory.\"\"\"\n    if isinstance(dir_map, str):\n        directory_path = dir_map\n        if not os.path.isdir(directory_path): return None\n    else:\n        directory_path = dir_map.get(tomo_id)\n        if not directory_path or not os.path.isdir(directory_path): return None\n\n    slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.jpg')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.png')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.tif')))\n    if not slice_files: return None\n\n    def get_slice_number(filepath):\n        try:\n            filename = os.path.basename(filepath)\n            num_str = filename.split('slice_')[1].split('.')[0]\n            return int(num_str)\n        except: return -1\n\n    valid_slice_files = [(f, get_slice_number(f)) for f in slice_files]\n    valid_slice_files = [(f, num) for f, num in valid_slice_files if num != -1]\n    if not valid_slice_files: return None\n    valid_slice_files.sort(key=lambda item: item[1])\n    slice_files_sorted = [f for f, num in valid_slice_files]\n\n    try:\n        first_slice = cv2.imread(slice_files_sorted[0], cv2.IMREAD_GRAYSCALE)\n        if first_slice is None: raise IOError(f\"cv2.imread failed for {slice_files_sorted[0]}\")\n        height, width = first_slice.shape\n        num_slices = len(slice_files_sorted)\n        dtype = first_slice.dtype\n    except Exception as e:\n        print(f\"Error reading first slice {slice_files_sorted[0]}: {e}\"); return None\n\n    tomogram_data = np.zeros((num_slices, height, width), dtype=dtype)\n    tomogram_data[0, :, :] = first_slice\n\n    for i in range(1, num_slices):\n        try:\n            slice_img = cv2.imread(slice_files_sorted[i], cv2.IMREAD_GRAYSCALE)\n            if slice_img is None: raise IOError(f\"cv2.imread failed for {slice_files_sorted[i]}\")\n            if slice_img.shape != (height, width): return None\n            tomogram_data[i, :, :] = slice_img\n        except Exception as e:\n            print(f\"Error reading slice {slice_files_sorted[i]}: {e}\"); return None\n    return tomogram_data\n\n# --- Plotting Function ---\ndef plot_slices_with_label(tomo_data, voxel_spacing, label_coords_A=None, title=\"Tomogram Slices\"):\n    \"\"\"Plots central slices along each axis, optionally marking the label(s).\"\"\"\n    if tomo_data is None: print(\"No data to plot.\"); return\n    if voxel_spacing <= 0: print(\"Invalid voxel spacing.\"); return\n    shape = tomo_data.shape\n    if not (len(shape) == 3 and all(s > 0 for s in shape)): print(f\"Invalid shape {shape}.\"); return\n\n    center_z, center_y, center_x = shape[0] // 2, shape[1] // 2, shape[2] // 2\n    fig, axes = plt.subplots(1, 3, figsize=(15, 5)); fig.suptitle(title, fontsize=16)\n    try:\n        axes[0].imshow(tomo_data[center_z, :, :], cmap='gray'); axes[0].set_title(f'Z Slice (Z={center_z})')\n        axes[1].imshow(tomo_data[:, center_y, :], cmap='gray', aspect=shape[0]/shape[2] if shape[2]>0 else 1); axes[1].set_title(f'Y Slice (Y={center_y})')\n        axes[2].imshow(tomo_data[:, :, center_x], cmap='gray', aspect=shape[0]/shape[1] if shape[1]>0 else 1); axes[2].set_title(f'X Slice (X={center_x})')\n    except IndexError as e: print(f\"Error plotting slices: {e}\"); plt.close(fig); return\n\n    if label_coords_A is not None and len(label_coords_A) > 0:\n        is_multiple = isinstance(label_coords_A[0], (list, np.ndarray))\n        if not is_multiple: label_coords_A = [label_coords_A]\n        plotted_legend = False\n        for i, motor_A in enumerate(label_coords_A):\n            if len(motor_A) < 3 or motor_A[0] == NO_MOTOR_COORD: continue\n            lz_A, ly_A, lx_A = motor_A\n            lz_px, ly_px, lx_px = int(lz_A / voxel_spacing), int(ly_A / voxel_spacing), int(lx_A / voxel_spacing)\n            if not (0 <= lz_px < shape[0] and 0 <= ly_px < shape[1] and 0 <= lx_px < shape[2]): continue\n            label_text = f'Motor {i+1}' if len(label_coords_A) > 1 else 'Motor'\n            z_thresh=max(5, shape[0]*0.05); y_thresh=max(5, shape[1]*0.05); x_thresh=max(5, shape[2]*0.05)\n            try:\n                lbl = None\n                if not plotted_legend: lbl = label_text; plotted_legend=True\n                if abs(lz_px - center_z) < z_thresh: axes[0].plot(lx_px, ly_px, 'ro', ms=8, label=lbl)\n                if abs(ly_px - center_y) < y_thresh: axes[1].plot(lx_px, lz_px, 'ro', ms=8)\n                if abs(lx_px - center_x) < x_thresh: axes[2].plot(ly_px, lz_px, 'ro', ms=8)\n                if lbl: lbl=None\n            except Exception as plot_e: print(f\"Error plotting marker: {plot_e}\")\n    if plotted_legend: axes[0].legend()\n    plt.tight_layout(rect=[0, 0.03, 1, 0.95]); plt.show()\n\n# --- Visualize a sample tomogram ---\n# (Keep this visualization, it's quick and useful)\nif unique_train_ids_with_data:\n    sample_tomo_id = unique_train_ids_with_data[0]\n    print(f\"\\nVisualizing sample tomogram: {sample_tomo_id}\")\n    sample_data = load_tomogram(sample_tomo_id, train_dir_map)\n    sample_labels = train_labels_df_filtered[train_labels_df_filtered['tomo_id'] == sample_tomo_id]\n    if sample_data is not None and not sample_labels.empty:\n        motor_locations_A = sample_labels[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].values.tolist()\n        voxel_spacing = sample_labels.iloc[0]['Voxel spacing']\n        print(f\"Sample Tomogram Shape: {sample_data.shape}, Voxel Spacing: {voxel_spacing} A/px\")\n        print(f\"Motor Location(s) (Angstroms): {motor_locations_A}\")\n        plot_slices_with_label(sample_data, voxel_spacing, motor_locations_A, title=f\"{sample_tomo_id}\")\n        plt.figure(figsize=(10, 4)); plt.hist(sample_data.flatten(), bins=100, color='blue', alpha=0.7); plt.title(f'Intensity Distribution {sample_tomo_id}'); plt.xlabel('Intensity'); plt.ylabel('Frequency'); plt.grid(True, alpha=0.3); plt.show()\n        del sample_data, sample_labels; gc.collect()\n    else: print(f\"Failed to load/find labels for sample {sample_tomo_id}\")\nelse: print(\"\\nSkipping Visualization: No valid training tomograms found.\")\n\n# %% [markdown]\n# ### 3.1 Visualizing Physics-Inspired Concepts\n# (Skip this visualization as it relies on patch extraction and calculations that mirror training data gen)\n\n# %%\n# --- Patch Extraction Function (Keep definition for prediction step) ---\ndef get_local_patch(tomo_data, center_px, patch_size_px=64):\n    \"\"\"Extracts a 3D patch centered at center_px. Patch size is in pixels.\"\"\"\n    if tomo_data is None or center_px is None: return None, None\n    center_px = np.round(center_px).astype(int); z, y, x = center_px\n    shape = tomo_data.shape\n    if not (0 <= z < shape[0] and 0 <= y < shape[1] and 0 <= x < shape[2]): return None, None\n    half_size = patch_size_px // 2\n    z_start, z_end = max(0, z - half_size), min(shape[0], z + half_size)\n    y_start, y_end = max(0, y - half_size), min(shape[1], y + half_size)\n    x_start, x_end = max(0, x - half_size), min(shape[2], x + half_size)\n    dest_z_start = half_size - (z - z_start); dest_z_end = dest_z_start + (z_end - z_start)\n    dest_y_start = half_size - (y - y_start); dest_y_end = dest_y_start + (y_end - y_start)\n    dest_x_start = half_size - (x - x_start); dest_x_end = dest_x_start + (x_end - x_start)\n    patch = np.zeros((patch_size_px, patch_size_px, patch_size_px), dtype=tomo_data.dtype)\n    try:\n        src_slice = (slice(z_start, z_end), slice(y_start, y_end), slice(x_start, x_end))\n        dest_slice = (slice(dest_z_start, dest_z_end), slice(dest_y_start, dest_y_end), slice(dest_x_start, dest_x_end))\n        patch[dest_slice] = tomo_data[src_slice]\n        patch_center_coords_in_patch = np.round([z - z_start + dest_z_start, y - y_start + dest_y_start, x - x_start + dest_x_start]).astype(int)\n        if not (0 <= patch_center_coords_in_patch[0] < patch_size_px and 0 <= patch_center_coords_in_patch[1] < patch_size_px and 0 <= patch_center_coords_in_patch[2] < patch_size_px):\n             patch_center_coords_in_patch = [patch_size_px // 2] * 3\n        return patch, patch_center_coords_in_patch\n    except (ValueError, IndexError) as e: return None, None\n\n# --- Visualize physics concepts ---\nprint(\"\\nSkipping Physics Concept Visualization (requires patch processing).\")\n# (The code for visualization is removed here)\n\n# %% [markdown]\n# # 4. Data Preprocessing & Feature Engineering\n# (Keep function definitions as they are needed for prediction, even if models aren't trained)\n\n# %%\n# --- Simple Preprocessing ---\ndef preprocess_tomogram(data):\n    \"\"\"Normalizes and smooths the tomogram data.\"\"\"\n    if data is None: return None\n    if not np.issubdtype(data.dtype, np.floating): data = data.astype(np.float32)\n    mean, std = np.mean(data), np.std(data)\n    if std > 1e-6: normalized_data = (data - mean) / std\n    else: normalized_data = data - mean\n    smoothed_data = gaussian(normalized_data, sigma=1.5, mode='reflect', preserve_range=True, truncate=4.0)\n    return smoothed_data\n\n# --- Feature Extraction (Keep definition) ---\ndef extract_features_patch(patch):\n    \"\"\"Extracts a feature vector from a 3D patch.\"\"\"\n    if patch is None or patch.size == 0: return None\n    if not np.issubdtype(patch.dtype, np.floating): patch = patch.astype(np.float32)\n    features = []\n    EXPECTED_N_FEATURES = 15\n    try:\n        features.extend([np.mean(patch), np.std(patch), np.median(patch), np.min(patch), np.max(patch)])\n        patch_smooth_grad = gaussian(patch, sigma=1.0, mode='reflect', preserve_range=True, truncate=4.0)\n        patch_smooth_st = gaussian(patch, sigma=1.5, mode='reflect', preserve_range=True, truncate=4.0)\n        patch_smooth_hess = gaussian(patch, sigma=2.0, mode='reflect', preserve_range=True, truncate=4.0)\n        try:\n            grad_z, grad_y, grad_x = np.gradient(patch_smooth_grad)\n            grad_mag = np.sqrt(grad_z**2 + grad_y**2 + grad_x**2)\n            features.extend([np.mean(grad_mag), np.std(grad_mag), np.max(grad_mag)])\n        except Exception: features.extend([0.0] * 3)\n        try:\n            shape = patch_smooth_st.shape; slice_start = max(0, shape[0]//5); slice_end = max(slice_start+1, 4*shape[0]//5)\n            center_slice = slice(slice_start, slice_end)\n            if not (0<=center_slice.start<shape[0] and 0<center_slice.stop<=shape[0] and center_slice.start<center_slice.stop): center_patch_smooth = patch_smooth_st\n            else: center_patch_smooth = patch_smooth_st[center_slice, center_slice, center_slice]\n            if center_patch_smooth.size > 0:\n                S_elems_c = structure_tensor(center_patch_smooth, sigma=1.5, mode='reflect')\n                eigvals_S_c = structure_tensor_eigenvalues(S_elems_c); eigvals_S_c = np.sort(eigvals_S_c, axis=0)\n                l3, l2, l1 = eigvals_S_c[0], eigvals_S_c[1], eigvals_S_c[2]\n                den = l1 + l2 + l3 + 1e-9; coh = np.where(den > 1e-8, (l1 - l3) / den, 0.0)\n                features.extend([np.mean(l1), np.mean(l3), np.mean(coh), np.max(coh)])\n            else: features.extend([0.0] * 4)\n        except Exception: features.extend([0.0] * 4)\n        try:\n            shape_h = patch_smooth_hess.shape; slice_start_h = max(0, shape_h[0]//5); slice_end_h = max(slice_start_h+1, 4*shape_h[0]//5)\n            center_slice_h = slice(slice_start_h, slice_end_h)\n            if not (0<=center_slice_h.start<shape_h[0] and 0<center_slice_h.stop<=shape_h[0] and center_slice_h.start<center_slice_h.stop): center_patch_smooth_h = patch_smooth_hess\n            else: center_patch_smooth_h = patch_smooth_hess[center_slice_h, center_slice_h, center_slice_h]\n            if center_patch_smooth_h.size > 0:\n                h_matrix = hessian_matrix(center_patch_smooth_h, sigma=2.0, mode='reflect', use_gaussian_derivatives=False)\n                eigvals_H_c = hessian_matrix_eigvals(h_matrix); eigvals_H_c = np.sort(eigvals_H_c, axis=0)\n                h1, h2, h3 = eigvals_H_c[0], eigvals_H_c[1], eigvals_H_c[2]\n                features.extend([np.mean(h1), np.mean(h3), np.std(h1)])\n            else: features.extend([0.0] * 3)\n        except Exception: features.extend([0.0] * 3)\n    except Exception as outer_e: print(f\"Outer FE error: {outer_e}\"); return None\n    if len(features) != EXPECTED_N_FEATURES:\n        if len(features) < EXPECTED_N_FEATURES: features.extend([0.0] * (EXPECTED_N_FEATURES - len(features)))\n        else: features = features[:EXPECTED_N_FEATURES]\n        if len(features) != EXPECTED_N_FEATURES: return None\n    features_arr = np.array(features, dtype=np.float32)\n    if np.any(np.isnan(features_arr)) or np.any(np.isinf(features_arr)):\n        features_arr = np.nan_to_num(features_arr, nan=0.0, posinf=0.0, neginf=0.0)\n    if len(features_arr) != EXPECTED_N_FEATURES: return None\n    return features_arr\n\n\n# %% [markdown]\n# # 5. Model Development (Example: Patch Classification + Regression)\n#\n# <<< SKIPPING Training Data Generation >>>\n\n# %%\n# --- Generate Training Data ---\nPATCH_SIZE_PX = 64 # Keep constants defined\nEXPECTED_N_FEATURES = 15\n\nprint(\"\\nSKIPPING Training Data Generation as requested.\")\n\n# Define variables that would have been created as None or empty\nfeatures_list = []; labels_clf_list = []; labels_reg_list = []\nfailed_tomos = []\nclf_model = None # Crucial: Ensure models are None\nreg_model = None # Crucial: Ensure models are None\nX = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_clf = np.array([], dtype=np.int32)\ny_reg_all = np.array([]).reshape(0, 3)\nX_reg = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_reg = np.array([]).reshape(0, 3)\nX_train_clf = X_val_clf = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_train_clf = y_val_clf = np.array([], dtype=np.int32)\nX_train_reg = X_val_reg = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_train_reg = y_val_reg = np.array([]).reshape(0, 3)\n\ngc.collect() # Clean up memory in case any large objects were created before skip\n\n\n# %% [markdown]\n# ### 5.1 Train Models\n#\n# <<< SKIPPING Model Training >>>\n\n# %%\n# --- Split data and Train ---\nprint(\"\\nSKIPPING Model Training as training data generation was skipped.\")\n\n# Models remain None as set in the previous step\nprint(f\"\\nModels after training attempts:\")\nprint(f\"  Classifier model: {'Available' if clf_model is not None else 'Not Available'}\")\nprint(f\"  Regressor model: {'Available' if reg_model is not None else 'Not Available'}\")\n\n\n# %% [markdown]\n# # 6. Evaluation Metric Implementation\n# (Keep definition, although it won't be used effectively without trained models/predictions)\n\n# %%\ndef calculate_fbeta_and_distance(y_true_df, y_pred_df, voxel_spacing_map, default_voxel_spacing, beta=F_BETA, threshold_a=DISTANCE_THRESHOLD_A):\n    \"\"\"\n    Calculates the competition metric (F-beta score based on distance).\n    Handles multiple ground truth motors per tomogram.\n    \"\"\"\n    tp = 0; fp = 0; fn = 0\n    no_motor = float(NO_MOTOR_COORD)\n\n    if y_pred_df.empty: pred_map = {}\n    else:\n        y_pred_df_unique = y_pred_df.drop_duplicates(subset='tomo_id', keep='first')\n        pred_map = y_pred_df_unique.set_index('tomo_id')[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].apply(lambda row: row.tolist(), axis=1).to_dict()\n\n    if y_true_df.empty: true_grouped = {}; total_true_motors = 0\n    else:\n        y_true_motors_only = y_true_df[y_true_df['Motor axis 0'] != no_motor].copy()\n        true_grouped = y_true_motors_only.groupby('tomo_id')\n        total_true_motors = len(y_true_motors_only)\n\n    tp_final = 0; fp_final = 0\n    matched_true_motor_indices = {}\n\n    for tomo_id, pred_coords_a in pred_map.items():\n        if pred_coords_a[0] == no_motor: continue\n        true_motors_a = []; true_indices = []; is_true_motor_present = False\n        if tomo_id in true_grouped.groups:\n             true_df_tomo = true_grouped.get_group(tomo_id)\n             true_motors_a = true_df_tomo[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].values\n             true_indices = true_df_tomo.index.tolist(); is_true_motor_present = True\n        if not is_true_motor_present: fp_final += 1; continue\n        pred_vec = np.array(pred_coords_a)\n        found_match = False; best_match_true_idx = -1; min_dist = float('inf')\n        for i, true_motor_a in enumerate(true_motors_a):\n            dist = np.linalg.norm(true_motor_a - pred_vec)\n            original_idx = true_indices[i]\n            if dist <= threshold_a:\n                 already_matched = tomo_id in matched_true_motor_indices and original_idx in matched_true_motor_indices[tomo_id]\n                 if not already_matched and dist < min_dist:\n                      min_dist = dist; best_match_true_idx = original_idx; found_match = True\n        if found_match:\n            tp_final += 1\n            if tomo_id not in matched_true_motor_indices: matched_true_motor_indices[tomo_id] = set()\n            matched_true_motor_indices[tomo_id].add(best_match_true_idx)\n        else: fp_final += 1\n    fn_final = total_true_motors - tp_final\n    f_beta_numerator = (1 + beta**2) * tp_final\n    f_beta_denominator = (1 + beta**2) * tp_final + (beta**2 * fn_final) + fp_final\n    if f_beta_denominator == 0: f_beta_score = 1.0 if total_true_motors == 0 else 0.0\n    else: f_beta_score = f_beta_numerator / f_beta_denominator\n    stats = {'TP': tp_final, 'FP': fp_final, 'FN': fn_final, 'Total True': total_true_motors}\n    # Print score even if it's based on default predictions\n    print(f\"Evaluation Stats (based on default predictions): TP={tp_final}, FP={fp_final}, FN={fn_final} (Total True={total_true_motors}), F{beta}_Score={f_beta_score:.4f}\")\n    return f_beta_score, stats\n\n\n# %% [markdown]\n# # 7. Prediction Pipeline on Specific Tomograms\n#\n# Apply the (non-existent) models ONLY to the requested test set IDs: `tomo_003acc`, `tomo_00e047`, `tomo_00e463`.\n# This will load the data but predict NO_MOTOR_COORD because models are None.\n\n# %%\n# --- Prediction Function (Definition unchanged, behavior changes as models are None) ---\ndef predict_motor_location(tomo_data, voxel_spacing, clf_model, reg_model, patch_size_px=PATCH_SIZE_PX, n_samples=1000, prob_threshold=0.6):\n    \"\"\"Predicts motor location by sampling patches.\"\"\"\n    # This check is now crucial and will always be true in this script version\n    if clf_model is None:\n        # No need to print warning every time here, we know it's None\n        return [NO_MOTOR_COORD] * 3\n    # The rest of the function will not execute if clf_model is None\n\n    # --- Original function logic (will not run) ---\n    if tomo_data is None or voxel_spacing <= 0: return [NO_MOTOR_COORD] * 3\n    pred_start_time = time.time()\n    preprocessed_data = preprocess_tomogram(tomo_data)\n    if preprocessed_data is None: return [NO_MOTOR_COORD] * 3\n    tomo_shape = preprocessed_data.shape\n    candidate_features = []; candidate_centers_px = []\n    n_samples_to_use = n_samples if n_samples > 0 else 2500\n    for i in range(n_samples_to_use):\n        center_px = np.random.randint(0, tomo_shape, size=3)\n        patch, _ = get_local_patch(preprocessed_data, center_px, patch_size_px=patch_size_px)\n        if patch is not None:\n            features = extract_features_patch(patch)\n            if features is not None: candidate_features.append(features); candidate_centers_px.append(center_px)\n    if not candidate_features: return [NO_MOTOR_COORD] * 3\n    candidate_features_arr = np.array(candidate_features, dtype=np.float32)\n    try:\n        if candidate_features_arr.shape[0] == 0: return [NO_MOTOR_COORD] * 3\n        probs = clf_model.predict_proba(candidate_features_arr)[:, 1]\n    except Exception as e: print(f\"Error predict_proba: {e}\"); return [NO_MOTOR_COORD] * 3\n    if len(probs) == 0: return [NO_MOTOR_COORD] * 3\n    max_prob_idx = np.argmax(probs); max_prob = probs[max_prob_idx]; best_center_px = candidate_centers_px[max_prob_idx]\n    if max_prob >= prob_threshold:\n        best_features = candidate_features_arr[max_prob_idx:max_prob_idx+1]; predicted_offset_px = np.zeros(3)\n        if reg_model is not None:\n            try:\n                 if best_features.shape[0] > 0: predicted_offset_px = reg_model.predict(best_features)[0]\n            except Exception as e: print(f\"Warn reg_predict: {e}\")\n        final_predicted_coords_px = best_center_px + predicted_offset_px\n        final_predicted_coords_a = final_predicted_coords_px * voxel_spacing\n        max_coords_a = (np.array(tomo_shape) - 1) * voxel_spacing\n        final_predicted_coords_a = np.clip(final_predicted_coords_a, 0, max_coords_a); return final_predicted_coords_a.tolist()\n    else: return [NO_MOTOR_COORD] * 3\n    # --- End of original function logic ---\n\n# --- Generate Predictions ONLY for Specific Tomogram IDs ---\nprint(\"\\nGenerating DEFAULT Predictions for SPECIFIC Tomograms (Models were not trained)...\")\npredictions = []\ntest_predict_start_time = time.time()\n\nids_to_predict = ['tomo_003acc', 'tomo_00e047', 'tomo_00e463']\nprint(f\"Processing ONLY the following {len(ids_to_predict)} tomograms: {ids_to_predict}\")\n\nif not ids_to_predict:\n     print(\"No specific tomograms defined to process.\")\n# No need to check clf_model here, predict_motor_location handles it\nelse:\n    if voxel_spacing_map:\n        valid_spacings = [v for v in voxel_spacing_map.values() if v > 0]\n        avg_voxel_spacing = np.mean(valid_spacings) if valid_spacings else DEFAULT_VOXEL_SPACING_A\n    else: avg_voxel_spacing = DEFAULT_VOXEL_SPACING_A\n    print(f\"Using average voxel spacing (fallback): {avg_voxel_spacing:.4f} A/px\")\n\n    for tomo_idx, tomo_id in enumerate(ids_to_predict):\n        print(f\"Generating default prediction for {tomo_idx+1}/{len(ids_to_predict)}: {tomo_id}...\")\n        tomo_dir_path = test_dir_map.get(tomo_id)\n        source_dir = \"test\"\n        if not tomo_dir_path:\n            tomo_dir_path = train_dir_map.get(tomo_id)\n            source_dir = \"train\"\n        if not tomo_dir_path:\n            print(f\"  Skipping {tomo_id}: Directory not found.\")\n            pred_coords = [NO_MOTOR_COORD] * 3\n            predictions.append({'tomo_id': tomo_id,'Motor axis 0': pred_coords[0],'Motor axis 1': pred_coords[1],'Motor axis 2': pred_coords[2]})\n            continue\n\n        # Load data to get shape info, even though model won't use features\n        tomo_data = load_tomogram(tomo_id, {tomo_id: tomo_dir_path})\n        current_voxel_spacing = voxel_spacing_map.get(tomo_id, avg_voxel_spacing)\n\n        # Call prediction function - it will return default due to clf_model being None\n        pred_coords = predict_motor_location(\n            tomo_data, current_voxel_spacing, clf_model, reg_model, # Pass None models\n            patch_size_px=PATCH_SIZE_PX, n_samples=10, prob_threshold=0.99 # Params don't matter much here\n        )\n        print(f\"  Predicted Coords (Default) for {tomo_id}: {[f'{c:.2f}' for c in pred_coords]}\")\n\n        predictions.append({'tomo_id': tomo_id,'Motor axis 0': pred_coords[0],'Motor axis 1': pred_coords[1],'Motor axis 2': pred_coords[2]})\n        if tomo_data is not None: del tomo_data; gc.collect()\n\n    print(f\"\\nFinished default predictions. Total time: {time.time() - test_predict_start_time:.2f} seconds.\")\n\n# Create dataframe from the default predictions\nif predictions:\n    submission_df = pd.DataFrame(predictions)\n    print(\"\\nDefault Predictions for the specified IDs:\")\n    print(submission_df)\nelse:\n    print(\"\\nNo default predictions were generated (check specific IDs).\")\n    submission_df = pd.DataFrame(columns=['tomo_id', 'Motor axis 0', 'Motor axis 1', 'Motor axis 2'])\n\n\n# %% [markdown]\n# # 8. Create Submission File\n#\n# Generate the submission file template, ensuring all original test IDs are present.\n\n# %%\nprint(\"\\nCreating Submission File Template...\")\nif 'submission_df' in locals():\n    try:\n        sample_sub = pd.read_csv(SAMPLE_SUB_PATH)\n        expected_test_ids = test_ids_found # Use all IDs found in the test directory\n\n        if not expected_test_ids:\n             print(\"Warning: No test IDs were found in the test directory structure.\")\n             expected_test_ids = sample_sub['tomo_id'].tolist()\n             print(f\"Using {len(expected_test_ids)} IDs from sample submission as base.\")\n        else:\n             print(f\"Generating submission structure based on {len(expected_test_ids)} test IDs found in directory.\")\n\n        sub_df_final = pd.DataFrame({'tomo_id': expected_test_ids})\n\n        # Merge the default predictions we generated for the *specific* IDs\n        if not submission_df.empty:\n             sub_df_final = sub_df_final.merge(submission_df, on='tomo_id', how='left')\n        else:\n             # If even the default predictions weren't made, fill everything\n             for col in ['Motor axis 0', 'Motor axis 1', 'Motor axis 2']: sub_df_final[col] = NO_MOTOR_COORD\n\n        # Fill any remaining NaNs (for test IDs not in the specific prediction list) with default\n        sub_df_final.fillna(NO_MOTOR_COORD, inplace=True)\n\n        # Ensure column order and types match sample\n        final_columns = sample_sub.columns.tolist()\n        for col in final_columns:\n             if col not in sub_df_final.columns:\n                  print(f\"Warning: Column '{col}' from sample not in generated df. Adding default.\")\n                  sub_df_final[col] = NO_MOTOR_COORD\n        sub_df_final = sub_df_final[final_columns]\n        for col in ['Motor axis 0', 'Motor axis 1', 'Motor axis 2']:\n            sub_df_final[col] = pd.to_numeric(sub_df_final[col], errors='coerce').fillna(NO_MOTOR_COORD).astype(float)\n\n        sub_df_final.to_csv('submission.csv', index=False)\n        print(\"Submission file template created: submission.csv\")\n        print(\"\\nFinal Submission Head:\")\n        print(sub_df_final.head())\n        print(f\"Submission shape: {sub_df_final.shape}\")\n\n        if len(sub_df_final) != len(expected_test_ids):\n            print(f\"CRITICAL WARNING: Submission row count ({len(sub_df_final)}) \"\n                  f\"mismatches expected test IDs ({len(expected_test_ids)})!\")\n        else:\n             print(\"Submission row count matches expected test IDs.\")\n\n    except FileNotFoundError:\n        print(f\"Error: Sample submission file not found at {SAMPLE_SUB_PATH}. Cannot verify format.\")\n        # Try to save anyway if df exists\n        if 'sub_df_final' in locals():\n             sub_df_final.to_csv('submission.csv', index=False)\n             print(\"Submission file template created (format not verified): submission.csv\")\n    except Exception as e:\n        print(f\"Error creating final submission file: {e}\")\nelse:\n    print(\"No prediction DataFrame ('submission_df') was generated. Cannot create submission template.\")\n\n\n# %% [markdown]\n# # 9. Discussion and Future Work\n#\n# *   **Data Loading & Handling:** Code structure maintained for loading data, labels, and finding directories.\n# *   **Feature Engineering:** Utility functions (`preprocess_tomogram`, `extract_features_patch`) remain defined but were not used for training.\n# *   **Modeling:** **Training data generation and model training sections were entirely skipped.** `clf_model` and `reg_model` were explicitly set to `None`.\n# *   **Targeted Prediction:** The prediction loop for specific IDs (`tomo_003acc`, `tomo_00e047`, `tomo_00e463`) was executed. However, since the models were `None`, the `predict_motor_location` function returned the default `NO_MOTOR_COORD` for these IDs.\n# *   **Submission Generation:** A `submission.csv` file was successfully created containing *all* test IDs found in the `test` directory. The specifically targeted IDs have `NO_MOTOR_COORD` because no model was trained, and all other test IDs were also filled with `NO_MOTOR_COORD`. This serves as a valid submission template.\n# *   **Limitations:** This script produces a baseline submission with no actual detection. Its primary purpose is to verify the data loading, directory handling, and submission file formatting logic.\n# *   **Future Improvements:** To get actual predictions, the skipping of Sections 5 and 5.1 must be removed, and sufficient training data (likely more than the 20 used in the previous version) needs to be processed to train meaningful models. Then, the prediction pipeline in Section 7 will use the trained models. All other future improvements mentioned previously (CNNs, FCNs, augmentation, etc.) still apply for building a competitive model.\n\n# %%\nprint(\"Script finished (Training and Model Fitting Skipped). Submission template generated.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T22:48:35.522190Z","iopub.execute_input":"2025-04-16T22:48:35.522433Z","iopub.status.idle":"2025-04-16T22:49:16.813623Z","shell.execute_reply.started":"2025-04-16T22:48:35.522418Z","shell.execute_reply":"2025-04-16T22:49:16.813051Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -*- coding: utf-8 -*-\n\"\"\"\nKaggle Notebook: BYU Flagellar Motor Localization (Revised for Directory Structure & Multi-Motor)\n\nGoal: Detect the presence and (x, y, z) coordinates of bacterial\n      flagellar motors in 3D cryo-ET tomograms stored as slice directories.\n\nInspired Concepts: Analogies from Non-Linear Electrodynamics / String Theory\n    - Tomogram Intensity: Scalar field (like potential)\n    - Gradient (∇I): Vector field (analogous to Electric Field E)\n    - Structure Tensor (∇I ⊗ ∇I averaged): Symmetric Tensor (analogous to Stress Tensor T_µν or Metric Perturbation S_µν)\n    - Hessian (∇²I): Tensor capturing curvature (related to field changes)\n    - Invariants: Rotationally invariant features derived from gradient/structure tensor/Hessian (analogous to Lorentz invariants x, y)\n    - Anisotropic Diffusion/Filtering: Guided propagation of information (analogous to wave propagation influenced by background fields/metrics)\n    - Multi-Metric Idea: Using different analysis scales/features (like g_µν vs G_µν) - e.g., global search vs. local refinement.\n\nMODIFICATION: This version SKIPS training data generation and model training entirely.\n            It focuses prediction ONLY on specified IDs ['tomo_003acc', 'tomo_00e047', 'tomo_00e463']\n            using NO trained model (will predict NO_MOTOR_COORD).\n            The final output is a submission file template for ALL test IDs found.\n\"\"\"\n\n# %% [markdown]\n# # 1. Setup and Imports\n#\n# Load necessary libraries and define constants.\n\n# %%\n# Attempt installation only if needed (e.g., in Kaggle environment)\ntry:\n    import cv2\n    import mrcfile # Though likely unused for loading now\n    print(\"Libraries opencv-python, mrcfile seem available.\")\nexcept ImportError:\n    print(\"Installing required libraries: opencv-python, mrcfile\")\n    # Use %pip instead of !pip in notebooks for better integration\n    %pip install opencv-python mrcfile --quiet\n    import cv2\n    import mrcfile\n    print(\"Libraries installed.\")\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport cv2 # Import OpenCV\n# import mrcfile # Keep if other parts might use it, but not for loading slices\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom scipy import ndimage as ndi\nfrom skimage.feature import structure_tensor, structure_tensor_eigenvalues, hessian_matrix, hessian_matrix_eigvals\nfrom skimage.filters import gaussian, median\nfrom skimage.measure import regionprops, label\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\nfrom sklearn.metrics import fbeta_score, mean_squared_error, make_scorer\nimport gc # Garbage collection\nimport time # For timing\n\n# Optional: Deep Learning Libraries (if used)\n# import torch\n# import torch.nn as nn\n# from torch.utils.data import Dataset, DataLoader\n# import monai # Medical Imaging AI library\n\n# Constants\nCOMPETITION_DIR = '/kaggle/input/byu-locating-bacterial-flagellar-motors-2025'\nTRAIN_DIR = os.path.join(COMPETITION_DIR, 'train')\nTEST_DIR = os.path.join(COMPETITION_DIR, 'test')\nTRAIN_LABELS_PATH = os.path.join(COMPETITION_DIR, 'train_labels.csv')\nSAMPLE_SUB_PATH = os.path.join(COMPETITION_DIR, 'sample_submission.csv')\n\n# Default Voxel spacing (used for test set or if lookup fails)\nDEFAULT_VOXEL_SPACING_A = 10.0\nprint(f\"Using default voxel spacing: {DEFAULT_VOXEL_SPACING_A} Angstroms/pixel (will try to read from labels for training)\")\n\n# Evaluation Threshold (Angstroms)\nDISTANCE_THRESHOLD_A = 1000.0\n# DISTANCE_THRESHOLD_PX will be calculated dynamically based on actual voxel spacing\n\n# F-beta Score Beta value\nF_BETA = 2.0\n\n# For submission: Coordinate value indicating no motor found\nNO_MOTOR_COORD = -1.0 # Ensure float\n\n# Seed for reproducibility\nSEED = 42\nnp.random.seed(SEED)\n\n# %% [markdown]\n# # 2. Load Data and Metadata\n#\n# Load the training labels and find the data directories. Check for consistency.\n\n# %%\nstart_time = time.time()\nprint(\"Loading labels...\")\ntry:\n    train_labels_df = pd.read_csv(TRAIN_LABELS_PATH)\n    print(f\"Training labels shape: {train_labels_df.shape}\")\n    print(train_labels_df.head())\nexcept FileNotFoundError:\n    print(f\"ERROR: Training labels file not found at {TRAIN_LABELS_PATH}\")\n    train_labels_df = pd.DataFrame() # Empty dataframe\n\n# --- Debug: List directories to confirm paths ---\nprint(\"\\n--- Directory Listing ---\")\nprint(f\"Listing {COMPETITION_DIR}:\")\ntry:\n    print(os.listdir(COMPETITION_DIR))\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {COMPETITION_DIR}\")\n\nprint(f\"\\nListing {TRAIN_DIR}:\")\ntry:\n    train_contents = os.listdir(TRAIN_DIR)\n    print(f\"  Found {len(train_contents)} items in train dir. First 10:\")\n    print(train_contents[:10])\n    if len(train_contents) > 10: print(\"  ...\")\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {TRAIN_DIR}\")\n    train_contents = []\n\nprint(f\"\\nListing {TEST_DIR}:\")\ntry:\n    test_contents = os.listdir(TEST_DIR)\n    print(f\"  Found {len(test_contents)} items in test dir. First 10:\")\n    print(test_contents[:10])\n    if len(test_contents) > 10: print(\"  ...\")\nexcept FileNotFoundError:\n    print(f\"  Error: Directory not found: {TEST_DIR}\")\n    test_contents = []\nprint(\"--- End Directory Listing ---\\n\")\n\n\n# --- Find Tomogram Directories ---\nprint(\"Finding tomogram directories...\")\ntrain_dirs = sorted(glob.glob(os.path.join(TRAIN_DIR, 'tomo_*')))\ntest_dirs = sorted(glob.glob(os.path.join(TEST_DIR, 'tomo_*')))\n\n# Extract UNIQUE tomo_ids present in the directories found\ntrain_ids_found = sorted([os.path.basename(d) for d in train_dirs if os.path.isdir(d)]) # Ensure it's a directory\ntest_ids_found = sorted([os.path.basename(d) for d in test_dirs if os.path.isdir(d)]) # Ensure it's a directory\n\nprint(f\"Found {len(train_ids_found)} potential training tomogram directories.\")\nprint(f\"First 5 found: {train_ids_found[:5]}\")\nprint(f\"Found {len(test_ids_found)} potential testing tomogram directories.\")\nprint(f\"First 5 found: {test_ids_found[:5]}\")\n\n# --- Filter Labels based on Found Directories ---\nif not train_labels_df.empty:\n    # Check if 'tomo_id' column exists\n    if 'tomo_id' not in train_labels_df.columns:\n        print(\"ERROR: 'tomo_id' column not found in labels CSV. Cannot proceed.\")\n        train_labels_df_filtered = pd.DataFrame()\n        unique_train_ids_with_data = []\n        voxel_spacing_map = {}\n    else:\n        train_labels_df_filtered = train_labels_df[train_labels_df['tomo_id'].isin(train_ids_found)].copy()\n        print(f\"\\nOriginal label rows: {len(train_labels_df)}\")\n        print(f\"Label rows after filtering by found train directories: {len(train_labels_df_filtered)}\")\n\n        # Get unique tomogram IDs that have both labels AND a found directory\n        unique_train_ids_with_data = sorted(train_labels_df_filtered['tomo_id'].unique())\n        print(f\"Number of unique train tomograms with labels AND data: {len(unique_train_ids_with_data)}\")\n\n        # Create map for Voxel Spacing from labels (use first entry per tomo_id)\n        if 'Voxel spacing' in train_labels_df_filtered.columns and not train_labels_df_filtered.empty:\n            voxel_spacing_map = train_labels_df_filtered.drop_duplicates(subset='tomo_id').set_index('tomo_id')['Voxel spacing'].to_dict()\n            print(f\"Example voxel spacing: {list(voxel_spacing_map.items())[:5]}\")\n        else:\n             print(\"WARNING: 'Voxel spacing' column not found or no valid labels. Cannot create voxel spacing map.\")\n             voxel_spacing_map = {}\n\nelse:\n    print(\"WARNING: Training labels dataframe is empty. Cannot filter or create maps.\")\n    train_labels_df_filtered = pd.DataFrame()\n    unique_train_ids_with_data = []\n    voxel_spacing_map = {}\n\n# Map tomo_id to DIRECTORY path (even if labels are missing)\ntrain_dir_map = {os.path.basename(d): d for d in train_dirs if os.path.isdir(d)}\ntest_dir_map = {os.path.basename(d): d for d in test_dirs if os.path.isdir(d)}\n\n\nprint(f\"Data loading setup took {time.time() - start_time:.2f} seconds.\")\n\nif not unique_train_ids_with_data and not train_labels_df.empty:\n    print(\"\\n\\nCRITICAL WARNING: No intersection between labels and found training directories. Training cannot proceed.\")\nelif not train_ids_found and not test_ids_found:\n     print(\"\\n\\nCRITICAL WARNING: No training or testing directories found. Check data paths.\")\n\n# %% [markdown]\n# # 3. Exploratory Data Analysis (EDA) & Concept Visualization\n#\n# Understand the data distribution, visualize tomograms and labels, and see how our physics-inspired concepts manifest.\n# (Keeping these utility functions defined, even if visualization is partially skipped)\n\n# %%\n# --- Tomogram Loading Function ---\ndef load_tomogram(tomo_id, dir_map):\n    \"\"\"Loads a tomogram by reading and stacking slices from a directory.\"\"\"\n    if isinstance(dir_map, str):\n        directory_path = dir_map\n        if not os.path.isdir(directory_path): return None\n    else:\n        directory_path = dir_map.get(tomo_id)\n        if not directory_path or not os.path.isdir(directory_path): return None\n\n    slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.jpg')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.png')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.tif')))\n    if not slice_files: return None\n\n    def get_slice_number(filepath):\n        try:\n            filename = os.path.basename(filepath)\n            num_str = filename.split('slice_')[1].split('.')[0]\n            return int(num_str)\n        except: return -1\n\n    valid_slice_files = [(f, get_slice_number(f)) for f in slice_files]\n    valid_slice_files = [(f, num) for f, num in valid_slice_files if num != -1]\n    if not valid_slice_files: return None\n    valid_slice_files.sort(key=lambda item: item[1])\n    slice_files_sorted = [f for f, num in valid_slice_files]\n\n    try:\n        first_slice = cv2.imread(slice_files_sorted[0], cv2.IMREAD_GRAYSCALE)\n        if first_slice is None: raise IOError(f\"cv2.imread failed for {slice_files_sorted[0]}\")\n        height, width = first_slice.shape\n        num_slices = len(slice_files_sorted)\n        dtype = first_slice.dtype\n    except Exception as e:\n        print(f\"Error reading first slice {slice_files_sorted[0]}: {e}\"); return None\n\n    tomogram_data = np.zeros((num_slices, height, width), dtype=dtype)\n    tomogram_data[0, :, :] = first_slice\n\n    for i in range(1, num_slices):\n        try:\n            slice_img = cv2.imread(slice_files_sorted[i], cv2.IMREAD_GRAYSCALE)\n            if slice_img is None: raise IOError(f\"cv2.imread failed for {slice_files_sorted[i]}\")\n            if slice_img.shape != (height, width): return None\n            tomogram_data[i, :, :] = slice_img\n        except Exception as e:\n            print(f\"Error reading slice {slice_files_sorted[i]}: {e}\"); return None\n    return tomogram_data\n\n# --- Plotting Function ---\ndef plot_slices_with_label(tomo_data, voxel_spacing, label_coords_A=None, title=\"Tomogram Slices\"):\n    \"\"\"Plots central slices along each axis, optionally marking the label(s).\"\"\"\n    if tomo_data is None: print(\"No data to plot.\"); return\n    if voxel_spacing <= 0: print(\"Invalid voxel spacing.\"); return\n    shape = tomo_data.shape\n    if not (len(shape) == 3 and all(s > 0 for s in shape)): print(f\"Invalid shape {shape}.\"); return\n\n    center_z, center_y, center_x = shape[0] // 2, shape[1] // 2, shape[2] // 2\n    fig, axes = plt.subplots(1, 3, figsize=(15, 5)); fig.suptitle(title, fontsize=16)\n    try:\n        axes[0].imshow(tomo_data[center_z, :, :], cmap='gray'); axes[0].set_title(f'Z Slice (Z={center_z})')\n        axes[1].imshow(tomo_data[:, center_y, :], cmap='gray', aspect=shape[0]/shape[2] if shape[2]>0 else 1); axes[1].set_title(f'Y Slice (Y={center_y})')\n        axes[2].imshow(tomo_data[:, :, center_x], cmap='gray', aspect=shape[0]/shape[1] if shape[1]>0 else 1); axes[2].set_title(f'X Slice (X={center_x})')\n    except IndexError as e: print(f\"Error plotting slices: {e}\"); plt.close(fig); return\n\n    if label_coords_A is not None and len(label_coords_A) > 0:\n        is_multiple = isinstance(label_coords_A[0], (list, np.ndarray))\n        if not is_multiple: label_coords_A = [label_coords_A]\n        plotted_legend = False\n        for i, motor_A in enumerate(label_coords_A):\n            if len(motor_A) < 3 or motor_A[0] == NO_MOTOR_COORD: continue\n            lz_A, ly_A, lx_A = motor_A\n            lz_px, ly_px, lx_px = int(lz_A / voxel_spacing), int(ly_A / voxel_spacing), int(lx_A / voxel_spacing)\n            if not (0 <= lz_px < shape[0] and 0 <= ly_px < shape[1] and 0 <= lx_px < shape[2]): continue\n            label_text = f'Motor {i+1}' if len(label_coords_A) > 1 else 'Motor'\n            z_thresh=max(5, shape[0]*0.05); y_thresh=max(5, shape[1]*0.05); x_thresh=max(5, shape[2]*0.05)\n            try:\n                lbl = None\n                if not plotted_legend: lbl = label_text; plotted_legend=True\n                if abs(lz_px - center_z) < z_thresh: axes[0].plot(lx_px, ly_px, 'ro', ms=8, label=lbl)\n                if abs(ly_px - center_y) < y_thresh: axes[1].plot(lx_px, lz_px, 'ro', ms=8)\n                if abs(lx_px - center_x) < x_thresh: axes[2].plot(ly_px, lz_px, 'ro', ms=8)\n                if lbl: lbl=None\n            except Exception as plot_e: print(f\"Error plotting marker: {plot_e}\")\n    if plotted_legend: axes[0].legend()\n    plt.tight_layout(rect=[0, 0.03, 1, 0.95]); plt.show()\n\n# --- Visualize a sample tomogram ---\n# (Keep this visualization, it's quick and useful)\nif unique_train_ids_with_data:\n    sample_tomo_id = unique_train_ids_with_data[0]\n    print(f\"\\nVisualizing sample tomogram: {sample_tomo_id}\")\n    sample_data = load_tomogram(sample_tomo_id, train_dir_map)\n    sample_labels = train_labels_df_filtered[train_labels_df_filtered['tomo_id'] == sample_tomo_id]\n    if sample_data is not None and not sample_labels.empty:\n        motor_locations_A = sample_labels[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].values.tolist()\n        voxel_spacing = sample_labels.iloc[0]['Voxel spacing']\n        print(f\"Sample Tomogram Shape: {sample_data.shape}, Voxel Spacing: {voxel_spacing} A/px\")\n        print(f\"Motor Location(s) (Angstroms): {motor_locations_A}\")\n        plot_slices_with_label(sample_data, voxel_spacing, motor_locations_A, title=f\"{sample_tomo_id}\")\n        plt.figure(figsize=(10, 4)); plt.hist(sample_data.flatten(), bins=100, color='blue', alpha=0.7); plt.title(f'Intensity Distribution {sample_tomo_id}'); plt.xlabel('Intensity'); plt.ylabel('Frequency'); plt.grid(True, alpha=0.3); plt.show()\n        del sample_data, sample_labels; gc.collect()\n    else: print(f\"Failed to load/find labels for sample {sample_tomo_id}\")\nelse: print(\"\\nSkipping Visualization: No valid training tomograms found.\")\n\n# %% [markdown]\n# ### 3.1 Visualizing Physics-Inspired Concepts\n# (Skip this visualization as it relies on patch extraction and calculations that mirror training data gen)\n\n# %%\n# --- Patch Extraction Function (Keep definition for prediction step) ---\ndef get_local_patch(tomo_data, center_px, patch_size_px=64):\n    \"\"\"Extracts a 3D patch centered at center_px. Patch size is in pixels.\"\"\"\n    if tomo_data is None or center_px is None: return None, None\n    center_px = np.round(center_px).astype(int); z, y, x = center_px\n    shape = tomo_data.shape\n    if not (0 <= z < shape[0] and 0 <= y < shape[1] and 0 <= x < shape[2]): return None, None\n    half_size = patch_size_px // 2\n    z_start, z_end = max(0, z - half_size), min(shape[0], z + half_size)\n    y_start, y_end = max(0, y - half_size), min(shape[1], y + half_size)\n    x_start, x_end = max(0, x - half_size), min(shape[2], x + half_size)\n    dest_z_start = half_size - (z - z_start); dest_z_end = dest_z_start + (z_end - z_start)\n    dest_y_start = half_size - (y - y_start); dest_y_end = dest_y_start + (y_end - y_start)\n    dest_x_start = half_size - (x - x_start); dest_x_end = dest_x_start + (x_end - x_start)\n    patch = np.zeros((patch_size_px, patch_size_px, patch_size_px), dtype=tomo_data.dtype)\n    try:\n        src_slice = (slice(z_start, z_end), slice(y_start, y_end), slice(x_start, x_end))\n        dest_slice = (slice(dest_z_start, dest_z_end), slice(dest_y_start, dest_y_end), slice(dest_x_start, dest_x_end))\n        patch[dest_slice] = tomo_data[src_slice]\n        patch_center_coords_in_patch = np.round([z - z_start + dest_z_start, y - y_start + dest_y_start, x - x_start + dest_x_start]).astype(int)\n        if not (0 <= patch_center_coords_in_patch[0] < patch_size_px and 0 <= patch_center_coords_in_patch[1] < patch_size_px and 0 <= patch_center_coords_in_patch[2] < patch_size_px):\n             patch_center_coords_in_patch = [patch_size_px // 2] * 3\n        return patch, patch_center_coords_in_patch\n    except (ValueError, IndexError) as e: return None, None\n\n# --- Visualize physics concepts ---\nprint(\"\\nSkipping Physics Concept Visualization (requires patch processing).\")\n# (The code for visualization is removed here)\n\n# %% [markdown]\n# # 4. Data Preprocessing & Feature Engineering\n# (Keep function definitions as they are needed for prediction, even if models aren't trained)\n\n# %%\n# --- Simple Preprocessing ---\ndef preprocess_tomogram(data):\n    \"\"\"Normalizes and smooths the tomogram data.\"\"\"\n    if data is None: return None\n    if not np.issubdtype(data.dtype, np.floating): data = data.astype(np.float32)\n    mean, std = np.mean(data), np.std(data)\n    if std > 1e-6: normalized_data = (data - mean) / std\n    else: normalized_data = data - mean\n    smoothed_data = gaussian(normalized_data, sigma=1.5, mode='reflect', preserve_range=True, truncate=4.0)\n    return smoothed_data\n\n# --- Feature Extraction (Keep definition) ---\ndef extract_features_patch(patch):\n    \"\"\"Extracts a feature vector from a 3D patch.\"\"\"\n    if patch is None or patch.size == 0: return None\n    if not np.issubdtype(patch.dtype, np.floating): patch = patch.astype(np.float32)\n    features = []\n    EXPECTED_N_FEATURES = 15\n    try:\n        features.extend([np.mean(patch), np.std(patch), np.median(patch), np.min(patch), np.max(patch)])\n        patch_smooth_grad = gaussian(patch, sigma=1.0, mode='reflect', preserve_range=True, truncate=4.0)\n        patch_smooth_st = gaussian(patch, sigma=1.5, mode='reflect', preserve_range=True, truncate=4.0)\n        patch_smooth_hess = gaussian(patch, sigma=2.0, mode='reflect', preserve_range=True, truncate=4.0)\n        try:\n            grad_z, grad_y, grad_x = np.gradient(patch_smooth_grad)\n            grad_mag = np.sqrt(grad_z**2 + grad_y**2 + grad_x**2)\n            features.extend([np.mean(grad_mag), np.std(grad_mag), np.max(grad_mag)])\n        except Exception: features.extend([0.0] * 3)\n        try:\n            shape = patch_smooth_st.shape; slice_start = max(0, shape[0]//5); slice_end = max(slice_start+1, 4*shape[0]//5)\n            center_slice = slice(slice_start, slice_end)\n            if not (0<=center_slice.start<shape[0] and 0<center_slice.stop<=shape[0] and center_slice.start<center_slice.stop): center_patch_smooth = patch_smooth_st\n            else: center_patch_smooth = patch_smooth_st[center_slice, center_slice, center_slice]\n            if center_patch_smooth.size > 0:\n                S_elems_c = structure_tensor(center_patch_smooth, sigma=1.5, mode='reflect')\n                eigvals_S_c = structure_tensor_eigenvalues(S_elems_c); eigvals_S_c = np.sort(eigvals_S_c, axis=0)\n                l3, l2, l1 = eigvals_S_c[0], eigvals_S_c[1], eigvals_S_c[2]\n                den = l1 + l2 + l3 + 1e-9; coh = np.where(den > 1e-8, (l1 - l3) / den, 0.0)\n                features.extend([np.mean(l1), np.mean(l3), np.mean(coh), np.max(coh)])\n            else: features.extend([0.0] * 4)\n        except Exception: features.extend([0.0] * 4)\n        try:\n            shape_h = patch_smooth_hess.shape; slice_start_h = max(0, shape_h[0]//5); slice_end_h = max(slice_start_h+1, 4*shape_h[0]//5)\n            center_slice_h = slice(slice_start_h, slice_end_h)\n            if not (0<=center_slice_h.start<shape_h[0] and 0<center_slice_h.stop<=shape_h[0] and center_slice_h.start<center_slice_h.stop): center_patch_smooth_h = patch_smooth_hess\n            else: center_patch_smooth_h = patch_smooth_hess[center_slice_h, center_slice_h, center_slice_h]\n            if center_patch_smooth_h.size > 0:\n                h_matrix = hessian_matrix(center_patch_smooth_h, sigma=2.0, mode='reflect', use_gaussian_derivatives=False)\n                eigvals_H_c = hessian_matrix_eigvals(h_matrix); eigvals_H_c = np.sort(eigvals_H_c, axis=0)\n                h1, h2, h3 = eigvals_H_c[0], eigvals_H_c[1], eigvals_H_c[2]\n                features.extend([np.mean(h1), np.mean(h3), np.std(h1)])\n            else: features.extend([0.0] * 3)\n        except Exception: features.extend([0.0] * 3)\n    except Exception as outer_e: print(f\"Outer FE error: {outer_e}\"); return None\n    if len(features) != EXPECTED_N_FEATURES:\n        if len(features) < EXPECTED_N_FEATURES: features.extend([0.0] * (EXPECTED_N_FEATURES - len(features)))\n        else: features = features[:EXPECTED_N_FEATURES]\n        if len(features) != EXPECTED_N_FEATURES: return None\n    features_arr = np.array(features, dtype=np.float32)\n    if np.any(np.isnan(features_arr)) or np.any(np.isinf(features_arr)):\n        features_arr = np.nan_to_num(features_arr, nan=0.0, posinf=0.0, neginf=0.0)\n    if len(features_arr) != EXPECTED_N_FEATURES: return None\n    return features_arr\n\n\n# %% [markdown]\n# # 5. Model Development (Example: Patch Classification + Regression)\n#\n# <<< SKIPPING Training Data Generation >>>\n\n# %%\n# --- Generate Training Data ---\nPATCH_SIZE_PX = 64 # Keep constants defined\nEXPECTED_N_FEATURES = 15\n\nprint(\"\\nSKIPPING Training Data Generation as requested.\")\n\n# Define variables that would have been created as None or empty\nfeatures_list = []; labels_clf_list = []; labels_reg_list = []\nfailed_tomos = []\nclf_model = None # Crucial: Ensure models are None\nreg_model = None # Crucial: Ensure models are None\nX = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_clf = np.array([], dtype=np.int32)\ny_reg_all = np.array([]).reshape(0, 3)\nX_reg = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_reg = np.array([]).reshape(0, 3)\nX_train_clf = X_val_clf = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_train_clf = y_val_clf = np.array([], dtype=np.int32)\nX_train_reg = X_val_reg = np.array([]).reshape(0, EXPECTED_N_FEATURES)\ny_train_reg = y_val_reg = np.array([]).reshape(0, 3)\n\ngc.collect() # Clean up memory in case any large objects were created before skip\n\n\n# %% [markdown]\n# ### 5.1 Train Models\n#\n# <<< SKIPPING Model Training >>>\n\n# %%\n# --- Split data and Train ---\nprint(\"\\nSKIPPING Model Training as training data generation was skipped.\")\n\n# Models remain None as set in the previous step\nprint(f\"\\nModels after training attempts:\")\nprint(f\"  Classifier model: {'Available' if clf_model is not None else 'Not Available'}\")\nprint(f\"  Regressor model: {'Available' if reg_model is not None else 'Not Available'}\")\n\n\n# %% [markdown]\n# # 6. Evaluation Metric Implementation\n# (Keep definition, although it won't be used effectively without trained models/predictions)\n\n# %%\ndef calculate_fbeta_and_distance(y_true_df, y_pred_df, voxel_spacing_map, default_voxel_spacing, beta=F_BETA, threshold_a=DISTANCE_THRESHOLD_A):\n    \"\"\"\n    Calculates the competition metric (F-beta score based on distance).\n    Handles multiple ground truth motors per tomogram.\n    \"\"\"\n    tp = 0; fp = 0; fn = 0\n    no_motor = float(NO_MOTOR_COORD)\n\n    if y_pred_df.empty: pred_map = {}\n    else:\n        y_pred_df_unique = y_pred_df.drop_duplicates(subset='tomo_id', keep='first')\n        pred_map = y_pred_df_unique.set_index('tomo_id')[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].apply(lambda row: row.tolist(), axis=1).to_dict()\n\n    if y_true_df.empty: true_grouped = {}; total_true_motors = 0\n    else:\n        y_true_motors_only = y_true_df[y_true_df['Motor axis 0'] != no_motor].copy()\n        true_grouped = y_true_motors_only.groupby('tomo_id')\n        total_true_motors = len(y_true_motors_only)\n\n    tp_final = 0; fp_final = 0\n    matched_true_motor_indices = {}\n\n    for tomo_id, pred_coords_a in pred_map.items():\n        if pred_coords_a[0] == no_motor: continue\n        true_motors_a = []; true_indices = []; is_true_motor_present = False\n        if tomo_id in true_grouped.groups:\n             true_df_tomo = true_grouped.get_group(tomo_id)\n             true_motors_a = true_df_tomo[['Motor axis 0', 'Motor axis 1', 'Motor axis 2']].values\n             true_indices = true_df_tomo.index.tolist(); is_true_motor_present = True\n        if not is_true_motor_present: fp_final += 1; continue\n        pred_vec = np.array(pred_coords_a)\n        found_match = False; best_match_true_idx = -1; min_dist = float('inf')\n        for i, true_motor_a in enumerate(true_motors_a):\n            dist = np.linalg.norm(true_motor_a - pred_vec)\n            original_idx = true_indices[i]\n            if dist <= threshold_a:\n                 already_matched = tomo_id in matched_true_motor_indices and original_idx in matched_true_motor_indices[tomo_id]\n                 if not already_matched and dist < min_dist:\n                      min_dist = dist; best_match_true_idx = original_idx; found_match = True\n        if found_match:\n            tp_final += 1\n            if tomo_id not in matched_true_motor_indices: matched_true_motor_indices[tomo_id] = set()\n            matched_true_motor_indices[tomo_id].add(best_match_true_idx)\n        else: fp_final += 1\n    fn_final = total_true_motors - tp_final\n    f_beta_numerator = (1 + beta**2) * tp_final\n    f_beta_denominator = (1 + beta**2) * tp_final + (beta**2 * fn_final) + fp_final\n    if f_beta_denominator == 0: f_beta_score = 1.0 if total_true_motors == 0 else 0.0\n    else: f_beta_score = f_beta_numerator / f_beta_denominator\n    stats = {'TP': tp_final, 'FP': fp_final, 'FN': fn_final, 'Total True': total_true_motors}\n    # Print score even if it's based on default predictions\n    print(f\"Evaluation Stats (based on default predictions): TP={tp_final}, FP={fp_final}, FN={fn_final} (Total True={total_true_motors}), F{beta}_Score={f_beta_score:.4f}\")\n    return f_beta_score, stats\n\n\n# %% [markdown]\n# # 7. Prediction Pipeline on Specific Tomograms\n#\n# Apply the (non-existent) models ONLY to the requested test set IDs: `tomo_003acc`, `tomo_00e047`, `tomo_00e463`.\n# This will load the data but predict NO_MOTOR_COORD because models are None.\n\n# %%\n# --- Prediction Function (Definition unchanged, behavior changes as models are None) ---\ndef predict_motor_location(tomo_data, voxel_spacing, clf_model, reg_model, patch_size_px=PATCH_SIZE_PX, n_samples=1000, prob_threshold=0.6):\n    \"\"\"Predicts motor location by sampling patches.\"\"\"\n    # This check is now crucial and will always be true in this script version\n    if clf_model is None:\n        # No need to print warning every time here, we know it's None\n        return [NO_MOTOR_COORD] * 3\n    # The rest of the function will not execute if clf_model is None\n\n    # --- Original function logic (will not run) ---\n    if tomo_data is None or voxel_spacing <= 0: return [NO_MOTOR_COORD] * 3\n    pred_start_time = time.time()\n    preprocessed_data = preprocess_tomogram(tomo_data)\n    if preprocessed_data is None: return [NO_MOTOR_COORD] * 3\n    tomo_shape = preprocessed_data.shape\n    candidate_features = []; candidate_centers_px = []\n    n_samples_to_use = n_samples if n_samples > 0 else 2500\n    for i in range(n_samples_to_use):\n        center_px = np.random.randint(0, tomo_shape, size=3)\n        patch, _ = get_local_patch(preprocessed_data, center_px, patch_size_px=patch_size_px)\n        if patch is not None:\n            features = extract_features_patch(patch)\n            if features is not None: candidate_features.append(features); candidate_centers_px.append(center_px)\n    if not candidate_features: return [NO_MOTOR_COORD] * 3\n    candidate_features_arr = np.array(candidate_features, dtype=np.float32)\n    try:\n        if candidate_features_arr.shape[0] == 0: return [NO_MOTOR_COORD] * 3\n        probs = clf_model.predict_proba(candidate_features_arr)[:, 1]\n    except Exception as e: print(f\"Error predict_proba: {e}\"); return [NO_MOTOR_COORD] * 3\n    if len(probs) == 0: return [NO_MOTOR_COORD] * 3\n    max_prob_idx = np.argmax(probs); max_prob = probs[max_prob_idx]; best_center_px = candidate_centers_px[max_prob_idx]\n    if max_prob >= prob_threshold:\n        best_features = candidate_features_arr[max_prob_idx:max_prob_idx+1]; predicted_offset_px = np.zeros(3)\n        if reg_model is not None:\n            try:\n                 if best_features.shape[0] > 0: predicted_offset_px = reg_model.predict(best_features)[0]\n            except Exception as e: print(f\"Warn reg_predict: {e}\")\n        final_predicted_coords_px = best_center_px + predicted_offset_px\n        final_predicted_coords_a = final_predicted_coords_px * voxel_spacing\n        max_coords_a = (np.array(tomo_shape) - 1) * voxel_spacing\n        final_predicted_coords_a = np.clip(final_predicted_coords_a, 0, max_coords_a); return final_predicted_coords_a.tolist()\n    else: return [NO_MOTOR_COORD] * 3\n    # --- End of original function logic ---\n\n# --- Generate Predictions ONLY for Specific Tomogram IDs ---\nprint(\"\\nGenerating DEFAULT Predictions for SPECIFIC Tomograms (Models were not trained)...\")\npredictions = []\ntest_predict_start_time = time.time()\n\nids_to_predict = ['tomo_003acc', 'tomo_00e047', 'tomo_00e463']\nprint(f\"Processing ONLY the following {len(ids_to_predict)} tomograms: {ids_to_predict}\")\n\nif not ids_to_predict:\n     print(\"No specific tomograms defined to process.\")\n# No need to check clf_model here, predict_motor_location handles it\nelse:\n    if voxel_spacing_map:\n        valid_spacings = [v for v in voxel_spacing_map.values() if v > 0]\n        avg_voxel_spacing = np.mean(valid_spacings) if valid_spacings else DEFAULT_VOXEL_SPACING_A\n    else: avg_voxel_spacing = DEFAULT_VOXEL_SPACING_A\n    print(f\"Using average voxel spacing (fallback): {avg_voxel_spacing:.4f} A/px\")\n\n    for tomo_idx, tomo_id in enumerate(ids_to_predict):\n        print(f\"Generating default prediction for {tomo_idx+1}/{len(ids_to_predict)}: {tomo_id}...\")\n        tomo_dir_path = test_dir_map.get(tomo_id)\n        source_dir = \"test\"\n        if not tomo_dir_path:\n            tomo_dir_path = train_dir_map.get(tomo_id)\n            source_dir = \"train\"\n        if not tomo_dir_path:\n            print(f\"  Skipping {tomo_id}: Directory not found.\")\n            pred_coords = [NO_MOTOR_COORD] * 3\n            predictions.append({'tomo_id': tomo_id,'Motor axis 0': pred_coords[0],'Motor axis 1': pred_coords[1],'Motor axis 2': pred_coords[2]})\n            continue\n\n        # Load data to get shape info, even though model won't use features\n        tomo_data = load_tomogram(tomo_id, {tomo_id: tomo_dir_path})\n        current_voxel_spacing = voxel_spacing_map.get(tomo_id, avg_voxel_spacing)\n\n        # Call prediction function - it will return default due to clf_model being None\n        pred_coords = predict_motor_location(\n            tomo_data, current_voxel_spacing, clf_model, reg_model, # Pass None models\n            patch_size_px=PATCH_SIZE_PX, n_samples=10, prob_threshold=0.99 # Params don't matter much here\n        )\n        print(f\"  Predicted Coords (Default) for {tomo_id}: {[f'{c:.2f}' for c in pred_coords]}\")\n\n        predictions.append({'tomo_id': tomo_id,'Motor axis 0': pred_coords[0],'Motor axis 1': pred_coords[1],'Motor axis 2': pred_coords[2]})\n        if tomo_data is not None: del tomo_data; gc.collect()\n\n    print(f\"\\nFinished default predictions. Total time: {time.time() - test_predict_start_time:.2f} seconds.\")\n\n# Create dataframe from the default predictions\nif predictions:\n    submission_df = pd.DataFrame(predictions)\n    print(\"\\nDefault Predictions for the specified IDs:\")\n    print(submission_df)\nelse:\n    print(\"\\nNo default predictions were generated (check specific IDs).\")\n    submission_df = pd.DataFrame(columns=['tomo_id', 'Motor axis 0', 'Motor axis 1', 'Motor axis 2'])\n\n\n# %% [markdown]\n# # 8. Create Submission File\n#\n# Generate the submission file template, ensuring all original test IDs are present.\n\n# %%\nprint(\"\\nCreating Submission File Template...\")\nif 'submission_df' in locals():\n    try:\n        sample_sub = pd.read_csv(SAMPLE_SUB_PATH)\n        expected_test_ids = test_ids_found # Use all IDs found in the test directory\n\n        if not expected_test_ids:\n             print(\"Warning: No test IDs were found in the test directory structure.\")\n             expected_test_ids = sample_sub['tomo_id'].tolist()\n             print(f\"Using {len(expected_test_ids)} IDs from sample submission as base.\")\n        else:\n             print(f\"Generating submission structure based on {len(expected_test_ids)} test IDs found in directory.\")\n\n        sub_df_final = pd.DataFrame({'tomo_id': expected_test_ids})\n\n        # Merge the default predictions we generated for the *specific* IDs\n        if not submission_df.empty:\n             sub_df_final = sub_df_final.merge(submission_df, on='tomo_id', how='left')\n        else:\n             # If even the default predictions weren't made, fill everything\n             for col in ['Motor axis 0', 'Motor axis 1', 'Motor axis 2']: sub_df_final[col] = NO_MOTOR_COORD\n\n        # Fill any remaining NaNs (for test IDs not in the specific prediction list) with default\n        sub_df_final.fillna(NO_MOTOR_COORD, inplace=True)\n\n        # Ensure column order and types match sample\n        final_columns = sample_sub.columns.tolist()\n        for col in final_columns:\n             if col not in sub_df_final.columns:\n                  print(f\"Warning: Column '{col}' from sample not in generated df. Adding default.\")\n                  sub_df_final[col] = NO_MOTOR_COORD\n        sub_df_final = sub_df_final[final_columns]\n        for col in ['Motor axis 0', 'Motor axis 1', 'Motor axis 2']:\n            sub_df_final[col] = pd.to_numeric(sub_df_final[col], errors='coerce').fillna(NO_MOTOR_COORD).astype(float)\n\n        sub_df_final.to_csv('submission.csv', index=False)\n        print(\"Submission file template created: submission.csv\")\n        print(\"\\nFinal Submission Head:\")\n        print(sub_df_final.head())\n        print(f\"Submission shape: {sub_df_final.shape}\")\n\n        if len(sub_df_final) != len(expected_test_ids):\n            print(f\"CRITICAL WARNING: Submission row count ({len(sub_df_final)}) \"\n                  f\"mismatches expected test IDs ({len(expected_test_ids)})!\")\n        else:\n             print(\"Submission row count matches expected test IDs.\")\n\n    except FileNotFoundError:\n        print(f\"Error: Sample submission file not found at {SAMPLE_SUB_PATH}. Cannot verify format.\")\n        # Try to save anyway if df exists\n        if 'sub_df_final' in locals():\n             sub_df_final.to_csv('submission.csv', index=False)\n             print(\"Submission file template created (format not verified): submission.csv\")\n    except Exception as e:\n        print(f\"Error creating final submission file: {e}\")\nelse:\n    print(\"No prediction DataFrame ('submission_df') was generated. Cannot create submission template.\")\n\n\n# %% [markdown]\n# # 9. Discussion and Future Work\n#\n# *   **Data Loading & Handling:** Code structure maintained for loading data, labels, and finding directories.\n# *   **Feature Engineering:** Utility functions (`preprocess_tomogram`, `extract_features_patch`) remain defined but were not used for training.\n# *   **Modeling:** **Training data generation and model training sections were entirely skipped.** `clf_model` and `reg_model` were explicitly set to `None`.\n# *   **Targeted Prediction:** The prediction loop for specific IDs (`tomo_003acc`, `tomo_00e047`, `tomo_00e463`) was executed. However, since the models were `None`, the `predict_motor_location` function returned the default `NO_MOTOR_COORD` for these IDs.\n# *   **Submission Generation:** A `submission.csv` file was successfully created containing *all* test IDs found in the `test` directory. The specifically targeted IDs have `NO_MOTOR_COORD` because no model was trained, and all other test IDs were also filled with `NO_MOTOR_COORD`. This serves as a valid submission template.\n# *   **Limitations:** This script produces a baseline submission with no actual detection. Its primary purpose is to verify the data loading, directory handling, and submission file formatting logic.\n# *   **Future Improvements:** To get actual predictions, the skipping of Sections 5 and 5.1 must be removed, and sufficient training data (likely more than the 20 used in the previous version) needs to be processed to train meaningful models. Then, the prediction pipeline in Section 7 will use the trained models. All other future improvements mentioned previously (CNNs, FCNs, augmentation, etc.) still apply for building a competitive model.\n\n# %%\nprint(\"Script finished (Training and Model Fitting Skipped). Submission template generated.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T22:48:35.522190Z","iopub.execute_input":"2025-04-16T22:48:35.522433Z","iopub.status.idle":"2025-04-16T22:49:16.813623Z","shell.execute_reply.started":"2025-04-16T22:48:35.522418Z","shell.execute_reply":"2025-04-16T22:49:16.813051Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -*- coding: utf-8 -*-\n\"\"\"\nKaggle Notebook: Calculate and Visualize Physics-Inspired Features for Specific Tomograms\n\nGoal: To load specific tomograms ('tomo_003acc', 'tomo_00e047', 'tomo_01a877'),\n      extract a patch around their geometric center (as no motor is labeled),\n      calculate the physics-inspired features\n      (Smoothed, Gradient Mag, ST Coherence, ST λ_max, Hessian λ_min, Blobness),\n      and display these calculated features in a 6-panel plot for each tomogram.\n\"\"\"\n\n# %% [markdown]\n# # 1. Setup and Imports\n#\n# Load necessary libraries for data loading, calculation, and visualization.\n\n# %%\n# Attempt installation only if needed (e.g., in Kaggle environment)\ntry:\n    import cv2\n    import mrcfile # Although not used for loading slices here\n    print(\"Libraries opencv-python, mrcfile seem available.\")\nexcept ImportError:\n    print(\"Installing required libraries: opencv-python, mrcfile\")\n    %pip install opencv-python mrcfile --quiet\n    import cv2\n    # import mrcfile # Not strictly needed if only using cv2 for slices\n    print(\"Libraries installed.\")\n\nimport os\nimport glob\nimport numpy as np\nimport pandas as pd\nimport cv2 # Import OpenCV\nimport matplotlib.pyplot as plt\nfrom scipy import ndimage as ndi\nfrom skimage.feature import structure_tensor, structure_tensor_eigenvalues, hessian_matrix, hessian_matrix_eigvals\nfrom skimage.filters import gaussian\nimport gc # Garbage collection\nimport time # For timing\n\n# Constants\nCOMPETITION_DIR = '/kaggle/input/byu-locating-bacterial-flagellar-motors-2025'\nTRAIN_DIR = os.path.join(COMPETITION_DIR, 'train')\nTEST_DIR = os.path.join(COMPETITION_DIR, 'test') # Include test dir as some IDs might be there\nTRAIN_LABELS_PATH = os.path.join(COMPETITION_DIR, 'train_labels.csv')\n\n# Default Voxel spacing (fallback if lookup fails)\nDEFAULT_VOXEL_SPACING_A = 10.0\n\n# Target Tomogram IDs (those with no labeled motors)\nTARGET_TOMO_IDS = ['tomo_003acc', 'tomo_00e047', 'tomo_01a877']\n\n# Visualization Parameters\nPATCH_SIZE_PX_VIS = 96        # Patch size for visualization\nSIGMA_SMOOTH = 1.5            # Sigma for initial smoothing\nSIGMA_ST = 2.0                # Sigma for Structure Tensor calculations\nSIGMA_HESSIAN = 3.0           # Sigma for Hessian calculations\n\nprint(f\"Target Tomograms for Visualization: {TARGET_TOMO_IDS}\")\nprint(f\"Patch Size: {PATCH_SIZE_PX_VIS}x{PATCH_SIZE_PX_VIS}x{PATCH_SIZE_PX_VIS}\")\n\n# %% [markdown]\n# # 2. Load Metadata and Prepare Directory Maps\n#\n# Load the labels to get voxel spacing (if available) for our target tomograms. Create maps for both train and test directories.\n\n# %%\nprint(\"Loading labels (primarily for voxel spacing)...\")\nvoxel_spacing_map = {}\ntry:\n    train_labels_df = pd.read_csv(TRAIN_LABELS_PATH)\n    print(f\"Training labels shape: {train_labels_df.shape}\")\n    # Create map for Voxel Spacing from labels (use first entry per tomo_id)\n    if 'Voxel spacing' in train_labels_df.columns and not train_labels_df.empty:\n        voxel_spacing_map = train_labels_df.drop_duplicates(subset='tomo_id').set_index('tomo_id')['Voxel spacing'].to_dict()\n        print(f\"Loaded voxel spacing for {len(voxel_spacing_map)} tomograms.\")\n    else:\n        print(\"WARNING: 'Voxel spacing' column not found or labels empty. Will use default spacing.\")\n\nexcept FileNotFoundError:\n    print(f\"ERROR: Training labels file not found at {TRAIN_LABELS_PATH}. Will use default spacing.\")\n    train_labels_df = pd.DataFrame()\n\n# --- Find Tomogram Directories ---\nprint(\"\\nFinding tomogram directories...\")\ntrain_dirs = sorted(glob.glob(os.path.join(TRAIN_DIR, 'tomo_*')))\ntest_dirs = sorted(glob.glob(os.path.join(TEST_DIR, 'tomo_*')))\n\ntrain_dir_map = {os.path.basename(d): d for d in train_dirs if os.path.isdir(d)}\ntest_dir_map = {os.path.basename(d): d for d in test_dirs if os.path.isdir(d)}\n\n# Combine maps for easier lookup, prioritizing test if duplicates exist (unlikely)\nfull_dir_map = {**train_dir_map, **test_dir_map}\n\nprint(f\"Found {len(train_dir_map)} potential training directories.\")\nprint(f\"Found {len(test_dir_map)} potential testing directories.\")\nprint(f\"Total unique directories mapped: {len(full_dir_map)}\")\n\n\n# %% [markdown]\n# # 3. Utility Functions (Load Tomogram, Extract Patch)\n\n# %%\n# --- Tomogram Loading Function ---\ndef load_tomogram(tomo_id, dir_map):\n    \"\"\"Loads a tomogram by reading and stacking slices from a directory.\"\"\"\n    directory_path = dir_map.get(tomo_id)\n    if not directory_path or not os.path.isdir(directory_path):\n        print(f\"Warning: Tomogram directory for {tomo_id} not found in provided map.\")\n        return None\n    slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.jpg')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.png')))\n    if not slice_files: slice_files = sorted(glob.glob(os.path.join(directory_path, 'slice_*.tif')))\n    if not slice_files: print(f\"Warning: No slice files found in {directory_path}.\"); return None\n    def get_slice_number(filepath):\n        try: return int(os.path.basename(filepath).split('slice_')[1].split('.')[0])\n        except: return -1\n    valid_slice_files = [(f, get_slice_number(f)) for f in slice_files if get_slice_number(f) != -1]\n    if not valid_slice_files: print(f\"Warning: No valid slice numbers found in {directory_path}.\"); return None\n    valid_slice_files.sort(key=lambda item: item[1]); slice_files_sorted = [f for f, num in valid_slice_files]\n    try:\n        first_slice = cv2.imread(slice_files_sorted[0], cv2.IMREAD_GRAYSCALE)\n        if first_slice is None: raise IOError(\"imread failed\")\n        height, width = first_slice.shape\n        num_slices = len(slice_files_sorted); dtype = first_slice.dtype\n    except Exception as e: print(f\"Error reading first slice {slice_files_sorted[0]}: {e}\"); return None\n    tomogram_data = np.zeros((num_slices, height, width), dtype=dtype); tomogram_data[0, :, :] = first_slice\n    for i in range(1, num_slices):\n        try:\n            slice_img = cv2.imread(slice_files_sorted[i], cv2.IMREAD_GRAYSCALE)\n            if slice_img is None: raise IOError(\"imread failed\")\n            if slice_img.shape != (height, width): print(\"Shape mismatch\"); return None\n            tomogram_data[i, :, :] = slice_img\n        except Exception as e: print(f\"Error reading slice {slice_files_sorted[i]}: {e}\"); return None\n    print(f\"  Successfully loaded {tomo_id}, shape: {tomogram_data.shape}\")\n    return tomogram_data\n\n# --- Patch Extraction Function ---\ndef get_local_patch(tomo_data, center_px, patch_size_px=64):\n    \"\"\"Extracts a 3D patch centered at center_px. Patch size is in pixels.\"\"\"\n    if tomo_data is None or center_px is None: return None, None\n    center_px = np.round(center_px).astype(int); z, y, x = center_px\n    shape = tomo_data.shape\n    if not (0 <= z < shape[0] and 0 <= y < shape[1] and 0 <= x < shape[2]):\n        print(f\"Warning: Center px {center_px} out of bounds {shape}. Cannot extract patch.\")\n        return None, None\n    half_size = patch_size_px // 2\n    z_start, z_end = max(0, z - half_size), min(shape[0], z + half_size)\n    y_start, y_end = max(0, y - half_size), min(shape[1], y + half_size)\n    x_start, x_end = max(0, x - half_size), min(shape[2], x + half_size)\n    dest_z_start = half_size - (z - z_start); dest_z_end = dest_z_start + (z_end - z_start)\n    dest_y_start = half_size - (y - y_start); dest_y_end = dest_y_start + (y_end - y_start)\n    dest_x_start = half_size - (x - x_start); dest_x_end = dest_x_start + (x_end - x_start)\n    patch = np.zeros((patch_size_px, patch_size_px, patch_size_px), dtype=tomo_data.dtype)\n    try:\n        src_slice = (slice(z_start, z_end), slice(y_start, y_end), slice(x_start, x_end))\n        dest_slice = (slice(dest_z_start, dest_z_end), slice(dest_y_start, dest_y_end), slice(dest_x_start, dest_x_end))\n        patch[dest_slice] = tomo_data[src_slice]\n        # Calculate center coordinates *within the patch* (should be near half_size)\n        patch_center_coords = np.round([z - z_start + dest_z_start, y - y_start + dest_y_start, x - x_start + dest_x_start]).astype(int)\n        if not (0 <= patch_center_coords[0] < patch_size_px and 0 <= patch_center_coords[1] < patch_size_px and 0 <= patch_center_coords[2] < patch_size_px):\n             patch_center_coords = [patch_size_px // 2] * 3 # Fallback\n        print(f\"  Extracted patch shape: {patch.shape}, Center within patch: {patch_center_coords}\")\n        return patch, patch_center_coords\n    except (ValueError, IndexError) as e: print(f\"Patch extraction error: {e}\"); return None, None\n\n\n# %% [markdown]\n# # 4. Process Each Target Tomogram\n#\n# Loop through the specified IDs, load data, extract a central patch, calculate features, and visualize.\n\n# %%\nfor tomo_id in TARGET_TOMO_IDS:\n    print(f\"\\n--- Processing Tomogram: {tomo_id} ---\")\n\n    # --- Load Tomogram Data ---\n    tomo_data = load_tomogram(tomo_id, full_dir_map)\n    if tomo_data is None:\n        print(f\"  Skipping {tomo_id} due to loading error.\")\n        continue\n\n    # --- Get Geometric Center and Voxel Spacing ---\n    tomo_shape = tomo_data.shape\n    geometric_center_px = np.array([tomo_shape[0] // 2, tomo_shape[1] // 2, tomo_shape[2] // 2])\n    voxel_spacing = voxel_spacing_map.get(tomo_id, DEFAULT_VOXEL_SPACING_A)\n    print(f\"  Geometric Center (px): {geometric_center_px}\")\n    print(f\"  Voxel Spacing (A/px): {voxel_spacing}\")\n\n    # --- Extract Central Patch ---\n    print(f\"  Extracting {PATCH_SIZE_PX_VIS}x{PATCH_SIZE_PX_VIS} patch around geometric center...\")\n    central_patch, patch_center_coords = get_local_patch(tomo_data, geometric_center_px, patch_size_px=PATCH_SIZE_PX_VIS)\n\n    # Free memory of full tomogram\n    del tomo_data; gc.collect()\n\n    if central_patch is None:\n        print(f\"  Skipping {tomo_id} because patch extraction failed.\")\n        continue\n\n    # --- Calculate Features ---\n    patch_smooth = grad_mag = st_coherence = st_lambda_max = hessian_lambda_min = blobness = None\n    calculation_successful = False\n    print(\"\\n  Calculating features on the central patch...\")\n    calc_start_time = time.time()\n    try:\n        if not np.issubdtype(central_patch.dtype, np.floating):\n            patch_float = central_patch.astype(np.float32)\n        else:\n            patch_float = central_patch\n\n        patch_smooth = gaussian(patch_float, sigma=SIGMA_SMOOTH, mode='reflect', preserve_range=True, truncate=4.0)\n        grad_z, grad_y, grad_x = np.gradient(patch_smooth)\n        grad_mag = np.sqrt(grad_z**2 + grad_y**2 + grad_x**2)\n\n        S_elems = structure_tensor(patch_smooth, sigma=SIGMA_ST, mode='reflect')\n        eigvals_S = structure_tensor_eigenvalues(S_elems)\n        eigvals_S = np.sort(eigvals_S, axis=0)\n        lambda3, lambda2, lambda1 = eigvals_S[0], eigvals_S[1], eigvals_S[2]\n        denominator = lambda1 + lambda2 + lambda3 + 1e-9\n        st_coherence = np.where(denominator > 1e-8, (lambda1 - lambda3) / denominator, 0.0)\n        st_lambda_max = lambda1\n\n        h_matrix = hessian_matrix(patch_smooth, sigma=SIGMA_HESSIAN, mode='reflect', use_gaussian_derivatives=False)\n        eigvals_H = hessian_matrix_eigvals(h_matrix)\n        eigvals_H = np.sort(eigvals_H, axis=0)\n        h_lambda1, h_lambda2, h_lambda3 = eigvals_H[0], eigvals_H[1], eigvals_H[2]\n        hessian_lambda_min = h_lambda1\n        blobness = np.abs(h_lambda1 * h_lambda2 * h_lambda3) * (h_lambda1 < 0) * (h_lambda2 < 0) * (h_lambda3 < 0)\n\n        calculation_successful = True\n        print(f\"  Feature calculation took {time.time() - calc_start_time:.2f} seconds.\")\n\n    except Exception as e:\n        print(f\"  ERROR during feature calculation for {tomo_id}: {e}\")\n        calculation_successful = False\n\n    # --- Visualize Features ---\n    if calculation_successful and patch_center_coords is not None:\n        print(\"\\n  Generating visualization...\")\n\n        # Use the Z-coordinate of the center *within the patch* for slicing\n        center_z_in_patch = patch_center_coords[0]\n        center_y_in_patch = patch_center_coords[1]\n        center_x_in_patch = patch_center_coords[2]\n\n        fig, axes = plt.subplots(2, 3, figsize=(18, 11))\n        fig.suptitle(f\"Calculated Features for {tomo_id} (Central Patch, Z-Slice={center_z_in_patch})\", fontsize=16)\n\n        # Panel 1: Smoothed\n        ax = axes[0, 0]; im = ax.imshow(patch_smooth[center_z_in_patch, :, :], cmap='gray'); ax.set_title(f'Smoothed (Z={center_z_in_patch})'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n        # No marker plotted\n\n        # Panel 2: Gradient Magnitude\n        ax = axes[0, 1]; im = ax.imshow(grad_mag[center_z_in_patch, :, :], cmap='magma'); ax.set_title('Gradient Mag'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n\n        # Panel 3: ST Coherence\n        ax = axes[0, 2]; im = ax.imshow(st_coherence[center_z_in_patch, :, :], cmap='viridis', vmin=0, vmax=1); ax.set_title('ST Coherence'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n\n        # Panel 4: ST Lambda_max\n        ax = axes[1, 0]; im = ax.imshow(st_lambda_max[center_z_in_patch, :, :], cmap='plasma'); ax.set_title('ST $\\lambda_{max}$'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n\n        # Panel 5: Hessian Lambda_min\n        ax = axes[1, 1]; vmin_h = np.percentile(hessian_lambda_min, 1); vmax_h = np.percentile(hessian_lambda_min, 99); im = ax.imshow(hessian_lambda_min[center_z_in_patch, :, :], cmap='coolwarm', vmin=vmin_h, vmax=vmax_h); ax.set_title('Hessian $\\lambda_{min}$'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n\n        # Panel 6: Blobness\n        ax = axes[1, 2]; vmax_b = np.percentile(blobness[blobness > 0], 99) if np.any(blobness > 0) else 1.0; im = ax.imshow(blobness[center_z_in_patch, :, :], cmap='cubehelix', vmin=0, vmax=vmax_b); ax.set_title('Blobness'); ax.axis('off'); plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n\n        plt.tight_layout(rect=[0, 0.03, 1, 0.95])\n        plt.show()\n\n    elif central_patch is None:\n        print(f\"  Could not visualize features for {tomo_id} because patch extraction failed earlier.\")\n    else:\n        print(f\"  Could not visualize features for {tomo_id} because calculation failed.\")\n\n    # Clean up memory before next loop iteration\n    del central_patch, patch_smooth, grad_mag, st_coherence, st_lambda_max, hessian_lambda_min, blobness\n    gc.collect()\n\n\n# %% [markdown]\n# # 5. Discussion\n#\n# This notebook processed the tomograms `tomo_003acc`, `tomo_00e047`, and `tomo_01a877`. Since these tomograms did not have labeled motor coordinates in the provided `train_labels.csv` file, the analysis focused on a patch extracted from the geometric center of each volume.\n#\n# For each tomogram, the six physics-inspired features were calculated on this central patch and visualized. The plots show the characteristics of the cellular environment at the center of these specific tomograms according to each feature:\n#\n# *   **Smoothed:** Shows the general density and texture at the center.\n# *   **Gradient Mag:** Highlights edges or textures present in the central region.\n# *   **ST Coherence & λ_max:** Reveal the degree of orientation and anisotropy of structures at the center.\n# *   **Hessian λ_min & Blobness:** Indicate the local curvature and presence of any blob-like structures near the center.\n#\n# Unlike visualizations centered on a known motor, these plots provide insight into the \"background\" characteristics captured by these features in tomograms potentially lacking the target structure (or where it wasn't labeled). This can be useful for understanding feature responses away from the target.\n\n# %%\nprint(\"Notebook execution finished.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T22:59:22.856186Z","iopub.execute_input":"2025-04-16T22:59:22.857273Z","iopub.status.idle":"2025-04-16T22:59:48.669695Z","shell.execute_reply.started":"2025-04-16T22:59:22.857233Z","shell.execute_reply":"2025-04-16T22:59:48.669069Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import math # For sqrt symbol, though not strictly needed for text output\n\n# Define a helper function for clarity (optional)\ndef print_concept(title, concept, equation, usability, equation2=None):\n    \"\"\"Prints the concept, equation, and usability in a standard format.\"\"\"\n    print(f\"\\n--- {title} ---\")\n    print(f\"Concept: {concept}\")\n    print(f\"General Equation: {equation}\")\n    if equation2:\n        print(f\"Derived/Related Equation: {equation2}\")\n    print(f\"Usability: {usability}\")\n\n# --- Print Explanations ---\n\nprint(\"=== Conceptual Overview of Physics-Inspired Tomogram Features ===\")\n\nprint_concept(\n    title=\"1. Intensity: The Base Scalar Field (Φ)\",\n    concept=\"The fundamental data assigning a density/intensity value to every 3D point.\",\n    equation=\"I_s(x, y, z)  (Intensity, typically after smoothing)\",\n    usability=\"Provides the raw density map. Smoothing creates a stable base for derivatives.\"\n)\n\nprint_concept(\n    title=\"2. Gradient: Rate and Direction of Change (E)\",\n    concept=\"A vector field pointing in the direction of steepest intensity increase.\",\n    equation=\"∇I_s = (∂I_s/∂x, ∂I_s/∂y, ∂I_s/∂z)\",\n    equation2=\"|∇I_s| = sqrt[ (∂I_s/∂x)² + (∂I_s/∂y)² + (∂I_s/∂z)² ]  (Magnitude)\",\n    usability=\"Magnitude (|∇I_s|) detects edges and textures. Direction indicates orientation perpendicular to iso-surfaces.\"\n)\n\nprint_concept(\n    title=\"3. Structure Tensor: Local Orientation and Anisotropy (T_μν / S_μν)\",\n    concept=\"Tensor summarizing predominant gradient orientations via local averaging.\",\n    # Using <>_sigma notation for averaging over scale sigma\n    equation=\"ST = <∇I_s ⊗ ∇I_s^T>_σ  (Locally averaged outer product of gradient)\",\n    # Explain eigenvalues and coherence conceptually\n    equation2=\"Derived features based on Eigenvalues (λ₁≥λ₂≥λ₃≥0): λ_max=λ₁, Coherence C ≈ (λ₁-λ₃)/(λ₁+λ₂+λ₃)\",\n    usability=\"Describes local geometry: discriminates lines vs. planes vs. isotropic regions using eigenvalues. Coherence measures degree/strength of orientation.\"\n)\n\nprint_concept(\n    title=\"4. Hessian Matrix: Local Curvature\",\n    concept=\"Matrix of second-order partial derivatives describing the local curvature of the intensity landscape.\",\n    equation=\"H_ij = ∂²I_s / ∂xᵢ ∂xⱼ\",\n    # Explain eigenvalues and blobness conceptually\n    equation2=\"Derived features based on Eigenvalues (h₁, h₂, h₃): λ_min=h₁, Blobness B ≈ |h₁h₂h₃| if all h<0\",\n    usability=\"Classifies local shape (blob, tube, sheet) via eigenvalue signs/magnitudes. λ_min helps find centers of dark structures. Blobness specifically highlights blob-like regions.\"\n)\n\nprint(\"\\n==============================================================\")\nprint(\"Notebook execution finished. Conceptual explanations printed.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-16T23:08:08.078202Z","iopub.execute_input":"2025-04-16T23:08:08.078506Z","iopub.status.idle":"2025-04-16T23:08:08.086295Z","shell.execute_reply.started":"2025-04-16T23:08:08.078483Z","shell.execute_reply":"2025-04-16T23:08:08.085426Z"}},"outputs":[],"execution_count":null}]}