Multi Epoch Processing¶

Initialization¶

Let's first set up the python environment and define the configuration file

In [2]:

# Import required standard modules
import shutil
import sys
from pathlib import Path

import numpy as np

# Import required icepy4d4D modules
from icepy4d import core as icepy4d_core
from icepy4d.core import Epoch, Epoches, EpochDataMap
from icepy4d import matching
from icepy4d import sfm
from icepy4d import io
from icepy4d import utils as icepy4d_utils
from icepy4d.metashape import metashape as MS
from icepy4d.utils import initialization as inizialization

# Import required standard modules import shutil import sys from pathlib import Path import numpy as np # Import required icepy4d4D modules from icepy4d import core as icepy4d_core from icepy4d.core import Epoch, Epoches, EpochDataMap from icepy4d import matching from icepy4d import sfm from icepy4d import io from icepy4d import utils as icepy4d_utils from icepy4d.metashape import metashape as MS from icepy4d.utils import initialization as inizialization

In [3]:

# Define the path to the configuration file
CFG_FILE = "config/config_2022.yaml"

# Define the path to the configuration file CFG_FILE = "config/config_2022.yaml"

Inizialize all the required variables

In [4]:

# Parse the configuration file
cfg_file = Path(CFG_FILE)
cfg = inizialization.parse_cfg(cfg_file)

# Initialize the timer and logger
timer_global = icepy4d_utils.AverageTimer()
logger = icepy4d_utils.get_logger()

# Get the list of cameras from the configuration file
cams = cfg.cams

# Get the list of images from the configuration file
images, epoch_dict = inizialization.initialize_image_ds(cfg)

# Initialize an empty Epoches object to store the results of each epoch
epoches = Epoches(starting_epoch=cfg.proc.epoch_to_process[0])

# Parse the configuration file cfg_file = Path(CFG_FILE) cfg = inizialization.parse_cfg(cfg_file) # Initialize the timer and logger timer_global = icepy4d_utils.AverageTimer() logger = icepy4d_utils.get_logger() # Get the list of cameras from the configuration file cams = cfg.cams # Get the list of images from the configuration file images, epoch_dict = inizialization.initialize_image_ds(cfg) # Initialize an empty Epoches object to store the results of each epoch epoches = Epoches(starting_epoch=cfg.proc.epoch_to_process[0])

================================================================
ICEpy4D
Image-based Continuos monitoring of glaciers' Evolution with low-cost stereo-cameras and Deep Learning photogrammetry
2023 - Francesco Ioli - francesco.ioli@polimi.it
================================================================

An exception has occurred, use %tb to see the full traceback.

SystemExit: Configuration file does not exist! Aborting.

/home/francesco/miniconda3/envs/icepy4d/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3516: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.
  warn("To exit: use 'exit', 'quit', or Ctrl-D.", stacklevel=1)

ICEpy4D main structures¶

For the processing, you have to initialize all the required variables. This procedure is the same also for multi-epoch processing.

EpochDataMap¶

The EpochDataMap class is a critical structure for multi-epoch processing in ICEpy4D. It helps organize and manage data for each epoch, including timestamps, associated images, and time differences. EpochDataMap is a dictionary that contains the information the timestamp (i.e., date and time) of each epoch and the images belonging to that epoch. The purpose of this class is to give a east-to-use tool for associating the timestamp of each epoch to the images belonging to that epoch.

The EpochDataMap allows for automatically define the association of the images to the epochs, based on the timestamp of the images. The images are associated to the epoch with the closest timestamp. The EpochDataMap first select the master camera (by default, the first camera). Then, for each image of the master camera, it looks for the closest image taken by all the other cameras (i.e., named 'slave cameras') based on their timestamp. As the cameras may not be perfectly syncronized, a time tolerance can be defined to allow for a maximum time difference between the images of different cameras (by default, 180 seconds). If the time difference between a slave image and the master image is larger than the time tolerance, the slave image is not associated to the current epoch. It may be possible that one epoch contains only the image of the master camera (e.g., due to a disruption of the other cameras during that day), therefore a minimum number of images can be defined for an epoch to be included in the EpochDataMap (by default, 2 images). If less than the minimum number of images are associated to an epoch, the epoch is discarded.

