{"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":"code","source":"import numpy as np\nfrom dataclasses import dataclass\nfrom time import time","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-04-11T17:10:31.618334Z","iopub.execute_input":"2023-04-11T17:10:31.618764Z","iopub.status.idle":"2023-04-11T17:10:31.624886Z","shell.execute_reply.started":"2023-04-11T17:10:31.618726Z","shell.execute_reply":"2023-04-11T17:10:31.623294Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Conveniency functions.\n\ndef arr_to_str(a):\n    return ';'.join([str(x) for x in a.reshape(-1)])","metadata":{"execution":{"iopub.status.busy":"2023-04-11T17:10:31.819452Z","iopub.execute_input":"2023-04-11T17:10:31.819889Z","iopub.status.idle":"2023-04-11T17:10:31.825531Z","shell.execute_reply.started":"2023-04-11T17:10:31.819844Z","shell.execute_reply":"2023-04-11T17:10:31.824250Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Evaluation metric.\n\n@dataclass\nclass Camera:\n    rotmat: np.array\n    tvec: np.array\n\ndef quaternion_from_matrix(matrix):\n    M = np.array(matrix, dtype=np.float64, copy=False)[:4, :4]\n    m00 = M[0, 0]\n    m01 = M[0, 1]\n    m02 = M[0, 2]\n    m10 = M[1, 0]\n    m11 = M[1, 1]\n    m12 = M[1, 2]\n    m20 = M[2, 0]\n    m21 = M[2, 1]\n    m22 = M[2, 2]\n\n    # Symmetric matrix K.\n    K = np.array([[m00 - m11 - m22, 0.0, 0.0, 0.0],\n                  [m01 + m10, m11 - m00 - m22, 0.0, 0.0],\n                  [m02 + m20, m12 + m21, m22 - m00 - m11, 0.0],\n                  [m21 - m12, m02 - m20, m10 - m01, m00 + m11 + m22]])\n    K /= 3.0\n\n    # Quaternion is eigenvector of K that corresponds to largest eigenvalue.\n    w, V = np.linalg.eigh(K)\n    q = V[[3, 0, 1, 2], np.argmax(w)]\n\n    if q[0] < 0.0:\n        np.negative(q, q)\n    return q\n\ndef evaluate_R_t(R_gt, t_gt, R, t, eps=1e-15):\n    t = t.flatten()\n    t_gt = t_gt.flatten()\n\n    q_gt = quaternion_from_matrix(R_gt)\n    q = quaternion_from_matrix(R)\n    q = q / (np.linalg.norm(q) + eps)\n    q_gt = q_gt / (np.linalg.norm(q_gt) + eps)\n    loss_q = np.maximum(eps, (1.0 - np.sum(q * q_gt)**2))\n    err_q = np.arccos(1 - 2 * loss_q)\n\n    GT_SCALE = np.linalg.norm(t_gt)\n    t = GT_SCALE * (t / (np.linalg.norm(t) + eps))\n    err_t = min(np.linalg.norm(t_gt - t), np.linalg.norm(t_gt + t))\n    \n    return np.degrees(err_q), err_t\n\ndef compute_dR_dT(R1, T1, R2, T2):\n    '''Given absolute (R, T) pairs for two cameras, compute the relative pose difference, from the first.'''\n    \n    dR = np.dot(R2, R1.T)\n    dT = T2 - np.dot(dR, T1)\n    return dR, dT\n\ndef compute_mAA(err_q, err_t, ths_q, ths_t):\n    '''Compute the mean average accuracy over a set of thresholds. Additionally returns the metric only over rotation and translation.'''\n\n    acc, acc_q, acc_t = [], [], []\n    for th_q, th_t in zip(ths_q, ths_t):\n        cur_acc_q = (err_q <= th_q)\n        cur_acc_t = (err_t <= th_t)\n        cur_acc = cur_acc_q & cur_acc_t\n        \n        acc.append(cur_acc.astype(np.float32).mean())\n        acc_q.append(cur_acc_q.astype(np.float32).mean())\n        acc_t.append(cur_acc_t.astype(np.float32).mean())\n    return np.array(acc), np.array(acc_q), np.array(acc_t)\n\ndef dict_from_csv(csv_path, has_header):\n    csv_dict = {}\n    with open(csv_path, 'r') as f:\n        for i, l in enumerate(f):\n            if has_header and i == 0:\n                continue\n            if l:\n                image, dataset, scene, R_str, T_str = l.strip().split(',')\n                R = np.fromstring(R_str.strip(), sep=';').reshape(3, 3)\n                T = np.fromstring(T_str.strip(), sep=';')\n                if dataset not in csv_dict:\n                    csv_dict[dataset] = {}\n                if scene not in csv_dict[dataset]:\n                    csv_dict[dataset][scene] = {}\n                csv_dict[dataset][scene][image] = Camera(rotmat=R, tvec=T)\n    return csv_dict\n\ndef eval_submission(submission_csv_path, ground_truth_csv_path, rotation_thresholds_degrees_dict, translation_thresholds_meters_dict, verbose=False):\n    '''Compute final metric given submission and ground truth files. Thresholds are specified per dataset.'''\n\n    submission_dict = dict_from_csv(submission_csv_path, has_header=True)\n    gt_dict = dict_from_csv(ground_truth_csv_path, has_header=True)\n\n    # Check that all necessary keys exist in the submission file\n    for dataset in gt_dict:\n        assert dataset in submission_dict, f'Unknown dataset: {dataset}'\n        for scene in gt_dict[dataset]:\n            assert scene in submission_dict[dataset], f'Unknown scene: {dataset}->{scene}'\n            for image in gt_dict[dataset][scene]:\n                assert image in submission_dict[dataset][scene], f'Unknown image: {dataset}->{scene}->{image}'\n\n    # Iterate over all the scenes\n    if verbose:\n        t = time()\n        print('*** METRICS ***')\n\n    metrics_per_dataset = []\n    for dataset in gt_dict:\n        metrics_per_scene = []\n        for scene in gt_dict[dataset]:\n            err_q_all = []\n            err_t_all = []\n            images = [camera for camera in gt_dict[dataset][scene]]\n            # Process all pairs in a scene\n            for i in range(len(images)):\n                for j in range(i + 1, len(images)):\n                    gt_i = gt_dict[dataset][scene][images[i]]\n                    gt_j = gt_dict[dataset][scene][images[j]]\n                    dR_gt, dT_gt = compute_dR_dT(gt_i.rotmat, gt_i.tvec, gt_j.rotmat, gt_j.tvec)\n\n                    pred_i = submission_dict[dataset][scene][images[i]]\n                    pred_j = submission_dict[dataset][scene][images[j]]\n                    dR_pred, dT_pred = compute_dR_dT(pred_i.rotmat, pred_i.tvec, pred_j.rotmat, pred_j.tvec)\n\n                    err_q, err_t = evaluate_R_t(dR_gt, dT_gt, dR_pred, dT_pred)\n                    err_q_all.append(err_q)\n                    err_t_all.append(err_t)\n\n            mAA, mAA_q, mAA_t = compute_mAA(err_q=err_q_all,\n                                            err_t=err_t_all,\n                                            ths_q=rotation_thresholds_degrees_dict[(dataset, scene)],\n                                            ths_t=translation_thresholds_meters_dict[(dataset, scene)])\n            if verbose:\n                print(f'{dataset} / {scene} ({len(images)} images, {len(err_q_all)} pairs) -> mAA={np.mean(mAA):.06f}, mAA_q={np.mean(mAA_q):.06f}, mAA_t={np.mean(mAA_t):.06f}')\n            metrics_per_scene.append(np.mean(mAA))\n\n        metrics_per_dataset.append(np.mean(metrics_per_scene))\n        if verbose:\n            print(f'{dataset} -> mAA={np.mean(metrics_per_scene):.06f}')\n            print()\n\n    if verbose:\n        print(f'Final metric -> mAA={np.mean(metrics_per_dataset):.06f} (t: {time() - t} sec.)')\n        print()\n\n    return np.mean(metrics_per_dataset)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T17:11:58.758679Z","iopub.execute_input":"2023-04-11T17:11:58.759336Z","iopub.status.idle":"2023-04-11T17:11:58.794212Z","shell.execute_reply.started":"2023-04-11T17:11:58.759284Z","shell.execute_reply":"2023-04-11T17:11:58.792957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Set rotation thresholds per scene.\n\nrotation_thresholds_degrees_dict = {\n    **{('haiper', scene): np.linspace(1, 10, 10) for scene in ['bike', 'chairs', 'fountain']},\n    **{('heritage', scene): np.linspace(1, 10, 10) for scene in ['cyprus', 'dioscuri']},\n    **{('heritage', 'wall'): np.linspace(0.2, 10, 10)},\n    **{('urban', 'kyiv-puppet-theater'): np.linspace(1, 10, 10)},\n}\n\ntranslation_thresholds_meters_dict = {\n    **{('haiper', scene): np.geomspace(0.05, 0.5, 10) for scene in ['bike', 'chairs', 'fountain']},\n    **{('heritage', scene): np.geomspace(0.1, 2, 10) for scene in ['cyprus', 'dioscuri']},\n    **{('heritage', 'wall'): np.geomspace(0.05, 1, 10)},\n    **{('urban', 'kyiv-puppet-theater'): np.geomspace(0.5, 5, 10)},\n}","metadata":{"execution":{"iopub.status.busy":"2023-04-11T17:41:33.577788Z","iopub.execute_input":"2023-04-11T17:41:33.578177Z","iopub.status.idle":"2023-04-11T17:41:33.590549Z","shell.execute_reply.started":"2023-04-11T17:41:33.578142Z","shell.execute_reply":"2023-04-11T17:41:33.589151Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate and evaluate a random submission.\n\nsrc = '/kaggle/input/image-matching-challenge-2023'\n\n# Note that the fields were reordered between dataset versions so you need to be careful.\n# Here we regenerate the ground truth file using the submission format, which lists the image path first.\nwith open(f'{src}/train/train_labels.csv', 'r') as fr, open('submission.csv', 'w') as fw:\n    for i, l in enumerate(fr):\n        if i == 0:\n            fw.write('image_path,dataset,scene,rotation_matrix,translation_vector\\n')\n        else:\n            dataset, scene, image, _, _ = l.strip().split(',')\n            R = np.random.rand(9)\n            T = np.random.rand(3)\n            fw.write(f'{image},{dataset},{scene},{arr_to_str(R)},{arr_to_str(T)}\\n')\n\nwith open(f'{src}/train/train_labels.csv', 'r') as fr, open('ground_truth.csv', 'w') as fw:\n    for i, l in enumerate(fr):\n        if i == 0:\n            fw.write('image_path,dataset,scene,rotation_matrix,translation_vector\\n')\n        else:\n            dataset, scene, image, R, T = l.strip().split(',')\n            fw.write(f'{image},{dataset},{scene},{R},{T}\\n')\n\neval_submission(submission_csv_path='submission.csv',\n                ground_truth_csv_path='ground_truth.csv',\n                rotation_thresholds_degrees_dict=rotation_thresholds_degrees_dict,\n                translation_thresholds_meters_dict=translation_thresholds_meters_dict,\n                verbose=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T17:41:35.355490Z","iopub.execute_input":"2023-04-11T17:41:35.355899Z","iopub.status.idle":"2023-04-11T17:41:38.075517Z","shell.execute_reply.started":"2023-04-11T17:41:35.355861Z","shell.execute_reply":"2023-04-11T17:41:38.074393Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now evaluate a perfect submission.\n\neval_submission(submission_csv_path='ground_truth.csv',\n                ground_truth_csv_path='ground_truth.csv',\n                rotation_thresholds_degrees_dict=rotation_thresholds_degrees_dict,\n                translation_thresholds_meters_dict=translation_thresholds_meters_dict,\n                verbose=True)","metadata":{"execution":{"iopub.status.busy":"2023-04-11T17:42:02.345139Z","iopub.execute_input":"2023-04-11T17:42:02.345545Z","iopub.status.idle":"2023-04-11T17:42:05.043693Z","shell.execute_reply.started":"2023-04-11T17:42:02.345512Z","shell.execute_reply":"2023-04-11T17:42:05.042393Z"},"trusted":true},"execution_count":null,"outputs":[]}]}