diff --git a/.docker/Dockerfile b/.docker/Dockerfile
index 59be9bd..2b2bac0 100644
--- a/.docker/Dockerfile
+++ b/.docker/Dockerfile
@@ -122,9 +122,23 @@ RUN git clone https://github.com/colmap/colmap.git --branch 3.8 ${HOME_DIR}/colm
 ####################
 
 FROM colmap as nerfstudio
+ARG NERFSTUDIO_VERSION=v1.1.4
 RUN pip install --no-cache-dir --upgrade pip
 RUN pip install --no-cache-dir torch==2.1.2+cu118 torchvision==0.16.2+cu118 'numpy<2.0.0' --extra-index-url https://download.pytorch.org/whl/cu118
 RUN TCNN_CUDA_ARCHITECTURES="${CUDA_ARCHITECTURES}" pip install --no-cache-dir "git+https://github.com/NVlabs/tiny-cuda-nn.git@b3473c81396fe927293bdfd5a6be32df8769927c#subdirectory=bindings/torch"
+RUN git clone https://github.com/nerfstudio-project/nerfstudio.git --branch ${NERFSTUDIO_VERSION} ${HOME_DIR}/nerfstudio &&\
+    cd ${HOME_DIR}/nerfstudio &&\
+    pip install -e .
+
+FROM colmap as sdfstudio
+ARG SDFSTUDIO_COMMIT=370902a
+RUN pip install --no-cache-dir --upgrade pip
+RUN pip install torch==1.12.1+cu113 torchvision==0.13.1+cu113 -f https://download.pytorch.org/whl/torch_stable.html
+RUN TCNN_CUDA_ARCHITECTURES="${CUDA_ARCHITECTURES}" pip install --no-cache-dir "git+https://github.com/NVlabs/tiny-cuda-nn.git@b3473c81396fe927293bdfd5a6be32df8769927c#subdirectory=bindings/torch"
+RUN git clone https://github.com/autonomousvision/sdfstudio.git ${HOME_DIR}/sdfstudio &&\
+    cd ${HOME_DIR}/sdfstudio &&\
+    git checkout ${SDFSTUDIO_COMMIT} &&\
+    pip install -e .
 
 ####################
 # Deployment image #
diff --git a/config/recon_benchmark.yaml b/config/recon_benchmark.yaml
new file mode 100644
index 0000000..fd130eb
--- /dev/null
+++ b/config/recon_benchmark.yaml
@@ -0,0 +1,13 @@
+reconstruction_benchmark:
+  project_folder: "/home/docker_dev/oxford_spires_dataset/data/2024-03-13-observatory-quarter-01"
+  gt_folder: "/home/docker_dev/oxford_spires_dataset/data/ground_truth_cloud/observatory-quarter"
+
+  run_gt_cloud_processing: True
+  run_lidar_cloud_processing: True
+  run_lidar_cloud_evaluation: True
+  run_colmap: True
+  run_colmap_sim3: True
+  run_mvs: True
+  run_mvs_evaluation: True
+  run_nerfstudio: True
+  run_novel_view_synthesis_only: True
\ No newline at end of file
diff --git a/oxford_spires_utils/eval.py b/oxford_spires_utils/eval.py
index d3c8806..86bdab2 100644
--- a/oxford_spires_utils/eval.py
+++ b/oxford_spires_utils/eval.py
@@ -1,3 +1,4 @@
+import logging
 from pathlib import Path
 
 import matplotlib.pyplot as plt
@@ -6,10 +7,12 @@
 from matplotlib.colors import LinearSegmentedColormap
 from scipy.spatial import cKDTree as KDTree
 
+logger = logging.getLogger(__name__)
+
 
 def compute_p2p_distance(query_cloud: np.ndarray, reference_cloud: np.ndarray):
     ref_kd_tree = KDTree(reference_cloud)
-    distances, _ = ref_kd_tree.query(query_cloud)
+    distances, _ = ref_kd_tree.query(query_cloud, workers=-1)
     return distances
 
 
@@ -21,26 +24,63 @@ def get_recon_metrics(
 ):
     assert isinstance(input_cloud, np.ndarray) and isinstance(gt_cloud, np.ndarray)
     assert input_cloud.shape[1] == 3 and gt_cloud.shape[1] == 3
-    print("Computing Accuracy and Precision ...")
+    logger.info(f"Computing Accuracy and Precision ({precision_threshold}) ...")
     distances = compute_p2p_distance(input_cloud, gt_cloud)
     accuracy = np.mean(distances)
     precision = np.sum(distances < precision_threshold) / len(distances)
 
-    print("Computing Completeness and Recall ...")
+    logger.info(f"Computing Completeness and Recall ({recall_threshold}) ...")
     distances = compute_p2p_distance(gt_cloud, input_cloud)
     completeness = np.mean(distances)
     recall = np.sum(distances < recall_threshold) / len(distances)
+    f1_score = 2 * (precision * recall) / (precision + recall)
 