1. The EpochDataMap structure

   The key of the dictionary is the epoch number (an integer number, starting from 0). The values of the dictionary are the EpochData class are again dictionaries that contains the timestamp of the epoch, the images associated at that epoch and the time difference between the timestamp of each camera and that of the master camera (by default, the first camera). For instance, a value of the EpochDataMap dictionary is the following:

   epoch_data_map[0] = {
       'timestamp': datetime.datetime(2023, 6, 23, 9, 59, 58), # timestamp of the first epoch
       'images': { 
           'p1': Image data/img/p1/p1_20230623_095958_IMG_0845.JPG,
           'p2': Image data/img/p2/p2_20230623_100000_IMG_0582.JPG
           }, # images of the first epoch
       'dt': {'p1': datetime.timedelta(0), 'p2': datetime.timedelta(seconds=2)} # time difference between each camera and the master camera
   }
   

   For accessing the data inside the dictionary for each epoch, you can use the dot notation for getting the timestamp, the image dictionary and the time differences.

   epoch_map[0].timestamp 
   epoch_map[0].images 
   epoch_map[0].dt
   

   The epoch timestamp (epoch_map[n].timestamp) is taken from the timestamp of the master camera, and it is stored as a datetime object.

   epoch_map[0].timestamp = datetime.datetime(2023, 6, 23, 9, 59, 58)
   

   The images of each epoch is again a dictionary with the camera names (i.e., the name of the folder containing the image sequence) as keys and Image objects as values.

   epoch_map[0].images['p1'] = Image("data/img/p1/p1_20230623_095958_IMG_0845.JPG")
   epoch_map[0].images['p2'] = Image("data/img/p2/p2_20230623_100000_IMG_0582.JPG")
   

2. How to inizialize a EpochDataMap

   Initialize the EpochDataMap object by providing the image directory and optional parameters like the time tolerance (maximum time difference allowed between images from different cameras) and minimum number of images required for an epoch to be included.

   epoch_map = EpochDataMap('path_to_image_directory', time_tolerance_sec=180, min_images=2)
   

   If not specified, the time tolerance is set to 180 seconds and the minimum number of images is set to 2.

In [6]:

# Build the EpochDataMap object find pairs of coheval images for each epoch
epoch_map = EpochDataMap(cfg.paths.image_dir, time_tolerance_sec=1200)

print(f"Epoch 0  Timestamp: {epoch_map[0].timestamp}")
print(f"\t Images: {epoch_map[0].images}")
print(f"\t Delta t from master camera (Cam 0): {epoch_map[0].dt}")

# Build the EpochDataMap object find pairs of coheval images for each epoch epoch_map = EpochDataMap(cfg.paths.image_dir, time_tolerance_sec=1200) print(f"Epoch 0 Timestamp: {epoch_map[0].timestamp}") print(f"\t Images: {epoch_map[0].images}") print(f"\t Delta t from master camera (Cam 0): {epoch_map[0].dt}")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/home/francesco/phd/icepy4d/notebooks/mutlitemporal_workflow.ipynb Cell 8 line 2
      <a href='vscode-notebook-cell:/home/francesco/phd/icepy4d/notebooks/mutlitemporal_workflow.ipynb#X44sZmlsZQ%3D%3D?line=0'>1</a> # Build the EpochDataMap object find pairs of coheval images for each epoch
----> <a href='vscode-notebook-cell:/home/francesco/phd/icepy4d/notebooks/mutlitemporal_workflow.ipynb#X44sZmlsZQ%3D%3D?line=1'>2</a> epoch_map = EpochDataMap(cfg.paths.image_dir, time_tolerance_sec=1200)
      <a href='vscode-notebook-cell:/home/francesco/phd/icepy4d/notebooks/mutlitemporal_workflow.ipynb#X44sZmlsZQ%3D%3D?line=3'>4</a> print(f"Epoch 0  Timestamp: {epoch_map[0].timestamp}")
      <a href='vscode-notebook-cell:/home/francesco/phd/icepy4d/notebooks/mutlitemporal_workflow.ipynb#X44sZmlsZQ%3D%3D?line=4'>5</a> print(f"\t Images: {epoch_map[0].images}")

NameError: name 'cfg' is not defined

In [ ]:

# Set id to process
ep = cfg.proc.epoch_to_process[0]

# Define paths to the epoch directory
epoch_name = epoch_map.get_timestamp_str(ep)
epochdir = cfg.paths.results_dir / epoch_name

# Build a dictionary of images containing the name of the cameras as keys and the image path as values
im_epoch = epoch_map.get_images(ep)

# Set id to process ep = cfg.proc.epoch_to_process[0] # Define paths to the epoch directory epoch_name = epoch_map.get_timestamp_str(ep) epochdir = cfg.paths.results_dir / epoch_name # Build a dictionary of images containing the name of the cameras as keys and the image path as values im_epoch = epoch_map.get_images(ep)

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 processing for all the epoches is then iterated in a big loop.

Load or create a new Epoch object¶

In [4]:

