{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Image matching Challenge\n\nAs a robotics engineer 3D reconstruction is an interesting problem.\nEspecially one of my interest, as the perception side is my focus of work","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"mode = \"development\"\nmode = \"submission\"","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:02.998045Z","iopub.execute_input":"2023-05-16T19:35:02.998479Z","iopub.status.idle":"2023-05-16T19:35:03.035036Z","shell.execute_reply.started":"2023-05-16T19:35:02.998440Z","shell.execute_reply":"2023-05-16T19:35:03.033573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if mode != \"submission\":\n    !pip install open3d","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:03.037400Z","iopub.execute_input":"2023-05-16T19:35:03.037894Z","iopub.status.idle":"2023-05-16T19:35:03.044060Z","shell.execute_reply.started":"2023-05-16T19:35:03.037852Z","shell.execute_reply":"2023-05-16T19:35:03.042739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport cv2\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nimport pycolmap\nif mode != \"submission\":\n    import open3d as o3d\nfrom tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:03.046196Z","iopub.execute_input":"2023-05-16T19:35:03.046678Z","iopub.status.idle":"2023-05-16T19:35:05.073948Z","shell.execute_reply.started":"2023-05-16T19:35:03.046638Z","shell.execute_reply":"2023-05-16T19:35:05.072355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Discover a single object\nLets start by discovering the images and 3D reconstruction of the haiper bike","metadata":{}},{"cell_type":"code","source":"category = \"phototourism\"\nobject_name = \"st_peters_square\"\npath_object = f\"/kaggle/input/image-matching-challenge-2023/train/{category}/{object_name}/images/\"\nimage_files = os.listdir(path_object)\nlen(image_files)","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:05.078162Z","iopub.execute_input":"2023-05-16T19:35:05.080864Z","iopub.status.idle":"2023-05-16T19:35:05.230673Z","shell.execute_reply.started":"2023-05-16T19:35:05.080798Z","shell.execute_reply":"2023-05-16T19:35:05.229451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## looking at the images","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(16,12))\nf, axarr = plt.subplots(3, 5) \n\ni = 0\nj = 0\n\nfor file in image_files[:15]:\n    im = plt.imread(path_object+file, format='jpeg')\n    dim = (im.shape[1] // 2, im.shape[0] // 2)\n    im = cv2.resize(im, dim, interpolation = cv2.INTER_AREA)\n    axarr[i,j].imshow(im)\n    \n    i += 1\n    if i == 3:\n        j += 1\n        i = 0\n","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:05.234541Z","iopub.execute_input":"2023-05-16T19:35:05.235097Z","iopub.status.idle":"2023-05-16T19:35:08.052861Z","shell.execute_reply.started":"2023-05-16T19:35:05.235065Z","shell.execute_reply":"2023-05-16T19:35:08.051183Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## Select two similair images\nimages_compare = [image_files[3], image_files[12]]","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.054533Z","iopub.execute_input":"2023-05-16T19:35:08.054943Z","iopub.status.idle":"2023-05-16T19:35:08.062872Z","shell.execute_reply.started":"2023-05-16T19:35:08.054910Z","shell.execute_reply":"2023-05-16T19:35:08.060821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## orientation of every image","metadata":{}},{"cell_type":"markdown","source":"The SfM files are produced with COLMAP so the files are in this order\n\nLets use pycolmap to extract features and summarize them\n\nhttps://colmap.github.io/index.html#","metadata":{}},{"cell_type":"code","source":"if mode != \"submission\":\n    reconstruction_dir = f\"/kaggle/input/image-matching-challenge-2023/train/{category}/{object_name}/sfm\"\n    images_dir = f\"/kaggle/input/image-matching-challenge-2023/train/{category}/{object_name}/images\"\n\n    reconstruction = pycolmap.Reconstruction(reconstruction_dir)\n    print(reconstruction.summary())\n","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.064735Z","iopub.execute_input":"2023-05-16T19:35:08.065790Z","iopub.status.idle":"2023-05-16T19:35:08.085667Z","shell.execute_reply.started":"2023-05-16T19:35:08.065745Z","shell.execute_reply":"2023-05-16T19:35:08.084740Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if mode != \"submission\":\n\n    for image_id, image in reconstruction.images.items():\n        print(image.summary())\n        print(image.rotation_matrix())\n        break\n\n    def get_images_from_reconstruction(reconstruction):\n        images = []\n        for image_id, image in reconstruction.images.items():\n            images.append([image_id, image.camera_id, image.name, image.rotation_matrix(), image.tvec])\n\n        df = pd.DataFrame(np.array(images), columns=[\"image id\", \"camera_id\", \"image name\", \"rotation_matrix\", \"transformation\"])\n        return df\n\n\n    df_images = get_images_from_reconstruction(reconstruction)\n    df_images","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.087005Z","iopub.execute_input":"2023-05-16T19:35:08.088179Z","iopub.status.idle":"2023-05-16T19:35:08.099215Z","shell.execute_reply.started":"2023-05-16T19:35:08.088139Z","shell.execute_reply":"2023-05-16T19:35:08.098247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if mode != \"submission\":\n    for point3D_id, point3D in reconstruction.points3D.items():\n        print(point3D.summary())\n        break\n\n    def get_points_from_reconstruction(reconstruction):\n        points = []\n        colors = []\n        for point3D_id, point3D in reconstruction.points3D.items():\n            points.append(point3D.xyz)\n            colors.append(point3D.color)\n\n        points = np.array(points)\n        colors = np.array(colors)\n        return points, colors\n\n    points, colors = get_points_from_reconstruction(reconstruction)\n\n\n    pcd = o3d.geometry.PointCloud()\n    pcd.points = o3d.utility.Vector3dVector(points)\n    pcd.colors = o3d.utility.Vector3dVector(colors)\n    o3d.visualization.draw_plotly([pcd])","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.100851Z","iopub.execute_input":"2023-05-16T19:35:08.101951Z","iopub.status.idle":"2023-05-16T19:35:08.116958Z","shell.execute_reply.started":"2023-05-16T19:35:08.101915Z","shell.execute_reply":"2023-05-16T19:35:08.116065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if mode != \"submission\":\n    i = 0\n    for camera_id, camera in reconstruction.cameras.items():\n        print(camera.summary())\n        break\n\n    def get_cameras_from_reconstruction(reconstruction):\n        cameras = []\n        for camera_id, camera in reconstruction.cameras.items():\n            cameras.append([camera_id, camera.model_name, camera.width, camera.height, camera.params[0], camera.params[1], camera.params[2], camera.params[3]])\n\n        df = pd.DataFrame(np.array(cameras), columns=[\"camera_id\", \"camera model\", \"width\", \"height\", \"Focal length\", \"Optical center X\", \"Optical center Y\", \"K?\"])\n        return df\n\n\n    df_cameras = get_cameras_from_reconstruction(reconstruction)\n    df_cameras","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.121862Z","iopub.execute_input":"2023-05-16T19:35:08.122449Z","iopub.status.idle":"2023-05-16T19:35:08.132332Z","shell.execute_reply.started":"2023-05-16T19:35:08.122411Z","shell.execute_reply":"2023-05-16T19:35:08.130840Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we have a function to extract the data from a reconstruction for images, points and cameras.","metadata":{}},{"cell_type":"markdown","source":"# Structure from Motion (SfM) methods\nColmap is the most well know method for SfM, lets see if we can reproduce some results from the training data.\n\npycolmap will be used to create an understanding of the selected object\n\nFirst SIFT features are extracted for the images\n","metadata":{}},{"cell_type":"code","source":"from PIL import Image, ImageOps\n\ndescriptor_list = []\nkeypoint_list = []\nimgs = []\nmax_features = 512\n\nfor file in images_compare:\n    im = Image.open(path_object+file).convert('RGB')\n    img = ImageOps.grayscale(im)\n    img = np.array(img).astype(np.float) / 255.\n\n    # Optional parameters:\n    # - options: dict or pycolmap.SiftExtractionOptions\n    # - device: default pycolmap.Device.auto uses the GPU if available\n    sift = pycolmap.Sift()\n\n    # Parameters:\n    # - image: HxW float array\n    keypoints, scores, descriptors = sift.extract(img)\n    scores_idx = np.argsort(scores)\n    descriptor_list.append(descriptors[scores_idx[-(max_features):]])\n    keypoint_list.append(keypoints[scores_idx[-(max_features):]])\n\n    im = np.array(im)\n    im_append = im.copy()\n    imgs.append(im_append)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:08.134440Z","iopub.execute_input":"2023-05-16T19:35:08.134884Z","iopub.status.idle":"2023-05-16T19:35:09.457155Z","shell.execute_reply.started":"2023-05-16T19:35:08.134844Z","shell.execute_reply":"2023-05-16T19:35:09.456036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Calculating the loss between selected images","metadata":{}},{"cell_type":"code","source":"print(descriptor_list[0].shape)\nprint(descriptor_list[1].shape)\n\n\nscore_matrix = np.zeros((descriptor_list[0].shape[0], descriptor_list[1].shape[0]))\n\nmatching_keypoints = []\nerrors = []\nerror_threshold = 0.1\n\nfrom tqdm import tqdm\nfor i in tqdm(range(score_matrix.shape[0])):\n    for j in range(score_matrix.shape[1]):\n        if abs(np.sum((descriptor_list[0][i] - descriptor_list[1][j])**2)) < error_threshold:\n            matching_keypoints.append([i,j])\n\nprint(len(matching_keypoints), score_matrix.shape[0]*score_matrix.shape[1])","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:09.458865Z","iopub.execute_input":"2023-05-16T19:35:09.459348Z","iopub.status.idle":"2023-05-16T19:35:13.357114Z","shell.execute_reply.started":"2023-05-16T19:35:09.459304Z","shell.execute_reply":"2023-05-16T19:35:13.355971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## plotting the similair pixels\n","metadata":{}},{"cell_type":"code","source":"# Radius of circle\nradius = 5\n# Line thickness of 2 px\nthickness = 2\n\ncolormaps = mpl.colormaps['tab20c']\nnorm = mpl.colors.Normalize(vmin=0, vmax=(len(matching_keypoints)-1))\n\nif imgs[0].shape[0] > imgs[1].shape[0]:\n    add = (imgs[0].shape[0] - imgs[1].shape[0])\n    pad = np.zeros((add, imgs[1].shape[1], 3))*255\n    imgs[1] = np.append(imgs[1], pad, axis=0)\nelif imgs[1].shape[0] > imgs[0].shape[0]:\n    add = (imgs[1].shape[0] - imgs[0].shape[0])\n    pad = np.zeros((add, imgs[0].shape[1], 3))*255\n    imgs[0] = np.append(imgs[0], pad, axis=0)\n\n\ncombined_image = np.append(imgs[0],imgs[1], axis=1).astype(np.int32)\nwidth = imgs[0].shape[1]\n\ni=0\nfor point in matching_keypoints:\n    im1_point_id = descriptor_list[0][point[0]]\n    im2_point_id = descriptor_list[1][point[1]]\n    im1_point = np.array([int(keypoint_list[0][point[0]][0]), int(keypoint_list[0][point[0]][1])])\n    im2_point = np.array([int(keypoint_list[1][point[1]][0]), int(keypoint_list[1][point[1]][1])])\n    im2_point[0] += width\n\n    color = colormaps(norm(i))\n    color = (int(color[0]*255), int(color[1]*255), int(color[2]*255))\n\n    cv2.circle(combined_image, (im1_point[0], im1_point[1]), radius, color, thickness)\n    cv2.circle(combined_image, (im2_point[0], im2_point[1]), radius, color, thickness)\n    cv2.line(combined_image, (im1_point[0], im1_point[1]), (int(im2_point[0]), int(im2_point[1])), color, thickness)\n    i += 1\n\nplt.figure(figsize=(16,8))\nplt.imshow(combined_image)","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:13.358759Z","iopub.execute_input":"2023-05-16T19:35:13.362860Z","iopub.status.idle":"2023-05-16T19:35:14.443739Z","shell.execute_reply.started":"2023-05-16T19:35:13.362821Z","shell.execute_reply":"2023-05-16T19:35:14.442481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3d points estimation\n\nread the paper: https://demuc.de/colmap/","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission Code\n\nAlready setup the submission code, so after we can focus on the COLMAP method.","metadata":{}},{"cell_type":"code","source":"sample_submission = pd.read_csv(\"/kaggle/input/image-matching-challenge-2023/sample_submission.csv\")\nsample_submission","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.445289Z","iopub.execute_input":"2023-05-16T19:35:14.445848Z","iopub.status.idle":"2023-05-16T19:35:14.495455Z","shell.execute_reply.started":"2023-05-16T19:35:14.445809Z","shell.execute_reply":"2023-05-16T19:35:14.494349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_submission[\"dataset/scene\"] = sample_submission[\"dataset\"] + \"/\" + sample_submission[\"scene\"]\n\nall_scenes = sample_submission[\"dataset/scene\"].unique()\nall_scenes","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.496705Z","iopub.execute_input":"2023-05-16T19:35:14.497010Z","iopub.status.idle":"2023-05-16T19:35:14.510306Z","shell.execute_reply.started":"2023-05-16T19:35:14.496984Z","shell.execute_reply":"2023-05-16T19:35:14.508920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.spatial.transform import Rotation as R\n\ndef get_rotation_matrix(roll, pitch, yaw):\n    rotation_matrix = R.from_euler('zyx', [yaw, pitch, roll], degrees=True).as_matrix()\n    rotation_matrix_string = \"\"\n    for i in range(3):\n        for j in range(3):\n            # LAST EDIT, NO semi colon in the last number, test tomorrow 16 May 2023\n            if j == 2 and i == 2:\n                rotation_matrix_string += f\"{rotation_matrix[i,j]}\"\n            else:\n                rotation_matrix_string += f\"{rotation_matrix[i,j]};\"\n    \n    return rotation_matrix_string\n    ","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.511869Z","iopub.execute_input":"2023-05-16T19:35:14.512568Z","iopub.status.idle":"2023-05-16T19:35:14.520936Z","shell.execute_reply.started":"2023-05-16T19:35:14.512536Z","shell.execute_reply":"2023-05-16T19:35:14.519727Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\n\ndef get_all_images(scene_directory):\n    sub_df = sample_submission[sample_submission[\"dataset/scene\"] == scene_directory]\n    images = sub_df[\"image_path\"].values\n    return images\n\ndef analyze_scenes(scene_directory):\n    images = get_all_images(scene_directory)\n    scene_dict = {\n        \"image_path\": [],\n        \"dataset\": [],\n        \"scene\": [],\n        \"rotation_matrix\": [],\n        \"translation_vector\": [],\n    }\n    dataset = scene_directory.split(\"/\")[-2]\n    scene = scene_directory.split(\"/\")[-1]\n    \n    translation_vector = \"0.1;0.2;0.3\"\n    \n    for image in images:\n        roll = random.randint(-10, 10)        \n        pitch = random.randint(-30, 30)\n        yaw = random.randint(-180, 180)\n        rotation_matrix_string = get_rotation_matrix(roll, pitch, yaw)\n        \n        scene_dict[\"image_path\"].append(image)\n        scene_dict[\"dataset\"].append(dataset)\n        scene_dict[\"scene\"].append(scene)\n        scene_dict[\"rotation_matrix\"].append(rotation_matrix_string)\n        scene_dict[\"translation_vector\"].append(translation_vector)\n       \n    df = pd.DataFrame(scene_dict)\n    return df\n        ","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.522558Z","iopub.execute_input":"2023-05-16T19:35:14.522873Z","iopub.status.idle":"2023-05-16T19:35:14.535428Z","shell.execute_reply.started":"2023-05-16T19:35:14.522845Z","shell.execute_reply":"2023-05-16T19:35:14.534028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df = pd.DataFrame(columns=[\"image_path\", \"dataset\", \"scene\", \"rotation_matrix\", \"translation_vector\"])\nprint(submission_df)\n\nfor scene in tqdm(all_scenes):\n    scene_dict = analyze_scenes(scene)\n    submission_df = pd.concat([submission_df, scene_dict])\n    \nsubmission_df","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.537081Z","iopub.execute_input":"2023-05-16T19:35:14.537448Z","iopub.status.idle":"2023-05-16T19:35:14.575224Z","shell.execute_reply.started":"2023-05-16T19:35:14.537420Z","shell.execute_reply":"2023-05-16T19:35:14.574131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2023-05-16T19:35:14.576607Z","iopub.execute_input":"2023-05-16T19:35:14.577447Z","iopub.status.idle":"2023-05-16T19:35:14.588332Z","shell.execute_reply.started":"2023-05-16T19:35:14.577409Z","shell.execute_reply":"2023-05-16T19:35:14.587099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}