#!/usr/bin/env python
# coding: utf8
#
# Copyright (c) 2020 Centre National d'Etudes Spatiales (CNES).
#
# This file is part of CARS
# (see https://github.com/CNES/cars).
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# pylint: disable=too-many-lines
"""
this module contains the LineOfSightIntersection application class.
"""
# Standard imports
import logging
import os
from typing import Dict, Tuple
# Third party imports
import numpy as np
import xarray as xr
from json_checker import Checker
from shareloc.geofunctions.rectification_grid import RectificationGrid
# CARS imports
import cars.orchestrator.orchestrator as ocht
from cars.applications import application_constants
from cars.applications.grid_generation import grid_generation_algo
from cars.applications.triangulation import (
pc_transform,
triangulation_algo,
triangulation_constants,
)
from cars.applications.triangulation import (
triangulation_wrappers as triangulation_wrap,
)
from cars.applications.triangulation.abstract_triangulation_app import (
Triangulation,
)
from cars.core import constants as cst
from cars.core import inputs, projection, tiling
from cars.core.utils import safe_makedirs
from cars.data_structures import cars_dataset
from cars.orchestrator.cluster.log_wrapper import cars_profile
from cars.pipelines.parameters import sensor_inputs_constants as sens_cst
[docs]
class LineOfSightIntersection(
Triangulation, short_name="line_of_sight_intersection"
):
"""
Triangulation
"""
def __init__(self, conf=None):
"""
Init function of Triangulation
:param conf: configuration for triangulation
:return: an application_to_use object
"""
super().__init__(conf=conf)
# check conf
self.used_method = self.used_config["method"]
self.snap_to_img1 = self.used_config["snap_to_img1"]
# Saving files
self.save_intermediate_data = self.used_config["save_intermediate_data"]
self.z_inf_sup_bh_approximation = self.used_config[
"z_inf_sup_bh_approximation"
]
# global value for left image to check if snap_to_img1 can
# be applied : Need than same application object is run
# for all pairs
self.ref_left_img = None
# Init orchestrator
self.orchestrator = None
[docs]
def check_conf(self, conf):
"""
Check configuration
:param conf: configuration to check
:type conf: dict
:return: overloaded configuration
:rtype: dict
"""
# init conf
if conf is not None:
overloaded_conf = conf.copy()
else:
conf = {}
overloaded_conf = {}
# Overload conf
overloaded_conf["method"] = conf.get(
"method", "line_of_sight_intersection"
)
overloaded_conf["snap_to_img1"] = conf.get("snap_to_img1", False)
# Saving files
overloaded_conf["save_intermediate_data"] = conf.get(
"save_intermediate_data", False
)
overloaded_conf["z_inf_sup_bh_approximation"] = conf.get(
"z_inf_sup_bh_approximation", True
)
triangulation_schema = {
"method": str,
"snap_to_img1": bool,
"save_intermediate_data": bool,
"z_inf_sup_bh_approximation": bool,
}
# Check conf
checker = Checker(triangulation_schema)
checker.validate(overloaded_conf)
return overloaded_conf
# pylint: disable=too-many-positional-arguments
@cars_profile(name="Save triangulation output", interval=0.5)
def save_triangulation_output( # noqa: C901 function is too complex
self,
epi_or_sens_point_cloud,
sensor_image_left,
output_dir,
dump_dir=None,
performance_maps_to_generate=None,
save_output_coordinates=True,
save_output_color=True,
save_output_classification=False,
save_output_filling=False,
save_output_performance_map=False,
save_output_ambiguity=False,
save_output_edges=False,
save_residues=False,
):
"""
Save the triangulation output. The different TIFs composing the depth
map are written to the output directory. Auxiliary products can be
requested or not using the parameters. A dump directory can also be
provided to write any additionnal files that have not been written
to the output directory (because they are not part of the depth map
definition, or because they have not been requested).
:param epipolar_point_cloud: tiled epipolar left image
:type epipolar_point_cloud: CarsDataset
:param sensor_image_left: tiled sensor left image
Dict Must contain keys : "image", "texture", "geomodel",
"no_data", "mask". Paths must be absolutes
:type sensor_image_left: CarsDataset
:param output_dir: directory to write triangulation output depth
map.
:type output_dir: None or str
:param dump_dir: folder used as dump directory for current pair, None to
deactivate intermediate data writing
:type dump_dir: str
:param performance_maps_to_generate: None or list containing
"risk" or "intervals"
:type performance_maps_to_generate: None or list containing
"risk" or "intervals"
:param save_output_coordinates: Save X, Y and Z coords in output_dir
:type save_output_coordinates: bool
:param save_output_color: Save color depth map in output_dir
:type save_output_color: bool
:param save_output_classification: Save classification depth map in
output_dir
:type save_output_classification: bool
:param save_output_filling: Save filling depth map in output_dir
:type save_output_filling: bool
:param save_output_performance_map: Save performance map in output_dir
:type save_output_performance_map: bool
:param save_output_ambiguity: Save ambiguity in output_dir
:type save_output_ambiguity: bool
:param save_output_edges: Save edges in output_dir
:type save_output_edges: bool
"""
if dump_dir:
dump_dir = os.path.join(dump_dir, "tif")
safe_makedirs(dump_dir)
# Propagate color type in output file
color_type = inputs.rasterio_get_image_type(
sensor_image_left[sens_cst.INPUT_IMG]["bands"]["b0"]["path"]
)
if output_dir is None:
output_dir = dump_dir
if output_dir is not None:
output_dir = os.path.join(output_dir, "tif")
os.makedirs(output_dir, exist_ok=True)
if save_output_coordinates or dump_dir:
coords_output_dir = (
output_dir if save_output_coordinates else dump_dir
)
self.orchestrator.add_to_save_lists(
os.path.join(coords_output_dir, "X.tif"),
cst.X,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_x",
dtype=np.float64,
)
self.orchestrator.add_to_save_lists(
os.path.join(coords_output_dir, "Y.tif"),
cst.Y,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_y",
dtype=np.float64,
)
self.orchestrator.add_to_save_lists(
os.path.join(coords_output_dir, "Z.tif"),
cst.Z,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_z",
dtype=np.float64,
)
if save_residues or dump_dir:
residues_output_dir = output_dir if save_residues else dump_dir
self.orchestrator.add_to_save_lists(
os.path.join(residues_output_dir, "intersections_residues.tif"),
cst.INDEX_DEPTH_MAP_INTERSECTIONS_RESIDUES,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_intersections_residues",
dtype=np.float64,
)
if save_output_color or dump_dir:
color_output_dir = output_dir if save_output_color else dump_dir
self.orchestrator.add_to_save_lists(
os.path.join(color_output_dir, "image.tif"),
cst.EPI_TEXTURE,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_color",
dtype=color_type,
)
if save_output_performance_map or dump_dir:
map_output_dir = (
output_dir if save_output_performance_map else dump_dir
)
performance_map_tail_path = "performance_map.tif"
if "intervals" in performance_maps_to_generate:
tail_path = performance_map_tail_path
if len(performance_maps_to_generate) == 2:
tail_path = "performance_map_from_intervals.tif"
self.orchestrator.add_to_save_lists(
os.path.join(map_output_dir, tail_path),
cst.POINT_CLOUD_PERFORMANCE_MAP_FROM_INTERVALS,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_performance_map_from_intervals",
optional_data=True,
dtype=np.float64,
)
if "risk" in performance_maps_to_generate:
tail_path = performance_map_tail_path
if len(performance_maps_to_generate) == 2:
tail_path = "performance_map_from_risk.tif"
self.orchestrator.add_to_save_lists(
os.path.join(map_output_dir, tail_path),
cst.POINT_CLOUD_PERFORMANCE_MAP_FROM_RISK,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_performance_map_from_risk",
optional_data=True,
dtype=np.float64,
)
if save_output_ambiguity or dump_dir:
map_output_dir = output_dir if save_output_ambiguity else dump_dir
self.orchestrator.add_to_save_lists(
os.path.join(map_output_dir, "ambiguity.tif"),
cst.EPI_AMBIGUITY,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_ambiguity",
optional_data=True,
dtype=np.float64,
)
if save_output_classification or dump_dir:
classif_output_dir = (
output_dir if save_output_classification else dump_dir
)
self.orchestrator.add_to_save_lists(
os.path.join(classif_output_dir, "classification.tif"),
cst.EPI_CLASSIFICATION,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_classification",
optional_data=True,
dtype=np.uint8,
)
if save_output_filling or dump_dir:
filling_output_dir = output_dir if save_output_filling else dump_dir
self.orchestrator.add_to_save_lists(
os.path.join(filling_output_dir, "filling.tif"),
cst.EPI_FILLING,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_filling",
optional_data=True,
dtype=np.uint8,
nodata=255,
)
if save_output_edges or dump_dir:
edges_output_dir = output_dir if save_output_edges else dump_dir
self.orchestrator.add_to_save_lists(
os.path.join(edges_output_dir, "edges_mask.tif"),
cst.EPI_EDGES_MASK,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_edges_mask",
dtype=np.uint8,
optional_data=True,
)
self.orchestrator.add_to_save_lists(
os.path.join(edges_output_dir, "edges_depth_map.tif"),
cst.EPI_EDGES_DEPTH_MAP,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_edges_depth_map",
dtype=np.float32,
optional_data=True,
)
self.orchestrator.add_to_save_lists(
os.path.join(edges_output_dir, "edges_normals.tif"),
cst.EPI_EDGES_NORMALS,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_edges_normals",
dtype=np.float32,
optional_data=True,
)
self.orchestrator.add_to_save_lists(
os.path.join(edges_output_dir, "edges_tile_id.tif"),
cst.EPI_EDGES_TILE_ID,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_edges_tile_id",
dtype=np.uint8,
optional_data=True,
)
if dump_dir and performance_maps_to_generate is not None:
if "intervals" in performance_maps_to_generate:
self.orchestrator.add_to_save_lists(
os.path.join(dump_dir, "Z_inf_from_intervals.tif"),
cst.POINT_CLOUD_LAYER_INF_FROM_INTERVALS,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_z_inf_from_intervals",
)
self.orchestrator.add_to_save_lists(
os.path.join(dump_dir, "Z_sup_from_intervals.tif"),
cst.POINT_CLOUD_LAYER_SUP_FROM_INTERVALS,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_z_sup_from_intervals",
)
if "risk" in performance_maps_to_generate:
self.orchestrator.add_to_save_lists(
os.path.join(dump_dir, "Z_inf_from_risk.tif"),
cst.POINT_CLOUD_LAYER_INF_FROM_RISK,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_z_inf_from_risk",
)
self.orchestrator.add_to_save_lists(
os.path.join(dump_dir, "Z_sup_from_risk.tif"),
cst.POINT_CLOUD_LAYER_SUP_FROM_RISK,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_z_sup_from_risk",
)
if dump_dir:
self.orchestrator.add_to_save_lists(
os.path.join(dump_dir, "corr_mask.tif"),
cst.POINT_CLOUD_CORR_MSK,
epi_or_sens_point_cloud,
cars_ds_name="depth_map_corr_msk",
optional_data=True,
)
@cars_profile(name="Fill index", interval=0.5)
def fill_index( # pylint: disable=too-many-positional-arguments
self,
save_output_coordinates=True,
save_output_color=True,
save_output_classification=False,
save_output_filling=False,
save_output_performance_map=False,
save_output_ambiguity=False,
save_output_edges=False,
save_residues=False,
pair_key="PAIR_0",
):
"""
Fill depth map index for current pair, according to which product
should be saved
:param save_output_coordinates: Save X, Y and Z coords in output_dir
:type save_output_coordinates: bool
:param save_output_color: Save color depth map in output_dir
:type save_output_color: bool
:param save_output_classification: Save classification depth map in
output_dir
:type save_output_classification: bool
:param save_output_filling: Save filling depth map in output_dir
:type save_output_filling: bool
:param save_output_performance_map: Save performance map in output_dir
:type save_output_performance_map: bool
:param save_output_ambiguity: Save ambiguity in output_dir
:type save_output_ambiguity: bool
:param pair_key: name of the current pair
:type pair_key: str
"""
# index file for this depth map
index = {}
if save_output_coordinates:
index[cst.INDEX_DEPTH_MAP_X] = os.path.join(pair_key, "X.tif")
index[cst.INDEX_DEPTH_MAP_Y] = os.path.join(pair_key, "Y.tif")
index[cst.INDEX_DEPTH_MAP_Z] = os.path.join(pair_key, "Z.tif")
if save_output_color:
index[cst.INDEX_DEPTH_MAP_COLOR] = os.path.join(
pair_key, "image.tif"
)
if save_residues:
index[cst.INDEX_DEPTH_MAP_INTERSECTIONS_RESIDUES] = os.path.join(
pair_key, "intersections_residues.tif"
)
if save_output_performance_map:
index[cst.INDEX_DEPTH_MAP_PERFORMANCE_MAP] = os.path.join(
pair_key, "performance_map.tif"
)
if save_output_ambiguity:
index[cst.INDEX_DEPTH_MAP_AMBIGUITY] = os.path.join(
pair_key, "ambiguity.tif"
)
if save_output_classification:
index[cst.INDEX_DEPTH_MAP_CLASSIFICATION] = os.path.join(
pair_key, "classification.tif"
)
if save_output_filling:
index[cst.INDEX_DEPTH_MAP_FILLING] = os.path.join(
pair_key, "filling.tif"
)
if save_output_edges:
index[cst.INDEX_DEPTH_MAP_EDGES_MASK] = os.path.join(
pair_key, "edges_mask.tif"
)
index[cst.INDEX_DEPTH_MAP_EDGES_NORMALS] = os.path.join(
pair_key, "edges_normals.tif"
)
index[cst.INDEX_DEPTH_MAP_EDGES_DEPTH_MAP] = os.path.join(
pair_key, "edges_depth_map.tif"
)
index[cst.INDEX_DEPTH_MAP_EDGES_TILE_ID] = os.path.join(
pair_key, "edges_tile_id.tif"
)
# update orchestrator index if it has been filled
if index:
# Add epsg code (always lon/lat in triangulation)
index[cst.INDEX_DEPTH_MAP_EPSG] = 4326
self.orchestrator.update_index({"point_cloud": {pair_key: index}})
[docs]
def create_point_cloud_directories(
self, pair_dump_dir, point_cloud_dir, point_cloud
):
"""
Set and create directories for point cloud disk output (laz and csv)
The function return None path if the point cloud should not be saved
:param pair_dump_dir: folder used as dump directory for current pair
:type pair_dump_dir: str
:param point_cloud_dir: folder used for laz official product directory
:type point_cloud_dir: str
:param point_cloud: input point cloud (for orchestrator registration)
:type point_cloud: Dataset
"""
csv_pc_dir_name = None
if self.save_intermediate_data:
csv_pc_dir_name = os.path.join(pair_dump_dir, "csv")
safe_makedirs(csv_pc_dir_name)
self.orchestrator.add_to_compute_lists(
point_cloud, cars_ds_name="point_cloud_csv"
)
laz_pc_dir_name = None
if self.save_intermediate_data or point_cloud_dir is not None:
if point_cloud_dir is not None:
laz_pc_dir_name = os.path.join(point_cloud_dir, "laz")
os.makedirs(laz_pc_dir_name, exist_ok=True)
else:
laz_pc_dir_name = os.path.join(pair_dump_dir, "laz")
safe_makedirs(laz_pc_dir_name)
self.orchestrator.add_to_compute_lists(
point_cloud, cars_ds_name="point_cloud_laz"
)
return csv_pc_dir_name, laz_pc_dir_name
[docs]
def run( # pylint: disable=too-many-positional-arguments # noqa: C901
self,
sensor_image_left,
sensor_images_right,
geometry_plugin,
sensor_matches=None,
grid_left=None,
grid_right=None,
epipolar_disparity_map=None,
epipolar_image=None,
epsg=None,
denoising_overload_fun=None,
source_pc_names=None,
orchestrator=None,
pair_dump_dir=None,
pair_key="PAIR_0",
uncorrected_grid_right=None,
geoid_path=None,
cloud_id=None,
performance_maps_param=None,
point_cloud_format="laz",
point_cloud_dir=None,
save_output_coordinates=False,
save_output_color=False,
save_output_classification=False,
save_output_filling=False,
save_output_performance_map=False,
save_output_ambiguity=False,
save_output_edges=False,
save_residues=False,
):
"""
Run Triangulation application.
Created left and right CarsDataset filled with xarray.Dataset,
corresponding to 3D point clouds, stored on epipolar geometry grid.
:param sensor_image_left: tiled sensor left image
Dict Must contain keys : "image", "texture", "geomodel",
"no_data", "mask". Paths must be absolutes
:type sensor_image_left: CarsDataset
:param sensor_images_right: tiled sensor right image
Dict Must contain keys : "image", "texture", "geomodel",
"no_data", "mask". Paths must be absolutes
:type sensor_images_right: CarsDataset
:param grid_left: left epipolar grid. Grid dict contains :
- "grid_spacing", "grid_origin",\
"epipolar_size_x", epipolar_size_y", "epipolar_origin_x",\
"epipolar_origin_y","epipolar_spacing_x",\
"epipolar_spacing", "disp_to_alt_ratio", "path"
:type grid_left: dict
:param grid_right: right epipolar grid. Grid dict contains :
- "grid_spacing", "grid_origin",\
"epipolar_size_x", epipolar_size_y", "epipolar_origin_x",\
"epipolar_origin_y","epipolar_spacing_x",\
"epipolar_spacing", "disp_to_alt_ratio", "path"
:type grid_right: dict
:param epipolar_disparity_map: tiled left disparity map or \
sparse matches:
- if CarsDataset is instance of "arrays", CarsDataset contains:
- N x M Delayed tiles \
Each tile will be a future xarray Dataset containing:
- data with keys : "disp", "disp_msk"
- attrs with keys: profile, window, overlaps
- attributes containing:"largest_epipolar_region"\
"opt_epipolar_tile_size"
- if CarsDataset is instance of "points", CarsDataset contains:
- N x M Delayed tiles \
Each tile will be a future pandas DataFrame containing:
- data : (L, 4) shape matches
- attributes containing:"disp_lower_bound","disp_upper_bound",\
"elevation_delta_lower_bound","elevation_delta_upper_bound"
:type epipolar_disparity_map: CarsDataset
:param epipolar_image: tiled epipolar left image
:type epipolar_image: CarsDataset
:param denoising_overload_fun: function to overload dataset
:type denoising_overload_fun: fun
:param source_pc_names: source pc names
:type source_pc_names: list[str]
:param orchestrator: orchestrator used
:param pair_dump_dir: folder used as dump directory for current pair
:type pair_dump_dir: str
:param pair_key: pair key id
:type pair_key: str
:param uncorrected_grid_right: not corrected right epipolar grid
used if self.snap_to_img1
:type uncorrected_grid_right: CarsDataset
:param geoid_path: geoid path
:type geoid_path: str
:param performance_maps_param: parameters used
to generate performance map
:type performance_maps_param: dict or None
:param depth_map_dir: directory to write triangulation output depth
map.
:type depth_map_dir: None or str
:param save_output_coordinates: Save X, Y, Z coords in depth_map_dir
:type save_output_coordinates: bool
:param save_output_color: Save color depth map in depth_map_dir
:type save_output_color: bool
:param save_output_classification: Save classification depth map in
depth_map_dir
:type save_output_classification: bool
:param save_output_filling: Save filling depth map in depth_map_dir
:type save_output_filling: bool
:param save_output_performance_map: Save performance map in
depth_map_dir
:type save_output_performance_map: bool
:param save_output_ambiguity: Save ambiguity in
depth_map_dir
:type save_output_ambiguity: bool
:param save_output_edges: Save edges in
depth_map_dir
:type save_output_edges: bool
:param save_residues: Save intersections residues in depth_map_dir
:type save_residues: bool
:return: point cloud \
The CarsDataset contains:
- N x M Delayed tiles \
Each tile will be a future xarray Dataset containing:
- data : with keys : "x", "y", "z", "corr_msk"\
optional: "texture", "msk", "z_inf", "z_sup"
- attrs with keys: "margins", "epi_full_size", "epsg"
- attributes containing: "disp_lower_bound", "disp_upper_bound", \
"elevation_delta_lower_bound","elevation_delta_upper_bound"
:rtype: Tuple(CarsDataset, CarsDataset)
"""
if not isinstance(sensor_images_right, list):
sensor_images_right = [sensor_images_right]
# Check parameters, to diferentiate sensor or epipolar triangulation
epi_parameters = [
epipolar_disparity_map,
epipolar_image,
grid_left,
grid_right,
uncorrected_grid_right,
]
sensor_parameters = [sensor_matches]
check_epi = all(x is None for x in epi_parameters) or all(
x is not None for x in epi_parameters
)
check_sensor = all(x is None for x in sensor_parameters) or all(
x is not None for x in sensor_parameters
)
is_epipolar = check_epi and sensor_matches is None
if is_epipolar:
matches_cars_dss = [epipolar_disparity_map]
else:
matches_cars_dss = sensor_matches
if not check_epi or not check_sensor:
raise RuntimeError(
"For triangulation on epipolar geometry, all these parameters "
"must be provided: epipolar_disparity_map, epipolar_image,"
"uncorrected_grid_right "
"grid_left and grid_right. \n"
"For triangulation on sensor geometry N views, all "
"these parameters "
"must be provided: sensor matches"
)
# Default orchestrator
if orchestrator is None:
# Create default sequential orchestrator for current application
# be awere, no out_json will be shared between orchestrators
# No files saved
self.orchestrator = ocht.Orchestrator(
orchestrator_conf={"mode": "sequential"}
)
else:
self.orchestrator = orchestrator
if matches_cars_dss[0].dataset_type not in ("arrays", "points"):
raise RuntimeError(
"Triangulation application doesn't support this input "
"data format"
)
if isinstance(point_cloud_format, str):
point_cloud_format = [point_cloud_format]
sensor1 = sensor_image_left[sens_cst.INPUT_IMG]
sensors2 = [sens[sens_cst.INPUT_IMG] for sens in sensor_images_right]
geomodel1 = sensor_image_left[sens_cst.INPUT_GEO_MODEL]
geomodels2 = [
sens[sens_cst.INPUT_GEO_MODEL] for sens in sensor_images_right
]
if is_epipolar:
# Interpolate grid
interpolated_grid_left = RectificationGrid(
grid_left["path"],
interpolator=geometry_plugin.interpolator,
)
interpolated_grid_right = RectificationGrid(
grid_right["path"],
interpolator=geometry_plugin.interpolator,
)
broadcasted_interpolated_grid_left = (
self.orchestrator.cluster.scatter(interpolated_grid_left)
)
broadcasted_interpolated_grid_right = (
self.orchestrator.cluster.scatter(interpolated_grid_right)
)
disp_to_alt_ratio = grid_left["disp_to_alt_ratio"]
# Transform
disp_min_tiling = epipolar_disparity_map.attributes[
"disp_min_tiling"
]
disp_max_tiling = epipolar_disparity_map.attributes[
"disp_max_tiling"
]
# change to N+1 M+1 dimension, fitting to tiling
(
disp_min_tiling,
disp_max_tiling,
) = tiling.transform_disp_range_grid_to_two_layers(
disp_min_tiling, disp_max_tiling
)
(
epipolar_grid_min,
epipolar_grid_max,
) = grid_generation_algo.compute_epipolar_grid_min_max(
geometry_plugin,
tiling.transform_four_layers_to_two_layers_grid(
epipolar_image.tiling_grid
),
sensor1,
sensors2[0],
geomodel1,
geomodels2[0],
grid_left,
grid_right,
epsg,
disp_min_tiling,
disp_max_tiling,
)
largest_epipolar_region = epipolar_image.attributes[
"largest_epipolar_region"
]
color_type = epipolar_image.attributes["image_type"]
opt_epipolar_tile_size = epipolar_image.attributes["tile_width"]
else:
broadcasted_interpolated_grid_left = None
broadcasted_interpolated_grid_right = None
disp_to_alt_ratio = (
None # TODO evolve for performance map generation
)
epipolar_grid_min = None
epipolar_grid_max = None
largest_epipolar_region = None
color_type = None
opt_epipolar_tile_size = None
if source_pc_names is None:
source_pc_names = ["PAIR_0"]
if pair_dump_dir is None:
pair_dump_dir = os.path.join(self.orchestrator.out_dir, "tmp")
# Get local conf left image for this in_json iteration
conf_left_img = sensor_image_left[sens_cst.INPUT_IMG]
# Check left image and raise a warning
# if different left images are used along with snap_to_img1 mode
if self.ref_left_img is None:
self.ref_left_img = conf_left_img
else:
if self.snap_to_img1 and self.ref_left_img != conf_left_img:
logging.warning(
"snap_to_left_image mode is used but inputs "
"have different images as their "
"left image in pair. This may result in "
"increasing registration discrepancies between pairs"
)
# Add log about geoid
if geoid_path is not None:
alt_reference = "geoid"
else:
alt_reference = "ellipsoid"
# Add infos to orchestrator.out_json
updating_dict = {
application_constants.APPLICATION_TAG: {
triangulation_constants.TRIANGULATION_RUN_TAG: {
pair_key: {
triangulation_constants.ALT_REFERENCE_TAG: (
alt_reference
)
},
}
}
}
self.orchestrator.update_out_info(updating_dict)
if self.snap_to_img1:
grid_right = uncorrected_grid_right
if grid_right is None:
logging.error(
"Uncorrected grid was not given in order "
"to snap it to img1"
)
# Compute disp_min and disp_max location for epipolar grid
# update attributes for corresponding tiles in cloud fusion
# TODO remove with refactoring
pc_attributes = {
"used_epsg_for_terrain_grid": epsg,
"epipolar_grid_min": epipolar_grid_min,
"epipolar_grid_max": epipolar_grid_max,
"largest_epipolar_region": largest_epipolar_region,
"source_pc_names": source_pc_names,
"source_pc_name": pair_key,
"color_type": color_type,
"opt_epipolar_tile_size": opt_epipolar_tile_size,
}
if geoid_path:
pc_attributes["geoid"] = (geoid_path,)
if matches_cars_dss[0].dataset_type != "arrays":
raise RuntimeError(
"Triangulation application doesn't support this input "
"data format"
)
# Check performance_maps_parameters
performance_maps_to_generate = None
if performance_maps_param is not None:
if "performance_map_method" not in performance_maps_param:
raise RuntimeError("No performance_map_method specified")
performance_maps_to_generate = performance_maps_param[
"performance_map_method"
]
if "perf_ambiguity_threshold" not in performance_maps_param:
raise RuntimeError("No perf_ambiguity_threshold specified")
# Create CarsDataset
# Epipolar_point_cloud or sensor depth map
sens_or_epi_depth_map = cars_dataset.CarsDataset(
matches_cars_dss[0].dataset_type,
name="triangulation_" + pair_key,
)
sens_or_epi_depth_map.create_empty_copy(matches_cars_dss[0])
# Update attributes to get epipolar info
sens_or_epi_depth_map.attributes.update(pc_attributes)
# Save objects
# Save as depth map
self.save_triangulation_output(
sens_or_epi_depth_map,
sensor_image_left,
point_cloud_dir,
pair_dump_dir if self.save_intermediate_data else None,
performance_maps_to_generate,
save_output_coordinates,
save_output_color,
save_output_classification,
save_output_filling,
save_output_performance_map,
save_output_ambiguity,
save_output_edges,
save_residues,
)
self.fill_index(
save_output_coordinates,
save_output_color,
save_output_classification,
save_output_filling,
save_output_performance_map,
save_output_ambiguity,
save_output_edges,
save_residues,
pair_key,
)
# Save as point cloud
point_cloud = cars_dataset.CarsDataset(
"points",
name="triangulation_flatten_" + pair_key,
)
point_cloud.create_empty_copy(sens_or_epi_depth_map)
point_cloud.attributes = sens_or_epi_depth_map.attributes
csv_pc_dir_name = None
laz_pc_dir_name = None
if "laz" in point_cloud_format:
csv_pc_dir_name, laz_pc_dir_name = (
self.create_point_cloud_directories(
pair_dump_dir, point_cloud_dir, point_cloud
)
)
# Get saving infos in order to save tiles when they are computed
[saving_info_point_cloud] = self.orchestrator.get_saving_infos(
[sens_or_epi_depth_map]
)
saving_info_flatten = None
if self.save_intermediate_data or point_cloud_dir is not None:
[saving_info_flatten] = self.orchestrator.get_saving_infos(
[point_cloud]
)
# Generate Point clouds
# initialize empty index file for point cloud product if official
# product is requested
pc_index = None
if point_cloud_dir:
pc_index = {}
for col in range(sens_or_epi_depth_map.shape[1]):
for row in range(sens_or_epi_depth_map.shape[0]):
if matches_cars_dss[0][row, col] is not None:
# update saving infos for potential replacement
full_saving_info_epipolar = ocht.update_saving_infos(
saving_info_point_cloud, row=row, col=col
)
full_saving_info_flatten = None
if saving_info_flatten is not None:
full_saving_info_flatten = ocht.update_saving_infos(
saving_info_flatten, row=row, col=col
)
csv_pc_file_name, laz_pc_file_name = (
triangulation_wrap.generate_point_cloud_file_names(
csv_pc_dir_name,
laz_pc_dir_name,
row,
col,
pc_index,
pair_key,
)
)
if is_epipolar:
wrapper_kwargs = {
"grid1": broadcasted_interpolated_grid_left,
"grid2": broadcasted_interpolated_grid_right,
"disparity_object": epipolar_disparity_map[
row, col
],
"disp_to_alt_ratio": disp_to_alt_ratio,
}
else:
sensor_matches = [
cars_ds[row, col] for cars_ds in matches_cars_dss
]
wrapper_kwargs = {
"sensor_matches": sensor_matches,
}
# Compute points
(
sens_or_epi_depth_map[row][col],
point_cloud[row][col],
) = self.orchestrator.cluster.create_task(
triangulation_wrapper, nout=2
)(
is_epipolar,
sensor1,
sensors2,
geomodel1,
geomodels2,
geometry_plugin,
epsg,
z_inf_sup_bh_approximation=(
self.z_inf_sup_bh_approximation
),
geoid_path=geoid_path,
denoising_overload_fun=denoising_overload_fun,
cloud_id=cloud_id,
performance_maps_to_generate=(
performance_maps_to_generate
),
performance_maps_parameters=performance_maps_param,
point_cloud_csv_file_name=csv_pc_file_name,
point_cloud_laz_file_name=laz_pc_file_name,
saving_info_epipolar=full_saving_info_epipolar,
saving_info_flatten=full_saving_info_flatten,
**wrapper_kwargs,
)
# update point cloud index
if point_cloud_dir:
self.orchestrator.update_index(pc_index)
return sens_or_epi_depth_map
[docs]
def triangulation_wrapper( # noqa: C901 function is too complex
# pylint: disable=too-many-positional-arguments
is_epipolar,
sensor1,
sensors2,
geomodel1,
geomodels2,
geometry_plugin,
epsg,
grid1=None,
grid2=None,
disparity_object=None,
sensor_matches=None,
z_inf_sup_bh_approximation=True,
disp_to_alt_ratio=None,
geoid_path=None,
denoising_overload_fun=None,
cloud_id=None,
performance_maps_to_generate=None,
performance_maps_parameters=None,
point_cloud_csv_file_name=None,
point_cloud_laz_file_name=None,
saving_info_epipolar=None,
saving_info_flatten=None,
) -> Dict[str, Tuple[xr.Dataset, xr.Dataset]]:
"""
Compute point clouds from image objects and disparity objects.
:param disparity_object: Left disparity map dataset with :
- cst_disp.MAP
- cst_disp.VALID
- cst.EPI_TEXTURE
:type disparity_object: xr.Dataset
:param sensor1: path to left sensor image
:type sensor1: str
:param sensor2: path to right sensor image
:type sensor2: str
:param geomodel1: path and attributes for left geomodel
:type geomodel1: dict
:param geomodel2: path and attributes for right geomodel
:type geomodel2: dict
:param grid1: dataset of the reference image grid file
:type grid1: CarsDataset
:param grid2: dataset of the secondary image grid file
:type grid2: CarsDataset
:param geometry_plugin: geometry plugin to use
:type geometry_plugin: AbstractGeometry
:param geoid_path: Geoid used for altimetric reference. Defaults to None
for using ellipsoid as altimetric reference.
:type geoid_path: str
:param performance_maps_to_generate: None or list containing
"risk" or "intervals"
:param performance_maps_parameters: parameters used to
generate performance map
:type performance_maps_parameters: dict or None
:param denoising_overload_fun: function to overload dataset
:type denoising_overload_fun: fun
:return: Left disparity object
Returned object is composed of :
- dataset with :
- cst.X
- cst.Y
- cst.Z
- cst.EPI_TEXTURE
- cst.Z_INF (optional)
- cst.Z_SUP (optional)
"""
if is_epipolar:
# Get disparity maps
matches_dataset = disparity_object
# Triangulate
if not isinstance(matches_dataset, xr.Dataset):
logging.error("Disp ref not xarray Dataset")
raise TypeError("Disp ref not xarray Dataset")
# Triangulate epipolar dense disparities
points = triangulation_algo.triangulate(
geometry_plugin,
sensor1,
sensors2,
geomodel1,
geomodels2,
grid1=grid1,
grid2=grid2,
disp_ref=matches_dataset,
)
else:
matches_dataset = sensor_matches[0]
# Triangulate sensor dense matches
points = triangulation_algo.triangulate(
geometry_plugin,
sensor1,
sensors2,
geomodel1,
geomodels2,
sensor_matches=sensor_matches,
)
if performance_maps_to_generate is not None and is_epipolar:
for perf_method in performance_maps_to_generate:
# Generate keys to use
if perf_method == "risk":
# From pandora loader overload
disp_keys = [
"confidence_from_disp_inf_from_risk.cars_2",
"confidence_from_disp_sup_from_risk.cars_2",
]
cars_inf_key = cst.POINT_CLOUD_LAYER_SUP_FROM_RISK
cars_sup_key = cst.POINT_CLOUD_LAYER_INF_FROM_RISK
cars_perf_key = cst.POINT_CLOUD_PERFORMANCE_MAP_FROM_RISK
use_ambiguity = True
else:
# From pandora loader overload
disp_keys = [
"confidence_from_interval_bounds_inf.cars_3",
"confidence_from_interval_bounds_sup.cars_3",
]
cars_inf_key = cst.POINT_CLOUD_LAYER_SUP_FROM_INTERVALS
cars_sup_key = cst.POINT_CLOUD_LAYER_INF_FROM_INTERVALS
cars_perf_key = cst.POINT_CLOUD_PERFORMANCE_MAP_FROM_INTERVALS
use_ambiguity = False
ambiguity_map = None
perf_ambiguity_threshold = None
if use_ambiguity:
ambiguity_map = matches_dataset["ambiguity"]
perf_ambiguity_threshold = performance_maps_parameters[
"perf_ambiguity_threshold"
]
if z_inf_sup_bh_approximation:
z_points = points[cst.STEREO_REF][cst.Z]
disp_min = matches_dataset[disp_keys[0]]
disp_max = matches_dataset[disp_keys[1]]
points_inf = (
-disp_to_alt_ratio * (disp_min - matches_dataset["disp"])
+ z_points.values
)
points_sup = (
-disp_to_alt_ratio * (disp_max - matches_dataset["disp"])
+ z_points.values
)
points[cst.STEREO_REF][cars_inf_key] = (
("row", "col"),
points_inf.data,
)
points[cst.STEREO_REF][cars_sup_key] = (
("row", "col"),
points_sup.data,
)
else:
# triangulate disp inf and supp
points_inf = triangulation_algo.triangulate(
geometry_plugin,
sensor1,
sensors2,
geomodel1,
geomodels2,
grid1,
grid2,
matches_dataset,
disp_key=disp_keys[0],
)
points_sup = triangulation_algo.triangulate(
geometry_plugin,
sensor1,
sensors2,
geomodel1,
geomodels2,
grid1,
grid2,
matches_dataset,
disp_key=disp_keys[1],
)
points[cst.STEREO_REF][cars_inf_key] = points_inf[
cst.STEREO_REF
][cst.Z]
points[cst.STEREO_REF][cars_sup_key] = points_sup[
cst.STEREO_REF
][cst.Z]
# Generate performance map
points[cst.STEREO_REF][cars_perf_key] = (
triangulation_wrap.compute_performance_map(
points[cst.STEREO_REF][cst.Z],
points[cst.STEREO_REF][cars_inf_key],
points[cst.STEREO_REF][cars_sup_key],
ambiguity_map=ambiguity_map,
perf_ambiguity_threshold=perf_ambiguity_threshold,
)
)
if geoid_path is not None: # if user pass a geoid, use it as alt reference
for key, point in points.items():
points[key] = triangulation_wrap.geoid_offset(point, geoid_path)
# Fill datasets
pc_dataset = points[cst.STEREO_REF]
pc_dataset.attrs["cloud_id"] = cloud_id
# Overload dataset with denoising fun
if denoising_overload_fun is not None:
if isinstance(pc_dataset, xr.Dataset):
denoising_overload_fun(
pc_dataset,
sensor1,
sensors2,
geomodel1,
geomodels2,
grid1,
grid2,
geometry_plugin,
matches_dataset,
)
else:
raise RuntimeError("wrong pc type for denoising func")
attributes = {
cst.CROPPED_DISPARITY_RANGE: (
ocht.get_disparity_range_cropped(matches_dataset)
)
}
cars_dataset.fill_dataset(
pc_dataset,
saving_info=saving_info_epipolar,
window=cars_dataset.get_window_dataset(matches_dataset),
profile=cars_dataset.get_profile_rasterio(matches_dataset),
attributes=attributes,
overlaps=cars_dataset.get_overlaps_dataset(matches_dataset),
)
# Flatten point cloud to save it as LAZ
flatten_pc_dataset = None
if point_cloud_csv_file_name or point_cloud_laz_file_name:
# Convert epipolar array into point cloud
flatten_pc_dataset, cloud_epsg = (
pc_transform.depth_map_dataset_to_dataframe(pc_dataset, epsg)
)
# Convert to wanted epsg
if epsg is not None and cloud_epsg != epsg:
projection.point_cloud_conversion_dataframe(
flatten_pc_dataset, cloud_epsg, epsg
)
cloud_epsg = epsg
# Fill attributes for LAZ saving
color_type = pc_transform.get_color_type([pc_dataset])
attributes = {
"epsg": cloud_epsg,
"color_type": color_type,
}
cars_dataset.fill_dataframe(
flatten_pc_dataset,
saving_info=saving_info_flatten,
attributes=attributes,
)
# Save point cloud in worker
if point_cloud_csv_file_name:
cars_dataset.run_save_points(
flatten_pc_dataset,
point_cloud_csv_file_name,
overwrite=True,
point_cloud_format="csv",
overwrite_file_name=False,
)
if point_cloud_laz_file_name:
cars_dataset.run_save_points(
flatten_pc_dataset,
point_cloud_laz_file_name,
overwrite=True,
point_cloud_format="laz",
overwrite_file_name=False,
)
return pc_dataset, flatten_pc_dataset
# pylint: disable=too-many-positional-arguments
[docs]
def triangulation_wrapper_matches(
matches,
sensor1,
sensor2,
geomodel1,
geomodel2,
grid1,
grid2,
geometry_plugin,
full_saving_info_matches,
epsg,
):
"""
:param matches: matches
:type matches: cars_dataset
:param sensor1: path to left sensor image
:type sensor1: str
:param sensor2: path to right sensor image
:type sensor2: str
:param interpolated_grid1: rectification grid left
:type interpolated_grid1: shareloc.rectificationGrid
:param interpolated_grid2: rectification grid right
:type interpolated_grid2: shareloc.rectificationGrid
:param geometry_plugin: geometry plugin to use
:type geometry_plugin: AbstractGeometry
:param epsg: ground epsg
:type epsg: int
"""
# Ignore this tile if there are not match
if matches is None or matches.empty:
return None
epipolar_matches = matches.to_numpy()
sensor_matches = geometry_plugin.matches_to_sensor_coords(
grid1,
grid2,
epipolar_matches,
cst.MATCHES_MODE,
)
sensor_matches = np.concatenate(sensor_matches, axis=1)
matches = np.concatenate(
[
epipolar_matches,
sensor_matches,
],
axis=1,
)
# Triangulate matches
triangulated_matches = triangulation_algo.triangulate_sparse_matches(
sensor1,
sensor2,
geomodel1,
geomodel2,
grid1,
grid2,
matches,
geometry_plugin,
epsg,
)
cars_dataset.fill_dataframe(
triangulated_matches,
saving_info=full_saving_info_matches,
attributes={"epsg": triangulated_matches.attrs["epsg"]},
)
return triangulated_matches