# Initialize a timer to measure the processing time
timer = icepy4d_utils.AverageTimer()

# Get epoch id to process
ep = cfg.proc.epoch_to_process[0]

# Define paths to the epoch directory
epoch_name = epoch_map.get_timestamp(ep)
epochdir = cfg.paths.results_dir / epoch_name

# Load an existing epoch or create a new one
if cfg.proc.load_existing_results:
    try:
        # Load existing epcoh from pickle file
        epoch = Epoch.read_pickle(epochdir / f"{epoch_name}.pickle")

    except:
        logger.error(
            f"Unable to load epoch {epoch_name} from pickle file. Creating new epoch..."
        )
        epoch = inizialization.initialize_epoch(
            cfg=cfg, images=images, epoch_id=ep, epoch_dir=epochdir
        )

else:
    # Create new epoch object
    epoch = inizialization.initialize_epoch(
        cfg=cfg, images=images, epoch_id=ep, epoch_dir=epochdir
    )

# Initialize a timer to measure the processing time timer = icepy4d_utils.AverageTimer() # Get epoch id to process ep = cfg.proc.epoch_to_process[0] # Define paths to the epoch directory epoch_name = epoch_map.get_timestamp(ep) epochdir = cfg.paths.results_dir / epoch_name # Load an existing epoch or create a new one if cfg.proc.load_existing_results: try: # Load existing epcoh from pickle file epoch = Epoch.read_pickle(epochdir / f"{epoch_name}.pickle") except: logger.error( f"Unable to load epoch {epoch_name} from pickle file. Creating new epoch..." ) epoch = inizialization.initialize_epoch( cfg=cfg, images=images, epoch_id=ep, epoch_dir=epochdir ) else: # Create new epoch object epoch = inizialization.initialize_epoch( cfg=cfg, images=images, epoch_id=ep, epoch_dir=epochdir )

Feature matching with SuperGlue¶

In [5]:

# Define matching parameters
matching_quality = matching.Quality.HIGH
tile_selection = matching.TileSelection.PRESELECTION
tiling_grid = [4, 3]
tiling_overlap = 200
geometric_verification = matching.GeometricVerification.PYDEGENSAC
geometric_verification_threshold = 1
geometric_verification_confidence = 0.9999
match_dir = epoch.epoch_dir / "matching"

# Create a new matcher object
matcher = matching.SuperGlueMatcher(cfg.matching)
matcher.match(
    epoch.images[cams[0]].value,
    epoch.images[cams[1]].value,
    quality=matching_quality,
    tile_selection=tile_selection,
    grid=tiling_grid,
    overlap=tiling_overlap,
    do_viz_matches=True,
    do_viz_tiles=False,
    save_dir=match_dir,
    geometric_verification=geometric_verification,
    threshold=geometric_verification_threshold,
    confidence=geometric_verification_confidence,
)
timer.update("matching")

# Define matching parameters matching_quality = matching.Quality.HIGH tile_selection = matching.TileSelection.PRESELECTION tiling_grid = [4, 3] tiling_overlap = 200 geometric_verification = matching.GeometricVerification.PYDEGENSAC geometric_verification_threshold = 1 geometric_verification_confidence = 0.9999 match_dir = epoch.epoch_dir / "matching" # Create a new matcher object matcher = matching.SuperGlueMatcher(cfg.matching) matcher.match( epoch.images[cams[0]].value, epoch.images[cams[1]].value, quality=matching_quality, tile_selection=tile_selection, grid=tiling_grid, overlap=tiling_overlap, do_viz_matches=True, do_viz_tiles=False, save_dir=match_dir, geometric_verification=geometric_verification, threshold=geometric_verification_threshold, confidence=geometric_verification_confidence, ) timer.update("matching")

