diff --git a/s2p/__init__.py b/s2p/__init__.py index 4467f16a..2fc807f0 100644 --- a/s2p/__init__.py +++ b/s2p/__init__.py @@ -129,7 +129,7 @@ def pointing_correction(tile, i, cfg): A, m = pointing_accuracy.compute_correction( img1, img2, rpc1, rpc2, x, y, w, h, method, cfg['sift_match_thresh'], cfg['max_pointing_error'], cfg['matching_method'], - cfg['min_value'], cfg['max_value'], cfg['confidence_threshold'] + cfg['min_value'], cfg['max_value'], cfg['confidence_threshold'], cfg ) if A is not None: # A is the correction matrix @@ -142,7 +142,7 @@ def pointing_correction(tile, i, cfg): visualisation.plot_matches(img1, img2, rpc1, rpc2, m, os.path.join(out_dir, 'sift_matches_pointing.png'), - x, y, w, h) + x, y, w, h, cfg) def global_pointing_correction(cfg, tiles): @@ -212,10 +212,14 @@ def rectification_pair(tile, i, cfg): H1, H2, disp_min, disp_max = rectification.rectify_pair(img1, img2, rpc1, rpc2, x, y, w, h, - rect1, rect2, A, m, + rect1, rect2, + cfg, + A=A, + sift_matches=m, method=cfg['rectification_method'], hmargin=cfg['horizontal_margin'], - vmargin=cfg['vertical_margin']) + vmargin=cfg['vertical_margin'], + ) np.savetxt(os.path.join(out_dir, 'H_ref.txt'), H1, fmt='%12.6f') np.savetxt(os.path.join(out_dir, 'H_sec.txt'), H2, fmt='%12.6f') @@ -246,7 +250,7 @@ def stereo_matching(tile, i, cfg): disp_min, disp_max = np.loadtxt(os.path.join(out_dir, 'disp_min_max.txt')) block_matching.compute_disparity_map(rect1, rect2, disp, mask, - cfg['matching_algorithm'], disp_min, + cfg['matching_algorithm'], cfg, disp_min, disp_max, timeout=cfg['mgm_timeout'], max_disp_range=cfg['max_disp_range']) @@ -335,7 +339,7 @@ def disparity_to_ply(tile, cfg): with rasterio.open(os.path.join(out_dir, 'pair_1', 'rectified_ref.tif')) as f: ww, hh = f.width, f.height - colors = common.tmpfile(".tif") + colors = common.tmpfile(cfg['temporary_dir'], ".tif") success_state = common.image_apply_homography(colors, cfg['images'][0]['clr'], np.loadtxt(H_ref), ww, hh) if success_state is False: @@ -450,7 +454,7 @@ def heights_fusion(cfg, tile): # merge the height maps (applying mean offset to register) fusion.merge_n(os.path.join(tile_dir, 'height_map.tif'), height_maps, global_mean_heights, averaging=cfg['fusion_operator'], - threshold=cfg['fusion_thresh']) + threshold=cfg['fusion_thresh'], debug=cfg['debug']) if cfg['clean_intermediate']: for f in height_maps: @@ -729,7 +733,7 @@ def main(user_cfg, start_from=0, merge_matches=False): common.print_elapsed_time() # cleanup - common.garbage_cleanup() + common.garbage_cleanup(cfg['clean_tmp']) common.print_elapsed_time(since_first_call=True) diff --git a/s2p/block_matching.py b/s2p/block_matching.py index b61baa73..388db431 100644 --- a/s2p/block_matching.py +++ b/s2p/block_matching.py @@ -8,14 +8,13 @@ import rasterio from s2p import common -from s2p.config import cfg class MaxDisparityRangeError(Exception): pass -def create_rejection_mask(disp, im1, im2, mask): +def create_rejection_mask(disp, im1, im2, mask, temporary_dir): """ Create rejection mask (0 means rejected, 1 means accepted) Keep only the points that are matched and present in both input images @@ -25,14 +24,14 @@ def create_rejection_mask(disp, im1, im2, mask): im1, im2: rectified stereo pair mask: path to the output rejection mask """ - tmp1 = common.tmpfile('.tif') - tmp2 = common.tmpfile('.tif') + tmp1 = common.tmpfile(temporary_dir, '.tif') + tmp2 = common.tmpfile(temporary_dir, '.tif') common.run(["plambda", disp, "x 0 join", "-o", tmp1]) common.run(["backflow", tmp1, im2, tmp2]) common.run(["plambda", disp, im1, tmp2, "x isfinite y isfinite z isfinite and and vmul", "-o", mask]) -def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, +def compute_disparity_map(im1, im2, disp, mask, algo, cfg, disp_min=None, disp_max=None, timeout=600, max_disp_range=None, extra_params=''): """ @@ -124,14 +123,14 @@ def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, win = 3 # matched block size. It must be a positive odd number lr = 1 # maximum difference allowed in the left-right disparity check - cost = common.tmpfile('.tif') + cost = common.tmpfile(cfg['temporary_dir'], '.tif') common.run('sgbm {} {} {} {} {} {} {} {} {} {}'.format(im1, im2, disp, cost, disp_min, disp_max, win, p1, p2, lr)) - create_rejection_mask(disp, im1, im2, mask) + create_rejection_mask(disp, im1, im2, mask, cfg['temporary_dir']) if algo == 'tvl1': tvl1 = 'callTVL1.sh' @@ -185,14 +184,14 @@ def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, timeout=timeout, ) - create_rejection_mask(disp, im1, im2, mask) + create_rejection_mask(disp, im1, im2, mask, cfg['temporary_dir']) if algo == 'mgm_multi_lsd': ref = im1 sec = im2 - wref = common.tmpfile('.tif') - wsec = common.tmpfile('.tif') + wref = common.tmpfile(cfg['temporary_dir'], '.tif') + wsec = common.tmpfile(cfg['temporary_dir'], '.tif') # TODO TUNE LSD PARAMETERS TO HANDLE DIRECTLY 12 bits images? # image dependent weights based on lsd segments with rasterio.open(ref, "r") as f: @@ -263,7 +262,7 @@ def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, timeout=timeout, ) - create_rejection_mask(disp, im1, im2, mask) + create_rejection_mask(disp, im1, im2, mask, cfg['temporary_dir']) if algo == 'mgm_multi': @@ -307,7 +306,7 @@ def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, timeout=timeout, ) - create_rejection_mask(disp, im1, im2, mask) + create_rejection_mask(disp, im1, im2, mask, cfg['temporary_dir']) if (algo == 'micmac'): # add micmac binaries to the PATH environment variable diff --git a/s2p/common.py b/s2p/common.py index c1bf3a4b..c4628e8f 100644 --- a/s2p/common.py +++ b/s2p/common.py @@ -15,8 +15,6 @@ import rasterio -from s2p.config import cfg - # silent rasterio NotGeoreferencedWarning warnings.filterwarnings("ignore", @@ -38,18 +36,18 @@ def remove(target): except OSError: pass -def garbage_cleanup(): +def garbage_cleanup(clean_tmp=True): """ Removes all the files listed in the global variable 'garbage'. """ - if cfg['clean_tmp']: + if clean_tmp: while garbage: remove(garbage.pop()) -def tmpfile(ext=''): +def tmpfile(temporary_dir, ext=''): """ - Creates a temporary file in the cfg['temporary_dir'] directory. + Creates a temporary file in the temporary_dir directory. Args: ext: desired file extension. The dot has to be included. @@ -61,7 +59,7 @@ def tmpfile(ext=''): at the end of the pipeline. """ fd, out = tempfile.mkstemp(suffix=ext, prefix='s2p_', - dir=os.path.expandvars(cfg['temporary_dir'])) + dir=os.path.expandvars(temporary_dir)) os.close(fd) # http://www.logilab.org/blogentry/17873 garbage.append(out) return out diff --git a/s2p/fusion.py b/s2p/fusion.py index 17662baf..64c901f3 100644 --- a/s2p/fusion.py +++ b/s2p/fusion.py @@ -9,7 +9,6 @@ import numpy as np import rasterio -from s2p.config import cfg from s2p import common @@ -22,7 +21,7 @@ def average_if_close(x, threshold): return np.nanmedian(x) -def merge_n(output, inputs, offsets, averaging='average_if_close', threshold=1): +def merge_n(output, inputs, offsets, averaging='average_if_close', threshold=1, debug=False): """ Merge n images of equal sizes by taking the median/mean/min/max pixelwise. @@ -47,7 +46,7 @@ def merge_n(output, inputs, offsets, averaging='average_if_close', threshold=1): for i, img in enumerate(inputs): with rasterio.open(img, 'r') as f: x[:, :, i] = f.read(1) - offsets[i] - if cfg['debug']: + if debug: common.rasterio_write('{}_registered.tif'.format(os.path.splitext(img)[0]), x[:, :, i] + np.mean(offsets)) diff --git a/s2p/initialization.py b/s2p/initialization.py index a85072d0..b76e820f 100644 --- a/s2p/initialization.py +++ b/s2p/initialization.py @@ -100,12 +100,10 @@ def check_parameters(d): def build_cfg(cfg, user_cfg): """ - Populate a dictionary containing the s2p parameters from a user config file. - - This dictionary is contained in the global variable 'cfg' of the config - module. + Populate a default cfg dictionary with the s2p parameters from a user config file. Args: + cfg: default config dictionary user_cfg: user config dictionary """ # check that all the mandatory arguments are defined @@ -312,7 +310,7 @@ def is_tile_all_nodata(path:str, window:rasterio.windows.Window): return False -def is_this_tile_useful(x, y, w, h, images, images_sizes, border_margin): +def is_this_tile_useful(x, y, w, h, images, images_sizes, border_margin, cfg): """ Check if a tile contains valid pixels. @@ -336,7 +334,7 @@ def is_this_tile_useful(x, y, w, h, images, images_sizes, border_margin): rpc = images[0]['rpcm'] for img, size in zip(images[1:], images_sizes[1:]): - coords = rpc_utils.corresponding_roi(rpc, img['rpcm'], x, y, w, h) + coords = rpc_utils.corresponding_roi(rpc, img['rpcm'], x, y, w, h, cfg) if rectangles_intersect(coords, (0, 0, size[1], size[0])): break # the tile is partly contained else: # we've reached the end of the loop hence the tile is not contained @@ -346,7 +344,7 @@ def is_this_tile_useful(x, y, w, h, images, images_sizes, border_margin): cld_msk = images[0]['cld'] wat_msk = images[0]['wat'] mask = masking.image_tile_mask(x, y, w, h, roi_msk, cld_msk, wat_msk, - images_sizes[0], border_margin=border_margin) + images_sizes[0], border_margin=border_margin, temporary_dir=cfg['temporary_dir']) if not mask.any(): return False, None return True, mask @@ -388,6 +386,7 @@ def tiles_full_info(cfg, tw, th, tiles_txt, create_masks=False): cfg['images'], images_sizes, cfg['border_margin'], + cfg, tilewise=False, timeout=cfg['timeout'], ) diff --git a/s2p/masking.py b/s2p/masking.py index 67a42676..9e23743d 100644 --- a/s2p/masking.py +++ b/s2p/masking.py @@ -16,7 +16,7 @@ def image_tile_mask(x, y, w, h, roi_gml=None, cld_gml=None, raster_mask=None, - img_shape=None, border_margin=10): + img_shape=None, border_margin=10, temporary_dir=None): """ Compute a validity mask for an image tile from vector/raster image masks. @@ -44,7 +44,7 @@ def image_tile_mask(x, y, w, h, roi_gml=None, cld_gml=None, raster_mask=None, mask = np.ones((h, w), dtype=bool) if roi_gml is not None: # image domain mask (polygons) - tmp = common.tmpfile('.png') + tmp = common.tmpfile(temporary_dir, '.png') subprocess.check_call('cldmask %d %d -h "%s" %s %s' % (w, h, hij, roi_gml, tmp), shell=True) @@ -55,7 +55,7 @@ def image_tile_mask(x, y, w, h, roi_gml=None, cld_gml=None, raster_mask=None, return mask if cld_gml is not None: # cloud mask (polygons) - tmp = common.tmpfile('.png') + tmp = common.tmpfile(temporary_dir, '.png') subprocess.check_call('cldmask %d %d -h "%s" %s %s' % (w, h, hij, cld_gml, tmp), shell=True) diff --git a/s2p/pointing_accuracy.py b/s2p/pointing_accuracy.py index 1b1b8344..0c825df7 100644 --- a/s2p/pointing_accuracy.py +++ b/s2p/pointing_accuracy.py @@ -10,7 +10,6 @@ from s2p import sift from s2p import rpc_utils from s2p import estimation -from s2p.config import cfg def error_vectors(m, F, ind='ref'): @@ -60,7 +59,7 @@ def error_vectors(m, F, ind='ref'): return np.vstack((np.multiply(a, l[:, 0]), np.multiply(a, l[:, 1]))).T -def local_translation(r1, r2, x, y, w, h, m): +def local_translation(r1, r2, x, y, w, h, m, cfg): """ Estimates the optimal translation to minimise the relative pointing error on a given tile. @@ -80,7 +79,7 @@ def local_translation(r1, r2, x, y, w, h, m): """ # estimate the affine fundamental matrix between the two views n = cfg['n_gcp_per_axis'] - rpc_matches = rpc_utils.matches_from_rpc(r1, r2, x, y, w, h, n) + rpc_matches = rpc_utils.matches_from_rpc(r1, r2, x, y, w, h, n, cfg) F = estimation.affine_fundamental_matrix(rpc_matches) # compute the error vectors @@ -101,7 +100,7 @@ def local_translation(r1, r2, x, y, w, h, m): def compute_correction(img1, img2, rpc1, rpc2, x, y, w, h, method, sift_thresh, epipolar_threshold, - matching_method, min_value, max_value, confidence_threshold): + matching_method, min_value, max_value, confidence_threshold, cfg): """ Computes pointing correction matrix for specific ROI @@ -126,11 +125,11 @@ def compute_correction(img1, img2, rpc1, rpc2, x, y, w, h, raise ValueError(f"width or height <= 0 for:\n{img1}\n{img2}\nx={x}, y={y}, w={w}, h={h}. Try a different" f"tilesize or different ROI.") m = sift.matches_on_rpc_roi(img1, img2, rpc1, rpc2, x, y, w, h, - method, sift_thresh, epipolar_threshold, + method, sift_thresh, epipolar_threshold, cfg, matching_method, min_value, max_value, confidence_threshold) if m is not None: - A = local_translation(rpc1, rpc2, x, y, w, h, m) + A = local_translation(rpc1, rpc2, x, y, w, h, m, cfg) else: A = None diff --git a/s2p/rectification.py b/s2p/rectification.py index 49aa65ea..fad27e3e 100644 --- a/s2p/rectification.py +++ b/s2p/rectification.py @@ -13,7 +13,6 @@ from s2p import evaluation from s2p import common from s2p import visualisation -from s2p.config import cfg class NoRectificationMatchesError(Exception): @@ -50,7 +49,7 @@ def filter_matches_epipolar_constraint(F, matches, thresh): return np.array(out) -def register_horizontally_shear(matches, H1, H2): +def register_horizontally_shear(matches, H1, H2, debug=False): """ Adjust rectifying homographies with tilt, shear and translation to reduce the disparity range. @@ -73,7 +72,7 @@ def register_horizontally_shear(matches, H1, H2): x2 = p2[:, 0] y2 = p2[:, 1] - if cfg['debug']: + if debug: print("Residual vertical disparities: max, min, mean. Should be zero") print(np.max(y2 - y1), np.min(y2 - y1), np.mean(y2 - y1)) @@ -86,7 +85,7 @@ def register_horizontally_shear(matches, H1, H2): return np.dot(np.array([[a, b, c], [0, 1, 0], [0, 0, 1]]), H2) -def register_horizontally_translation(matches, H1, H2, flag='center'): +def register_horizontally_translation(matches, H1, H2, flag='center', debug=False): """ Adjust rectifying homographies with a translation to modify the disparity range. @@ -117,7 +116,7 @@ def register_horizontally_translation(matches, H1, H2, flag='center'): y2 = p2[:, 1] # for debug, print the vertical disparities. Should be zero. - if cfg['debug']: + if debug: print("Residual vertical disparities: max, min, mean. Should be zero") print(np.max(y2 - y1), np.min(y2 - y1), np.mean(y2 - y1)) @@ -134,7 +133,7 @@ def register_horizontally_translation(matches, H1, H2, flag='center'): return np.dot(common.matrix_translation(-t, 0), H2) -def disparity_range_from_matches(matches, H1, H2, w, h): +def disparity_range_from_matches(matches, H1, H2, w, h, disp_range_extra_margin=0.2): """ Compute the disparity range of a ROI from a list of point matches. @@ -158,12 +157,12 @@ def disparity_range_from_matches(matches, H1, H2, w, h): disp_max = np.ceil(np.max(x2 - x1)) # add a security margin to the disparity range - disp_min -= (disp_max - disp_min) * cfg['disp_range_extra_margin'] - disp_max += (disp_max - disp_min) * cfg['disp_range_extra_margin'] + disp_min -= (disp_max - disp_min) * disp_range_extra_margin + disp_max += (disp_max - disp_min) * disp_range_extra_margin return disp_min, disp_max -def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, A=None): +def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, cfg, A=None): """ Compute the disparity range of a ROI from a list of point matches. @@ -186,7 +185,7 @@ def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, A=None): if cfg['disp_range_method'] in ['exogenous', 'wider_sift_exogenous']: exogenous_disp = rpc_utils.exogenous_disp_range_estimation(rpc1, rpc2, x, y, w, h, - H1, H2, A, + H1, H2, cfg, A, cfg['disp_range_exogenous_high_margin'], cfg['disp_range_exogenous_low_margin']) @@ -195,7 +194,7 @@ def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, A=None): # compute SIFT disparity range if needed if cfg['disp_range_method'] in ['sift', 'wider_sift_exogenous']: if matches is not None and len(matches) >= 2: - sift_disp = disparity_range_from_matches(matches, H1, H2, w, h) + sift_disp = disparity_range_from_matches(matches, H1, H2, w, h, cfg['disp_range_extra_margin']) else: sift_disp = None print("SIFT disparity range:", sift_disp) @@ -239,7 +238,7 @@ def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, A=None): return disp -def rectification_homographies(matches, x, y, w, h): +def rectification_homographies(matches, x, y, w, h, debug=False): """ Computes rectifying homographies from point matches for a given ROI. @@ -262,9 +261,9 @@ def rectification_homographies(matches, x, y, w, h): F = estimation.affine_fundamental_matrix(matches) # compute rectifying similarities - S1, S2 = estimation.rectifying_similarities_from_affine_fundamental_matrix(F, cfg['debug']) + S1, S2 = estimation.rectifying_similarities_from_affine_fundamental_matrix(F, debug) - if cfg['debug']: + if debug: y1 = common.points_apply_homography(S1, matches[:, :2])[:, 1] y2 = common.points_apply_homography(S2, matches[:, 2:])[:, 1] err = np.abs(y1 - y2) @@ -278,8 +277,7 @@ def rectification_homographies(matches, x, y, w, h): return np.dot(T, S1), np.dot(T, S2), F -def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matches=None, - method='rpc', hmargin=0, vmargin=0): +def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, cfg, A=None, sift_matches=None, method='rpc', hmargin=0, vmargin=0): """ Rectify a ROI in a pair of images. @@ -308,8 +306,7 @@ def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matc # compute real or virtual matches if method == 'rpc': # find virtual matches from RPC camera models - matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, - cfg['n_gcp_per_axis']) + matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, cfg['n_gcp_per_axis'], cfg) # correct second image coordinates with the pointing correction matrix if A is not None: @@ -326,8 +323,7 @@ def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matc # If there are still no matches, raise the no matches error. if matches is None or len(matches) < 4 or len(matches.shape) != 2: # find virtual matches from RPC camera models - matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, - cfg['n_gcp_per_axis']) + matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, cfg['n_gcp_per_axis'], cfg) # correct second image coordinates with the pointing correction matrix if A is not None: @@ -335,29 +331,28 @@ def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matc matches[:, 2:]) # compute rectifying homographies - H1, H2, F = rectification_homographies(matches, x, y, w, h) + H1, H2, F = rectification_homographies(matches, x, y, w, h, debug=cfg['debug']) if cfg['register_with_shear']: # compose H2 with a horizontal shear to reduce the disparity range - a = np.mean(rpc_utils.altitude_range(rpc1, x, y, w, h)) + a = np.mean(rpc_utils.altitude_range(rpc1, x, y, w, h, cfg)) lon, lat, alt = rpc_utils.ground_control_points(rpc1, x, y, w, h, a, a, 4) x1, y1 = rpc1.projection(lon, lat, alt) x2, y2 = rpc2.projection(lon, lat, alt) m = np.vstack([x1, y1, x2, y2]).T m = np.vstack(list({tuple(row) for row in m})) # remove duplicates due to no alt range - H2 = register_horizontally_shear(m, H1, H2) + H2 = register_horizontally_shear(m, H1, H2, debug=cfg['debug']) # compose H2 with a horizontal translation to center disp range around 0 if sift_matches is not None: - sift_matches = filter_matches_epipolar_constraint(F, sift_matches, - cfg['epipolar_thresh']) + sift_matches = filter_matches_epipolar_constraint(F, sift_matches, cfg['epipolar_thresh']) if len(sift_matches) < 1: warnings.warn( "Need at least one sift match for the horizontal registration", category=NoHorizontalRegistrationWarning, ) else: - H2 = register_horizontally_translation(sift_matches, H1, H2) + H2 = register_horizontally_translation(sift_matches, H1, H2, debug=cfg['debug']) # compute disparity range if cfg['debug']: @@ -367,9 +362,10 @@ def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matc visualisation.plot_matches(im1, im2, rpc1, rpc2, sift_matches, os.path.join(out_dir, 'sift_matches_disp.png'), - x, y, w, h) - disp_m, disp_M = disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, - sift_matches, A) + x, y, w, h, cfg) + disp_m, disp_M = disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, + sift_matches, cfg, A=A + ) # recompute hmargin and homographies hmargin = int(np.ceil(max([hmargin, np.fabs(disp_m), np.fabs(disp_M)]))) diff --git a/s2p/rpc_utils.py b/s2p/rpc_utils.py index f93b14a4..0cc524dd 100644 --- a/s2p/rpc_utils.py +++ b/s2p/rpc_utils.py @@ -12,7 +12,6 @@ from s2p import geographiclib from s2p import common -from s2p.config import cfg warnings.filterwarnings("ignore", category=rasterio.errors.NotGeoreferencedWarning) @@ -89,7 +88,7 @@ def altitude_range_coarse(rpc, scale_factor=1): return m, M -def min_max_heights_from_bbx(im, lon_m, lon_M, lat_m, lat_M, rpc): +def min_max_heights_from_bbx(im, lon_m, lon_M, lat_m, lat_M, rpc, exogenous_dem_geoid_mode=True, rpc_alt_range_scale_factor=1): """ Compute min, max heights from bounding box @@ -139,7 +138,7 @@ def min_max_heights_from_bbx(im, lon_m, lon_M, lat_m, lat_M, rpc): hmin = np.nanmin(array) hmax = np.nanmax(array) - if cfg['exogenous_dem_geoid_mode'] is True: + if exogenous_dem_geoid_mode is True: offset = geographiclib.geoid_to_ellipsoid((lat_m + lat_M)/2, (lon_m + lon_M)/2, 0) hmin += offset hmax += offset @@ -147,10 +146,10 @@ def min_max_heights_from_bbx(im, lon_m, lon_M, lat_m, lat_M, rpc): else: print("WARNING: rpc_utils.min_max_heights_from_bbx: access window out of range") print("returning coarse range from rpc") - return altitude_range_coarse(rpc, cfg['rpc_alt_range_scale_factor']) + return altitude_range_coarse(rpc, rpc_alt_range_scale_factor) -def altitude_range(rpc, x, y, w, h, margin_top=0, margin_bottom=0): +def altitude_range(rpc, x, y, w, h, cfg, margin_top=0, margin_bottom=0): """ Computes an altitude range using the exogenous dem. @@ -180,7 +179,8 @@ def altitude_range(rpc, x, y, w, h, margin_top=0, margin_bottom=0): # compute heights on this bounding box if cfg['exogenous_dem'] is not None: h_m, h_M = min_max_heights_from_bbx(cfg['exogenous_dem'], - lon_m, lon_M, lat_m, lat_M, rpc) + lon_m, lon_M, lat_m, lat_M, rpc, + cfg['exogenous_dem_geoid_mode'], cfg['rpc_alt_range_scale_factor']) h_m += margin_bottom h_M += margin_top elif cfg['use_srtm']: @@ -192,6 +192,8 @@ def altitude_range(rpc, x, y, w, h, margin_top=0, margin_bottom=0): h_m = min(alts) + margin_bottom h_M = max(alts) + margin_top else: + print(f"WARNING: Using coarse altitude range because use_srtm=False and no exogenous DEM is provided") + print(f"This may lead to poor results") h_m, h_M = altitude_range_coarse(rpc, cfg['rpc_alt_range_scale_factor']) return h_m, h_M @@ -316,7 +318,7 @@ def ground_control_points(rpc, x, y, w, h, m, M, n): return lon, lat, alt -def corresponding_roi(rpc1, rpc2, x, y, w, h): +def corresponding_roi(rpc1, rpc2, x, y, w, h, cfg): """ Uses RPC functions to determine the region of im2 associated to the specified ROI of im1. @@ -338,7 +340,7 @@ def corresponding_roi(rpc1, rpc2, x, y, w, h): rpc1 = rpcm.RPCModel(rpc1) if not isinstance(rpc2, rpcm.RPCModel): rpc2 = rpcm.RPCModel(rpc2) - m, M = altitude_range(rpc1, x, y, w, h, 0, 0) + m, M = altitude_range(rpc1, x, y, w, h, cfg, 0, 0) # build an array with vertices of the 3D ROI, obtained as {2D ROI} x [m, M] a = np.array([x, x, x, x, x+w, x+w, x+w, x+w]) @@ -353,7 +355,7 @@ def corresponding_roi(rpc1, rpc2, x, y, w, h): return np.round(out) -def matches_from_rpc(rpc1, rpc2, x, y, w, h, n): +def matches_from_rpc(rpc1, rpc2, x, y, w, h, n, cfg): """ Uses RPC functions to generate matches between two Pleiades images. @@ -368,7 +370,7 @@ def matches_from_rpc(rpc1, rpc2, x, y, w, h, n): Returns: an array of matches, one per line, expressed as x1, y1, x2, y2. """ - m, M = altitude_range(rpc1, x, y, w, h, 100, -100) + m, M = altitude_range(rpc1, x, y, w, h, cfg, 100, -100) lon, lat, alt = ground_control_points(rpc1, x, y, w, h, m, M, n) x1, y1 = rpc1.projection(lon, lat, alt) x2, y2 = rpc2.projection(lon, lat, alt) @@ -411,7 +413,7 @@ def alt_to_disp(rpc1, rpc2, x, y, alt, H1, H2, A=None): return disp -def exogenous_disp_range_estimation(rpc1, rpc2, x, y, w, h, H1, H2, A=None, +def exogenous_disp_range_estimation(rpc1, rpc2, x, y, w, h, H1, H2, cfg, A=None, margin_top=0, margin_bottom=0): """ Args: @@ -436,7 +438,7 @@ def exogenous_disp_range_estimation(rpc1, rpc2, x, y, w, h, H1, H2, A=None, if cfg['exogenous_dem'] is None: return - m, M = altitude_range(rpc1, x, y, w, h, margin_top, margin_bottom) + m, M = altitude_range(rpc1, x, y, w, h, cfg, margin_top, margin_bottom) return altitude_range_to_disp_range(m, M, rpc1, rpc2, x, y, w, h, H1, H2, A, margin_top, margin_bottom) diff --git a/s2p/sift.py b/s2p/sift.py index 12aa06e2..c7973cfe 100644 --- a/s2p/sift.py +++ b/s2p/sift.py @@ -254,7 +254,7 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold, def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h, - method, sift_thresh, epipolar_threshold, matching_method="sift", + method, sift_thresh, epipolar_threshold, cfg, matching_method="sift", min_value=200, max_value=3000, confidence_threshold=0.5): """ Compute a list of SIFT matches between two images on a given roi. @@ -273,10 +273,10 @@ def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h, contains one pair of points, ordered as x1 y1 x2 y2. The coordinate system is that of the full images. """ - x2, y2, w2, h2 = rpc_utils.corresponding_roi(rpc1, rpc2, x, y, w, h) + x2, y2, w2, h2 = rpc_utils.corresponding_roi(rpc1, rpc2, x, y, w, h, cfg) # estimate an approximate affine fundamental matrix from the rpcs - rpc_matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, 5) + rpc_matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, 5, cfg) F = estimation.affine_fundamental_matrix(rpc_matches) if matching_method == "loftr": diff --git a/s2p/triangulation.py b/s2p/triangulation.py index a04d8b2b..cd10f439 100644 --- a/s2p/triangulation.py +++ b/s2p/triangulation.py @@ -11,7 +11,6 @@ import rasterio from s2p import common -from s2p.config import cfg from s2p import ply from s2p import geographiclib diff --git a/s2p/visualisation.py b/s2p/visualisation.py index 7c6b1c98..ef36269a 100644 --- a/s2p/visualisation.py +++ b/s2p/visualisation.py @@ -107,7 +107,7 @@ def plot_matches_low_level(crop1, crop2, matches, outfile, max_matches=100): common.rasterio_write(outfile, out) -def plot_matches(img1, img2, rpc1, rpc2, matches, outfile, x, y, w, h): +def plot_matches(img1, img2, rpc1, rpc2, matches, outfile, x, y, w, h, cfg): """ Plot keypoint matches on images corresponding ROIs. @@ -129,7 +129,7 @@ def plot_matches(img1, img2, rpc1, rpc2, matches, outfile, x, y, w, h): return x1, y1, w1, h1 = x, y, w, h - x2, y2, w2, h2 = map(int, rpc_utils.corresponding_roi(rpc1, rpc2, x1, y1, w1, h1)) + x2, y2, w2, h2 = map(int, rpc_utils.corresponding_roi(rpc1, rpc2, x1, y1, w1, h1, cfg)) # read the crops with rasterio.open(img1, "r") as f: diff --git a/setup.py b/setup.py index 0c500cb6..dd4d3056 100644 --- a/setup.py +++ b/setup.py @@ -63,7 +63,7 @@ def finalize_options(self): } setup(name="s2p", - version="1.6.1", + version="1.6.2", description="Satellite Stereo Pipeline.", long_description=readme(), long_description_content_type='text/markdown', diff --git a/tests/block_matching_test.py b/tests/block_matching_test.py index f6814efc..bac3e23b 100644 --- a/tests/block_matching_test.py +++ b/tests/block_matching_test.py @@ -4,6 +4,7 @@ import pytest import s2p +from s2p.config import cfg from tests_utils import data_path @@ -17,7 +18,7 @@ def test_compute_disparity_map_timeout(timeout=1): with pytest.raises(subprocess.TimeoutExpired): s2p.block_matching.compute_disparity_map(img, img, disp, mask, - "mgm_multi", -100, 100, + "mgm_multi", cfg, -100, 100, timeout) @@ -32,5 +33,5 @@ def test_compute_disparity_map_max_disp_range(max_disp_range=10): with pytest.raises(s2p.block_matching.MaxDisparityRangeError): s2p.block_matching.compute_disparity_map(img, img, disp, mask, - "mgm_multi", -100, 100, + "mgm_multi", cfg, -100, 100, max_disp_range=max_disp_range) diff --git a/tests/rectification_test.py b/tests/rectification_test.py index f54675c4..d7e00340 100644 --- a/tests/rectification_test.py +++ b/tests/rectification_test.py @@ -5,6 +5,7 @@ from rpcm import rpc_from_geotiff import s2p +from s2p.config import cfg from tests_utils import data_path @@ -55,6 +56,7 @@ def test_rectify_pair_with_matches(tmp_path, matches, images): h=200, out1=str(tmp_path / 'out1.tiff'), out2=str(tmp_path / 'out2.tiff'), + cfg=cfg, sift_matches=matches, method='sift', ) diff --git a/tests/rpc_utils_test.py b/tests/rpc_utils_test.py index 42996efa..5742f1c5 100644 --- a/tests/rpc_utils_test.py +++ b/tests/rpc_utils_test.py @@ -7,6 +7,7 @@ from s2p import rpc_utils from tests_utils import data_path +from s2p.config import cfg def test_matches_from_rpc(): @@ -16,7 +17,7 @@ def test_matches_from_rpc(): r1 = rpcm.rpc_from_geotiff(data_path(os.path.join('input_pair', 'img_01.tif'))) r2 = rpcm.rpc_from_geotiff(data_path(os.path.join('input_pair', 'img_02.tif'))) - test_matches = rpc_utils.matches_from_rpc(r1, r2, 100, 100, 200, 200, 5) + test_matches = rpc_utils.matches_from_rpc(r1, r2, 100, 100, 200, 200, 5, cfg) expected_matches = np.loadtxt(data_path(os.path.join('expected_output', 'units', 'unit_matches_from_rpc.txt'))) diff --git a/tests/sift_test.py b/tests/sift_test.py index f9605b56..adb1b297 100644 --- a/tests/sift_test.py +++ b/tests/sift_test.py @@ -7,6 +7,7 @@ from s2p import sift from tests_utils import data_path +from s2p.config import cfg def test_image_keypoints(): @@ -43,7 +44,7 @@ def test_matches_on_rpc_roi(): rpc2 = rpcm.rpc_from_geotiff(img2) computed = sift.matches_on_rpc_roi( img1, img2, rpc1, rpc2, 100, 100, 200, 200, - method='relative', sift_thresh=0.6, epipolar_threshold=10 + method='relative', sift_thresh=0.6, epipolar_threshold=10, cfg=cfg ) expected = np.loadtxt(data_path('expected_output/units/matches_on_rpc_roi.txt')) np.testing.assert_allclose(computed, expected, rtol=0.01, atol=0.1, diff --git a/utils/svg_tilemap.py b/utils/svg_tilemap.py index c9659825..50a28599 100755 --- a/utils/svg_tilemap.py +++ b/utils/svg_tilemap.py @@ -90,7 +90,7 @@ def main(user_cfg): # cleanup - common.garbage_cleanup() + common.garbage_cleanup(cfg['clean_temp']) if __name__ == '__main__':