Stereo reconstruction¶
This notebook introduce you to the procedure of building a 3D reconstruction from two stereo cameras. The procedure is as follows:
- Create a
Epochobject that stores all the relevant information for the 3D reconstruction. This includes:- Load the images from the two cameras taken at the same time
- Create the camera objects given the pre-calibrated intrinsic parameters
- Load information about the targets
- Detect and match features on the two images
- Perform a relative orientation between the two cameras based on the matched tie points (optionally, you can also set the scale of the reconstruction, e.g., by giving the camera baseline)
- Triangulate the tie points into the object space
- Perform an absolute orientation to refine the camera positions and the 3D points
- (optional) Perform a bundle adjustment to refine the camera positions and the 3D points by using Agisoft Metashape
Initialization¶
Let's first set up the python environment by importing the required libraries.
%load_ext autoreload
%autoreload 2
# Import required standard modules
from pathlib import Path
import numpy as np
# Import required icepy4d4D modules
from icepy4d import core as icecore
from icepy4d.core import Epoch
from icepy4d import matching
from icepy4d import sfm
from icepy4d import io
from icepy4d import utils
from icepy4d.metashape import metashape as MS
from icepy4d.utils import initialization
Jupyter environment detected. Enabling Open3D WebVisualizer. [Open3D INFO] WebRTC GUI backend enabled. [Open3D INFO] WebRTCWindowSystem: HTTP handshake server disabled.
First, you have to define the path to the configuration file (.yaml file).
This file contains all the paths and parameters needed to run the code.
See the config.yaml file in the nootebook folder for an example and refer to the documentation for how to prepare all the data for ICEpy4D.
Additionally, you can setup a logger for the code to print out some information and a timer to measure the runtime of the code.
# Parse the configuration file
CFG_FILE = "config.yaml"
# Parse the configuration file
cfg_file = Path(CFG_FILE)
cfg = initialization.parse_cfg(cfg_file, ignore_errors=True)
# Initialize the logger
logger = utils.get_logger()
# Initialize a timer to measure the processing time
timer = utils.AverageTimer()
# Get the list of cameras from the configuration file
cams = cfg.cams
================================================================ ICEpy4D Image-based Continuos monitoring of glaciers' Evolution with low-cost stereo-cameras and Deep Learning photogrammetry 2023 - Francesco Ioli - francesco.ioli@polimi.it ================================================================ 2023-10-03 10:41:41 | [INFO ] Configuration file: config.yaml 2023-10-03 10:41:41 | [INFO ] Epoch_to_process set to a pair of values. Expanding it for a range of epoches from epoch 0 to 158.
The Epoch class¶
An Epoch object is the main object that stores all the information about the 3D reconstruction. It is initialized by defining the timestamp of the epoch and, optionally, by storing the other information about the epoch (e.g., images, cameras, tie points etc...).
Let's first define a dictionary that contains the information about all the images that we are going to process. This dictionary has the camera name as key and an Image object as value. The Image object is an ICEpy4D object that manage the image in a lazy way (i.e., it doesn't load the image into memory, but read only image metadata). An image object is initialized by passing the path to the image to the constructor, and it automatically reads the image timestamp and other exif metadata.
image = Image('path_to_image')
To get the image timestamp, you can access the timestamp attribute of the image object.
image.timestamp
# Build a dictionary of images containing the name of the cameras as keys and Image objects as values
im_epoch = {
cams[0]: icecore.Image("../data/img/p1/IMG_2637.jpg"),
cams[1]: icecore.Image("../data/img/p2/IMG_1112.jpg")
}
# Get epoch timestamp as the timestamp of the first image and define epoch directory
epoch_timestamp = im_epoch[cams[0]].timestamp
epochdir = cfg.paths.results_dir / epoch_timestamp
We now have to load/define all the other information for the 3D reconstruction. In particular, this includes the camera objects (by loading their pre-calibrated intrinsics orientation) and the GCPs (by loading their coordinates in the object space and their projections in the image space).
# Load cameras
cams_ep = {}
for cam in cams:
calib = icecore.Calibration(cfg.paths.calibration_dir / f"{cam}.txt")
cams_ep[cam] = calib.to_camera()
# Load targets
target_paths = [
cfg.georef.target_dir
/ (im_epoch[cam].stem + cfg.georef.target_file_ext)
for cam in cams
]
targ_ep = icecore.Targets(
im_file_path=target_paths,
obj_file_path=cfg.georef.target_dir
/ cfg.georef.target_world_file,
)
# Create empty features
feat_ep = {cam: icecore.Features() for cam in cams}
Now we can create an Epoch object by passing all the information that we have previously defined:
# Create the epoch object
epoch = Epoch(
timestamp=epoch_timestamp,
images=im_epoch,
cameras=cams_ep,
features=feat_ep,
targets=targ_ep,
epoch_dir=epochdir,
)
print(f"Epoch: {epoch}")
Epoch: 2022-05-01_14-01-15
The information stored inside the epoch object can be easily accessed as attributes of the object. For example, to get the camera objects, you can access the cameras attribute of the epoch object.
epoch.timestamp # Gives the timestamp of the epoch
epoch.cameras # Gives a dictionary with the camera name as key and the camera object as value
epoch.images # GIves a dictionary with the camera name as key and the image object as value
epoch.features # Gives a dictionary with the camera name as key and the features object as value
epoch.targets # Gives a Target object containing both the coordinates in the object space and the projections in the image space
print(epoch.timestamp)
print(epoch.cameras)
print(epoch.images)
print(epoch.features)
print(epoch.targets)
2022-05-01 14:01:15
{'p1': Camera (f=6621.743457206283, img_size=(6012.0, 4008.0), 'p2': Camera (f=9267.892627662095, img_size=(6012.0, 4008.0)}
{'p1': Image ../data/img/p1/IMG_2637.jpg, 'p2': Image ../data/img/p2/IMG_1112.jpg}
{'p1': Features with 0 features, 'p2': Features with 0 features}
<icepy4d.core.targets.Targets object at 0x7f915b22b4f0>
Stereo Processing¶
The stereo processing is carried out for each epoch in order to find matched features, estimating camera pose, and triangulating the 3D points. The output of this step is a set of 3D points and their corresponding descriptors.
The same procedure is then iterated in a big loop for all the epoches in a multitemporal processing (refer to the multitemporal_workflow.ipynb notebook).
Feature matching with LightGlue¶
Wide-baseline feature matching is performed using the LightGlue algorithm.
Refer to the matching.ipynb notebook for more details about the matching process and explanation of the parameters.
matcher = matching.LightGlueMatcher()
matcher.match(
epoch.images[cams[0]].value,
epoch.images[cams[1]].value,
quality=matching.Quality.HIGH,
tile_selection= matching.TileSelection.PRESELECTION,
grid=[2, 3],
overlap=200,
origin=[0, 0],
do_viz_matches=True,
do_viz_tiles=True,
min_matches_per_tile = 3,
max_keypoints = 8196,
save_dir=epoch.epoch_dir / "matching",
geometric_verification=matching.GeometricVerification.PYDEGENSAC,
threshold=2,
confidence=0.9999,
)
timer.update("matching")
2023-10-03 10:30:44 | [INFO ] Running inference on device cuda
2023-10-03 10:30:45 | [INFO ] Matching by tiles... 2023-10-03 10:30:45 | [INFO ] Matching tiles by preselection tile selection 2023-10-03 10:30:46 | [INFO ] - Matching tile pair (1, 1) 2023-10-03 10:30:48 | [INFO ] - Matching tile pair (1, 4) 2023-10-03 10:30:49 | [INFO ] - Matching tile pair (2, 4) 2023-10-03 10:30:50 | [INFO ] - Matching tile pair (2, 5) 2023-10-03 10:30:52 | [INFO ] - Matching tile pair (3, 3) 2023-10-03 10:30:53 | [INFO ] - Matching tile pair (4, 3) 2023-10-03 10:30:54 | [INFO ] - Matching tile pair (4, 4) 2023-10-03 10:30:56 | [INFO ] - Matching tile pair (5, 4) 2023-10-03 10:30:57 | [INFO ] - Matching tile pair (5, 5) 2023-10-03 10:30:58 | [INFO ] Restoring full image coordinates of matches... 2023-10-03 10:30:58 | [INFO ] Matching by tile completed. 2023-10-03 10:30:58 | [INFO ] Matching done! 2023-10-03 10:30:58 | [INFO ] Performing geometric verification... 2023-10-03 10:30:59 | [INFO ] Pydegensac found 1061 inliers (36.49%) 2023-10-03 10:30:59 | [INFO ] Geometric verification done. 2023-10-03 10:31:00 | [INFO ] [Timer] | [Matching] preselection=1.626, matching=12.153, geometric_verification=0.146, Function match took 15.1186 seconds
You can now extract the matched features from the Matcher object and save them in the current Epoch object
# Define a dictionary with empty Features objects for each camera, which will be filled with the matched keypoints, descriptors and scores
f = {cam: icecore.Features() for cam in cams}
# Stack matched keypoints, descriptors and scores into Features objects
f[cams[0]].append_features_from_numpy(
x=matcher.mkpts0[:, 0],
y=matcher.mkpts0[:, 1],
descr=matcher.descriptors0,
scores=matcher.scores0,
)
f[cams[1]].append_features_from_numpy(
x=matcher.mkpts1[:, 0],
y=matcher.mkpts1[:, 1],
descr=matcher.descriptors1,
scores=matcher.scores1,
)
# Store the dictionary with the features in the Epoch object
epoch.features = f
3D Scene reconstruction¶
Relative orientation¶
First, perform Relative orientation of the two cameras by using the matched features and the a-priori camera interior orientation.
To perform the relative orientation, you have to define a RelativeOrientation object by passing first a list containing the two camera objects and then a list containing the matched features on each image. The matched features are Nx2 numpy arrary containing the x-y pixel coordinates of the matched features.
relative_orientation = RelativeOrientation([camera1, camera2], [features1, features2])
To get the pixel coordinates of the matched features as numpy arrays you can use the kpts_to_numpy() method of a Features object (that is now stored into the current Epoch object).
epoch.features[cams[0]].kpts_to_numpy()
The relative orientation is then performed by calling the estimate_pose method of the object.
You can pass some additional parameters such as the camera baseline (in meters) to scale the reconstruction.
relative_orientation.estimate_pose(scale_factor=camera_baseline)
# Compute the camera baseline from a-priori camera positions
baseline = np.linalg.norm(
cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1]
)
# Initialize RelativeOrientation class with a list containing the two
# cameras and a list contaning the matched features location on each camera.
relative_ori = sfm.RelativeOrientation(
[epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
[
epoch.features[cams[0]].kpts_to_numpy(),
epoch.features[cams[1]].kpts_to_numpy(),
],
)
# Estimate the relative orientation
relative_ori.estimate_pose(
threshold=cfg.matching.pydegensac_threshold,
confidence=0.999999,
scale_factor=baseline,
)
# Store result in camera 1 object
epoch.cameras[cams[1]] = relative_ori.cameras[1]
timer.update("relative orientation")
2023-10-03 10:33:11 | [INFO ] Relative Orientation - valid points: 963/1061 2023-10-03 10:33:11 | [INFO ] Relative orientation Succeded.
Triangulation¶
You can now triangulate the tie points (i.e., the matched features) into the object space.
Similarly as before, you have to define a Triangulation object by passing first a list containing the two camera objects and then a list containing the matched features on each image. The matched features are Nx2 numpy arrary containing the x-y pixel coordinates of the matched features.
triangulation = Triangulation([camera1, camera2], [features1, features2])
The triangulation is then performed by calling the triangulate method of the object.
triangulation.triangulate()
You can decide if you want to compute the point colors by interpolating them from one of the two images. If so, you have to pass the index of the image that you want to use for the color interpolation. For example, to interpolate colors from image 1, you can do:
triangulation.triangulate(compute_colors=True, image=epoch.images[cams[1]].value, cam_id=1)
triang = sfm.Triangulate(
[epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
[
epoch.features[cams[0]].kpts_to_numpy(),
epoch.features[cams[1]].kpts_to_numpy(),
],
)
points3d = triang.triangulate_two_views(
compute_colors=True, image=epoch.images[cams[1]].value, cam_id=1
)
# Update timer
timer.update("triangulation")
2023-10-03 10:40:02 | [INFO ] Point triangulation succeded: 1.0. 2023-10-03 10:40:02 | [INFO ] Point colors interpolated
Absolute orientation¶
Now, you can perform an absolute orientation of the current solution (i.e., cameras' exterior orientation and 3D points) by using the ground control points.
The coordinates of the two cameras are used as additional ground control points for estimating a Helmert transformation.
You need first to extract the image and object coordinates of the targets from the Target object, stored into the current Epoch. Alternatively, you can also define manually the coordinates of the targets as numpy arrays.
To extract the image coordinates of the targets, you can use the get_image_coor_by_label() method of the Target object, by passing the list of the target labels that you want to extract and the camera id they are referred to.
epoch.targets.get_image_coor_by_label(cfg.georef.targets_to_use, cam_id=id)
To get the targets object coordinates, you can use the get_object_coor_by_label() method of the Target object, by passing the list of the target labels that you want to extract.
epoch.targets.get_object_coor_by_label(cfg.georef.targets_to_use)
Eventually, you need to have 3 numpy array, all with the number of rows (the number of targets).
# Extract the image coordinates of the targets from the Targets object
image_coords = [
epoch.targets.get_image_coor_by_label(cfg.georef.targets_to_use, cam_id=id)[0] for id, cam in enumerate(cams)
]
print(f"Targets coordinates on image 0:\n{image_coords[0]}")
print(f"Targets coordinates on image 1:\n{image_coords[1]}")
obj_coords = epoch.targets.get_object_coor_by_label(cfg.georef.targets_to_use)[0]
print(f"Targets coordinates in object space:\n{obj_coords}")
Targets coordinates on image 0: [[4002.709 3543.0627] [1611.3804 1420.486 ] [4671.8179 3465.3052]] Targets coordinates on image 1: [[1003.5037 3859.1558] [5694.3535 620.8673] [1565.0667 3927.6331]] Targets coordinates in object space: [[ 49.6488 192.0875 71.7466] [-532.7409 391.02 238.8015] [ 51.1682 210.4649 70.9032]]
You can now perform the absolute orientation in a similar way as before:
# Initialize AbsoluteOrientation object with a list containing the two
abs_ori = sfm.Absolute_orientation(
(epoch.cameras[cams[0]], epoch.cameras[cams[1]]),
points3d_final=obj_coords,
image_points=image_coords,
camera_centers_world=cfg.georef.camera_centers_world,
)
# Estimate the absolute orientation transformation
T = abs_ori.estimate_transformation_linear(estimate_scale=True)
# Transform the 3D points
points3d = abs_ori.apply_transformation(points3d=points3d)
2023-10-03 10:49:29 | [INFO ] Point triangulation succeded: 1.0.
You can now save the estimated camera positions and the 3D points into the current Epoch object. The 3D coordinates of the points in the object space can be saved as a ICEpy4D Points object, as follows:
# Store the absolute orientation transformation in the camera objects
for i, cam in enumerate(cams):
epoch.cameras[cam] = abs_ori.cameras[i]
# Convert the 3D points to an icepy4d Points object
pts = icecore.Points()
pts.append_points_from_numpy(
points3d,
track_ids=epoch.features[cams[0]].get_track_ids(),
colors=triang.colors,
)
# Store the points in the Epoch object
epoch.points = pts
# Update timer
timer.update("absolute orientation")
You can save the current Epoch object as a pickle file in the previously defined epoch directory.
# Save epoch as a pickle object
if epoch.save_pickle(f"{epoch.epoch_dir}/{epoch}.pickle"):
logger.info(f"{epoch} saved successfully")
else:
logger.error(f"Unable to save {epoch}")
2023-10-03 10:50:49 | [INFO ] 2022-05-01_14-01-15 saved successfully