2023-09-03 13:19:12 | [INFO    ] Running inference on device cuda
Loaded SuperPoint model
Loaded SuperGlue model ("outdoor" weights)
2023-09-03 13:19:13 | [INFO    ] Matching by tiles...
2023-09-03 13:19:13 | [INFO    ] Matching tiles by preselection tile selection
2023-09-03 13:19:14 | [INFO    ] Matching completed.
2023-09-03 13:19:14 | [INFO    ]  - Matching tile pair (3, 2)
2023-09-03 13:19:16 | [INFO    ]  - Matching tile pair (4, 7)
2023-09-03 13:19:18 | [INFO    ]  - Matching tile pair (5, 7)
2023-09-03 13:19:20 | [INFO    ]  - Matching tile pair (5, 8)
2023-09-03 13:19:23 | [INFO    ]  - Matching tile pair (6, 6)
2023-09-03 13:19:25 | [INFO    ]  - Matching tile pair (6, 9)
2023-09-03 13:19:27 | [INFO    ]  - Matching tile pair (7, 6)
2023-09-03 13:19:29 | [INFO    ]  - Matching tile pair (7, 7)
2023-09-03 13:19:31 | [INFO    ]  - Matching tile pair (7, 9)
2023-09-03 13:19:33 | [INFO    ]  - Matching tile pair (7, 10)
2023-09-03 13:19:36 | [INFO    ]  - Matching tile pair (8, 7)
2023-09-03 13:19:38 | [INFO    ]  - Matching tile pair (8, 8)
2023-09-03 13:19:40 | [INFO    ]  - Matching tile pair (8, 10)
2023-09-03 13:19:42 | [INFO    ]  - Matching tile pair (8, 11)
2023-09-03 13:19:44 | [INFO    ]  - Matching tile pair (9, 9)
2023-09-03 13:19:46 | [INFO    ]  - Matching tile pair (10, 9)
2023-09-03 13:19:49 | [INFO    ]  - Matching tile pair (10, 10)
2023-09-03 13:19:51 | [INFO    ]  - Matching tile pair (11, 10)
2023-09-03 13:19:53 | [INFO    ] Restoring full image coordinates of matches...
2023-09-03 13:19:53 | [INFO    ] Matching by tile completed.
2023-09-03 13:19:53 | [INFO    ] Matching done!
2023-09-03 13:19:53 | [INFO    ] Performing geometric verification...
2023-09-03 13:19:53 | [INFO    ] Pydegensac found 2012 inliers (36.58%)
2023-09-03 13:19:53 | [INFO    ] Geometric verification done.
2023-09-03 13:19:55 | [INFO    ] [Timer] | [Matching] preselection=1.014, matching=39.122, geometric_verification=0.449, 
Function match took 41.6525 seconds

Extract the matched features from the Matcher object and save them in the current Epoch object

In [ ]:

# Define a dictionary with empty Features objects for each camera, which will be filled with the matched keypoints, descriptors and scores
f = {cam: icepy4d_classes.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

# Define a dictionary with empty Features objects for each camera, which will be filled with the matched keypoints, descriptors and scores f = {cam: icepy4d_classes.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

Scene reconstruction¶

First, perform Relative orientation of the two cameras by using the matched features and the a-priori camera interior orientation.

In [7]:

# 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(),
    ],
)
relative_ori.estimate_pose(
    threshold=cfg.matching.pydegensac_threshold,
    confidence=0.999999,
    scale_factor=np.linalg.norm(
        cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1]
    ),
)
# Store result in camera 1 object
epoch.cameras[cams[1]] = relative_ori.cameras[1]

cfg.georef.camera_centers_world

# 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(), ], ) relative_ori.estimate_pose( threshold=cfg.matching.pydegensac_threshold, confidence=0.999999, scale_factor=np.linalg.norm( cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1] ), ) # Store result in camera 1 object epoch.cameras[cams[1]] = relative_ori.cameras[1] cfg.georef.camera_centers_world

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[7], line 10
      1 # Initialize RelativeOrientation class with a list containing the two
      2 # cameras and a list contaning the matched features location on each camera.
      3 relative_ori = sfm.RelativeOrientation(
      4     [epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
      5     [
   (...)
      8     ],
      9 )
---> 10 relative_ori.estimate_pose(
     11     threshold=cfg.matching.pydegensac_threshold,
     12     confidence=0.999999,
     13     scale_factor=np.linalg.norm(
     14         cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1]
     15     ),
     16 )
     17 # Store result in camera 1 object
     18 epoch.cameras[cams[1]] = relative_ori.cameras[1]

File ~/phd/icepy4d/src/icepy4d/sfm/two_view_geometry.py:81, in RelativeOrientation.estimate_pose(self, threshold, confidence, scale_factor)
     75 assert self.cameras[0].extrinsics is not None, print(
     76     "Extrinsics matrix is not available for camera 0. Please, compute it before running RelativeOrientation estimation."
     77 )
     79 # Estimate Realtive Pose with Essential Matrix
     80 # R, t make up a tuple that performs a change of basis from the first camera's coordinate system to the second camera's coordinate system.
---> 81 R, t, valid = estimate_pose(
     82     self.features[0],
     83     self.features[1],
     84     self.cameras[0].K,
     85     self.cameras[1].K,
     86     thresh=threshold,
     87     conf=confidence,
     88 )
     89 logging.info(f"Relative Orientation - valid points: {valid.sum()}/{len(valid)}")
     91 # If the scaling factor is given, scale the stereo model

TypeError: cannot unpack non-iterable NoneType object

