{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":71885,"databundleVersionId":8143495,"sourceType":"competition"},{"sourceId":1650695,"sourceType":"datasetVersion","datasetId":976194},{"sourceId":3149286,"sourceType":"datasetVersion","datasetId":1916691},{"sourceId":6204030,"sourceType":"datasetVersion","datasetId":3277821},{"sourceId":7884485,"sourceType":"datasetVersion","datasetId":4628051},{"sourceId":7884725,"sourceType":"datasetVersion","datasetId":4628331},{"sourceId":4534,"sourceType":"modelInstanceVersion","modelInstanceId":3326,"modelId":986},{"sourceId":17191,"sourceType":"modelInstanceVersion","modelInstanceId":14317,"modelId":21716},{"sourceId":17555,"sourceType":"modelInstanceVersion","modelInstanceId":14611,"modelId":22086},{"sourceId":41690,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":35049,"modelId":49313}],"dockerImageVersionId":30665,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Constructing precise 3D maps from sets of images in diverse scenarios and environments including Structure from Motion (SfM), Multi-View Stereo (MVS)","metadata":{}},{"cell_type":"code","source":"# import cv2\n# import numpy as np\n# import os\n\n# def load_images_from_folder(folder):\n# images = []\n# for filename in os.listdir(folder):\n# img = cv2.imread(os.path.join(folder, filename))\n# if img is not None:\n# images.append(img)\n# return images\n\n# folder_path = '/kaggle/input/image-matching-challenge-2024/train/multi-temporal-temple-baalshamin/images/'\n# images = load_images_from_folder(folder_path)\n\n\n# def detect_and_match_features(images):\n# sift = cv2.SIFT_create()\n# keypoints_list = []\n# descriptors_list = []\n\n# for img in images:\n# keypoints, descriptors = sift.detectAndCompute(img, None)\n# keypoints_list.append(keypoints)\n# descriptors_list.append(descriptors)\n\n# return keypoints_list, descriptors_list\n\n# def match_features(descriptors_list):\n# bf = cv2.BFMatcher()\n# matches = []\n# for i in range(len(descriptors_list) - 1):\n# matches.append(bf.knnMatch(descriptors_list[i], descriptors_list[i + 1], k=2))\n# return matches\n\n# keypoints_list, descriptors_list = detect_and_match_features(images)\n# matches = match_features(descriptors_list)\n\n\n\n# def estimate_camera_pose(matches, keypoints_list, K):\n# poses = [np.eye(4)]\n# for i, match in enumerate(matches):\n# good_matches = [m for m, n in match if m.distance < 0.75 * n.distance]\n# if len(good_matches) < 5:\n# print(f\"Not enough good matches between image {i} and image {i+1}\")\n# continue\n# src_pts = np.float32([keypoints_list[i][m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n# dst_pts = np.float32([keypoints_list[i + 1][m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n# print(f\"src_pts shape: {src_pts.shape}, dst_pts shape: {dst_pts.shape}\")\n# E, mask = cv2.findEssentialMat(src_pts, dst_pts, K)\n# if E is None or E.shape[0] != 3 or E.shape[1] != 3:\n# print(f\"Invalid Essential Matrix between image {i} and image {i+1}\")\n# continue\n# _, R, t, _ = cv2.recoverPose(E, src_pts, dst_pts, K)\n# pose = np.eye(4)\n# pose[:3, :3] = R\n# pose[:3, 3] = t.squeeze()\n# poses.append(pose @ poses[-1])\n# return poses\n\n# fx, fy, cx, cy = 1000, 1000, 320, 240 # Example values\n# K = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1]])\n\n# poses = estimate_camera_pose(matches, keypoints_list, K)\n\n\n# import g2o\n\n# class BundleAdjustment:\n# def __init__(self):\n# self.optimizer = g2o.SparseOptimizer()\n# solver = g2o.BlockSolverSE3(g2o.LinearSolverCSparseSE3())\n# self.optimizer.set_algorithm(g2o.OptimizationAlgorithmLevenberg(solver))\n# self.vertices = []\n# self.edges = []\n\n# def add_pose(self, pose):\n# vertex = g2o.VertexSE3()\n# vertex.set_id(len(self.vertices))\n# vertex.set_estimate(g2o.SE3Quat(pose[:3, :3], pose[:3, 3]))\n# self.optimizer.add_vertex(vertex)\n# self.vertices.append(vertex)\n\n# def add_edge(self, v0_id, v1_id, measurement, information=np.identity(6)):\n# edge = g2o.EdgeSE3()\n# edge.set_vertex(0, self.optimizer.vertex(v0_id))\n# edge.set_vertex(1, self.optimizer.vertex(v1_id))\n# edge.set_measurement(g2o.SE3Quat(measurement[:3, :3], measurement[:3, 3]))\n# edge.set_information(information)\n# self.optimizer.add_edge(edge)\n# self.edges.append(edge)\n\n# def optimize(self, max_iterations=10):\n# self.optimizer.initialize_optimization()\n# self.optimizer.optimize(max_iterations)\n\n# bundle_adjustment = BundleAdjustment()\n\n# # Add poses\n# for pose in poses:\n# bundle_adjustment.add_pose(pose)\n\n# # Add edges\n# for i in range(len(poses) - 1):\n# relative_pose = np.linalg.inv(poses[i]) @ poses[i + 1]\n# bundle_adjustment.add_edge(i, i + 1, relative_pose)\n\n# bundle_adjustment.optimize()\n\n\n# def plot_poses(poses):\n# fig = plt.figure()\n# ax = fig.add_subplot(111, projection='3d')\n# for pose in poses:\n# ax.scatter(pose[0, 3], pose[1, 3], pose[2, 3])\n# ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[0, 0], pose[0, 1], pose[0, 2], length=0.1, color='r')\n# ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[1, 0], pose[1, 1], pose[1, 2], length=0.1, color='g')\n# ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[2, 0], pose[2, 1], pose[2, 2], length=0.1, color='b')\n# ax.set_xlabel('X')\n# ax.set_ylabel('Y')\n# ax.set_zlabel('Z')\n# plt.show()\n\n# optimized_poses = [vertex.estimate().matrix() for vertex in bundle_adjustment.vertices]\n# plot_poses(optimized_poses)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-07-05T01:48:16.583967Z","iopub.execute_input":"2024-07-05T01:48:16.584781Z","iopub.status.idle":"2024-07-05T01:48:16.594957Z","shell.execute_reply.started":"2024-07-05T01:48:16.584748Z","shell.execute_reply":"2024-07-05T01:48:16.594027Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# OpenCV's built-in functions\n- Constructing precise 3D maps using sets of images with OpenCV for feature detection, matching, and camera pose estimation:\n - Detect and Match Features: Features are detected and matched using SIFT.\n - Estimate Camera Poses: Estimate the camera poses using essential matrices and pose recovery.\n - Visualize the Results:visualize camera poses in a 3D plot.","metadata":{}},{"cell_type":"code","source":"import cv2\nimport numpy as np\nimport os\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\n\ndef load_images_from_folder(folder):\n images = []\n for filename in os.listdir(folder):\n img = cv2.imread(os.path.join(folder, filename))\n if img is not None:\n images.append(img)\n return images\n\nfolder_path = '/kaggle/input/image-matching-challenge-2024/train/dioscuri/images/'\nimages = load_images_from_folder(folder_path)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T03:44:42.051635Z","iopub.execute_input":"2024-06-02T03:44:42.052539Z","iopub.status.idle":"2024-06-02T03:44:44.760157Z","shell.execute_reply.started":"2024-06-02T03:44:42.052501Z","shell.execute_reply":"2024-06-02T03:44:44.758637Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def detect_and_match_features(images):\n sift = cv2.SIFT_create()\n keypoints_list = []\n descriptors_list = []\n\n for img in images:\n keypoints, descriptors = sift.detectAndCompute(img, None)\n keypoints_list.append(keypoints)\n descriptors_list.append(descriptors)\n\n return keypoints_list, descriptors_list\n\ndef match_features(descriptors_list):\n bf = cv2.BFMatcher()\n matches = []\n for i in range(len(descriptors_list) - 1):\n matches.append(bf.knnMatch(descriptors_list[i], descriptors_list[i + 1], k=2))\n return matches\n\nkeypoints_list, descriptors_list = detect_and_match_features(images)\nmatches = match_features(descriptors_list)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T03:44:53.998300Z","iopub.execute_input":"2024-06-02T03:44:53.999203Z","iopub.status.idle":"2024-06-02T03:45:32.652068Z","shell.execute_reply.started":"2024-06-02T03:44:53.999154Z","shell.execute_reply":"2024-06-02T03:45:32.651167Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def estimate_camera_pose(matches, keypoints_list, K):\n poses = [np.eye(4)]\n for i, match in enumerate(matches):\n good_matches = [m for m, n in match if m.distance < 0.75 * n.distance]\n if len(good_matches) < 5:\n print(f\"Not enough good matches between image {i} and image {i+1}\")\n continue\n src_pts = np.float32([keypoints_list[i][m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n dst_pts = np.float32([keypoints_list[i + 1][m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n print(f\"src_pts shape: {src_pts.shape}, dst_pts shape: {dst_pts.shape}\")\n E, mask = cv2.findEssentialMat(src_pts, dst_pts, K)\n if E is None or E.shape[0] != 3 or E.shape[1] != 3:\n print(f\"Invalid Essential Matrix between image {i} and image {i+1}\")\n continue\n _, R, t, _ = cv2.recoverPose(E, src_pts, dst_pts, K)\n pose = np.eye(4)\n pose[:3, :3] = R\n pose[:3, 3] = t.squeeze()\n poses.append(pose @ poses[-1])\n return poses\n\nfx, fy, cx, cy = 1000, 1000, 320, 240 # Example values\nK = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1]])\n\nposes = estimate_camera_pose(matches, keypoints_list, K)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T03:46:19.527539Z","iopub.execute_input":"2024-06-02T03:46:19.528447Z","iopub.status.idle":"2024-06-02T03:46:33.652294Z","shell.execute_reply.started":"2024-06-02T03:46:19.528392Z","shell.execute_reply":"2024-06-02T03:46:33.651255Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_poses(poses):\n fig = plt.figure()\n ax = fig.add_subplot(111, projection='3d')\n for pose in poses:\n ax.scatter(pose[0, 3], pose[1, 3], pose[2, 3])\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[0, 0], pose[0, 1], pose[0, 2], length=0.1, color='r')\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[1, 0], pose[1, 1], pose[1, 2], length=0.1, color='g')\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[2, 0], pose[2, 1], pose[2, 2], length=0.1, color='b')\n ax.set_xlabel('X')\n ax.set_ylabel('Y')\n ax.set_zlabel('Z')\n plt.show()\n\nplot_poses(poses)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T03:46:42.627731Z","iopub.execute_input":"2024-06-02T03:46:42.628271Z","iopub.status.idle":"2024-06-02T03:46:44.449134Z","shell.execute_reply.started":"2024-06-02T03:46:42.628230Z","shell.execute_reply":"2024-06-02T03:46:44.447711Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# import tensorflow as tf\n# import tensorflow_datasets as tfds\n# import cv2\n# import numpy as np\n# import matplotlib.pyplot as plt\n# from mpl_toolkits.mplot3d import Axes3D\n\n# # Oxford IIIT Pet dataset\n# dataset, info = tfds.load('oxford_iiit_pet', with_info=True)\n# train_dataset = dataset['train']\n\n# def extract_images(dataset, num_images=10):\n# images = []\n# for example in tfds.as_numpy(dataset.take(num_images)):\n# img = example['image']\n# img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)\n# images.append(img)\n# return images\n\n# num_images = 10 \n# images = extract_images(train_dataset, num_images)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T04:18:22.952527Z","iopub.execute_input":"2024-06-02T04:18:22.953260Z","iopub.status.idle":"2024-06-02T04:18:22.958568Z","shell.execute_reply.started":"2024-06-02T04:18:22.953229Z","shell.execute_reply":"2024-06-02T04:18:22.957196Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# In-Built Dataset \n - ImageCollection feature from skimage to load a sequence of images. Images from the data module in skimage or cv2's sample data. ","metadata":{}},{"cell_type":"code","source":"import cv2\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage import data\nfrom skimage.color import rgb2gray\nfrom skimage.transform import rotate, AffineTransform, warp\n\nimage = rgb2gray(data.astronaut())\n\n# sequence of images by applying transformations\nimages = [image]\ntransformations = [\n AffineTransform(translation=(10, 0)),\n AffineTransform(translation=(-10, 0)),\n AffineTransform(translation=(0, 10)),\n AffineTransform(translation=(0, -10)),\n rotate(image, angle=15, resize=False),\n rotate(image, angle=-15, resize=False)\n]\n\nfor tform in transformations:\n if isinstance(tform, AffineTransform):\n images.append(warp(image, tform.inverse, output_shape=image.shape))\n else:\n images.append(tform)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T04:18:28.109178Z","iopub.execute_input":"2024-06-02T04:18:28.109578Z","iopub.status.idle":"2024-06-02T04:18:28.677812Z","shell.execute_reply.started":"2024-06-02T04:18:28.109549Z","shell.execute_reply":"2024-06-02T04:18:28.676131Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def detect_and_match_features(images):\n sift = cv2.SIFT_create()\n keypoints_list = []\n descriptors_list = []\n\n for img in images:\n img = (img * 255).astype(np.uint8) # Convert to uint8 for SIFT\n keypoints, descriptors = sift.detectAndCompute(img, None)\n keypoints_list.append(keypoints)\n descriptors_list.append(descriptors)\n\n return keypoints_list, descriptors_list\n\ndef match_features(descriptors_list):\n bf = cv2.BFMatcher()\n matches = []\n for i in range(len(descriptors_list) - 1):\n matches.append(bf.knnMatch(descriptors_list[i], descriptors_list[i + 1], k=2))\n return matches\n\nkeypoints_list, descriptors_list = detect_and_match_features(images)\nmatches = match_features(descriptors_list)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T04:18:47.227303Z","iopub.execute_input":"2024-06-02T04:18:47.228160Z","iopub.status.idle":"2024-06-02T04:18:47.927343Z","shell.execute_reply.started":"2024-06-02T04:18:47.228128Z","shell.execute_reply":"2024-06-02T04:18:47.926461Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def estimate_camera_pose(matches, keypoints_list, K):\n poses = [np.eye(4)]\n for i, match in enumerate(matches):\n good_matches = [m for m, n in match if m.distance < 0.75 * n.distance]\n if len(good_matches) < 5:\n print(f\"Not enough good matches between image {i} and image {i+1}\")\n continue\n src_pts = np.float32([keypoints_list[i][m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n dst_pts = np.float32([keypoints_list[i + 1][m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)\n print(f\"src_pts shape: {src_pts.shape}, dst_pts shape: {dst_pts.shape}\")\n E, mask = cv2.findEssentialMat(src_pts, dst_pts, K)\n if E is None or E.shape[0] != 3 or E.shape[1] != 3:\n print(f\"Invalid Essential Matrix between image {i} and image {i+1}\")\n continue\n _, R, t, _ = cv2.recoverPose(E, src_pts, dst_pts, K)\n pose = np.eye(4)\n pose[:3, :3] = R\n pose[:3, 3] = t.squeeze()\n poses.append(pose @ poses[-1])\n return poses\n\n# Replace these with example intrinsic parameters\nfx, fy, cx, cy = 1000, 1000, image.shape[1] // 2, image.shape[0] // 2\nK = np.array([[fx, 0, cx], [0, fy, cy], [0, 0, 1]])\n\nposes = estimate_camera_pose(matches, keypoints_list, K)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T04:18:58.087508Z","iopub.execute_input":"2024-06-02T04:18:58.087911Z","iopub.status.idle":"2024-06-02T04:18:58.166039Z","shell.execute_reply.started":"2024-06-02T04:18:58.087879Z","shell.execute_reply":"2024-06-02T04:18:58.164821Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_poses(poses):\n fig = plt.figure()\n ax = fig.add_subplot(111, projection='3d')\n for pose in poses:\n ax.scatter(pose[0, 3], pose[1, 3], pose[2, 3])\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[0, 0], pose[0, 1], pose[0, 2], length=0.1, color='r')\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[1, 0], pose[1, 1], pose[1, 2], length=0.1, color='g')\n ax.quiver(pose[0, 3], pose[1, 3], pose[2, 3], pose[2, 0], pose[2, 1], pose[2, 2], length=0.1, color='b')\n ax.set_xlabel('X')\n ax.set_ylabel('Y')\n ax.set_zlabel('Z')\n plt.show()\n\nplot_poses(poses)\n","metadata":{"execution":{"iopub.status.busy":"2024-06-02T04:19:09.277117Z","iopub.execute_input":"2024-06-02T04:19:09.277960Z","iopub.status.idle":"2024-06-02T04:19:09.640303Z","shell.execute_reply.started":"2024-06-02T04:19:09.277926Z","shell.execute_reply":"2024-06-02T04:19:09.638946Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # !pip install open3d\n\n# import open3d as o3d\n\n# def create_camera_actor(pose, scale=0.1):\n# camera = o3d.geometry.LineSet()\n \n# # camera frame in homogeneous coordinates\n# camera_frame = np.array([\n# [0, 0, 0, 1],\n# [1, 0, 0, 1],\n# [0, 1, 0, 1],\n# [0, 0, 1, 1],\n# [1, 1, 1, 1]\n# ]).T * scale\n \n# # Transform camera frame by the pose\n# transformed_frame = (pose @ camera_frame).T[:, :3]\n \n# points = [\n# transformed_frame[0],\n# transformed_frame[1],\n# transformed_frame[2],\n# transformed_frame[3]\n# ]\n \n# lines = [\n# [0, 1],\n# [0, 2],\n# [0, 3],\n# [1, 2],\n# [2, 3],\n# [3, 1]\n# ]\n \n# colors = [[1, 0, 0] for _ in range(len(lines))]\n \n# camera.points = o3d.utility.Vector3dVector(points)\n# camera.lines = o3d.utility.Vector2iVector(lines)\n# camera.colors = o3d.utility.Vector3dVector(colors)\n \n# return camera\n\n# def plot_poses_with_open3d(poses):\n# vis = o3d.visualization.Visualizer()\n# vis.create_window()\n \n# for pose in poses:\n# camera_actor = create_camera_actor(pose)\n# vis.add_geometry(camera_actor)\n \n# vis.run()\n# vis.destroy_window()\n\n# plot_poses_with_open3d(poses)\n","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-06-02T04:27:10.018717Z","iopub.execute_input":"2024-06-02T04:27:10.019922Z","iopub.status.idle":"2024-06-02T04:27:10.026579Z","shell.execute_reply.started":"2024-06-02T04:27:10.019875Z","shell.execute_reply":"2024-06-02T04:27:10.025094Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # url check\n# import urllib.request\n# import time\n\n# def fetch_url(url, retries=5, delay=2):\n# for attempt in range(retries):\n# try:\n# response = urllib.request.urlopen(url)\n# return response.read()\n# except urllib.error.URLError as e:\n# print(f\"Attempt {attempt + 1} failed: {e.reason}\")\n# time.sleep(delay)\n# raise Exception(f\"Failed to fetch URL after {retries} attempts\")\n\n# url = \"http://example.com\"\n# try:\n# content = fetch_url(url)\n# print(\"Successfully fetched URL\")\n# except Exception as e:\n# print(e)\n","metadata":{"execution":{"iopub.status.busy":"2024-07-05T01:54:19.194489Z","iopub.execute_input":"2024-07-05T01:54:19.195301Z","iopub.status.idle":"2024-07-05T01:54:19.199618Z","shell.execute_reply.started":"2024-07-05T01:54:19.195271Z","shell.execute_reply":"2024-07-05T01:54:19.198530Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import torch\nimport torchvision.models as models\n\n# full pre-trained VGG16 model\nmodel = models.vgg16(pretrained=False)\n\n# pre-trained weights\nmodel_weights_path = \"/kaggle/input/vgg16-397923af/pytorch/vgg16-397923af/1/vgg16-397923af.pth\"\nmodel.load_state_dict(torch.load(model_weights_path))\n\n# Extract features part of the VGG16 model\nmodel = model.features\nmodel.eval()\n","metadata":{"execution":{"iopub.status.busy":"2024-07-05T01:56:00.235013Z","iopub.execute_input":"2024-07-05T01:56:00.235396Z","iopub.status.idle":"2024-07-05T01:56:03.088036Z","shell.execute_reply.started":"2024-07-05T01:56:00.235364Z","shell.execute_reply":"2024-07-05T01:56:03.087070Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Match sections of one image to another by creating patches, extracting features using a pre-trained convolutional neural network (VGG16), and comparing these features to find the best match.","metadata":{}},{"cell_type":"code","source":"import torch\nimport torchvision.transforms as transforms\nfrom PIL import Image\nimport torchvision.models as models\nimport torch.nn.functional as F\n\ndef load_image(image_path):\n image = Image.open(image_path)\n return image\n\n\ndef preprocess_image(image):\n preprocess = transforms.Compose([\n transforms.Resize((224, 224)),\n transforms.ToTensor(),\n transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),\n ])\n image = preprocess(image)\n return image\n\n# patches from images\ndef create_patches(image, patch_size):\n patches = []\n _, h, w = image.size()\n for i in range(0, h, patch_size):\n for j in range(0, w, patch_size):\n patch = image[:, i:i+patch_size, j:j+patch_size]\n if patch.size(1) == patch_size and patch.size(2) == patch_size:\n patches.append(patch)\n return patches\n\n# Extract features using pre-trained CNN\ndef extract_features(patch, model):\n with torch.no_grad():\n patch = patch.unsqueeze(0) # batch dimension\n features = model(patch)\n return features\n\n# Comparing features to find best match\ndef compare_features(features1, features2):\n similarity = []\n for f1 in features1:\n for f2 in features2:\n f1 = f1.view(f1.size(0), -1) # Flatten the features\n f2 = f2.view(f2.size(0), -1)\n sim = F.cosine_similarity(f1, f2)\n similarity.append(sim.item())\n return similarity\n\n# images\noriginal_image_path = \"/kaggle/input/image-matching-challenge-2024/train/church/images/00003.png\"\ntarget_image_path = \"/kaggle/input/image-matching-challenge-2024/train/multi-temporal-temple-baalshamin/images/070313_img_6187_syria_palmyra_roman_ruins_999x.png\"\n\n\noriginal_image = load_image(original_image_path)\ntarget_image = load_image(target_image_path)\n\noriginal_image = preprocess_image(original_image)\ntarget_image = preprocess_image(target_image)\n\n# patches from images\npatch_size = 64 # Define the size of each patch\noriginal_patches = create_patches(original_image, patch_size)\ntarget_patches = create_patches(target_image, patch_size)\n\n# full pre-trained VGG16 model\nmodel = models.vgg16(pretrained=False)\n\n# pre-trained weights\nmodel_weights_path = \"/kaggle/input/vgg16-397923af/pytorch/vgg16-397923af/1/vgg16-397923af.pth\"\nmodel.load_state_dict(torch.load(model_weights_path))\n\n# Extract features part of the VGG16 model\nmodel = model.features\nmodel.eval()\n\n# Extract features from patches\noriginal_features = [extract_features(patch, model) for patch in original_patches]\ntarget_features = [extract_features(patch, model) for patch in target_patches]\n\n# best match\nsimilarities = compare_features(original_features, target_features)\nbest_match = max(similarities)\nprint(f\"Best Match Similarity: {best_match}\")\n","metadata":{"execution":{"iopub.status.busy":"2024-07-05T01:57:45.267266Z","iopub.execute_input":"2024-07-05T01:57:45.268140Z","iopub.status.idle":"2024-07-05T01:57:48.713192Z","shell.execute_reply.started":"2024-07-05T01:57:45.268104Z","shell.execute_reply":"2024-07-05T01:57:48.712219Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Patch visualization and aggregation of top-N similarities","metadata":{}},{"cell_type":"code","source":"import torch\nimport torchvision.transforms as transforms\nfrom PIL import Image\nimport torchvision.models as models\nimport torch.nn.functional as F\nimport matplotlib.pyplot as plt\n\n\ndef load_image(image_path):\n image = Image.open(image_path)\n return image\n\n\ndef preprocess_image(image):\n preprocess = transforms.Compose([\n transforms.Resize((224, 224)),\n transforms.ToTensor(),\n transforms.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),\n ])\n image = preprocess(image)\n return image\n\n# patches from images\ndef create_patches(image, patch_size):\n patches = []\n _, h, w = image.size()\n for i in range(0, h, patch_size):\n for j in range(0, w, patch_size):\n patch = image[:, i:i+patch_size, j:j+patch_size]\n if patch.size(1) == patch_size and patch.size(2) == patch_size:\n patches.append(patch)\n return patches\n\n\ndef extract_features(patch, model):\n with torch.no_grad():\n patch = patch.unsqueeze(0) # Batch dimension\n features = model(patch)\n return features\n\n# Comparing features to find best match\ndef compare_features(features1, features2):\n similarity = []\n for f1 in features1:\n for f2 in features2:\n f1 = f1.view(f1.size(0), -1) # Flatten features\n f2 = f2.view(f2.size(0), -1)\n sim = F.cosine_similarity(f1, f2)\n similarity.append(sim.item())\n return similarity\n\n\noriginal_image_path = \"/kaggle/input/image-matching-challenge-2024/train/church/images/00003.png\"\ntarget_image_path = \"/kaggle/input/image-matching-challenge-2024/train/multi-temporal-temple-baalshamin/images/070313_img_6187_syria_palmyra_roman_ruins_999x.png\"\n\n\noriginal_image = load_image(original_image_path)\ntarget_image = load_image(target_image_path)\n\noriginal_image = preprocess_image(original_image)\ntarget_image = preprocess_image(target_image)\n\n# Patches from images\npatch_size = 64 # Define the size of each patch\noriginal_patches = create_patches(original_image, patch_size)\ntarget_patches = create_patches(target_image, patch_size)\n\n# VGG16 model\nmodel = models.vgg16(pretrained=False)\n\n# pre-trained weights\nmodel_weights_path = \"/kaggle/input/vgg16-397923af/pytorch/vgg16-397923af/1/vgg16-397923af.pth\"\nmodel.load_state_dict(torch.load(model_weights_path))\n\n# features part of VGG16 model\nmodel = model.features\nmodel.eval()\n\n# features from patches\noriginal_features = [extract_features(patch, model) for patch in original_patches]\ntarget_features = [extract_features(patch, model) for patch in target_patches]\n\n# Compare features to find best match\nsimilarities = compare_features(original_features, target_features)\n\n# Top-N matches\nN = 5\ntop_n_similarities = sorted(similarities, reverse=True)[:N]\naverage_similarity = sum(top_n_similarities) / N\n\nprint(f\"Best Match Similarity: {max(similarities)}\")\nprint(f\"Average Top-{N} Match Similarity: {average_similarity}\")\n\n\ndef show_patches(patches, title, nrows=1):\n fig, axes = plt.subplots(nrows, len(patches)//nrows, figsize=(15, 3 * nrows))\n for ax, patch in zip(axes.flatten(), patches):\n patch = patch.permute(1, 2, 0) # Rearrange dimensions for plotting\n patch = patch * torch.tensor([0.229, 0.224, 0.225]) + torch.tensor([0.485, 0.456, 0.406]) # Denormalize\n patch = patch.clamp(0, 1) # Ensuring pixel values are between 0 and 1\n ax.imshow(patch)\n ax.axis('off')\n fig.suptitle(title)\n plt.tight_layout(rect=[0, 0.03, 1, 0.95])\n plt.show()\n\n# original patches\nshow_patches(original_patches[:10], \"Original Image Patches\", nrows=2)\n\n# target patches\nshow_patches(target_patches[:10], \"Target Image Patches\", nrows=2)\n","metadata":{"execution":{"iopub.status.busy":"2024-07-05T02:15:26.679939Z","iopub.execute_input":"2024-07-05T02:15:26.680310Z","iopub.status.idle":"2024-07-05T02:15:30.943155Z","shell.execute_reply.started":"2024-07-05T02:15:26.680282Z","shell.execute_reply":"2024-07-05T02:15:30.942179Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"#
📸 Image Matching Challenge - 📊 Understanding the baseline
\n

