Merge branch 'reentrance-dsm-phases' into 'master'

feat: add use of phased dsm

Closes #967

See merge request 3d/cars-park/cars!795

cars/applications/rasterization/rasterization_tools.py

@@ -65,7 +65,6 @@ def compute_xy_starts_and_sizes(
 
     # Clamp to a regular grid
     x_start = np.floor(xmin / resolution) * resolution
-    x_size = int(1 + np.floor((xmax - x_start) / resolution))
 
     # Derive ystart
     ymin = np.nanmin(cloud[cst.Y].values)
@@ -74,6 +73,8 @@ def compute_xy_starts_and_sizes(
 
     # Clamp to a regular grid
     y_start = np.ceil(ymax / resolution) * resolution
+
+    x_size = int(1 + np.floor((xmax - x_start) / resolution))
     y_size = int(1 + np.floor((y_start - ymin) / resolution))
 
     return x_start, y_start, x_size, y_size
@@ -202,6 +203,25 @@ def substring_in_list(src_list, substring):
     return len(res) > 0
 
 
+def phased_dsm(start: float, phase: float, resolution: float):
+    """
+    Phased the dsm
+
+    :param start: start of the roi
+    :param phase: the point for phasing
+    :param resolution: resolution of the dsm
+    """
+
+    div = np.abs(start - phase) / resolution
+
+    if phase > start:
+        start = phase - resolution * np.floor(div)
+    else:
+        start = resolution * np.floor(div) + phase
+
+    return start
+
+
 def find_indexes_in_point_cloud(
     cloud: pandas.DataFrame, tag: str, list_computed_layers: List[str] = None
 ) -> List[str]:
@@ -845,19 +865,21 @@ def update_data(
 
     :return: updated current data
     """
+
     new_data = current_data
+    data = old_data
     if old_data is not None:
         old_data = np.squeeze(old_data)
         old_weights = np.squeeze(old_weights)
         shape = old_data.shape
-        if len(old_data.shape) == 3:
+        if len(data.shape) == 3 and data.shape[0] > 1:
             old_weights = np.repeat(
                 np.expand_dims(old_weights, axis=0), old_data.shape[0], axis=0
             )
 
         current_data = np.squeeze(current_data)
         weights = np.squeeze(weights)
-        if len(current_data.shape) == 3:
+        if len(new_data.shape) == 3 and new_data.shape[0] > 1:
             weights = np.repeat(
                 np.expand_dims(weights, axis=0), current_data.shape[0], axis=0
             )
@@ -867,7 +889,9 @@ def update_data(
         old_valid = old_weights != 0
 
         both_valid = np.logical_and(current_valid, old_valid)
+
         total_weights = np.zeros(shape)
+
         total_weights[both_valid] = (
             weights[both_valid] + old_weights[both_valid]
         )

cars/applications/rasterization/simple_gaussian.py

@@ -217,6 +217,7 @@ def run(  # noqa: C901 function is too complex
         filling_file_name=None,
         color_dtype=None,
         dump_dir=None,
+        phasing=None,
     ):
         """
         Run PointCloudRasterisation application.
@@ -271,6 +272,8 @@ def run(  # noqa: C901 function is too complex
         :type color_dtype: str (numpy type)
         :param dump_dir: directory used for outputs with no associated filename
         :type dump_dir: str
+        :param phasing: if activated, we phase the dsm on this point
+        :type phasing: dict
 
         :return: raster DSM. CarsDataset contains:
 
@@ -360,6 +363,21 @@ def run(  # noqa: C901 function is too complex
 
         terrain_raster.generate_none_tiles()
 
+        if phasing is not None:
+            res = resolution
+            x_phase = phasing["point"][0]
+            y_phase = phasing["point"][1]
+
+            for index, value in enumerate(bounds):
+                if index in (0, 2):
+                    bounds[index] = rasterization_step.phased_dsm(
+                        value, x_phase, res
+                    )
+                else:
+                    bounds[index] = rasterization_step.phased_dsm(
+                        value, y_phase, res
+                    )
+
         # Derive output image files parameters to pass to rasterio
         xsize, ysize = tiling.roi_to_start_and_size(bounds, resolution)[2:]
         logging.info("DSM output image size: {}x{} pixels".format(xsize, ysize))
@@ -1000,6 +1018,7 @@ def rasterization_wrapper(
     xstart, ystart, xsize, ysize = tiling.roi_to_start_and_size(
         terrain_region, resolution
     )
+
     if window is None:
         transform = rio.Affine(*profile["transform"][0:6])
         row_pix_pos, col_pix_pos = rio.transform.AffineTransformer(

cars/core/constants.py

@@ -130,3 +130,26 @@
 INDEX_DEPTH_MAP_PERFORMANCE_MAP = "performance_map"
 INDEX_DEPTH_MAP_FILLING = "filling"
 INDEX_DEPTH_MAP_EPSG = "epsg"
+
+# dsms inputs index
+DSM_CLASSIF = "classification"
+DSM_ALT = "dsm"
+DSM_ALT_INF = "dsm_inf"
+DSM_ALT_SUP = "dsm_sup"
+DSM_WEIGHTS_SUM = "weights"
+DSM_MSK = "mask"
+DSM_NB_PTS = "dsm_n_pts"
+DSM_NB_PTS_IN_CELL = "dsm_pts_in_cell"
+DSM_MEAN = "dsm_mean"
+DSM_STD_DEV = "dsm_std"
+DSM_INF_MEAN = "dsm_inf_mean"
+DSM_INF_STD = "dsm_inf_std"
+DSM_SUP_MEAN = "dsm_sup_mean"
+DSM_SUP_STD = "dsm_sup_std"
+DSM_CONFIDENCE_AMBIGUITY = "confidence_from_ambiguity"
+DSM_CONFIDENCE_RISK_MIN = "confidence_from_risk_min"
+DSM_CONFIDENCE_RISK_MAX = "confidence_from_risk_max"
+DSM_PERFORMANCE_MAP = "performance_map"
+DSM_SOURCE_PC = "source_pc"
+DSM_FILLING = "filling"
+DSM_COLOR = "color"

cars/core/inputs.py

@@ -182,6 +182,28 @@ def rasterio_get_size(raster_file: str) -> Tuple[int, int]:
         return (descriptor.width, descriptor.height)
 
 
+def rasterio_get_nodata(raster_file: str) -> Tuple[int, int]:
+    """
+    Get the no data value
+
+    :param raster_file: Image file
+    :return: the no data value
+    """
+    with rio.open(raster_file, "r") as descriptor:
+        return descriptor.nodata
+
+
+def rasterio_get_dtype(raster_file: str) -> Tuple[int, int]:
+    """
+    Get the dtype of an image (file)
+
+    :param raster_file: Image file
+    :return: The dtype
+    """
+    with rio.open(raster_file, "r") as descriptor:
+        return descriptor.dtypes[0]
+
+
 def rasterio_get_pixel_points(raster_file: str, terrain_points) -> list:
     """
     Get pixel point coordinates of terrain points