In [16]:

camera_baseline = np.linalg.norm(
        cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1]
    )
image = epoch.images[cams[0]].value

camera_baseline = np.linalg.norm( cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1] ) image = epoch.images[cams[0]].value

In [17]:

# Relative orientation
feature0 = epoch.features[cams[0]].kpts_to_numpy()
feature1 = epoch.features[cams[1]].kpts_to_numpy()

relative_ori = sfm.RelativeOrientation(
    [epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
    [feature0, feature1],
)
relative_ori.estimate_pose(scale_factor=camera_baseline)
epoch.cameras[cams[1]] = relative_ori.cameras[1]

# Triangulation
triang = sfm.Triangulate(
    [epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
    [feature0, feature1],
)
points3d = triang.triangulate_two_views(
    compute_colors=True, image=image, cam_id=1
)

# Relative orientation feature0 = epoch.features[cams[0]].kpts_to_numpy() feature1 = epoch.features[cams[1]].kpts_to_numpy() relative_ori = sfm.RelativeOrientation( [epoch.cameras[cams[0]], epoch.cameras[cams[1]]], [feature0, feature1], ) relative_ori.estimate_pose(scale_factor=camera_baseline) epoch.cameras[cams[1]] = relative_ori.cameras[1] # Triangulation triang = sfm.Triangulate( [epoch.cameras[cams[0]], epoch.cameras[cams[1]]], [feature0, feature1], ) points3d = triang.triangulate_two_views( compute_colors=True, image=image, cam_id=1 )

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[17], line 9
      3 feature1 = epoch.features[cams[1]].kpts_to_numpy()
      5 relative_ori = sfm.RelativeOrientation(
      6     [epoch.cameras[cams[0]], epoch.cameras[cams[1]]],
      7     [feature0, feature1],
      8 )
----> 9 relative_ori.estimate_pose(
     10     scale_factor=camera_baseline,
     11 )
     12 epoch.cameras[cams[1]] = relative_ori.cameras[1]
     14 # Traingulation

File ~/phd/icepy4d/src/icepy4d/sfm/two_view_geometry.py:81, in RelativeOrientation.estimate_pose(self, threshold, confidence, scale_factor)
     75 assert self.cameras[0].extrinsics is not None, print(
     76     "Extrinsics matrix is not available for camera 0. Please, compute it before running RelativeOrientation estimation."
     77 )
     79 # Estimate Realtive Pose with Essential Matrix
     80 # R, t make up a tuple that performs a change of basis from the first camera's coordinate system to the second camera's coordinate system.
---> 81 R, t, valid = estimate_pose(
     82     self.features[0],
     83     self.features[1],
     84     self.cameras[0].K,
     85     self.cameras[1].K,
     86     thresh=threshold,
     87     conf=confidence,
     88 )
     89 logging.info(f"Relative Orientation - valid points: {valid.sum()}/{len(valid)}")
     91 # If the scaling factor is given, scale the stereo model

TypeError: cannot unpack non-iterable NoneType object

Triangulate points into the object space

In [14]:

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=images[cams[1]].read_image(ep).value, cam_id=1
)

# Update timer
timer.update("triangulation")

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=images[cams[1]].read_image(ep).value, cam_id=1 ) # Update timer timer.update("triangulation")

2023-08-29 17:34:37 | [INFO    ] Point triangulation succeded: 1.0.
2023-08-29 17:34:37 | [INFO    ] Point colors interpolated

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.

In [13]:

# Get targets available in all cameras. The Labels of valid targets are returned as second element by the get_image_coor_by_label() method
valid_targets = epoch.targets.get_image_coor_by_label(
    cfg.georef.targets_to_use, cam_id=0
)[1]

# Check if the same targets are available in all cameras
for id in range(1, len(cams)):
    assert (
        valid_targets
        == epoch.targets.get_image_coor_by_label(
            cfg.georef.targets_to_use, cam_id=id
        )[1]
    ), f"""epoch {ep} - {epoch_map.get_timestamp(ep)}: 
    Different targets found in image {id} - {images[cams[id]][ep]}"""

# Check if there are enough targets
assert len(valid_targets) > 1, f"Not enough targets found in epoch {ep}"

# If not all the targets defined in the config file are found, log a warning and use only the valid targets
if valid_targets != cfg.georef.targets_to_use:
    logger.warning(f"Not all targets found. Using onlys {valid_targets}")

# Get image and object coordinates of valid targets
image_coords = [
    epoch.targets.get_image_coor_by_label(valid_targets, cam_id=id)[0]
    for id, cam in enumerate(cams)
]
obj_coords = epoch.targets.get_object_coor_by_label(valid_targets)[0]

