{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":91498,"databundleVersionId":11655853,"sourceType":"competition"}],"dockerImageVersionId":31193,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"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\n#for 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-12-08T13:09:48.771115Z","iopub.execute_input":"2025-12-08T13:09:48.771684Z","iopub.status.idle":"2025-12-08T13:09:49.041855Z","shell.execute_reply.started":"2025-12-08T13:09:48.771658Z","shell.execute_reply":"2025-12-08T13:09:49.041198Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Import Libraries and Setup","metadata":{}},{"cell_type":"code","source":"import os\nimport sys\nimport numpy as np\nimport cv2\nimport torch\nimport torch.nn.functional as F\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# Check GPU availability\nprint(f\"PyTorch version: {torch.__version__}\")\nprint(f\"CUDA available: {torch.cuda.is_available()}\")\nif torch.cuda.is_available():\n    print(f\"GPU: {torch.cuda.get_device_name(0)}\")\n    print(f\"GPU Memory: {torch.cuda.get_device_properties(0).total_memory / 1e9:.2f} GB\")\n\n# Set device\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\nprint(f\"\\nUsing device: {device}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:10:50.972033Z","iopub.execute_input":"2025-12-08T13:10:50.972957Z","iopub.status.idle":"2025-12-08T13:10:51.102717Z","shell.execute_reply.started":"2025-12-08T13:10:50.972926Z","shell.execute_reply":"2025-12-08T13:10:51.101851Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install --no-deps git+https://github.com/cvg/LightGlue.git","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:27:12.016154Z","iopub.execute_input":"2025-12-08T13:27:12.016852Z","iopub.status.idle":"2025-12-08T13:27:24.032910Z","shell.execute_reply.started":"2025-12-08T13:27:12.016823Z","shell.execute_reply":"2025-12-08T13:27:24.032037Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import cv2\nimport torch\nfrom lightglue import match_pair\nfrom lightglue import ALIKED, LightGlue\nfrom lightglue.utils import load_image, rbd\nfrom kornia.feature import LightGlue","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:28:11.707965Z","iopub.execute_input":"2025-12-08T13:28:11.708490Z","iopub.status.idle":"2025-12-08T13:28:11.712430Z","shell.execute_reply.started":"2025-12-08T13:28:11.708466Z","shell.execute_reply":"2025-12-08T13:28:11.711857Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Configuration","metadata":{}},{"cell_type":"code","source":"class Config:\n    # Paths - MODIFY THESE ACCORDING TO YOUR KAGGLE SETUP\n    DATA_PATH = \"/kaggle/input/image-matching-challenge-2025\"\n    TRAIN_PATH = f\"{DATA_PATH}/train\"\n    TEST_PATH = f\"{DATA_PATH}/test\"\n    OUTPUT_PATH = \"/kaggle/working/output\"\n    \n    # Which dataset to use\n    USE_TRAIN = True  # Set False to use test data\n    PROCESS_ALL_SCENES = True  # Set False to process single scene\n    SCENE_NAME = None  # Set to specific scene name or None to auto-detect\n    \n    # Pipeline parameters\n    TOP_K_SIMILAR = 15  # Number of most similar images to match per image\n    MIN_MATCHES = 15  # Minimum matches required for a valid pair\n    MATCH_CONFIDENCE = 0.2  # LightGlue confidence threshold\n    MAX_IMAGES = 50  # Maximum images to process per scene (set None for all)\n    \n    # Visualization\n    DISPLAY_MATCHES = True\n    NUM_MATCH_EXAMPLES = 2  # Per scene\n\n\nconfig = Config()\n# Create output directory\nos.makedirs(config.OUTPUT_PATH, exist_ok=True)\nprint(f\"Configuration loaded. Output directory: {config.OUTPUT_PATH}\")\nprint(f\"Process all scenes: {config.PROCESS_ALL_SCENES}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:24:35.439270Z","iopub.execute_input":"2025-12-08T13:24:35.439954Z","iopub.status.idle":"2025-12-08T13:24:35.445672Z","shell.execute_reply.started":"2025-12-08T13:24:35.439928Z","shell.execute_reply":"2025-12-08T13:24:35.445063Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Config:\n    # Paths - MODIFY THESE ACCORDING TO YOUR KAGGLE SETUP\n    DATA_PATH = \"/kaggle/input/image-matching-challenge-2025\"\n    TRAIN_PATH = f\"{DATA_PATH}/train\"\n    TEST_PATH = f\"{DATA_PATH}/test\"\n    OUTPUT_PATH = \"/kaggle/working/output\"\n    \n    # Which dataset to use\n    USE_TRAIN = True  # Set False to use test data\n    PROCESS_ALL_SCENES = True  # Set False to process single scene\n    SCENE_NAME = None  # Set to specific scene name or None to auto-detect\n    \n    # Pipeline parameters\n    TOP_K_SIMILAR = 15  # Number of most similar images to match per image\n    MIN_MATCHES = 15  # Minimum matches required for a valid pair\n    MATCH_CONFIDENCE = 0.2  # LightGlue confidence threshold\n    MAX_IMAGES = 50  # Maximum images to process per scene (set None for all)\n    \n    # Visualization\n    DISPLAY_MATCHES = True\n    NUM_MATCH_EXAMPLES = 2  # Per scene\n    \n    # COLMAP settings\n    USE_COLMAP = True  # Use actual COLMAP for reconstruction\n    COLMAP_VOCAB_TREE = \"/kaggle/input/colmap-vocab-tree/vocab_tree_flickr100K_words256K.bin\"  # Optional\n    \nconfig = Config()\n\n# Create output directory\nos.makedirs(config.OUTPUT_PATH, exist_ok=True)\nprint(f\"Configuration loaded. Output directory: {config.OUTPUT_PATH}\")\nprint(f\"Process all scenes: {config.PROCESS_ALL_SCENES}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:27:56.518677Z","iopub.execute_input":"2025-12-08T16:27:56.518976Z","iopub.status.idle":"2025-12-08T16:27:56.525673Z","shell.execute_reply.started":"2025-12-08T16:27:56.518956Z","shell.execute_reply":"2025-12-08T16:27:56.525022Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load and Prepare Images","metadata":{}},{"cell_type":"code","source":"def load_images_from_scene(scene_path, max_images=None):\n    \"\"\"Load images from a scene directory\"\"\"\n    image_files = []\n    valid_extensions = {'.jpg', '.jpeg', '.png', '.JPG', '.JPEG', '.PNG'}\n    \n    scene_path = Path(scene_path)\n    if not scene_path.exists():\n        raise FileNotFoundError(f\"Scene path does not exist: {scene_path}\")\n    \n    # Find all image files\n    for ext in valid_extensions:\n        image_files.extend(list(scene_path.glob(f\"*{ext}\")))\n    \n    image_files = sorted(image_files)\n    \n    if max_images:\n        image_files = image_files[:max_images]\n    \n    if len(image_files) == 0:\n        raise ValueError(f\"No images found in {scene_path}\")\n    \n    print(f\"  Found {len(image_files)} images in scene\")\n    \n    # Load images\n    images = []\n    image_names = []\n    \n    for img_path in image_files:\n        try:\n            img = cv2.imread(str(img_path))\n            if img is not None:\n                img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n                images.append(img)\n                image_names.append(img_path.name)\n        except Exception as e:\n            print(f\"  Error loading {img_path}: {e}\")\n    \n    print(f\"  Successfully loaded {len(images)} images\")\n    return images, image_names, image_files\n\n# Get all scenes to process\nbase_path = config.TRAIN_PATH if config.USE_TRAIN else config.TEST_PATH\n\nif config.PROCESS_ALL_SCENES:\n    # Get all scenes\n    all_scenes = sorted([d for d in Path(base_path).iterdir() if d.is_dir()])\n    if len(all_scenes) == 0:\n        raise ValueError(f\"No scenes found in {base_path}\")\n    print(f\"Found {len(all_scenes)} scenes to process:\")\n    for scene in all_scenes:\n        print(f\"  - {scene.name}\")\nelse:\n    # Single scene\n    if config.SCENE_NAME is None:\n        scenes = [d for d in Path(base_path).iterdir() if d.is_dir()]\n        if len(scenes) == 0:\n            raise ValueError(f\"No scenes found in {base_path}\")\n        all_scenes = [scenes[0]]\n        print(f\"Auto-detected scene: {all_scenes[0].name}\")\n    else:\n        all_scenes = [Path(base_path) / config.SCENE_NAME]\n\n# Store all scene data\nscenes_data = {}\n\nprint(f\"\\n{'='*70}\")\nprint(\"LOADING ALL SCENES\")\nprint(f\"{'='*70}\\n\")\n\nfor scene_path in all_scenes:\n    scene_name = scene_path.name\n    print(f\"Loading scene: {scene_name}\")\n    \n    try:\n        images, image_names, image_paths = load_images_from_scene(\n            scene_path, \n            max_images=config.MAX_IMAGES\n        )\n        \n        scenes_data[scene_name] = {\n            'images': images,\n            'image_names': image_names,\n            'image_paths': image_paths,\n            'scene_path': scene_path\n        }\n        \n        print(f\"  ✓ Loaded {len(images)} images from '{scene_name}'\\n\")\n    \n    except Exception as e:\n        print(f\"  ✗ Error loading scene '{scene_name}': {e}\\n\")\n        continue\n\nprint(f\"\\n✓ Successfully loaded {len(scenes_data)} scenes\")\nprint(f\"  Total images: {sum(len(data['images']) for data in scenes_data.values())}\")\n\n# Display sample from first scene\nif scenes_data:\n    first_scene = list(scenes_data.keys())[0]\n    sample_images = scenes_data[first_scene]['images']\n    sample_names = scenes_data[first_scene]['image_names']\n    \n    fig, axes = plt.subplots(1, min(4, len(sample_images)), figsize=(15, 4))\n    if len(sample_images) == 1:\n        axes = [axes]\n    for i, ax in enumerate(axes):\n        if i < len(sample_images):\n            ax.imshow(sample_images[i])\n            ax.set_title(f\"{sample_names[i]}\")\n            ax.axis('off')\n    plt.suptitle(f\"Sample from scene: {first_scene}\")\n    plt.tight_layout()\n    plt.savefig(f\"{config.OUTPUT_PATH}/sample_images.png\", dpi=150, bbox_inches='tight')\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:25:36.459073Z","iopub.execute_input":"2025-12-08T13:25:36.459459Z","iopub.status.idle":"2025-12-08T13:26:41.476559Z","shell.execute_reply.started":"2025-12-08T13:25:36.459423Z","shell.execute_reply":"2025-12-08T13:26:41.475652Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Process Each Scene","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"PROCESSING ALL SCENES\")\nprint(\"=\"*70 + \"\\n\")\n\n# Initialize models once (reuse across scenes)\nprint(\"Loading models...\")\n\n# Load DINOv2\ntry:\n    dinov2_model = torch.hub.load('facebookresearch/dinov2', 'dinov2_vitb14')\n    dinov2_model = dinov2_model.to(device).eval()\n    print(\"✓ DINOv2 loaded\")\nexcept Exception as e:\n    print(f\"✗ Error loading DINOv2: {e}\")\n    dinov2_model = None\n\naliked = ALIKED(\n        max_num_keypoints=2048,\n        detection_threshold=0.01,\n        resize=1024\n    ).to(device).eval()\n    \nprint(\"✓ ALIKED loaded successfully\")\n\nlightglue = LightGlue(\n        features='aliked',\n        depth_confidence=0.95,\n        width_confidence=0.95\n    ).to(device).eval()\n    \nprint(\"✓ LightGlue loaded successfully\")\n\nprint()\n\ndef preprocess_for_dinov2(image, size=224):\n    h, w = image.shape[:2]\n    if h > w:\n        new_h, new_w = size, int(w * size / h)\n    else:\n        new_h, new_w = int(h * size / w), size\n    resized = cv2.resize(image, (new_w, new_h))\n    pad_h, pad_w = size - new_h, size - new_w\n    padded = np.pad(resized, ((pad_h//2, pad_h-pad_h//2), (pad_w//2, pad_w-pad_w//2), (0, 0)), mode='constant')\n    img_tensor = torch.from_numpy(padded).float() / 255.0\n    img_tensor = img_tensor.permute(2, 0, 1)\n    mean = torch.tensor([0.485, 0.456, 0.406]).view(3, 1, 1)\n    std = torch.tensor([0.229, 0.224, 0.225]).view(3, 1, 1)\n    return (img_tensor - mean) / std\n\ndef extract_features_opencv(image, max_features=2048):\n    gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)\n    try:\n        detector = cv2.SIFT_create(nfeatures=max_features)\n    except:\n        detector = cv2.ORB_create(nfeatures=max_features)\n    keypoints, descriptors = detector.detectAndCompute(gray, None)\n    if descriptors is None:\n        return np.array([]), np.array([])\n    return np.array([kp.pt for kp in keypoints]), descriptors\n\ndef match_features_opencv(desc1, desc2):\n    if len(desc1) == 0 or len(desc2) == 0:\n        return np.array([])\n    bf = cv2.BFMatcher(cv2.NORM_L2 if desc1.dtype == np.float32 else cv2.NORM_HAMMING, crossCheck=False)\n    matches = bf.knnMatch(desc1, desc2, k=2)\n    good_matches = []\n    for pair in matches:\n        if len(pair) == 2:\n            m, n = pair\n            if m.distance < 0.75 * n.distance:\n                good_matches.append([m.queryIdx, m.trainIdx])\n    return np.array(good_matches)\n\ndef match_image_pair(img1, img2):\n    if aliked is not None and lightglue is not None:\n        try:\n            with torch.no_grad():\n                img1_t = torch.from_numpy(img1).float().permute(2, 0, 1).unsqueeze(0).to(device) / 255.0\n                img2_t = torch.from_numpy(img2).float().permute(2, 0, 1).unsqueeze(0).to(device) / 255.0\n                feats1 = aliked(img1_t)\n                feats2 = aliked(img2_t)\n                matches_dict = lightglue({'image0': feats1, 'image1': feats2})\n                matches = matches_dict['matches0'][0].cpu().numpy()\n                valid = matches > -1\n                kpts0 = feats1['keypoints'][0].cpu().numpy()\n                kpts1 = feats2['keypoints'][0].cpu().numpy()\n                return kpts0[valid], kpts1[matches[valid].astype(int)]\n        except Exception as e:\n            pass\n    kpts1, desc1 = extract_features_opencv(img1)\n    kpts2, desc2 = extract_features_opencv(img2)\n    matches = match_features_opencv(desc1, desc2)\n    if len(matches) == 0:\n        return np.array([]), np.array([])\n    return kpts1[matches[:, 0]], kpts2[matches[:, 1]]\n\n# Process each scene\nall_results = {}\n\nfor scene_idx, (scene_name, scene_data) in enumerate(scenes_data.items()):\n    print(f\"\\n{'='*70}\")\n    print(f\"SCENE {scene_idx+1}/{len(scenes_data)}: {scene_name}\")\n    print(f\"{'='*70}\\n\")\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    \n    # Extract embeddings\n    if dinov2_model is not None:\n        print(\"Extracting DINOv2 embeddings...\")\n        embeddings = []\n        with torch.no_grad():\n            for i in range(0, len(images), 8):\n                batch = torch.stack([preprocess_for_dinov2(img) for img in images[i:i+8]]).to(device)\n                embeddings.append(dinov2_model(batch).cpu())\n                print(f\"  {min(i+8, len(images))}/{len(images)}\", end='\\r')\n        embeddings = F.normalize(torch.cat(embeddings, dim=0), p=2, dim=1)\n        similarity_matrix = torch.mm(embeddings, embeddings.t()).numpy()\n        print(f\"\\n✓ Extracted embeddings\")\n    else:\n        print(\"⚠ Skipping DINOv2, using sequential pairs\")\n        similarity_matrix = np.eye(len(images))\n    \n    # Create pairs\n    image_pairs = []\n    for i in range(len(images)):\n        similarities = similarity_matrix[i].copy()\n        similarities[i] = -1\n        top_k = np.argsort(similarities)[::-1][:config.TOP_K_SIMILAR]\n        for j in top_k:\n            if j > i:\n                image_pairs.append((i, j))\n    print(f\"✓ Created {len(image_pairs)} pairs\")\n    \n    # Match pairs\n    print(\"Matching pairs...\")\n    all_matches = {}\n    for idx, (i, j) in enumerate(image_pairs):\n        try:\n            mkpts0, mkpts1 = match_image_pair(images[i], images[j])\n            if len(mkpts0) >= config.MIN_MATCHES:\n                all_matches[(i, j)] = {'mkpts0': mkpts0, 'mkpts1': mkpts1, 'num_matches': len(mkpts0)}\n            if (idx+1) % 20 == 0:\n                print(f\"  {idx+1}/{len(image_pairs)}\", end='\\r')\n        except:\n            continue\n    print(f\"\\n✓ Matched {len(all_matches)} pairs\")\n    \n    # Geometric verification\n    print(\"Geometric verification...\")\n    verified_matches = {}\n    for (i, j), match_data in all_matches.items():\n        mkpts0, mkpts1 = match_data['mkpts0'], match_data['mkpts1']\n        if len(mkpts0) < 8:\n            continue\n        try:\n            E, mask = cv2.findEssentialMat(mkpts0, mkpts1, focal=1.0, pp=(0., 0.), method=cv2.RANSAC, prob=0.999, threshold=1.0)\n            if E is not None and mask is not None:\n                inliers = mask.ravel() == 1\n                if np.sum(inliers) >= config.MIN_MATCHES:\n                    verified_matches[(i, j)] = {\n                        'mkpts0': mkpts0[inliers],\n                        'mkpts1': mkpts1[inliers],\n                        'num_inliers': np.sum(inliers)\n                    }\n        except:\n            continue\n    print(f\"✓ Verified {len(verified_matches)} pairs\")\n    \n    all_results[scene_name] = {\n        'verified_matches': verified_matches,\n        'similarity_matrix': similarity_matrix,\n        'images': images,\n        'image_names': image_names,\n        'scene_path': scene_data['scene_path']\n    }\n    \n    torch.cuda.empty_cache()\n\nprint(f\"\\n✓ Processed all {len(all_results)} scenes\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T13:29:47.854082Z","iopub.execute_input":"2025-12-08T13:29:47.854395Z","iopub.status.idle":"2025-12-08T15:52:20.773350Z","shell.execute_reply.started":"2025-12-08T13:29:47.854375Z","shell.execute_reply":"2025-12-08T15:52:20.772312Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Extract DinoV2 Embeddings","metadata":{}},{"cell_type":"code","source":"def preprocess_for_dinov2(image, size=224):\n    \"\"\"Preprocess image for DINOv2\"\"\"\n    # Resize\n    h, w = image.shape[:2]\n    if h > w:\n        new_h, new_w = size, int(w * size / h)\n    else:\n        new_h, new_w = int(h * size / w), size\n    \n    resized = cv2.resize(image, (new_w, new_h))\n    \n    # Pad to square\n    pad_h = size - new_h\n    pad_w = size - new_w\n    padded = np.pad(resized, ((pad_h//2, pad_h - pad_h//2), \n                               (pad_w//2, pad_w - pad_w//2), \n                               (0, 0)), mode='constant')\n    \n    # Normalize\n    img_tensor = torch.from_numpy(padded).float() / 255.0\n    img_tensor = img_tensor.permute(2, 0, 1)\n    \n    # ImageNet normalization\n    mean = torch.tensor([0.485, 0.456, 0.406]).view(3, 1, 1)\n    std = torch.tensor([0.229, 0.224, 0.225]).view(3, 1, 1)\n    img_tensor = (img_tensor - mean) / std\n    \n    return img_tensor\n\nprint(\"Extracting DINOv2 embeddings...\")\n\nembeddings = []\nbatch_size = 8\n\nwith torch.no_grad():\n    for i in range(0, len(images), batch_size):\n        batch_images = images[i:i+batch_size]\n        batch_tensors = torch.stack([\n            preprocess_for_dinov2(img) for img in batch_images\n        ]).to(device)\n        \n        batch_embeddings = dinov2_model(batch_tensors)\n        embeddings.append(batch_embeddings.cpu())\n        \n        print(f\"Processed {min(i+batch_size, len(images))}/{len(images)} images\", end='\\r')\n\nembeddings = torch.cat(embeddings, dim=0)\nembeddings = F.normalize(embeddings, p=2, dim=1)\n\nprint(f\"\\n✓ Extracted embeddings: {embeddings.shape}\")\n\n# Free up memory\ndel dinov2_model\ntorch.cuda.empty_cache()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T15:52:26.812027Z","iopub.execute_input":"2025-12-08T15:52:26.812837Z","iopub.status.idle":"2025-12-08T15:52:27.731130Z","shell.execute_reply.started":"2025-12-08T15:52:26.812811Z","shell.execute_reply":"2025-12-08T15:52:27.730461Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Compute Image Similarity and Create Pairs","metadata":{}},{"cell_type":"code","source":"print(\"Computing image similarity matrix...\")\n\n# Compute similarity matrix\nsimilarity_matrix = torch.mm(embeddings, embeddings.t()).numpy()\n\n# Create image pairs based on similarity\nimage_pairs = []\npair_similarities = []\n\nfor i in range(len(images)):\n    # Get top-k most similar images (excluding self)\n    similarities = similarity_matrix[i].copy()\n    similarities[i] = -1  # Exclude self\n    \n    top_k_indices = np.argsort(similarities)[::-1][:config.TOP_K_SIMILAR]\n    \n    for j in top_k_indices:\n        if j > i:  # Avoid duplicates\n            image_pairs.append((i, j))\n            pair_similarities.append(similarities[j])\n\nprint(f\"✓ Created {len(image_pairs)} image pairs for matching\")\n\n# Visualize similarity matrix\nplt.figure(figsize=(10, 8))\nplt.imshow(similarity_matrix, cmap='viridis', aspect='auto')\nplt.colorbar(label='Cosine Similarity')\nplt.title('Image Similarity Matrix (DINOv2)')\nplt.xlabel('Image Index')\nplt.ylabel('Image Index')\nplt.tight_layout()\nplt.savefig(f\"{config.OUTPUT_PATH}/similarity_matrix.png\", dpi=150, bbox_inches='tight')\nplt.show()\n\n# Show top pairs\nprint(\"\\nTop 5 most similar image pairs:\")\nsorted_pairs = sorted(zip(image_pairs, pair_similarities), key=lambda x: x[1], reverse=True)\nfor (i, j), sim in sorted_pairs[:5]:\n    print(f\"  {image_names[i]} <-> {image_names[j]} (similarity: {sim:.3f})\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T15:52:36.787743Z","iopub.execute_input":"2025-12-08T15:52:36.788504Z","iopub.status.idle":"2025-12-08T15:52:37.568027Z","shell.execute_reply.started":"2025-12-08T15:52:36.788478Z","shell.execute_reply":"2025-12-08T15:52:37.567163Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature Detection and Matching","metadata":{}},{"cell_type":"code","source":"def extract_features_opencv(image, max_features=2048):\n    \"\"\"Fallback feature extraction using OpenCV\"\"\"\n    gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)\n    \n    # Try SIFT first, fallback to ORB\n    try:\n        detector = cv2.SIFT_create(nfeatures=max_features)\n    except:\n        detector = cv2.ORB_create(nfeatures=max_features)\n    \n    keypoints, descriptors = detector.detectAndCompute(gray, None)\n    \n    if descriptors is None:\n        return np.array([]), np.array([])\n    \n    kpts = np.array([kp.pt for kp in keypoints])\n    return kpts, descriptors\n\ndef match_features_opencv(desc1, desc2):\n    \"\"\"Fallback matching using OpenCV\"\"\"\n    if len(desc1) == 0 or len(desc2) == 0:\n        return np.array([])\n    \n    # Use BFMatcher with cross-check\n    if desc1.dtype == np.float32:\n        bf = cv2.BFMatcher(cv2.NORM_L2, crossCheck=False)\n    else:\n        bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=False)\n    \n    matches = bf.knnMatch(desc1, desc2, k=2)\n    \n    # Lowe's ratio test\n    good_matches = []\n    for pair in matches:\n        if len(pair) == 2:\n            m, n = pair\n            if m.distance < 0.75 * n.distance:\n                good_matches.append([m.queryIdx, m.trainIdx])\n    \n    return np.array(good_matches)\n\ndef match_image_pair(img1, img2, use_aliked=True):\n    \"\"\"Match a pair of images\"\"\"\n    if use_aliked and aliked is not None and lightglue is not None:\n        # Use ALIKED + LightGlue\n        try:\n            with torch.no_grad():\n                # Prepare images\n                h1, w1 = img1.shape[:2]\n                h2, w2 = img2.shape[:2]\n                \n                img1_t = torch.from_numpy(img1).float().permute(2, 0, 1).unsqueeze(0) / 255.0\n                img2_t = torch.from_numpy(img2).float().permute(2, 0, 1).unsqueeze(0) / 255.0\n                \n                img1_t = img1_t.to(device)\n                img2_t = img2_t.to(device)\n                \n                # Extract features\n                feats1 = aliked(img1_t)\n                feats2 = aliked(img2_t)\n                \n                # Match with LightGlue\n                matches_dict = lightglue({'image0': feats1, 'image1': feats2})\n                \n                # Extract matched keypoints - CORRECT indexing for Kornia ALIKED\n                # Check the actual structure first\n                if 'matches0' in matches_dict:\n                    matches = matches_dict['matches0']\n                    # Handle batch dimension properly\n                    if matches.dim() == 2:  # Already unbatched or shape [N]\n                        matches = matches.squeeze() if matches.dim() > 1 else matches\n                    elif matches.dim() == 1:  # Already 1D\n                        pass\n                    else:  # Has batch dimension\n                        matches = matches[0]\n                    \n                    matches = matches.cpu().numpy()\n                    valid = matches > -1\n                    \n                    # Get keypoints - handle different possible structures\n                    if 'keypoints' in feats1:\n                        kpts0 = feats1['keypoints']\n                        if kpts0.dim() == 3:  # [B, N, 2]\n                            kpts0 = kpts0[0]\n                        kpts0 = kpts0.cpu().numpy()\n                    else:\n                        raise KeyError(\"No keypoints in features\")\n                    \n                    if 'keypoints' in feats2:\n                        kpts1 = feats2['keypoints']\n                        if kpts1.dim() == 3:  # [B, N, 2]\n                            kpts1 = kpts1[0]\n                        kpts1 = kpts1.cpu().numpy()\n                    else:\n                        raise KeyError(\"No keypoints in features\")\n                    \n                    # Get matched points\n                    mkpts0 = kpts0[valid]\n                    mkpts1 = kpts1[matches[valid].astype(int)]\n                    \n                    return mkpts0, mkpts1\n                else:\n                    raise KeyError(\"No matches0 in output\")\n        \n        except Exception as e:\n            # Silently fall back to OpenCV - comment this line to debug\n            # print(f\"\\nALIKED failed, using OpenCV fallback: {e}\")\n            pass\n    \n    # Use OpenCV fallback\n    kpts1, desc1 = extract_features_opencv(img1)\n    kpts2, desc2 = extract_features_opencv(img2)\n    \n    matches = match_features_opencv(desc1, desc2)\n    \n    if len(matches) == 0:\n        return np.array([]), np.array([])\n    \n    mkpts0 = kpts1[matches[:, 0]]\n    mkpts1 = kpts2[matches[:, 1]]\n    \n    return mkpts0, mkpts1\n\nprint(\"Matching image pairs...\")\n\nall_matches = {}\nmatch_counts = []\n\nfor idx, (i, j) in enumerate(image_pairs):\n    try:\n        mkpts0, mkpts1 = match_image_pair(images[i], images[j])\n        \n        if len(mkpts0) >= config.MIN_MATCHES:\n            all_matches[(i, j)] = {\n                'mkpts0': mkpts0,\n                'mkpts1': mkpts1,\n                'num_matches': len(mkpts0)\n            }\n            match_counts.append(len(mkpts0))\n        \n        if (idx + 1) % 10 == 0:\n            print(f\"Matched {idx+1}/{len(image_pairs)} pairs\", end='\\r')\n            \n    except Exception as e:\n        print(f\"\\nError matching pair ({i}, {j}): {e}\")\n        continue\n\nprint(f\"\\n✓ Successfully matched {len(all_matches)} pairs\")\nprint(f\"  Average matches per pair: {np.mean(match_counts):.1f}\")\nprint(f\"  Max matches: {np.max(match_counts) if match_counts else 0}\")\nprint(f\"  Min matches: {np.min(match_counts) if match_counts else 0}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:01:58.704876Z","iopub.execute_input":"2025-12-08T16:01:58.705521Z","iopub.status.idle":"2025-12-08T16:06:49.653885Z","shell.execute_reply.started":"2025-12-08T16:01:58.705498Z","shell.execute_reply":"2025-12-08T16:06:49.653127Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualize Feature Matches","metadata":{}},{"cell_type":"code","source":"def draw_matches(img1, img2, mkpts0, mkpts1, num_display=50):\n    \"\"\"Draw matches between two images\"\"\"\n    h1, w1 = img1.shape[:2]\n    h2, w2 = img2.shape[:2]\n    \n    # Create side-by-side image\n    h = max(h1, h2)\n    w = w1 + w2\n    canvas = np.zeros((h, w, 3), dtype=np.uint8)\n    canvas[:h1, :w1] = img1\n    canvas[:h2, w1:w1+w2] = img2\n    \n    # Randomly sample matches to display\n    if len(mkpts0) > num_display:\n        indices = np.random.choice(len(mkpts0), num_display, replace=False)\n        mkpts0_display = mkpts0[indices]\n        mkpts1_display = mkpts1[indices]\n    else:\n        mkpts0_display = mkpts0\n        mkpts1_display = mkpts1\n    \n    # Draw matches\n    for pt1, pt2 in zip(mkpts0_display, mkpts1_display):\n        pt1 = tuple(pt1.astype(int))\n        pt2 = tuple((pt2 + np.array([w1, 0])).astype(int))\n        \n        color = tuple(np.random.randint(0, 255, 3).tolist())\n        cv2.circle(canvas, pt1, 3, color, -1)\n        cv2.circle(canvas, pt2, 3, color, -1)\n        cv2.line(canvas, pt1, pt2, color, 1)\n    \n    return canvas\n\nif config.DISPLAY_MATCHES and len(all_matches) > 0:\n    print(\"\\nVisualizing matches...\")\n    \n    # Get top matches by count\n    sorted_matches = sorted(all_matches.items(), \n                           key=lambda x: x[1]['num_matches'], \n                           reverse=True)\n    \n    num_to_display = min(config.NUM_MATCH_EXAMPLES, len(sorted_matches))\n    \n    fig, axes = plt.subplots(num_to_display, 1, figsize=(20, 6*num_to_display))\n    if num_to_display == 1:\n        axes = [axes]\n    \n    for idx, ((i, j), match_data) in enumerate(sorted_matches[:num_to_display]):\n        match_img = draw_matches(\n            images[i], images[j],\n            match_data['mkpts0'], match_data['mkpts1']\n        )\n        \n        axes[idx].imshow(match_img)\n        axes[idx].set_title(\n            f\"{image_names[i]} ↔ {image_names[j]} \"\n            f\"({match_data['num_matches']} matches)\",\n            fontsize=14\n        )\n        axes[idx].axis('off')\n    \n    plt.tight_layout()\n    plt.savefig(f\"{config.OUTPUT_PATH}/feature_matches.png\", dpi=150, bbox_inches='tight')\n    plt.show()\n    \n    print(f\"✓ Visualized {num_to_display} match examples\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:08:17.503385Z","iopub.execute_input":"2025-12-08T16:08:17.503734Z","iopub.status.idle":"2025-12-08T16:08:22.827845Z","shell.execute_reply.started":"2025-12-08T16:08:17.503709Z","shell.execute_reply":"2025-12-08T16:08:22.826648Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Geometric Verification with RANSAC","metadata":{}},{"cell_type":"code","source":"print(\"Performing geometric verification...\")\n\nverified_matches = {}\n\nfor (i, j), match_data in all_matches.items():\n    mkpts0 = match_data['mkpts0']\n    mkpts1 = match_data['mkpts1']\n    \n    if len(mkpts0) < 8:\n        continue\n    \n    try:\n        # Estimate essential matrix\n        E, inlier_mask = cv2.findEssentialMat(\n            mkpts0, mkpts1,\n            focal=1.0,\n            pp=(0., 0.),\n            method=cv2.RANSAC,\n            prob=0.999,\n            threshold=1.0\n        )\n        \n        if E is not None and inlier_mask is not None:\n            inliers = inlier_mask.ravel() == 1\n            num_inliers = np.sum(inliers)\n            \n            if num_inliers >= config.MIN_MATCHES:\n                verified_matches[(i, j)] = {\n                    'mkpts0': mkpts0[inliers],\n                    'mkpts1': mkpts1[inliers],\n                    'num_inliers': num_inliers,\n                    'E': E\n                }\n    \n    except Exception as e:\n        continue\n\nprint(f\"✓ Verified {len(verified_matches)} pairs with geometric consistency\")\n\nif len(verified_matches) == 0:\n    print(\"⚠ Warning: No geometrically verified matches found!\")\n    print(\"  The reconstruction may fail. Consider:\")\n    print(\"  - Reducing MIN_MATCHES threshold\")\n    print(\"  - Increasing TOP_K_SIMILAR\")\n    print(\"  - Using more/different images\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:08:39.807283Z","iopub.execute_input":"2025-12-08T16:08:39.807972Z","iopub.status.idle":"2025-12-08T16:09:13.705739Z","shell.execute_reply.started":"2025-12-08T16:08:39.807949Z","shell.execute_reply":"2025-12-08T16:09:13.705066Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Prepare Data for COLMAP","metadata":{}},{"cell_type":"code","source":"print(\"Preparing data for 3D reconstruction...\")\n\n# Create COLMAP workspace\ncolmap_workspace = Path(config.OUTPUT_PATH) / \"colmap_workspace\"\ncolmap_workspace.mkdir(exist_ok=True)\n\nimages_dir = colmap_workspace / \"images\"\nimages_dir.mkdir(exist_ok=True)\n\n# Copy/save images to workspace\nfor idx, (img, name) in enumerate(zip(images, image_names)):\n    img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n    cv2.imwrite(str(images_dir / name), img_bgr)\n\n# Save matches to text file for COLMAP\nmatches_file = colmap_workspace / \"matches.txt\"\nwith open(matches_file, 'w') as f:\n    for (i, j), match_data in verified_matches.items():\n        f.write(f\"{image_names[i]} {image_names[j]}\\n\")\n        mkpts0 = match_data['mkpts0']\n        mkpts1 = match_data['mkpts1']\n        for pt0, pt1 in zip(mkpts0, mkpts1):\n            f.write(f\"{pt0[0]:.2f} {pt0[1]:.2f} {pt1[0]:.2f} {pt1[1]:.2f}\\n\")\n\nprint(f\"✓ Prepared {len(images)} images and {len(verified_matches)} match pairs\")\nprint(f\"  Workspace: {colmap_workspace}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:10:22.704823Z","iopub.execute_input":"2025-12-08T16:10:22.705101Z","iopub.status.idle":"2025-12-08T16:10:26.391780Z","shell.execute_reply.started":"2025-12-08T16:10:22.705080Z","shell.execute_reply":"2025-12-08T16:10:26.391051Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Run COLMAP Reconstruction","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"RUNNING COLMAP 3D RECONSTRUCTION\")\nprint(\"=\"*70 + \"\\n\")\n\n# Check if COLMAP is available\ntry:\n    import subprocess\n    result = subprocess.run(['colmap', '--help'], capture_output=True, timeout=5)\n    colmap_available = result.returncode == 0\n    print(\"✓ COLMAP is available\")\nexcept:\n    colmap_available = False\n    print(\"⚠ COLMAP not available, will use simple triangulation\")\n\ndef run_colmap_reconstruction(scene_name, scene_data, output_base):\n    \"\"\"Run COLMAP reconstruction for a scene\"\"\"\n    \n    # Create workspace\n    workspace = Path(output_base) / scene_name\n    workspace.mkdir(parents=True, exist_ok=True)\n    \n    images_dir = workspace / \"images\"\n    database_path = workspace / \"database.db\"\n    sparse_dir = workspace / \"sparse\"\n    sparse_dir.mkdir(exist_ok=True)\n    \n    images_dir.mkdir(exist_ok=True)\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    verified_matches = scene_data['verified_matches']\n    \n    print(f\"  Scene: {scene_name}\")\n    print(f\"    Images: {len(images)}\")\n    print(f\"    Verified pairs: {len(verified_matches)}\")\n    \n    # Save images\n    for img, name in zip(images, image_names):\n        img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n        cv2.imwrite(str(images_dir / name), img_bgr)\n    \n    # Write matches in COLMAP format\n    matches_import_dir = workspace / \"matches_import\"\n    matches_import_dir.mkdir(exist_ok=True)\n    \n    # Write image pairs\n    pairs_file = matches_import_dir / \"image_pairs.txt\"\n    with open(pairs_file, 'w') as f:\n        for (i, j) in verified_matches.keys():\n            f.write(f\"{image_names[i]} {image_names[j]}\\n\")\n    \n    # Write matches\n    matches_file = matches_import_dir / \"matches.txt\"\n    with open(matches_file, 'w') as f:\n        for (i, j), match_data in verified_matches.items():\n            f.write(f\"{image_names[i]} {image_names[j]}\\n\")\n            mkpts0 = match_data['mkpts0']\n            mkpts1 = match_data['mkpts1']\n            for pt0, pt1 in zip(mkpts0, mkpts1):\n                f.write(f\"{pt0[0]} {pt0[1]} {pt1[0]} {pt1[1]}\\n\")\n    \n    if not colmap_available:\n        print(f\"    ⚠ COLMAP not available, skipping reconstruction\")\n        return None\n    \n    try:\n        # Feature extraction\n        print(f\"    Running feature extraction...\")\n        subprocess.run([\n            'colmap', 'feature_extractor',\n            '--database_path', str(database_path),\n            '--image_path', str(images_dir),\n            '--ImageReader.single_camera', '1',\n            '--ImageReader.camera_model', 'SIMPLE_RADIAL',\n            '--SiftExtraction.max_num_features', '8192'\n        ], check=True, capture_output=True, timeout=300)\n        \n        # Feature matching\n        print(f\"    Running feature matching...\")\n        subprocess.run([\n            'colmap', 'exhaustive_matcher',\n            '--database_path', str(database_path),\n            '--SiftMatching.guided_matching', '1'\n        ], check=True, capture_output=True, timeout=300)\n        \n        # Mapper (reconstruction)\n        print(f\"    Running mapper...\")\n        subprocess.run([\n            'colmap', 'mapper',\n            '--database_path', str(database_path),\n            '--image_path', str(images_dir),\n            '--output_path', str(sparse_dir),\n            '--Mapper.ba_refine_focal_length', '1',\n            '--Mapper.ba_refine_extra_params', '1',\n            '--Mapper.min_num_matches', str(config.MIN_MATCHES)\n        ], check=True, capture_output=True, timeout=600)\n        \n        print(f\"    ✓ COLMAP reconstruction complete\")\n        \n        # Find the reconstruction (usually in sparse/0)\n        recon_dir = sparse_dir / \"0\"\n        if recon_dir.exists():\n            return recon_dir\n        else:\n            print(f\"    ⚠ No reconstruction found\")\n            return None\n            \n    except subprocess.TimeoutExpired:\n        print(f\"    ✗ COLMAP timeout\")\n        return None\n    except subprocess.CalledProcessError as e:\n        print(f\"    ✗ COLMAP error: {e}\")\n        return None\n    except Exception as e:\n        print(f\"    ✗ Error: {e}\")\n        return None\n\n# Run COLMAP for each scene\ncolmap_results = {}\n\nif config.USE_COLMAP and colmap_available:\n    for scene_name, scene_data in all_results.items():\n        recon_dir = run_colmap_reconstruction(\n            scene_name, \n            scene_data, \n            config.OUTPUT_PATH + \"/colmap_scenes\"\n        )\n        colmap_results[scene_name] = recon_dir\n        print()\nelse:\n    print(\"Skipping COLMAP reconstruction (not available or disabled)\")\n\nprint(f\"\\n✓ Completed reconstructions for {len(colmap_results)} scenes\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:28:06.361527Z","iopub.execute_input":"2025-12-08T16:28:06.362250Z","iopub.status.idle":"2025-12-08T16:28:06.377273Z","shell.execute_reply.started":"2025-12-08T16:28:06.362226Z","shell.execute_reply":"2025-12-08T16:28:06.376626Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load and Visualize COLMAP Results","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"LOADING COLMAP RECONSTRUCTIONS\")\nprint(\"=\"*70 + \"\\n\")\n\ndef read_colmap_points3D(points3D_path):\n    \"\"\"Read COLMAP points3D.txt file\"\"\"\n    points = []\n    colors = []\n    \n    with open(points3D_path, 'r') as f:\n        for line in f:\n            if line.startswith('#'):\n                continue\n            parts = line.strip().split()\n            if len(parts) >= 7:\n                # Format: POINT3D_ID, X, Y, Z, R, G, B, ERROR, TRACK[]\n                x, y, z = float(parts[1]), float(parts[2]), float(parts[3])\n                r, g, b = int(parts[4]), int(parts[5]), int(parts[6])\n                points.append([x, y, z])\n                colors.append([r, g, b])\n    \n    return np.array(points), np.array(colors)\n\ndef read_colmap_binary_points3D(points3D_path):\n    \"\"\"Read COLMAP points3D.bin file\"\"\"\n    try:\n        import struct\n        points = []\n        colors = []\n        \n        with open(points3D_path, 'rb') as f:\n            num_points = struct.unpack('Q', f.read(8))[0]\n            for _ in range(num_points):\n                point3D_id = struct.unpack('Q', f.read(8))[0]\n                xyz = struct.unpack('ddd', f.read(24))\n                rgb = struct.unpack('BBB', f.read(3))\n                error = struct.unpack('d', f.read(8))[0]\n                track_length = struct.unpack('Q', f.read(8))[0]\n                f.read(8 * track_length)  # Skip track elements\n                \n                points.append(xyz)\n                colors.append(rgb)\n        \n        return np.array(points), np.array(colors)\n    except:\n        return np.array([]), np.array([])\n\n# Load all COLMAP reconstructions\nreconstructions = {}\n\nfor scene_name, recon_dir in colmap_results.items():\n    if recon_dir is None or not Path(recon_dir).exists():\n        continue\n    \n    print(f\"Loading reconstruction: {scene_name}\")\n    \n    # Try binary format first\n    points3D_bin = Path(recon_dir) / \"points3D.bin\"\n    points3D_txt = Path(recon_dir) / \"points3D.txt\"\n    \n    if points3D_bin.exists():\n        points, colors = read_colmap_binary_points3D(points3D_bin)\n    elif points3D_txt.exists():\n        points, colors = read_colmap_points3D(points3D_txt)\n    else:\n        print(f\"  ⚠ No points3D file found\")\n        continue\n    \n    if len(points) > 0:\n        # Remove outliers\n        mean = np.mean(points, axis=0)\n        std = np.std(points, axis=0)\n        mask = np.all(np.abs(points - mean) < 3 * std, axis=1)\n        points = points[mask]\n        colors = colors[mask]\n        \n        reconstructions[scene_name] = {\n            'points': points,\n            'colors': colors,\n            'recon_dir': recon_dir\n        }\n        print(f\"  ✓ Loaded {len(points)} 3D points\")\n    else:\n        print(f\"  ⚠ No points loaded\")\n\nprint(f\"\\n✓ Loaded {len(reconstructions)} reconstructions\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:29:56.026172Z","iopub.execute_input":"2025-12-08T16:29:56.027045Z","iopub.status.idle":"2025-12-08T16:29:56.040340Z","shell.execute_reply.started":"2025-12-08T16:29:56.026989Z","shell.execute_reply":"2025-12-08T16:29:56.039685Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualize All Reconstructions","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"CREATING VISUALIZATIONS\")\nprint(\"=\"*70 + \"\\n\")\n\nimport plotly.graph_objects as go\nfrom mpl_toolkits.mplot3d import Axes3D\n\n# Create visualizations for each reconstruction\nfor scene_name, recon_data in reconstructions.items():\n    print(f\"Visualizing: {scene_name}\")\n    \n    points = recon_data['points']\n    colors = recon_data['colors']\n    \n    if len(points) == 0:\n        continue\n    \n    # Create output directory for this scene\n    scene_output = Path(config.OUTPUT_PATH) / \"visualizations\" / scene_name\n    scene_output.mkdir(parents=True, exist_ok=True)\n    \n    # Plotly 3D interactive visualization\n    fig = go.Figure(data=[\n        go.Scatter3d(\n            x=points[:, 0],\n            y=points[:, 1],\n            z=points[:, 2],\n            mode='markers',\n            marker=dict(\n                size=1,\n                color=colors if len(colors) > 0 else 'blue',\n                opacity=0.8\n            ),\n            name='3D Points'\n        )\n    ])\n    \n    fig.update_layout(\n        title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n        scene=dict(\n            xaxis_title='X',\n            yaxis_title='Y',\n            zaxis_title='Z',\n            aspectmode='data'\n        ),\n        width=1000,\n        height=800\n    )\n    \n    html_path = scene_output / \"reconstruction_3d.html\"\n    fig.write_html(str(html_path))\n    print(f\"  Saved interactive: {html_path}\")\n    \n    # Matplotlib 3D plot\n    fig = plt.figure(figsize=(15, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    # Subsample for faster rendering\n    subsample = min(10000, len(points))\n    indices = np.random.choice(len(points), subsample, replace=False)\n    \n    ax.scatter(\n        points[indices, 0],\n        points[indices, 1],\n        points[indices, 2],\n        c=colors[indices] / 255.0 if len(colors) > 0 else 'blue',\n        s=1,\n        alpha=0.6\n    )\n    \n    ax.set_xlabel('X')\n    ax.set_ylabel('Y')\n    ax.set_zlabel('Z')\n    ax.set_title(f'3D Reconstruction - {scene_name}\\n{len(points)} points')\n    \n    png_path = scene_output / \"reconstruction_3d.png\"\n    plt.savefig(str(png_path), dpi=150, bbox_inches='tight')\n    plt.close()\n    print(f\"  Saved image: {png_path}\")\n    \n    # Show first reconstruction\n    if scene_name == list(reconstructions.keys())[0]:\n        fig = go.Figure(data=[\n            go.Scatter3d(\n                x=points[:, 0],\n                y=points[:, 1],\n                z=points[:, 2],\n                mode='markers',\n                marker=dict(size=1, color=colors if len(colors) > 0 else 'blue', opacity=0.8),\n                name='3D Points'\n            )\n        ])\n        fig.update_layout(\n            title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n            scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z', aspectmode='data'),\n            width=1000, height=800\n        )\n        fig.show()\n    \n    print()\n\nprint(f\"✓ Created visualizations for {len(reconstructions)} scenes\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:30:35.258093Z","iopub.execute_input":"2025-12-08T16:30:35.258747Z","iopub.status.idle":"2025-12-08T16:30:35.322099Z","shell.execute_reply.started":"2025-12-08T16:30:35.258708Z","shell.execute_reply":"2025-12-08T16:30:35.321480Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Summary and Statistics","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"FINAL PIPELINE SUMMARY\")\nprint(\"=\"*70 + \"\\n\")\n\ntotal_images = sum(len(data['images']) for data in scenes_data.values())\ntotal_matches = sum(len(data['verified_matches']) for data in all_results.values())\ntotal_points = sum(len(data['points']) for data in reconstructions.values())\n\nprint(f\"{'Metric':<40} {'Value':>20}\")\nprint(\"-\" * 62)\nprint(f\"{'Total scenes processed':<40} {len(scenes_data):>20}\")\nprint(f\"{'Total images loaded':<40} {total_images:>20}\")\nprint(f\"{'Total verified match pairs':<40} {total_matches:>20}\")\nprint(f\"{'Successful reconstructions':<40} {len(reconstructions):>20}\")\nprint(f\"{'Total 3D points reconstructed':<40} {total_points:>20,}\")\nprint()\n\nprint(\"Per-Scene Statistics:\")\nprint(\"-\" * 62)\nprint(f\"{'Scene Name':<30} {'Images':>10} {'Matches':>10} {'3D Points':>10}\")\nprint(\"-\" * 62)\n\nfor scene_name in scenes_data.keys():\n    num_images = len(scenes_data[scene_name]['images'])\n    num_matches = len(all_results[scene_name]['verified_matches'])\n    num_points = len(reconstructions[scene_name]['points']) if scene_name in reconstructions else 0\n    print(f\"{scene_name:<30} {num_images:>10} {num_matches:>10} {num_points:>10,}\")\n\nprint()\nprint(\"Output Locations:\")\nprint(f\"  Main output: {config.OUTPUT_PATH}\")\nprint(f\"  COLMAP workspaces: {config.OUTPUT_PATH}/colmap_scenes/\")\nprint(f\"  Visualizations: {config.OUTPUT_PATH}/visualizations/\")\nprint()\n\n# Create summary visualization\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Plot 1: Images per scene\nscene_names = list(scenes_data.keys())\nimage_counts = [len(scenes_data[s]['images']) for s in scene_names]\naxes[0, 0].barh(range(len(scene_names)), image_counts, color='steelblue')\naxes[0, 0].set_yticks(range(len(scene_names)))\naxes[0, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 0].set_xlabel('Number of Images')\naxes[0, 0].set_title('Images per Scene')\naxes[0, 0].grid(axis='x', alpha=0.3)\n\n# Plot 2: Verified matches per scene\nmatch_counts = [len(all_results[s]['verified_matches']) for s in scene_names]\naxes[0, 1].barh(range(len(scene_names)), match_counts, color='forestgreen')\naxes[0, 1].set_yticks(range(len(scene_names)))\naxes[0, 1].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 1].set_xlabel('Number of Verified Matches')\naxes[0, 1].set_title('Verified Matches per Scene')\naxes[0, 1].grid(axis='x', alpha=0.3)\n\n# Plot 3: 3D points per scene\npoint_counts = [len(reconstructions[s]['points']) if s in reconstructions else 0 for s in scene_names]\naxes[1, 0].barh(range(len(scene_names)), point_counts, color='coral')\naxes[1, 0].set_yticks(range(len(scene_names)))\naxes[1, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[1, 0].set_xlabel('Number of 3D Points')\naxes[1, 0].set_title('3D Points per Scene')\naxes[1, 0].grid(axis='x', alpha=0.3)\n\n# Plot 4: Pipeline overview\npipeline_stages = ['Images\\nLoaded', 'Pairs\\nCreated', 'Pairs\\nMatched', 'Pairs\\nVerified', '3D Points\\nReconstructed']\npipeline_counts = [\n    total_images,\n    sum(len(all_results[s]['verified_matches']) * 2 for s in scene_names),  # Approximate pairs created\n    total_matches,\n    total_matches,\n    total_points\n]\naxes[1, 1].plot(pipeline_stages, pipeline_counts, marker='o', linewidth=2, markersize=10, color='darkviolet')\naxes[1, 1].set_ylabel('Count')\naxes[1, 1].set_title('Pipeline Flow')\naxes[1, 1].grid(alpha=0.3)\naxes[1, 1].tick_params(axis='x', rotation=0)\n\nplt.tight_layout()\nplt.savefig(f\"{config.OUTPUT_PATH}/pipeline_summary.png\", dpi=150, bbox_inches='tight')\nplt.show()\n\nprint(\"✓ Pipeline complete!\")\nprint(f\"\\n{'='*70}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:31:12.032627Z","iopub.execute_input":"2025-12-08T16:31:12.032913Z","iopub.status.idle":"2025-12-08T16:31:13.777332Z","shell.execute_reply.started":"2025-12-08T16:31:12.032894Z","shell.execute_reply":"2025-12-08T16:31:13.776661Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Run COLMAP/pycolmap Reconstruction","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"RUNNING 3D RECONSTRUCTION\")\nprint(\"=\"*70 + \"\\n\")\n\n# Check if pycolmap is available\ntry:\n    import pycolmap\n    print(\"✓ pycolmap is available\")\n    use_pycolmap = True\nexcept ImportError:\n    print(\"⚠ pycolmap not available, will use simple triangulation\")\n    use_pycolmap = False\n\n\ndef reconstruct_with_pycolmap(scene_name, scene_data, output_base):\n    \"\"\"Reconstruct using pycolmap (Python-based)\"\"\"\n    import pycolmap\n    \n    workspace = Path(output_base) / scene_name\n    workspace.mkdir(parents=True, exist_ok=True)\n    \n    images_dir = workspace / \"images\"\n    images_dir.mkdir(exist_ok=True)\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    verified_matches = scene_data['verified_matches']\n    \n    print(f\"  Scene: {scene_name}\")\n    print(f\"    Images: {len(images)}\")\n    print(f\"    Verified pairs: {len(verified_matches)}\")\n    \n    # Save images to disk\n    for img, name in zip(images, image_names):\n        img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n        cv2.imwrite(str(images_dir / name), img_bgr)\n    \n    # Create pycolmap reconstruction\n    output_path = workspace / \"sparse\"\n    output_path.mkdir(exist_ok=True)\n    \n    try:\n        # Setup reconstruction\n        reconstruction = pycolmap.Reconstruction()\n        \n        # Add camera (assume all images use same camera with approximate parameters)\n        camera = pycolmap.Camera(\n            model='SIMPLE_RADIAL',\n            width=images[0].shape[1],\n            height=images[0].shape[0],\n            params=[max(images[0].shape[:2]), images[0].shape[1]/2, images[0].shape[0]/2, 0]\n        )\n        camera_id = reconstruction.add_camera(camera)\n        \n        # Add images\n        image_id_map = {}\n        for idx, name in enumerate(image_names):\n            image = pycolmap.Image(\n                id=idx + 1,\n                name=name,\n                camera_id=camera_id\n            )\n            image_id = reconstruction.add_image(image)\n            image_id_map[idx] = image_id\n        \n        # Add keypoints and matches\n        for (i, j), match_data in verified_matches.items():\n            img_id1 = image_id_map[i]\n            img_id2 = image_id_map[j]\n            \n            # Add keypoints to images\n            mkpts0 = match_data['mkpts0']\n            mkpts1 = match_data['mkpts1']\n            \n            # This is simplified - in production you'd need proper feature management\n            # For now, we'll use a simpler triangulation approach\n        \n        print(f\"    ⚠ pycolmap setup complete but needs full implementation\")\n        print(f\"    Falling back to simple triangulation...\")\n        return None\n        \n    except Exception as e:\n        print(f\"    ✗ pycolmap error: {e}\")\n        return None\n\ndef reconstruct_simple_triangulation(scene_name, scene_data, output_base):\n    \"\"\"Simple reconstruction using triangulation (fallback)\"\"\"\n    \n    workspace = Path(output_base) / scene_name\n    workspace.mkdir(parents=True, exist_ok=True)\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    verified_matches = scene_data['verified_matches']\n    \n    print(f\"  Scene: {scene_name}\")\n    print(f\"    Images: {len(images)}\")\n    print(f\"    Verified pairs: {len(verified_matches)}\")\n    \n    if len(verified_matches) == 0:\n        print(f\"    ⚠ No matches to reconstruct\")\n        return None\n    \n    # Simple triangulation approach\n    all_points_3d = []\n    all_colors = []\n    \n    # Process multiple pairs\n    for pair_idx, ((i, j), match_data) in enumerate(list(verified_matches.items())[:10]):\n        try:\n            mkpts0 = match_data['mkpts0']\n            mkpts1 = match_data['mkpts1']\n            \n            if len(mkpts0) < 8:\n                continue\n            \n            # Estimate essential matrix\n            E, mask = cv2.findEssentialMat(\n                mkpts0, mkpts1,\n                focal=max(images[i].shape[:2]),\n                pp=(images[i].shape[1]/2, images[i].shape[0]/2),\n                method=cv2.RANSAC,\n                prob=0.999,\n                threshold=1.0\n            )\n            \n            if E is None:\n                continue\n            \n            # Recover pose\n            _, R, t, mask_pose = cv2.recoverPose(E, mkpts0, mkpts1)\n            \n            # Create projection matrices\n            P1 = np.hstack([np.eye(3), np.zeros((3, 1))])\n            P2 = np.hstack([R, t])\n            \n            # Triangulate points\n            points_4d = cv2.triangulatePoints(P1, P2, mkpts0.T, mkpts1.T)\n            points_3d = (points_4d[:3] / points_4d[3]).T\n            \n            # Get colors from first image\n            colors = []\n            for pt in mkpts0.astype(int):\n                y, x = np.clip(pt[1], 0, images[i].shape[0]-1), np.clip(pt[0], 0, images[i].shape[1]-1)\n                colors.append(images[i][y, x])\n            \n            all_points_3d.append(points_3d)\n            all_colors.append(np.array(colors))\n            \n        except Exception as e:\n            continue\n    \n    if len(all_points_3d) == 0:\n        print(f\"    ⚠ No 3D points reconstructed\")\n        return None\n    \n    # Combine all points\n    points_3d = np.vstack(all_points_3d)\n    colors = np.vstack(all_colors)\n    \n    # Remove outliers (points too far from median)\n    median = np.median(points_3d, axis=0)\n    distances = np.linalg.norm(points_3d - median, axis=1)\n    threshold = np.percentile(distances, 95)  # Keep 95% of points\n    mask = distances < threshold\n    \n    points_3d = points_3d[mask]\n    colors = colors[mask]\n    \n    print(f\"    ✓ Reconstructed {len(points_3d)} 3D points\")\n    \n    # Save as simple format\n    output_file = workspace / \"points3D.txt\"\n    with open(output_file, 'w') as f:\n        f.write(\"# 3D point list with RGB colors\\n\")\n        f.write(\"# Format: X Y Z R G B\\n\")\n        for pt, col in zip(points_3d, colors):\n            f.write(f\"{pt[0]:.6f} {pt[1]:.6f} {pt[2]:.6f} {int(col[0])} {int(col[1])} {int(col[2])}\\n\")\n    \n    return {\n        'points': points_3d,\n        'colors': colors,\n        'output_file': output_file,\n        'workspace': workspace\n    }\n\n# Run reconstruction for each scene\nreconstructions = {}\n\nfor scene_name, scene_data in all_results.items():\n    print(f\"\\n{'─'*60}\")\n    \n    # Try pycolmap first, fallback to simple triangulation\n    if use_pycolmap and config.USE_COLMAP:\n        recon = reconstruct_with_pycolmap(\n            scene_name, \n            scene_data, \n            config.OUTPUT_PATH + \"/reconstructions\"\n        )\n    else:\n        recon = None\n    \n    # Fallback to simple triangulation\n    if recon is None:\n        recon = reconstruct_simple_triangulation(\n            scene_name,\n            scene_data,\n            config.OUTPUT_PATH + \"/reconstructions\"\n        )\n    \n    if recon is not None:\n        reconstructions[scene_name] = recon\n\nprint(f\"\\n{'='*70}\")\nprint(f\"✓ Completed {len(reconstructions)}/{len(all_results)} reconstructions\")\nprint(f\"{'='*70}\\n\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:44:38.284813Z","iopub.execute_input":"2025-12-08T16:44:38.285457Z","iopub.status.idle":"2025-12-08T16:44:39.266318Z","shell.execute_reply.started":"2025-12-08T16:44:38.285415Z","shell.execute_reply":"2025-12-08T16:44:39.265634Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load Reconstructions","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"LOADING RECONSTRUCTIONS\")\nprint(\"=\"*70 + \"\\n\")\n\n# Reconstructions are already loaded from Cell 15!\n# They're in the 'reconstructions' dictionary\n\nif len(reconstructions) > 0:\n    print(f\"✓ Loaded {len(reconstructions)} reconstructions\\n\")\n    \n    for scene_name, recon_data in reconstructions.items():\n        num_points = len(recon_data['points'])\n        print(f\"  {scene_name}: {num_points:,} 3D points\")\nelse:\n    print(\"⚠ No reconstructions available\")\n    print(\"  This might be because:\")\n    print(\"  - No verified matches were found\")\n    print(\"  - Triangulation failed\")\n    print(\"  - Try adjusting MIN_MATCHES or TOP_K_SIMILAR in Cell 2\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:45:40.066221Z","iopub.execute_input":"2025-12-08T16:45:40.066785Z","iopub.status.idle":"2025-12-08T16:45:40.073452Z","shell.execute_reply.started":"2025-12-08T16:45:40.066758Z","shell.execute_reply":"2025-12-08T16:45:40.072468Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualize All Reconstructions","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"CREATING VISUALIZATIONS\")\nprint(\"=\"*70 + \"\\n\")\n\nimport plotly.graph_objects as go\nfrom mpl_toolkits.mplot3d import Axes3D\n\n# Create visualizations for each reconstruction\nfor scene_name, recon_data in reconstructions.items():\n    print(f\"Visualizing: {scene_name}\")\n    \n    points = recon_data['points']\n    colors = recon_data['colors']\n    \n    if len(points) == 0:\n        continue\n    \n    # Create output directory for this scene\n    scene_output = Path(config.OUTPUT_PATH) / \"visualizations\" / scene_name\n    scene_output.mkdir(parents=True, exist_ok=True)\n    \n    # Plotly 3D interactive visualization\n    fig = go.Figure(data=[\n        go.Scatter3d(\n            x=points[:, 0],\n            y=points[:, 1],\n            z=points[:, 2],\n            mode='markers',\n            marker=dict(\n                size=1,\n                color=colors if len(colors) > 0 else 'blue',\n                opacity=0.8\n            ),\n            name='3D Points'\n        )\n    ])\n    \n    fig.update_layout(\n        title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n        scene=dict(\n            xaxis_title='X',\n            yaxis_title='Y',\n            zaxis_title='Z',\n            aspectmode='data'\n        ),\n        width=1000,\n        height=800\n    )\n    \n    html_path = scene_output / \"reconstruction_3d.html\"\n    fig.write_html(str(html_path))\n    print(f\"  Saved interactive: {html_path}\")\n    \n    # Matplotlib 3D plot\n    fig = plt.figure(figsize=(15, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    # Subsample for faster rendering\n    subsample = min(10000, len(points))\n    indices = np.random.choice(len(points), subsample, replace=False)\n    \n    ax.scatter(\n        points[indices, 0],\n        points[indices, 1],\n        points[indices, 2],\n        c=colors[indices] / 255.0 if len(colors) > 0 else 'blue',\n        s=1,\n        alpha=0.6\n    )\n    \n    ax.set_xlabel('X')\n    ax.set_ylabel('Y')\n    ax.set_zlabel('Z')\n    ax.set_title(f'3D Reconstruction - {scene_name}\\n{len(points)} points')\n    \n    png_path = scene_output / \"reconstruction_3d.png\"\n    plt.savefig(str(png_path), dpi=150, bbox_inches='tight')\n    plt.close()\n    print(f\"  Saved image: {png_path}\")\n    \n    # Show first reconstruction\n    if scene_name == list(reconstructions.keys())[0]:\n        fig = go.Figure(data=[\n            go.Scatter3d(\n                x=points[:, 0],\n                y=points[:, 1],\n                z=points[:, 2],\n                mode='markers',\n                marker=dict(size=1, color=colors if len(colors) > 0 else 'blue', opacity=0.8),\n                name='3D Points'\n            )\n        ])\n        fig.update_layout(\n            title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n            scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z', aspectmode='data'),\n            width=1000, height=800\n        )\n        fig.show()\n    \n    print()\n\nprint(f\"✓ Created visualizations for {len(reconstructions)} scenes\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:47:01.941688Z","iopub.execute_input":"2025-12-08T16:47:01.941962Z","iopub.status.idle":"2025-12-08T16:47:08.374419Z","shell.execute_reply.started":"2025-12-08T16:47:01.941945Z","shell.execute_reply":"2025-12-08T16:47:08.373695Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Summary and Statistics","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"FINAL PIPELINE SUMMARY\")\nprint(\"=\"*70 + \"\\n\")\n\ntotal_images = sum(len(data['images']) for data in scenes_data.values())\ntotal_matches = sum(len(data['verified_matches']) for data in all_results.values())\ntotal_points = sum(len(data['points']) for data in reconstructions.values())\n\nprint(f\"{'Metric':<40} {'Value':>20}\")\nprint(\"-\" * 62)\nprint(f\"{'Total scenes processed':<40} {len(scenes_data):>20}\")\nprint(f\"{'Total images loaded':<40} {total_images:>20}\")\nprint(f\"{'Total verified match pairs':<40} {total_matches:>20}\")\nprint(f\"{'Successful reconstructions':<40} {len(reconstructions):>20}\")\nprint(f\"{'Total 3D points reconstructed':<40} {total_points:>20,}\")\nprint()\n\nprint(\"Per-Scene Statistics:\")\nprint(\"-\" * 62)\nprint(f\"{'Scene Name':<30} {'Images':>10} {'Matches':>10} {'3D Points':>10}\")\nprint(\"-\" * 62)\n\nfor scene_name in scenes_data.keys():\n    num_images = len(scenes_data[scene_name]['images'])\n    num_matches = len(all_results[scene_name]['verified_matches'])\n    num_points = len(reconstructions[scene_name]['points']) if scene_name in reconstructions else 0\n    print(f\"{scene_name:<30} {num_images:>10} {num_matches:>10} {num_points:>10,}\")\n\nprint()\nprint(\"Output Locations:\")\nprint(f\"  Main output: {config.OUTPUT_PATH}\")\nprint(f\"  COLMAP workspaces: {config.OUTPUT_PATH}/colmap_scenes/\")\nprint(f\"  Visualizations: {config.OUTPUT_PATH}/visualizations/\")\nprint()\n\n# Create summary visualization\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Plot 1: Images per scene\nscene_names = list(scenes_data.keys())\nimage_counts = [len(scenes_data[s]['images']) for s in scene_names]\naxes[0, 0].barh(range(len(scene_names)), image_counts, color='steelblue')\naxes[0, 0].set_yticks(range(len(scene_names)))\naxes[0, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 0].set_xlabel('Number of Images')\naxes[0, 0].set_title('Images per Scene')\naxes[0, 0].grid(axis='x', alpha=0.3)\n\n# Plot 2: Verified matches per scene\nmatch_counts = [len(all_results[s]['verified_matches']) for s in scene_names]\naxes[0, 1].barh(range(len(scene_names)), match_counts, color='forestgreen')\naxes[0, 1].set_yticks(range(len(scene_names)))\naxes[0, 1].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 1].set_xlabel('Number of Verified Matches')\naxes[0, 1].set_title('Verified Matches per Scene')\naxes[0, 1].grid(axis='x', alpha=0.3)\n\n# Plot 3: 3D points per scene\npoint_counts = [len(reconstructions[s]['points']) if s in reconstructions else 0 for s in scene_names]\naxes[1, 0].barh(range(len(scene_names)), point_counts, color='coral')\naxes[1, 0].set_yticks(range(len(scene_names)))\naxes[1, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[1, 0].set_xlabel('Number of 3D Points')\naxes[1, 0].set_title('3D Points per Scene')\naxes[1, 0].grid(axis='x', alpha=0.3)\n\n# Plot 4: Pipeline overview\npipeline_stages = ['Images\\nLoaded', 'Pairs\\nCreated', 'Pairs\\nMatched', 'Pairs\\nVerified', '3D Points\\nReconstructed']\npipeline_counts = [\n    total_images,\n    sum(len(all_results[s]['verified_matches']) * 2 for s in scene_names),  # Approximate pairs created\n    total_matches,\n    total_matches,\n    total_points\n]\naxes[1, 1].plot(pipeline_stages, pipeline_counts, marker='o', linewidth=2, markersize=10, color='darkviolet')\naxes[1, 1].set_ylabel('Count')\naxes[1, 1].set_title('Pipeline Flow')\naxes[1, 1].grid(alpha=0.3)\naxes[1, 1].tick_params(axis='x', rotation=0)\n\nplt.tight_layout()\nplt.savefig(f\"{config.OUTPUT_PATH}/pipeline_summary.png\", dpi=150, bbox_inches='tight')\nplt.show()\n\nprint(\"✓ Pipeline complete!\")\nprint(f\"\\n{'='*70}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:52:12.576150Z","iopub.execute_input":"2025-12-08T16:52:12.576743Z","iopub.status.idle":"2025-12-08T16:52:14.196488Z","shell.execute_reply.started":"2025-12-08T16:52:12.576718Z","shell.execute_reply":"2025-12-08T16:52:14.195781Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# PyColMap","metadata":{}},{"cell_type":"code","source":"! pip install pycolmap","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:55:23.332684Z","iopub.execute_input":"2025-12-08T16:55:23.333514Z","iopub.status.idle":"2025-12-08T16:55:27.839789Z","shell.execute_reply.started":"2025-12-08T16:55:23.333487Z","shell.execute_reply":"2025-12-08T16:55:27.838910Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Prepare Data for COLMAP","metadata":{}},{"cell_type":"code","source":"print(\"Preparing data for 3D reconstruction...\")\n\n# Create COLMAP workspace\ncolmap_workspace = Path(config.OUTPUT_PATH) / \"colmap_workspace\"\ncolmap_workspace.mkdir(exist_ok=True)\n\nimages_dir = colmap_workspace / \"images\"\nimages_dir.mkdir(exist_ok=True)\n\n# Copy/save images to workspace\nfor idx, (img, name) in enumerate(zip(images, image_names)):\n    img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n    cv2.imwrite(str(images_dir / name), img_bgr)\n\n# Save matches to text file for COLMAP\nmatches_file = colmap_workspace / \"matches.txt\"\nwith open(matches_file, 'w') as f:\n    for (i, j), match_data in verified_matches.items():\n        f.write(f\"{image_names[i]} {image_names[j]}\\n\")\n        mkpts0 = match_data['mkpts0']\n        mkpts1 = match_data['mkpts1']\n        for pt0, pt1 in zip(mkpts0, mkpts1):\n            f.write(f\"{pt0[0]:.2f} {pt0[1]:.2f} {pt1[0]:.2f} {pt1[1]:.2f}\\n\")\n\nprint(f\"✓ Prepared {len(images)} images and {len(verified_matches)} match pairs\")\nprint(f\"  Workspace: {colmap_workspace}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:56:12.170562Z","iopub.execute_input":"2025-12-08T16:56:12.171267Z","iopub.status.idle":"2025-12-08T16:56:15.870070Z","shell.execute_reply.started":"2025-12-08T16:56:12.171236Z","shell.execute_reply":"2025-12-08T16:56:15.869292Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Run COLMAP/pycolmap Reconstruction","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"RUNNING 3D RECONSTRUCTION\")\nprint(\"=\"*70 + \"\\n\")\n\n# Check if pycolmap is available\ntry:\n    import pycolmap\n    print(\"✓ pycolmap is available\")\n    use_pycolmap = True\nexcept ImportError:\n    print(\"⚠ pycolmap not available, will use simple triangulation\")\n    use_pycolmap = False\n\n\ndef reconstruct_with_pycolmap(scene_name, scene_data, output_base):\n    \"\"\"Reconstruct using pycolmap (Python-based)\"\"\"\n    import pycolmap\n    \n    workspace = Path(output_base) / scene_name\n    workspace.mkdir(parents=True, exist_ok=True)\n    \n    images_dir = workspace / \"images\"\n    images_dir.mkdir(exist_ok=True)\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    verified_matches = scene_data['verified_matches']\n    \n    print(f\"  Scene: {scene_name}\")\n    print(f\"    Images: {len(images)}\")\n    print(f\"    Verified pairs: {len(verified_matches)}\")\n    \n    # Save images to disk\n    for img, name in zip(images, image_names):\n        img_bgr = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n        cv2.imwrite(str(images_dir / name), img_bgr)\n    \n    # Create pycolmap reconstruction\n    output_path = workspace / \"sparse\"\n    output_path.mkdir(exist_ok=True)\n    \n    try:\n        # Setup reconstruction\n        reconstruction = pycolmap.Reconstruction()\n        \n        # Add camera (assume all images use same camera with approximate parameters)\n        camera = pycolmap.Camera(\n            model='SIMPLE_RADIAL',\n            width=images[0].shape[1],\n            height=images[0].shape[0],\n            params=[max(images[0].shape[:2]), images[0].shape[1]/2, images[0].shape[0]/2, 0]\n        )\n        camera_id = reconstruction.add_camera(camera)\n        \n        # Add images\n        image_id_map = {}\n        for idx, name in enumerate(image_names):\n            image = pycolmap.Image(\n                id=idx + 1,\n                name=name,\n                camera_id=camera_id\n            )\n            image_id = reconstruction.add_image(image)\n            image_id_map[idx] = image_id\n        \n        # Add keypoints and matches\n        for (i, j), match_data in verified_matches.items():\n            img_id1 = image_id_map[i]\n            img_id2 = image_id_map[j]\n            \n            # Add keypoints to images\n            mkpts0 = match_data['mkpts0']\n            mkpts1 = match_data['mkpts1']\n            \n            # This is simplified - in production you'd need proper feature management\n            # For now, we'll use a simpler triangulation approach\n        \n        print(f\"    ⚠ pycolmap setup complete but needs full implementation\")\n        print(f\"    Falling back to simple triangulation...\")\n        return None\n        \n    except Exception as e:\n        print(f\"    ✗ pycolmap error: {e}\")\n        return None\n\ndef reconstruct_simple_triangulation(scene_name, scene_data, output_base):\n    \"\"\"Simple reconstruction using triangulation (fallback)\"\"\"\n    \n    workspace = Path(output_base) / scene_name\n    workspace.mkdir(parents=True, exist_ok=True)\n    \n    images = scene_data['images']\n    image_names = scene_data['image_names']\n    verified_matches = scene_data['verified_matches']\n    \n    print(f\"  Scene: {scene_name}\")\n    print(f\"    Images: {len(images)}\")\n    print(f\"    Verified pairs: {len(verified_matches)}\")\n    \n    if len(verified_matches) == 0:\n        print(f\"    ⚠ No matches to reconstruct\")\n        return None\n    \n    # Simple triangulation approach\n    all_points_3d = []\n    all_colors = []\n    \n    # Process multiple pairs\n    for pair_idx, ((i, j), match_data) in enumerate(list(verified_matches.items())[:10]):\n        try:\n            mkpts0 = match_data['mkpts0']\n            mkpts1 = match_data['mkpts1']\n            \n            if len(mkpts0) < 8:\n                continue\n            \n            # Estimate essential matrix\n            E, mask = cv2.findEssentialMat(\n                mkpts0, mkpts1,\n                focal=max(images[i].shape[:2]),\n                pp=(images[i].shape[1]/2, images[i].shape[0]/2),\n                method=cv2.RANSAC,\n                prob=0.999,\n                threshold=1.0\n            )\n            \n            if E is None:\n                continue\n            \n            # Recover pose\n            _, R, t, mask_pose = cv2.recoverPose(E, mkpts0, mkpts1)\n            \n            # Create projection matrices\n            P1 = np.hstack([np.eye(3), np.zeros((3, 1))])\n            P2 = np.hstack([R, t])\n            \n            # Triangulate points\n            points_4d = cv2.triangulatePoints(P1, P2, mkpts0.T, mkpts1.T)\n            points_3d = (points_4d[:3] / points_4d[3]).T\n            \n            # Get colors from first image\n            colors = []\n            for pt in mkpts0.astype(int):\n                y, x = np.clip(pt[1], 0, images[i].shape[0]-1), np.clip(pt[0], 0, images[i].shape[1]-1)\n                colors.append(images[i][y, x])\n            \n            all_points_3d.append(points_3d)\n            all_colors.append(np.array(colors))\n            \n        except Exception as e:\n            continue\n    \n    if len(all_points_3d) == 0:\n        print(f\"    ⚠ No 3D points reconstructed\")\n        return None\n    \n    # Combine all points\n    points_3d = np.vstack(all_points_3d)\n    colors = np.vstack(all_colors)\n    \n    # Remove outliers (points too far from median)\n    median = np.median(points_3d, axis=0)\n    distances = np.linalg.norm(points_3d - median, axis=1)\n    threshold = np.percentile(distances, 95)  # Keep 95% of points\n    mask = distances < threshold\n    \n    points_3d = points_3d[mask]\n    colors = colors[mask]\n    \n    print(f\"    ✓ Reconstructed {len(points_3d)} 3D points\")\n    \n    # Save as simple format\n    output_file = workspace / \"points3D.txt\"\n    with open(output_file, 'w') as f:\n        f.write(\"# 3D point list with RGB colors\\n\")\n        f.write(\"# Format: X Y Z R G B\\n\")\n        for pt, col in zip(points_3d, colors):\n            f.write(f\"{pt[0]:.6f} {pt[1]:.6f} {pt[2]:.6f} {int(col[0])} {int(col[1])} {int(col[2])}\\n\")\n    \n    return {\n        'points': points_3d,\n        'colors': colors,\n        'output_file': output_file,\n        'workspace': workspace\n    }\n\n# Run reconstruction for each scene\nreconstructions = {}\n\nfor scene_name, scene_data in all_results.items():\n    print(f\"\\n{'─'*60}\")\n    \n    # Try pycolmap first, fallback to simple triangulation\n    if use_pycolmap and config.USE_COLMAP:\n        recon = reconstruct_with_pycolmap(\n            scene_name, \n            scene_data, \n            config.OUTPUT_PATH + \"/reconstructions\"\n        )\n    else:\n        recon = None\n    \n    # Fallback to simple triangulation\n    if recon is None:\n        recon = reconstruct_simple_triangulation(\n            scene_name,\n            scene_data,\n            config.OUTPUT_PATH + \"/reconstructions\"\n        )\n    \n    if recon is not None:\n        reconstructions[scene_name] = recon\n\nprint(f\"\\n{'='*70}\")\nprint(f\"✓ Completed {len(reconstructions)}/{len(all_results)} reconstructions\")\nprint(f\"{'='*70}\\n\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T16:58:54.884074Z","iopub.execute_input":"2025-12-08T16:58:54.884747Z","iopub.status.idle":"2025-12-08T17:00:14.842357Z","shell.execute_reply.started":"2025-12-08T16:58:54.884722Z","shell.execute_reply":"2025-12-08T17:00:14.841441Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load Reconstructions","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"LOADING RECONSTRUCTIONS\")\nprint(\"=\"*70 + \"\\n\")\n\n# Reconstructions are already loaded from Cell 15!\n# They're in the 'reconstructions' dictionary\n\nif len(reconstructions) > 0:\n    print(f\"✓ Loaded {len(reconstructions)} reconstructions\\n\")\n    \n    for scene_name, recon_data in reconstructions.items():\n        num_points = len(recon_data['points'])\n        print(f\"  {scene_name}: {num_points:,} 3D points\")\nelse:\n    print(\"⚠ No reconstructions available\")\n    print(\"  This might be because:\")\n    print(\"  - No verified matches were found\")\n    print(\"  - Triangulation failed\")\n    print(\"  - Try adjusting MIN_MATCHES or TOP_K_SIMILAR in Cell 2\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T17:02:58.619520Z","iopub.execute_input":"2025-12-08T17:02:58.620309Z","iopub.status.idle":"2025-12-08T17:02:58.625780Z","shell.execute_reply.started":"2025-12-08T17:02:58.620284Z","shell.execute_reply":"2025-12-08T17:02:58.624975Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Visualize All Reconstructions","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"CREATING VISUALIZATIONS\")\nprint(\"=\"*70 + \"\\n\")\n\nimport plotly.graph_objects as go\nfrom mpl_toolkits.mplot3d import Axes3D\n\n# Create visualizations for each reconstruction\nfor scene_name, recon_data in reconstructions.items():\n    print(f\"Visualizing: {scene_name}\")\n    \n    points = recon_data['points']\n    colors = recon_data['colors']\n    \n    if len(points) == 0:\n        continue\n    \n    # Create output directory for this scene\n    scene_output = Path(config.OUTPUT_PATH) / \"pycolmap_visualizations\" / scene_name\n    scene_output.mkdir(parents=True, exist_ok=True)\n    \n    # Plotly 3D interactive visualization\n    fig = go.Figure(data=[\n        go.Scatter3d(\n            x=points[:, 0],\n            y=points[:, 1],\n            z=points[:, 2],\n            mode='markers',\n            marker=dict(\n                size=1,\n                color=colors if len(colors) > 0 else 'blue',\n                opacity=0.8\n            ),\n            name='3D Points'\n        )\n    ])\n    \n    fig.update_layout(\n        title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n        scene=dict(\n            xaxis_title='X',\n            yaxis_title='Y',\n            zaxis_title='Z',\n            aspectmode='data'\n        ),\n        width=1000,\n        height=800\n    )\n    \n    html_path = scene_output / \"reconstruction_3d.html\"\n    fig.write_html(str(html_path))\n    print(f\"  Saved interactive: {html_path}\")\n    \n    # Matplotlib 3D plot\n    fig = plt.figure(figsize=(15, 10))\n    ax = fig.add_subplot(111, projection='3d')\n    \n    # Subsample for faster rendering\n    subsample = min(10000, len(points))\n    indices = np.random.choice(len(points), subsample, replace=False)\n    \n    ax.scatter(\n        points[indices, 0],\n        points[indices, 1],\n        points[indices, 2],\n        c=colors[indices] / 255.0 if len(colors) > 0 else 'blue',\n        s=1,\n        alpha=0.6\n    )\n    \n    ax.set_xlabel('X')\n    ax.set_ylabel('Y')\n    ax.set_zlabel('Z')\n    ax.set_title(f'3D Reconstruction - {scene_name}\\n{len(points)} points')\n    \n    png_path = scene_output / \"reconstruction_3d.png\"\n    plt.savefig(str(png_path), dpi=150, bbox_inches='tight')\n    plt.close()\n    print(f\"  Saved image: {png_path}\")\n    \n    # Show first reconstruction\n    if scene_name == list(reconstructions.keys())[0]:\n        fig = go.Figure(data=[\n            go.Scatter3d(\n                x=points[:, 0],\n                y=points[:, 1],\n                z=points[:, 2],\n                mode='markers',\n                marker=dict(size=1, color=colors if len(colors) > 0 else 'blue', opacity=0.8),\n                name='3D Points'\n            )\n        ])\n        fig.update_layout(\n            title=f'3D Reconstruction - {scene_name}<br>{len(points)} points',\n            scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z', aspectmode='data'),\n            width=1000, height=800\n        )\n        fig.show()\n    \n    print()\n\nprint(f\"✓ Created visualizations for {len(reconstructions)} scenes\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T17:03:56.473038Z","iopub.execute_input":"2025-12-08T17:03:56.473311Z","iopub.status.idle":"2025-12-08T17:04:01.102096Z","shell.execute_reply.started":"2025-12-08T17:03:56.473289Z","shell.execute_reply":"2025-12-08T17:04:01.101232Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Summary and Statistics","metadata":{}},{"cell_type":"code","source":"print(\"\\n\" + \"=\"*70)\nprint(\"FINAL PIPELINE SUMMARY\")\nprint(\"=\"*70 + \"\\n\")\n\ntotal_images = sum(len(data['images']) for data in scenes_data.values())\ntotal_matches = sum(len(data['verified_matches']) for data in all_results.values())\ntotal_points = sum(len(data['points']) for data in reconstructions.values())\n\nprint(f\"{'Metric':<40} {'Value':>20}\")\nprint(\"-\" * 62)\nprint(f\"{'Total scenes processed':<40} {len(scenes_data):>20}\")\nprint(f\"{'Total images loaded':<40} {total_images:>20}\")\nprint(f\"{'Total verified match pairs':<40} {total_matches:>20}\")\nprint(f\"{'Successful reconstructions':<40} {len(reconstructions):>20}\")\nprint(f\"{'Total 3D points reconstructed':<40} {total_points:>20,}\")\nprint()\n\nprint(\"Per-Scene Statistics:\")\nprint(\"-\" * 62)\nprint(f\"{'Scene Name':<30} {'Images':>10} {'Matches':>10} {'3D Points':>10}\")\nprint(\"-\" * 62)\n\nfor scene_name in scenes_data.keys():\n    num_images = len(scenes_data[scene_name]['images'])\n    num_matches = len(all_results[scene_name]['verified_matches'])\n    num_points = len(reconstructions[scene_name]['points']) if scene_name in reconstructions else 0\n    print(f\"{scene_name:<30} {num_images:>10} {num_matches:>10} {num_points:>10,}\")\n\nprint()\nprint(\"Output Locations:\")\nprint(f\"  Main output: {config.OUTPUT_PATH}\")\nprint(f\"  COLMAP workspaces: {config.OUTPUT_PATH}/colmap_scenes/\")\nprint(f\"  Visualizations: {config.OUTPUT_PATH}/visualizations/\")\nprint()\n\n# Create summary visualization\nfig, axes = plt.subplots(2, 2, figsize=(16, 12))\n\n# Plot 1: Images per scene\nscene_names = list(scenes_data.keys())\nimage_counts = [len(scenes_data[s]['images']) for s in scene_names]\naxes[0, 0].barh(range(len(scene_names)), image_counts, color='steelblue')\naxes[0, 0].set_yticks(range(len(scene_names)))\naxes[0, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 0].set_xlabel('Number of Images')\naxes[0, 0].set_title('Images per Scene')\naxes[0, 0].grid(axis='x', alpha=0.3)\n\n# Plot 2: Verified matches per scene\nmatch_counts = [len(all_results[s]['verified_matches']) for s in scene_names]\naxes[0, 1].barh(range(len(scene_names)), match_counts, color='forestgreen')\naxes[0, 1].set_yticks(range(len(scene_names)))\naxes[0, 1].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[0, 1].set_xlabel('Number of Verified Matches')\naxes[0, 1].set_title('Verified Matches per Scene')\naxes[0, 1].grid(axis='x', alpha=0.3)\n\n# Plot 3: 3D points per scene\npoint_counts = [len(reconstructions[s]['points']) if s in reconstructions else 0 for s in scene_names]\naxes[1, 0].barh(range(len(scene_names)), point_counts, color='coral')\naxes[1, 0].set_yticks(range(len(scene_names)))\naxes[1, 0].set_yticklabels([s[:20] for s in scene_names], fontsize=8)\naxes[1, 0].set_xlabel('Number of 3D Points')\naxes[1, 0].set_title('3D Points per Scene')\naxes[1, 0].grid(axis='x', alpha=0.3)\n\n# Plot 4: Pipeline overview\npipeline_stages = ['Images\\nLoaded', 'Pairs\\nCreated', 'Pairs\\nMatched', 'Pairs\\nVerified', '3D Points\\nReconstructed']\npipeline_counts = [\n    total_images,\n    sum(len(all_results[s]['verified_matches']) * 2 for s in scene_names),  # Approximate pairs created\n    total_matches,\n    total_matches,\n    total_points\n]\naxes[1, 1].plot(pipeline_stages, pipeline_counts, marker='o', linewidth=2, markersize=10, color='darkviolet')\naxes[1, 1].set_ylabel('Count')\naxes[1, 1].set_title('Pipeline Flow')\naxes[1, 1].grid(alpha=0.3)\naxes[1, 1].tick_params(axis='x', rotation=0)\n\nplt.tight_layout()\nplt.savefig(f\"{config.OUTPUT_PATH}/pipeline_summary.png\", dpi=150, bbox_inches='tight')\nplt.show()\n\nprint(\"✓ Pipeline complete!\")\nprint(f\"\\n{'='*70}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-08T17:05:42.067933Z","iopub.execute_input":"2025-12-08T17:05:42.068775Z","iopub.status.idle":"2025-12-08T17:05:43.710710Z","shell.execute_reply.started":"2025-12-08T17:05:42.068747Z","shell.execute_reply":"2025-12-08T17:05:43.709851Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}