{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":91498,"databundleVersionId":11655853,"isSourceIdPinned":false,"sourceType":"competition"},{"sourceId":11924468,"sourceType":"datasetVersion","datasetId":6988459},{"sourceId":4537,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":3329,"modelId":986}],"dockerImageVersionId":31041,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Import Libraries","metadata":{}},{"cell_type":"code","source":"!pip install faiss-cpu -q\n!pip install git+https://github.com/cvg/LightGlue.git -q\n!pip install pycolmap -q","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:17:26.497883Z","iopub.execute_input":"2025-06-23T07:17:26.498371Z","iopub.status.idle":"2025-06-23T07:19:04.706755Z","shell.execute_reply.started":"2025-06-23T07:17:26.498345Z","shell.execute_reply":"2025-06-23T07:19:04.705970Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os\nimport matplotlib.pyplot as plt\nfrom collections import defaultdict\nfrom tqdm import tqdm\nfrom itertools import combinations\nimport h5py\nimport dataclasses\n\nimport cv2\nfrom PIL import Image\n\nimport torch\nfrom transformers import AutoImageProcessor, AutoModel\nfrom torch.utils.data import Dataset, DataLoader\nfrom PIL import Image\nfrom torchvision import transforms\nfrom typing import Optional, Tuple\nimport torch.nn.functional as F\n\nimport faiss\nimport networkx as nx\nfrom community import community_louvain\n\nfrom lightglue import LightGlue, SuperPoint, ALIKED\nfrom lightglue.utils import rbd\nfrom lightglue import viz2d\nimport pycolmap\nfrom torchvision import transforms as T\n\nimport sys, warnings\nsys.path.append(\"/kaggle/input/imc25-utils\")\n\nfrom database import *\nfrom h5_to_db import *\nimport metric\nimport shutil\n\nfrom sklearn.manifold import TSNE\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom typing import Callable, List, Optional, Tuple, Union\nfrom pathlib import Path\n\n# import open3d as o3d\n# from open3d import JVisualizer\nfrom mpl_toolkits.mplot3d import Axes3D\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:04.708591Z","iopub.execute_input":"2025-06-23T07:19:04.708816Z","iopub.status.idle":"2025-06-23T07:19:29.903121Z","shell.execute_reply.started":"2025-06-23T07:19:04.708793Z","shell.execute_reply":"2025-06-23T07:19:29.902555Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Proposed Pipeline\n\n### Step 1: Global Clustering with DINO + FAISS + Louvain\n\n* Used DINOv2 Vision Transformers to extract global descriptors.\n* Construct a similarity graph via FAISS (cosine similarity).\n* Apply Louvain community detection to segment images into scene-level clusters.\n\n\n### Step 2: Local Matching with SuperPoint + LightGlue\n\n* Within each cluster, extract keypoints using SuperPoint.\n* Match image pairs using LightGlue with geometric filtering.\n* Store keypoints and inliers in HDF5 format for downstream Structure from Motion(SfM).\n\n\n### Step 3: 3D Reconstruction with COLMAP + PyCOLMAP\n* For each cluster, import images and feature matches into COLMAP.\n* Run incremental SfM to recover camera poses and sparse 3D points.\n* Run incremental SfM to reconstruct poses and sparse 3D points.\n* Reproject 3D points to generate sparse depth maps per image.\n* Store metadata in a Predictions dataclass: pose, intrinsics, depth maps, and 2D projections.\n\n\n### Step 4: Visualization \n* Visualize SuperPoint keypoints and LightGlue matches per image pair.\n* Overlay sparse depth maps on input images to highlight scene geometry.\n* Plot epipolar lines using predicted camera intrinsics and relative pose (from Step 3).  ","metadata":{}},{"cell_type":"markdown","source":"# STEP 1","metadata":{}},{"cell_type":"code","source":"# to extract  global feature using DINO Vision Transformers (ViTs)\nall_features = []\nimage_paths = []\n\ndef global_feature_extractor(\n    device, \n    Dataloader,\n    verbose = True,\n    pooling = 'mac',\n    model_path = '/kaggle/input/dinov2/pytorch/giant/1'):\n    \n    assert pooling in [\"mac\", \"cls\"], \"Pooling must be 'mac' or 'cls'\"\n    loop = tqdm(Dataloader, desc=\"Extracting DINO features\", disable=not verbose)\n    \n    processor = AutoImageProcessor.from_pretrained(model_path,use_fast=True)\n    dino_model = AutoModel.from_pretrained(model_path)\n    dino_model = dino_model.eval().to(device)\n    \n    with torch.inference_mode():\n        feature_list = []\n        for images, metadata in loop:\n            inputs = processor(images=images, return_tensors=\"pt\",do_rescale=False).to(device)\n            outputs = dino_model(**inputs)\n            if pooling == 'mac':\n                # last_hidden_state shape: (batch_size, sequence_length, hidden_dim)\n                vec = F.normalize(outputs.last_hidden_state[:,1:].max(dim=1)[0], dim=1, p=2) #all\n            else:\n                vec = F.normalize(outputs.last_hidden_state[:,0], dim=1, p=2) # batch, hidden_dim\n            feature_list.append(vec.detach().cpu())\n            image_paths.extend(metadata['image_path'])\n            \n        all_feature_list = torch.cat(feature_list, dim=0)\n        \n    return all_feature_list.numpy(), image_paths","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.903788Z","iopub.execute_input":"2025-06-23T07:19:29.904259Z","iopub.status.idle":"2025-06-23T07:19:29.910690Z","shell.execute_reply.started":"2025-06-23T07:19:29.904240Z","shell.execute_reply":"2025-06-23T07:19:29.909986Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# to build FAISS knn graph from global features\ndef build_knn_graph(features, data_dir, img_paths, k=10):\n    \"\"\"\n    Build a kNN graph using cosine similarity and per-scene threshold filtering.\n    \n    features: np.ndarray of global descriptors\n    img_paths: List of paths aligned with features\n    scene_thresholds: Dict[scene_name] -> float [mid-value]\n    k: number of neighbors\n\n    Returns:\n        List of (i, j, sim) edges that pass scene threshold check\n    \"\"\"\n    def load_threshold():\n        \"\"\"\n        return mid threshold value per scene in a Dict\n        \"\"\"\n        threshold = os.path.join(data_dir, \"train_thresholds.csv\")\n        thresh_df = pd.read_csv(threshold)\n        tdict = dict()\n        for _, row in thresh_df.iterrows():\n            tdict[row[\"scene\"]] = float(row['thresholds'].split(';')[-1])\n            \n        return tdict\n    \n    \n    features = features.astype(np.float32)\n    features /= np.linalg.norm(features, axis=1, keepdims=True)\n\n    index = faiss.IndexFlatIP(features.shape[1])  # inner product = cosine sim if normalized\n    index.add(features)\n    distances, indices = index.search(features, k + 1)  # +1 because self-match at index 0\n\n    paths = [os.path.dirname(p).split(\"/\")[-1] for p in img_paths]\n    og_scene_threshold = load_threshold()\n\n\n    edges = []\n    num_points = features.shape[0]\n    for i in range(num_points):\n        for j in range(1, k + 1):  # skip self-match\n            neighbor = indices[i][j]\n            sim = distances[i][j]\n\n            scene_1 = paths[i]\n            scene_2 = paths[neighbor]\n\n            if scene_1 == scene_2:\n                threshold = og_scene_threshold.get(scene_1, 0.1)\n                if sim >= threshold:\n                    edges.append((i, neighbor, sim))\n                else:\n                    continue\n\n    return edges","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.912456Z","iopub.execute_input":"2025-06-23T07:19:29.912802Z","iopub.status.idle":"2025-06-23T07:19:29.941658Z","shell.execute_reply.started":"2025-06-23T07:19:29.912768Z","shell.execute_reply":"2025-06-23T07:19:29.941044Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def build_graph_from_edges(edges):\n    \"\"\"\n    Builds a weighted undirected NetworkX graph from edge list (i, j, weight)\n    \"\"\"\n    G = nx.Graph()\n    for i, j, sim in edges:\n        G.add_edge(i, j, weight=sim)\n    return G\n\ndef louvian_cluster(graph):\n    return community_louvain.best_partition(graph,weight = 'weight')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.942296Z","iopub.execute_input":"2025-06-23T07:19:29.942521Z","iopub.status.idle":"2025-06-23T07:19:29.961060Z","shell.execute_reply.started":"2025-06-23T07:19:29.942498Z","shell.execute_reply":"2025-06-23T07:19:29.960467Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# STEP 2","metadata":{}},{"cell_type":"code","source":"# superpoint + lightglue -> match_pair [reformat]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.961832Z","iopub.execute_input":"2025-06-23T07:19:29.962135Z","iopub.status.idle":"2025-06-23T07:19:29.979905Z","shell.execute_reply.started":"2025-06-23T07:19:29.962115Z","shell.execute_reply":"2025-06-23T07:19:29.979353Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# img_load -> from lightglue util.py\n\ndef read_image(path, grayscale: bool = False) -> np.ndarray:\n    \"\"\"Read an image from path as RGB or grayscale\"\"\"\n    if not Path(path).exists():\n        raise FileNotFoundError(f\"No image at path {path}.\")\n    mode = cv2.IMREAD_GRAYSCALE if grayscale else cv2.IMREAD_COLOR\n    image = cv2.imread(str(path), mode)\n    if image is None:\n        raise IOError(f\"Could not read image at {path}.\")\n    if not grayscale:\n        image = image[..., ::-1]\n    return image\n\ndef resize_image(\n    image: np.ndarray,\n    size: Union[List[int], int],\n    fn: str = \"max\",\n    interp: Optional[str] = \"area\",\n) -> np.ndarray:\n    \"\"\"Resize an image to a fixed size, or according to max or min edge.\"\"\"\n    h, w = image.shape[:2]\n\n    fn = {\"max\": max, \"min\": min}[fn]\n    if isinstance(size, int):\n        scale = size / fn(h, w)\n        h_new, w_new = int(round(h * scale)), int(round(w * scale))\n        scale = (w_new / w, h_new / h)\n    elif isinstance(size, (tuple, list)):\n        h_new, w_new = size\n        scale = (w_new / w, h_new / h)\n    else:\n        raise ValueError(f\"Incorrect new size: {size}\")\n    mode = {\n        \"linear\": cv2.INTER_LINEAR,\n        \"cubic\": cv2.INTER_CUBIC,\n        \"nearest\": cv2.INTER_NEAREST,\n        \"area\": cv2.INTER_AREA,\n    }[interp]\n    return cv2.resize(image, (w_new, h_new), interpolation=mode), scale\n\n\ndef numpy_image_to_torch(image: np.ndarray) -> torch.Tensor:\n    \"\"\"Normalize the image tensor and reorder the dimensions.\"\"\"\n    if image.ndim == 3:\n        image = image.transpose((2, 0, 1))  # HxWxC to CxHxW\n    elif image.ndim == 2:\n        image = image[None]  # add channel axis\n    else:\n        raise ValueError(f\"Not an image: {image.shape}\")\n    return torch.tensor(image / 255.0, dtype=torch.float)\n\ndef load_image(path , resize: int = None, **kwargs) -> torch.Tensor:\n    image = read_image(path)\n    if resize is not None:\n        image, _ = resize_image(image, resize, **kwargs)\n    return numpy_image_to_torch(image)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.980596Z","iopub.execute_input":"2025-06-23T07:19:29.980823Z","iopub.status.idle":"2025-06-23T07:19:29.995827Z","shell.execute_reply.started":"2025-06-23T07:19:29.980800Z","shell.execute_reply":"2025-06-23T07:19:29.995100Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def match_pair_lightglue_superpoint(\n    image_fnames,\n    image_pairs,\n    device,\n    min_matches=20\n):\n    os.makedirs(feature_dir, exist_ok=True)\n    \n    with h5py.File(f'{feature_dir}/keypoints.h5', 'a') as f_kp, \\\n        h5py.File(f'{feature_dir}/descriptors.h5', 'a') as f_desc, \\\n        h5py.File(f'{feature_dir}/matches.h5', 'a') as f_match:\n        \n        extractor = SuperPoint(max_num_keypoints=4096).eval().to(device)\n        matcher = LightGlue(features='superpoint').eval().to(device)\n\n        for path1, path2 in tqdm(image_pairs, desc=\"Extracting & Matching\"):\n            key1 = os.path.basename(path1)\n            key2 = os.path.basename(path2)\n\n              \n            keyA, keyB = sorted([key1, key2])\n            pathA, pathB = (path1, path2) if key1 == keyA else (path2, path1)\n\n            if keyA in f_match and keyB in f_match[keyA]:\n                continue\n\n            try:\n                imgA = load_image(pathA)\n                imgB = load_image(pathB)\n\n                featsA = extractor.extract(imgA.to(device))\n                featsB = extractor.extract(imgB.to(device))\n\n                matchesAB = matcher({\"image0\": featsA, \"image1\": featsB})\n                featsA, featsB, matchesAB = [\n                    rbd(x) for x in [featsA, featsB, matchesAB]]# remove batch_dimension\n\n                matches = matchesAB[\"matches\"]\n\n                if matches.shape[0] < min_matches:\n                    continue\n\n                f_match.require_group(keyA).create_dataset(keyB, data=matches.cpu().numpy().astype(np.int32))\n\n                if keyA not in f_kp:\n                    f_kp.create_dataset(keyA, data=featsA[\"keypoints\"].cpu().numpy())\n                    f_desc.create_dataset(keyA, data=featsA[\"descriptors\"].cpu().numpy())\n\n                if keyB not in f_kp:\n                    f_kp.create_dataset(keyB, data=featsB[\"keypoints\"].cpu().numpy())\n                    f_desc.create_dataset(keyB, data=featsB[\"descriptors\"].cpu().numpy())\n\n            except Exception as e:\n                print(f\"[ERROR] Matching failed for {key1} ↔ {key2}: {e}\")\n                continue\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:29.996591Z","iopub.execute_input":"2025-06-23T07:19:29.996844Z","iopub.status.idle":"2025-06-23T07:19:30.272460Z","shell.execute_reply.started":"2025-06-23T07:19:29.996819Z","shell.execute_reply":"2025-06-23T07:19:30.271641Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# STEP 3","metadata":{}},{"cell_type":"code","source":"# COLMAP & pycolmap","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.273282Z","iopub.execute_input":"2025-06-23T07:19:30.273566Z","iopub.status.idle":"2025-06-23T07:19:30.290592Z","shell.execute_reply.started":"2025-06-23T07:19:30.273544Z","shell.execute_reply":"2025-06-23T07:19:30.289829Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def run_colmap_reconstruction(\n    cluster_imgs,\n    feature_dir,\n    cluster_idx,\n    output_root\n):\n    \"\"\"\n    Performs COLMAP-based SfM for a given cluster.\n    Steps:\n      - Creates a working directory\n      - Imports images, keypoints, matches into COLMAP\n      - Runs incremental mapping\n      - Returns pycolmap Reconstruction class objects \n    \"\"\"\n\n    os.makedirs(output_root, exist_ok=True)\n\n    db_path = os.path.join(output_root, \"colmap.db\")\n    if os.path.exists(db_path):\n        os.remove(db_path)\n\n    base_cluster_idic = {os.path.basename(p): p for p in cluster_imgs}\n    \n\n    def add_keypoints_flat(db, feature_dir, img_root, img_ext, camera_model, single_camera):\n        f_kp = h5py.File(os.path.join(feature_dir, \"keypoints.h5\"), \"r\")\n        fname_to_id = {}\n        camera_id = None\n        \n        for img_name in tqdm(list(f_kp.keys()), desc=\"Add_Keypoints\"):\n            if img_name in base_cluster_idic.keys():\n                full_img_path = base_cluster_idic.get(img_name)\n                scene_img_path = os.path.relpath(full_img_path, img_root)\n                \n                if not os.path.isfile(full_img_path):\n                    continue\n            \n                keypoints = f_kp[img_name][()]\n                \n                if len(keypoints) == 0:\n                    print(f\"for image {img_name} 0 keypoints\")\n                    continue\n                \n                if camera_id is None or not single_camera:\n                    camera_id = create_camera(db, full_img_path, camera_model)\n                    \n                image_id = db.add_image(scene_img_path, camera_id)\n                db.add_keypoints(image_id, keypoints)\n                fname_to_id[img_name] = image_id\n    \n        return fname_to_id\n\n\n    def add_matches_flat(db, feature_dir, fname_to_id):\n\n        f_match = h5py.File(os.path.join(feature_dir, \"matches.h5\"), \"r\")\n        added = set()\n        n_total  = 0\n        valid_pairs = 0\n    \n        # Estimate total number of match entries for tqdm\n        n_keys = len(f_match.keys())\n        n_total = (n_keys * (n_keys - 1)) // 2\n\n\n        with tqdm(total=n_total, desc=\"Importing Matches\") as pbar:         \n            for key1 in f_match.keys():\n                key_group = f_match[key1] \n                for key2 in key_group.keys():\n                    \n                    if key1 not in fname_to_id and key2 not in fname_to_id:\n                        continue\n\n                    id_1 = fname_to_id[key1]\n                    id_2 = fname_to_id[key2]\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                        pbar.update(1)\n                        continue\n                        \n                    matches = key_group[key2][()] \n                                        \n                    matches = np.unique(matches.astype(np.uint32), axis=0)\n                    if matches.shape[0] == 0 or matches.ndim != 2 or matches.shape[1]!=2:\n                        pbar.update(1)\n                        continue\n                    db.add_matches(id_1, id_2, matches.astype(np.uint32))\n                    db.add_two_view_geometry(id_1, id_2, matches)\n                    added.add(pair_id)\n                    pbar.update(1)\n                    valid_pairs += 1\n        \n        print(f\"[DEBUG] Total valid pairs added to DB: {valid_pairs}\")\n                        \n\n    \n    def import_into_colmap(img_list, feature_dir, database_path):\n        \n        db = COLMAPDatabase.connect(database_path)\n        db.create_tables()\n        \n        fname_to_id = add_keypoints_flat(\n            db,\n            feature_dir,\n            img_root=\"/kaggle/input/image-matching-challenge-2025/train\",\n            img_ext=\".png\", camera_model='pinhole', single_camera=False)\n    \n        add_matches_flat(\n            db,\n            feature_dir,\n            fname_to_id\n        )\n    \n        db.commit()\n        return\n    \n    import_into_colmap(cluster_imgs, feature_dir, db_path)\n    image_dir = \"/kaggle/input/image-matching-challenge-2025/train\"\n\n    # RANSAC - pycolmap ver\n    # pycolmap.match_exhaustive(db_path)\n        \n    \n    maps = dict()\n    mapper_options = pycolmap.IncrementalPipelineOptions()\n    mapper_options.num_threads = 2\n    \n    scene_size = len(cluster_imgs)\n    if scene_size <= 25:\n        mapper_options.min_model_size = 3\n        mapper_options.ba_refine_focal_length = True\n        mapper_options.ba_refine_principal_point = False\n        mapper_options.ba_refine_extra_params = False\n        mapper_options.ba_local_num_images = 4\n        mapper_options.ba_local_max_num_iterations = 15\n        \n    else:\n        mapper_options.min_model_size = 8\n        mapper_options.max_num_models = 25\n        mapper_options.ba_refine_focal_length = True\n        mapper_options.ba_refine_principal_point = True\n        mapper_options.ba_refine_extra_params = True\n        mapper_options.ba_local_num_images = 8\n        mapper_options.ba_global_max_num_iterations = 80\n        mapper_options.ba_local_max_num_iterations = 25\n        mapper_options.ba_global_function_tolerance = 1e-5\n        mapper_options.ba_local_function_tolerance = 1e-5\n\n    \n    maps = pycolmap.incremental_mapping(\n        database_path=db_path,\n        image_path=image_dir,\n        output_path=output_root,\n        options=mapper_options\n    )\n\n\n    return maps","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.292933Z","iopub.execute_input":"2025-06-23T07:19:30.293616Z","iopub.status.idle":"2025-06-23T07:19:30.308577Z","shell.execute_reply.started":"2025-06-23T07:19:30.293596Z","shell.execute_reply":"2025-06-23T07:19:30.308005Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# STEP 4","metadata":{}},{"cell_type":"markdown","source":"## keypoints and matches visualisation\n","metadata":{}},{"cell_type":"code","source":"def keypoint_match_visualiser(pose_img1, pose_img2):\n             \n    image0 = load_image(os.path.join(image_dir, \"train\", pose_img1))\n    image1 = load_image(os.path.join(image_dir, \"train\", pose_img2))\n\n    extractor = SuperPoint(max_num_keypoints=4096).eval().to(device)\n    matcher = LightGlue(features='superpoint').eval().to(device)\n    \n    feats0 = extractor.extract(image0.to(device))\n    feats1 = extractor.extract(image1.to(device))\n    matches01 = matcher({\"image0\": feats0, \"image1\": feats1})\n    feats0, feats1, matches01 = [\n        rbd(x) for x in [feats0, feats1, matches01]\n    ]  \n    \n    kpts0, kpts1, matches = feats0[\"keypoints\"], feats1[\"keypoints\"], matches01[\"matches\"]\n    m_kpts0, m_kpts1 = kpts0[matches[..., 0]], kpts1[matches[..., 1]]\n    \n    axes = viz2d.plot_images([image0, image1])\n    viz2d.plot_matches(m_kpts0, m_kpts1, color=\"lime\", lw=0.2)\n    viz2d.add_text(0, f'Stop after {matches01[\"stop\"]} layers', fs=20)\n    \n    kpc0, kpc1 = viz2d.cm_prune(matches01[\"prune0\"]), viz2d.cm_prune(matches01[\"prune1\"])\n    viz2d.plot_images([image0, image1])\n    viz2d.plot_keypoints([kpts0, kpts1], colors=[kpc0, kpc1], ps=10)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.309347Z","iopub.execute_input":"2025-06-23T07:19:30.309601Z","iopub.status.idle":"2025-06-23T07:19:30.328000Z","shell.execute_reply.started":"2025-06-23T07:19:30.309578Z","shell.execute_reply":"2025-06-23T07:19:30.327230Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Depth Map & Points\n","metadata":{}},{"cell_type":"code","source":"def sparse_depth_map(image, camera, points3D, image_shape):\n    depth_map = np.zeros(image_shape, dtype=np.float32)\n    projected = []\n\n    K = camera.calibration_matrix()\n    R = np.asarray(image.cam_from_world.rotation.matrix()).reshape(3, 3)\n    t = np.asarray(image.cam_from_world.translation).reshape(3)\n\n    for pt in points3D.values():\n        image_ids_in_track = [te.image_id for te in pt.track.elements]\n        if image.image_id not in image_ids_in_track:\n            continue\n\n        X = np.asarray(pt.xyz).reshape(3)\n        X_cam = R @ X + t\n\n        if X_cam[2] <= 0:\n            continue\n\n        x_proj = K @ X_cam\n        u = int(round(x_proj[0] / x_proj[2]))\n        v = int(round(x_proj[1] / x_proj[2]))\n        \n        if 0 <= u < image_shape[1] and 0 <= v < image_shape[0]:\n            depth = X_cam[2]            \n            depth_map[v,u] = depth\n            projected.append([u, v, depth])\n            \n    return depth_map, np.array(projected)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.328734Z","iopub.execute_input":"2025-06-23T07:19:30.328986Z","iopub.status.idle":"2025-06-23T07:19:30.346661Z","shell.execute_reply.started":"2025-06-23T07:19:30.328963Z","shell.execute_reply":"2025-06-23T07:19:30.346070Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Epipolar on Depth","metadata":{}},{"cell_type":"code","source":"def compute_fundamental_matrix(K_A, R_A, t_A, K_B, R_B, t_B):\n    R_rel = R_B @ R_A.T\n    t_rel = t_B - R_rel @ t_A\n    t_x = np.array([\n        [0, -t_rel[2], t_rel[1]],\n        [t_rel[2], 0, -t_rel[0]],\n        [-t_rel[1], t_rel[0], 0]\n    ])\n    E = t_x @ R_rel\n    F = np.linalg.inv(K_B).T @ E @ np.linalg.inv(K_A)\n    return F\n\n\ndef epipolar(model, predA, predB, max_points):\n    imgA = predA.datapath\n    imgB = predB.datapath\n    camA = model.cameras[predA.camera_id]\n    camB = model.cameras[predB.camera_id]\n\n    K_A = camA.calibration_matrix()\n    K_B = camB.calibration_matrix()\n    R_A, t_A = predA.rotation, predA.translation\n    R_B, t_B = predB.rotation, predB.translation\n\n    F = compute_fundamental_matrix(K_A, R_A, t_A, K_B, R_B, t_B)\n\n    ptsA, ptsB = [], []\n    for pt3D in model.points3D.values():\n        track = {e.image_id: e.point2D_idx for e in pt3D.track.elements}\n        if predA.image_id in track and predB.image_id in track:\n            ptA = model.images[predA.image_id].points2D[track[predA.image_id]].xy\n            ptB = model.images[predB.image_id].points2D[track[predB.image_id]].xy\n            ptsA.append(ptA)\n            ptsB.append(ptB)\n\n    if len(ptsA) < 8:\n        print(\"Insufficient shared 3D points.\")\n        return\n\n    ptsA, ptsB = np.array(ptsA[:max_points]), np.array(ptsB[:max_points])\n    linesA = cv2.computeCorrespondEpilines(ptsB.reshape(-1, 1, 2), 2, F).reshape(-1, 3)\n    linesB = cv2.computeCorrespondEpilines(ptsA.reshape(-1, 1, 2), 1, F).reshape(-1, 3)\n\n    imgA_rgb = np.array(Image.open(imgA).convert(\"RGB\"))\n    imgB_rgb = np.array(Image.open(imgB).convert(\"RGB\"))\n    imgA_draw = imgA_rgb.copy()\n    imgB_draw = imgB_rgb.copy()\n\n    for r, pt in zip(linesA, ptsA):\n        color = tuple(np.random.randint(64, 255, 3).tolist())\n        x0, y0 = 0, int(-r[2]/r[1])\n        x1, y1 = imgA_draw.shape[1], int(-(r[2] + r[0]*x0) / r[1])\n        cv2.line(imgA_draw, (x0, y0), (x1, y1), color, 7)\n        cv2.circle(imgA_draw, tuple(pt.astype(int)), 6, color, thickness=15)\n\n    for r, pt in zip(linesB, ptsB):\n        color = tuple(np.random.randint(64, 255, 3).tolist())\n        x0, y0 = 0, int(-r[2]/r[1])\n        x1, y1 = imgB_draw.shape[1], int(-(r[2] + r[0]*x0) / r[1])\n        cv2.line(imgB_draw, (x0, y0), (x1, y1), color, 7)\n        cv2.circle(imgB_draw, tuple(pt.astype(int)), 6, color, thickness=15)\n\n    fig, axs = plt.subplots(1, 2, figsize=(16, 8))\n    axs[0].imshow(imgA_draw)\n    axs[0].set_title(f\"Epipolar Lines on {predA.filename}\")\n    axs[0].axis(\"off\")\n    axs[1].imshow(imgB_draw)\n    axs[1].set_title(f\"Epipolar Lines on {predB.filename}\")\n    axs[1].axis(\"off\")\n    plt.tight_layout()\n    plt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:55:34.611700Z","iopub.execute_input":"2025-06-23T08:55:34.612011Z","iopub.status.idle":"2025-06-23T08:55:34.626027Z","shell.execute_reply.started":"2025-06-23T08:55:34.611993Z","shell.execute_reply":"2025-06-23T08:55:34.625419Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# DATALOADER","metadata":{}},{"cell_type":"code","source":"image_transform = transforms.Compose([\n    transforms.Resize((224, 224)),\n    transforms.ToTensor(),\n])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.368006Z","iopub.execute_input":"2025-06-23T07:19:30.368470Z","iopub.status.idle":"2025-06-23T07:19:30.385741Z","shell.execute_reply.started":"2025-06-23T07:19:30.368447Z","shell.execute_reply":"2025-06-23T07:19:30.384966Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class ImageMatching():\n    def __init__(\n        self,\n        csv_path: str,\n        image_dir: str,\n        to_train: bool = True,\n        datasets_filter: Optional[list] = None,\n        image_transform: Optional[transforms.Compose] = None\n        ):\n        self.data = []\n        self.is_train = to_train\n        self.data_dir = image_dir\n        self.image_transform = image_transform \n        df = pd.read_csv(csv_path)\n\n        for _, row in df.iterrows():\n            if row['dataset'] in datasets_filter:\n                entry = {\n                    'dataset': row['dataset'],\n                    'scene': row['scene'],\n                    'filename': row['image'],\n                    'image_id': row.get('image_id', None),  \n                    'rotation': row.get('rotation_matrix', None) if to_train else None,\n                    'translation': row.get('translation_vector', None) if to_train else None,\n                }\n                self.data.append(entry)\n            \n\n    def __len__(self):\n        return len(self.data)\n\n    def __getitem__(self, idx:int):\n        entry = self.data[idx]\n        base_dir = 'train' if self.is_train else 'test'\n        image_path = os.path.join(self.data_dir, base_dir, entry['dataset'], entry['filename'])\n        img = Image.open(image_path).convert(\"RGB\")\n        img = self.image_transform(img)\n\n        metadata = {\n            'image_path': image_path,\n            'dataset': entry['dataset'],\n            'scene': entry['scene'],\n            'filename': entry['filename'],\n        }\n\n        if not self.is_train:\n            metadata['image_id'] = entry['image_id']\n                \n        else:\n            metadata['rotation_matrix'] = entry['rotation']\n            metadata['translation_vector'] = entry['translation']\n        \n        return img, metadata        ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.386623Z","iopub.execute_input":"2025-06-23T07:19:30.386867Z","iopub.status.idle":"2025-06-23T07:19:30.400498Z","shell.execute_reply.started":"2025-06-23T07:19:30.386845Z","shell.execute_reply":"2025-06-23T07:19:30.399859Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"@dataclasses.dataclass\nclass Predictions:\n    image_id: str\n    camera_id:int\n    datapath: str \n    filename: str\n    cluster_index: Optional[int] = None\n    colmap_model_id: Optional[int] = None\n    rotation: Optional[np.ndarray] = None\n    translation: Optional[np.ndarray] = None\n    intrinsics: Optional[np.ndarray] = None  \n    image_shape: Optional[Tuple[int, int]] = None  # (H, W)\n    depth_map: Optional[np.ndarray] = None  # sparse depth map\n    depth_points: Optional[np.ndarray] = None  # (N, 3) [u, v, z]\n\n\nimage_dir = \"/kaggle/input/image-matching-challenge-2025/\"\nres_dir = \"/kaggle/working/result\"\nto_train = True\ndevice = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\ntrain_csv = os.path.join(image_dir,'train_labels.csv')\nfeature_dir = os.path.join(res_dir, 'featureout')\ncolmap_output_dir = os.path.join(res_dir, \"output/colmap\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.401114Z","iopub.execute_input":"2025-06-23T07:19:30.401305Z","iopub.status.idle":"2025-06-23T07:19:30.420470Z","shell.execute_reply.started":"2025-06-23T07:19:30.401291Z","shell.execute_reply":"2025-06-23T07:19:30.419737Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Flow","metadata":{}},{"cell_type":"code","source":"def generate_image_pairs(clusters):\n    \"\"\"\n    Returns a set of unique image pairs using canonical (sorted) \n    tuple to prevent (i,j) vs (j,i) duplication.\n    \"\"\"\n    unique_pairs = set()\n    for cluster_id, img_list in clusters.items():\n        print(f\"Processing cluster {cluster_id} with {len(img_list)} images\")\n        if len(img_list) < 2:\n            continue\n        for img1, img2 in combinations(img_list, 2):\n            pair = tuple(sorted((img1, img2)))  \n            unique_pairs.add(pair)\n    return list(unique_pairs)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.439910Z","iopub.execute_input":"2025-06-23T07:19:30.440117Z","iopub.status.idle":"2025-06-23T07:19:30.453548Z","shell.execute_reply.started":"2025-06-23T07:19:30.440104Z","shell.execute_reply":"2025-06-23T07:19:30.452879Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"dataset = ImageMatching(\n            train_csv,\n            image_dir,\n            to_train,\n            datasets_filter=['imc2023_haiper','pt_brandenburg_british_buckingham','amy_gardens'], #'imc2023_theather_imc2024_church','ETs','amy_gardens', 'pt_brandenburg_british_buckingham\n            image_transform= image_transform)\n\ntrain_Dataloader = DataLoader(dataset, batch_size=16, shuffle=False)\n\n# STEP 1\nfeature, img_paths = global_feature_extractor(device, train_Dataloader)\nedges = build_knn_graph(feature, image_dir, img_paths, k=15)\nGraph = build_graph_from_edges(edges)\ncluster_map = louvian_cluster(Graph)\n\nclusters = defaultdict(list)\nfor idx, cluster_id in cluster_map.items():\n    clusters[cluster_id].append(img_paths[idx])\n\n# STEP 2\n\nimage_pairs = generate_image_pairs(clusters)\n\nmatch_pair_lightglue_superpoint(\n    image_paths,\n    image_pairs,\n    device,\n    min_matches=20\n)\n\n# STEP 3\n\n#intrinsics = [2048, 2048, 1024, 1024] # [fx,fy, cx, cy] approximates \n\nfull_img_dir = {}\nmodel_maps = {}\ncolmap_results = {}\n\nfor cluster_idx, cluster_imgpath in tqdm(clusters.items(), desc=\"Running COLMAP per cluster\"):\n    \n    pose_dict = defaultdict(dict)\n    colmap_preds = []\n    output_dir = os.path.join(colmap_output_dir, f\"cluster_{cluster_idx}\")\n    \n    maps = run_colmap_reconstruction(\n    cluster_imgpath,\n    feature_dir,\n    cluster_idx,\n    output_dir\n    )    \n\n    for i in cluster_imgpath:\n        scene_img = \"/\".join(i.split(\"/\")[-2:])\n        full_img_dir[scene_img] = i\n        \n    \n    for model_id, model in tqdm(maps.items(), desc = \"Pose Extraction\"):\n        count = 0\n        \n        for image_id, image in model.images.items():\n            pose_dict[image.name] = {\n                \"rotation\": image.cam_from_world.rotation.matrix(),\n                \"translation\": image.cam_from_world.translation\n            }\n\n            cam = model.cameras[image.camera_id]\n            img_shape = (cam.height, cam.width)\n            depth_map, projected_pts = sparse_depth_map(\n                    image, cam, model.points3D, img_shape)\n            colmap_preds.append(Predictions(\n                image_id=image.image_id,\n                camera_id= image.camera_id,\n                datapath=full_img_dir[image.name],\n                filename=image.name,\n                cluster_index=cluster_idx,\n                colmap_model_id = model_id,\n                rotation=image.cam_from_world.rotation.matrix(),\n                translation=image.cam_from_world.translation,\n                intrinsics=cam.calibration_matrix(),\n                image_shape=img_shape,\n                depth_map=depth_map,\n                depth_points=projected_pts                \n                ))\n                \n    colmap_results[cluster_idx] = {\n        \"models\" : maps,\n        \"colmap_res\" : colmap_preds,\n        \"pose_rt\" : pose_dict \n    }\n\n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T07:19:30.454290Z","iopub.execute_input":"2025-06-23T07:19:30.454564Z","iopub.status.idle":"2025-06-23T08:29:56.594720Z","shell.execute_reply.started":"2025-06-23T07:19:30.454546Z","shell.execute_reply":"2025-06-23T08:29:56.594076Z"},"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\nfig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\n\nfor idx,_ in clusters.items():\n    poses = colmap_results[idx][\"pose_rt\"]\n    for img_name, pose in poses.items():\n        c = -pose[\"rotation\"].T @ pose[\"translation\"]\n        ax.scatter(c[0], c[1], c[2], c='r')\n\nax.set_title(\"Camera Centers from COLMAP\")\nax.set_xlabel(\"X\")\nax.set_ylabel(\"Y\")\nax.set_zlabel(\"Z\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:29:56.595724Z","iopub.execute_input":"2025-06-23T08:29:56.595964Z","iopub.status.idle":"2025-06-23T08:30:00.225576Z","shell.execute_reply.started":"2025-06-23T08:29:56.595946Z","shell.execute_reply":"2025-06-23T08:30:00.224814Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# from cluster 0\npose_img1 = 'imc2023_haiper/fountain_image_012.png' \npose_img2 = 'imc2023_haiper/fountain_image_025.png'\n\nkeypoint_match_visualiser(pose_img1, pose_img2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:48:04.303957Z","iopub.execute_input":"2025-06-23T08:48:04.304489Z","iopub.status.idle":"2025-06-23T08:48:08.085469Z","shell.execute_reply.started":"2025-06-23T08:48:04.304464Z","shell.execute_reply":"2025-06-23T08:48:08.084761Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# cluster 1\npose1 = 'imc2023_haiper/bike_image_115.png'\npose2 = 'imc2023_haiper/bike_image_119.png'\nkeypoint_match_visualiser(pose1, pose2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:30:05.916380Z","iopub.execute_input":"2025-06-23T08:30:05.916619Z","iopub.status.idle":"2025-06-23T08:30:09.978075Z","shell.execute_reply.started":"2025-06-23T08:30:05.916603Z","shell.execute_reply":"2025-06-23T08:30:09.977302Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# cluster 5\npose_img1 = 'imc2023_theather_imc2024_church/kyiv_puppet_theater_img_20220127_165703.png'\npose_img2 = 'imc2023_theather_imc2024_church/kyiv_puppet_theater_img_20220127_165709.png'\nkeypoint_match_visualiser(pose_img1, pose_img2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:30:09.978886Z","iopub.execute_input":"2025-06-23T08:30:09.979092Z","iopub.status.idle":"2025-06-23T08:30:13.845433Z","shell.execute_reply.started":"2025-06-23T08:30:09.979075Z","shell.execute_reply":"2025-06-23T08:30:13.844725Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#cluster 5\npose_img1 = 'imc2023_theather_imc2024_church/church_00004.png'\npose_img2 = 'imc2023_theather_imc2024_church/church_00005.png'\nkeypoint_match_visualiser(pose_img1, pose_img2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:30:13.846363Z","iopub.execute_input":"2025-06-23T08:30:13.846629Z","iopub.status.idle":"2025-06-23T08:30:16.454525Z","shell.execute_reply.started":"2025-06-23T08:30:13.846610Z","shell.execute_reply":"2025-06-23T08:30:16.453799Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pose_img1 = 'pt_brandenburg_british_buckingham/british_museum_04790098_4924719784.png'\npose_img2 = 'pt_brandenburg_british_buckingham/british_museum_09218102_3818070255.png'\nkeypoint_match_visualiser(pose_img1, pose_img2)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:30:16.455422Z","iopub.execute_input":"2025-06-23T08:30:16.455690Z","iopub.status.idle":"2025-06-23T08:30:19.644276Z","shell.execute_reply.started":"2025-06-23T08:30:16.455669Z","shell.execute_reply":"2025-06-23T08:30:19.643171Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Depth Points Visualization\n> * ### The depth point overlays visualize sparse 3D points, triangulated by COLMAP using geometrically verified 2D matches obtained from LightGlue + SuperPoint. These points are projected back onto the original image planes and color-coded by depth.\n> * ### Closer points → cool colors (purple/blue)\n> * ### Farther points → warm colors (red/yellow)","metadata":{}},{"cell_type":"code","source":"# STEP 4\n# Depth Points visualisation\n\nimg_count = 0\n\nfor idx, _ in clusters.items():\n    colmap_preds = colmap_results[idx][\"colmap_res\"]\n    \n    for i in range(0, len(colmap_preds) - 1, 2): \n        preds_pair = colmap_preds[i:i+2]\n        fig, axs = plt.subplots(1, 2, figsize=(16, 8))\n\n        for j, pred in enumerate(preds_pair):\n            rgb = np.array(Image.open(pred.datapath).convert(\"RGB\")).astype(np.float32) / 255.0\n            depth_overlay = np.zeros_like(rgb)\n\n            if pred.depth_points.size == 0:\n                continue\n\n            depths = pred.depth_points[:, 2]\n            z_min, z_max = depths.min(), depths.max()\n\n            for x, y, z in pred.depth_points:\n                x, y = int(round(x)), int(round(y))\n                if 0 <= y < rgb.shape[0] and 0 <= x < rgb.shape[1]:\n                    norm_z = (z - z_min) / (z_max - z_min + 1e-8)\n                    color = plt.cm.plasma(norm_z)[:3]\n                    cv2.circle(depth_overlay, (x, y), radius=3, color=color[::-1], thickness=12)\n\n            blended = 0.6 * rgb + 0.4 * depth_overlay\n            axs[j].imshow(blended)\n            axs[j].set_title(f\"Depth Points: {pred.filename}\")\n            axs[j].axis(\"off\")\n\n        plt.tight_layout()\n        plt.show()\n        \n        img_count += 2\n        if img_count >= 4:\n            break\n\n            \n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:53:10.993256Z","iopub.execute_input":"2025-06-23T08:53:10.994144Z","iopub.status.idle":"2025-06-23T08:53:22.667010Z","shell.execute_reply.started":"2025-06-23T08:53:10.994119Z","shell.execute_reply":"2025-06-23T08:53:22.666250Z"},"collapsed":true,"jupyter":{"outputs_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Epipolar Geometry Insight\n> ### In epipolar geometry, each matched point between two images defines a line of sight extending in 3D space. The epipolar line in one image indicates where the corresponding point should lie in the other image, assuming a consistent 3D structure and camera motion. When multiple views are available (i.e., different camera orientations), the intersection of these lines in 3D space reveals the true depth of the scene point—this is the principle behind triangulation in Structure-from-Motion (SfM). Thus, the colored points lying along corresponding epipolar lines across image pairs indicate geometrically verified matches, which COLMAP uses to build both the 3D point cloud and sparse depth maps.","metadata":{}},{"cell_type":"code","source":"\n# epipolar on depth plot\nimg_nameA = 'imc2023_haiper/fountain_image_012.png' \nimg_nameB = 'imc2023_haiper/fountain_image_025.png'\nfor idx,_ in clusters.items():\n    model = colmap_results[idx][\"models\"]\n    colmap_preds = colmap_results[idx][\"colmap_res\"]\n    predA = predB = None\n    for pred in colmap_preds:\n        if pred.filename == img_nameA:\n            predA = pred\n            img_model = model[predA.colmap_model_id]\n        if pred.filename == img_nameB:\n            predB = pred\n            \n        if predA and predB and img_model:            \n            print(img_model)\n            epipolar(img_model, predA, predB, max_points=80)\n            break\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-23T08:55:45.604436Z","iopub.execute_input":"2025-06-23T08:55:45.605201Z","iopub.status.idle":"2025-06-23T08:55:46.945633Z","shell.execute_reply.started":"2025-06-23T08:55:45.605164Z","shell.execute_reply":"2025-06-23T08:55:46.944623Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"* ### Visualized epipolar lines across the image pair using only shared 3D points.\n* ### Each color corresponds to a point match (derived from a 3D point seen in both images).\n* ### The epipolar line drawn on the second image shows where that 3D point should appear, given the camera motion.\n* ### The fact that the point lies on the line confirms geometric correctness.\n* ### This is also a consistency check for triangulated 3D structure and camera motion.","metadata":{}}]}