# Perform absolute orientation
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,
)
T = abs_ori.estimate_transformation_linear(estimate_scale=True)
points3d = abs_ori.apply_transformation(points3d=points3d)
for i, cam in enumerate(cams):
    epoch.cameras[cam] = abs_ori.cameras[i]

# Convert the 3D points to an icepy4d Points object
pts = icepy4d_classes.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")

# Get targets available in all cameras. The Labels of valid targets are returned as second element by the get_image_coor_by_label() method valid_targets = epoch.targets.get_image_coor_by_label( cfg.georef.targets_to_use, cam_id=0 )[1] # Check if the same targets are available in all cameras for id in range(1, len(cams)): assert ( valid_targets == epoch.targets.get_image_coor_by_label( cfg.georef.targets_to_use, cam_id=id )[1] ), f"""epoch {ep} - {epoch_map.get_timestamp(ep)}: Different targets found in image {id} - {images[cams[id]][ep]}""" # Check if there are enough targets assert len(valid_targets) > 1, f"Not enough targets found in epoch {ep}" # If not all the targets defined in the config file are found, log a warning and use only the valid targets if valid_targets != cfg.georef.targets_to_use: logger.warning(f"Not all targets found. Using onlys {valid_targets}") # Get image and object coordinates of valid targets image_coords = [ epoch.targets.get_image_coor_by_label(valid_targets, cam_id=id)[0] for id, cam in enumerate(cams) ] obj_coords = epoch.targets.get_object_coor_by_label(valid_targets)[0] # Perform absolute orientation 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, ) T = abs_ori.estimate_transformation_linear(estimate_scale=True) points3d = abs_ori.apply_transformation(points3d=points3d) for i, cam in enumerate(cams): epoch.cameras[cam] = abs_ori.cameras[i] # Convert the 3D points to an icepy4d Points object pts = icepy4d_classes.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")

2023-08-29 17:34:18 | [WARNING ] Warning: target T2 is not present on camera 0.
2023-08-29 17:34:18 | [WARNING ] Warning: target F10_2 is not present on camera 0.
2023-08-29 17:34:18 | [WARNING ] Warning: target T2 is not present on camera 1.
2023-08-29 17:34:18 | [WARNING ] Warning: target F10_2 is not present on camera 1.
2023-08-29 17:34:18 | [WARNING ] Not all targets found. Using onlys ['F2', 'F12', 'F13']
2023-08-29 17:34:18 | [INFO    ] Point triangulation succeded: 1.0.

Save the current Epoch object as a pickle file.

In [9]:

# 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}")

# 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-08-29 18:32:30 | [INFO    ] 2022-05-01_14:01:15 saved successfully

Big loop over the epoches¶

Stack all the processing of a single epoch into a function and iterate over all the epoches

In [7]:

# Define processing for single epoch
def process_epoch(epoch, cfg, timer) -> Epoch:

    cams = cfg.cams
    epochdir = epoch.epoch_dir
    match_dir = epochdir / "matching"
   
    # Matching
    matching_quality = matching.Quality.HIGH
    tile_selection = matching.TileSelection.PRESELECTION
    tiling_grid = [4, 3]
    tiling_overlap = 200
    geometric_verification = matching.GeometricVerification.PYDEGENSAC
    geometric_verification_threshold = 1
    geometric_verification_confidence = 0.9999
    matcher = matching.SuperGlueMatcher(cfg.matching)    
    matcher.match(
        epoch.images[cams[0]].value,
        epoch.images[cams[1]].value,
        quality=matching_quality,
        tile_selection=tile_selection,
        grid=tiling_grid,
        overlap=tiling_overlap,
        do_viz_matches=True,
        do_viz_tiles=False,
        save_dir=match_dir,
        geometric_verification=geometric_verification,
        threshold=geometric_verification_threshold,
        confidence=geometric_verification_confidence,
    )
    f = {cam: icepy4d_classes.Features() for cam in cams}
    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,
    )
    epoch.features = f
    timer.update("matching")
    
    # Relative orientation
    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(),
    ],
    )
    relative_ori.estimate_pose(
    threshold=cfg.matching.pydegensac_threshold,
    confidence=0.999999,
    scale_factor=np.linalg.norm(
        cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1]
    ),
    )
    epoch.cameras[cams[1]] = relative_ori.cameras[1]
    timer.update("relative orientation")

    # Triangulation
    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=images[cams[1]].read_image(ep).value, cam_id=1
    )
    timer.update("triangulation")    

    # Absolute orientation
    valid_targets = epoch.targets.get_image_coor_by_label(
        cfg.georef.targets_to_use, cam_id=0
    )[1]
    for id in range(1, len(cams)):
        assert (
            valid_targets
            == epoch.targets.get_image_coor_by_label(
                cfg.georef.targets_to_use, cam_id=id
            )[1]
        ), f"""epoch {ep} - {epoch_map.get_timestamp(ep)}: 
        Different targets found in image {id} - {images[cams[id]][ep]}"""
    assert len(valid_targets) > 1, f"Not enough targets found in epoch {ep}"
    if valid_targets != cfg.georef.targets_to_use:
        logger.warning(f"Not all targets found. Using onlys {valid_targets}")

    image_coords = [
        epoch.targets.get_image_coor_by_label(valid_targets, cam_id=id)[0]
        for id, cam in enumerate(cams)
    ]
    obj_coords = epoch.targets.get_object_coor_by_label(valid_targets)[0]

    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,
    )
    T = abs_ori.estimate_transformation_linear(estimate_scale=True)
    points3d = abs_ori.apply_transformation(points3d=points3d)
    for i, cam in enumerate(cams):
        epoch.cameras[cam] = abs_ori.cameras[i]

    pts = icepy4d_classes.Points()
    pts.append_points_from_numpy(
        points3d,
        track_ids=epoch.features[cams[0]].get_track_ids(),
        colors=triang.colors,
    )
    epoch.points = pts
    timer.update("absolute orientation") 

    # 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}")

    return epoch