-    print("Done!")
-    return {
+    results = {
         "accuracy": accuracy,
         "precision": precision,
         "completeness": completeness,
         "recall": recall,
+        "f1_score": f1_score,
     }
+    return results
 
 
-def save_error_cloud(input_cloud: np.ndarray, reference_cloud: np.ndarray, save_path, cmap="bgyr"):
+def get_recon_metrics_multi_thresholds(
+    input_cloud: np.ndarray, gt_cloud: np.ndarray, thresholds: list = [0.02, 0.05, 0.1], max_distance=2.0
+):
+    assert isinstance(input_cloud, np.ndarray) and isinstance(gt_cloud, np.ndarray)
+    assert input_cloud.shape[1] == 3 and gt_cloud.shape[1] == 3
+    results = []
+
+    logger.info("Computing Accuracy and Precision ...")
+    input_to_gt_dist = compute_p2p_distance(input_cloud, gt_cloud)
+    input_to_gt_dist = input_to_gt_dist[input_to_gt_dist <= max_distance]
+    accuracy = np.mean(input_to_gt_dist)
+
+    logger.info("Computing Completeness and Recall ...")
+    gt_to_input_dist = compute_p2p_distance(gt_cloud, input_cloud)
+    gt_to_input_dist = gt_to_input_dist[gt_to_input_dist <= max_distance]
+    completeness = np.mean(gt_to_input_dist)
+
+    logger.info(f"Accuracy: {accuracy:.4f}, Completeness: {completeness:.4f}")
+    results.append({"accuracy": accuracy, "completeness": completeness})
+    for threshold in thresholds:
+        precision = np.sum(input_to_gt_dist < threshold) / len(input_to_gt_dist)
+        recall = np.sum(gt_to_input_dist < threshold) / len(gt_to_input_dist)
+        f1_score = 2 * (precision * recall) / (precision + recall)
+        results.append(
+            {
+                "threshold": threshold,
+                "precision": precision,
+                "recall": recall,
+                "f1_score": f1_score,
+            }
+        )
+        logger.info(f"threshold {threshold} m, precision: {precision:.4f}, recall: {recall:.4f}, f1: {f1_score:.4f}")
+    return results
+
+
+def save_error_cloud(input_cloud: np.ndarray, reference_cloud: np.ndarray, save_path, cmap="bgyr", max_distance=2.0):
     def get_BGYR_colourmap():
         colours = [
             (0, 0, 255),  # Blue
@@ -56,6 +96,8 @@ def get_BGYR_colourmap():
         return cmap
 
     distances = compute_p2p_distance(input_cloud, reference_cloud)
+    input_cloud = input_cloud[distances <= max_distance]
+    distances = distances[distances <= max_distance]
     distances = np.clip(distances, 0, 1)
     if cmap == "bgyr":
         cmap = get_BGYR_colourmap()
@@ -67,7 +109,7 @@ def get_BGYR_colourmap():
     test_cloud.points = o3d.utility.Vector3dVector(input_cloud)
     test_cloud.colors = o3d.utility.Vector3dVector(distances_cmap[:, :3])
     o3d.io.write_point_cloud(save_path, test_cloud)
-    print(f"Error cloud saved to {save_path}")
+    logger.info(f"diff cloud saved to {save_path}")
 
 
 if __name__ == "__main__":
diff --git a/oxford_spires_utils/point_cloud.py b/oxford_spires_utils/point_cloud.py
index e194f25..59782b1 100644
--- a/oxford_spires_utils/point_cloud.py
+++ b/oxford_spires_utils/point_cloud.py
@@ -1,3 +1,4 @@
+import logging
 from pathlib import Path
 
 import numpy as np
@@ -9,6 +10,8 @@
 
 from oxford_spires_utils.se3 import is_se3_matrix, xyz_quat_xyzw_to_se3_matrix
 
+logger = logging.getLogger(__name__)
+
 
 def transform_3d_cloud(cloud_np, transform_matrix):
     """Apply a transformation to the point cloud."""
@@ -174,3 +177,13 @@ def convert_e57_to_pcd(e57_file_path, pcd_file_path, check_output=True, pcd_lib=
         if has_colour:
             colours_np = np.array(saved_cloud.colors)
             assert np.allclose(colours_np, colours / 255, rtol=1e-5, atol=1e-8)
+
+
+def transform_cloud_with_se3(cloud_file, se3_matrix, output_cloud_file):
+    assert str(cloud_file).endswith(".pcd") or str(cloud_file).endswith(".ply")
+    assert is_se3_matrix(se3_matrix)[0], is_se3_matrix(se3_matrix)[1]
+    assert str(output_cloud_file).endswith(".pcd") or str(output_cloud_file).endswith(".ply")
+    cloud = o3d.io.read_point_cloud(str(cloud_file))
+    cloud.transform(se3_matrix)
+    o3d.io.write_point_cloud(str(output_cloud_file), cloud)
+    logger.info(f"Transformed point cloud with SE(3) and saved as {output_cloud_file}")
diff --git a/oxford_spires_utils/trajectory/file_interfaces/nerf.py b/oxford_spires_utils/trajectory/file_interfaces/nerf.py
index ae4d3c1..dd0f360 100644
--- a/oxford_spires_utils/trajectory/file_interfaces/nerf.py
+++ b/oxford_spires_utils/trajectory/file_interfaces/nerf.py
@@ -30,7 +30,7 @@ def __init__(self, file_path, nerf_reader_valid_folder_path="", nerf_reader_sort
         self.valid_folder_path = nerf_reader_valid_folder_path
         self.sort_timestamp = nerf_reader_sort_timestamp
 
-    def read_file(self):
+    def read_file(self, has_timestamp=True):
         """
         Read NeRF trajectory file (transforms.json)
         @return: PosePath3D from evo
@@ -44,21 +44,24 @@ def read_file(self):
             if self.valid_folder_path != "":
                 if not frame["file_path"].startswith(self.valid_folder_path):
                     continue
-            t_float128 = NeRFTrajUtils.get_t_float128_from_fname(frame["file_path"])
             T = np.array(frame["transform_matrix"])
             assert T.shape == (4, 4)
             assert np.allclose(T[3, :], np.array([0, 0, 0, 1]))
-            timestamps.append(t_float128)
             poses_se3.append(T)
+            if has_timestamp:
+                t_float128 = NeRFTrajUtils.get_t_float128_from_fname(frame["file_path"])
+                timestamps.append(t_float128)
 
-        timestamps = np.array(timestamps)
         poses_se3 = np.array(poses_se3)
-        if self.sort_timestamp:
-            sort_idx = np.argsort(timestamps)
-            timestamps = timestamps[sort_idx]
-            poses_se3 = poses_se3[sort_idx]
-
-        return evo.core.trajectory.PoseTrajectory3D(poses_se3=poses_se3, timestamps=timestamps)
+        if has_timestamp:
+            timestamps = np.array(timestamps)
+            if self.sort_timestamp:
+                sort_idx = np.argsort(timestamps)
+                timestamps = timestamps[sort_idx]
+                poses_se3 = poses_se3[sort_idx]
+            return evo.core.trajectory.PoseTrajectory3D(poses_se3=poses_se3, timestamps=timestamps)
+
+        return evo.core.trajectory.PosePath3D(poses_se3=poses_se3)
 
 
 class NeRFTrajWriter(BasicTrajWriter):
diff --git a/oxford_spires_utils/trajectory/nerf_json_handler.py b/oxford_spires_utils/trajectory/nerf_json_handler.py
new file mode 100644
index 0000000..dac6e5a
--- /dev/null
+++ b/oxford_spires_utils/trajectory/nerf_json_handler.py
@@ -0,0 +1,212 @@
+import json
+import shutil
+from pathlib import Path
+
+
+class NeRFJsonHandler:
+    # remember to create a copy of the json when removing,
+    # otherwise frames will be skipped
+    def __init__(self, input_json_path) -> None:
+        self.load_json(input_json_path)
+
+    def load_json(self, json_path):
+        with open(json_path, "r") as f:
+            self.traj = json.load(f)
+
+    def save_json(self, json_path):
+        with open(json_path, "w") as f:
+            json.dump(self.traj, f, indent=4)
+
+    def get_n_frames(self):
+        return len(self.traj["frames"])
+
+    def sort_frames(self):
+        self.traj["frames"].sort(key=lambda x: x["file_path"])
+
+    def sync_with_folder(self, folder_path, valid_ext=".jpg"):
+        # get all files in subfolder
+        ref_files = list(Path(folder_path).glob("**/*" + valid_ext))
+        print(f"{len(self.traj['frames'])} files in json")
+
+        count = 0
+        frames_copy = self.traj["frames"].copy()
+        for frame in self.traj["frames"]:
+            file_path = frame["file_path"]
+            exist = [file_path for ref_file in ref_files if file_path == ref_file.__str__()[-len(file_path) :]]
+            if len(exist) == 0:
+                # print(f"file {file_path} not exist, remove it from json")
+                frames_copy.remove(frame)
+                count += 1
+        print(f"{len(ref_files)} files in reference folder")
+        print(f"removed {count} files, {len(frames_copy)} left")
+        self.traj["frames"] = frames_copy
+
+    def remove_folder(self, folder_path):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            file_path = Path(frame["file_path"])
+            if file_path.parent == Path(folder_path):
+                self.traj["frames"].remove(frame)
+                print(f"removed {file_path} from json")
+        print(f"filter_folder {len(frames_copy)} files in json, {len(self.traj['frames'])} left")
+
+    def keep_folder(self, folder_path):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            file_path = Path(frame["file_path"])
+            if file_path.parent != Path(folder_path):
+                self.traj["frames"].remove(frame)
+                print(f"removed {file_path} from json")
+        print(f"filter_folder {len(frames_copy)} files in json, {len(self.traj['frames'])} left")
+
+    def rename_filename(self, old_folder=None, new_folder=None, prefix="", suffix="", base_folder=None):
+        for frame in self.traj["frames"]:
+            file_path = Path(frame["file_path"])
+            if old_folder is not None and new_folder is not None:
+                assert str(file_path).startswith(old_folder), f"{file_path} does not start with {old_folder}"
+                new_file_path = Path(str(file_path).replace(old_folder, new_folder))
+            new_file_path = str(new_file_path.parent / (prefix + new_file_path.stem + suffix + new_file_path.suffix))
+            frame["file_path"] = new_file_path
+            if base_folder is not None:
+                abs_old_file = Path(base_folder) / file_path
+                assert abs_old_file.exists(), f"{abs_old_file} not exist"
+                abs_new_file = Path(base_folder) / new_file_path
+                abs_new_file.parent.mkdir(parents=True, exist_ok=True)
+                shutil.copy(abs_old_file, abs_new_file)
+
+    def keep_timestamp_only(self, start_time, end_time):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            file_path = Path(frame["file_path"])
+            timestamp = float(file_path.stem)
+            if timestamp < start_time or timestamp > end_time:
+                self.traj["frames"].remove(frame)
+                print(f"removed {file_path} from json")
+        print(f"keep_timestamp_only {len(frames_copy)} files in json, {len(self.traj['frames'])} left")
+
+    def remove_timestamp(self, start_time, end_time):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            file_path = Path(frame["file_path"])
+            timestamp = float(file_path.stem)
+            if timestamp >= start_time and timestamp <= end_time:
+                self.traj["frames"].remove(frame)
+                print(f"removed {file_path} from json")
+        print(f"remove_timestamp {len(frames_copy)} files in json, {len(self.traj['frames'])} left")
+
+    def skip_frames(self, skip):
+        frames_copy = self.traj["frames"].copy()
+        # sort
+        # frames_copy.sort(key=lambda x: float(Path(x["file_path"]).stem))
+
+        self.traj["frames"] = frames_copy[::skip]
+        print(f"Skipping: {len(frames_copy)} files in json, {len(self.traj['frames'])} left")
+
+    def remove_intrinsics(self):
+        # frames_copy = self.traj["frames"].copy()
+        for frame in self.traj["frames"]:
+            frame.pop("fl_x")
+            frame.pop("fl_y")
+            frame.pop("cx")
+            frame.pop("cy")
+            frame.pop("k1")
+            frame.pop("k2")
+            frame.pop("k3")
+            frame.pop("k4")
+            frame.pop("h")
+            frame.pop("w")
+
+    def add_depth(self, depth_folder=None):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            if depth_folder is None:
+                # simply add a path to the depth file modified from the image file path
+                depth_file_path = frame["file_path"].replace("images", "depth").replace(".jpg", ".png")
+                frame["depth_file_path"] = depth_file_path
+            else:
+                # add if it exists, otherwise remove
+                depth_folder_stem = Path(depth_folder).stem
+                depth_file_path = frame["file_path"].replace("images", depth_folder_stem).replace(".jpg", ".png")
+                depth_file_path_full = Path(depth_folder).parent / depth_file_path
+                if depth_file_path_full.exists():
+                    frame["depth_file_path"] = depth_file_path.__str__()
+                else:
+                    print(f"{depth_file_path_full} not exist")
+                    self.traj["frames"].remove(frame)
+
+    def add_normal(self, normal_folder=None):
+        frames_copy = self.traj["frames"].copy()
+        for frame in frames_copy:
+            if normal_folder is None:
+                # simply add a path to the depth file modified from the image file path
+                normal_file_path = frame["file_path"].replace("images", "normal").replace(".jpg", ".png")
+                frame["normal_file_path"] = normal_file_path
+            else:
+                # only add if it exists, otherwise remove
+                normal_folder_stem = Path(normal_folder).stem
+                normal_file_path = frame["file_path"].replace("images", normal_folder_stem).replace(".jpg", ".png")
+                normal_file_path_full = Path(normal_folder).parent / normal_file_path
+                if normal_file_path_full.exists():
+                    frame["normal_file_path"] = normal_file_path.__str__()
+                else:
+                    print(f"{normal_file_path_full} not exist")
+                    self.traj["frames"].remove(frame)
+
+    def add_mask(self):
+        for frame in self.traj["frames"]:
+            frame["mask_path"] = frame["file_path"].replace("images", "masks")
+
+    def get_clouds_in_json(self, cloud_dir, output_dir):
+        cloud_dir = Path(cloud_dir)
+        output_dir = Path(output_dir)
+        if output_dir.exists():
+            shutil.rmtree(output_dir)
+        output_dir.mkdir(parents=True, exist_ok=True)
+        for frame in self.traj["frames"]:
+            # find / in file_path
+            img_path = Path(frame["file_path"])
+
+            cloud_path = cloud_dir / img_path.parent.name / img_path.name.replace(".jpg", ".pcd")
+            output_path = output_dir / img_path.parent.name / img_path.name.replace(".jpg", ".pcd")
+            # remove if exist
+            output_path.parent.mkdir(parents=True, exist_ok=True)
+            if cloud_path.exists():
+                shutil.copy(cloud_path, output_dir / cloud_path.parent.name / cloud_path.name)
+            else:
+                test_path = cloud_path.parent / (cloud_path.stem[:-1] + cloud_path.suffix)
+                if test_path.exists():
+                    shutil.copy(test_path, output_dir / cloud_path.parent.name / test_path.name)
+                else:
+                    print(f"{cloud_path} not exist")
+
+    def get_images_in_json(self, image_dir, output_dir):
+        image_dir = Path(image_dir)
+        output_dir = Path(output_dir)
+        if output_dir.exists():
+            shutil.rmtree(output_dir)
+        output_dir.mkdir(parents=True, exist_ok=True)
+        for frame in self.traj["frames"]:
+            # find / in file_path
+            img_path = image_dir / frame["file_path"]
+            assert img_path.exists(), f"{img_path} not exist"
+
+            output_path = output_dir / img_path.parent.name / img_path.name
+            output_path.parent.mkdir(parents=True, exist_ok=True)
+            if img_path.exists():
+                shutil.copy(img_path, output_dir / img_path.parent.name / img_path.name)
+            else:
+                print(f"{img_path} not exist")
+
+    def update_hw(self):
+        for frame in self.traj["frames"]:
+            frame["h"] = 935
+            # frame["w"] = w
+
+    def write_pose_cloud(self, output_file):
+        import open3d as o3d
+
+        output_cloud = o3d.geometry.PointCloud()
+        for frame in self.traj["frames"]:
+            xyz = get_xyz(frame)
+            output_cloud.points.append(xyz)
+        o3d.io.write_point_cloud(str(output_file), output_cloud)
diff --git a/requirements.txt b/requirements.txt
index d5488ee..1eb008f 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,7 +4,7 @@ scipy>=1.10.1
 pytest>=8.0.0
 pypcd4>=1.1.0
 pye57>=0.4.13
-nerfstudio==1.1.4
+# nerfstudio==1.1.4
 evo>=1.29.0
 pytransform3d>=3.5.0
 huggingface_hub>=0.25.1
\ No newline at end of file
diff --git a/scripts/dataset_download.py b/scripts/dataset_download.py
index 3cf1491..c8b86b2 100644
--- a/scripts/dataset_download.py
+++ b/scripts/dataset_download.py
@@ -9,10 +9,20 @@
 repo_id = "ori-drs/oxford_spires_dataset"
 download_sequences = True
 download_ground_truth = True
-dataset_sequences = ["2024-03-13-observatory-quarter-01", "2024-03-14-blenheim-05"]
+dataset_sequences = [
+    "2024-03-12-keble-college-04",
+    "2024-03-13-observatory-quarter-01",
+    "2024-03-14-blenheim-palace-05",
+    "2024-03-18-christ-church-02",
+]
 file_lists = ["images.zip", "lidar_slam.zip", "T_gt_lidar.txt"]
 
-ground_truth_lists = ["observatory-quarter", "blenheim-palace"]
+ground_truth_lists = [
+    "blenheim-palace",
+    "christ-church",
+    "keble-college",
+    "observatory-quarter",
+]
 ground_truth_lists = [f"ground_truth_cloud/{site}" for site in ground_truth_lists]
 cloud_file = "individual_cloud_e57.zip"
 
diff --git a/scripts/reconstruction_benchmark/lidar_cloud_eval.py b/scripts/reconstruction_benchmark/lidar_cloud_eval.py
deleted file mode 100644
index dc023da..0000000
--- a/scripts/reconstruction_benchmark/lidar_cloud_eval.py
+++ /dev/null
@@ -1,34 +0,0 @@
-from pathlib import Path
-
-import numpy as np
-import open3d as o3d
-
-from oxford_spires_utils.eval import get_recon_metrics, save_error_cloud
-from oxford_spires_utils.point_cloud import merge_downsample_vilens_slam_clouds
-from spires_cpp import convertOctreeToPointCloud, processPCDFolder, removeUnknownPoints
-
-
-def evaluate_lidar_cloud(
-    project_folder,
-    lidar_cloud_folder_path,
-    gt_octree_path,
-    gt_cloud_path,
-    octomap_resolution=0.1,
-    downsample_voxel_size=0.05,
-):
-    input_cloud_bt_path = Path(project_folder) / "input_cloud.bt"
-    processPCDFolder(str(lidar_cloud_folder_path), octomap_resolution, str(input_cloud_bt_path))
-
-    input_cloud_free_path = str(Path(input_cloud_bt_path).with_name(f"{Path(input_cloud_bt_path).stem}_free.pcd"))
-    input_cloud_occ_path = str(Path(input_cloud_bt_path).with_name(f"{Path(input_cloud_bt_path).stem}_occ.pcd"))
-    convertOctreeToPointCloud(str(input_cloud_bt_path), str(input_cloud_free_path), str(input_cloud_occ_path))
-
-    input_cloud_merged_path = Path(project_folder) / "input_cloud_merged.pcd"
-    _ = merge_downsample_vilens_slam_clouds(lidar_cloud_folder_path, downsample_voxel_size, input_cloud_merged_path)
-    input_cloud_filtered_path = Path(project_folder) / "input_cloud_merged_filtered.pcd"
-    removeUnknownPoints(str(input_cloud_merged_path), str(gt_octree_path), str(input_cloud_filtered_path))
-    input_cloud_np = np.asarray(o3d.io.read_point_cloud(str(input_cloud_filtered_path)).points)
-    gt_cloud_np = np.asarray(o3d.io.read_point_cloud(str(gt_cloud_path)).points)
-
-    print(get_recon_metrics(input_cloud_np, gt_cloud_np))
-    save_error_cloud(input_cloud_np, gt_cloud_np, str(Path(project_folder) / "input_error.pcd"))
diff --git a/scripts/reconstruction_benchmark/main.py b/scripts/reconstruction_benchmark/main.py
index 6322be3..dfb6258 100644
--- a/scripts/reconstruction_benchmark/main.py
+++ b/scripts/reconstruction_benchmark/main.py
@@ -1,23 +1,28 @@
+import argparse
+import csv
 import logging
+import shutil
 from copy import deepcopy
 from datetime import datetime
 from pathlib import Path
 
 import numpy as np
+import open3d as o3d
 import yaml
-from lidar_cloud_eval import evaluate_lidar_cloud
 from mvs import rescale_openmvs_cloud, run_openmvs
 from nerf import create_nerfstudio_dir, generate_nerfstudio_config, run_nerfstudio
-from sfm import rescale_colmap_json, run_colmap
+from sfm import export_json, rescale_colmap_json, run_colmap
 
 from oxford_spires_utils.bash_command import print_with_colour
-from oxford_spires_utils.point_cloud import merge_downsample_vilens_slam_clouds
+from oxford_spires_utils.eval import get_recon_metrics_multi_thresholds, save_error_cloud
+from oxford_spires_utils.point_cloud import merge_downsample_vilens_slam_clouds, transform_cloud_with_se3
+from oxford_spires_utils.se3 import is_se3_matrix
 from oxford_spires_utils.sensor import Sensor
 from oxford_spires_utils.trajectory.align import align
 from oxford_spires_utils.trajectory.file_interfaces import NeRFTrajReader, VilensSlamTrajReader
 from oxford_spires_utils.trajectory.utils import pose_to_ply
 from oxford_spires_utils.utils import convert_e57_folder_to_pcd_folder, transform_pcd_folder
-from spires_cpp import convertOctreeToPointCloud, processPCDFolder
+from spires_cpp import convertOctreeToPointCloud, processPCDFolder, removeUnknownPoints
 
 logger = logging.getLogger(__name__)
 
@@ -37,87 +42,138 @@ def setup_logging():
 
 
 class ReconstructionBenchmark:
-    def __init__(self, project_folder, sensor):
+    def __init__(self, project_folder, gt_folder, sensor):
         self.project_folder = Path(project_folder)
         self.project_folder.mkdir(parents=True, exist_ok=True)
+        self.gt_folder = Path(gt_folder)
         self.sensor = sensor
         self.camera_for_alignment = "cam_front"
         self.image_folder = self.project_folder / "images"
-        self.gt_individual_folder = self.project_folder / "gt_clouds"
-        self.individual_clouds_folder = self.project_folder / "lidar_clouds"
+        self.individual_clouds_folder = self.project_folder / "lidar_slam" / "individual_clouds"
+        self.lidar_slam_traj_file = self.project_folder / "lidar_slam" / "slam_poses.csv"
         self.output_folder = self.project_folder / "outputs"
         self.lidar_output_folder = self.output_folder / "lidar"
         self.lidar_output_folder.mkdir(exist_ok=True, parents=True)
         self.colmap_output_folder = self.output_folder / "colmap"
         self.colmap_output_folder.mkdir(exist_ok=True, parents=True)
+        self.recon_benchmark_dir = self.output_folder / "recon_benchmark"
+        self.recon_benchmark_dir.mkdir(exist_ok=True, parents=True)
         # TODO: check lidar cloud folder has viewpoints and is pcd, check gt folder is pcd, check image folder is jpg/png
         self.octomap_resolution = 0.1
-        self.cloud_downsample_voxel_size = 0.05
-        self.gt_octree_path = self.output_folder / "gt_cloud.bt"
-        self.gt_cloud_merged_path = self.output_folder / "gt_cloud_merged.pcd"
+        self.cloud_downsample_voxel_size = 0.01
+        self.gt_octree_path = self.recon_benchmark_dir / "gt_cloud.bt"
+        self.gt_cloud_merged_path = self.recon_benchmark_dir / "gt_cloud_merged.pcd"
+        self.gt_cloud_individual_e57_folder = self.gt_folder / "individual_cloud_e57"
+        self.gt_cloud_individual_pcd_folder = self.gt_folder / "individual_cloud_pcd"
 
-        self.colmap_sparse_folder = self.colmap_output_folder / "sparse" / "0"
+        self.colmap_sparse_folder = self.colmap_output_folder / "sparse"
+        self.colmap_sparse_0_folder = self.colmap_sparse_folder / "0"
         self.openmvs_bin = "/usr/local/bin/OpenMVS"
         self.mvs_output_folder = self.output_folder / "mvs"
         self.mvs_output_folder.mkdir(exist_ok=True, parents=True)
-        self.mvs_max_image_size = 600
+        self.scaled_mvs_cloud_gt_frame_file = self.recon_benchmark_dir / "OpenMVS_dense_cloud_gt_frame.pcd"
+        self.colmap_undistort_max_image_size = -1
 
         self.ns_data_dir = self.output_folder / "nerfstudio" / self.project_folder.name
         self.metric_json_filename = "transforms_metric.json"
-        self.ns_model_dir = self.ns_data_dir / "trained_models"
         logger.info(f"Project folder: {self.project_folder}")
 
+        self.lidar_cloud_merged_path = self.recon_benchmark_dir / "lidar_cloud_merged.pcd"
+        self.lidar_occ_benchmark_file = self.recon_benchmark_dir / "lidar_occ.pcd"
+
     def process_gt_cloud(self):
-        print_with_colour("Creating Octree and merged cloud from ground truth clouds")
-        processPCDFolder(str(self.gt_individual_folder), self.octomap_resolution, str(self.gt_octree_path))
+        logger.info("Converting ground truth clouds from e57 to pcd")
+        convert_e57_folder_to_pcd_folder(self.gt_cloud_individual_e57_folder, self.gt_cloud_individual_pcd_folder)
+        logger.info("Creating Octree and merged cloud from ground truth clouds")
+        processPCDFolder(str(self.gt_cloud_individual_pcd_folder), self.octomap_resolution, str(self.gt_octree_path))
         gt_cloud_free_path = str(Path(self.gt_octree_path).with_name(f"{Path(self.gt_octree_path).stem}_free.pcd"))
         gt_cloud_occ_path = str(Path(self.gt_octree_path).with_name(f"{Path(self.gt_octree_path).stem}_occ.pcd"))
         convertOctreeToPointCloud(str(self.gt_octree_path), str(gt_cloud_free_path), str(gt_cloud_occ_path))
+        logger.info("Merging and downsampling ground truth clouds")
         _ = merge_downsample_vilens_slam_clouds(
-            self.gt_individual_folder, self.cloud_downsample_voxel_size, self.gt_cloud_merged_path
-        )
-
-    def evaluate_lidar_clouds(self):
-        evaluate_lidar_cloud(
-            self.lidar_output_folder,
-            self.individual_clouds_folder,
-            self.gt_octree_path,
-            self.gt_cloud_merged_path,
-            self.octomap_resolution,
+            self.gt_cloud_individual_pcd_folder, self.cloud_downsample_voxel_size, self.gt_cloud_merged_path
         )
 
-    def tranform_lidar_clouds(self, transform_matrix_path=None):
+    def load_lidar_gt_transform(self, transform_matrix_path=None):
         if transform_matrix_path is None:
             transform_matrix_path = self.project_folder / "T_gt_lidar.txt"
+        logger.info(f"Loading transform matrix from {transform_matrix_path}")
         assert transform_matrix_path.exists(), f"Transform matrix not found at {transform_matrix_path}"
-        transform_matrix = np.loadtxt(transform_matrix_path)
-        new_individual_clouds_folder = self.project_folder / "lidar_clouds_transformed"
-        transform_pcd_folder(self.individual_clouds_folder, new_individual_clouds_folder, transform_matrix)
+        self.transform_matrix = np.loadtxt(transform_matrix_path)
+        assert is_se3_matrix(self.transform_matrix)[0], is_se3_matrix(self.transform_matrix)[1]
+
+    def process_lidar_clouds(self):
+        logger.info("Transforming lidar clouds to the same frame as the ground truth clouds")
+        new_individual_clouds_folder = self.lidar_output_folder / "lidar_clouds_transformed"
+        transform_pcd_folder(self.individual_clouds_folder, new_individual_clouds_folder, self.transform_matrix)
         self.individual_clouds_folder = new_individual_clouds_folder
+        logger.info("Creating Octree from transformed lidar clouds")
+        lidar_cloud_octomap_file = self.lidar_output_folder / "lidar_cloud.bt"
+        processPCDFolder(str(self.individual_clouds_folder), self.octomap_resolution, str(lidar_cloud_octomap_file))
+        logger.info("Converting Octree to point cloud")
+        lidar_cloud_free_path = Path(lidar_cloud_octomap_file).with_name(
+            f"{Path(lidar_cloud_octomap_file).stem}_free.pcd"
+        )
+        lidar_cloud_occ_path = Path(lidar_cloud_octomap_file).with_name(
+            f"{Path(lidar_cloud_octomap_file).stem}_occ.pcd"
+        )
+        convertOctreeToPointCloud(str(lidar_cloud_octomap_file), str(lidar_cloud_free_path), str(lidar_cloud_occ_path))
+        shutil.copy(lidar_cloud_occ_path, self.lidar_occ_benchmark_file)
+        logger.info("Merging and downsampling lidar clouds")
+        _ = merge_downsample_vilens_slam_clouds(
+            self.individual_clouds_folder, self.cloud_downsample_voxel_size, self.lidar_cloud_merged_path
+        )
 
-    def run_colmap(self):
-        run_colmap(self.image_folder, self.colmap_output_folder)
+    def run_colmap(self, matcher="vocab_tree_matcher"):
+        camera_model = "OPENCV_FISHEYE"
+        run_colmap(
+            self.image_folder,
+            self.colmap_output_folder,
+            matcher=matcher,
+            camera_model=camera_model,
+            max_image_size=self.colmap_undistort_max_image_size,
+            loop_detection_period=5,
+        )
+        export_json(
+            self.colmap_sparse_0_folder,
+            json_file_name="transforms.json",
+            output_dir=self.colmap_output_folder,
+        )
+        export_json(
+            input_bin_dir=self.colmap_output_folder / "dense" / "sparse",
+            json_file_name="transforms.json",
+            output_dir=self.colmap_output_folder / "dense",
+            db_file=self.colmap_output_folder / "database.db",
+        )
         create_nerfstudio_dir(self.colmap_output_folder, self.ns_data_dir, self.image_folder)
+        create_nerfstudio_dir(
+            self.colmap_output_folder / "dense",
+            self.ns_data_dir.with_name(self.ns_data_dir.name + "_undistorted"),
+            self.ns_data_dir / "dense" / self.image_folder.name,
+        )
 
     def run_openmvs(self):
         # check if multiple sparse folders exist
         num_sparse_folders = len(list(self.colmap_sparse_folder.glob("*")))
         if num_sparse_folders > 1:
             print_with_colour(f"Multiple sparse folders found in {self.colmap_output_folder}. Using the first one.")
-        run_openmvs(
+        assert self.transform_matrix is not None, "Ground truth to lidar transform not loaded"
+        mvs_cloud_file = run_openmvs(
             self.image_folder,
             self.colmap_output_folder,
-            self.colmap_sparse_folder,
+            self.colmap_sparse_0_folder,
             self.mvs_output_folder,
-            self.mvs_max_image_size,
             self.openmvs_bin,
         )
+        scaled_mvs_cloud_file = self.mvs_output_folder / "OpenMVS_dense_cloud_metric.pcd"
+        rescale_openmvs_cloud(mvs_cloud_file, T_lidar_colmap, scaled_mvs_cloud_file)
+        transform_cloud_with_se3(scaled_mvs_cloud_file, self.transform_matrix, self.scaled_mvs_cloud_gt_frame_file)
 
-    def compute_sim3(self):
-        lidar_slam_traj_file = self.project_folder / "slam_poses_robotics.csv"
+    def compute_colmap_sim3(self):
         colmap_traj_file = self.colmap_output_folder / "transforms.json"
         rescaled_colmap_traj_file = self.colmap_output_folder / self.metric_json_filename  # TODO refactor
-        lidar_slam_traj = VilensSlamTrajReader(lidar_slam_traj_file).read_file()
+        lidar_slam_traj = VilensSlamTrajReader(self.lidar_slam_traj_file).read_file()
+        pose_to_ply(lidar_slam_traj, self.colmap_output_folder / "lidar_slam_traj.ply", [1.0, 0.0, 0.0])
         camera_alignment = self.sensor.get_camera(self.camera_for_alignment)
         valid_folder_path = "images/" + Sensor.convert_camera_topic_to_folder_name(camera_alignment.topic)
         logger.info(f'Loading only "{self.camera_for_alignment}" with directory "{valid_folder_path}" from json file')
@@ -131,40 +187,120 @@ def compute_sim3(self):
         lidar_slam_traj_cam_frame = deepcopy(lidar_slam_traj)
         lidar_slam_traj_cam_frame.transform(T_base_cam, right_mul=True)
         # T_lidar_colmap = align(lidar_slam_traj, colmap_traj_single_cam, self.colmap_output_folder)
-        T_lidar_colmap = align(lidar_slam_traj_cam_frame, colmap_traj_single_cam, self.colmap_output_folder)
-        rescale_colmap_json(colmap_traj_file, T_lidar_colmap, rescaled_colmap_traj_file)
-        mvs_cloud_file = self.mvs_output_folder / "scene_dense_nerf_world.ply"
-        scaled_mvs_cloud_file = self.mvs_output_folder / "scene_dense_nerf_world_scaled.ply"
-        rescale_openmvs_cloud(mvs_cloud_file, T_lidar_colmap, scaled_mvs_cloud_file)
+        self.T_lidar_colmap = align(lidar_slam_traj_cam_frame, colmap_traj_single_cam, self.colmap_output_folder)
+        rescale_colmap_json(colmap_traj_file, self.T_lidar_colmap, rescaled_colmap_traj_file)
         rescaled_colmap_traj = NeRFTrajReader(rescaled_colmap_traj_file).read_file()
         pose_to_ply(rescaled_colmap_traj, self.colmap_output_folder / "rescaled_colmap_traj.ply", [0.0, 1.0, 0.0])
-        pose_to_ply(lidar_slam_traj, self.colmap_output_folder / "lidar_slam_traj.ply", [1.0, 0.0, 0.0])
+
         ns_metric_json_file = self.ns_data_dir / self.metric_json_filename
-        if not ns_metric_json_file.exists():
-            ns_metric_json_file.symlink_to(rescaled_colmap_traj_file)  # TODO remove old ones?
+        ns_metric_json_file.unlink(missing_ok=True)
+        ns_metric_json_file.symlink_to(rescaled_colmap_traj_file)
+
+    def run_nerfstudio(
+        self,
+        method="nerfacto",
+        ns_data_dir=None,
+        json_filename="transforms_metric.json",
+        eval_mode="fraction",
+        export_cloud=True,
+    ):
+        ns_data_dir = self.ns_data_dir if ns_data_dir is None else Path(ns_data_dir)
+        ns_model_dir = ns_data_dir / "trained_models"
+        assert ns_data_dir.exists(), f"nerfstudio directory not found at {ns_data_dir}"
+        ns_config, ns_data_config = generate_nerfstudio_config(
+            method, ns_data_dir / json_filename, ns_model_dir, eval_mode=eval_mode
+        )
+        metric_cloud_file = run_nerfstudio(ns_config, ns_data_config, export_cloud)
+        if metric_cloud_file is not None:
+            metric_cloud_gt_frame = self.recon_benchmark_dir / (metric_cloud_file.stem + "_gt_frame.pcd")
+            transform_cloud_with_se3(metric_cloud_file, self.transform_matrix, metric_cloud_gt_frame)
+            self.evaluate_reconstruction(metric_cloud_gt_frame)
 
-    def run_nerfstudio(self, method="nerfacto", json_filename="transforms_metric.json"):
-        assert self.ns_data_dir.exists(), f"nerfstudio directory not found at {self.ns_data_dir}"
-        ns_config = generate_nerfstudio_config(method, self.ns_data_dir / json_filename, self.ns_model_dir)
-        run_nerfstudio(ns_config)
+    def evaluate_reconstruction(self, input_cloud_path, results_dir=None):
+        assert input_cloud_path.exists(), f"Input cloud not found at {input_cloud_path}"
+        assert Path(input_cloud_path).suffix == ".pcd", "Input cloud must be a pcd file"
+        assert self.gt_octree_path.exists(), f"Ground truth octree not found at {self.gt_octree_path}"
+        recon_thresholds = [0.03, 0.05, 0.1, 0.2]
+        results_dir = self.recon_benchmark_dir if results_dir is None else Path(results_dir)
+        input_cloud = o3d.io.read_point_cloud(str(input_cloud_path))
+        logger.info(f"Cropping input cloud using bounding box from {self.gt_cloud_merged_path}")
+        gt_bbox = o3d.io.read_point_cloud(str(self.gt_cloud_merged_path)).get_axis_aligned_bounding_box()
+        input_cloud = input_cloud.crop(gt_bbox)
+        logger.info(f"Downsampling filtered cloud to {self.cloud_downsample_voxel_size} m")
+        input_cloud_downsampled = input_cloud.voxel_down_sample(voxel_size=self.cloud_downsample_voxel_size)
+        input_cloud_downsampled_path = results_dir / (input_cloud_path.stem + "_downsampled.pcd")
+        o3d.io.write_point_cloud(str(input_cloud_downsampled_path), input_cloud_downsampled)
+
+        logger.info(f"Removing unknown points using {self.gt_octree_path}")
+        filtered_input_cloud_path = results_dir / f"{Path(input_cloud_downsampled_path).stem}_filtered.pcd"
+        removeUnknownPoints(str(input_cloud_downsampled_path), str(self.gt_octree_path), str(filtered_input_cloud_path))
+        input_cloud_np = np.asarray(o3d.io.read_point_cloud(str(filtered_input_cloud_path)).points)
+        gt_cloud_np = np.asarray(o3d.io.read_point_cloud(str(self.gt_cloud_merged_path)).points)
+        recon_metrics = get_recon_metrics_multi_thresholds(input_cloud_np, gt_cloud_np, thresholds=recon_thresholds)
+        error_csv_file = results_dir / f"{Path(input_cloud_path).stem}_metrics.csv"
+        with open(error_csv_file, mode="w", newline="") as f:
+            writer = csv.DictWriter(f, fieldnames=recon_metrics[1].keys())
+            writer.writeheader()
+            for metric in recon_metrics[1:]:
+                writer.writerow(metric)
+            writer = csv.DictWriter(f, fieldnames=recon_metrics[0].keys())
+            writer.writeheader()
+            writer.writerow(recon_metrics[0])
+
+        error_cloud_file = results_dir / f"{Path(input_cloud_path).stem}_error.pcd"
+        save_error_cloud(input_cloud_np, gt_cloud_np, str(error_cloud_file))
+
+
+def get_args():
+    parser = argparse.ArgumentParser(description="Reconstruction Benchmark")
+    default_recon_config_file = Path(__file__).parent.parent.parent / "config" / "recon_benchmark.yaml"
+    parser.add_argument(
+        "--config-file",
+        type=str,
+        default=str(default_recon_config_file),
+    )
+    return parser.parse_args()
 
 
 if __name__ == "__main__":
     setup_logging()
+    recon_config_file = get_args().config_file
     logger.info("Starting Reconstruction Benchmark")
+    with open(recon_config_file, "r") as f:
+        recon_config = yaml.safe_load(f)["reconstruction_benchmark"]
     with open(Path(__file__).parent.parent.parent / "config" / "sensor.yaml", "r") as f:
         sensor_config = yaml.safe_load(f)["sensor"]
         sensor = Sensor(**sensor_config)
-    gt_cloud_folder_e57_path = "/home/oxford_spires_dataset/data/2024-03-13-maths_1/gt_individual_e57"
-    gt_cloud_folder_pcd_path = "/home/oxford_spires_dataset/data/2024-03-13-maths_1/gt_clouds"
-    convert_e57_folder_to_pcd_folder(gt_cloud_folder_e57_path, gt_cloud_folder_pcd_path)
-    project_folder = "/home/oxford_spires_dataset/data/2024-03-13-observatory-quarter-01"
-    recon_benchmark = ReconstructionBenchmark(project_folder, sensor)
-    recon_benchmark.process_gt_cloud()
-    recon_benchmark.tranform_lidar_clouds()
-    recon_benchmark.evaluate_lidar_clouds()
-    recon_benchmark.run_colmap()
-    recon_benchmark.run_openmvs()
-    recon_benchmark.compute_sim3()
-    recon_benchmark.run_nerfstudio("nerfacto", json_filename="transforms_metric.json")
-    recon_benchmark.run_nerfstudio("splatfacto")
+
+    recon_benchmark = ReconstructionBenchmark(recon_config["project_folder"], recon_config["gt_folder"], sensor)
+    recon_benchmark.load_lidar_gt_transform()
+    if recon_config["run_gt_cloud_processing"]:
+        recon_benchmark.process_gt_cloud()
+    if recon_config["run_lidar_cloud_processing"]:
+        recon_benchmark.process_lidar_clouds()
+    if recon_config["run_lidar_cloud_evaluation"]:
+        recon_benchmark.evaluate_reconstruction(recon_benchmark.lidar_cloud_merged_path)
+    if recon_config["run_colmap"]:
+        recon_benchmark.run_colmap("sequential_matcher")
+    if recon_config["run_colmap_sim3"]:
+        recon_benchmark.compute_colmap_sim3()
+    if recon_config["run_mvs"]:
+        recon_benchmark.run_openmvs()
+    if recon_config["run_mvs_evaluation"]:
+        recon_benchmark.evaluate_reconstruction(recon_benchmark.scaled_mvs_cloud_gt_frame_file)
+    if recon_config["run_nerfstudio"]:
+        if recon_config["run_novel_view_synthesis_only"]:
+            # fmt: off
+            undistorted_ns_dir = recon_benchmark.ns_data_dir.with_name(recon_benchmark.ns_data_dir.name + "_undistorted") # noqa: F401
+            recon_benchmark.run_nerfstudio("nerfacto", json_filename="transforms_train.json", eval_mode="fraction", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("nerfacto", json_filename="transforms_train_eval.json", eval_mode="filename", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("nerfacto-big", json_filename="transforms_train.json", eval_mode="fraction", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("nerfacto-big", json_filename="transforms_train_eval.json", eval_mode="filename", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("splatfacto", json_filename="transforms_train.json", eval_mode="fraction", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("splatfacto", json_filename="transforms_train_eval.json", eval_mode="filename", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("splatfacto-big", json_filename="transforms_train.json", eval_mode="fraction", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            recon_benchmark.run_nerfstudio("splatfacto-big", json_filename="transforms_train_eval.json", eval_mode="filename", ns_data_dir=undistorted_ns_dir, export_cloud=False)# noqa
+            # fmt: on
+        else:
+            recon_benchmark.run_nerfstudio("nerfacto", json_filename="transforms_metric.json")
+            recon_benchmark.run_nerfstudio("splatfacto")
diff --git a/scripts/reconstruction_benchmark/mvs.py b/scripts/reconstruction_benchmark/mvs.py
index 3be9d09..9a3aade 100644
--- a/scripts/reconstruction_benchmark/mvs.py
+++ b/scripts/reconstruction_benchmark/mvs.py
@@ -52,24 +52,11 @@ def run_colmap_mvs(image_path, colmap_output_path, sparse_folder, max_image_size
     run_command(colmap_delauany_mesh_filter_cmd, print_command=True)
 
 
-def run_openmvs(
-    image_path, colmap_output_path, sparse_folder, mvs_dir, max_image_size, openmvs_bin="/usr/local/bin/OpenMVS"
-):
+def run_openmvs(image_path, colmap_output_path, sparse_folder, mvs_dir, openmvs_bin="/usr/local/bin/OpenMVS"):
     logger.info(f"Running OpenMVS; img_path {image_path}; output: {mvs_dir}")
     colmap_output_path = Path(colmap_output_path)
     mvs_dir.mkdir(parents=True, exist_ok=True)
 
-    colmap_image_undistorter_cmd = [
-        "colmap image_undistorter",
-        f"--image_path {image_path}",
-        f"--input_path {sparse_folder}",
-        f"--output_path {colmap_output_path/'dense'}",
-        "--output_type COLMAP",
-        f"--max_image_size {max_image_size}",
-    ]
-    colmap_image_undistorter_cmd = " ".join(colmap_image_undistorter_cmd)
-    run_command(colmap_image_undistorter_cmd, print_command=True)
-
     # Export to openMVS
     export_cmd = [
         f"{openmvs_bin}/InterfaceCOLMAP",
@@ -95,8 +82,10 @@ def run_openmvs(
         logger.error(f"Failed to generate dense point cloud at {output_ply_file}")
     dense_ply = o3d.io.read_point_cloud(str(output_ply_file))
     dense_ply.transform(colmap_to_nerf_world_transform)
-    o3d.io.write_point_cloud(str(output_ply_file.with_name("scene_dense_nerf_world.ply")), dense_ply)
+    output_file = output_ply_file.with_name("scene_dense_nerf_world.pcd")
+    o3d.io.write_point_cloud(str(output_file), dense_ply)
     logger.info("Transformed MVS point cloud to the world frame defined by the nerf convention")
+    return output_file
     # Reconstruct the mesh
     # reconstruct_cmd = [f"{openmvs_bin}/ReconstructMesh", "scene_dense.mvs", "-p scene_dense.ply", f"-w {mvs_dir}"]
     # run_command(" ".join(reconstruct_cmd), print_command=True)
diff --git a/scripts/reconstruction_benchmark/nerf.py b/scripts/reconstruction_benchmark/nerf.py
index 505750e..d8a68c2 100644
--- a/scripts/reconstruction_benchmark/nerf.py
+++ b/scripts/reconstruction_benchmark/nerf.py
@@ -1,53 +1,67 @@
 import json
+import logging
+import os
 import sys
 from pathlib import Path
 
 import numpy as np
 import open3d as o3d
+from nerfstudio.scripts.eval import entrypoint as eval_entrypoint
 from nerfstudio.scripts.exporter import entrypoint as exporter_entrypoint
 from nerfstudio.scripts.train import entrypoint as train_entrypoint
 
 from oxford_spires_utils.bash_command import print_with_colour
 
+logger = logging.getLogger(__name__)
+
 
 def generate_nerfstudio_config(
-    method, data_dir, output_dir, iterations=30000, eval_step=500, vis="wandb", cam_opt_mode="off"
+    method, data_dir, output_dir, iterations=5000, eval_step=500, vis="wandb", cam_opt_mode="off", eval_mode="fraction"
 ):
+    exp_name = Path(data_dir).name if Path(data_dir).is_dir() else Path(data_dir).parent.name
     ns_config = {
         "method": method,
-        "data": str(data_dir),
+        "experiment-name": str(exp_name),
         "output-dir": str(output_dir),
         "vis": vis,
         "max-num-iterations": iterations,
         "pipeline.model.camera-optimizer.mode": cam_opt_mode,
         "steps-per-eval-image": eval_step,
     }
-    return ns_config
+    ns_data_config = {
+        "dataparser": "nerfstudio-data",
+        "data": str(data_dir),
+        "eval-mode": eval_mode,
+    }
+    return ns_config, ns_data_config
 
 
 def create_nerfstudio_dir(colmap_dir, ns_dir, image_dir):
-    ns_dir = Path(ns_dir)
-    colmap_dir = Path(colmap_dir)
-    image_dir = Path(image_dir)
+    ns_dir = Path(ns_dir).resolve()
+    colmap_dir = Path(colmap_dir).resolve()
+    image_dir = Path(image_dir).resolve()
     # Ensure ns_dir exists
     ns_dir.mkdir(parents=True, exist_ok=True)
 
     # Symlink image_dir to ns_dir
     image_symlink = ns_dir / image_dir.name
     if not image_symlink.exists():
-        image_symlink.symlink_to(image_dir)
+        relative_image_dir = Path(os.path.relpath(str(image_dir.parent), str(ns_dir))) / image_dir.name
+        image_symlink.symlink_to(relative_image_dir)
 
     # Symlink contents of colmap_dir to ns_dir
     for item in colmap_dir.iterdir():
         item_symlink = ns_dir / item.name
         if not item_symlink.exists():
-            item_symlink.symlink_to(item)
+            relative_item = Path(os.path.relpath(str(colmap_dir), str(ns_dir))) / item.name
+            item_symlink.symlink_to(relative_item)
 
 
-def update_argv(nerfstudio_config):
-    assert sys.argv[0].endswith(".py") and len(sys.argv) == 1, "No args should be provided for the script"
+def update_argv(nerfstudio_config, follow_up=False):
+    if not follow_up:
+        assert sys.argv[0].endswith(".py") and len(sys.argv) == 1, "No args should be provided for the script"
     for k, v in nerfstudio_config.items():
-        if k == "method":
+        if k in ("method", "dataparser"):
             sys.argv.append(f"{v}")
         else:
             sys.argv.append(f"--{k}")
@@ -55,23 +69,51 @@ def update_argv(nerfstudio_config):
     print_with_colour(" ".join(sys.argv))
 
 
-def run_nerfstudio(ns_config):
+def run_nerfstudio(ns_config, ns_data_config, export_cloud=True):
+    logger.info(f"Running '{ns_config['method']}' on {ns_data_config['data']}")
+    logging.disable(logging.DEBUG)
     update_argv(ns_config)
+    update_argv(ns_data_config, follow_up=True)
     train_entrypoint()
     sys.argv = [sys.argv[0]]
-    ns_data = Path(ns_config["data"])
-    folder_name = ns_data.name if ns_data.is_dir() else ns_data.parent.name
-    output_log_dir = Path(ns_config["output-dir"]) / folder_name / ns_config["method"]
-    lastest_output_folder = sorted([x for x in output_log_dir.glob("*") if x.is_dir()])[-1]
-    latest_output_config = lastest_output_folder / "config.yml"
+    folder_name = ns_config["experiment-name"]
+    # rename nerfacto-big or nerfacto-huge to nerfacto, splatfacto-big to splatfacto
+    method_dir_name = ns_config["method"].replace("-big", "").replace("-huge", "")
+    output_log_dir = Path(ns_config["output-dir"]) / folder_name / method_dir_name
+    latest_output_folder = sorted([x for x in output_log_dir.glob("*") if x.is_dir()])[-1]
+    latest_output_config = latest_output_folder / "config.yml"
+
+    # evaluate renders
+    logger.info(f"Evaluating {latest_output_folder.name} from {latest_output_folder}")
+    render_dir = latest_output_folder / "renders"
+    run_nerfstudio_eval(latest_output_config, render_dir)
+    logging.disable(logging.NOTSET)
+
+    # export cloud
+    if not export_cloud:
+        return None
     export_method = "gaussian-splat" if ns_config["method"] == "splatfacto" else "pointcloud"
     output_cloud_file = run_nerfstudio_exporter(latest_output_config, export_method)
-    ns_se3, scale_matrix = load_ns_transform(lastest_output_folder)
+    ns_se3, scale_matrix = load_ns_transform(latest_output_folder)
     cloud = o3d.io.read_point_cloud(str(output_cloud_file))
-
     cloud.transform(scale_matrix)
     cloud.transform(np.linalg.inv(ns_se3))
-    o3d.io.write_point_cloud(str(output_cloud_file.with_name("input_scale.ply")), cloud)
+    final_metric_cloud_file = output_cloud_file.with_name(f'{ns_config["method"]}_cloud_metric.pcd')
+    o3d.io.write_point_cloud(str(final_metric_cloud_file), cloud)
+    return final_metric_cloud_file
+
+
+def run_nerfstudio_eval(config_file, render_dir):
+    output_eval_file = config_file.parent / "eval_results.json"
+    eval_config = {
+        "load-config": config_file,
+        "output-path": output_eval_file,
+        "render-output-path": render_dir,
+    }
+    update_argv(eval_config)
+    eval_entrypoint()
+    logger.info(f"Nerfstudio eval results\n{json.load(output_eval_file.open())}")
+    sys.argv = [sys.argv[0]]
 
 
 def run_nerfstudio_exporter(config_file, export_method):
@@ -79,6 +121,7 @@ def run_nerfstudio_exporter(config_file, export_method):
         "method": export_method,
         "load-config": config_file,
         "output-dir": config_file.parent,
+        "num-points": 100000000,
     }
     if export_method == "pointcloud":
         exporter_config["normal-method"] = "open3d"
diff --git a/scripts/reconstruction_benchmark/nerf_json.py b/scripts/reconstruction_benchmark/nerf_json.py
new file mode 100644
index 0000000..d552f8b
--- /dev/null
+++ b/scripts/reconstruction_benchmark/nerf_json.py
@@ -0,0 +1,103 @@
+from copy import deepcopy
+from pathlib import Path
+
+from nerf import create_nerfstudio_dir
+
+from oxford_spires_utils.trajectory.file_interfaces.nerf import NeRFTrajReader
+from oxford_spires_utils.trajectory.nerf_json_handler import NeRFJsonHandler
+from oxford_spires_utils.trajectory.utils import pose_to_ply
+
+
+def select_json_with_time_range(json_file, start_time, end_time, save_path):
+    nerf_json_handler = NeRFJsonHandler(json_file)
+    nerf_json_handler.sort_frames()
+    nerf_json_handler.keep_timestamp_only(start_time, end_time)
+    nerf_json_handler.save_json(save_path)
+
+    nerf_traj = NeRFTrajReader(save_path)
+    nerf_pose = nerf_traj.read_file()
+    pose_ply_file = save_path.with_suffix(".ply")
+    pose_to_ply(nerf_pose, pose_ply_file)
+
+
+def select_json_with_folder(json_file, img_folder, save_path):
+    nerf_json_handler = NeRFJsonHandler(json_file)
+    nerf_json_handler.sort_frames()
+    nerf_json_handler.sync_with_folder(img_folder)
+    nerf_json_handler.save_json(save_path)
+
+
+def split_json_every_n(json_file, n, save_path_removed, save_path_kept):
+    nerf_json_handler = NeRFJsonHandler(json_file)
+    nerf_json_handler.sort_frames()
+    train_json = deepcopy(nerf_json_handler)
+    removed_frames = nerf_json_handler.skip_frames(8, return_removed=True)
+    train_json.traj["frames"] = removed_frames
+    train_json.save_json(save_path_removed)
+    nerf_json_handler.save_json(save_path_kept)
+
+
+def merge_json_files(json_train_file, json_eval_file, image_dir, new_image_dir, merged_json_file):
+    json_train = NeRFJsonHandler(json_train_file)
+    json_eval = NeRFJsonHandler(json_eval_file)
+    json_train.rename_filename(
+        old_folder="images", new_folder=new_image_dir.stem, prefix="train_", base_folder=str(Path(image_dir).parent)
+    )
+    json_train.save_json(str(merged_json_file.with_name("transforms_train.json")))
+    json_eval.rename_filename(
+        old_folder="images", new_folder=new_image_dir.stem, prefix="eval_", base_folder=str(Path(image_dir).parent)
+    )
+    # merge
+    assert json_train.traj.keys() == json_eval.traj.keys()
+    new_json = deepcopy(json_train)
+    new_json.traj["frames"] += json_eval.traj["frames"]
+    new_json.save_json(str(merged_json_file))
+
+
+# dataset_folder = "/home/docker_dev/oxford_spires_dataset/data/roq-full"
+
+# train_name = "seq_1_fountain"
+# train_start_time = 1710338123.042936934
+# train_end_time = 1710338186.342030327
+
+# eval_name = "seq_1_fountain_back"
+# eval_start_time = 1710338353.039451086
+# eval_end_time = 1710338386.638942551
+
+dataset_folder = "/home/docker_dev/oxford_spires_dataset/data/2024-03-14-blenheim-palace-all"
+train_name = "seq_5"
+eval_name = "seq_1"
+
+use_undistorted_image = True
+
+dataset_folder = Path(dataset_folder).resolve()
+colmap_folder = dataset_folder / "outputs/colmap"
+if use_undistorted_image:
+    colmap_folder = colmap_folder / "dense"
+json_file = colmap_folder / "transforms.json"
+train_image_folder = dataset_folder / "train_val_image" / "train"
+eval_image_folder = dataset_folder / "train_val_image" / "eval"
+
+assert (train_image_folder / "images").exists(), f"{train_image_folder/'images'} does not exist"
+assert (eval_image_folder / "images").exists(), f"{eval_image_folder/'images'} does not exist"
+
+train_save_path = Path(json_file).parent / f"{train_name}.json"
+eval_save_path = Path(json_file).parent / f"{eval_name}.json"
+select_json_with_folder(json_file, train_image_folder, train_save_path)
+select_json_with_folder(json_file, eval_image_folder, eval_save_path)
+# select_json_with_time_range(json_file, train_start_time, train_end_time, train_save_path)
+# select_json_with_time_range(json_file, eval_start_time, eval_end_time, eval_save_path)
+
+
+# split_json_every_n(train_save_path, n=8, save_path_kept = train_save_path, save_path_removed
+image_dir = dataset_folder / "images" if not use_undistorted_image else colmap_folder / "images"
+json_train_file = colmap_folder / (train_name + ".json")
+json_eval_file = colmap_folder / (eval_name + ".json")
+new_image_dir = image_dir.parent / "images_train_eval"
+merged_json_file = colmap_folder / "transforms_train_eval.json"
+merge_json_files(json_train_file, json_eval_file, image_dir, new_image_dir, merged_json_file)
+# create json with the new train/eval prefix
+
+ns_dir = dataset_folder / "outputs" / "nerfstudio" / (dataset_folder.stem + "_undistorted")
+print(f"Creating NeRF Studio directory at {ns_dir}")
+create_nerfstudio_dir(colmap_folder, ns_dir, new_image_dir)
diff --git a/scripts/reconstruction_benchmark/sfm.py b/scripts/reconstruction_benchmark/sfm.py
index a84355b..beb0704 100644
--- a/scripts/reconstruction_benchmark/sfm.py
+++ b/scripts/reconstruction_benchmark/sfm.py
@@ -49,7 +49,15 @@ def get_vocab_tree(image_num) -> Path:
     return vocab_tree_filename
 
 
-def run_colmap(image_path, output_path, camera_model="OPENCV_FISHEYE"):
+def run_colmap(
+    image_path,
+    output_path,
+    camera_model="OPENCV_FISHEYE",
+    matcher="vocab_tree_matcher",
+    loop_detection_period=10,
+    sift_max_num_features=8192,
+    max_image_size=1000,
+):
     logger.debug(f"Running colmap; img_path {image_path}; output: {output_path}, {camera_model}")
     assert camera_model in camera_model_list, f"{camera_model} not supported. Supported models: {camera_model_list}"
     database_path = output_path / "database.db"
@@ -63,23 +71,30 @@ def run_colmap(image_path, output_path, camera_model="OPENCV_FISHEYE"):
         f"--database_path {database_path}",
         "--ImageReader.single_camera_per_folder 1",
         f"--ImageReader.camera_model {camera_model}",
+        f"--SiftExtraction.max_num_features {sift_max_num_features}",
     ]
     colmap_feature_extractor_cmd = " ".join(colmap_feature_extractor_cmd)
     logger.info(f"Running {colmap_feature_extractor_cmd}")
     run_command(colmap_feature_extractor_cmd, print_command=False)
-    colmap_db = COLMAPDatabase.connect(database_path)
-    total_image_num = colmap_db.execute("SELECT COUNT(*) FROM images").fetchone()[0]
-    logger.debug(f"Total number of images in COLMAP database: {total_image_num}")
 
     image_num = len(list(image_path.rglob("*")))
-    colmap_vocab_tree_matcher_cmd = [
-        "colmap vocab_tree_matcher",
+    colmap_matcher_cmd = [
+        f"colmap {matcher}",
         f"--database_path {database_path}",
-        f"--VocabTreeMatching.vocab_tree_path {get_vocab_tree(image_num)}",
     ]
-    colmap_vocab_tree_matcher_cmd = " ".join(colmap_vocab_tree_matcher_cmd)
-    logger.info(f"Running {colmap_vocab_tree_matcher_cmd}")
-    run_command(colmap_vocab_tree_matcher_cmd, print_command=False)
+    if matcher == "vocab_tree_matcher":
+        colmap_matcher_cmd.append(f"--VocabTreeMatching.vocab_tree_path {get_vocab_tree(image_num)}")
+    elif matcher == "sequential_matcher":
+        colmap_matcher_cmd.append("--SequentialMatching.loop_detection 1")
+        colmap_matcher_cmd.append(f"--SequentialMatching.vocab_tree_path {get_vocab_tree(image_num)}")
+        colmap_matcher_cmd.append(f"--SequentialMatching.loop_detection_period {loop_detection_period}")
+    else:
+        raise ValueError(
+            f"matcher {matcher} not supported. Supported matchers: ['vocab_tree_matcher', 'sequential_matcher']"
+        )
+    colmap_matcher_cmd = " ".join(colmap_matcher_cmd)
+    logger.info(f"Running {colmap_matcher_cmd}")
+    run_command(colmap_matcher_cmd, print_command=False)
 
     mapper_ba_global_function_tolerance = 1e-5
     colmap_mapper_cmd = [
@@ -103,15 +118,20 @@ def run_colmap(image_path, output_path, camera_model="OPENCV_FISHEYE"):
     logger.info(f"Running {colmap_ba_cmd}")
     run_command(colmap_ba_cmd, print_command=False)
 
+    colmap_image_undistorter_cmd = [
+        "colmap image_undistorter",
+        f"--image_path {image_path}",
+        f"--input_path {sparse_0_path}",
+        f"--output_path {output_path/'dense'}",
+        "--output_type COLMAP",
+        f"--max_image_size {max_image_size}",
+    ]
+    colmap_image_undistorter_cmd = " ".join(colmap_image_undistorter_cmd)
+    logger.info(f"Running {colmap_image_undistorter_cmd}")
+    run_command(colmap_image_undistorter_cmd, print_command=False)
+
     # from nerfstudio.process_data.colmap_utils import colmap_to_json
     # num_image_matched = colmap_to_json(recon_dir=sparse_0_path, output_dir=output_path)
-    logger.info("Exporting COLMAP to json file")
-    num_frame_matched = export_json(
-        sparse_0_path, json_file_name="transforms.json", output_dir=output_path, camera_model=camera_model
-    )
-    logger.info(
-        f"COLMAP matched {num_frame_matched} / {total_image_num} images {num_frame_matched / total_image_num * 100:.2f}%"
-    )
 
 
 def rescale_colmap_json(json_file, sim3_matrix, output_file):
@@ -138,11 +158,13 @@ def rescale_colmap_json(json_file, sim3_matrix, output_file):
         json.dump(data, f, indent=2)
 
 
-def export_json(input_bin_dir=None, json_file_name="transforms.json", output_dir=None, camera_model="OPENCV_FISHEYE"):
+def export_json(input_bin_dir, json_file_name="transforms.json", output_dir=None, db_file=None):
+    logger.info("Exporting COLMAP to json file")
     camera_mask_path = None
     input_bin_dir = Path(input_bin_dir)
     cameras_path = input_bin_dir / "cameras.bin"
     images_path = input_bin_dir / "images.bin"
+    database_path = output_dir / "database.db" if db_file is None else Path(db_file)
     output_dir = input_bin_dir if output_dir is None else Path(output_dir)
 
     cameras = read_cameras_binary(cameras_path)
@@ -158,7 +180,7 @@ def export_json(input_bin_dir=None, json_file_name="transforms.json", output_dir
         c2w = np.linalg.inv(w2c)  # this is the coordinate for openMVS
         c2w = get_nerf_pose(c2w)
 
-        frame = generate_json_camera_data(camera, camera_model)
+        frame = generate_json_camera_data(camera)
         frame["file_path"] = Path(f"./images/{im_data.name}").as_posix()  # assume images not in image path in colmap
         frame["transform_matrix"] = c2w.tolist()
         frame["colmap_img_id"] = img_id
@@ -166,20 +188,28 @@ def export_json(input_bin_dir=None, json_file_name="transforms.json", output_dir
             frame["mask_path"] = camera_mask_path.relative_to(camera_mask_path.parent.parent).as_posix()
 
         frames.append(frame)
+    frames = sorted(frames, key=lambda x: x["file_path"])
 
     out = {}
-    out["camera_model"] = camera_model
+    out["camera_model"] = camera.model
     out["frames"] = frames
+    num_frame_matched = len(frames)
+
+    colmap_db = COLMAPDatabase.connect(database_path)
+    total_image_num = colmap_db.execute("SELECT COUNT(*) FROM images").fetchone()[0]
+    logger.info(
+        f"COLMAP matched {num_frame_matched} / {total_image_num} images {num_frame_matched / total_image_num * 100:.2f}%"
+    )
 
     # Save for scale adjustment later
     assert json_file_name[-5:] == ".json"
     with open(output_dir / json_file_name, "w", encoding="utf-8") as f:
         json.dump(out, f, indent=4)
-    return len(frames)
 
 
-def generate_json_camera_data(camera, camera_model):
-    assert camera_model in ["OPENCV_FISHEYE", "OPENCV"]
+def generate_json_camera_data(camera):
+    camera_model = camera.model
+    assert camera_model in ["OPENCV_FISHEYE", "OPENCV", "PINHOLE"]
     data = {
         "fl_x": float(camera.params[0]),
         "fl_y": float(camera.params[1]),