{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","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"},{"sourceId":7884485,"sourceType":"datasetVersion","datasetId":4628051},{"sourceId":11924468,"sourceType":"datasetVersion","datasetId":6988459},{"sourceId":234271505,"sourceType":"kernelVersion"},{"sourceId":237741314,"sourceType":"kernelVersion"},{"sourceId":4534,"sourceType":"modelInstanceVersion","modelInstanceId":3326,"modelId":986},{"sourceId":4535,"sourceType":"modelInstanceVersion","modelInstanceId":3327,"modelId":986},{"sourceId":4537,"sourceType":"modelInstanceVersion","modelInstanceId":3329,"modelId":986},{"sourceId":17191,"sourceType":"modelInstanceVersion","modelInstanceId":14317,"modelId":21716},{"sourceId":17555,"sourceType":"modelInstanceVersion","modelInstanceId":14611,"modelId":22086},{"sourceId":17578,"sourceType":"modelInstanceVersion","modelInstanceId":14634,"modelId":22086},{"sourceId":17579,"sourceType":"modelInstanceVersion","modelInstanceId":14635,"modelId":22086},{"sourceId":677159,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":513463,"modelId":528097}],"dockerImageVersionId":30919,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Dependencies","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"# install colmap(CPU)\n!cd /kaggle/input/pkg-colmap/colmap_offline && dpkg -i ./*.deb\n\n# test\n!colmap -h","metadata":{"trusted":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2025-12-09T03:24:28.178096Z","iopub.execute_input":"2025-12-09T03:24:28.178376Z","iopub.status.idle":"2025-12-09T03:24:42.451635Z","shell.execute_reply.started":"2025-12-09T03:24:28.178323Z","shell.execute_reply":"2025-12-09T03:24:42.450823Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# IMPORTANT \n#Install dependencies and copy model weights to run the notebook without internet access when submitting to the competition.\n\n!pip install --no-index /kaggle/input/imc2024-packages-lightglue-rerun-kornia/* --no-deps\n!mkdir -p /root/.cache/torch/hub/checkpoints\n!cp /kaggle/input/aliked/pytorch/aliked-n16/1/aliked-n16.pth /root/.cache/torch/hub/checkpoints/\n!cp /kaggle/input/aliked/pytorch/aliked-n32/1/aliked-n32.pth /root/.cache/torch/hub/checkpoints/\n!cp /kaggle/input/aliked/pytorch/aliked-n16rot/1/aliked-n16rot.pth /root/.cache/torch/hub/checkpoints/\n!cp /kaggle/input/lightglue/pytorch/aliked/1/aliked_lightglue.pth /root/.cache/torch/hub/checkpoints/\n!cp /kaggle/input/lightglue/pytorch/aliked/1/aliked_lightglue.pth /root/.cache/torch/hub/checkpoints/aliked_lightglue_v0-1_arxiv-pth","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:32.713026Z","iopub.execute_input":"2025-04-17T13:00:32.713318Z","iopub.status.idle":"2025-04-17T13:00:38.467498Z","shell.execute_reply.started":"2025-04-17T13:00:32.713288Z","shell.execute_reply":"2025-04-17T13:00:38.466379Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import sys\nimport os, glob\nfrom tqdm import tqdm\nfrom fastprogress import progress_bar\nfrom time import time, sleep\nimport gc\nimport numpy as np\nimport h5py\nimport dataclasses\nimport pandas as pd\nfrom IPython.display import clear_output\nfrom collections import defaultdict\nfrom copy import deepcopy\nfrom PIL import Image\nimport networkx as nx\n\nimport cv2\nimport torch\nimport torch.nn.functional as F\nimport kornia as K\nimport kornia.feature as KF\n\nimport torch\nfrom lightglue import match_pair\nfrom lightglue import ALIKED, LightGlue\nfrom lightglue.utils import load_image, rbd\nfrom transformers import AutoImageProcessor, AutoModel\n\n# IMPORTANT Utilities: importing data into colmap and competition metric\nimport pycolmap\nsys.path.append('/kaggle/input/imc25-utils')\nfrom database import *\nfrom h5_to_db import *\nimport metric\n# Do not forget to select an accelerator on the sidebar to the right.\n#device = K.utils.get_cuda_device_if_available(0)\n#print(f'{device=}')\n\nimport concurrent.futures","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:38.468464Z","iopub.execute_input":"2025-04-17T13:00:38.468714Z","iopub.status.idle":"2025-04-17T13:00:59.413357Z","shell.execute_reply.started":"2025-04-17T13:00:38.468691Z","shell.execute_reply":"2025-04-17T13:00:59.412645Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Configurations","metadata":{}},{"cell_type":"code","source":"class Param_AlikedLightGlue:\n    def __init__(\n        self,\n        min_matches = 15,\n        max_num_keypoints = 4096,\n        image_size = 1024,\n    ):\n        self.min_matches = min_matches\n        self.max_num_keypoints = max_num_keypoints\n        self.image_size = image_size\n        \nclass CONFIG:\n    # Image pairs\n    sim_th = 2.0\n    min_pairs = 20\n    exhaustive_if_less = 20\n    GLOBAL_TOPK    = 150 # 50 # 近傍候補数\n    N_KEYPOINTS    = 2048          # 1画像あたり上限\n    RATIO_THR      = 1.0               # Lowe ratio\n    MATCH_THRESH   = 15                # good‑matches 閾値\n\n    # Image Matching\n    params_alikedlg = [\n        #Param_AlikedLightGlue( min_matches = 15, max_num_keypoints = 4096, image_size = 512 ),\n        #Param_AlikedLightGlue( min_matches = 15, max_num_keypoints = 4096, image_size = 1024 ),\n        #Param_AlikedLightGlue( min_matches = 15, max_num_keypoints = 8192, image_size = 1536 ),\n        Param_AlikedLightGlue( min_matches = 50, max_num_keypoints = 8192, image_size = 1536 ),\n        #Param_AlikedLightGlue( min_matches = 50, max_num_keypoints = 4096*3, image_size = 1536 ),\n    ]\n\n    roma_max_num_keypoints = 512\n    roma_image_size = 140\n    roma_batch_size = 20\n\n    # SfM\n    num_pallalel_sfm = 4","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.586021Z","iopub.execute_input":"2025-04-17T13:00:59.586364Z","iopub.status.idle":"2025-04-17T13:00:59.590351Z","shell.execute_reply.started":"2025-04-17T13:00:59.58633Z","shell.execute_reply":"2025-04-17T13:00:59.589535Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"device0=torch.device('cuda:0')\ndevice1=torch.device('cuda:1')\n\ndef switch_gpu(device, device_index):\n    if device is None:\n        return device0, 0\n    elif device == device0:\n        return device1, 1\n    else:\n        return device0, 0","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.591024Z","iopub.execute_input":"2025-04-17T13:00:59.591286Z","iopub.status.idle":"2025-04-17T13:00:59.602005Z","shell.execute_reply.started":"2025-04-17T13:00:59.591267Z","shell.execute_reply":"2025-04-17T13:00:59.601445Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"# Image Orientation Correction (Using Pre-trained Orientation ResNet)\"\"\"\n\nimport torchvision.models as models\nfrom PIL import Image, ImageOps\nimport torch.nn.functional as F\n\nclass OrientationCorrector:\n    def __init__(self, device):\n        self.device = device\n        print(\"Loading orientation detection model...\")\n        \n        # Load ResNeXt50_32x4d architecture\n        self.model = models.resnext50_32x4d(pretrained=False)\n        \n        # Modify final layer for 4-class orientation (0°, 90°, 180°, 270°)\n        num_ftrs = self.model.fc.in_features\n        self.model.fc = torch.nn.Linear(num_ftrs, 4)\n        \n        # Load the orientation-specific weights\n        weights_path = '/kaggle/input/resnet/pytorch/default/1/2020-11-16_resnext50_32x4d.pth'\n        try:\n            state_dict = torch.load(weights_path, map_location=device)\n            self.model.load_state_dict(state_dict, strict=False)\n            print(f\"✓ Loaded orientation detection weights from {weights_path}\")\n        except Exception as e:\n            print(f\"Error loading weights: {e}\")\n            raise\n        \n        self.model = self.model.to(device)\n        self.model.eval()\n        \n        # Standard ImageNet normalization (typically used for orientation models)\n        self.transform = transforms.Compose([\n            transforms.Resize((224, 224)),\n            transforms.ToTensor(),\n            transforms.Normalize([0.485, 0.456, 0.406], [0.229, 0.224, 0.225])\n        ])\n        \n        print(\"Orientation detection model loaded successfully\")\n    \n    def get_exif_orientation(self, img_pil):\n        \"\"\"Extract orientation from EXIF data as backup\"\"\"\n        try:\n            exif = img_pil.getexif()\n            if exif is not None:\n                orientation = exif.get(274, 1)  # 274 is the EXIF orientation tag\n                return orientation\n        except:\n            pass\n        return 1  # Normal orientation\n    \n    def predict_orientation(self, img_pil):\n        \"\"\"\n        Use trained model to predict orientation.\n        Returns: rotation class (0, 1, 2, 3) and confidence\n        \"\"\"\n        img_tensor = self.transform(img_pil).unsqueeze(0).to(self.device)\n        \n        with torch.no_grad():\n            output = self.model(img_tensor)\n            probs = F.softmax(output, dim=1)\n            confidence, predicted = torch.max(probs, 1)\n        \n        return predicted.item(), confidence.item()\n    \n    def correct_orientation(self, img_path):\n        \"\"\"\n        Detect and correct image orientation.\n        Returns: corrected PIL Image, rotation angle applied, method used, confidence\n        \"\"\"\n        img_pil = Image.open(img_path)\n        \n        # Convert to RGB if needed\n        if img_pil.mode != 'RGB':\n            img_pil = img_pil.convert('RGB')\n        \n        # Method 1: Check EXIF first (fast and reliable when available)\n        exif_orientation = self.get_exif_orientation(img_pil)\n        \n        if exif_orientation in [3, 6, 8]:\n            # Use EXIF-based correction\n            img_corrected = ImageOps.exif_transpose(Image.open(img_path))\n            if img_corrected.mode != 'RGB':\n                img_corrected = img_corrected.convert('RGB')\n            \n            angle = 0\n            if exif_orientation == 3:\n                angle = 180\n            elif exif_orientation == 6:\n                angle = 90\n            elif exif_orientation == 8:\n                angle = 270\n            \n            return img_corrected, angle, 'exif', 1.0\n        \n        # Method 2: Use trained orientation detection model\n        rotation_class, confidence = self.predict_orientation(img_pil)\n        \n        # Map class to degrees: 0=0°, 1=90°, 2=180°, 3=270°\n        rotation_angles = {0: 0, 1: 90, 2: 180, 3: 270}\n        angle = rotation_angles[rotation_class]\n        \n        # Apply rotation if needed\n        if angle != 0:\n            if angle == 90:\n                img_pil = img_pil.rotate(-90, expand=True)\n            elif angle == 180:\n                img_pil = img_pil.rotate(180, expand=True)\n            elif angle == 270:\n                img_pil = img_pil.rotate(90, expand=True)\n        \n        return img_pil, angle, 'model', confidence\n    \n    def correct_and_save_batch(self, img_paths, output_dir):\n        \"\"\"\n        Process a batch of images and save corrected versions.\n        Returns: dict mapping original path to correction info\n        \"\"\"\n        os.makedirs(output_dir, exist_ok=True)\n        corrections = {}\n        \n        stats = {\n            'exif': 0, \n            'model': 0, \n            'failed': 0,\n            'by_angle': {0: 0, 90: 0, 180: 0, 270: 0}\n        }\n        \n        for img_path in tqdm(img_paths, desc=\"Correcting orientations\"):\n            try:\n                img_pil, angle, method, confidence = self.correct_orientation(img_path)\n                \n                # Save corrected image\n                filename = os.path.basename(img_path)\n                corrected_path = os.path.join(output_dir, filename)\n                img_pil.save(corrected_path, quality=95)\n                \n                corrections[img_path] = {\n                    'corrected_path': corrected_path,\n                    'angle': angle,\n                    'was_corrected': angle != 0,\n                    'method': method,\n                    'confidence': confidence\n                }\n                \n                stats[method] += 1\n                stats['by_angle'][angle] += 1\n                \n            except Exception as e:\n                print(f\"Failed to correct {img_path}: {e}\")\n                # Copy original if correction fails\n                import shutil\n                corrected_path = os.path.join(output_dir, os.path.basename(img_path))\n                try:\n                    shutil.copy2(img_path, corrected_path)\n                except:\n                    # If copy fails, just use original path\n                    corrected_path = img_path\n                \n                corrections[img_path] = {\n                    'corrected_path': corrected_path,\n                    'angle': 0,\n                    'was_corrected': False,\n                    'method': 'failed',\n                    'confidence': 0.0\n                }\n                stats['failed'] += 1\n        \n        # Print detailed statistics\n        total = len(img_paths)\n        corrected = sum(1 for c in corrections.values() if c['was_corrected'])\n        \n        print(f\"\\n{'='*60}\")\n        print(f\"Orientation Correction Statistics:\")\n        print(f\"{'='*60}\")\n        print(f\"Total images processed: {total}\")\n        print(f\"Images corrected: {corrected} ({100*corrected/total:.1f}%)\")\n        print(f\"\\nCorrection method:\")\n        print(f\"  - EXIF metadata: {stats['exif']}\")\n        print(f\"  - Model prediction: {stats['model']}\")\n        print(f\"  - Failed: {stats['failed']}\")\n        print(f\"\\nRotations applied:\")\n        print(f\"  - 0° (no rotation): {stats['by_angle'][0]}\")\n        print(f\"  - 90° clockwise: {stats['by_angle'][90]}\")\n        print(f\"  - 180°: {stats['by_angle'][180]}\")\n        print(f\"  - 270° clockwise: {stats['by_angle'][270]}\")\n        \n        # Show confidence statistics for model predictions\n        model_confidences = [c['confidence'] for c in corrections.values() if c['method'] == 'model']\n        if model_confidences:\n            avg_conf = sum(model_confidences) / len(model_confidences)\n            min_conf = min(model_confidences)\n            max_conf = max(model_confidences)\n            print(f\"\\nModel prediction confidence:\")\n            print(f\"  - Average: {avg_conf:.3f}\")\n            print(f\"  - Min: {min_conf:.3f}, Max: {max_conf:.3f}\")\n        print(f\"{'='*60}\\n\")\n        \n        return corrections\n    \n    def cleanup(self):\n        \"\"\"Free up model memory\"\"\"\n        if hasattr(self, 'model'):\n            del self.model\n            torch.cuda.empty_cache()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Colmap Utilities","metadata":{}},{"cell_type":"code","source":"def import_into_colmap(img_dir, feature_dir ='.featureout', database_path = 'colmap.db'):\n    db = COLMAPDatabase.connect(database_path)\n    db.create_tables()\n    single_camera = False\n    fname_to_id = add_keypoints(db, feature_dir, img_dir, '', 'simple-radial', single_camera)\n    add_matches(\n        db,\n        feature_dir,\n        fname_to_id,\n    )\n    db.commit()\n    return","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.604375Z","iopub.execute_input":"2025-04-17T13:00:59.604574Z","iopub.status.idle":"2025-04-17T13:00:59.613341Z","shell.execute_reply.started":"2025-04-17T13:00:59.604556Z","shell.execute_reply":"2025-04-17T13:00:59.612634Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def add_keypoints_with_scene(db, scene, h5_path, image_path, img_ext, camera_model, single_camera = True):\n    keypoint_f = h5py.File(os.path.join(h5_path, 'keypoints.h5'), 'r')\n\n    camera_id = None\n    fname_to_id = {}\n    for filename in tqdm(list(keypoint_f.keys())):\n        if not filename in scene:\n            continue\n        \n        keypoints = keypoint_f[filename][()]\n\n        fname_with_ext = filename# + img_ext\n        path = os.path.join(image_path, fname_with_ext)\n        if not os.path.isfile(path):\n            raise IOError(f'Invalid image path {path}')\n\n        if camera_id is None or not single_camera:\n            camera_id = create_camera(db, path, camera_model)\n        image_id = db.add_image(fname_with_ext, camera_id)\n        fname_to_id[filename] = image_id\n\n        db.add_keypoints(image_id, keypoints)\n\n    return fname_to_id\n\ndef add_matches_with_scene(db, scene, h5_path, fname_to_id):\n    match_file = h5py.File(os.path.join(h5_path, 'matches.h5'), 'r')\n    \n    added = set()\n    n_keys = len(match_file.keys())\n    n_total = (n_keys * (n_keys - 1)) // 2\n\n    with tqdm(total=n_total) as pbar:\n        for key_1 in match_file.keys():\n            if not key_1 in scene:\n                continue\n\n            group = match_file[key_1]\n            for key_2 in group.keys():\n                if not key_2 in scene:\n                    continue\n                    \n                id_1 = fname_to_id[key_1]\n                id_2 = fname_to_id[key_2]\n\n                pair_id = image_ids_to_pair_id(id_1, id_2)\n                if pair_id in added:\n                    warnings.warn(f'Pair {pair_id} ({id_1}, {id_2}) already added!')\n                    continue\n            \n                matches = group[key_2][()]\n                db.add_matches(id_1, id_2, matches)\n\n                added.add(pair_id)\n\n                pbar.update(1)\n\ndef import_into_colmap_with_scene(img_dir, scene, feature_dir ='.featureout', database_path = 'colmap.db'):\n    db = COLMAPDatabase.connect(database_path)\n    db.create_tables()\n    single_camera = False\n    fname_to_id = add_keypoints(db, feature_dir, img_dir, '', 'simple-radial', single_camera)\n    #fname_to_id = add_keypoints_with_scene(db, scene, feature_dir, img_dir, '', 'simple-radial', single_camera)\n    add_matches_with_scene(\n        db,\n        scene,\n        feature_dir,\n        fname_to_id,\n    )\n    db.commit()\n    return","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Pairs Reranker","metadata":{}},{"cell_type":"code","source":"import os, torch, urllib.request, tempfile\nfrom PIL import Image\nfrom tqdm import tqdm\nfrom torchvision import transforms\n\n\n# ---------- Global embeddings ---------- #\ndef _extract_dino_embeddings(fnames, device = torch.device('cuda')):\n    processor = AutoImageProcessor.from_pretrained('/kaggle/input/dinov2/pytorch/base/1')\n    model = AutoModel.from_pretrained('/kaggle/input/dinov2/pytorch/base/1')\n    model = model.eval()\n    model = model.to(device)\n    global_descs_dinov2 = []\n    for i, img_fname_full in tqdm(enumerate(fnames),total= len(fnames)):\n        key = os.path.splitext(os.path.basename(img_fname_full))[0]\n        timg = load_torch_image(img_fname_full)\n        with torch.inference_mode():\n            inputs = processor(images=timg, return_tensors=\"pt\", do_rescale=False).to(device)\n            outputs = model(**inputs)\n            dino_mac = F.normalize(outputs.last_hidden_state[:,1:].max(dim=1)[0], dim=1, p=2)\n        global_descs_dinov2.append(dino_mac.detach().cpu())\n    global_descs_dinov2 = torch.cat(global_descs_dinov2, dim=0)\n    return global_descs_dinov2\n\ndef _build_topk_lists(global_feats, device):\n    \"\"\"\n    global_feats : (N, D)  L2‑normed\n    returns      : list of length‑≤(N‑1) neighbor indices for each image\n    \"\"\"\n    g   = global_feats.to(device)          # (N,D)\n    sim = g @ g.T                          # cosine similarity\n    sim.fill_diagonal_(-1)\n\n    N   = sim.size(0)\n    k   = min(CONFIG.GLOBAL_TOPK, N - 1)   # ★ ここが修正点\n    k   = max(k, 1)                        # 万一 N==2 でも k=1 を確保\n\n    topk = torch.topk(sim, k, dim=1).indices.cpu()\n    return [row.tolist() for row in topk]\n\n\n# ---------- ALIKED local features ---------- #\ndef _get_aliked_model(device, num_features, resize_to=1024, detection_threshold=0.01):\n    dtype = torch.float32 # ALIKED has issues with float16\n    extractor = ALIKED(\n        model_name=\"aliked-n16\",\n        max_num_keypoints=num_features,\n        detection_threshold=detection_threshold, \n        resize=resize_to\n    ).eval().to(device, dtype)\n    extractor.preprocess_conf[\"resize\"] = resize_to\n    return extractor\n    \n@torch.no_grad()\ndef _extract_local_features(img_paths, model, device):\n    dtype = torch.float32\n    descs = []\n    for p in tqdm(img_paths, desc=\"ALIKED\"):\n        image0 = load_torch_image(p, device=device).to(dtype)\n        h, w = image0.shape[2], image0.shape[3]\n        feats0 = model.extract(image0)  # auto-resize the image, disable with resize=None\n        d = feats0['descriptors'].reshape(-1, 128).detach()\n        descs.append(torch.nn.functional.normalize(d, dim=1).half().cpu())\n    return descs\n\n\n# ---------- good‑match カウント (Torch) ---------- #\n@torch.no_grad()\ndef _mutual_nn_score(desc1, desc2, device):\n    if desc1.size(0) == 0 or desc2.size(0) == 0:\n        return 0\n    d12 = torch.cdist(desc1.to(device), desc2.to(device), p=2)\n    min_val, _ = torch.min(d12, 1)\n    n_matches = np.sum( min_val.cpu().numpy() < (CONFIG.RATIO_THR ** 2) )\n    return n_matches\n\ndef _build_final_pairs(topk_lists, descs, device):\n    pairs = []\n    for i, nbrs in enumerate(tqdm(topk_lists, desc=\"Local verify\")):\n        for j in nbrs:\n            if i < j:\n                score = _mutual_nn_score(descs[i], descs[j], device)\n                if score >= CONFIG.MATCH_THRESH:\n                    pairs.append((i, j))\n                    #pairs.append((i, j, score))\n    pairs = sorted(list(set(pairs)))\n    return pairs\n\n\n# =========================================================\ndef get_image_pairs(fnames, device=torch.device(\"cuda\")):\n    \"\"\"\n    fnames: list[str] - 画像ファイルパスのリスト\n    device: torch.device - 計算に使用する GPU / CPU\n    returns: list[tuple[int,int]] - マッチングすべき画像 index ペア\n    \"\"\"\n    assert len(fnames) > 1, \"fnames must contain at least two images\"\n\n    # 1) DINO global features\n    global_feats = _extract_dino_embeddings(fnames, device)\n    topk_lists   = _build_topk_lists(global_feats, device)\n    cnt = 0\n    for topk_list in topk_lists:\n        cnt += len(topk_list)\n    print(cnt)\n\n    # 2) ALIKED local descriptors\n    aliked = _get_aliked_model(device, CONFIG.N_KEYPOINTS)\n    descs  = _extract_local_features(fnames, aliked, device)\n\n    # 3) local verification\n    pairs  = _build_final_pairs(topk_lists, descs, device)\n\n    del aliked\n    torch.cuda.empty_cache()\n    gc.collect()\n    return pairs\n# =========================================================","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Image pairs","metadata":{}},{"cell_type":"code","source":"def load_torch_image(fname, device=torch.device('cuda')):\n    img = K.io.load_image(fname, K.io.ImageLoadType.RGB32, device=device)[None, ...]\n    return img\n\n# Must Use efficientnet global descriptor to get matching shortlists.\ndef get_global_desc(fnames, device = torch.device('cuda')):\n    processor = AutoImageProcessor.from_pretrained('/kaggle/input/dinov2/pytorch/base/1')\n    model = AutoModel.from_pretrained('/kaggle/input/dinov2/pytorch/base/1')\n    model = model.eval()\n    model = model.to(device)\n    global_descs_dinov2 = []\n    for i, img_fname_full in tqdm(enumerate(fnames),total= len(fnames)):\n        key = os.path.splitext(os.path.basename(img_fname_full))[0]\n        timg = load_torch_image(img_fname_full)\n        with torch.inference_mode():\n            inputs = processor(images=timg, return_tensors=\"pt\", do_rescale=False).to(device)\n            outputs = model(**inputs)\n            dino_mac = F.normalize(outputs.last_hidden_state[:,1:].max(dim=1)[0], dim=1, p=2)\n        global_descs_dinov2.append(dino_mac.detach().cpu())\n    global_descs_dinov2 = torch.cat(global_descs_dinov2, dim=0)\n    return global_descs_dinov2\n\ndef get_img_pairs_exhaustive(img_fnames):\n    index_pairs = []\n    for i in range(len(img_fnames)):\n        for j in range(i+1, len(img_fnames)):\n            index_pairs.append((i,j))\n    return index_pairs\n\ndef get_image_pairs_shortlist(fnames,\n                              sim_th = 0.6, # should be strict\n                              min_pairs = 20,\n                              exhaustive_if_less = 20,\n                              device=torch.device('cuda')):\n    num_imgs = len(fnames)\n    if num_imgs <= exhaustive_if_less:\n        return get_img_pairs_exhaustive(fnames)\n    matching_list = get_image_pairs(fnames, device)\n    return matching_list\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.61466Z","iopub.execute_input":"2025-04-17T13:00:59.614904Z","iopub.status.idle":"2025-04-17T13:00:59.628751Z","shell.execute_reply.started":"2025-04-17T13:00:59.614883Z","shell.execute_reply":"2025-04-17T13:00:59.628036Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Aliked+LightGlue","metadata":{}},{"cell_type":"code","source":"def convert_coord(r, w, h, rotk):\n    if rotk == 0:\n        return r\n    elif rotk == 1:\n        rx = w-1-r[:, 1]\n        ry = r[:, 0]\n        return torch.concat([rx[None], ry[None]], dim=0).T # np.array([rx, ry]).T\n    elif rotk == 2:\n        rx = w-1-r[:, 0]\n        ry = h-1-r[:, 1]\n        return torch.concat([rx[None], ry[None]], dim=0).T # np.array([rx, ry]).T\n    elif rotk == 3:\n        rx = r[:, 1]\n        ry = h-1-r[:, 0]\n        return torch.concat([rx[None], ry[None]], dim=0).T # np.array([rx, ry]).T\n\ndef matching_aliked_lightglue_rot(\n    img_fnames,\n    index_pairs,\n    rot,\n    saved_file,\n    feature_dir = '.featureout',\n    num_features = 4096,\n    resize_to = 1024,\n    device=torch.device('cuda'),\n    min_matches=15,\n    verbose=True,\n):\n    if not os.path.isdir(feature_dir):\n        os.makedirs(feature_dir)\n\n    #####################################################\n    # Extract keypoints and descriptions\n    #####################################################\n    dtype = torch.float32 # ALIKED has issues with float16\n    extractor = ALIKED(\n        model_name=\"aliked-n16\",\n        max_num_keypoints=num_features,\n        detection_threshold=0.01, #0.001, \n        resize=resize_to\n    ).eval().to(device, dtype)\n    print(\"aliked> image size =\", extractor.preprocess_conf[\"resize\"], \"-->\", resize_to )\n    extractor.preprocess_conf[\"resize\"] = resize_to\n    if not os.path.isdir(feature_dir):\n        os.makedirs(feature_dir)\n    \n    dict_kpts_cuda = {}\n    dict_descs_cuda = {}\n    for img_path in img_fnames:\n        img_fname = img_path.split('/')[-1]\n        key = img_fname\n\n        with torch.inference_mode():\n            rot_k = 0\n            image0 = load_torch_image(img_path, device=device).to(dtype)\n            h, w = image0.shape[2], image0.shape[3]\n            feats0 = extractor.extract(image0)  # auto-resize the image, disable with resize=None\n            kpts = feats0['keypoints'].reshape(-1, 2).detach()\n            descs = feats0['descriptors'].reshape(len(kpts), -1).detach()\n            dict_kpts_cuda[f\"{key}_{rot_k}\"] = kpts\n            dict_descs_cuda[f\"{key}_{rot_k}\"] = descs\n            if verbose:\n                print(f\"aliked_rot> rot_k={rot_k}, kpts.shape={kpts.shape}, descs.shape={descs.shape}\")\n\n        if rot != 0:\n            with torch.inference_mode():\n                rot_k = rot\n                image0 = load_torch_image(img_path, device=device).to(dtype)\n                h, w = image0.shape[2], image0.shape[3]\n                image1 = torch.rot90(image0, rot, [2, 3])\n                feats0 = extractor.extract(image1)  # auto-resize the image, disable with resize=None\n                kpts = feats0['keypoints'].reshape(-1, 2).detach()\n                descs = feats0['descriptors'].reshape(len(kpts), -1).detach()\n                kpts = convert_coord(kpts, w, h, rot_k)\n                dict_kpts_cuda[f\"{key}_{rot_k}\"] = kpts\n                dict_descs_cuda[f\"{key}_{rot_k}\"] = descs\n                if verbose:\n                    print(f\"aliked_rot> rot_k={rot_k}, kpts.shape={kpts.shape}, descs.shape={descs.shape}\")\n    del extractor\n    gc.collect()\n\n    #####################################################\n    # Matching keypoints\n    #####################################################\n    lg_matcher = KF.LightGlueMatcher(\n        \"aliked\", {\n            \"width_confidence\": -1,\n            \"depth_confidence\": -1,\n            #\"filter_threshold\": 0.5,\n            \"mp\": True if 'cuda' in str(device) else False\n        }\n    ).eval().to(device).half()\n    \n    cnt_pairs = 0\n    with h5py.File(saved_file, mode='w') as f_match:\n        for pair_idx in tqdm(index_pairs):\n            idx1, idx2 = pair_idx\n            fname1, fname2 = img_fnames[idx1], img_fnames[idx2]\n            if (\"outliers\" in fname1) or ( \"outliers\" in fname2 ):\n                continue\n            key1, key2 = fname1.split('/')[-1], fname2.split('/')[-1]\n            \n            kp1 = dict_kpts_cuda[f\"{key1}_0\"]\n            desc1 = dict_descs_cuda[f\"{key1}_0\"]\n                        \n            kp2 = dict_kpts_cuda[f\"{key2}_{rot}\"]\n            desc2 = dict_descs_cuda[f\"{key2}_{rot}\"]\n            with torch.inference_mode():\n                dists, idxs = lg_matcher(desc1.half(),\n                                     desc2.half(),\n                                     KF.laf_from_center_scale_ori(kp1.half()[None]),\n                                     KF.laf_from_center_scale_ori(kp2.half()[None]))\n            #print( f\"dists.shape={dists.shape}\")\n            #print( f\"type(dists)={type(dists)}\")\n            #print( f\"dists.mean()={dists.mean()}, {dists.max()}, {dists.min()}\")\n            if len(idxs)  == 0:\n                continue\n            kp1 = kp1[idxs[:,0], :].cpu().numpy().reshape(-1, 2).astype(np.float32)\n            kp2 = kp2[idxs[:,1], :].cpu().numpy().reshape(-1, 2).astype(np.float32)\n            confs = (1.0-dists).cpu().numpy().reshape(-1, 1).astype(np.float32)\n            n_matches = kp1.shape[0]\n            group  = f_match.require_group(key1)\n            if n_matches >= min_matches:\n                matches = np.concatenate([kp1, kp2, confs], axis=1)\n                #print(f\"@@@@ matches.shape = {matches.shape}\")\n                group.create_dataset(key2, data=matches)\n                cnt_pairs+=1\n                if verbose:\n                    print (f'aliked_rot> {key1}-{key2}@rot={rot}: {n_matches} matches @ {cnt_pairs}th pair')            \n            else:\n                if verbose:\n                    print (f'aliked_rot> {key1}-{key2}: {n_matches} matches --> skipped')\n    del lg_matcher\n    torch.cuda.empty_cache()\n    gc.collect()\n    return\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.646519Z","iopub.execute_input":"2025-04-17T13:00:59.64679Z","iopub.status.idle":"2025-04-17T13:00:59.662745Z","shell.execute_reply.started":"2025-04-17T13:00:59.646761Z","shell.execute_reply":"2025-04-17T13:00:59.661951Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Keypoints merger & Matching Filter","metadata":{}},{"cell_type":"code","source":"def get_unique_idxs(A, dim=0):\n    # https://stackoverflow.com/questions/72001505/how-to-get-unique-elements-and-their-firstly-appeared-indices-of-a-pytorch-tenso\n    unique, idx, counts = torch.unique(A, dim=dim, sorted=True, return_inverse=True, return_counts=True)\n    _, ind_sorted = torch.sort(idx, stable=True)\n    cum_sum = counts.cumsum(0)\n    cum_sum = torch.cat((torch.tensor([0],device=cum_sum.device), cum_sum[:-1]))\n    first_indices = ind_sorted[cum_sum]\n    return first_indices\n\ndef get_keypoint_from_h5(fp, key1, key2):\n    rc = -1\n    try:\n        kpts = np.array(fp[key1][key2])\n        rc = 0\n        return (rc, kpts)\n    except:\n        return (rc, None)\n\ndef get_keypoint_from_multi_h5(fps, key1, key2):\n    list_mkpts = []\n    for fp in fps:\n        rc, mkpts = get_keypoint_from_h5(fp, key1, key2)\n        if mkpts is not None:\n            list_mkpts.append(mkpts)\n    if len(list_mkpts) > 0:\n        for data in list_mkpts:\n            print(f\"--> \", data.shape)\n        list_mkpts = np.concatenate(list_mkpts, axis=0)\n    else:\n        list_mkpts = None\n    return list_mkpts\n\ndef matches_merger(\n    img_fnames,\n    index_pairs,\n    files_keypoints,\n    save_file,\n    feature_dir = 'featureout',\n    filter_FundamentalMatrix = False,\n    filter_iterations = 10,\n    filter_threshold = 8,\n    min_matches=15,\n):\n    print( files_keypoints )\n    # open h5 files\n    fps = [ h5py.File(file, mode=\"r\") for file in files_keypoints ]\n    \n    with h5py.File(save_file, mode='w') as f_match:\n        counter = 0\n        for pair_idx in progress_bar(index_pairs):\n            idx1, idx2 = pair_idx\n            fname1, fname2 = img_fnames[idx1], img_fnames[idx2]\n            key1, key2 = fname1.split('/')[-1], fname2.split('/')[-1]\n\n            # extract keypoints\n            mkpts = get_keypoint_from_multi_h5(fps, key1, key2)\n            if mkpts is None:\n                print(f\"skipped key1={key1}, key2={key2}\")\n                continue\n\n            ori_size = mkpts.shape[0]\n            if mkpts.shape[0] < min_matches:\n                continue\n            \n            if filter_FundamentalMatrix:\n                store_inliers = { idx:0 for idx in range(mkpts.shape[0]) }\n                idxs = np.array(range(mkpts.shape[0]))\n                for iter in range(filter_iterations):\n                    try:\n                        Fm, inliers = cv2.findFundamentalMat(\n                            mkpts[:,:2], mkpts[:,2:4], cv2.USAC_MAGSAC, 0.15, 0.9999, 20000)\n                        if Fm is not None:\n                            inliers = inliers > 0\n                            inlier_idxs = idxs[inliers[:, 0]]\n                            #print(inliers.shape, inlier_idxs[:5])\n                            for idx in inlier_idxs:\n                                store_inliers[idx] += 1\n                    except:\n                        print(f\"Failed to cv2.findFundamentalMat. mkpts.shape={mkpts.shape}\")\n                inliers = np.array([ count for (idx, count) in store_inliers.items() ]) >= filter_threshold\n                mkpts = mkpts[inliers]\n                if mkpts.shape[0] < 15:\n                    print(f\"skipped key1={key1}, key2={key2}: mkpts.shape={mkpts.shape} after filtered.\")\n                    continue\n                #print(f\"filter_FundamentalMatrix: {len(store_inliers)} matches --> {mkpts.shape[0]} matches\")\n            \n            \n            print (f'{key1}-{key2}: {ori_size} --> {mkpts.shape[0]} matches')            \n            # regist tmp file\n            group  = f_match.require_group(key1)\n            group.create_dataset(key2, data=mkpts[:, :4])\n            counter += 1\n    print( f\"Ensembled pairs : {counter} pairs\" )\n    for fp in fps:\n        fp.close()\n\ndef get_adjacent_edges_with_weights(G, node):\n    \"\"\"\n    指定されたノードに隣接するエッジとその重みをリストで返す関数\n    \n    :param G: networkxのグラフオブジェクト\n    :param node: 対象となるノード\n    :return: (隣接ノード, 重み)のタプルのリスト、重みの降順でソートされる\n    \"\"\"\n    adjacent_edges = [(neighbor, G[node][neighbor]['weight']) for neighbor in G.neighbors(node)]\n    return sorted(adjacent_edges, key=lambda x: x[1], reverse=True)\n\ndef filter_mkpts(\n    file_mkpts, \n    save_file,\n    th_num_pairs,\n):\n    sleep(10)\n    G = nx.Graph()\n    pairs = []\n    with h5py.File(save_file, mode='w') as f_match, h5py.File(file_mkpts, \"r\") as fp:\n        for key1 in fp.keys():\n            for key2 in fp[key1].keys():\n                mkpts = fp[key1][key2]\n                G.add_edge(key1, key2, weight=mkpts.shape[0])\n                pairs.append([key1, key2, mkpts])\n\n        for (key1, key2, mkpts) in pairs:\n            neibors1 = get_adjacent_edges_with_weights(G, key1)\n            neibors2 = get_adjacent_edges_with_weights(G, key2)\n            if (key2 in [ val[0] for val in neibors1[:th_num_pairs] ]) or (key1 in [ val[0] for val in neibors2[:th_num_pairs] ]):\n                group  = f_match.require_group(key1)\n                group.create_dataset(key2, data=mkpts)\n\ndef keypoints_merger(\n    img_fnames,\n    index_pairs,\n    files_keypoints,\n    feature_dir = 'featureout',\n    filter_FundamentalMatrix = False,\n    filter_iterations = 10,\n    filter_threshold = 8,\n):\n    print(f\"files_keypoints = {files_keypoints}\")\n    save_file0 = f'{feature_dir}/merge_tmp0.h5'\n    !rm -rf {save_file0}\n    matches_merger(\n        img_fnames,\n        index_pairs,\n        files_keypoints,\n        save_file0,\n        feature_dir = feature_dir,\n        filter_FundamentalMatrix = filter_FundamentalMatrix,\n        filter_iterations = filter_iterations,\n        filter_threshold = filter_threshold,\n    )\n\n    th_num_pairs = 6\n    save_file = f'{feature_dir}/merge_tmp.h5'\n    !rm -rf {save_file}\n    filter_mkpts(\n        save_file0, \n        save_file,\n        th_num_pairs,\n    )\n    \n    # Let's find unique loftr pixels and group them together.\n    kpts = defaultdict(list)\n    match_indexes = defaultdict(dict)\n    total_kpts=defaultdict(int)\n    with h5py.File(save_file, mode='r') as f_match:\n        for k1 in f_match.keys():\n            group  = f_match[k1]\n            for k2 in group.keys():\n                matches = group[k2][...]\n                total_kpts[k1]\n                kpts[k1].append(matches[:, :2])\n                kpts[k2].append(matches[:, 2:])\n                current_match = torch.arange(len(matches)).reshape(-1, 1).repeat(1, 2)\n                current_match[:, 0]+=total_kpts[k1]\n                current_match[:, 1]+=total_kpts[k2]\n                total_kpts[k1]+=len(matches)\n                total_kpts[k2]+=len(matches)\n                match_indexes[k1][k2]=current_match\n\n    for k in kpts.keys():\n        kpts[k] = np.round(np.concatenate(kpts[k], axis=0))\n    unique_kpts = {}\n    unique_match_idxs = {}\n    out_match = defaultdict(dict)\n    for k in kpts.keys():\n        uniq_kps, uniq_reverse_idxs = torch.unique(torch.from_numpy(kpts[k]),dim=0, return_inverse=True)\n        unique_match_idxs[k] = uniq_reverse_idxs\n        unique_kpts[k] = uniq_kps.numpy()\n    for k1, group in match_indexes.items():\n        for k2, m in group.items():\n            m2 = deepcopy(m)\n            m2[:,0] = unique_match_idxs[k1][m2[:,0]]\n            m2[:,1] = unique_match_idxs[k2][m2[:,1]]\n            mkpts = np.concatenate([unique_kpts[k1][ m2[:,0]],\n                                    unique_kpts[k2][  m2[:,1]],\n                                   ],\n                                   axis=1)\n            unique_idxs_current = get_unique_idxs(torch.from_numpy(mkpts), dim=0)\n            m2_semiclean = m2[unique_idxs_current]\n            unique_idxs_current1 = get_unique_idxs(m2_semiclean[:, 0], dim=0)\n            m2_semiclean = m2_semiclean[unique_idxs_current1]\n            unique_idxs_current2 = get_unique_idxs(m2_semiclean[:, 1], dim=0)\n            m2_semiclean2 = m2_semiclean[unique_idxs_current2]\n            out_match[k1][k2] = m2_semiclean2.numpy()\n    with h5py.File(f'{feature_dir}/keypoints.h5', mode='w') as f_kp:\n        for k, kpts1 in unique_kpts.items():\n            f_kp[k] = kpts1\n    \n    with h5py.File(f'{feature_dir}/matches.h5', mode='w') as f_match:\n        for k1, gr in out_match.items():\n            group  = f_match.require_group(k1)\n            for k2, match in gr.items():\n                group[k2] = match\n    return","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.663535Z","iopub.execute_input":"2025-04-17T13:00:59.663808Z","iopub.status.idle":"2025-04-17T13:00:59.692456Z","shell.execute_reply.started":"2025-04-17T13:00:59.663788Z","shell.execute_reply":"2025-04-17T13:00:59.691754Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Pipeline: Image matching","metadata":{}},{"cell_type":"code","source":"def exec_feature_extraction(params):\n    images = params.images\n    feature_dir = params.feature_dir\n    device = params.device\n\n    # **NEW: Correct image orientations first**\n    print(f\"\\n{'='*60}\")\n    print(f\"STEP 1: Orientation Detection and Correction\")\n    print(f\"{'='*60}\")\n    \n    t_orient = time()\n    orientation_corrector = OrientationCorrector(device)\n    corrected_dir = os.path.join(feature_dir, 'corrected_images')\n    corrections = orientation_corrector.correct_and_save_batch(images, corrected_dir)\n    \n    # Update image paths to use corrected versions\n    corrected_images = [corrections[img]['corrected_path'] for img in images]\n    params.corrected_images = corrected_images\n    params.orientation_corrections = corrections\n    \n    # Cleanup orientation model to free memory\n    orientation_corrector.cleanup()\n    del orientation_corrector\n    torch.cuda.empty_cache()\n    gc.collect()\n    \n    params.laptime_orientation = time() - t_orient\n    print(f\"Orientation correction completed in {params.laptime_orientation:.2f} sec\\n\")\n    \n    # **Continue with existing pipeline using corrected images**\n    print(f\"{'='*60}\")\n    print(f\"STEP 2: Image Pair Selection\")\n    print(f\"{'='*60}\")\n    \n    t = time()\n    index_pairs = get_image_pairs_shortlist(\n        corrected_images,  # **CHANGED: use corrected images**\n        sim_th = CONFIG.sim_th,\n        min_pairs = CONFIG.min_pairs,\n        exhaustive_if_less = CONFIG.exhaustive_if_less,\n        device=device\n    )\n    params.laptime_image_pairs = time() - t\n    params.image_pairs = index_pairs\n    print (f'Shortlisting. Number of pairs to match: {len(index_pairs)}. Done in {params.laptime_image_pairs:.4f} sec\\n')\n    gc.collect()\n\n    print(f\"{'='*60}\")\n    print(f\"STEP 3: Feature Extraction and Matching\")\n    print(f\"{'='*60}\")\n    \n    t = time()\n    files_keypoints = []\n    for i_param, imc_param in enumerate(CONFIG.params_alikedlg):\n        for rot in range(1):\n            file_keypoints = f'{feature_dir}/aliked_lightglue_param{i_param}_rot{rot}_mkpts.h5'\n            max_num_keypoints = imc_param.max_num_keypoints\n            resize_to = imc_param.image_size\n            min_matches = imc_param.min_matches\n            matching_aliked_lightglue_rot(\n                corrected_images,  # **CHANGED: use corrected images**\n                index_pairs,\n                rot,\n                file_keypoints,\n                feature_dir,\n                max_num_keypoints,\n                resize_to=resize_to,\n                device=device,\n                min_matches=min_matches,\n                verbose=True,\n            )\n            files_keypoints.append( file_keypoints )\n            gc.collect()\n    params.laptime_lightglue_4rots = time() - t\n    print(f'Features matched in {params.laptime_lightglue_4rots:.4f} sec\\n')\n\n    # ... rest of your code remains the same ...\n\n   \n\n    # RoMa\n    \"\"\"\n    t = time()\n    file_keypoints = f'{feature_dir}/roma_mkpts.h5'\n    matching_roma(\n        images,\n        index_pairs,\n        file_keypoints,\n        feature_dir,\n        num_features=CONFIG.roma_max_num_keypoints,\n        resize_to=CONFIG.roma_image_size,\n        device=device,\n        min_matches=CONFIG.min_matches,\n        verbose=True,\n    )\n    files_keypoints.append( file_keypoints )\n    gc.collect()\n    params.laptime_roma = time() - t\n    print(f'Features matched in {params.laptime_roma:.4f} sec')\n    \"\"\"\n    \n    # Merge mkpts and create keypoints/matches\n    sleep(10)\n    keypoints_merger(\n        images,\n        index_pairs,\n        files_keypoints,\n        feature_dir,\n        filter_FundamentalMatrix = False,\n        filter_iterations = 10,\n        filter_threshold = 8,\n    )\n    \n    gc.collect()\n    return params","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.746932Z","iopub.execute_input":"2025-04-17T13:00:59.747163Z","iopub.status.idle":"2025-04-17T13:00:59.753465Z","shell.execute_reply.started":"2025-04-17T13:00:59.747143Z","shell.execute_reply":"2025-04-17T13:00:59.752642Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# SfM Maps Merger","metadata":{}},{"cell_type":"code","source":"# Merge maps of colmap\nimport argparse, shutil, subprocess, tempfile, sys\nfrom pathlib import Path\n\nimport networkx as nx\nimport pycolmap  # pip install pycolmap\n\n# ---------- util ----------------------------------------------------------------\ndef read_image_names(model_dir: Path):\n    \"\"\"\n    pycolmap 0.6.x では Reconstruction() で読み込む。\n    返り値: 画像ファイル名(set[str])\n    \"\"\"\n    rec = pycolmap.Reconstruction(str(model_dir))   # ← これで images.bin をパース\n    return {img.name for img in rec.images.values()}\n\n#def read_image_names(model_dir: Path):\n#    images = pycolmap.read_images_binary(model_dir / \"images.bin\")\n#    return {img.name for img in images.values()}\n\ndef build_overlap_graph(image_sets, min_common):\n    G = nx.Graph()\n    G.add_nodes_from(range(len(image_sets)))\n    for i in range(len(image_sets)):\n        for j in range(i + 1, len(image_sets)):\n            if len(image_sets[i] & image_sets[j]) >= min_common:\n                G.add_edge(i, j)\n    return G\n\ndef run_cmd(cmd):\n    \"\"\"stdout+stderr をストリーム出力しつつ実行\"\"\"\n    proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)\n    for line in proc.stdout:\n        print(line, end=\"\")\n    proc.wait()\n    if proc.returncode != 0:\n        raise subprocess.CalledProcessError(proc.returncode, cmd)\n\n# ---------- COLMAP wrappers ------------------------------------------------------\ndef merge_two_models(colmap_bin, ref_model, src_model, out_dir):\n    cmd = [\n        colmap_bin, \"model_merger\",\n        \"--input_path1\", str(ref_model),\n        \"--input_path2\", str(src_model),\n        \"--output_path\", str(out_dir),\n    ]\n    try:\n        run_cmd(cmd)\n        return True\n    except subprocess.CalledProcessError:\n        return False\n\ndef run_bundle_adjuster(colmap_bin, model_dir, num_threads=-1):\n    \"\"\"\n    実行例:\n        colmap bundle_adjuster --input_path scene_000 --output_path scene_000 \\\n            --BundleAdjustment.refine_principal_point 1\n    \"\"\"\n    cmd = [\n        colmap_bin, \"bundle_adjuster\",\n        \"--input_path\", str(model_dir),\n        \"--output_path\", str(model_dir),\n        \"--BundleAdjustment.refine_focal_length\", \"1\",\n        \"--BundleAdjustment.refine_principal_point\", \"1\",\n        \"--BundleAdjustment.refine_extra_params\", \"1\",\n    ]\n    if num_threads > 0:\n        cmd += [\"--Mapper.num_threads\", str(num_threads)]\n    run_cmd(cmd)     # raise on error\n\n# ---------- main logic -----------------------------------------------------------\ndef cluster_and_merge(\n    model_dirs,\n    out_root,\n    colmap_bin=\"colmap\",\n    min_common=10,\n    run_bundle_adjustment=True,\n    keep_tmp=False,\n):\n    out_root = Path(out_root)\n    out_root.mkdir(parents=True, exist_ok=True)\n\n    # 1) 画像集合\n    image_sets = [read_image_names(Path(m)) for m in model_dirs]\n\n    # 2) 重複グラフ → 連結成分\n    clusters = list(nx.connected_components(build_overlap_graph(image_sets, min_common)))\n    print(f\"num of clusters = {len(clusters)}\")\n    \n    merged_paths = []\n    for cidx, comp in enumerate(clusters):\n        comp = list(comp)\n        if len(comp) == 1:\n            src = Path(model_dirs[comp[0]])\n            dst = out_root / f\"scene_{cidx:03d}\"\n            shutil.copytree(src, dst, dirs_exist_ok=True)\n            merged_paths.append(dst)\n            if run_bundle_adjustment:\n                run_bundle_adjuster(colmap_bin, dst)\n            continue\n\n        # 3) 画像数の多い順に\n        comp_sorted = sorted(comp, key=lambda i: len(image_sets[i]), reverse=True)\n        work_dir = Path(model_dirs[comp_sorted[0]])\n\n        # 4) 順次マージ\n        failed_sources = []\n        for midx in comp_sorted[1:]:\n            tmp_dir = Path(tempfile.mkdtemp(prefix=f\"merge_s{cidx}_\"))\n            ok = merge_two_models(colmap_bin, work_dir, Path(model_dirs[midx]), tmp_dir)\n            if ok:\n                work_dir = tmp_dir\n            else:\n                failed_sources.append(midx)\n\n        # 5) 出力 & BA\n        final_dir = out_root / f\"scene_{cidx:03d}\"\n        shutil.move(str(work_dir), final_dir)\n        merged_paths.append(final_dir)\n\n        if run_bundle_adjustment:\n            run_bundle_adjuster(colmap_bin, final_dir)\n\n        # 6) マージ失敗分を別シーンとしてコピー\n        for fidx in failed_sources:\n            dst = out_root / f\"scene_{cidx}_{fidx}\"\n            shutil.copytree(model_dirs[fidx], dst, dirs_exist_ok=True)\n            merged_paths.append(dst)\n            if run_bundle_adjustment:\n                run_bundle_adjuster(colmap_bin, dst)\n\n    return merged_paths","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Update submission values","metadata":{}},{"cell_type":"code","source":"def update_prediction(params, result_map_dirs, retry):\n    # Parse\n    feature_dir = params.feature_dir\n    images_dir = params.images_dir\n    filename_to_index = params.filename_to_index\n    predictions = params.predictions\n    dataset = params.dataset\n\n    maps = {}\n    for idx, result_map_dir in enumerate(result_map_dirs):\n        maps[idx] = pycolmap.Reconstruction(result_map_dir)\n    \n    print (\"Counting map size...\")\n    list_num_images = []\n    if isinstance(maps, dict):\n        for idx1, rec in maps.items():\n            list_num_images.append( len(rec.images) )\n    list_num_images = np.array(list_num_images)\n    print(f\"list_num_images = {list_num_images}\")\n    if params.list_model_size is None:\n        params.list_model_size = [list_num_images]\n    else:\n        params.list_model_size.append( list_num_images )\n    sort_idx = np.argsort( np.array(list_num_images) )\n    \n    registered = 0\n    for map_index, cur_idx in enumerate(sort_idx):\n        cur_map = maps[cur_idx]\n        cur_map_size = list_num_images[cur_idx]\n        if cur_map_size < 4:\n            continue\n        for index, image in cur_map.images.items():\n            prediction_index = filename_to_index[image.name]\n            if (cur_map_size > predictions[prediction_index].map_size) and (cur_map_size > 5):\n                predictions[prediction_index].map_size = cur_map_size\n                predictions[prediction_index].cluster_index = map_index\n                predictions[prediction_index].retry_index = retry\n                predictions[prediction_index].rotation = deepcopy(image.cam_from_world.rotation.matrix())\n                predictions[prediction_index].translation = deepcopy(image.cam_from_world.translation)\n                registered += 1\n    del maps\n    gc.collect()\n    return params","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Clustering","metadata":{}},{"cell_type":"code","source":"import networkx as nx\nfrom networkx.algorithms.community import louvain_communities\nfrom sklearn.cluster import SpectralClustering\nfrom sklearn.cluster import AgglomerativeClustering\nfrom typing import List, Tuple, Iterable\n\ndef get_network_from_matches_h5( file, images, num_max_keypoints = 8192, th_matches=150):\n    image_to_index = {file:i for i,file in enumerate(images)}\n    index_to_image = {i:file for i,file in enumerate(images)}\n\n    edges = []\n    with h5py.File(file, \"r\") as f_mat:\n        for key1 in f_mat.keys():\n            for key2 in f_mat[key1].keys():\n                if f_mat[ key1 ][ key2 ].shape[0] >=th_matches:\n                    edges.append( (key1, key2, f_mat[key1][key2].shape[0] / num_max_keypoints) )\n    G = nx.Graph()\n    G.add_weighted_edges_from(edges, weight=\"weight\")\n    return G, image_to_index, index_to_image\n\ndef get_components(\n    G: nx.Graph | nx.DiGraph,\n    *,\n    directed_mode: str = \"auto\"  # \"auto\" | \"weak\" | \"strong\"\n) -> Tuple[int, List[List]]:\n    # --- 無向グラフ ----------------------------------------------------\n    if not G.is_directed():\n        comp_iter: Iterable[set] = nx.connected_components(G)\n\n    # --- 有向グラフ ----------------------------------------------------\n    else:\n        if directed_mode == \"strong\":\n            comp_iter = nx.strongly_connected_components(G)\n        elif directed_mode in {\"weak\", \"auto\"}:\n            comp_iter = nx.weakly_connected_components(G)\n        else:\n            raise ValueError(\"directed_mode must be 'auto', 'weak', or 'strong'\")\n\n    components = [sorted(list(c)) for c in comp_iter]\n    n_components = len(components)\n    return n_components, components\n\ndef get_scenes_from_graph(G, random_state=42, th=0.2):\n    communities = louvain_communities(\n        G,\n        weight=\"weight\",\n        resolution=th,\n        seed=random_state\n    )\n    \n    scenes = []\n    print(f\"#clusters = {len(communities)}\")\n    for cid, members in enumerate(communities):\n        scene = set(members)\n        print(f\"Cluster {cid}: {len(members)} | {sorted(members)}\")\n        if len(scene) >= 5:\n            scenes.append( scene )\n        print(\"=\"*50)\n    return scenes","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Pipeline: Multi-Scene SfM","metadata":{}},{"cell_type":"code","source":"def exec_reconstruction(params):\n    feature_dir = params.feature_dir\n    params.laptime_calc_fmat = []\n    params.laptime_sfm = []\n\n    # Clustering\n    feature_dir = params.feature_dir\n    fnames = [ p.filename for p in params.predictions ]\n\n    matches_file = f'{feature_dir}/matches.h5'\n    scenes = []\n    num_retry = 4\n    for _ in range(num_retry):\n        G, images_to_index, index_to_image = get_network_from_matches_h5( \n            matches_file, fnames, num_max_keypoints = 8192, th_matches=150,\n        )\n        n_components, components = get_components(G)\n        for component in components:\n            if len(component) > 4:\n                scenes.append( set(component) )\n\n    # SfM for each scne cluster\n    for retry, scene in enumerate(scenes):\n        params = reconstruction_single(params, scene, retry, 2)\n\n    # merge models\n    input_map_dirs = sorted( glob.glob( f\"{feature_dir}/colmap_rec_*/*/cameras.bin\"))\n    input_map_dirs = [ \n        os.path.dirname(x)\n        for x in input_map_dirs\n        if len(pycolmap.Reconstruction( os.path.dirname(x) ).images) > 8\n    ]\n    input_map_dirs = [ \n        x \n        for x in input_map_dirs\n        if pycolmap.Reconstruction(x).compute_mean_reprojection_error() < 2.0\n    ]\n\n    if len(input_map_dirs) > 1:\n        output_dir = f\"{feature_dir}/colmap_rec_merged/\"\n        COLMAP_BIN = \"/usr/bin/colmap\"\n        MERGE_MIN_COMMON = 6\n        MERGE_RUN_BA = True\n        result_map_dirs = cluster_and_merge(\n            input_map_dirs,\n            output_dir,\n            colmap_bin=COLMAP_BIN,\n            min_common=MERGE_MIN_COMMON,\n            run_bundle_adjustment=MERGE_RUN_BA,\n            keep_tmp=False,\n        )\n    \n        # update motion info\n        update_prediction(params, result_map_dirs, len(scenes) )\n    clear_output(wait=False)\n    return params\n\ndef reconstruction_single(params, scene, retry, min_model_size):\n    # Parse\n    feature_dir = params.feature_dir\n    images_dir = params.images_dir\n    filename_to_index = params.filename_to_index\n    predictions = params.predictions\n\n    # Execute process\n    database_path = os.path.join(feature_dir, f'colmap_{retry}.db')\n    if os.path.isfile(database_path):\n        os.remove(database_path)\n    gc.collect()\n    sleep(1)\n    #import_into_colmap(images_dir, feature_dir=feature_dir, database_path=database_path)\n    import_into_colmap_with_scene(images_dir, scene, feature_dir=feature_dir, database_path=database_path)\n    output_path = f'{feature_dir}/colmap_rec_{retry}'\n    \n    t = time()\n    pycolmap.match_exhaustive(database_path)\n    params.laptime_calc_fmat = time() - t\n    print(f'Ran RANSAC in {params.laptime_calc_fmat:.4f} sec')\n    \n    # By default colmap does not generate a reconstruction if less than 10 images are registered.\n    # Lower it to 3.\n    mapper_options = pycolmap.IncrementalPipelineOptions()\n    mapper_options.min_model_size = min_model_size #25 # 30 # 20 # 10 # 3\n    mapper_options.max_num_models = 15 # 25 # 25\n\n    # Strict initial pairs\n    #mapper_options.mapper.init_min_num_inliers = 200 # default:100\n    #mapper_options.mapper.init_max_error = 2         # default:4\n    #mapper_options.mapper.init_min_tri_angle = 10    # default:16\n\n    ##\n    #mapper_options.max_model_overlap = 30 # 20\n    #mapper_options.min_num_matches   = 50\n    ##\n    mapper_options.init_num_trials = 1000 # 200\n    \n    os.makedirs(output_path, exist_ok=True)\n    t = time()\n    maps = pycolmap.incremental_mapping(\n        database_path=database_path, \n        image_path=images_dir,\n        output_path=output_path,\n        options=mapper_options)\n    sleep(1)\n    params.laptime_sfm = time() - t\n    print(f'Reconstruction done in  {params.laptime_sfm:.4f} sec')\n    print(maps)\n\n    clear_output(wait=False)\n\n    #print(\"merge maps...\")\n    #print(f\"before: len(maps) = \", len(maps))\n    #maps = merge_models_dict(maps)\n    #print(f\"after: len(maps) = \", len(maps))\n\n    print (\"Counting map size...\")\n    list_num_images = []\n    if isinstance(maps, dict):\n        for idx1, rec in maps.items():\n            list_num_images.append( len(rec.images) )\n    list_num_images = np.array(list_num_images)\n    print(f\"list_num_images = {list_num_images}\")\n    if params.list_model_size is None:\n        params.list_model_size = [list_num_images]\n    else:\n        params.list_model_size.append( list_num_images )\n    sort_idx = np.argsort( np.array(list_num_images) )\n\n    registered = 0\n    for map_index, cur_idx in enumerate(sort_idx):\n        cur_map = maps[cur_idx]\n        cur_map_size = list_num_images[cur_idx]\n        for index, image in cur_map.images.items():\n            prediction_index = filename_to_index[image.name]\n            if (cur_map_size > predictions[prediction_index].map_size) and (cur_map_size > 5):\n                predictions[prediction_index].map_size = cur_map_size\n                predictions[prediction_index].cluster_index = map_index\n                predictions[prediction_index].retry_index = retry\n                predictions[prediction_index].rotation = deepcopy(image.cam_from_world.rotation.matrix())\n                predictions[prediction_index].translation = deepcopy(image.cam_from_world.translation)\n                registered += 1\n    gc.collect()\n\n    return params","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.754294Z","iopub.execute_input":"2025-04-17T13:00:59.754565Z","iopub.status.idle":"2025-04-17T13:00:59.767964Z","shell.execute_reply.started":"2025-04-17T13:00:59.754538Z","shell.execute_reply":"2025-04-17T13:00:59.767342Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Main","metadata":{}},{"cell_type":"code","source":"# Collect vital info from the dataset\n\n\n@dataclasses.dataclass\nclass Prediction:\n    image_id: str | None  # A unique identifier for the row -- unused otherwise. Used only on the hidden test set.\n    dataset: str\n    filename: str\n    map_size: int = -1\n    cluster_index: int | None = None\n    retry_index: int = -1\n    rotation: np.ndarray | None = None\n    translation: np.ndarray | None = None\n\n@dataclasses.dataclass\nclass DatasetParams:\n    dataset: str\n    feature_dir: str | None = None\n    images_dir: str | None = None\n    filename_to_index: dict | None = None\n    predictions: list | None = None\n    images: list | None = None\n    device: str | None = None\n    laptime_image_pairs: float = -1.0\n    laptime_lightglue_4rots: float = -1.0\n    laptime_roma: float = -1.0\n    laptime_calc_fmat: float = -1.0\n    laptime_sfm: float = -1.0\n    image_pairs: list | None = None\n    list_model_size: list | None = None\n\n    def summary(self):\n        print(\"[Summary]\")\n        print(f\"- Dataset                   : {self.dataset}\")\n        if self.images is not None:\n            print(f\"- images                    : {len(self.images)}\")\n        if self.image_pairs is not None:\n            print(f\"- pairs                     : {len(self.image_pairs)}\")\n        print(f\"- laptime_image_pairs       : {self.laptime_image_pairs:.2f}\")\n        print(f\"- laptime_lightglue_4rots   : {self.laptime_lightglue_4rots:.2f}\")\n        print(f\"- laptime_roma              : {self.laptime_roma:.2f}\")\n        print(f\"- laptime_calc_fmat         : {self.laptime_calc_fmat:.2f}\")\n        print(f\"- laptime_sfm               : {self.laptime_sfm:.2f}\")\n        print(f\"- list_model_size           : {self.list_model_size}\")\n\n# Set is_train=True to run the notebook on the training data.\n# Set is_train=False if submitting an entry to the competition (test data is hidden, and different from what you see on the \"test\" folder).\nis_train = False\ndata_dir = '/kaggle/input/image-matching-challenge-2025'\nworkdir = '/tmp/kaggle/working/result/'\n\nif is_train:\n    sample_submission_csv = os.path.join(data_dir, 'train_labels.csv')\nelse:\n    sample_submission_csv = os.path.join(data_dir, 'sample_submission.csv')\n\nsamples = {}\ndataset_logs = []\ncompetition_data = pd.read_csv(sample_submission_csv)\nif (not is_train) and (competition_data.shape[0] == 1945):\n    competition_data = competition_data[ competition_data[\"dataset\"].isin([\"ETs\", \"stairs\"]) ]\n    workdir = '/kaggle/working/result/'\n\nos.makedirs(workdir, exist_ok=True)\n\nfor _, row in competition_data.iterrows():\n    # Note: For the test data, the \"scene\" column has no meaning, and the rotation_matrix and translation_vector columns are random.\n    if row.dataset not in samples:\n        samples[row.dataset] = []\n    prediction = Prediction(\n        image_id=None if is_train else row.image_id,\n        dataset=row.dataset,\n        filename=row.image\n    )\n\n    #if len(samples[row.dataset]) < 10:\n    samples[row.dataset].append( prediction )\n\nfor dataset in samples:\n    print(f'Dataset \"{dataset}\" -> num_images={len(samples[dataset])}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.693268Z","iopub.execute_input":"2025-04-17T13:00:59.693498Z","iopub.status.idle":"2025-04-17T13:00:59.746332Z","shell.execute_reply.started":"2025-04-17T13:00:59.693479Z","shell.execute_reply":"2025-04-17T13:00:59.745721Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gc.collect()\n\nmax_images = None  # Used For debugging only. Set to None to disable.\ndatasets_to_process = None  # Not the best convention, but None means all datasets.\n\nif is_train:\n    # max_images = 5\n\n    # Note: When running on the training dataset, the notebook will hit the time limit and die. Use this filter to run on a few specific datasets.\n    datasets_to_process = [\n    \t# New data.\n    \t#'amy_gardens',\n    \t'ETs',\n    \t#'fbk_vineyard',\n    \t'stairs',\n    \t# Data from IMC 2023 and 2024.\n    \t# 'imc2024_dioscuri_baalshamin',\n    \t# 'imc2023_theather_imc2024_church',\n    \t# 'imc2023_heritage',\n    \t 'imc2023_haiper',\n    \t# 'imc2024_lizard_pond',\n    \t# Crowdsourced PhotoTourism data.\n    \t# 'pt_stpeters_stpauls',\n    \t# 'pt_brandenburg_british_buckingham',\n    \t# 'pt_piazzasanmarco_grandplace',\n    \t# 'pt_sacrecoeur_trevi_tajmahal',\n    ]\n\nmapping_result_strs = []\n\n\n# Enqeue feature extraction processing\nfutures_gpu0 = {}\nfutures_gpu1 = {}\ndevice = None\ndevice_index = None\nwith concurrent.futures.ThreadPoolExecutor(max_workers=1) as executors_gpu0, \\\n     concurrent.futures.ThreadPoolExecutor(max_workers=1) as executors_gpu1:\n\n    executors_gpus = [executors_gpu0, executors_gpu1]\n    futures_gpus = [futures_gpu0, futures_gpu1]\n    for dataset, predictions in samples.items():\n        if datasets_to_process and dataset not in datasets_to_process:\n            print(f'Skipping \"{dataset}\"')\n            continue\n        \n        images_dir = os.path.join(data_dir, 'train' if is_train else 'test', dataset)\n        images = [os.path.join(images_dir, p.filename) for p in predictions]\n        if max_images is not None:\n            images = images[:max_images]\n    \n        print(f'\\nProcessing dataset \"{dataset}\": {len(images)} images')\n    \n        filename_to_index = {p.filename: idx for idx, p in enumerate(predictions)}\n    \n        feature_dir = os.path.join(workdir, 'featureout', dataset)\n        os.makedirs(feature_dir, exist_ok=True)\n\n        # Switch device\n        device, device_index = switch_gpu(device, device_index)\n        \n        # Dataset parameter\n        dataset_params = DatasetParams(\n            dataset = dataset,\n            feature_dir = feature_dir,\n            images_dir = images_dir,\n            predictions = predictions,\n            filename_to_index = filename_to_index,\n            images = images,  # image filepaths\n            device = device,\n        )\n\n        # Enqueue feature extraction\n        futures_gpus[device_index][dataset] = executors_gpus[device_index].submit(\n            exec_feature_extraction, dataset_params,\n        )\n        #dataset_params = exec_feature_extraction(dataset_params)\n\n    # Enqeue reconstruction processing\n    futures_cpu  = {}\n    with concurrent.futures.ThreadPoolExecutor(max_workers=CONFIG.num_pallalel_sfm) as executors:\n        for dataset, predictions in samples.items():\n            if datasets_to_process and (dataset not in datasets_to_process):\n                print(f'Skipping \"{dataset}\"')\n                continue\n                \n            if dataset in futures_gpu0.keys():\n                future = futures_gpu0[dataset]\n            else:\n                future = futures_gpu1[dataset]\n    \n            # Wait for feature extraction\n            print(f\"waiting feature extraction at dataset = {dataset}\")\n            dataset_params = future.result()\n    \n            # if feature extraction is failed:\n            if dataset_params is None:\n                continue\n                \n            # Enqueue reconstruction\n            gc.collect()\n            futures_cpu[dataset] = executors.submit(\n                exec_reconstruction, dataset_params,\n            )\n\n        # Wait to reconstruction\n        for dataset, predictions in samples.items():\n            if datasets_to_process and (dataset not in datasets_to_process):\n                print(f'Skipping \"{dataset}\"')\n                continue\n            gc.collect()\n            dataset_params = futures_cpu[dataset].result()\n            samples[dataset] = dataset_params.predictions\n            dataset_logs.append( dataset_params )\n            gc.collect()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:00:59.76879Z","iopub.execute_input":"2025-04-17T13:00:59.768975Z","iopub.status.idle":"2025-04-17T13:02:03.861577Z","shell.execute_reply.started":"2025-04-17T13:00:59.768958Z","shell.execute_reply":"2025-04-17T13:02:03.860647Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Create submission","metadata":{}},{"cell_type":"code","source":"# Must Create a submission file.\n\narray_to_str = lambda array: ';'.join([f\"{x:.09f}\" for x in array])\nnone_to_str = lambda n: ';'.join(['nan'] * n)\n\nsubmission_file = '/kaggle/working/submission.csv'\nwith open(submission_file, 'w') as f:\n    if is_train:\n        f.write('dataset,scene,image,rotation_matrix,translation_vector\\n')\n        for dataset in samples:\n            for prediction in samples[dataset]:\n                cluster_name = 'outliers' if prediction.cluster_index is None else f'retry{prediction.retry_index}cluster{prediction.cluster_index}'\n                rotation = none_to_str(9) if prediction.rotation is None else array_to_str(prediction.rotation.flatten())\n                translation = none_to_str(3) if prediction.translation is None else array_to_str(prediction.translation)\n                f.write(f'{prediction.dataset},{cluster_name},{prediction.filename},{rotation},{translation}\\n')\n    else:\n        f.write('image_id,dataset,scene,image,rotation_matrix,translation_vector\\n')\n        for dataset in samples:\n            for prediction in samples[dataset]:\n                cluster_name = 'outliers' if prediction.cluster_index is None else f'retry{prediction.retry_index}cluster{prediction.cluster_index}'\n                rotation = none_to_str(9) if prediction.rotation is None else array_to_str(prediction.rotation.flatten())\n                translation = none_to_str(3) if prediction.translation is None else array_to_str(prediction.translation)\n                f.write(f'{prediction.image_id},{prediction.dataset},{cluster_name},{prediction.filename},{rotation},{translation}\\n')\n\n!head {submission_file}","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:02:03.862482Z","iopub.execute_input":"2025-04-17T13:02:03.862698Z","iopub.status.idle":"2025-04-17T13:02:04.148939Z","shell.execute_reply.started":"2025-04-17T13:02:03.862679Z","shell.execute_reply":"2025-04-17T13:02:04.147891Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Running logs","metadata":{}},{"cell_type":"code","source":"for dataset_params in dataset_logs:\n    print(\"=\" * 100)\n    dataset_params.summary()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:02:04.149913Z","iopub.execute_input":"2025-04-17T13:02:04.15018Z","iopub.status.idle":"2025-04-17T13:02:04.157583Z","shell.execute_reply.started":"2025-04-17T13:02:04.150156Z","shell.execute_reply":"2025-04-17T13:02:04.156729Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# CV scoring","metadata":{}},{"cell_type":"code","source":"# Definitely Compute results if running on the training set.\n# Do not do this when submitting a notebook for scoring. All you have to do is save your submission to /kaggle/working/submission.csv.\n\nif is_train:\n    t = time()\n    final_score, dataset_scores = metric.score(\n        gt_csv='/kaggle/input/image-matching-challenge-2025/train_labels.csv',\n        user_csv=submission_file,\n        thresholds_csv='/kaggle/input/image-matching-challenge-2025/train_thresholds.csv',\n        mask_csv=None if is_train else os.path.join(data_dir, 'mask.csv'),\n        inl_cf=0,\n        strict_cf=-1,\n        verbose=True,\n    )\n    print(f'Computed metric in: {time() - t:.02f} sec.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-04-17T13:02:04.158495Z","iopub.execute_input":"2025-04-17T13:02:04.158835Z","iopub.status.idle":"2025-04-17T13:02:04.171534Z","shell.execute_reply.started":"2025-04-17T13:02:04.158804Z","shell.execute_reply":"2025-04-17T13:02:04.170853Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#CELL 22: Visualization of Results\n\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import ConnectionPatch\nimport matplotlib.patches as mpatches\nfrom mpl_toolkits.mplot3d import Axes3D\nimport numpy as np\nimport cv2\nfrom PIL import Image\nimport h5py\nimport os\n\n# ============================================================================\n# 1. VISUALIZE IMAGE PAIRS WITH MATCHES\n# ============================================================================\n\ndef visualize_matches(img_path1, img_path2, matches_h5_path, dataset_name, \n                      feature_dir, max_matches=50, save_path=None):\n    \"\"\"\n    Visualize feature matches between two images\n    \n    Args:\n        img_path1, img_path2: Full paths to images\n        matches_h5_path: Path to matches.h5 file\n        dataset_name: Name of dataset\n        feature_dir: Feature directory path\n        max_matches: Maximum number of matches to display\n        save_path: If provided, saves figure to this path\n    \"\"\"\n    # Load images\n    img1 = cv2.imread(img_path1)\n    img2 = cv2.imread(img_path2)\n    img1 = cv2.cvtColor(img1, cv2.COLOR_BGR2RGB)\n    img2 = cv2.cvtColor(img2, cv2.COLOR_BGR2RGB)\n    \n    # Get filenames\n    fname1 = os.path.basename(img_path1)\n    fname2 = os.path.basename(img_path2)\n    \n    # Load keypoints and matches\n    kpts_h5 = h5py.File(os.path.join(feature_dir, 'keypoints.h5'), 'r')\n    matches_h5 = h5py.File(matches_h5_path, 'r')\n    \n    try:\n        kpts1 = kpts_h5[fname1][:]\n        kpts2 = kpts_h5[fname2][:]\n        match_indices = matches_h5[fname1][fname2][:]\n        \n        # Limit number of matches for visualization\n        if len(match_indices) > max_matches:\n            indices = np.random.choice(len(match_indices), max_matches, replace=False)\n            match_indices = match_indices[indices]\n        \n        # Create side-by-side visualization\n        h1, w1 = img1.shape[:2]\n        h2, w2 = img2.shape[:2]\n        max_h = max(h1, h2)\n        \n        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n        \n        ax1.imshow(img1)\n        ax2.imshow(img2)\n        \n        # Draw matches\n        colors = plt.cm.rainbow(np.linspace(0, 1, len(match_indices)))\n        \n        for idx, (m_idx1, m_idx2) in enumerate(match_indices):\n            kp1 = kpts1[int(m_idx1)]\n            kp2 = kpts2[int(m_idx2)]\n            \n            # Draw keypoints\n            ax1.plot(kp1[0], kp1[1], 'o', color=colors[idx], markersize=8)\n            ax2.plot(kp2[0], kp2[1], 'o', color=colors[idx], markersize=8)\n            \n            # Draw connection line\n            con = ConnectionPatch(\n                xyA=(kp1[0], kp1[1]), coordsA=ax1.transData,\n                xyB=(kp2[0], kp2[1]), coordsB=ax2.transData,\n                color=colors[idx], linewidth=1, alpha=0.5\n            )\n            fig.add_artist(con)\n        \n        ax1.set_title(f'{fname1}\\n{len(kpts1)} keypoints', fontsize=12)\n        ax2.set_title(f'{fname2}\\n{len(kpts2)} keypoints', fontsize=12)\n        ax1.axis('off')\n        ax2.axis('off')\n        \n        plt.suptitle(f'Dataset: {dataset_name}\\nMatches: {len(match_indices)}', \n                     fontsize=14, fontweight='bold')\n        plt.tight_layout()\n        \n        if save_path:\n            plt.savefig(save_path, dpi=150, bbox_inches='tight')\n            print(f\"Saved to {save_path}\")\n        \n        plt.show()\n        \n    except Exception as e:\n        print(f\"Error visualizing {fname1} - {fname2}: {e}\")\n    finally:\n        kpts_h5.close()\n        matches_h5.close()\n\n\n# ============================================================================\n# 2. VISUALIZE 3D RECONSTRUCTION (Point Cloud + Cameras)\n# ============================================================================\n\ndef visualize_3d_reconstruction(reconstruction_path, dataset_name, \n                                max_points=5000, save_path=None):\n    \"\"\"\n    Visualize 3D reconstruction with point cloud and camera positions\n    \n    Args:\n        reconstruction_path: Path to COLMAP reconstruction directory\n        dataset_name: Name of dataset\n        max_points: Maximum number of 3D points to display\n        save_path: If provided, saves figure to this path\n    \"\"\"\n    try:\n        rec = pycolmap.Reconstruction(reconstruction_path)\n        \n        # Extract 3D points\n        points_3d = []\n        colors = []\n        for point3D in rec.points3D.values():\n            points_3d.append(point3D.xyz)\n            colors.append(point3D.color / 255.0)  # Normalize to [0,1]\n        \n        points_3d = np.array(points_3d)\n        colors = np.array(colors)\n        \n        # Subsample if too many points\n        if len(points_3d) > max_points:\n            indices = np.random.choice(len(points_3d), max_points, replace=False)\n            points_3d = points_3d[indices]\n            colors = colors[indices]\n        \n        # Extract camera positions and orientations\n        camera_positions = []\n        camera_directions = []\n        for image in rec.images.values():\n            # Camera center in world coordinates\n            cam_center = -image.cam_from_world.rotation.matrix().T @ image.cam_from_world.translation\n            camera_positions.append(cam_center)\n            \n            # Camera viewing direction (z-axis of camera frame)\n            cam_direction = image.cam_from_world.rotation.matrix().T @ np.array([0, 0, 1])\n            camera_directions.append(cam_direction)\n        \n        camera_positions = np.array(camera_positions)\n        camera_directions = np.array(camera_directions)\n        \n        # Create 3D plot\n        fig = plt.figure(figsize=(16, 12))\n        ax = fig.add_subplot(111, projection='3d')\n        \n        # Plot 3D points\n        ax.scatter(points_3d[:, 0], points_3d[:, 1], points_3d[:, 2],\n                  c=colors, s=1, alpha=0.5, label='3D Points')\n        \n        # Plot camera positions\n        ax.scatter(camera_positions[:, 0], camera_positions[:, 1], camera_positions[:, 2],\n                  c='red', s=100, marker='^', label='Cameras', edgecolors='black', linewidth=1)\n        \n        # Plot camera viewing directions\n        for pos, direction in zip(camera_positions, camera_directions):\n            ax.quiver(pos[0], pos[1], pos[2],\n                     direction[0], direction[1], direction[2],\n                     length=1.0, color='blue', alpha=0.6, arrow_length_ratio=0.3)\n        \n        ax.set_xlabel('X', fontsize=12)\n        ax.set_ylabel('Y', fontsize=12)\n        ax.set_zlabel('Z', fontsize=12)\n        ax.set_title(f'3D Reconstruction: {dataset_name}\\n'\n                    f'{len(rec.images)} cameras, {len(rec.points3D)} 3D points',\n                    fontsize=14, fontweight='bold')\n        ax.legend(fontsize=10)\n        \n        # Equal aspect ratio\n        max_range = np.array([\n            points_3d[:, 0].max() - points_3d[:, 0].min(),\n            points_3d[:, 1].max() - points_3d[:, 1].min(),\n            points_3d[:, 2].max() - points_3d[:, 2].min()\n        ]).max() / 2.0\n        \n        mid_x = (points_3d[:, 0].max() + points_3d[:, 0].min()) * 0.5\n        mid_y = (points_3d[:, 1].max() + points_3d[:, 1].min()) * 0.5\n        mid_z = (points_3d[:, 2].max() + points_3d[:, 2].min()) * 0.5\n        \n        ax.set_xlim(mid_x - max_range, mid_x + max_range)\n        ax.set_ylim(mid_y - max_range, mid_y + max_range)\n        ax.set_zlim(mid_z - max_range, mid_z + max_range)\n        \n        if save_path:\n            plt.savefig(save_path, dpi=150, bbox_inches='tight')\n            print(f\"Saved to {save_path}\")\n        \n        plt.show()\n        \n    except Exception as e:\n        print(f\"Error visualizing reconstruction: {e}\")\n\n\n# ============================================================================\n# 3. VISUALIZE CAMERA POSES (Bird's Eye View)\n# ============================================================================\n\ndef visualize_camera_poses_2d(reconstruction_path, dataset_name, save_path=None):\n    \"\"\"\n    Visualize camera positions from bird's eye view (top-down)\n    \n    Args:\n        reconstruction_path: Path to COLMAP reconstruction directory\n        dataset_name: Name of dataset\n        save_path: If provided, saves figure to this path\n    \"\"\"\n    try:\n        rec = pycolmap.Reconstruction(reconstruction_path)\n        \n        camera_positions = []\n        camera_names = []\n        \n        for image in rec.images.values():\n            cam_center = -image.cam_from_world.rotation.matrix().T @ image.cam_from_world.translation\n            camera_positions.append(cam_center)\n            camera_names.append(image.name[:15])  # Truncate long names\n        \n        camera_positions = np.array(camera_positions)\n        \n        fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 8))\n        \n        # XY plane (top view)\n        ax1.scatter(camera_positions[:, 0], camera_positions[:, 1], \n                   c='red', s=100, alpha=0.6, edgecolors='black')\n        for i, name in enumerate(camera_names):\n            ax1.annotate(name, (camera_positions[i, 0], camera_positions[i, 1]),\n                        fontsize=8, alpha=0.7)\n        ax1.set_xlabel('X', fontsize=12)\n        ax1.set_ylabel('Y', fontsize=12)\n        ax1.set_title('Top View (XY plane)', fontsize=12)\n        ax1.grid(True, alpha=0.3)\n        ax1.set_aspect('equal')\n        \n        # XZ plane (side view)\n        ax2.scatter(camera_positions[:, 0], camera_positions[:, 2], \n                   c='blue', s=100, alpha=0.6, edgecolors='black')\n        for i, name in enumerate(camera_names):\n            ax2.annotate(name, (camera_positions[i, 0], camera_positions[i, 2]),\n                        fontsize=8, alpha=0.7)\n        ax2.set_xlabel('X', fontsize=12)\n        ax2.set_ylabel('Z', fontsize=12)\n        ax2.set_title('Side View (XZ plane)', fontsize=12)\n        ax2.grid(True, alpha=0.3)\n        ax2.set_aspect('equal')\n        \n        plt.suptitle(f'Camera Poses: {dataset_name}\\n{len(rec.images)} cameras',\n                    fontsize=14, fontweight='bold')\n        plt.tight_layout()\n        \n        if save_path:\n            plt.savefig(save_path, dpi=150, bbox_inches='tight')\n            print(f\"Saved to {save_path}\")\n        \n        plt.show()\n        \n    except Exception as e:\n        print(f\"Error visualizing camera poses: {e}\")\n\n\n# ============================================================================\n# 4. VISUALIZE MATCH STATISTICS\n# ============================================================================\n\ndef visualize_match_statistics(feature_dir, dataset_name, save_path=None):\n    \"\"\"\n    Visualize statistics about matches in the dataset\n    \n    Args:\n        feature_dir: Feature directory path\n        dataset_name: Name of dataset\n        save_path: If provided, saves figure to this path\n    \"\"\"\n    matches_path = os.path.join(feature_dir, 'matches.h5')\n    \n    try:\n        matches_h5 = h5py.File(matches_path, 'r')\n        \n        num_matches_list = []\n        pair_names = []\n        \n        for key1 in matches_h5.keys():\n            for key2 in matches_h5[key1].keys():\n                num_matches = len(matches_h5[key1][key2])\n                num_matches_list.append(num_matches)\n                pair_names.append(f\"{key1[:10]}-{key2[:10]}\")\n        \n        num_matches_list = np.array(num_matches_list)\n        \n        fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n        \n        # Histogram of matches\n        axes[0, 0].hist(num_matches_list, bins=50, color='skyblue', edgecolor='black')\n        axes[0, 0].set_xlabel('Number of Matches', fontsize=12)\n        axes[0, 0].set_ylabel('Frequency', fontsize=12)\n        axes[0, 0].set_title('Distribution of Match Counts', fontsize=12)\n        axes[0, 0].axvline(np.median(num_matches_list), color='red', \n                          linestyle='--', label=f'Median: {np.median(num_matches_list):.0f}')\n        axes[0, 0].legend()\n        axes[0, 0].grid(True, alpha=0.3)\n        \n        # Top pairs by match count\n        top_k = min(20, len(num_matches_list))\n        top_indices = np.argsort(num_matches_list)[-top_k:]\n        top_pairs = [pair_names[i] for i in top_indices]\n        top_counts = num_matches_list[top_indices]\n        \n        axes[0, 1].barh(range(top_k), top_counts, color='lightcoral')\n        axes[0, 1].set_yticks(range(top_k))\n        axes[0, 1].set_yticklabels(top_pairs, fontsize=8)\n        axes[0, 1].set_xlabel('Number of Matches', fontsize=12)\n        axes[0, 1].set_title(f'Top {top_k} Pairs by Match Count', fontsize=12)\n        axes[0, 1].grid(True, alpha=0.3, axis='x')\n        \n        # Cumulative distribution\n        sorted_matches = np.sort(num_matches_list)\n        cumulative = np.arange(1, len(sorted_matches) + 1) / len(sorted_matches)\n        axes[1, 0].plot(sorted_matches, cumulative, linewidth=2, color='green')\n        axes[1, 0].set_xlabel('Number of Matches', fontsize=12)\n        axes[1, 0].set_ylabel('Cumulative Probability', fontsize=12)\n        axes[1, 0].set_title('Cumulative Distribution', fontsize=12)\n        axes[1, 0].grid(True, alpha=0.3)\n        \n        # Statistics summary\n        stats_text = f\"\"\"\n        Statistics Summary:\n        \n        Total Pairs: {len(num_matches_list)}\n        \n        Mean Matches: {np.mean(num_matches_list):.1f}\n        Median Matches: {np.median(num_matches_list):.1f}\n        Std Dev: {np.std(num_matches_list):.1f}\n        \n        Min Matches: {np.min(num_matches_list):.0f}\n        Max Matches: {np.max(num_matches_list):.0f}\n        \n        25th Percentile: {np.percentile(num_matches_list, 25):.1f}\n        75th Percentile: {np.percentile(num_matches_list, 75):.1f}\n        \"\"\"\n        \n        axes[1, 1].text(0.1, 0.5, stats_text, fontsize=11, \n                       verticalalignment='center', family='monospace',\n                       bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))\n        axes[1, 1].axis('off')\n        \n        plt.suptitle(f'Match Statistics: {dataset_name}', \n                    fontsize=14, fontweight='bold')\n        plt.tight_layout()\n        \n        if save_path:\n            plt.savefig(save_path, dpi=150, bbox_inches='tight')\n            print(f\"Saved to {save_path}\")\n        \n        plt.show()\n        \n        matches_h5.close()\n        \n    except Exception as e:\n        print(f\"Error visualizing match statistics: {e}\")\n\n\n# ============================================================================\n# 5. VISUALIZE RECONSTRUCTION QUALITY (Reprojection Errors)\n# ============================================================================\n\ndef visualize_reconstruction_quality(reconstruction_path, dataset_name, save_path=None):\n    \"\"\"\n    Visualize reconstruction quality metrics\n    \n    Args:\n        reconstruction_path: Path to COLMAP reconstruction directory\n        dataset_name: Name of dataset\n        save_path: If provided, saves figure to this path\n    \"\"\"\n    try:\n        rec = pycolmap.Reconstruction(reconstruction_path)\n        \n        # Collect reprojection errors\n        reprojection_errors = []\n        track_lengths = []\n        \n        for point3D in rec.points3D.values():\n            reprojection_errors.append(point3D.error)\n            track_lengths.append(point3D.track.length())\n        \n        reprojection_errors = np.array(reprojection_errors)\n        track_lengths = np.array(track_lengths)\n        \n        fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n        \n        # Reprojection error histogram\n        axes[0, 0].hist(reprojection_errors, bins=50, color='salmon', edgecolor='black')\n        axes[0, 0].set_xlabel('Reprojection Error (pixels)', fontsize=12)\n        axes[0, 0].set_ylabel('Frequency', fontsize=12)\n        axes[0, 0].set_title('Distribution of Reprojection Errors', fontsize=12)\n        axes[0, 0].axvline(np.median(reprojection_errors), color='red', \n                          linestyle='--', label=f'Median: {np.median(reprojection_errors):.2f}px')\n        axes[0, 0].legend()\n        axes[0, 0].grid(True, alpha=0.3)\n        \n        # Track length histogram\n        axes[0, 1].hist(track_lengths, bins=50, color='lightgreen', edgecolor='black')\n        axes[0, 1].set_xlabel('Track Length (observations)', fontsize=12)\n        axes[0, 1].set_ylabel('Frequency', fontsize=12)\n        axes[0, 1].set_title('Distribution of Track Lengths', fontsize=12)\n        axes[0, 1].axvline(np.median(track_lengths), color='red', \n                          linestyle='--', label=f'Median: {np.median(track_lengths):.1f}')\n        axes[0, 1].legend()\n        axes[0, 1].grid(True, alpha=0.3)\n        \n        # Error vs Track Length scatter\n        axes[1, 0].scatter(track_lengths, reprojection_errors, alpha=0.3, s=10)\n        axes[1, 0].set_xlabel('Track Length', fontsize=12)\n        axes[1, 0].set_ylabel('Reprojection Error (pixels)', fontsize=12)\n        axes[1, 0].set_title('Error vs Track Length', fontsize=12)\n        axes[1, 0].grid(True, alpha=0.3)\n        \n        # Summary statistics\n        mean_error = rec.compute_mean_reprojection_error()\n        \n        stats_text = f\"\"\"\n        Quality Metrics:\n        \n        3D Points: {len(rec.points3D)}\n        Registered Images: {len(rec.images)}\n        \n        Mean Reprojection Error: {mean_error:.3f} px\n        Median Reprojection Error: {np.median(reprojection_errors):.3f} px\n        Std Reprojection Error: {np.std(reprojection_errors):.3f} px\n        \n        Mean Track Length: {np.mean(track_lengths):.2f}\n        Median Track Length: {np.median(track_lengths):.1f}\n        Max Track Length: {np.max(track_lengths):.0f}\n        \"\"\"\n        \n        axes[1, 1].text(0.1, 0.5, stats_text, fontsize=11,\n                       verticalalignment='center', family='monospace',\n                       bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))\n        axes[1, 1].axis('off')\n        \n        plt.suptitle(f'Reconstruction Quality: {dataset_name}',\n                    fontsize=14, fontweight='bold')\n        plt.tight_layout()\n        \n        if save_path:\n            plt.savefig(save_path, dpi=150, bbox_inches='tight')\n            print(f\"Saved to {save_path}\")\n        \n        plt.show()\n        \n    except Exception as e:\n        print(f\"Error visualizing reconstruction quality: {e}\")\n\n\n# ============================================================================\n# 6. MAIN VISUALIZATION FUNCTION - RUN THIS!\n# ============================================================================\n\ndef visualize_dataset_results(dataset_name, workdir='/tmp/kaggle/working/result/', \n                              data_dir='/kaggle/input/image-matching-challenge-2025',\n                              num_match_examples=3):\n    \"\"\"\n    Comprehensive visualization for a dataset\n    \n    Args:\n        dataset_name: Name of the dataset (e.g., 'ETs', 'stairs')\n        workdir: Working directory where results are stored\n        data_dir: Data directory\n        num_match_examples: Number of match pair examples to show\n    \"\"\"\n    print(f\"\\n{'='*80}\")\n    print(f\"VISUALIZING RESULTS FOR DATASET: {dataset_name}\")\n    print(f\"{'='*80}\\n\")\n    \n    feature_dir = os.path.join(workdir, 'featureout', dataset_name)\n    images_dir = os.path.join(data_dir, 'train', dataset_name)\n    \n    # Create output directory for visualizations\n    vis_dir = os.path.join(workdir, 'visualizations', dataset_name)\n    os.makedirs(vis_dir, exist_ok=True)\n    \n    # 1. Visualize match statistics\n    print(\"1. Generating match statistics...\")\n    visualize_match_statistics(\n        feature_dir, \n        dataset_name,\n        save_path=os.path.join(vis_dir, '1_match_statistics.png')\n    )\n    \n    # 2. Visualize some example matches\n    print(f\"\\n2. Generating {num_match_examples} example match visualizations...\")\n    matches_h5 = h5py.File(os.path.join(feature_dir, 'matches.h5'), 'r')\n    \n    # Get some good matches to visualize\n    match_pairs = []\n    for key1 in list(matches_h5.keys())[:10]:  # Check first 10 images\n        for key2 in list(matches_h5[key1].keys())[:3]:  # Get 3 pairs per image\n            num_matches = len(matches_h5[key1][key2])\n            match_pairs.append((key1, key2, num_matches))\n    \n    matches_h5.close()\n    \n    # Sort by number of matches and take top examples\n    match_pairs.sort(key=lambda x: x[2], reverse=True)\n    match_pairs = match_pairs[:num_match_examples]\n    \n    for idx, (fname1, fname2, num_matches) in enumerate(match_pairs):\n        print(f\"   Example {idx+1}: {fname1} - {fname2} ({num_matches} matches)\")\n        visualize_matches(\n            os.path.join(images_dir, fname1),\n            os.path.join(images_dir, fname2),\n            os.path.join(feature_dir, 'matches.h5'),\n            dataset_name,\n            feature_dir,\n            max_matches=100,\n            save_path=os.path.join(vis_dir, f'2_matches_example_{idx+1}.png')\n        )\n    \n    # 3. Find and visualize reconstructions\n    print(\"\\n3. Looking for reconstructions...\")\n    \n    # Check for merged reconstructions first\n    merged_rec_dirs = glob.glob(os.path.join(feature_dir, 'colmap_rec_merged', 'scene_*'))\n    \n    if not merged_rec_dirs:\n        # Check for individual reconstructions\n        rec_dirs = glob.glob(os.path.join(feature_dir, 'colmap_rec_*', '*'))\n        rec_dirs = [d for d in rec_dirs if os.path.isfile(os.path.join(d, 'cameras.bin'))]\n    else:\n        rec_dirs = merged_rec_dirs\n    \n    if rec_dirs:\n        # Sort by size (number of images)\n        rec_sizes = []\n        for rec_dir in rec_dirs:\n            try:\n                rec = pycolmap.Reconstruction(rec_dir)\n                rec_sizes.append(len(rec.images))\n            except:\n                rec_sizes.append(0)\n        \n        # Take largest reconstruction\n        largest_idx = np.argmax(rec_sizes)\n        largest_rec_dir = rec_dirs[largest_idx]\n        \n        print(f\"   Found {len(rec_dirs)} reconstructions\")\n        print(f\"   Visualizing largest: {largest_rec_dir}\")\n        print(f\"   ({rec_sizes[largest_idx]} images)\")\n        \n        # 3D reconstruction\n        print(\"\\n   a) Generating 3D reconstruction visualization...\")\n        visualize_3d_reconstruction(\n            largest_rec_dir,\n            dataset_name,\n            max_points=10000,\n            save_path=os.path.join(vis_dir, '3_reconstruction_3d.png')\n        )\n        \n        # Camera poses\n        print(\"   b) Generating camera pose visualization...\")\n        visualize_camera_poses_2d(\n            largest_rec_dir,\n            dataset_name,\n            save_path=os.path.join(vis_dir, '4_camera_poses.png')\n        )\n        \n        # Quality metrics\n        print(\"   c) Generating quality metrics visualization...\")\n        visualize_reconstruction_quality(\n            largest_rec_dir,\n            dataset_name,\n            save_path=os.path.join(vis_dir, '5_reconstruction_quality.png')\n        )\n    else:\n        print(\"   No reconstructions found!\")\n    \n    print(f\"\\n{'='*80}\")\n    print(f\"VISUALIZATION COMPLETE!\")\n    print(f\"Results saved to: {vis_dir}\")\n    print(f\"{'='*80}\\n\")\n\n\n# ============================================================================\n# USAGE EXAMPLES\n# ============================================================================\n\n# Example 1: Visualize a single dataset\nprint(\"To visualize results for a dataset, run:\")\nprint(\"visualize_dataset_results('ETs')\")\nprint(\"\\nOr for another dataset:\")\nprint(\"visualize_dataset_results('stairs')\")\n\n# Example 2: Visualize all processed datasets\nprint(\"\\n\" + \"=\"*80)\nprint(\"To visualize ALL processed datasets, run:\")\nprint(\"=\"*80)\nprint(\"\"\"\nfor dataset_name in samples.keys():\n    try:\n        visualize_dataset_results(dataset_name)\n    except Exception as e:\n        print(f\"Failed to visualize {dataset_name}: {e}\")\n\"\"\")\n\n# Example 3: Quick visualization for specific dataset\nprint(\"\\n\" + \"=\"*80)\nprint(\"QUICK START - Visualize your first dataset:\")\nprint(\"=\"*80)\nvisualize_dataset_results('ETs', num_match_examples=2","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}