# Define processing for single epoch def process_epoch(epoch, cfg, timer) -> Epoch: cams = cfg.cams epochdir = epoch.epoch_dir match_dir = epochdir / "matching" # Matching matching_quality = matching.Quality.HIGH tile_selection = matching.TileSelection.PRESELECTION tiling_grid = [4, 3] tiling_overlap = 200 geometric_verification = matching.GeometricVerification.PYDEGENSAC geometric_verification_threshold = 1 geometric_verification_confidence = 0.9999 matcher = matching.SuperGlueMatcher(cfg.matching) matcher.match( epoch.images[cams[0]].value, epoch.images[cams[1]].value, quality=matching_quality, tile_selection=tile_selection, grid=tiling_grid, overlap=tiling_overlap, do_viz_matches=True, do_viz_tiles=False, save_dir=match_dir, geometric_verification=geometric_verification, threshold=geometric_verification_threshold, confidence=geometric_verification_confidence, ) f = {cam: icepy4d_classes.Features() for cam in cams} 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, ) epoch.features = f timer.update("matching") # Relative orientation 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(), ], ) relative_ori.estimate_pose( threshold=cfg.matching.pydegensac_threshold, confidence=0.999999, scale_factor=np.linalg.norm( cfg.georef.camera_centers_world[0] - cfg.georef.camera_centers_world[1] ), ) epoch.cameras[cams[1]] = relative_ori.cameras[1] timer.update("relative orientation") # Triangulation 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=images[cams[1]].read_image(ep).value, cam_id=1 ) timer.update("triangulation") # Absolute orientation valid_targets = epoch.targets.get_image_coor_by_label( cfg.georef.targets_to_use, cam_id=0 )[1] for id in range(1, len(cams)): assert ( valid_targets == epoch.targets.get_image_coor_by_label( cfg.georef.targets_to_use, cam_id=id )[1] ), f"""epoch {ep} - {epoch_map.get_timestamp(ep)}: Different targets found in image {id} - {images[cams[id]][ep]}""" assert len(valid_targets) > 1, f"Not enough targets found in epoch {ep}" if valid_targets != cfg.georef.targets_to_use: logger.warning(f"Not all targets found. Using onlys {valid_targets}") image_coords = [ epoch.targets.get_image_coor_by_label(valid_targets, cam_id=id)[0] for id, cam in enumerate(cams) ] obj_coords = epoch.targets.get_object_coor_by_label(valid_targets)[0] 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, ) T = abs_ori.estimate_transformation_linear(estimate_scale=True) points3d = abs_ori.apply_transformation(points3d=points3d) for i, cam in enumerate(cams): epoch.cameras[cam] = abs_ori.cameras[i] pts = icepy4d_classes.Points() pts.append_points_from_numpy( points3d, track_ids=epoch.features[cams[0]].get_track_ids(), colors=triang.colors, ) epoch.points = pts timer.update("absolute orientation") # 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}") return epoch

In [8]:

epoch = process_epoch(epoch, cfg, timer)

epoch = process_epoch(epoch, cfg, timer)

