Merge branch '797-cars-up-down-sampled' into 'master'

Resolve "CARS avec données  up/down sampled" + grille rectification loc

Closes #818 and #797

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

cars/applications/dense_matching/census_mccnn_sgm.py

@@ -704,7 +704,7 @@ def generate_disparity_grids(  # noqa: C901
                 dem_max, lon_mean, lat_mean, points_cloud_conversion
             )
 
-            # sensors positions as index
+            # sensors physical positions
             (
                 ind_cols_sensor,
                 ind_rows_sensor,

cars/applications/grid_generation/grid_correction.py

@@ -218,23 +218,21 @@ def estimate_right_grid_correction(
     spacing = grid_right.attributes["grid_spacing"]
 
     # Form 3D array with grid positions
-    x_values_1d = np.linspace(
+    x_values_1d = np.arange(
         origin[0],
         origin[0] + right_grid.shape[0] * spacing[0],
-        right_grid.shape[0],
+        spacing[0],
     )
-    y_values_1d = np.linspace(
+    y_values_1d = np.arange(
         origin[1],
         origin[1] + right_grid.shape[1] * spacing[1],
-        right_grid.shape[1],
+        spacing[1],
     )
     x_values_2d, y_values_2d = np.meshgrid(y_values_1d, x_values_1d)
 
     # Compute corresponding point in sensor geometry (grid encodes (x_sensor -
     # x_epi,y_sensor - y__epi)
     source_points = right_grid
-    source_points[:, :, 0] += x_values_2d
-    source_points[:, :, 1] += y_values_2d
 
     # Extract matches for convenience
     matches_y1 = matches[:, 1]

cars/applications/resampling/resampling_tools.py

@@ -23,22 +23,21 @@
 contains functions used for epipolar resampling
 """
 
-import logging
-
 # Standard imports
 import math
 
 # Third party imports
 import numpy as np
 import rasterio as rio
 import resample as cresample
-from rasterio.windows import Window, bounds, from_bounds
+from rasterio.windows import Window, bounds
 
 from cars.conf import mask_cst as msk_cst
 
 # CARS imports
 from cars.core import constants as cst
 from cars.core import datasets, inputs, tiling
+from cars.core.geometry import abstract_geometry
 from cars.data_structures import cars_dataset
 
 
@@ -261,6 +260,7 @@ def resample_image(
     # Localize blocks of the tile to resample
     if step is None:
         step = region[2] - region[0]
+
     xmin_of_blocks = np.arange(region[0], region[2], step)
     xmax_of_blocks = np.append(
         np.arange(region[0] + step, region[2], step), region[2]
@@ -300,74 +300,48 @@ def resample_image(
             grid_as_array = grid_as_array.astype(np.float32)
             grid_as_array = grid_as_array.astype(np.float64)
 
-            # deformation to localization
-            grid_as_array[0, ...] += np.arange(
-                oversampling * grid_region[0],
-                oversampling * (grid_region[2] + 1),
-                step=oversampling,
-            )
-            grid_as_array[1, ...] += np.arange(
-                oversampling * grid_region[1],
-                oversampling * (grid_region[3] + 1),
-                step=oversampling,
-            ).T[..., np.newaxis]
-
             # get needed source bounding box
             left = math.floor(np.amin(grid_as_array[0, ...]))
             right = math.ceil(np.amax(grid_as_array[0, ...]))
             top = math.floor(np.amin(grid_as_array[1, ...]))
             bottom = math.ceil(np.amax(grid_as_array[1, ...]))
 
+            # transform xmin and xmax positions to index
+            (top, bottom, left, right) = (
+                abstract_geometry.min_max_to_index_min_max(
+                    left, right, top, bottom, img_reader.transform
+                )
+            )
+
             # filter margin for bicubic = 2
             filter_margin = 2
             top -= filter_margin
             bottom += filter_margin
             left -= filter_margin
             right += filter_margin
 
-            # extract src according to grid values
-            transform = img_reader.transform
-            res_x = int(transform[0] / abs(transform[0]))
-            res_y = int(transform[4] / abs(transform[4]))
-
-            (full_left, full_bottom, full_right, full_top) = img_reader.bounds
-
-            left, right, top, bottom = (
-                res_x * left,
-                res_x * right,
-                res_y * top,
-                res_y * bottom,
-            )
+            left, right = list(np.clip([left, right], 0, img_reader.shape[0]))
+            top, bottom = list(np.clip([top, bottom], 0, img_reader.shape[1]))
 
-            full_bounds_window = from_bounds(
-                full_left, full_bottom, full_right, full_top, transform
-            )
-            img_window = from_bounds(left, bottom, right, top, transform)
+            img_window = Window.from_slices([left, right], [top, bottom])
 
-            # Crop window to be in image
             in_sensor = True
-            try:
-                img_window = img_window.intersection(full_bounds_window)
-            except rio.errors.WindowError:
-                # Window not in sensor image
-                logging.debug("Window not in sensor image")
+            if right - left == 0 or bottom - top == 0:
                 in_sensor = False
 
             # round window
             img_window = img_window.round_offsets()
             img_window = img_window.round_lengths()
 
             # Compute offset
-            tile_bounds = bounds(img_window, transform)
-            tile_bounds_with_res = [
-                res_x * tile_bounds[0],
-                res_y * tile_bounds[1],
-                res_x * tile_bounds[2],
-                res_y * tile_bounds[3],
-            ]
+            transform = img_reader.transform
 
-            x_offset = min(tile_bounds_with_res[0], tile_bounds_with_res[2])
-            y_offset = min(tile_bounds_with_res[1], tile_bounds_with_res[3])
+            res_x = float(abs(transform[0]))
+            res_y = float(abs(transform[4]))
+            tile_bounds = list(bounds(img_window, transform))
+
+            x_offset = min(tile_bounds[0], tile_bounds[2])
+            y_offset = min(tile_bounds[1], tile_bounds[3])
 
             if in_sensor:
                 # Get sensor data
@@ -377,6 +351,10 @@ def resample_image(
                 grid_as_array[0, ...] -= x_offset
                 grid_as_array[1, ...] -= y_offset
 
+                # apply input resolution
+                grid_as_array[0, ...] /= res_x
+                grid_as_array[1, ...] /= res_y
+
                 block_resamp = cresample.grid(
                     img_as_array,
                     grid_as_array,

cars/core/geometry/abstract_geometry.py

@@ -32,6 +32,7 @@
 from scipy import interpolate
 from scipy.interpolate import LinearNDInterpolator
 from shapely.geometry import Polygon
+from shareloc import proj_utils
 
 from cars.core import constants as cst
 from cars.core import constants_disparity as cst_disp
@@ -363,11 +364,8 @@ def sensor_position_from_grid(
 
         # create regular grid points positions
         points = (cols, rows)
-        grid_row, grid_col = np.mgrid[
-            ori_row:last_row:step_row, ori_col:last_col:step_col
-        ]
-        sensor_row_positions = row_dep + grid_row
-        sensor_col_positions = col_dep + grid_col
+        sensor_row_positions = row_dep
+        sensor_col_positions = col_dep
 
         # interpolate sensor positions
         interp_row = interpolate.interpn(
@@ -644,3 +642,75 @@ def image_envelope(self, sensor, geomodel, shp=None):
             outputs.write_vector([poly_bb], shp, 4326, driver="ESRI Shapefile")
 
         return u_l, u_r, l_l, l_r
+
+
+def min_max_to_physical_min_max(xmin, xmax, ymin, ymax, transform):
+    """
+    Transform min max index to position min max
+
+    :param xmin: xmin
+    :type xmin: int
+    :param xmax: xmax
+    :type xmax: int
+    :param ymin: ymin
+    :type ymin: int
+    :param ymax: ymax
+    :type ymax: int
+    :param transform: transform
+    :type transform: Affine
+
+    :return: xmin, xmax, ymin, ymax
+    :rtype: list(int)
+    """
+
+    cols_ind = np.array([xmin, xmin, xmax, xmax])
+    rows_ind = np.array([ymin, ymax, ymin, ymax])
+
+    rows_pos, cols_pos = proj_utils.transform_index_to_physical_point(
+        transform,
+        rows_ind,
+        cols_ind,
+    )
+
+    return (
+        np.min(cols_pos),
+        np.max(cols_pos),
+        np.min(rows_pos),
+        np.max(rows_pos),
+    )
+
+
+def min_max_to_index_min_max(xmin, xmax, ymin, ymax, transform):
+    """
+    Transform min max position to index min max
+
+    :param xmin: xmin
+    :type xmin: int
+    :param xmax: xmax
+    :type xmax: int
+    :param ymin: ymin
+    :type ymin: int
+    :param ymax: ymax
+    :type ymax: int
+    :param transform: transform
+    :type transform: Affine
+
+    :return: xmin, xmax, ymin, ymax
+    :rtype: list(int)
+    """
+
+    cols_ind = np.array([xmin, xmin, xmax, xmax])
+    rows_ind = np.array([ymin, ymax, ymin, ymax])
+
+    rows_pos, cols_pos = proj_utils.transform_physical_point_to_index(
+        ~transform,
+        rows_ind,
+        cols_ind,
+    )
+
+    return (
+        np.min(cols_pos),
+        np.max(cols_pos),
+        np.min(rows_pos),
+        np.max(rows_pos),
+    )

cars/core/geometry/shareloc_geometry.py

@@ -25,7 +25,7 @@
 import logging
 from typing import List, Tuple, Union
 
-import bindings_cpp
+import bindings_cpp  # pylint: disable=E0401
 import numpy as np
 import rasterio as rio
 import shareloc.geofunctions.rectification as rectif
@@ -373,7 +373,6 @@ def generate_epipolar_grids(
             grid2,
             [epipolar_size_y, epipolar_size_x],
             alt_to_disp_ratio,
-            _,
         ) = rectif.compute_stereorectification_epipolar_grids(
             image1,
             shareloc_model1,
@@ -388,21 +387,20 @@ def generate_epipolar_grids(
         grid1 = grid1[:, :, 0:2][:, :, ::-1]
         grid2 = grid2[:, :, 0:2][:, :, ::-1]
 
-        # compute associated characteristics
+        epipolar_size_x = int(np.floor(epipolar_size_x))
+        epipolar_size_y = int(np.floor(epipolar_size_y))
+
+        origin = [0.0, 0.0]
+        spacing = [float(epipolar_step), float(epipolar_step)]
+
+        # alt_to_disp_ratio does not consider image resolution
         with rio.open(sensor1, "r") as rio_dst:
             pixel_size_x, pixel_size_y = (
                 rio_dst.transform[0],
                 rio_dst.transform[4],
             )
-
         mean_size = (abs(pixel_size_x) + abs(pixel_size_y)) / 2
-        epipolar_size_x = int(np.floor(epipolar_size_x * mean_size))
-        epipolar_size_y = int(np.floor(epipolar_size_y * mean_size))
-
-        origin = [0.0, 0.0]
-        spacing = [float(epipolar_step), float(epipolar_step)]
-
-        disp_to_alt_ratio = 1 / alt_to_disp_ratio
+        disp_to_alt_ratio = mean_size / alt_to_disp_ratio
 
         return (
             grid1,

cars/pipelines/sensor_to_dense_dsm/sensor_to_dense_dsm_pipeline.py

@@ -731,7 +731,6 @@ def run(self):  # noqa C901
                         pair_key=pair_key,
                         orchestrator=cars_orchestrator,
                     )
-
                     # Correct grid right
                     pairs[pair_key]["corrected_grid_right"] = (
                         grid_correction.correct_grid(

docs/source/howto.rst

@@ -207,6 +207,46 @@ Any OTB application can be ran in docker
     docker run  --entrypoint=/bin/bash  cnes/cars otbcli_BandMath -help
 
 
+.. _resample_image:
+
+Resample an image
+========================
+
+If you want to upscale or downscale the resolution of you input data, use rasterio:
+
+.. code-block:: python
+
+    import rasterio
+    from rasterio.enums import Resampling
+    # Get data
+    upscale_factor = 0.5
+    with rasterio.open("example.tif") as dataset:
+        # resample data to target shape
+        data = dataset.read(
+            out_shape=(
+                dataset.count,
+                int(dataset.height * upscale_factor),
+                int(dataset.width * upscale_factor)
+            ),
+            resampling=Resampling.bilinear
+        )
+        # scale image transform
+        transform = dataset.transform * dataset.transform.scale(
+            (dataset.width / data.shape[-1]),
+            (dataset.height / data.shape[-2])
+        )
+        profile = dataset.profile
+        # Save data
+        profile.update(
+            width=data.shape[-1],
+            height=data.shape[-2],
+            transform=transform
+        )
+        with rasterio.open('resampled_example.tif', 'w', **profile) as dst:
+            dst.write(data)
+
+
+
 Use CARS with Pleiades images ...
 ========================================
 

setup.cfg

@@ -81,7 +81,7 @@ install_requires =
     cars-rasterize==0.2.*
     cars-resample==0.1.*
     vlsift==0.1.*
-    shareloc==0.2.2
+    shareloc==0.2.3
 
 package_dir =
     . = cars