Reconstruct 3D scenes from 2D images over six different domains

\n\n***\n\nIn this notebook I explain the baseline solution provided by the organizers in [this notebook](https://www.kaggle.com/code/oldufo/imc-2024-submission-example).\n\nI have made the code a bit easier to read, adding comments and type annotations to make it easier for you to understand what is going on.\n\nHope you enjoy ❤️","metadata":{}},{"cell_type":"markdown","source":"### Matching many images across different viewpoints is a fundamental Computer Vision problem that is far from solved. ","metadata":{}},{"cell_type":"markdown","source":"# Structure from Motion","metadata":{}},{"cell_type":"markdown","source":"Structure from Motion (SfM) is the name given to the procedure of **reconstructing a 3D scene and simultaneously obtaining the camera poses of a camera w.r.t. the given scene**. This means that, as the name suggests, we are creating the entire rigid structure from a set of images with different view points (or equivalently a camera in motion).\n\nIn this competition, the important aspect of SfM we are interested in is *obtaining the camera poses* of where each image was taken, described by a rotation matrix and translation vector from the origin. These are the objects that will be scored in our submission!","metadata":{}},{"cell_type":"markdown","source":"
\"Structure
","metadata":{}},{"cell_type":"markdown","source":"# Baseline solution steps\nIn order to be able to estimate the camera poses, the solution provided by the organizers consists in the following steps:\n\n* [1. Find pairs of images that are similar](#1)\n* [2. Compute image keypoints](#2)\n* [3. Match keypoints between images](#3)\n* [4. Outlier detection with RANSAC](#4)\n* [5. Sparse reconstruction](#5)\n\nLet's understand how these steps are carried out","metadata":{}},{"cell_type":"markdown","source":"# Installing & importing relevant packages and models","metadata":{}},{"cell_type":"code","source":"!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/* /root/.cache/torch/hub/checkpoints/\n!cp /kaggle/input/lightglue/pytorch/aliked/1/* /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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-03-27T20:44:43.769283Z","iopub.execute_input":"2024-03-27T20:44:43.770095Z","iopub.status.idle":"2024-03-27T20:44:54.026101Z","shell.execute_reply.started":"2024-03-27T20:44:43.770064Z","shell.execute_reply":"2024-03-27T20:44:54.025Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# Finding image pairs\n\nTo find pairs of similar images, we use [DINOv2](https://arxiv.org/pdf/2304.07193.pdf) to obtain normalized image embeddings.\n\n
\"DINOv2
\nThen, we calculate the distances between all the embeddings, and only keep those below a given distance threshold. For images with less than a set minimum number of pairs, the closest ones are kept instead.","metadata":{}},{"cell_type":"code","source":"def embed_images(\n paths: list[Path],\n model_name: str,\n device: torch.device = torch.device(\"cpu\"),\n) -> T:\n \"\"\"Computes image embeddings.\n \n Returns a tensor of shape [len(filenames), output_dim]\n \"\"\"\n processor = AutoImageProcessor.from_pretrained(model_name)\n model = AutoModel.from_pretrained(model_name).eval().to(device)\n \n embeddings = []\n \n for i, path in tqdm(enumerate(paths), desc=\"Global descriptors\"):\n image = load_torch_image(path)\n \n with torch.inference_mode():\n inputs = processor(images=image, return_tensors=\"pt\", do_rescale=False).to(device)\n outputs = model(**inputs) # last_hidden_state and pooled\n \n # Max pooling over all the hidden states but the first (starting token)\n # To obtain a tensor of shape [1, output_dim]\n # We normalize so that distances are computed in a better fashion later\n embedding = F.normalize(outputs.last_hidden_state[:,1:].max(dim=1)[0], dim=-1, p=2)\n \n embeddings.append(embedding.detach().cpu())\n return torch.cat(embeddings, dim=0)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:19.965897Z","iopub.execute_input":"2024-03-27T19:26:19.966233Z","iopub.status.idle":"2024-03-27T19:26:20.278034Z","shell.execute_reply.started":"2024-03-27T19:26:19.9662Z","shell.execute_reply":"2024-03-27T19:26:20.276792Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_pairs_exhaustive(lst: list[Any]) -> list[tuple[int, int]]:\n \"\"\"Obtains all possible index pairs of a list\"\"\"\n return list(itertools.combinations(range(len(lst)), 2)) \n \ndef get_image_pairs(\n paths: list[Path],\n model_name: str,\n similarity_threshold: float = 0.6,\n tolerance: int = 1000,\n min_matches: int = 20,\n exhaustive_if_less: int = 20,\n p: float = 2.0,\n device: torch.device = torch.device(\"cpu\"),\n) -> list[tuple[int, int]]:\n \"\"\"Obtains pairs of similar images\"\"\"\n if len(paths) <= exhaustive_if_less:\n return get_pairs_exhaustive(paths)\n \n matches = []\n \n # Embed images and compute distances for filtering\n embeddings = embed_images(paths, model_name)\n distances = torch.cdist(embeddings, embeddings, p=p)\n \n # Remove pairs above similarity threshold (if enough)\n mask = distances <= similarity_threshold\n image_indices = np.arange(len(paths))\n \n for current_image_index in range(len(paths)):\n mask_row = mask[current_image_index]\n indices_to_match = image_indices[mask_row]\n \n # We don't have enough matches below the threshold, we pick most similar ones\n if len(indices_to_match) < min_matches:\n indices_to_match = np.argsort(distances[current_image_index])[:min_matches]\n \n for other_image_index in indices_to_match:\n # Skip an image matching itself\n if other_image_index == current_image_index:\n continue\n \n # We need to check if we are below a certain distance tolerance \n # since for images that don't have enough matches, we picked\n # the most similar ones (which could all still be very different \n # to the image we are analyzing)\n if distances[current_image_index, other_image_index] < tolerance:\n # Add the pair in a sorted manner to avoid redundancy\n matches.append(tuple(sorted((current_image_index, other_image_index.item()))))\n \n return sorted(list(set(matches)))","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.279657Z","iopub.execute_input":"2024-03-27T19:26:20.280099Z","iopub.status.idle":"2024-03-27T19:26:20.295316Z","shell.execute_reply.started":"2024-03-27T19:26:20.280059Z","shell.execute_reply":"2024-03-27T19:26:20.294217Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n images_list = list(Path(\"/kaggle/input/image-matching-challenge-2024/test/church/images/\").glob(\"*.png\"))[:10]\n index_pairs = get_image_pairs(images_list, \"/kaggle/input/dinov2/pytorch/base/1\")\n print(index_pairs)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.297963Z","iopub.execute_input":"2024-03-27T19:26:20.298297Z","iopub.status.idle":"2024-03-27T19:26:20.311085Z","shell.execute_reply.started":"2024-03-27T19:26:20.298262Z","shell.execute_reply":"2024-03-27T19:26:20.310156Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# Computing keypoints\n\nIn order to be able to know the position of each camera, we must be able to relate images to each other. For this, we extract relevant keypoints and compare pairs of image keypoints against each other. There are many ways to extract relevant keypoints, the most traditional one being [SIFT](https://en.wikipedia.org/wiki/Scale-invariant_feature_transform). However, newer and improved methods exist now, one of which is [ALIKED](https://arxiv.org/abs/2304.03608), the keypoint extraction method used in the solution.\n\n
\"ALIKED
\n\n\nLet's take a closer look at the keypoints that ALIKED extracts.","metadata":{}},{"cell_type":"code","source":"if DEBUG:\n dtype = torch.float32 # ALIKED has issues with float16\n\n extractor = ALIKED(\n max_num_keypoints=4096, \n detection_threshold=0.01, \n resize=1024\n ).eval().to(device, dtype)\n\n path = images_list[0]\n image = load_torch_image(path, device=device).to(dtype)\n features = extractor.extract(image)\n\n fig, ax = plt.subplots(1, 2, figsize=(10, 20))\n ax[0].imshow(image[0, ...].permute(1,2,0).cpu())\n ax[1].imshow(image[0, ...].permute(1,2,0).cpu())\n ax[1].scatter(features[\"keypoints\"][0, :, 0].cpu(), features[\"keypoints\"][0, :, 1].cpu(), s=0.5, c=\"red\")\n\n del extractor","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.312401Z","iopub.execute_input":"2024-03-27T19:26:20.312775Z","iopub.status.idle":"2024-03-27T19:26:20.323257Z","shell.execute_reply.started":"2024-03-27T19:26:20.312742Z","shell.execute_reply":"2024-03-27T19:26:20.322392Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def detect_keypoints(\n paths: list[Path],\n feature_dir: Path,\n num_features: int = 4096,\n resize_to: int = 1024,\n device: torch.device = torch.device(\"cpu\"),\n) -> None:\n \"\"\"Detects the keypoints in a list of images with ALIKED\n \n Stores them in feature_dir/keypoints.h5 and feature_dir/descriptors.h5\n to be used later with LightGlue\n \"\"\"\n dtype = torch.float32 # ALIKED has issues with float16\n \n extractor = ALIKED(\n max_num_keypoints=num_features, \n detection_threshold=0.01, \n resize=resize_to\n ).eval().to(device, dtype)\n \n feature_dir.mkdir(parents=True, exist_ok=True)\n \n with h5py.File(feature_dir / \"keypoints.h5\", mode=\"w\") as f_keypoints, \\\n h5py.File(feature_dir / \"descriptors.h5\", mode=\"w\") as f_descriptors:\n \n for path in tqdm(paths, desc=\"Computing keypoints\"):\n key = path.name\n \n with torch.inference_mode():\n image = load_torch_image(path, device=device).to(dtype)\n features = extractor.extract(image)\n \n f_keypoints[key] = features[\"keypoints\"].squeeze().detach().cpu().numpy()\n f_descriptors[key] = features[\"descriptors\"].squeeze().detach().cpu().numpy()","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.324258Z","iopub.execute_input":"2024-03-27T19:26:20.324524Z","iopub.status.idle":"2024-03-27T19:26:20.338576Z","shell.execute_reply.started":"2024-03-27T19:26:20.3245Z","shell.execute_reply":"2024-03-27T19:26:20.337641Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n feature_dir = Path(\"./sample_test_features\")\n detect_keypoints(images_list, feature_dir)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.339756Z","iopub.execute_input":"2024-03-27T19:26:20.340111Z","iopub.status.idle":"2024-03-27T19:26:20.34904Z","shell.execute_reply.started":"2024-03-27T19:26:20.340086Z","shell.execute_reply":"2024-03-27T19:26:20.3482Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# Match and compute keypoint distances\n\nNow that we have the relevant image pairs and keypoints, we can go ahead and compare the keypoints of the images in a pair to find a good relationship between them. This is done with [LightGlue](https://arxiv.org/abs/2306.13643), which matches the keypoints and their descriptors between two images.\n\n
\"LightGlue
","metadata":{}},{"cell_type":"code","source":"if DEBUG:\n matcher_params = {\n \"width_confidence\": -1,\n \"depth_confidence\": -1,\n \"mp\": True if 'cuda' in str(device) else False,\n }\n matcher = KF.LightGlueMatcher(\"aliked\", matcher_params).eval().to(device)\n\n with h5py.File(feature_dir / \"keypoints.h5\", mode=\"r\") as f_keypoints, \\\n h5py.File(feature_dir / \"descriptors.h5\", mode=\"r\") as f_descriptors:\n idx1, idx2 = index_pairs[0]\n key1, key2 = images_list[idx1].name, images_list[idx2].name\n\n keypoints1 = torch.from_numpy(f_keypoints[key1][...]).to(device)\n keypoints2 = torch.from_numpy(f_keypoints[key2][...]).to(device)\n print(\"Keypoints:\", keypoints1.shape, keypoints2.shape)\n descriptors1 = torch.from_numpy(f_descriptors[key1][...]).to(device)\n descriptors2 = torch.from_numpy(f_descriptors[key2][...]).to(device)\n print(\"Descriptors:\", descriptors1.shape, descriptors2.shape)\n\n with torch.inference_mode():\n distances, indices = matcher(\n descriptors1, \n descriptors2, \n KF.laf_from_center_scale_ori(keypoints1[None]),\n KF.laf_from_center_scale_ori(keypoints2[None]),\n )\n print(distances, indices)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.350103Z","iopub.execute_input":"2024-03-27T19:26:20.350361Z","iopub.status.idle":"2024-03-27T19:26:20.363754Z","shell.execute_reply.started":"2024-03-27T19:26:20.350339Z","shell.execute_reply":"2024-03-27T19:26:20.362738Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def keypoint_distances(\n paths: list[Path],\n index_pairs: list[tuple[int, int]],\n feature_dir: Path,\n min_matches: int = 15,\n verbose: bool = True,\n device: torch.device = torch.device(\"cpu\"),\n) -> None:\n \"\"\"Computes distances between keypoints of images.\n \n Stores output at feature_dir/matches.h5\n \"\"\"\n \n matcher_params = {\n \"width_confidence\": -1,\n \"depth_confidence\": -1,\n \"mp\": True if 'cuda' in str(device) else False,\n }\n matcher = KF.LightGlueMatcher(\"aliked\", matcher_params).eval().to(device)\n \n with h5py.File(feature_dir / \"keypoints.h5\", mode=\"r\") as f_keypoints, \\\n h5py.File(feature_dir / \"descriptors.h5\", mode=\"r\") as f_descriptors, \\\n h5py.File(feature_dir / \"matches.h5\", mode=\"w\") as f_matches:\n \n for idx1, idx2 in tqdm(index_pairs, desc=\"Computing keypoing distances\"):\n key1, key2 = paths[idx1].name, paths[idx2].name\n\n keypoints1 = torch.from_numpy(f_keypoints[key1][...]).to(device)\n keypoints2 = torch.from_numpy(f_keypoints[key2][...]).to(device)\n descriptors1 = torch.from_numpy(f_descriptors[key1][...]).to(device)\n descriptors2 = torch.from_numpy(f_descriptors[key2][...]).to(device)\n\n with torch.inference_mode():\n distances, indices = matcher(\n descriptors1, \n descriptors2, \n KF.laf_from_center_scale_ori(keypoints1[None]),\n KF.laf_from_center_scale_ori(keypoints2[None]),\n )\n\n # We have matches to consider\n n_matches = len(indices)\n if n_matches:\n if verbose:\n print(f\"{key1}-{key2}: {n_matches} matches\")\n # Store the matches in the group of one image\n if n_matches >= min_matches:\n group = f_matches.require_group(key1)\n group.create_dataset(key2, data=indices.detach().cpu().numpy().reshape(-1, 2))","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.364842Z","iopub.execute_input":"2024-03-27T19:26:20.365126Z","iopub.status.idle":"2024-03-27T19:26:20.378255Z","shell.execute_reply.started":"2024-03-27T19:26:20.365102Z","shell.execute_reply":"2024-03-27T19:26:20.377301Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n keypoint_distances(images_list, index_pairs, feature_dir, verbose=False)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.379532Z","iopub.execute_input":"2024-03-27T19:26:20.379881Z","iopub.status.idle":"2024-03-27T19:26:20.391846Z","shell.execute_reply.started":"2024-03-27T19:26:20.379837Z","shell.execute_reply":"2024-03-27T19:26:20.390916Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# RANSAC\nUp to now, we have matched keypoints and their descriptors extracted from pairs of images. This is described by a [fundamental matrix](https://en.wikipedia.org/wiki/Fundamental_matrix_(computer_vision)) denoted as $F$. In epipolar geometry, with homogeneous image coordinates, $x$ and $x′$, of corresponding points in a stereo image pair, $Fx$ describes a line (an epipolar line) on which the corresponding point $x′$ on the other image must lie. That means, for all pairs of corresponding points, $x'Fx = 0$ holds. This is known as epipolar constraint or correspondance condition (or Longuet-Higgins equation), and is solved via the [eight-point algorithm](https://en.wikipedia.org/wiki/Eight-point_algorithm).\n\n
\"Fundamental
\n\nSince the keypoint correspondences are computed using feature descriptors, the data is bound to be noisy and (in general) contains several outliers. Thus, to remove these outliers, we use a [RANSAC](https://en.wikipedia.org/wiki/Random_sample_consensus) algorithm to find the best possible fundamental matrix. So, out of all possibilities, the $F$ matrix with maximum number of inliers is chosen.\n\n
\"RANSAC\"
","metadata":{}},{"cell_type":"code","source":"def import_into_colmap(\n path: Path,\n feature_dir: Path,\n database_path: str = \"colmap.db\",\n) -> None:\n \"\"\"Adds keypoints into colmap\"\"\"\n db = COLMAPDatabase.connect(database_path)\n db.create_tables()\n single_camera = False\n fname_to_id = add_keypoints(db, feature_dir, path, \"\", \"simple-pinhole\", single_camera)\n add_matches(\n db,\n feature_dir,\n fname_to_id,\n )\n db.commit()","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.393044Z","iopub.execute_input":"2024-03-27T19:26:20.393344Z","iopub.status.idle":"2024-03-27T19:26:20.402371Z","shell.execute_reply.started":"2024-03-27T19:26:20.393321Z","shell.execute_reply":"2024-03-27T19:26:20.401664Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n database_path = \"colmap.db\"\n images_dir = images_list[0].parent\n import_into_colmap(\n images_dir, \n feature_dir, \n database_path,\n )\n\n # This does RANSAC\n pycolmap.match_exhaustive(database_path)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.40346Z","iopub.execute_input":"2024-03-27T19:26:20.403745Z","iopub.status.idle":"2024-03-27T19:26:20.415434Z","shell.execute_reply.started":"2024-03-27T19:26:20.403721Z","shell.execute_reply":"2024-03-27T19:26:20.414517Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n# Sparse Reconstruction\n\nNow we have similar image pairs, with matched keypoint descriptors, without outliers! All that is left is to construct the scene and obtain the camera positions. We do this with pycolmap, which offers an incremental reconstruction algorithm that starts from two pairs of images and continually adds more and more images to the scene, resulting in a reconstructed scene with camera information. We can then use the camera rotation and translation as our submission!","metadata":{}},{"cell_type":"code","source":"if DEBUG:\n mapper_options = pycolmap.IncrementalPipelineOptions()\n mapper_options.min_model_size = 3\n mapper_options.max_num_models = 2\n\n maps = pycolmap.incremental_mapping(\n database_path=database_path, \n image_path=images_dir,\n output_path=Path.cwd() / \"incremental_pipeline_outputs\", \n options=mapper_options,\n )","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-03-27T19:26:20.418947Z","iopub.execute_input":"2024-03-27T19:26:20.419258Z","iopub.status.idle":"2024-03-27T19:26:20.426206Z","shell.execute_reply.started":"2024-03-27T19:26:20.419234Z","shell.execute_reply":"2024-03-27T19:26:20.425348Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"if DEBUG:\n print(maps[0].summary())\n for k, im in maps[0].images.items():\n print(\"Rotation\", im.cam_from_world.rotation.matrix(), \"Translation:\", im.cam_from_world.translation, sep=\"\\n\")\n print()","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.427189Z","iopub.execute_input":"2024-03-27T19:26:20.427423Z","iopub.status.idle":"2024-03-27T19:26:20.437259Z","shell.execute_reply.started":"2024-03-27T19:26:20.427402Z","shell.execute_reply":"2024-03-27T19:26:20.436438Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Running everything","metadata":{}},{"cell_type":"code","source":"def parse_sample_submission(\n base_path: Path,\n) -> dict[dict[str, list[Path]]]:\n \"\"\"Construct a dict describing the test data as \n \n {\"dataset\": {\"scene\": []}}\n \"\"\"\n data_dict = {}\n with open(base_path / \"sample_submission.csv\", \"r\") as f:\n for i, l in enumerate(f):\n # Skip header\n if i == 0:\n print(\"header:\", l)\n\n if l and i > 0:\n image_path, dataset, scene, _, _ = l.strip().split(',')\n if dataset not in data_dict:\n data_dict[dataset] = {}\n if scene not in data_dict[dataset]:\n data_dict[dataset][scene] = []\n data_dict[dataset][scene].append(Path(base_path / image_path))\n\n for dataset in data_dict:\n for scene in data_dict[dataset]:\n print(f\"{dataset} / {scene} -> {len(data_dict[dataset][scene])} images\")\n\n return data_dict","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-03-27T19:26:20.438391Z","iopub.execute_input":"2024-03-27T19:26:20.438648Z","iopub.status.idle":"2024-03-27T19:26:20.44954Z","shell.execute_reply.started":"2024-03-27T19:26:20.438626Z","shell.execute_reply":"2024-03-27T19:26:20.44869Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def create_submission(\n results: dict,\n data_dict: dict[dict[str, list[Path]]],\n base_path: Path,\n) -> None:\n \"\"\"Prepares a submission file.\"\"\"\n \n with open(\"submission.csv\", \"w\") as f:\n f.write(\"image_path,dataset,scene,rotation_matrix,translation_vector\\n\")\n \n for dataset in data_dict:\n # Only write results for datasets with images that have results \n if dataset in results:\n res = results[dataset]\n else:\n res = {}\n \n # Same for scenes\n for scene in data_dict[dataset]:\n if scene in res:\n scene_res = res[scene]\n else:\n scene_res = {\"R\":{}, \"t\":{}}\n \n # Write the row with rotation and translation matrices\n for image in data_dict[dataset][scene]:\n if image in scene_res:\n print(image)\n R = scene_res[image][\"R\"].reshape(-1)\n T = scene_res[image][\"t\"].reshape(-1)\n else:\n R = np.eye(3).reshape(-1)\n T = np.zeros((3))\n image_path = str(image.relative_to(base_path))\n f.write(f\"{image_path},{dataset},{scene},{arr_to_str(R)},{arr_to_str(T)}\\n\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-03-27T19:26:20.450681Z","iopub.execute_input":"2024-03-27T19:26:20.450974Z","iopub.status.idle":"2024-03-27T19:26:20.464365Z","shell.execute_reply.started":"2024-03-27T19:26:20.450951Z","shell.execute_reply":"2024-03-27T19:26:20.463622Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Config:\n base_path: Path = Path(\"/kaggle/input/image-matching-challenge-2024\")\n feature_dir: Path = Path.cwd() / \".feature_outputs\"\n \n device: torch.device = K.utils.get_cuda_device_if_available(0)\n \n pair_matching_args = {\n \"model_name\": \"/kaggle/input/dinov2/pytorch/base/1\",\n \"similarity_threshold\": 0.3,\n \"tolerance\": 500,\n \"min_matches\": 50,\n \"exhaustive_if_less\": 50,\n \"p\": 2.0,\n }\n \n keypoint_detection_args = {\n \"num_features\": 4096,\n \"resize_to\": 1024,\n }\n \n keypoint_distances_args = {\n \"min_matches\": 15,\n \"verbose\": False,\n }\n \n colmap_mapper_options = {\n \"min_model_size\": 3, # By default colmap does not generate a reconstruction if less than 10 images are registered. Lower it to 3.\n \"max_num_models\": 2,\n }","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.465376Z","iopub.execute_input":"2024-03-27T19:26:20.466245Z","iopub.status.idle":"2024-03-27T19:26:20.478884Z","shell.execute_reply.started":"2024-03-27T19:26:20.466214Z","shell.execute_reply":"2024-03-27T19:26:20.478113Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def run_from_config(config: Config) -> None:\n results = {}\n \n data_dict = parse_sample_submission(config.base_path)\n datasets = list(data_dict.keys())\n \n for dataset in datasets:\n if dataset not in results:\n results[dataset] = {}\n \n for scene in data_dict[dataset]:\n images_dir = data_dict[dataset][scene][0].parent\n results[dataset][scene] = {}\n image_paths = data_dict[dataset][scene]\n print (f\"Got {len(image_paths)} images\")\n \n try:\n feature_dir = config.feature_dir / f\"{dataset}_{scene}\"\n feature_dir.mkdir(parents=True, exist_ok=True)\n database_path = feature_dir / \"colmap.db\"\n if database_path.exists():\n database_path.unlink()\n \n # 1. Get the pairs of images that are somewhat similar\n index_pairs = get_image_pairs(\n image_paths,\n **config.pair_matching_args,\n device=config.device,\n )\n gc.collect()\n \n # 2. Detect keypoints of all images\n detect_keypoints(\n image_paths,\n feature_dir,\n **config.keypoint_detection_args,\n device=device,\n )\n gc.collect()\n \n # 3. Match keypoints of pairs of similar images\n keypoint_distances(\n image_paths, \n index_pairs, \n feature_dir,\n **config.keypoint_distances_args,\n device=device,\n )\n gc.collect()\n \n sleep(1)\n \n # 4.1. Import keypoint distances of matches into colmap for RANSAC \n import_into_colmap(\n images_dir, \n feature_dir, \n database_path,\n )\n \n output_path = feature_dir / \"colmap_rec_aliked\"\n output_path.mkdir(parents=True, exist_ok=True)\n \n # 4.2. Compute RANSAC (detect match outliers)\n # By doing it exhaustively we guarantee we will find the best possible configuration\n pycolmap.match_exhaustive(database_path)\n \n mapper_options = pycolmap.IncrementalPipelineOptions(**config.colmap_mapper_options)\n \n # 5.1 Incrementally start reconstructing the scene (sparse reconstruction)\n # The process starts from a random pair of images and is incrementally extended by \n # registering new images and triangulating new points.\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 )\n \n print(maps)\n clear_output(wait=False)\n \n # 5.2. Look for the best reconstruction: The incremental mapping offered by \n # pycolmap attempts to reconstruct multiple models, we must pick the best one\n images_registered = 0\n best_idx = None\n \n print (\"Looking for the best reconstruction\")\n \n if isinstance(maps, dict):\n for idx1, rec in maps.items():\n print(idx1, rec.summary())\n try:\n if len(rec.images) > images_registered:\n images_registered = len(rec.images)\n best_idx = idx1\n except Exception:\n continue\n \n # Parse the reconstruction object to get the rotation matrix and translation vector\n # obtained for each image in the reconstruction\n if best_idx is not None:\n for k, im in maps[best_idx].images.items():\n key = config.base_path / \"test\" / scene / \"images\" / im.name\n results[dataset][scene][key] = {}\n results[dataset][scene][key][\"R\"] = deepcopy(im.cam_from_world.rotation.matrix())\n results[dataset][scene][key][\"t\"] = deepcopy(np.array(im.cam_from_world.translation))\n \n print(f\"Registered: {dataset} / {scene} -> {len(results[dataset][scene])} images\")\n print(f\"Total: {dataset} / {scene} -> {len(data_dict[dataset][scene])} images\")\n create_submission(results, data_dict, config.base_path)\n gc.collect()\n \n except Exception as e:\n print(e)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.480193Z","iopub.execute_input":"2024-03-27T19:26:20.480514Z","iopub.status.idle":"2024-03-27T19:26:20.499335Z","shell.execute_reply.started":"2024-03-27T19:26:20.480486Z","shell.execute_reply":"2024-03-27T19:26:20.498449Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"run_from_config(Config)","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:26:20.500314Z","iopub.execute_input":"2024-03-27T19:26:20.500575Z","iopub.status.idle":"2024-03-27T19:31:55.75141Z","shell.execute_reply.started":"2024-03-27T19:26:20.500542Z","shell.execute_reply":"2024-03-27T19:31:55.750638Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!cat submission.csv","metadata":{"execution":{"iopub.status.busy":"2024-03-27T19:31:55.752732Z","iopub.execute_input":"2024-03-27T19:31:55.753102Z","iopub.status.idle":"2024-03-27T19:31:56.760105Z","shell.execute_reply.started":"2024-03-27T19:31:55.753069Z","shell.execute_reply":"2024-03-27T19:31:56.758889Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# What to do next?\n\nHere are some ways in which you can explore potential improvements:\n\n- Using a different image embedding model to obtain the image pairs\n- Trying other approaches for keypoint extraction, such as SIFT or DISK\n- Leveraging the training data to train a better models for each dataset","metadata":{}}]}