2023-08-29 18:40:56 | [INFO    ] Running inference on device cuda
Loaded SuperPoint model
Loaded SuperGlue model ("outdoor" weights)
2023-08-29 18:40:57 | [INFO    ] Matching by tiles...
2023-08-29 18:40:57 | [INFO    ] Matching tiles by preselection tile selection
2023-08-29 18:40:57 | [INFO    ] Matching completed.
2023-08-29 18:40:57 | [INFO    ]  - Matching tile pair (3, 2)
2023-08-29 18:40:59 | [INFO    ]  - Matching tile pair (4, 7)
2023-08-29 18:41:02 | [INFO    ]  - Matching tile pair (5, 7)
2023-08-29 18:41:04 | [INFO    ]  - Matching tile pair (5, 8)
2023-08-29 18:41:06 | [INFO    ]  - Matching tile pair (6, 6)
2023-08-29 18:41:09 | [INFO    ]  - Matching tile pair (6, 9)
2023-08-29 18:41:11 | [INFO    ]  - Matching tile pair (7, 6)
2023-08-29 18:41:13 | [INFO    ]  - Matching tile pair (7, 7)
2023-08-29 18:41:16 | [INFO    ]  - Matching tile pair (7, 9)
2023-08-29 18:41:18 | [INFO    ]  - Matching tile pair (7, 10)
2023-08-29 18:41:20 | [INFO    ]  - Matching tile pair (8, 7)
2023-08-29 18:41:23 | [INFO    ]  - Matching tile pair (8, 8)
2023-08-29 18:41:25 | [INFO    ]  - Matching tile pair (8, 10)
2023-08-29 18:41:27 | [INFO    ]  - Matching tile pair (8, 11)
2023-08-29 18:41:30 | [INFO    ]  - Matching tile pair (9, 9)
2023-08-29 18:41:32 | [INFO    ]  - Matching tile pair (10, 9)
2023-08-29 18:41:34 | [INFO    ]  - Matching tile pair (10, 10)
2023-08-29 18:41:37 | [INFO    ]  - Matching tile pair (11, 10)
2023-08-29 18:41:39 | [INFO    ] Restoring full image coordinates of matches...
2023-08-29 18:41:39 | [INFO    ] Matching by tile completed.
2023-08-29 18:41:39 | [INFO    ] Matching done!
2023-08-29 18:41:39 | [INFO    ] Performing geometric verification...
2023-08-29 18:41:40 | [INFO    ] Pydegensac found 2011 inliers (36.56%)
2023-08-29 18:41:40 | [INFO    ] Geometric verification done.
2023-08-29 18:41:41 | [INFO    ] [Timer] | [Matching] preselection=0.314, matching=42.250, geometric_verification=0.390, 
Function match took 44.0152 seconds
2023-08-29 18:41:41 | [INFO    ] Relative Orientation - valid points: 1874/2011
2023-08-29 18:41:41 | [INFO    ] Relative orientation Succeded.
2023-08-29 18:41:41 | [INFO    ] Point triangulation succeded: 1.0.
2023-08-29 18:41:41 | [INFO    ] Point colors interpolated
2023-08-29 18:41:41 | [WARNING ] Warning: target T2 is not present on camera 0.
2023-08-29 18:41:41 | [WARNING ] Warning: target F10_2 is not present on camera 0.
2023-08-29 18:41:41 | [WARNING ] Warning: target T2 is not present on camera 1.
2023-08-29 18:41:41 | [WARNING ] Warning: target F10_2 is not present on camera 1.
2023-08-29 18:41:41 | [WARNING ] Not all targets found. Using onlys ['F2', 'F12', 'F13']
2023-08-29 18:41:41 | [INFO    ] Point triangulation succeded: 1.0.
2023-08-29 18:41:41 | [INFO    ] 2022-05-01_14:01:15 saved successfully

Warning: License lost, attempting to recover (license status -21)
Warning: Failed to recover license (attempt 1 of 5, license status -21)
Warning: Failed to recover license (attempt 2 of 5, license status -21)
Warning: Failed to recover license (attempt 3 of 5, license status -21)
Warning: Failed to recover license (attempt 4 of 5, license status -21)
Warning: Failed to recover license (attempt 5 of 5, license status -43)
License successfully recovered

In [ ]:

# Add epoch to epoches object
epoches.add_epoch(epoch)

# Add epoch to epoches object epoches.add_epoch(epoch)

In [ ]:

logger.info("------------------------------------------------------")
logger.info("Processing started:")
timer = icepy4d_utils.AverageTimer()
iter = 0  # necessary only for printing the number of processed iteration

logger.info("------------------------------------------------------") logger.info("Processing started:") timer = icepy4d_utils.AverageTimer() iter = 0 # necessary only for printing the number of processed iteration