diff --git a/s2p/.ipynb_checkpoints/__init__-checkpoint.py b/s2p/.ipynb_checkpoints/__init__-checkpoint.py new file mode 100644 index 00000000..c201002c --- /dev/null +++ b/s2p/.ipynb_checkpoints/__init__-checkpoint.py @@ -0,0 +1,778 @@ +#!/usr/bin/env python + +# s2p - Satellite Stereo Pipeline +# Copyright (C) 2015, Carlo de Franchis +# Copyright (C) 2015, Gabriele Facciolo +# Copyright (C) 2015, Enric Meinhardt +# Copyright (C) 2015, Julien Michel + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero General Public License as published +# by the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +import subprocess +import sys +import os.path +import json +import datetime +import multiprocessing +import glob +import os + +import numpy as np +import rasterio +import rasterio.merge +from plyflatten import plyflatten_from_plyfiles_list + + +from s2p.config import cfg +from s2p import common +from s2p import parallel +from s2p import geographiclib +from s2p import initialization +from s2p import pointing_accuracy +from s2p import rectification +from s2p import block_matching +from s2p import masking +from s2p import ply +from s2p import triangulation +from s2p import fusion +from s2p import visualisation +from s2p.image_coordinates_to_coordinates import matches_to_geojson + + +def remove_missing_tiles(tiles): + """Remove tiles where any of the rectified images is missing (also remove from neighborhood_dirs).""" + n = len(cfg['images']) + tiles_new = [] + deleted_tiles = [] + for tile in tiles: + paths = [] + for i in range(1, n): + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + paths += [ + os.path.join(out_dir, 'rectified_ref.tif'), + os.path.join(out_dir, 'rectified_sec.tif') + ] + missing = sum(not os.path.exists(p) for p in paths) + if missing > 0: + relative_path = os.path.join('../../..', initialization.get_tile_dir(*tile['coordinates'])) + print(f"WARNING: tile {relative_path} is missing {missing}/{len(paths)} " + "input files, removing tile...") + deleted_tiles.append(relative_path) + continue + tiles_new.append(tile) + + # Remove deleted tiles from neighborhood_dirs + for i in range(len(tiles_new)): + tiles_new[i]['neighborhood_dirs'] = [d for d in tiles_new[i]['neighborhood_dirs'] if d not in deleted_tiles] + + tiles_pairs = [(t, i) for i in range(1, n) for t in tiles_new] + return tiles_new, tiles_pairs + + +def check_missing_sift(tiles_pairs): + missing_sift = [] + with open(os.path.join(cfg["out_dir"], "missing_sift.txt"), "w") as f: + for tile, i in tiles_pairs: + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + path = os.path.join(out_dir, 'sift_matches.txt') + if not os.path.exists(path): + missing_sift.append(path) + f.write(path + "\n") + if len(missing_sift) > 0: + print(" --- ") + print(f"WARNING: missing {len(missing_sift)}/{len(tiles_pairs)} " + "SIFT matches, this may deteriorate output quality") + + +def merge_all_match_files(): + """ + Merges together all non-empty match files into one .txt file and saves them in the output directory. + """ + read_files = glob.glob("**/*sift_matches.txt", recursive=True) + with open(os.path.join(cfg["out_dir"], "merged_sift_matches.txt"), "w") as outfile: + for f in read_files: + if os.stat(f).st_size == 0: + continue + with open(f, "r") as infile: + outfile.write(infile.read()) + + +def pointing_correction(tile, i): + """ + Compute the translation that corrects the pointing error on a pair of tiles. + + Args: + tile: dictionary containing the information needed to process the tile + i: index of the processed pair + """ + x, y, w, h = tile['coordinates'] + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + img1 = cfg['images'][0]['img'] + rpc1 = cfg['images'][0]['rpcm'] + img2 = cfg['images'][i]['img'] + rpc2 = cfg['images'][i]['rpcm'] + + # correct pointing error + print('correcting pointing on tile {} {} pair {}...'.format(x, y, i)) + method = 'relative' if cfg['relative_sift_match_thresh'] is True else 'absolute' + 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'] + ) + + if A is not None: # A is the correction matrix + np.savetxt(os.path.join(out_dir, 'pointing.txt'), A, fmt='%6.3f') + if m is not None: # m is the list of sift matches + np.savetxt(os.path.join(out_dir, 'sift_matches.txt'), m, fmt='%9.3f') + np.savetxt(os.path.join(out_dir, 'center_keypts_sec.txt'), + np.mean(m[:, 2:], 0), fmt='%9.3f') + if cfg['debug']: + visualisation.plot_matches(img1, img2, rpc1, rpc2, m, + os.path.join(out_dir, + 'sift_matches_pointing.png'), + x, y, w, h) + + +def global_pointing_correction(tiles): + """ + Compute the global pointing corrections for each pair of images. + + Args: + tiles: list of tile dictionaries + """ + for i in range(1, len(cfg['images'])): + out = os.path.join(cfg['out_dir'], 'global_pointing_pair_%d.txt' % i) + l = [os.path.join(t['dir'], 'pair_%d' % i) for t in tiles] + np.savetxt(out, pointing_accuracy.global_from_local(l), + fmt='%12.6f') + if cfg['clean_intermediate']: + for d in l: + common.remove(os.path.join(d, 'center_keypts_sec.txt')) + + +def rectification_pair(tile, i): + """ + Rectify a pair of images on a given tile. + + Args: + tile: dictionary containing the information needed to process a tile. + i: index of the processed pair + """ + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + x, y, w, h = tile['coordinates'] + img1 = cfg['images'][0]['img'] + rpc1 = cfg['images'][0]['rpcm'] + img2 = cfg['images'][i]['img'] + rpc2 = cfg['images'][i]['rpcm'] + pointing = os.path.join(cfg['out_dir'], + 'global_pointing_pair_{}.txt'.format(i)) + + print('rectifying tile {} {} pair {}...'.format(x, y, i)) + try: + A = np.loadtxt(os.path.join(out_dir, 'pointing.txt')) + except IOError: + A = np.loadtxt(pointing) + try: + m = np.loadtxt(os.path.join(out_dir, 'sift_matches.txt')) + except IOError: + m = None + + cur_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + for n in tile['neighborhood_dirs']: + nei_dir = os.path.join(tile['dir'], n, 'pair_{}'.format(i)) + if os.path.exists(nei_dir) and not os.path.samefile(cur_dir, nei_dir): + sift_from_neighborhood = os.path.join(nei_dir, 'sift_matches.txt') + try: + m_n = np.loadtxt(sift_from_neighborhood) + # added sifts in the ellipse of semi axes : (3*w/4, 3*h/4) + m_n = m_n[np.where(np.linalg.norm([(m_n[:, 0] - (x + w/2)) / w, + (m_n[:, 1] - (y + h/2)) / h], + axis=0) < 3/4)] + if m is None: + m = m_n + else: + m = np.concatenate((m, m_n)) + except IOError: + print('%s does not exist' % sift_from_neighborhood) + + rect1 = os.path.join(out_dir, 'rectified_ref.tif') + rect2 = os.path.join(out_dir, 'rectified_sec.tif') + H1, H2, disp_min, disp_max = rectification.rectify_pair(img1, img2, + rpc1, rpc2, + x, y, w, h, + rect1, rect2, A, m, + method=cfg['rectification_method'], + hmargin=cfg['horizontal_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') + np.savetxt(os.path.join(out_dir, 'disp_min_max.txt'), [disp_min, disp_max], + fmt='%3.1f') + + if cfg['clean_intermediate']: + common.remove(os.path.join(out_dir, 'pointing.txt')) + common.remove(os.path.join(out_dir, 'sift_matches.txt')) + +def stereo_matching(tile, i): + """ + Compute the disparity of a pair of images on a given tile. + + Args: + tile: dictionary containing the information needed to process a tile. + i: index of the processed pair + """ + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + x, y = tile['coordinates'][:2] + + print('estimating disparity on tile {} {} pair {}...'.format(x, y, i)) + rect1 = os.path.join(out_dir, 'rectified_ref.tif') + rect2 = os.path.join(out_dir, 'rectified_sec.tif') + disp = os.path.join(out_dir, 'rectified_disp.tif') + mask = os.path.join(out_dir, 'rectified_mask.png') + + 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, + disp_max, timeout=cfg['mgm_timeout'], + max_disp_range=cfg['max_disp_range']) + + # add margin around masked pixels + masking.erosion(mask, mask, cfg['msk_erosion']) + + if cfg['clean_intermediate']: + if len(cfg['images']) > 2: + common.remove(rect1) + common.remove(rect2) + common.remove(os.path.join(out_dir, 'disp_min_max.txt')) + +def disparity_to_height(tile, i): + """ + Compute a height map from the disparity map of a pair of image tiles. + + Args: + tile: dictionary containing the information needed to process a tile. + i: index of the processed pair. + """ + out_dir = os.path.join(tile['dir'], 'pair_{}'.format(i)) + x, y, w, h = tile['coordinates'] + + print('triangulating tile {} {} pair {}...'.format(x, y, i)) + rpc1 = cfg['images'][0]['rpcm'] + rpc2 = cfg['images'][i]['rpcm'] + H_ref = np.loadtxt(os.path.join(out_dir, 'H_ref.txt')) + H_sec = np.loadtxt(os.path.join(out_dir, 'H_sec.txt')) + disp = os.path.join(out_dir, 'rectified_disp.tif') + mask = os.path.join(out_dir, 'rectified_mask.png') + mask_orig = os.path.join(tile['dir'], 'mask.png') + pointing = os.path.join(cfg['out_dir'], + 'global_pointing_pair_{}.txt'.format(i)) + + with rasterio.open(disp, 'r') as f: + disp_img = f.read().squeeze() + with rasterio.open(mask, 'r') as f: + mask_rect_img = f.read().squeeze() + with rasterio.open(mask_orig, 'r') as f: + mask_orig_img = f.read().squeeze() + height_map = triangulation.height_map(x, y, w, h, rpc1, rpc2, H_ref, H_sec, + disp_img, mask_rect_img, + mask_orig_img, + A=np.loadtxt(pointing)) + + # write height map to a file + common.rasterio_write(os.path.join(out_dir, 'height_map.tif'), height_map) + + if cfg['clean_intermediate']: + common.remove(H_ref) + common.remove(H_sec) + common.remove(disp) + common.remove(mask) + + +def disparity_to_ply(tile): + """ + Compute a point cloud from the disparity map of a pair of image tiles. + + This function is called by s2p.main only if there are two input images (not + three). + + Args: + tile: dictionary containing the information needed to process a tile. + """ + out_dir = tile['dir'] + ply_file = os.path.join(out_dir, 'cloud.ply') + x, y, w, h = tile['coordinates'] + rpc1 = cfg['images'][0]['rpcm'] + rpc2 = cfg['images'][1]['rpcm'] + + print('triangulating tile {} {}...'.format(x, y)) + H_ref = os.path.join(out_dir, 'pair_1', 'H_ref.txt') + H_sec = os.path.join(out_dir, 'pair_1', 'H_sec.txt') + pointing = os.path.join(cfg['out_dir'], 'global_pointing_pair_1.txt') + disp = os.path.join(out_dir, 'pair_1', 'rectified_disp.tif') + extra = os.path.join(out_dir, 'pair_1', 'rectified_disp_confidence.tif') + if not os.path.exists(extra): # confidence file not always generated + extra = '' + mask_rect = os.path.join(out_dir, 'pair_1', 'rectified_mask.png') + mask_orig = os.path.join(out_dir, 'mask.png') + + # prepare the image needed to colorize point cloud + if cfg['images'][0]['clr']: + # we want colors image and rectified_ref.tif to have the same size + 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") + success_state = common.image_apply_homography(colors, cfg['images'][0]['clr'], + np.loadtxt(H_ref), ww, hh) + if success_state is False: + # If homography fails the tile will be skipped (instead of crashing the entire process) + if cfg['clean_intermediate']: + common.remove(H_ref) + common.remove(H_sec) + common.remove(disp) + common.remove(mask_rect) + common.remove(mask_orig) + return None + with rasterio.open(colors, "r") as f: + colors = f.read() + + else: + with rasterio.open(os.path.join(out_dir, 'pair_1', 'rectified_ref.tif')) as f: + img = f.read() + colors = common.linear_stretching_and_quantization_8bit(img) + + # compute the point cloud + with rasterio.open(disp, 'r') as f: + disp_img = f.read().squeeze() + with rasterio.open(mask_rect, 'r') as f: + mask_rect_img = f.read().squeeze() + with rasterio.open(mask_orig, 'r') as f: + mask_orig_img = f.read().squeeze() + + out_crs = geographiclib.pyproj_crs(cfg['out_crs']) + xyz_array, err = triangulation.disp_to_xyz(rpc1, rpc2, + np.loadtxt(H_ref), np.loadtxt(H_sec), + disp_img, mask_rect_img, + img_bbx=(x, x+w, y, y+h), + mask_orig=mask_orig_img, + A=np.loadtxt(pointing), + out_crs=out_crs) + + # 3D filtering + r = cfg['3d_filtering_r'] + n = cfg['3d_filtering_n'] + if r and n: + triangulation.filter_xyz(xyz_array, r, n, cfg['gsd']) + + proj_com = "CRS {}".format(cfg['out_crs']) + triangulation.write_to_ply(ply_file, xyz_array, colors, proj_com, confidence=extra) + + if cfg['clean_intermediate']: + common.remove(H_ref) + common.remove(H_sec) + common.remove(disp) + common.remove(mask_rect) + common.remove(mask_orig) + common.remove(os.path.join(out_dir, 'pair_1', 'rectified_ref.tif')) + return tile + +def mean_heights(tile): + """ + """ + w, h = tile['coordinates'][2:] + n = len(cfg['images']) - 1 + maps = np.empty((h, w, n)) + for i in range(n): + try: + with rasterio.open(os.path.join(tile['dir'], 'pair_{}'.format(i + 1), + 'height_map.tif'), 'r') as f: + maps[:, :, i] = f.read(1) + except RuntimeError: # the file is not there + maps[:, :, i] *= np.nan + + validity_mask = maps.sum(axis=2) # sum to propagate nan values + validity_mask += 1 - validity_mask # 1 on valid pixels, and nan on invalid + + # save the n mean height values to a txt file in the tile directory + np.savetxt(os.path.join(tile['dir'], 'local_mean_heights.txt'), + [np.nanmean(validity_mask * maps[:, :, i]) for i in range(n)]) + + +def global_mean_heights(tiles): + """ + """ + local_mean_heights = [np.loadtxt(os.path.join(t['dir'], 'local_mean_heights.txt')) + for t in tiles] + global_mean_heights = np.nanmean(local_mean_heights, axis=0) + for i in range(len(cfg['images']) - 1): + np.savetxt(os.path.join(cfg['out_dir'], + 'global_mean_height_pair_{}.txt'.format(i+1)), + [global_mean_heights[i]]) + + +def heights_fusion(tile): + """ + Merge the height maps computed for each image pair and generate a ply cloud. + + Args: + tile: a dictionary that provides all you need to process a tile + """ + tile_dir = tile['dir'] + height_maps = [os.path.join(tile_dir, 'pair_%d' % (i + 1), 'height_map.tif') + for i in range(len(cfg['images']) - 1)] + + # remove spurious matches + if cfg['cargarse_basura']: + for img in height_maps: + common.cargarse_basura(img, img) + + # load global mean heights + global_mean_heights = [] + for i in range(len(cfg['images']) - 1): + x = np.loadtxt(os.path.join(cfg['out_dir'], + 'global_mean_height_pair_{}.txt'.format(i+1))) + global_mean_heights.append(x) + + # 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']) + + if cfg['clean_intermediate']: + for f in height_maps: + common.remove(f) + + +def heights_to_ply(tile): + """ + Generate a ply cloud. + + Args: + tile: a dictionary that provides all you need to process a tile + """ + # merge the n-1 height maps of the tile (n = nb of images) + heights_fusion(tile) + + # compute a ply from the merged height map + out_dir = tile['dir'] + x, y, w, h = tile['coordinates'] + plyfile = os.path.join(out_dir, 'cloud.ply') + height_map = os.path.join(out_dir, 'height_map.tif') + + if cfg['images'][0]['clr']: + with rasterio.open(cfg['images'][0]['clr'], "r") as f: + colors = f.read(window=((y, y + h), (x, x + w))) + else: + with rasterio.open(cfg['images'][0]['img'], "r") as f: + colors = f.read(window=((y, y + h), (x, x + w))) + + colors = common.linear_stretching_and_quantization_8bit(colors) + + out_crs = geographiclib.pyproj_crs(cfg['out_crs']) + xyz_array = triangulation.height_map_to_xyz(height_map, + cfg['images'][0]['rpcm'], x, y, + out_crs) + + # 3D filtering + r = cfg['3d_filtering_r'] + n = cfg['3d_filtering_n'] + if r and n: + triangulation.filter_xyz(xyz_array, r, n, cfg['gsd']) + + proj_com = "CRS {}".format(cfg['out_crs']) + triangulation.write_to_ply(plyfile, xyz_array, colors, proj_com) + + if cfg['clean_intermediate']: + common.remove(height_map) + common.remove(os.path.join(out_dir, 'mask.png')) + + +def plys_to_dsm(tile): + """ + Generates DSM from plyfiles (cloud.ply) + + Args: + tile: a dictionary that provides all you need to process a tile + """ + out_dsm = os.path.join(tile['dir'], 'dsm.tif') + out_conf = os.path.join(tile['dir'], 'confidence.tif') + r = cfg['dsm_resolution'] + + # compute the point cloud x, y bounds + points, _ = ply.read_3d_point_cloud_from_ply(os.path.join(tile['dir'], + 'cloud.ply')) + if len(points) == 0: + return + + xmin, ymin, *_ = np.min(points, axis=0) + xmax, ymax, *_ = np.max(points, axis=0) + + # compute xoff, yoff, xsize, ysize on a grid of unit r + xoff = np.floor(xmin / r) * r + xsize = int(1 + np.floor((xmax - xoff) / r)) + + yoff = np.ceil(ymax / r) * r + ysize = int(1 - np.floor((ymin - yoff) / r)) + + roi = xoff, yoff, xsize, ysize + + clouds = [os.path.join(tile['dir'], n_dir, 'cloud.ply') for n_dir in tile['neighborhood_dirs']] + raster, profile = plyflatten_from_plyfiles_list(clouds, + resolution=r, + roi=roi, + radius=cfg['dsm_radius'], + sigma=cfg['dsm_sigma']) + + # save output image with utm georeferencing + common.rasterio_write(out_dsm, raster[:, :, 0], profile=profile) + + # export confidence (optional) + # note that the plys are assumed to contain the fields: + # [x(float32), y(float32), z(float32), r(uint8), g(uint8), b(uint8), confidence(optional, float32)] + # so the raster has 4 or 5 columns: [z, r, g, b, confidence (optional)] + if raster.shape[-1] == 5: + common.rasterio_write(out_conf, raster[:, :, 4], profile=profile) + + +def global_dsm(tiles): + """ + Merge tilewise DSMs and confidence maps in a global DSM and confidence map. + """ + bounds = None + if "roi_geojson" in cfg: + ll_poly = geographiclib.read_lon_lat_poly_from_geojson(cfg["roi_geojson"]) + pyproj_crs = geographiclib.pyproj_crs(cfg["out_crs"]) + bounds = geographiclib.crs_bbx(ll_poly, pyproj_crs, + align=cfg["dsm_resolution"]) + + creation_options = {"tiled": True, + "blockxsize": 256, + "blockysize": 256, + "compress": "deflate", + "predictor": 2} + + dsms = [] + confidence_maps = [] + + for t in tiles: + + d = os.path.join(t["dir"], "dsm.tif") + if os.path.exists(d): + dsms.append(d) + + c = os.path.join(t["dir"], "confidence.tif") + if os.path.exists(c): + confidence_maps.append(c) + + if dsms: + rasterio.merge.merge(dsms, + bounds=bounds, + res=cfg["dsm_resolution"], + nodata=np.nan, + indexes=[1], + dst_path=os.path.join(cfg["out_dir"], "dsm.tif"), + dst_kwds=creation_options) + + if confidence_maps: + rasterio.merge.merge(confidence_maps, + bounds=bounds, + res=cfg["dsm_resolution"], + nodata=np.nan, + indexes=[1], + dst_path=os.path.join(cfg["out_dir"], "confidence.tif"), + dst_kwds=creation_options) + + +def main(user_cfg, start_from=0): + """ + Launch the s2p pipeline with the parameters given in a json file. + + Args: + user_cfg: user config dictionary + start_from: the step to start from (default: 0) + """ + common.print_elapsed_time.t0 = datetime.datetime.now() + initialization.build_cfg(user_cfg) + initialization.make_dirs() + + # multiprocessing setup + nb_workers = multiprocessing.cpu_count() # nb of available cores + if cfg['max_processes'] is not None: + nb_workers = cfg['max_processes'] + print(f"Running s2p using {nb_workers} workers.") + + tw, th = initialization.adjust_tile_size() + tiles_txt = os.path.join(cfg['out_dir'], 'tiles.txt') + tiles = initialization.tiles_full_info(tw, th, tiles_txt, create_masks=True) + if not tiles: + print('ERROR: the ROI is not seen in two images or is totally masked.') + sys.exit(1) + + if start_from > 0: + assert os.path.exists(tiles_txt), "start_from set to {} but tiles.txt is not found in '{}'. Make sure this is" \ + " the output directory of a previous run.".format(start_from, cfg['out_dir']) + else: + # initialisation: write the list of tilewise json files to outdir/tiles.txt + with open(tiles_txt, 'w') as f: + for t in tiles: + print(t['json'], file=f) + + n = len(cfg['images']) + tiles_pairs = [(t, i) for i in range(1, n) for t in tiles] + timeout = cfg['timeout'] + + # local-pointing step: + if start_from <= 1: + print('1) correcting pointing locally...') + parallel.launch_calls(pointing_correction, tiles_pairs, nb_workers, + timeout=timeout) + check_missing_sift(tiles_pairs) + + # global-pointing step: + if start_from <= 2: + print('2) correcting pointing globally...') + global_pointing_correction(tiles) + common.print_elapsed_time() + + # Create matches GeoJSON. + print("Creating matches GeoJSON") + merge_all_match_files() + matches_to_geojson(f"{cfg['out_dir']}/merged_sift_matches.txt", + cfg['images'][0]['rpcm'], + 10, + [0, 1], + f"{cfg['out_dir']}/matches.geojson" + ) + # rectification step: + if start_from <= 3: + print('3) rectifying tiles...') + parallel.launch_calls(rectification_pair, tiles_pairs, nb_workers, + timeout=timeout) + + # matching step: + if start_from <= 4: + # Check which tiles were rectified correctly, and skip tiles that have missing files + tiles, tiles_pairs = remove_missing_tiles(tiles) + + if cfg['max_processes_stereo_matching'] is not None: + nb_workers_stereo = cfg['max_processes_stereo_matching'] + else: + nb_workers_stereo = nb_workers + try: + + print(f'4) running stereo matching using {nb_workers_stereo} workers...') + parallel.launch_calls(stereo_matching, tiles_pairs, nb_workers_stereo, + timeout=timeout) + except subprocess.CalledProcessError as e: + print(f'ERROR: stereo matching failed. In case this is due too little RAM set ' + f'"max_processes_stereo_matching" to a lower value (currently set to: {nb_workers_stereo}).') + raise e + + if start_from <= 5: + if n > 2: + # disparity-to-height step: + print('5a) computing height maps...') + parallel.launch_calls(disparity_to_height, tiles_pairs, nb_workers, + timeout=timeout) + + print('5b) computing local pairwise height offsets...') + parallel.launch_calls(mean_heights, tiles, nb_workers, timeout=timeout) + + # global-mean-heights step: + print('5c) computing global pairwise height offsets...') + global_mean_heights(tiles) + + # heights-to-ply step: + print('5d) merging height maps and computing point clouds...') + parallel.launch_calls(heights_to_ply, tiles, nb_workers, + timeout=timeout) + else: + # triangulation step: + print('5) triangulating tiles...') + num_tiles = len(tiles) + tiles = parallel.launch_calls(disparity_to_ply, tiles, nb_workers, + timeout=timeout) + tiles = [t for t in tiles if t is not None] + if len(tiles) != num_tiles: + print(f'WARNING: {len(tiles)-num_tiles}/{num_tiles} tiles failed when applying homography ' + 'and will be skipped.') + + # local-dsm-rasterization step: + if start_from <= 6: + print('6) computing DSM by tile...') + parallel.launch_calls(plys_to_dsm, tiles, nb_workers, timeout=timeout) + + # global-dsm-rasterization step: + if start_from <= 7: + print('7) computing global DSM...') + global_dsm(tiles) + + common.print_elapsed_time() + + # cleanup + common.garbage_cleanup() + common.print_elapsed_time(since_first_call=True) + + +def make_path_relative_to_file(path, f): + return os.path.join(os.path.abspath(os.path.dirname(f)), path) + + +def read_tiles(tiles_file): + outdir = os.path.dirname(tiles_file) + + with open(tiles_file) as f: + tiles = f.readlines() + + # Strip trailing \n + tiles = list(map(str.strip, tiles)) + tiles = [os.path.join(outdir, t) for t in tiles] + + return tiles + + +def read_config_file(config_file): + """ + Read a json configuration file and interpret relative paths. + + If any input or output path is a relative path, it is interpreted as + relative to the config_file location (and not relative to the current + working directory). Absolute paths are left unchanged. + """ + with open(config_file, 'r') as f: + user_cfg = json.load(f) + + # output paths + if not os.path.isabs(user_cfg['out_dir']): + user_cfg['out_dir'] = make_path_relative_to_file(user_cfg['out_dir'], + config_file) + + # ROI path + k = "roi_geojson" + if k in user_cfg and isinstance(user_cfg[k], str) and not os.path.isabs(user_cfg[k]): + user_cfg[k] = make_path_relative_to_file(user_cfg[k], config_file) + + if 'exogenous_dem' in user_cfg and user_cfg['exogenous_dem'] is not None: + if not os.path.isabs(user_cfg['exogenous_dem']): + user_cfg['exogenous_dem'] = make_path_relative_to_file(user_cfg['exogenous_dem'], config_file) + + # input paths + for img in user_cfg['images']: + for d in ['img', 'rpc', 'clr', 'cld', 'roi', 'wat']: + if d in img and isinstance(img[d], str) and not os.path.isabs(img[d]): + img[d] = make_path_relative_to_file(img[d], config_file) + + return user_cfg diff --git a/s2p/.ipynb_checkpoints/cli-checkpoint.py b/s2p/.ipynb_checkpoints/cli-checkpoint.py new file mode 100644 index 00000000..a10765bc --- /dev/null +++ b/s2p/.ipynb_checkpoints/cli-checkpoint.py @@ -0,0 +1,29 @@ +import os +import shutil +import argparse + +import s2p + + +def main(): + """ + Command line parsing for s2p command line interface. + """ + parser = argparse.ArgumentParser(description=('S2P: Satellite Stereo ' + 'Pipeline')) + parser.add_argument('config', metavar='config.json', + help=('path to a json file containing the paths to ' + 'input and output files and the algorithm ' + 'parameters')) + parser.add_argument('--start_from', dest='start_from', type=int, + default=0, help="Restart the process from a given step in " + "case of an interruption or to try different parameters.") + args = parser.parse_args() + + user_cfg = s2p.read_config_file(args.config) + + s2p.main(user_cfg, start_from=args.start_from) + + # Backup input file for sanity check + if not args.config.startswith(os.path.abspath(s2p.cfg['out_dir']+os.sep)): + shutil.copy2(args.config,os.path.join(s2p.cfg['out_dir'],'config.json.orig')) diff --git a/s2p/.ipynb_checkpoints/common-checkpoint.py b/s2p/.ipynb_checkpoints/common-checkpoint.py new file mode 100644 index 00000000..cdb1c6ed --- /dev/null +++ b/s2p/.ipynb_checkpoints/common-checkpoint.py @@ -0,0 +1,355 @@ +# Copyright (C) 2015, Carlo de Franchis +# Copyright (C) 2015, Gabriele Facciolo +# Copyright (C) 2015, Enric Meinhardt +# Copyright (C) 2015, Julien Michel + + +import os +import sys +import errno +import datetime +import warnings +import tempfile +import subprocess +import numpy as np +import rasterio + + +from s2p.config import cfg + + +# silent rasterio NotGeoreferencedWarning +warnings.filterwarnings("ignore", + category=rasterio.errors.NotGeoreferencedWarning) + +# add the bin folder to system path +parent_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +bin_dir = os.path.join(parent_dir, 'bin') +os.environ['PATH'] = bin_dir + os.pathsep + os.environ['PATH'] + +# global variable +# list of intermediary files generated by s2p +garbage = list() + + +def remove(target): + try: + os.remove(target) + except OSError: + pass + +def garbage_cleanup(): + """ + Removes all the files listed in the global variable 'garbage'. + """ + if cfg['clean_tmp']: + while garbage: + remove(garbage.pop()) + + +def tmpfile(ext=''): + """ + Creates a temporary file in the cfg['temporary_dir'] directory. + + Args: + ext: desired file extension. The dot has to be included. + + Returns: + absolute path to the created file + + The path of the created file is added to the garbage list to allow cleaning + at the end of the pipeline. + """ + fd, out = tempfile.mkstemp(suffix=ext, prefix='s2p_', + dir=os.path.expandvars(cfg['temporary_dir'])) + os.close(fd) # http://www.logilab.org/blogentry/17873 + garbage.append(out) + return out + + +def run(cmd, env=os.environ, timeout=None, shell=False, raise_errors=True): + """ + Runs a shell command, and print it before running. + + Arguments: + cmd: list of a command and its arguments, or as a fallback, + a string to be passed to a shell that will be split into a list. + env (optional, default value is os.environ): dictionary containing the + environment variables + timeout (optional, int): time in seconds after which the function will + raise an error if the command hasn't returned + + TODO: remove the temporary `shell` argument once all commands use shell=False + shell (bool): run the command in a subshell. Defaults to False. + + Both stdout and stderr of the shell in which the command is run are those + of the parent process. + """ + print("\nRUN: %s" % cmd) + t = datetime.datetime.now() + if not isinstance(cmd, list) and not shell: + cmd = cmd.split() + try: + subprocess.run(cmd, shell=shell, stdout=sys.stdout, stderr=sys.stderr, + env=env, timeout=timeout, check=True) + except subprocess.CalledProcessError as e: + print("ERROR: command failed: %s" % cmd) + print("ERROR: return code: %s" % e.returncode) + print("ERROR: output: %s" % e.output) + if raise_errors: + raise e + return False + print(datetime.datetime.now() - t) + return True + +def mkdir_p(path): + """ + Create a directory without complaining if it already exists. + """ + try: + os.makedirs(path) + except OSError as exc: # requires Python > 2.5 + if exc.errno == errno.EEXIST and os.path.isdir(path): + pass + else: raise + + +def matrix_translation(x, y): + t = np.eye(3) + t[0, 2] = x + t[1, 2] = y + return t + + +def rio_read_as_array_with_nans(im): + """ + Read an image replacing gdal nodata value with np.nan + + Args: + im: path to the input image file + + Returns: + array: raster as numpy array + """ + with rasterio.open(im, 'r') as src: + array = src.read() + nodata_values = src.nodatavals + + for band, nodata in zip(array, nodata_values): + if nodata is not None: + band[band == nodata] = np.nan + + return array.squeeze() + + +def rasterio_write(path, array, profile={}, tags={}): + """ + Write a numpy array in a tiff or png file with rasterio. + + Args: + path (str): path to the output tiff/png file + array (numpy array): 2D or 3D array containing the image to write. + profile (dict): rasterio profile (ie dictionary of metadata) + tags (dict): dictionary with additional geotiff tags + """ + # determine the driver based on the file extension + extension = os.path.splitext(path)[1].lower() + if extension in ['.tif', '.tiff']: + driver = 'GTiff' + elif extension in ['.png']: + driver = 'png' + else: + raise NotImplementedError('format {} not supported'.format(extension)) + + # read image size and number of bands + array = np.atleast_3d(array) + height, width, nbands = array.shape + + # define image metadata dict + profile.update(driver=driver, count=nbands, width=width, height=height, + dtype=array.dtype) + + # write to file + with rasterio.Env(): + with rasterio.open(path, 'w', **profile) as dst: + dst.write(np.transpose(array, (2, 0, 1))) + dst.update_tags(**tags) + + +def image_apply_homography(out, im, H, w, h): + """ + Applies an homography to an image. + + Args: + out: path to the output image file + im: path to the input image file + H: numpy array containing the 3x3 homography matrix + w, h: dimensions (width and height) of the output image + + The output image is defined on the domain [0, w] x [0, h]. Its pixels + intensities are defined by out(x) = im(H^{-1}(x)). + + This function calls the homography binary, rewritten by Marc Lebrun and + Carlo de Franchis based on a code of Pascal Monasse refactored by Gabriele + Facciolo. + """ + # write the matrix to a string + hij = " ".join(str(x) for x in H.flatten()) + + # apply the homography + return run(["homography", im, "-h", hij, out, "%d" % w, "%d" % h], raise_errors=False) + + +def points_apply_homography(H, pts): + """ + Applies an homography to a list of 2D points. + + Args: + H: numpy array containing the 3x3 homography matrix + pts: numpy array containing the list of 2D points, one per line + + Returns: + a numpy array containing the list of transformed points, one per line + """ + # if the list of points is not a numpy array, convert it + if (type(pts) == list): + pts = np.array(pts) + + # convert the input points to homogeneous coordinates + if len(pts[0]) < 2: + raise ValueError( + "The input must be a numpy array" + "of 2D points, one point per line" + ) + pts = np.hstack((pts[:, 0:2], pts[:, 0:1]*0+1)) + + # apply the transformation + Hpts = (np.dot(H, pts.T)).T + + # normalize the homogeneous result and trim the extra dimension + Hpts = Hpts * (1.0 / np.tile( Hpts[:, 2], (3, 1)) ).T + return Hpts[:, 0:2] + + +def bounding_box2D(pts): + """ + bounding box for the points pts + """ + dim = len(pts[0]) # should be 2 + bb_min = [min([t[i] for t in pts]) for i in range(dim)] + bb_max = [max([t[i] for t in pts]) for i in range(dim)] + return bb_min[0], bb_min[1], bb_max[0] - bb_min[0], bb_max[1] - bb_min[1] + + +def crop_array(img, x, y, w, h, fill_value=0): + """ + Crop an image represented as an array. + + Args: + img (array): 2D input image + x, y (ints): coordinate of the top-left corner of the crop + w, h (ints): width and height of the crop + fill_value (img.dtype): constant value used for filling the crop + outside the input image domain + + Returns: + array with the cropped image + """ + crop = fill_value * np.ones((h, w), dtype=img.dtype) + + y0 = max(y, 0) + y1 = min(y + h, img.shape[0]) + x0 = max(x, 0) + x1 = min(x + w, img.shape[1]) + + if y0 < y1 and x0 < x1: # the requested crop overlaps the image + crop[y0 - y:y1 - y, x0 - x:x1 - x] = img[y0:y1, x0:x1] + + return crop + + +def run_binary_on_list_of_points(points, binary, option=None, env_var=None): + """ + Runs a binary that reads its input on stdin. + + Args: + points: numpy array containing all the input points, one per line + binary: path to the binary. It is supposed to write one output value on + stdout for each input point + option: optional option to pass to the binary + env_var (optional): environment variable that modifies the behaviour of + the binary. It is a tuple containing 2 strings, eg ('PATH', '/bin') + + Returns: + a numpy array containing all the output points, one per line. + """ + # send the input points to stdin + pts_file = tmpfile('.txt') + np.savetxt(pts_file, points, '%.18f') + p1 = subprocess.Popen(['cat', pts_file], stdout=subprocess.PIPE) + + # run the binary + env = os.environ.copy() + if env_var is not None: + env[env_var[0]] = env_var[1] + cmd = [binary] + if option is not None: + cmd.append(option) + p2 = subprocess.Popen(cmd, env=env, stdin=p1.stdout, + stdout=subprocess.PIPE) + + # recover output values + out = [] + for i in range(len(points)): + out.append([float(x) for x in p2.stdout.readline().split()]) + + return np.array(out) + + +def cargarse_basura(inputf, outputf): + se=5 + tmp1 = outputf + '1.tif' + tmp2 = outputf + '2.tif' + tmpM = outputf + 'M.tif' + run('morphoop %s min %d %s' % (inputf, se, tmpM)) + run('morphoop %s max %d %s' % (inputf, se, tmp1)) + run('morphoop %s max %d %s' % (inputf, se, tmpM)) + run('morphoop %s min %d %s' % (inputf, se, tmp2)) + run(["plambda", tmp1, tmp2, inputf, "x y - fabs %d > nan z if" % 5, "-o", tmpM]) + run('remove_small_cc %s %s %d %d' % (tmpM, outputf, 200, 5)) + run('rm -f %s %s %s' % (tmp1, tmp2, tmpM)) + + +def print_elapsed_time(since_first_call=False): + """ + Print the elapsed time since the last call or since the first call. + + Args: + since_first_call: + """ + t2 = datetime.datetime.now() + if since_first_call: + print("Total elapsed time:", t2 - print_elapsed_time.t0) + else: + try: + print("Elapsed time:", t2 - print_elapsed_time.t1) + except AttributeError: + print("Elapsed time:", t2 - print_elapsed_time.t0) + print_elapsed_time.t1 = t2 + print() + + +def linear_stretching_and_quantization_8bit(img, p=1): + """ + Simple 8-bit quantization with linear stretching. + + Args: + img (np.array): image to quantize + p (float): percentage of the darkest and brightest pixels to saturate, + from 0 to 100. + + Returns: + numpy array with the quantized uint8 image + """ + a, b = np.nanpercentile(img, (p, 100 - p)) + return np.round(255 * (np.clip(img, a, b) - a) / (b - a)).astype(np.uint8) diff --git a/s2p/.ipynb_checkpoints/config-checkpoint.py b/s2p/.ipynb_checkpoints/config-checkpoint.py new file mode 100644 index 00000000..6df14d70 --- /dev/null +++ b/s2p/.ipynb_checkpoints/config-checkpoint.py @@ -0,0 +1,189 @@ +# Copyright (C) 2013, Carlo de Franchis +# Copyright (C) 2013, Gabriele Facciolo +# Copyright (C) 2013, Enric Meinhardt +# Copyright (C) 2013, Julien Michel + +# This module contains a dictionary, cfg, containing all the parameters of the +# s2p pipeline. This dictionary is updated at runtime with parameters defined +# by the user in the config.json file. All the optional parameters (that the +# user is not forced to define in config.json) must be defined here, otherwise +# they won't have a default value. + +cfg = {} + +# path to output directory +cfg['out_dir'] = "s2p_output" + +# path to directory where (many) temporary files will be stored +cfg['temporary_dir'] = "s2p_tmp" + +# temporary files are erased when s2p terminates. Switch to False to keep them +cfg['clean_tmp'] = True + +# remove all generated files except from ply point clouds and tif raster dsm +cfg['clean_intermediate'] = False + +# switch to True if you want to process the whole image +cfg['full_img'] = False + +# s2p processes the images tile by tile. The tiles are squares cropped from the +# reference image. The width and height of the tiles are given by this param, in pixels. +cfg['tile_size'] = 800 + +# margins used to increase the footprint of the rectified tiles, to +# account for poor disparity estimation close to the borders +cfg['horizontal_margin'] = 50 # for regularity and occlusions +cfg['vertical_margin'] = 10 # for regularity + +# max number of processes launched in parallel. None means the number of available cores +cfg['max_processes'] = None + +# max number of processes launched in parallel for stereo_matching +# Uses the value of cfg['max_processes'] if None +cfg['max_processes_stereo_matching'] = None + +# max number of OMP threads used by programs compiled with openMP +cfg['omp_num_threads'] = 1 + +# timeout in seconds, after which a function that runs on a single tile is not +# waited for +cfg['timeout'] = 600 + +# debug mode (more verbose logs and intermediate results saved) +cfg['debug'] = False + +# resolution of the output digital surface model, in meters per pixel +cfg['dsm_resolution'] = 4 + +# radius to compute altitudes (and to interpolate the small holes) +cfg['dsm_radius'] = 0 + +# dsm_sigma controls the spread of the blob from each point for the dsm computation +# (dsm_resolution by default) +cfg['dsm_sigma'] = None + +# relative sift match threshold (else sift match threshold is absolute) +cfg['relative_sift_match_thresh'] = True + +# if cfg['relative_sift_match_thresh'] is True : +# sift threshold on the first over second best match ratio +# else (absolute) a reasonable value is between 200 and 300 (128-vectors SIFT descriptors) +cfg['sift_match_thresh'] = 0.6 + +# disp range expansion facto +cfg['disp_range_extra_margin'] = 0.2 + +# Maximum disparity range allowed in block matching +cfg['max_disp_range'] = None + +# estimate rectification homographies either blindly using the rpc data or from +# the images actual content thanks to sift matches +cfg['rectification_method'] = 'rpc' # either 'rpc' or 'sift' + +# register the rectified images with a shear estimated from the rpc data +cfg['register_with_shear'] = True + +# number of ground control points per axis in matches from rpc generation +cfg['n_gcp_per_axis'] = 5 + +# max distance allowed for a point to the epipolar line of its match +cfg['epipolar_thresh'] = 0.5 + +# maximal pointing error, in pixels +cfg['max_pointing_error'] = 10 + +# set these params if you want to impose the disparity range manually (cfg['disp_range_method'] == 'fixed_pixel_range') +cfg['disp_min'] = None +cfg['disp_max'] = None + +# set these params if you want to impose the altitude range manually (cfg['disp_range_method'] == 'fixed_altitude_range') +cfg['alt_min'] = None +cfg['alt_max'] = None + +# width of a stripe of pixels to be masked along the reference input image borders +cfg['border_margin'] = 10 + +# radius for erosion of valid disparity areas. Ignored if less than 2 +cfg['msk_erosion'] = 2 + +cfg['fusion_operator'] = 'average_if_close' + +# threshold (in meters) used for the fusion of two dems in triplet processing +cfg['fusion_thresh'] = 3 + +cfg['rpc_alt_range_scale_factor'] = 1 + +# method to compute the disparity range: "sift", "exogenous", "wider_sift_exogenous", "fixed_pixel_range", "fixed_altitude_range" +cfg['disp_range_method'] = "wider_sift_exogenous" +cfg['disp_range_exogenous_low_margin'] = -10 +cfg['disp_range_exogenous_high_margin'] = +100 + +# whether or not to use SRTM DEM (downloaded from internet) to estimate: +# - the average ground altitude (to project the input geographic AOI to the +# correct place in the input images) +# - a reasonable altitude range (to get a better rectification when +# "rectification_method" is set to "rpc") +cfg['use_srtm'] = False + +# exogenous dem. If set, it superseeds SRTM. +cfg['exogenous_dem'] = None +cfg['exogenous_dem_geoid_mode'] = True + +### stereo matching parameters + +# stereo matching algorithm: 'tvl1', 'msmw', 'hirschmuller08', +# hirschmuller08_laplacian', 'sgbm', 'mgm', 'mgm_multi' +cfg['matching_algorithm'] = 'mgm' + +# size of the Census NCC square windows used in mgm +cfg['census_ncc_win'] = 5 + +# MGM parameter: speckle filter minimum area (REMOVESMALLCC flag) +cfg['stereo_speckle_filter'] = 25 + +# MGM parameter: regularity (multiplies P1 and P2) +cfg['stereo_regularity_multiplier'] = 1.0 + +# MGM parameters: +# number of directions explored for regularization +cfg['mgm_nb_directions'] = 8 +# timeout in seconds, after which a running mgm process will be killed +cfg['mgm_timeout'] = 600 +# distance threshold (in pixels) for the left-right consistency test +cfg['mgm_leftright_threshold'] = 1.0 +# controls the mgm left-right consistency check. 0: disabled +# 1 (default): enabled at all scales +# 2: enables only at the last scale (faster) +cfg['mgm_leftright_control'] = 1 +# controls the mgm mindiff filter check. -1: disabled (default), produce denser maps +# 1: enabled, produce conservative results +cfg['mgm_mindiff_control'] = -1 + +# remove isolated 3d points in height maps +cfg['3d_filtering_r'] = None # radius in meters +cfg['3d_filtering_n'] = None # number of points + +# clean height maps outliers +cfg['cargarse_basura'] = True + +# Output coordinate reference system +# All formats accepted by `pyproj.CRS()` are allowed, for example: +# 32740 (int interpreted as an EPSG code), or +# "epsg:32740+5773" (authority string), or +# "+proj=utm +zone=40 +south +datum=WGS84 +units=m +vunits=m +no_defs +type=crs" (proj4 string) +# If None, the local UTM zone will be used +cfg['out_crs'] = None + +# If the out_crs is not set this parameter determines if the output CRS uses the EGM96 geoid vertical datum (if True) +# or the WGS84 ellipsoid vertical datum (if False). If out_crs is set, this parameter is ignored. +cfg['out_geoid'] = False + +# LoFTR parameters, points whose confidence is below this value are ignored. +cfg["confidence_threshold"] = 0.5 + +# Matching parameters +# Matching method to use, choose from ["sift", "loftr", "superglue", "all"] +cfg["matching_method"] = "sift" +# The min and max valeus to scale the image between 0 and 1 with, only used for non-SIFT methods. +cfg["max_value"] = 3000 +cfg["min_value"] = 200 diff --git a/s2p/.ipynb_checkpoints/estimation-checkpoint.py b/s2p/.ipynb_checkpoints/estimation-checkpoint.py new file mode 100644 index 00000000..743b10b4 --- /dev/null +++ b/s2p/.ipynb_checkpoints/estimation-checkpoint.py @@ -0,0 +1,227 @@ +# Copyright (C) 2015, Carlo de Franchis +# Copyright (C) 2015, Gabriele Facciolo +# Copyright (C) 2015, Enric Meinhardt +# Copyright (C) 2015, Julien Michel + +import numpy as np + + +def fundamental_matrix_cameras(P1, P2): + """ + Computes the fundamental matrix given the matrices of two cameras. + + Arguments: + P1, P2: 2D arrays of size 3x4 containing the camera matrices + + Returns: + the computed fundamental matrix, given by the formula 17.3 (p. 412) in + Hartley & Zisserman book (2nd ed.). + """ + X0 = P1[[1, 2], :] + X1 = P1[[2, 0], :] + X2 = P1[[0, 1], :] + Y0 = P2[[1, 2], :] + Y1 = P2[[2, 0], :] + Y2 = P2[[0, 1], :] + + F = np.zeros((3, 3)) + F[0, 0] = np.linalg.det(np.vstack([X0, Y0])) + F[0, 1] = np.linalg.det(np.vstack([X1, Y0])) + F[0, 2] = np.linalg.det(np.vstack([X2, Y0])) + F[1, 0] = np.linalg.det(np.vstack([X0, Y1])) + F[1, 1] = np.linalg.det(np.vstack([X1, Y1])) + F[1, 2] = np.linalg.det(np.vstack([X2, Y1])) + F[2, 0] = np.linalg.det(np.vstack([X0, Y2])) + F[2, 1] = np.linalg.det(np.vstack([X1, Y2])) + F[2, 2] = np.linalg.det(np.vstack([X2, Y2])) + + return F + + +def get_angle_from_cos_and_sin(c, s): + """ + Computes x in ]-pi, pi] such that cos(x) = c and sin(x) = s. + """ + if s >= 0: + return np.arccos(c) + else: + return -np.arccos(c) + + +def rectifying_similarities_from_affine_fundamental_matrix(F, debug=False): + """ + Computes two similarities from an affine fundamental matrix. + + Args: + F: 3x3 numpy array representing the input fundamental matrix + debug (optional, default is False): boolean flag to activate verbose + mode + + Returns: + S, S': two similarities such that, when used to resample the two images + related by the fundamental matrix, the resampled images are + stereo-rectified. + """ + # check that the input matrix is an affine fundamental matrix + assert(np.shape(F) == (3, 3)) + assert(np.linalg.matrix_rank(F) == 2) + np.testing.assert_allclose(F[:2, :2], np.zeros((2, 2))) + + # notations + a = F[0, 2] + b = F[1, 2] + c = F[2, 0] + d = F[2, 1] + e = F[2, 2] + + # rotations + r = np.sqrt(c*c + d*d) + s = np.sqrt(a*a + b*b) + R1 = 1 / r * np.array([[d, -c], [c, d]]) + R2 = 1 / s * np.array([[-b, a], [-a, -b]]) + + # zoom and translation + z = np.sqrt(r / s) + t = 0.5 * e / np.sqrt(r * s) + + if debug: + theta_1 = get_angle_from_cos_and_sin(d, c) + print("reference image:") + print("\trotation: %f deg" % np.rad2deg(theta_1)) + print("\tzoom: %f" % z) + print("\tvertical translation: %f" % t) + print() + theta_2 = get_angle_from_cos_and_sin(-b, -a) + print("secondary image:") + print("\trotation: %f deg" % np.rad2deg(theta_2)) + print("\tzoom: %f" % (1.0 / z)) + print("\tvertical translation: %f" % -t) + + # output similarities + S1 = np.zeros((3, 3)) + S1[0:2, 0:2] = z * R1 + S1[1, 2] = t + S1[2, 2] = 1 + + S2 = np.zeros((3, 3)) + S2[0:2, 0:2] = 1 / z * R2 + S2[1, 2] = -t + S2[2, 2] = 1 + + return S1, S2 + + +def affine_fundamental_matrix(matches): + """ + Estimates the affine fundamental matrix given a set of point correspondences + between two images. + + Arguments: + matches: 2D array of size Nx4 containing a list of pairs of matching + points. Each line is of the form x1, y1, x2, y2, where (x1, y1) is + the point in the first view while (x2, y2) is the matching point in + the second view. + + Returns: + the estimated affine fundamental matrix, given by the Gold Standard + algorithm, as described in Hartley & Zisserman book (see chap. 14). + """ + # revert the order of points to fit H&Z convention (see algo 14.1) + X = matches[:, [2, 3, 0, 1]] + + # compute the centroid + N = len(X) + XX = np.sum(X, axis=0) / N + + # compute the Nx4 matrix A + A = X - np.tile(XX, (N, 1)) + + # the solution is obtained as the singular vector corresponding to the + # smallest singular value of matrix A. See Hartley and Zissermann for + # details. + # It is the last line of matrix V (because np.linalg.svd returns V^T) + U, S, V = np.linalg.svd(A) + N = V[-1, :] + + # extract values and build F + F = np.zeros((3, 3)) + F[0, 2] = N[0] + F[1, 2] = N[1] + F[2, 0] = N[2] + F[2, 1] = N[3] + F[2, 2] = -np.dot(N, XX) + + return F + + +def affine_transformation(x, xx): + """ + Estimates an affine homography from a list of correspondences + + Args: + x: Nx2 numpy array, containing a list of points + xx: Nx2 numpy array, containing the list of corresponding points + + Returns: + 3x3 numpy array, representing (in homogeneous coordinates) an affine + homography that maps the points of x onto the points of xx. + + This function implements the Gold-Standard algorithm for estimating an + affine homography, described in Hartley & Zisserman page 130 (second + edition). + """ + # check that there are at least 3 points + if len(x) < 3: + print("ERROR: estimation.affine_transformation needs 3 correspondences") + return np.eye(3) + + # translate the input points so that the centroid is at the origin. + t = -np.mean(x, 0) + tt = -np.mean(xx, 0) + x = x + t + xx = xx + tt + + # compute the Nx4 matrix A + A = np.hstack((x, xx)) + + # two singular vectors corresponding to the two largest singular values of + # matrix A. See Hartley and Zissermann for details. These are the first + # two lines of matrix V (because np.linalg.svd returns V^T) + U, S, V = np.linalg.svd(A) + v1 = V[0, :] + v2 = V[1, :] + + # compute blocks B and C, then H + tmp = np.vstack((v1, v2)).T + assert(np.shape(tmp) == (4, 2)) + B = tmp[0:2, :] + C = tmp[2:4, :] + H = np.dot(C, np.linalg.inv(B)) + + # return A + A = np.eye(3) + A[0:2, 0:2] = H + A[0:2, 2] = np.dot(H, t) - tt + return A + + +def translation(x, xx): + """ + Estimates a planar translation from a list of correspondences + + Args: + x: Nx2 numpy array, containing a list of points + xx: Nx2 numpy array, containing the list of corresponding points + + Returns: + 3x3 numpy array, representing (in homogeneous coordinates) a planar + translation that maps the points of x onto the points of xx. + """ + # compute the mean of the displacement vectors + t = np.mean(xx - x, 0) + + # return A + A = np.eye(3) + A[0, 2] = t[0] + A[1, 2] = t[1] + return A diff --git a/s2p/.ipynb_checkpoints/image_coordinates_to_coordinates-checkpoint.py b/s2p/.ipynb_checkpoints/image_coordinates_to_coordinates-checkpoint.py new file mode 100644 index 00000000..c65b747f --- /dev/null +++ b/s2p/.ipynb_checkpoints/image_coordinates_to_coordinates-checkpoint.py @@ -0,0 +1,132 @@ +import glob +import os +import rasterio +import numpy as np +import geopandas as gpd +import shapely +from shapely.geometry import Point +import srtm4 +from rpcm import RPCModel, rpc_from_rpc_file +import geopy +import geopy.distance + + +def localize_row_col_geometry( + geometry: shapely.geometry.Point, + rpc_model: RPCModel, +): + """ + Project points in the row and column space to the geographic space, with altitude correction. + All the altitudes here are in meters. + Args: + geometry: The geometry to add altitudes to, should be a shapely point. + rpc_model: the RPCModel to base the projection on. + returns: + A tuple containing (new longitude, new latitude). + """ + points = [] + col_idx, row_idx = geometry.coords[0] + points.append((col_idx, row_idx)) + col_indices, row_indices = np.asarray(points).T + coords = [ + rpc_model.localization(col_idx, row_idx, 0) + for col_idx, row_idx in zip(col_indices, row_indices) + ] + corrected_lons, corrected_lats = [], [] + for (lon, lat), col_idx, row_idx in zip(coords, col_indices, row_indices): + + new_altitude = srtm4.srtm4(lon, lat) + new_lon, new_lat = 0, 0 + max_iterations = 10 + ground_dist_tol = 50 + ground_distance = ground_dist_tol + 1 + iteration_no = 0 + while (ground_distance > ground_dist_tol) and ( + iteration_no < max_iterations + ): + iteration_no += 1 + old_lon, old_lat = new_lon, new_lat + new_lon, new_lat = rpc_model.localization( + col_idx, row_idx, new_altitude + ) + new_altitude = srtm4.srtm4(new_lon, new_lat) + ground_distance = geopy.distance.distance( + (old_lat, old_lon), (new_lat, new_lon) + ).m + corrected_lons.append(new_lon) + corrected_lats.append(new_lat) + points.append( + ((lon, lat) for (lon, lat) in zip(corrected_lons, corrected_lats)) + ) + + corrected_lons.append(new_lon) + corrected_lats.append(new_lat) + points.append( + ((lon, lat) for (lon, lat) in zip(corrected_lons, corrected_lats)) + ) + return (new_lon, new_lat) + + +def read_and_prepare_matches(path_to_all_matches: str, match_columns: list=[0, 1]): + """ + Loads in the merged matches file and returns a geodataframe containing Point-like geometry. + + path_to_all_matches: The path to the text file containing all of the matches. + match_columns: The columns to be used as x, y for the matches geodataframe output. + """ + matches = np.loadtxt(path_to_all_matches) + matches = matches.astype(int) + matches = matches[:, match_columns] + matches_x = matches[:, 0] + matches_y = matches[:, 1] + matches_gdf = gpd.GeoDataFrame(geometry=gpd.points_from_xy(matches_x, matches_y)) + return matches_gdf + + +def correct_points(matches_gdf: gpd.GeoDataFrame, rpc_model_path: str, every_x_match: int=10): + """ + matches_gdf: A geodataframe containing point geometries of the matches created by s2p in image coordinates. + rpc_model_path: A path to the rpc model to be used to calculate real coordinates from image coordinates. + every_x_match: For decreasing the size of the matches geodataframe, takes every point, i.e 5 is every 5. + + return: two numpy arrays, (lons and lats). + """ + corrected_points = [] + if type(rpc_model_path) == str: + rpc_model = rpc_from_rpc_file(rpc_model_path) + else: + rpc_model = rpc_model_path + matches_gdf = matches_gdf["geometry"][::every_x_match] + + # Iterate over matches and correct them. + for i, point in enumerate(matches_gdf): + localized_point = localize_row_col_geometry(point, rpc_model) + corrected_points.append(list(localized_point)) + + # Prepare the lons and lats for returning + corrected_points = np.array(corrected_points) + lons = corrected_points[:, 0] + lats = corrected_points[:, 1] + + return lons, lats + + +def points_to_geojson(out_path: str, lons: np.ndarray, lats: np.ndarray): + """ + Saves the lons and lats from correct_points as a GeoJSON. + + out_path: The path to save the output GeoJSON to, + lons: The longitudes returned by correct_points. + lats: The latitudes returned by correct_points. + """ + gdf_corrected = gpd.GeoDataFrame(geometry=gpd.points_from_xy(lons, lats, crs="EPSG:4326")) + gdf_corrected.to_file(out_path, driver="GeoJSON") + + +def matches_to_geojson(path_to_matches: str, path_to_rpc: str, every_x_match: int, match_columns: list, out_path: str): + matches_gdf = read_and_prepare_matches(path_to_matches, match_columns) + lons, lats = correct_points(matches_gdf, path_to_rpc, every_x_match) + points_to_geojson(out_path, lons, lats) + + + \ No newline at end of file diff --git a/s2p/.ipynb_checkpoints/rectification-checkpoint.py b/s2p/.ipynb_checkpoints/rectification-checkpoint.py new file mode 100644 index 00000000..49aa65ea --- /dev/null +++ b/s2p/.ipynb_checkpoints/rectification-checkpoint.py @@ -0,0 +1,390 @@ +# Copyright (C) 2015, Carlo de Franchis +# Copyright (C) 2015, Gabriele Facciolo +# Copyright (C) 2015, Enric Meinhardt + + +import os +import warnings + +import numpy as np + +from s2p import rpc_utils +from s2p import estimation +from s2p import evaluation +from s2p import common +from s2p import visualisation +from s2p.config import cfg + + +class NoRectificationMatchesError(Exception): + pass + + +class NoHorizontalRegistrationWarning(Warning): + pass + + +def filter_matches_epipolar_constraint(F, matches, thresh): + """ + Discards matches that are not consistent with the epipolar constraint. + + Args: + F: fundamental matrix + matches: list of pairs of 2D points, stored as a Nx4 numpy array + thresh: maximum accepted distance between a point and its matched + epipolar line + + Returns: + the list of matches that satisfy the constraint. It is a sub-list of + the input list. + """ + out = [] + for match in matches: + x = np.array([match[0], match[1], 1]) + xx = np.array([match[2], match[3], 1]) + d1 = evaluation.distance_point_to_line(x, np.dot(F.T, xx)) + d2 = evaluation.distance_point_to_line(xx, np.dot(F, x)) + if max(d1, d2) < thresh: + out.append(match) + + return np.array(out) + + +def register_horizontally_shear(matches, H1, H2): + """ + Adjust rectifying homographies with tilt, shear and translation to reduce the disparity range. + + Args: + matches: list of pairs of 2D points, stored as a Nx4 numpy array + H1, H2: two homographies, stored as numpy 3x3 matrices + + Returns: + H2: corrected homography H2 + + The matches are provided in the original images coordinate system. By + transforming these coordinates with the provided homographies, we obtain + matches whose disparity is only along the x-axis. + """ + # transform the matches according to the homographies + p1 = common.points_apply_homography(H1, matches[:, :2]) + x1 = p1[:, 0] + y1 = p1[:, 1] + p2 = common.points_apply_homography(H2, matches[:, 2:]) + x2 = p2[:, 0] + y2 = p2[:, 1] + + if cfg['debug']: + print("Residual vertical disparities: max, min, mean. Should be zero") + print(np.max(y2 - y1), np.min(y2 - y1), np.mean(y2 - y1)) + + # we search the (a, b, c) vector that minimises \sum (x1 - (a*x2+b*y2+c))^2 + # it is a least squares minimisation problem + A = np.column_stack((x2, y2, y2*0+1)) + a, b, c = np.linalg.lstsq(A, x1, rcond=None)[0] + + # correct H2 with the estimated tilt, shear and translation + return np.dot(np.array([[a, b, c], [0, 1, 0], [0, 0, 1]]), H2) + + +def register_horizontally_translation(matches, H1, H2, flag='center'): + """ + Adjust rectifying homographies with a translation to modify the disparity range. + + Args: + matches: list of pairs of 2D points, stored as a Nx4 numpy array + H1, H2: two homographies, stored as numpy 3x3 matrices + flag: option needed to control how to modify the disparity range: + 'center': move the barycenter of disparities of matches to zero + 'positive': make all the disparities positive + 'negative': make all the disparities negative. Required for + Hirshmuller stereo (java) + + Returns: + H2: corrected homography H2 + + The matches are provided in the original images coordinate system. By + transforming these coordinates with the provided homographies, we obtain + matches whose disparity is only along the x-axis. The second homography H2 + is corrected with a horizontal translation to obtain the desired property + on the disparity range. + """ + # transform the matches according to the homographies + p1 = common.points_apply_homography(H1, matches[:, :2]) + x1 = p1[:, 0] + y1 = p1[:, 1] + p2 = common.points_apply_homography(H2, matches[:, 2:]) + x2 = p2[:, 0] + y2 = p2[:, 1] + + # for debug, print the vertical disparities. Should be zero. + if cfg['debug']: + print("Residual vertical disparities: max, min, mean. Should be zero") + print(np.max(y2 - y1), np.min(y2 - y1), np.mean(y2 - y1)) + + # compute the disparity offset according to selected option + t = 0 + if (flag == 'center'): + t = np.mean(x2 - x1) + if (flag == 'positive'): + t = np.min(x2 - x1) + if (flag == 'negative'): + t = np.max(x2 - x1) + + # correct H2 with a translation + return np.dot(common.matrix_translation(-t, 0), H2) + + +def disparity_range_from_matches(matches, H1, H2, w, h): + """ + Compute the disparity range of a ROI from a list of point matches. + + Args: + matches: Nx4 numpy array containing a list of matches, in the full + image coordinates frame, before rectification + w, h: width and height of the rectangular ROI in the first image. + H1, H2: two rectifying homographies, stored as numpy 3x3 matrices + + Returns: + disp_min, disp_max: horizontal disparity range + """ + # transform the matches according to the homographies + p1 = common.points_apply_homography(H1, matches[:, :2]) + x1 = p1[:, 0] + p2 = common.points_apply_homography(H2, matches[:, 2:]) + x2 = p2[:, 0] + + # compute the final disparity range + disp_min = np.floor(np.min(x2 - x1)) + 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'] + return disp_min, disp_max + + +def disparity_range(rpc1, rpc2, x, y, w, h, H1, H2, matches, A=None): + """ + Compute the disparity range of a ROI from a list of point matches. + + Args: + rpc1, rpc2 (rpcm.RPCModel): two RPC camera models + x, y, w, h (int): 4-tuple of integers defining the rectangular ROI in + the first image. (x, y) is the top-left corner, and (w, h) are the + dimensions of the rectangle. + H1, H2 (np.array): two rectifying homographies, stored as 3x3 arrays + matches (np.array): Nx4 array containing a list of sift matches, in the + full image coordinates frame + A (np.array): 3x3 array containing the pointing error correction for + im2. This matrix is usually estimated with the pointing_accuracy + module. + + Returns: + disp: 2-uple containing the horizontal disparity range + """ + # compute exogenous disparity range if needed + 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, + cfg['disp_range_exogenous_high_margin'], + cfg['disp_range_exogenous_low_margin']) + + print("exogenous disparity range:", exogenous_disp) + + # 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) + else: + sift_disp = None + print("SIFT disparity range:", sift_disp) + + # compute altitude range disparity if needed + if cfg['disp_range_method'] == 'fixed_altitude_range': + alt_disp = rpc_utils.altitude_range_to_disp_range(cfg['alt_min'], + cfg['alt_max'], + rpc1, rpc2, + x, y, w, h, H1, H2, A) + print("disparity range computed from fixed altitude range:", alt_disp) + + # now compute disparity range according to selected method + if cfg['disp_range_method'] == 'exogenous': + disp = exogenous_disp + + elif cfg['disp_range_method'] == 'sift': + disp = sift_disp + + elif cfg['disp_range_method'] == 'wider_sift_exogenous': + if sift_disp is not None and exogenous_disp is not None: + disp = min(exogenous_disp[0], sift_disp[0]), max(exogenous_disp[1], sift_disp[1]) + else: + disp = sift_disp or exogenous_disp + + elif cfg['disp_range_method'] == 'fixed_altitude_range': + disp = alt_disp + + elif cfg['disp_range_method'] == 'fixed_pixel_range': + disp = cfg['disp_min'], cfg['disp_max'] + + # default disparity range to return if everything else broke + if disp is None: + disp = -3, 3 + + # impose a minimal disparity range (TODO this is valid only with the + # 'center' flag for register_horizontally_translation) + disp = min(-3, disp[0]), max(3, disp[1]) + + print("Final disparity range:", disp) + return disp + + +def rectification_homographies(matches, x, y, w, h): + """ + Computes rectifying homographies from point matches for a given ROI. + + The affine fundamental matrix F is estimated with the gold-standard + algorithm, then two rectifying similarities (rotation, zoom, translation) + are computed directly from F. + + Args: + matches: numpy array of shape (n, 4) containing a list of 2D point + correspondences between the two images. + x, y, w, h: four integers defining the rectangular ROI in the first + image. (x, y) is the top-left corner, and (w, h) are the dimensions + of the rectangle. + Returns: + S1, S2, F: three numpy arrays of shape (3, 3) representing the + two rectifying similarities to be applied to the two images and the + corresponding affine fundamental matrix. + """ + # estimate the affine fundamental matrix with the Gold standard algorithm + F = estimation.affine_fundamental_matrix(matches) + + # compute rectifying similarities + S1, S2 = estimation.rectifying_similarities_from_affine_fundamental_matrix(F, cfg['debug']) + + if cfg['debug']: + y1 = common.points_apply_homography(S1, matches[:, :2])[:, 1] + y2 = common.points_apply_homography(S2, matches[:, 2:])[:, 1] + err = np.abs(y1 - y2) + print("max, min, mean rectification error on point matches: ", end=' ') + print(np.max(err), np.min(err), np.mean(err)) + + # pull back top-left corner of the ROI to the origin (plus margin) + pts = common.points_apply_homography(S1, [[x, y], [x+w, y], [x+w, y+h], [x, y+h]]) + x0, y0 = common.bounding_box2D(pts)[:2] + T = common.matrix_translation(-x0, -y0) + 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): + """ + Rectify a ROI in a pair of images. + + Args: + im1, im2: paths to two GeoTIFF image files + rpc1, rpc2: two instances of the rpcm.RPCModel class + x, y, w, h: four integers defining the rectangular ROI in the first + image. (x, y) is the top-left corner, and (w, h) are the dimensions + of the rectangle. + out1, out2: paths to the output rectified crops + A (optional): 3x3 numpy array containing the pointing error correction + for im2. This matrix is usually estimated with the pointing_accuracy + module. + sift_matches (optional): Nx4 numpy array containing a list of sift + matches, in the full image coordinates frame + method (default: 'rpc'): option to decide whether to use rpc of sift + matches for the fundamental matrix estimation. + {h,v}margin (optional): horizontal and vertical margins added on the + sides of the rectified images + + Returns: + H1, H2: Two 3x3 matrices representing the rectifying homographies that + have been applied to the two original (large) images. + disp_min, disp_max: horizontal disparity range + """ + # 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']) + + # correct second image coordinates with the pointing correction matrix + if A is not None: + matches[:, 2:] = common.points_apply_homography(np.linalg.inv(A), + matches[:, 2:]) + elif method == 'sift': + matches = sift_matches + + else: + raise Exception("Unknown value {} for argument 'method'".format(method)) + + + + # 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']) + + # correct second image coordinates with the pointing correction matrix + if A is not None: + matches[:, 2:] = common.points_apply_homography(np.linalg.inv(A), + matches[:, 2:]) + + # compute rectifying homographies + H1, H2, F = rectification_homographies(matches, x, y, w, h) + + 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)) + 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) + + # 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']) + 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) + + # compute disparity range + if cfg['debug']: + out_dir = os.path.dirname(out1) + np.savetxt(os.path.join(out_dir, 'sift_matches_disp.txt'), + sift_matches, fmt='%9.3f') + 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) + + # recompute hmargin and homographies + hmargin = int(np.ceil(max([hmargin, np.fabs(disp_m), np.fabs(disp_M)]))) + T = common.matrix_translation(hmargin, vmargin) + H1, H2 = np.dot(T, H1), np.dot(T, H2) + + # compute output images size + roi = [[x, y], [x+w, y], [x+w, y+h], [x, y+h]] + pts1 = common.points_apply_homography(H1, roi) + x0, y0, w0, h0 = common.bounding_box2D(pts1) + # check that the first homography maps the ROI in the positive quadrant + np.testing.assert_allclose(np.round([x0, y0]), [hmargin, vmargin], atol=.01) + + # apply homographies and do the crops + common.image_apply_homography(out1, im1, H1, w0 + 2*hmargin, h0 + 2*vmargin) + common.image_apply_homography(out2, im2, H2, w0 + 2*hmargin, h0 + 2*vmargin) + + return H1, H2, disp_m, disp_M diff --git a/s2p/.ipynb_checkpoints/sift-checkpoint.py b/s2p/.ipynb_checkpoints/sift-checkpoint.py new file mode 100644 index 00000000..12aa06e2 --- /dev/null +++ b/s2p/.ipynb_checkpoints/sift-checkpoint.py @@ -0,0 +1,378 @@ +# Copyright (C) 2015, Carlo de Franchis +# Copyright (C) 2015, Gabriele Facciolo +# Copyright (C) 2015, Enric Meinhardt +# Copyright (C) 2019, Julien Michel (CNES) + +import os +import ctypes +import warnings + +import numpy as np +import rasterio as rio +from numpy.ctypeslib import ndpointer +import cv2 +import ransac + +from s2p import rpc_utils +from s2p import estimation + +try: + from matching.matcher import get_keypoints_loftr +except ImportError as e: + print("LoFTR matching library could not be imported.") + print(e) + +try: + from matching.matcher import get_keypoints_superglue +except ImportError as e: + print("SuperGlue matching library could not be imported.") + print(e) + + +# Locate sift4ctypes library and raise an ImportError if it can not be +# found This call will raise an exception if library can not be found, +# at import time + +# TODO: This is kind of ugly. Cleaner way to do this is to update +# LD_LIBRARY_PATH, which we should do once we have a proper config file +here = os.path.dirname(os.path.abspath(__file__)) +sift4ctypes = os.path.join(os.path.dirname(here), 'lib', 'libsift4ctypes.so') +lib = ctypes.CDLL(sift4ctypes) + + +# Filter warnings from rasterio reading files wihtout georeferencing +warnings.filterwarnings("ignore", category=rio.errors.NotGeoreferencedWarning) + + +def keypoints_from_nparray(arr, thresh_dog=0.0133, nb_octaves=8, nb_scales=3, offset=None): + """ + Runs SIFT (the keypoints detection and description only, no matching) on an image stored in a 2D numpy array + + It uses Ives Rey Otero's implementation published in IPOL: + http://www.ipol.im/pub/pre/82/ + + Args: + arr: A 2D numpy array respresenting the input image + thresh_dog (optional): Threshold on gaussian derivative + nb_octaves (optional): Number of octaves + nb_scales (optional): Number of scales + offset (optional): offset to apply to sift position in case arr is an extract of a bigger image + + Returns: + A numpy array of shape (nb_points,132) containing for each row (y,x,scale,orientation, sift_descriptor) + """ + # retrieve numpy buffer dimensions + h, w = arr.shape + + # Set expected args and return types + lib.sift.argtypes = (ndpointer(dtype=ctypes.c_float, shape=(h, w)), ctypes.c_uint, ctypes.c_uint, ctypes.c_float, + ctypes.c_uint, ctypes.c_uint, ctypes.POINTER(ctypes.c_uint), ctypes.POINTER(ctypes.c_uint)) + lib.sift.restype = ctypes.POINTER(ctypes.c_float) + + # Create variables to be updated by function call + nb_points = ctypes.c_uint() + desc_size = ctypes.c_uint() + + # Call sift fonction from sift4ctypes.so + keypoints_ptr = lib.sift(arr.astype(np.float32), w, h, thresh_dog, + nb_octaves, nb_scales, ctypes.byref(desc_size), ctypes.byref(nb_points)) + + # Transform result into a numpy array + keypoints = np.asarray([keypoints_ptr[i] + for i in range(nb_points.value*desc_size.value)]) + + # Delete results to release memory + lib.delete_buffer.argtypes = (ctypes.POINTER(ctypes.c_float)), + lib.delete_buffer(keypoints_ptr) + + # Reshape keypoints array + keypoints = keypoints.reshape((nb_points.value, desc_size.value)) + + if offset is not None: + x, y = offset + keypoints[:, 0] += x + keypoints[:, 1] += y + + return keypoints + + +def image_keypoints(im, x, y, w, h, max_nb=None, thresh_dog=0.0133, nb_octaves=8, nb_scales=3): + """ + Runs SIFT (the keypoints detection and description only, no matching). + + It uses Ives Rey Otero's implementation published in IPOL: + http://www.ipol.im/pub/pre/82/ + + Args: + im (str): path to the input image + max_nb (optional): maximal number of keypoints. If more keypoints are + detected, those at smallest scales are discarded + + Returns: + numpy array of shape (n, 132) containing, on each row: (y, x, s, o, 128-descriptor) + """ + # Read file with rasterio + with rio.open(im) as ds: + # clip roi to stay inside the image boundaries + if x < 0: # if x is negative then replace it with 0 and reduce w + w += x + x = 0 + if y < 0: + h += y + y = 0 + # if extract not completely inside the full image then resize (w, h) + w = min(w, ds.width - x) + h = min(h, ds.height - y) + if w <= 0 or h <= 0: + return np.array([]) + in_buffer = ds.read(window=rio.windows.Window(x, y, w, h)) + + # Detect keypoints on first band + keypoints = keypoints_from_nparray(in_buffer[0], thresh_dog=thresh_dog, + nb_octaves=nb_octaves, + nb_scales=nb_scales, offset=(x, y)) + + # Limit number of keypoints if needed + if max_nb is not None: + keypoints = keypoints[:max_nb] + + return keypoints + + +def string_dump_of_keypoint_and_descriptor(k): + """ + Return a string representing a keypoint and its descriptor. + + Args: + k (array_like): list of 132 floats, the first four elements are the + keypoint (x, y, scale, orientation), the 128 following elements are + the coefficients of the SIFT descriptor and take integer values + between 0 and 255. + + Return: + string dump of the descriptor, such as for example + "342.254 003.570 0.91346 2.36788 000 001 005 000 000 000 028 029 179..." + """ + s = "{:8.3f} {:8.3f} {:7.3f} {: 5.3f} ".format(*k[:4]) + s += " ".join("{:3d}".format(int(x)) for x in k[4:]) + return s + + +def keypoints_match(k1, k2, method='relative', sift_thresh=0.6, F=None, + epipolar_threshold=10, model=None, ransac_max_err=0.3): + """ + Find matches among two lists of sift keypoints. + + Args: + k1 (array): numpy array of shape (n, 132), where each row represents a + sift keypoint with (y, x, scale, orientation, 128-descriptor) + k2 (array): numpy array of shape (m, 132), where each row represents a + sift keypoint + method (optional, default is 'relative'): flag ('relative' or + 'absolute') indicating whether to use absolute distance or relative + distance + sift_thresh (optional, default is 0.6): threshold for distance between SIFT + descriptors. These descriptors are 128-vectors, whose coefficients + range from 0 to 255, thus with absolute distance a reasonable value + for this threshold is between 200 and 300. With relative distance + (ie ratio between distance to nearest and distance to second + nearest), the commonly used value for the threshold is 0.6. + F (optional): affine fundamental matrix + epipolar_threshold (optional, default is 10): maximum distance allowed for + a point to the epipolar line of its match. + model (optional, default is None): model imposed by RANSAC when + searching the set of inliers. If None all matches are considered as + inliers. + ransac_max_err (float): maximum allowed epipolar error for + RANSAC inliers. Optional, default is 0.3. + + Returns: + if any, a numpy 2D array containing the list of inliers matches. + """ + # compute matches + matches = keypoints_match_from_nparray(k1, k2, method, sift_thresh, + epipolar_threshold, F) + + # filter matches with ransac + if model == 'fundamental' and len(matches) >= 7: + inliers = ransac.find_fundamental_matrix(matches, ntrials=1000, + max_err=ransac_max_err)[0] + matches = matches[inliers] + + return matches + + +def keypoints_match_from_nparray(k1, k2, method, sift_threshold, + epi_threshold=10, F=None): + """ + Wrapper for the sift keypoints matching function of libsift4ctypes.so. + """ + # Set expected args and return types + lib.matching.argtypes = (ndpointer(dtype=ctypes.c_float, shape=k1.shape), + ndpointer(dtype=ctypes.c_float, shape=k2.shape), + ctypes.c_uint, ctypes.c_uint, ctypes.c_uint, + ctypes.c_uint, ctypes.c_float, ctypes.c_float, + ndpointer(dtype=ctypes.c_double, shape=(5,)), + ctypes.c_bool, ctypes.c_bool, + ctypes.POINTER(ctypes.c_uint)) + lib.matching.restype = ctypes.POINTER(ctypes.c_float) + + # Get info of descriptor size + nb_sift_k1, descr = k1.shape + sift_offset = 4 + length_descr = descr - sift_offset + + # Transform information of method into boolean + use_relative_method = (method == 'relative') + + # Format fundamental matrix + use_fundamental_matrix = False + coeff_mat = np.zeros(5) + if F is not None: + coeff_mat = np.asarray([F[0, 2], F[1, 2], F[2, 0], F[2, 1], F[2, 2]]) + use_fundamental_matrix = True + + # Create variables to be updated by function call + nb_matches = ctypes.c_uint() + + # Call sift fonction from sift4ctypes.so + matches_ptr = lib.matching(k1.astype('float32'), k2.astype('float32'), + length_descr, sift_offset, len(k1), len(k2), + sift_threshold, epi_threshold, coeff_mat, + use_fundamental_matrix, use_relative_method, + ctypes.byref(nb_matches)) + + # Transform result into a numpy array + matches = np.asarray([matches_ptr[i] for i in range(nb_matches.value * 4)]) + + # Delete results to release memory + lib.delete_buffer.argtypes = ctypes.POINTER(ctypes.c_float), + lib.delete_buffer(matches_ptr) + + # Reshape keypoints array + return matches.reshape((nb_matches.value, 4)) + + +def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h, + method, sift_thresh, epipolar_threshold, 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. + The corresponding roi in the second image is determined using the rpc + functions. + Args: + im1, im2: paths to two large tif images + rpc1, rpc2: two instances of the rpcm.RPCModel class + x, y, w, h: four integers defining the rectangular ROI in the first + image. (x, y) is the top-left corner, and (w, h) are the dimensions + of the rectangle. + method, sift_thresh, epipolar_threshold: see docstring of + s2p.sift.keypoints_match() + Returns: + matches: 2D numpy array containing a list of matches. Each line + 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) + + # estimate an approximate affine fundamental matrix from the rpcs + rpc_matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, 5) + F = estimation.affine_fundamental_matrix(rpc_matches) + + if matching_method == "loftr": + p1, p2 = get_keypoints_loftr(im1, im2, min_value, max_value, confidence_threshold, + x, x2, y, y2, w, w2, h, h2, rpc_match=True) + if len(p1) < 10 or len(p2) < 10: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + return None + p1[:, 0] += x + p1[:, 1] += y + p2[:, 0] += x2 + p2[:, 1] += y2 + matches = np.hstack((p1, p2)) + + if len(matches) == 0 or matches.ndim != 2: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + return None + elif len(matches) > 10: + inliers = ransac.find_fundamental_matrix(matches, ntrials=1000, + max_err=0.3)[0] + matches = matches[inliers] + else: + print("WARNING: sift.matches_on_rpc_roi: not enough matches.") + return None + + + elif matching_method == "superglue": + p1, p2 = get_keypoints_superglue(im1, im2, min_value, max_value, x, x2, y, y2, w, w2, h, h2, rpc_match=True) + if len(p1) < 10 or len(p2) < 10: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + return None + + p1[:, 0] += x + p1[:, 1] += y + p2[:, 0] += x2 + p2[:, 1] += y2 + matches = np.hstack((p1, p2)) + + if matches.ndim != 2: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + return None + + elif matching_method == "sift": + # if less than 10 matches, lower thresh_dog. An alternative would be ASIFT + thresh_dog = 0.0133 + for _ in range(2): + p1 = image_keypoints(im1, x, y, w, h, thresh_dog=thresh_dog) + p2 = image_keypoints(im2, x2, y2, w2, h2, thresh_dog=thresh_dog) + matches = keypoints_match(p1, p2, method, sift_thresh, F, + epipolar_threshold=epipolar_threshold, + model='fundamental') + if matches is not None and matches.ndim == 2 and matches.shape[0] > 10: + break + thresh_dog /= 2.0 + else: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + return None + + + elif matching_method == "all": + p1, p2 = get_keypoints_loftr(im1, im2, min_value, max_value, confidence_threshold, + x, x2, y, y2, w, w2, h, h2, rpc_match=True) + p1[:, 0] += x + p1[:, 1] += y + p2[:, 0] += x2 + p2[:, 1] += y2 + matches = np.hstack((p1, p2)) + if len(matches) == 0 or matches.ndim != 2: + print("WARNING: sift.matches_on_rpc_roi: found no matches.") + matches = None + elif len(matches) > 10: + inliers = ransac.find_fundamental_matrix(matches, ntrials=1000, + max_err=0.3)[0] + matches = matches[inliers] + + p1_sg, p2_sg = get_keypoints_superglue(im1, im2, min_value, max_value, x, x2, y, y2, w, w2, h, h2, rpc_match=True) + p1_sg[:, 0] += x + p1_sg[:, 1] += y + p1_sg[:, 0] += x2 + p1_sg[:, 1] += y2 + matches_sg = np.hstack((p1_sg, p2_sg)) + + thresh_dog = 0.0133 + for _ in range(2): + p1_sift = image_keypoints(im1, x, y, w, h, thresh_dog=thresh_dog) + p2_sift = image_keypoints(im2, x2, y2, w2, h2, thresh_dog=thresh_dog) + matches_sift = keypoints_match(p1_sift, p2_sift, method, sift_thresh, F, + epipolar_threshold=epipolar_threshold, + model='fundamental') + if matches_sift is not None and matches_sift.ndim == 2 and matches_sift.shape[0] > 10: + break + thresh_dog /= 2.0 + + matches = np.vstack((matches, matches_sg, matches_sift)) + + if len(matches) < 10: + return None + + return matches diff --git a/s2p/aggregation.py b/s2p/aggregation.py new file mode 100644 index 00000000..b5c99f19 --- /dev/null +++ b/s2p/aggregation.py @@ -0,0 +1,116 @@ +import tensorflow as tf +import tensorflow.keras as keras + + +L2 = 1.0e-5 +alpha = 0.2 + + +def conv3d(filters, kernel_size, strides, padding): + return keras.layers.Conv3D(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + + +def conv3d_bn(filters, kernel_size, strides, padding, activation): + conv = keras.layers.Conv3D(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + use_bias=False, kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + bn = keras.layers.BatchNormalization() + leaky_relu = keras.layers.LeakyReLU(alpha=alpha) + + if activation: + return keras.Sequential([conv, bn, leaky_relu]) + else: + return keras.Sequential([conv, bn]) + + +def trans_conv3d_bn(filters, kernel_size, strides, padding, activation): + conv = keras.layers.Conv3DTranspose(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + use_bias=False, kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + bn = keras.layers.BatchNormalization() + leaky_relu = keras.layers.LeakyReLU(alpha=alpha) + + if activation: + return keras.Sequential([conv, bn, leaky_relu]) + else: + return keras.Sequential([conv, bn]) + + +class Hourglass(keras.Model): + def __init__(self, filters): + super(Hourglass, self).__init__() + + self.conv1 = conv3d_bn(filters, 3, 1, 'same', True) + self.conv2 = conv3d_bn(filters, 3, 1, 'same', True) + self.conv3 = conv3d_bn(2 * filters, 3, 2, 'same', True) + self.conv4 = conv3d_bn(2 * filters, 3, 1, 'same', True) + self.conv5 = conv3d_bn(2 * filters, 3, 2, 'same', True) + self.conv6 = conv3d_bn(2 * filters, 3, 1, 'same', True) + self.conv7 = trans_conv3d_bn(2 * filters, 4, 2, 'same', True) + self.conv8 = trans_conv3d_bn(filters, 4, 2, 'same', True) + + def call(self, inputs, training=None, mask=None): + x1 = self.conv1(inputs) + x1 = self.conv2(x1) + x2 = self.conv3(x1) + x2 = self.conv4(x2) + x3 = self.conv5(x2) + x3 = self.conv6(x3) + x4 = self.conv7(x3) + x4 += x2 + x5 = self.conv8(x4) + x5 += x1 + + return x5 # [N, D, H, W, C] + + +class FeatureFusion(keras.Model): + def __init__(self, units): + super(FeatureFusion, self).__init__() + + self.upsample = keras.layers.UpSampling3D(size=(2, 2, 2)) + self.add1 = keras.layers.Add() + self.avg_pool3d = keras.layers.GlobalAvgPool3D() + self.fc1 = keras.layers.Dense(units=units, use_bias=True, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + self.relu = keras.layers.Activation(activation='relu') + self.fc2 = keras.layers.Dense(units=units, use_bias=True, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + self.sigmoid = keras.layers.Activation(activation='sigmoid') + self.multiply1 = keras.layers.Multiply() + self.multiply2 = keras.layers.Multiply() + self.add2 = keras.layers.Add() + + def call(self, inputs, training=None, mask=None): + # inputs[0]: lower, inputs[1]: higher + assert len(inputs) == 2 + + x1 = self.upsample(inputs[0]) + x2 = self.add1([x1, inputs[1]]) + shape = x2.get_shape().as_list() + v = self.avg_pool3d(x2) + v = self.fc1(v) + v = self.relu(v) + v = self.fc2(v) + v = self.sigmoid(v) + v1 = 1.0 - v + v = tf.expand_dims(v, 1) + v = tf.expand_dims(v, 1) + v = tf.expand_dims(v, 1) + v = tf.tile(v, [1, shape[1], shape[2], shape[3], 1]) + v1 = tf.expand_dims(v1, 1) + v1 = tf.expand_dims(v1, 1) + v1 = tf.expand_dims(v1, 1) + v1 = tf.tile(v1, [1, shape[1], shape[2], shape[3], 1]) + x3 = self.multiply1([x1, v]) + x4 = self.multiply2([inputs[1], v1]) + x = self.add2([x3, x4]) + + return x diff --git a/s2p/block_matching.py b/s2p/block_matching.py index b61baa73..df5828da 100644 --- a/s2p/block_matching.py +++ b/s2p/block_matching.py @@ -8,6 +8,7 @@ import rasterio from s2p import common +from s2p.hmsmnet import HMSMNet from s2p.config import cfg @@ -87,6 +88,13 @@ def compute_disparity_map(im1, im2, disp, mask, algo, disp_min=None, env = os.environ.copy() env['OMP_NUM_THREADS'] = str(cfg['omp_num_threads']) + if algo == 'hmsmnet': + model = HMSMNet(1024, 1024, 1, -128, 64) + model.build_model() + model.load_weights("s2p/HMSM-Net.h5") + model.predict(im1, im2, disp) + create_rejection_mask(disp, im1, im2, mask) + # call the block_matching binary if algo == 'hirschmuller02': bm_binary = 'subpix.sh' diff --git a/s2p/computation.py b/s2p/computation.py new file mode 100644 index 00000000..ce72d808 --- /dev/null +++ b/s2p/computation.py @@ -0,0 +1,23 @@ +import tensorflow as tf +import tensorflow.keras as keras + + +class Estimation(keras.Model): + def __init__(self, min_disp=-112.0, max_disp=16.0): + super(Estimation, self).__init__() + self.min_disp = int(min_disp) + self.max_disp = int(max_disp) + self.conv = keras.layers.Conv3D(filters=1, kernel_size=3, + strides=1, padding='same', + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(1.0e-5)) + + def call(self, inputs, training=None, mask=None): + x = self.conv(inputs) # [N, D, H, W, 1] + x = tf.squeeze(x, -1) # [N, D, H, W] + x = tf.transpose(x, (0, 2, 3, 1)) # [N, H, W, D] + assert x.shape[-1] == self.max_disp - self.min_disp + candidates = tf.linspace(1.0 * self.min_disp, 1.0 * self.max_disp - 1.0, self.max_disp - self.min_disp) + probabilities = tf.math.softmax(-1.0 * x, -1) + disparities = tf.reduce_sum(candidates * probabilities, -1, True) + return disparities diff --git a/s2p/cost.py b/s2p/cost.py new file mode 100644 index 00000000..4b5037d5 --- /dev/null +++ b/s2p/cost.py @@ -0,0 +1,50 @@ +import tensorflow as tf +import tensorflow.keras as keras + + +class CostConcatenation(keras.Model): + def __init__(self, min_disp=-112.0, max_disp=16.0): + super(CostConcatenation, self).__init__() + self.min_disp = int(min_disp) + self.max_disp = int(max_disp) + + def call(self, inputs, training=None, mask=None): + assert len(inputs) == 2 + cost_volume = [] + for i in range(self.min_disp, self.max_disp): + if i < 0: + cost_volume.append(tf.pad( + tensor=tf.concat([inputs[0][:, :, :i, :], inputs[1][:, :, -i:, :]], -1), + paddings=[[0, 0], [0, 0], [0, -i], [0, 0]], mode='CONSTANT')) + elif i > 0: + cost_volume.append(tf.pad( + tensor=tf.concat([inputs[0][:, :, i:, :], inputs[1][:, :, :-i, :]], -1), + paddings=[[0, 0], [0, 0], [i, 0], [0, 0]], mode='CONSTANT')) + else: + cost_volume.append(tf.concat([inputs[0], inputs[1]], -1)) + cost_volume = tf.stack(cost_volume, 1) + return cost_volume + + +class CostDifference(keras.Model): + def __init__(self, min_disp=-112.0, max_disp=16.0): + super(CostDifference, self).__init__() + self.min_disp = int(min_disp) + self.max_disp = int(max_disp) + + def call(self, inputs, training=None, mask=None): + assert len(inputs) == 2 + cost_volume = [] + for i in range(self.min_disp, self.max_disp): + if i < 0: + cost_volume.append(tf.pad( + tensor=inputs[0][:, :, :i, :] - inputs[1][:, :, -i:, :], + paddings=[[0, 0], [0, 0], [0, -i], [0, 0]], mode='CONSTANT')) + elif i > 0: + cost_volume.append(tf.pad( + tensor=inputs[0][:, :, i:, :] - inputs[1][:, :, :-i, :], + paddings=[[0, 0], [0, 0], [i, 0], [0, 0]], mode='CONSTANT')) + else: + cost_volume.append(inputs[0] - inputs[1]) + cost_volume = tf.stack(cost_volume, 1) + return cost_volume diff --git a/s2p/data_reader.py b/s2p/data_reader.py new file mode 100644 index 00000000..26c2c099 --- /dev/null +++ b/s2p/data_reader.py @@ -0,0 +1,78 @@ +import cv2 +import random +import numpy as np +import scipy.signal as sig + + +kx = np.array([[-1, 0, 1]]) +ky = np.array([[-1], [0], [1]]) + + +def read_left(filename): + img = cv2.imread(filename, -1) + original_shape = img.shape + img = cv2.resize(img, (1024, 1024)) + img = (img - np.nanmean(img)) / np.nanstd(img) + dx = sig.convolve2d(img, kx, 'same') + dy = sig.convolve2d(img, ky, 'same') + img = np.expand_dims(img.astype('float32'), -1) + dx = np.expand_dims(dx.astype('float32'), -1) + dy = np.expand_dims(dy.astype('float32'), -1) + return img, dx, dy, original_shape + + +def read_right(filename): + img = cv2.imread(filename, -1) + original_shape = img.shape + img = cv2.resize(img, (1024, 1024)) + img = (img - np.nanmean(img)) / np.nanstd(img) + return np.expand_dims(img.astype('float32'), -1), original_shape + + +def read_disp(filename): + disp = cv2.imread(filename, cv2.IMREAD_UNCHANGED) + disp_16x = cv2.resize(disp, (64, 64)) / 16.0 + disp_8x = cv2.resize(disp, (128, 128)) / 8.0 + disp_4x = cv2.resize(disp, (256, 256)) / 4.0 + disp = np.expand_dims(disp, -1) + disp_16x = np.expand_dims(disp_16x, -1) + disp_8x = np.expand_dims(disp_8x, -1) + disp_4x = np.expand_dims(disp_4x, -1) + return disp_16x, disp_8x, disp_4x, disp + + +def read_batch(left_paths, right_paths, disp_paths): + lefts, dxs, dys, rights, d16s, d8s, d4s, ds = [], [], [], [], [], [], [], [] + for left_path, right_path, disp_path in zip(left_paths, right_paths, disp_paths): + left, dx, dy = read_left(left_path) + right = read_right(right_path) + d16, d8, d4, d = read_disp(disp_path) + lefts.append(left) + dxs.append(dx) + dys.append(dy) + rights.append(right) + d16s.append(d16) + d8s.append(d8) + d4s.append(d4) + ds.append(d) + return np.array(lefts), np.array(rights), np.array(dxs), np.array(dys),\ + np.array(d16s), np.array(d8s), np.array(d4s), np.array(ds) + + +def load_batch(all_left_paths, all_right_paths, all_disp_paths, batch_size=4, reshuffle=False): + assert len(all_left_paths) == len(all_disp_paths) + assert len(all_right_paths) == len(all_disp_paths) + + i = 0 + while True: + lefts, rights, dxs, dys, d16s, d8s, d4s, ds = read_batch( + left_paths=all_left_paths[i * batch_size:(i + 1) * batch_size], + right_paths=all_right_paths[i * batch_size:(i + 1) * batch_size], + disp_paths=all_disp_paths[i * batch_size:(i + 1) * batch_size]) + yield [lefts, rights, dxs, dys], [d16s, d8s, d4s, ds] + i = (i + 1) % (len(all_left_paths) // batch_size) + if reshuffle: + if i == 0: + paths = list(zip(all_left_paths, all_right_paths, all_disp_paths)) + random.shuffle(paths) + all_left_paths, all_right_paths, all_disp_paths = zip(*paths) diff --git a/s2p/evaluation_disparity.py b/s2p/evaluation_disparity.py new file mode 100644 index 00000000..5da89070 --- /dev/null +++ b/s2p/evaluation_disparity.py @@ -0,0 +1,70 @@ +import glob +import cv2 +import numpy as np + + +def compute_epe(est_path, gt_path, min_disp, max_disp): + est = cv2.imread(est_path, cv2.IMREAD_UNCHANGED) + gt = cv2.imread(gt_path, cv2.IMREAD_UNCHANGED) + + zeros = np.zeros_like(gt, 'int32') + ones = np.ones_like(gt, 'int32') + mask1 = np.where(gt >= max_disp, zeros, ones) + mask2 = np.where(gt < min_disp, zeros, ones) + mask = mask1 & mask2 + + error = np.sum(np.abs(est - gt) * mask) + nums = np.sum(mask) + epe = error / nums + + return error, nums, epe + + +def compute_d1(est_path, gt_path, min_disp, max_disp): + est = cv2.imread(est_path, cv2.IMREAD_UNCHANGED) + gt = cv2.imread(gt_path, cv2.IMREAD_UNCHANGED) + + zeros = np.zeros_like(gt, 'int32') + ones = np.ones_like(gt, 'int32') + mask1 = np.where(gt >= max_disp, zeros, ones) + mask2 = np.where(gt < min_disp, zeros, ones) + mask = mask1 & mask2 + + err_map = np.abs(est - gt) * mask + err_mask = err_map > 3 + err_disps = np.sum(err_mask.astype('float32')) + nums = np.sum(mask) + d1 = err_disps / nums + + return err_disps, nums, d1 + + +def evaluate(est_path, gt_path, min_disp, max_disp): + error, nums, epe = compute_epe(est_path, gt_path, min_disp, max_disp) + print('Sum of absolute error: %f, num of valid pixels: %d, end-point-error: %f' + % (error, int(nums), epe)) + err_disps, nums, d1 = compute_d1(est_path, gt_path, min_disp, max_disp) + print('Num of error disparities: %d, num of valid pixels: %d, d1: %f' + % (int(err_disps), int(nums), d1)) + + +def evaluate_all(est_dir, gt_dir, min_disp, max_disp): + est_paths = glob.glob(est_dir + '/*') + gt_paths = glob.glob(gt_dir + '/*') + est_paths.sort() + gt_paths.sort() + assert len(est_paths) == len(gt_paths) + + total_error, total_nums = 0, 0 + for est_path, gt_path in zip(est_paths, gt_paths): + error, nums, epe = compute_epe(est_path, gt_path, min_disp, max_disp) + total_error += error + total_nums += nums + print('\nEnd-point-error: %f' % (total_error / total_nums)) + + total_err_disps, total_nums = 0, 0 + for est_path, gt_path in zip(est_paths, gt_paths): + err_disps, nums, d1 = compute_d1(est_path, gt_path, min_disp, max_disp) + total_err_disps += err_disps + total_nums += nums + print('\nD1: %f' % (total_err_disps / total_nums)) diff --git a/s2p/feature.py b/s2p/feature.py new file mode 100644 index 00000000..ae5e8b0d --- /dev/null +++ b/s2p/feature.py @@ -0,0 +1,93 @@ +import tensorflow.keras as keras + + +L2 = 1.0e-5 + + +def conv2d(filters, kernel_size, strides, padding, dilation_rate): + return keras.layers.Conv2D(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + dilation_rate=dilation_rate, use_bias=True, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + + +def conv2d_bn(filters, kernel_size, strides, padding, dilation_rate, activation): + conv = keras.layers.Conv2D(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + dilation_rate=dilation_rate, use_bias=False, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + bn = keras.layers.BatchNormalization() + relu = keras.layers.ReLU() + + if activation: + return keras.Sequential([conv, bn, relu]) + else: + return keras.Sequential([conv, bn]) + + +def avg_pool(pool_size, filters): + pool = keras.layers.AvgPool2D(pool_size=pool_size) + conv = keras.layers.Conv2D(filters, kernel_size=1, strides=1, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + + return keras.Sequential([pool, conv]) + + +class BasicBlock(keras.Model): + def __init__(self, filters, dilation_rate): + super(BasicBlock, self).__init__() + + self.conv1 = conv2d_bn(filters, 3, 1, 'same', dilation_rate, True) + self.conv2 = conv2d_bn(filters, 3, 1, 'same', dilation_rate, False) + self.relu = keras.layers.ReLU() + + def call(self, inputs, training=None, mask=None): + x = self.conv1(inputs) + x = self.conv2(x) + x += inputs + x = self.relu(x) + + return x + + +def make_blocks(filters, dilation_rate, num): + blocks = keras.Sequential() + for i in range(num): + blocks.add(BasicBlock(filters, dilation_rate)) + + return blocks + + +class FeatureExtraction(keras.Model): + def __init__(self, filters): + super(FeatureExtraction, self).__init__() + + self.conv0_1 = conv2d_bn(filters, 5, 2, 'same', 1, True) + self.conv0_2 = conv2d_bn(2 * filters, 5, 2, 'same', 1, True) + + self.conv1_0 = make_blocks(2 * filters, 1, 4) + self.conv1_1 = make_blocks(2 * filters, 2, 2) + self.conv1_2 = make_blocks(2 * filters, 4, 2) + self.conv1_3 = make_blocks(2 * filters, 1, 2) + + self.branch0 = avg_pool(1, filters) + self.branch1 = avg_pool(2, filters) + self.branch2 = avg_pool(4, filters) + + def call(self, inputs, training=None, mask=None): + x = self.conv0_1(inputs) + x = self.conv0_2(x) + + x = self.conv1_0(x) + x = self.conv1_1(x) + x = self.conv1_2(x) + x = self.conv1_3(x) + + x0 = self.branch0(x) + x1 = self.branch1(x) + x2 = self.branch2(x) + + return [x0, x1, x2] # [1/4, 1/8, 1/16] diff --git a/s2p/hmsmnet.py b/s2p/hmsmnet.py new file mode 100644 index 00000000..6ea463a8 --- /dev/null +++ b/s2p/hmsmnet.py @@ -0,0 +1,119 @@ +import os +import glob +import time +import numpy as np +import tensorflow.keras as keras +from PIL import Image +from s2p.feature import FeatureExtraction +from s2p.cost import CostConcatenation +from s2p.aggregation import Hourglass, FeatureFusion +from s2p.computation import Estimation +from s2p.refinement import Refinement +from s2p.data_reader import read_left, read_right +from s2p.evaluation_disparity import evaluate_all +import cv2 + + +class HMSMNet: + def __init__(self, height, width, channel, min_disp, max_disp): + self.height = height + self.width = width + self.channel = channel + self.min_disp = int(min_disp) + self.max_disp = int(max_disp) + self.model = None + + def build_model(self): + left_image = keras.Input(shape=(self.height, self.width, self.channel)) + right_image = keras.Input(shape=(self.height, self.width, self.channel)) + gx = keras.Input(shape=(self.height, self.width, self.channel)) + gy = keras.Input(shape=(self.height, self.width, self.channel)) + + feature_extraction = FeatureExtraction(filters=16) + [l0, l1, l2] = feature_extraction(left_image) + [r0, r1, r2] = feature_extraction(right_image) + + cost0 = CostConcatenation(min_disp=self.min_disp // 4, max_disp=self.max_disp // 4) + cost1 = CostConcatenation(min_disp=self.min_disp // 8, max_disp=self.max_disp // 8) + cost2 = CostConcatenation(min_disp=self.min_disp // 16, max_disp=self.max_disp // 16) + cost_volume0 = cost0([l0, r0]) + cost_volume1 = cost1([l1, r1]) + cost_volume2 = cost2([l2, r2]) + + hourglass0 = Hourglass(filters=16) + hourglass1 = Hourglass(filters=16) + hourglass2 = Hourglass(filters=16) + agg_cost0 = hourglass0(cost_volume0) + agg_cost1 = hourglass1(cost_volume1) + agg_cost2 = hourglass2(cost_volume2) + + estimator2 = Estimation(min_disp=self.min_disp // 16, max_disp=self.max_disp // 16) + disparity2 = estimator2(agg_cost2) + + fusion1 = FeatureFusion(units=16) + fusion_cost1 = fusion1([agg_cost2, agg_cost1]) + hourglass3 = Hourglass(filters=16) + agg_fusion_cost1 = hourglass3(fusion_cost1) + + estimator1 = Estimation(min_disp=self.min_disp // 8, max_disp=self.max_disp // 8) + disparity1 = estimator1(agg_fusion_cost1) + + fusion2 = FeatureFusion(units=16) + fusion_cost2 = fusion2([agg_fusion_cost1, agg_cost0]) + hourglass4 = Hourglass(filters=16) + agg_fusion_cost2 = hourglass4(fusion_cost2) + + estimator0 = Estimation(min_disp=self.min_disp // 4, max_disp=self.max_disp // 4) + disparity0 = estimator0(agg_fusion_cost2) + + # refinement + refiner = Refinement(filters=32) + final_disp = refiner([disparity0, left_image, gx, gy]) + + self.model = keras.Model(inputs=[left_image, right_image, gx, gy], + outputs=[disparity2, disparity1, disparity0, final_disp]) + #self.model.summary() + + def load_weights(self, weights): + self.model.load_weights(weights) + + def predict(self, left_path, right_path, output_path): + left_image, gx, gy, left_shape = read_left(left_path) + left_image = np.nan_to_num(left_image) + right_image, right_shape = read_right(right_path) + right_image = np.nan_to_num(right_image) + left_image = np.expand_dims(left_image, 0) + gx = np.expand_dims(gx, 0) + gy = np.expand_dims(gy, 0) + right_image = np.expand_dims(right_image, 0) + disparity_increase = (left_shape[0] * left_shape[1]) / (1024*1024) + disparity = self.model.predict([left_image, right_image, gx, gy]) + disparity = disparity[-1][0, :, :, 0] + disparity *= disparity_increase + disparity[disparity == 0] = np.nan + disparity = cv2.resize(disparity, (left_shape[1], left_shape[0])) + disparity = Image.fromarray(disparity) + disparity.save(output_path) + + +if __name__ == '__main__': + # predict + weights = 'HMSMNet/HMSM-Net.h5' + net = HMSMNet(1024, 1024, 1, -128, 64) + folder = "/mnt/d/overstory/disparity_experiments/area_4" + net.build_model() + net.load_weights(weights) + + for tile in os.listdir(f"{folder}/tiles"): + for col in os.listdir(f"{folder}/tiles/{tile}"): + left_1 = f"{folder}/tiles/{tile}/{col}/pair_1/rectified_ref.tif" + left_2 = f"{folder}/tiles/{tile}/{col}/pair_2/rectified_ref.tif" + right_1 = f"{folder}/tiles/{tile}/{col}/pair_1/rectified_sec.tif" + right_2 = f"{folder}/tiles/{tile}/{col}/pair_2/rectified_sec.tif" + out_1 = f"{folder}/tiles/{tile}/{col}/pair_1/rectified_disp.tif" + out_2 = f"{folder}/tiles/{tile}/{col}/pair_2/rectified_disp.tif" + + + net.predict(right_1, left_1, out_1) + net.predict(right_2, left_2, out_2) + diff --git a/s2p/refinement.py b/s2p/refinement.py new file mode 100644 index 00000000..3eeed8a9 --- /dev/null +++ b/s2p/refinement.py @@ -0,0 +1,45 @@ +import tensorflow as tf +import tensorflow.keras as keras +from s2p.feature import conv2d, L2 + + +def conv_bn_act(filters, kernel_size, strides, padding, dilation_rate): + conv = keras.layers.Conv2D(filters=filters, kernel_size=kernel_size, + strides=strides, padding=padding, + dilation_rate=dilation_rate, use_bias=False, + kernel_initializer='he_normal', + kernel_regularizer=keras.regularizers.l2(L2)) + bn = keras.layers.BatchNormalization() + act = keras.layers.LeakyReLU() + + return keras.Sequential([conv, bn, act]) + + +class Refinement(keras.Model): + def __init__(self, filters): + super(Refinement, self).__init__() + self.conv1 = conv_bn_act(filters, 3, 1, 'same', 1) + self.conv2 = conv_bn_act(filters, 3, 1, 'same', 1) + self.conv3 = conv_bn_act(filters, 3, 1, 'same', 2) + self.conv4 = conv_bn_act(filters, 3, 1, 'same', 3) + self.conv5 = conv_bn_act(filters, 3, 1, 'same', 1) + self.conv6 = conv2d(1, 3, 1, 'same', 1) + + def call(self, inputs, training=None, mask=None): + # inputs: [disparity, rgb, gx, gy] + assert len(inputs) == 4 + + scale_factor = inputs[1].shape[1] / inputs[0].shape[1] + disp = tf.image.resize(inputs[0], [inputs[1].shape[1], inputs[1].shape[2]]) + disp = disp * scale_factor + + concat = tf.concat([disp, inputs[1], inputs[2], inputs[3]], -1) + delta = self.conv1(concat) + delta = self.conv2(delta) + delta = self.conv3(delta) + delta = self.conv4(delta) + delta = self.conv5(delta) + delta = self.conv6(delta) + disp_final = disp + delta + + return disp_final diff --git a/setup.py b/setup.py index b018fd86..254b3d47 100644 --- a/setup.py +++ b/setup.py @@ -56,7 +56,8 @@ def finalize_options(self): 'requests', 'opencv-python', 'geopandas', - 'geopy'] + 'geopy', + 'tensorflow'] extras_require = { "test": ["pytest", "pytest-cov", "psutil"], diff --git a/tests/data/input_pair/config.json b/tests/data/input_pair/config.json index ac0a8d6b..e0623e35 100644 --- a/tests/data/input_pair/config.json +++ b/tests/data/input_pair/config.json @@ -7,15 +7,18 @@ "roi" : { "x" : 150, "y" : 150, - "w" : 700, - "h" : 700 + "w" : 1000, + "h" : 1000 }, "horizontal_margin": 20, "vertical_margin": 5, - "tile_size" : 300, + "tile_size" : 1000, + "matching_algorithm": "hmsmnet", "disp_range_method" : "sift", "msk_erosion": 0, "dsm_resolution": 0.5, "3d_filtering_r": 5, - "3d_filtering_n": 50 + "3d_filtering_n": 50, + "debug": true, + "clean_intermediate": false } diff --git a/tests/data/input_triplet/config.json b/tests/data/input_triplet/config.json index 94624d6f..52c1bf90 100644 --- a/tests/data/input_triplet/config.json +++ b/tests/data/input_triplet/config.json @@ -1,20 +1,21 @@ { - "out_dir" : "../../testoutput/output_triplet", + "out_dir" : "testoutput/output_triplet_hmsmnet", "images" : [ - {"img" : "img_02.tif"}, {"img" : "img_01.tif"}, {"img" : "img_03.tif"} ], "roi" : { "x" : 150, "y" : 150, - "w" : 700, - "h" : 700 + "w" : 1024, + "h" : 1024 }, "horizontal_margin": 20, "vertical_margin": 5, - "tile_size" : 300, + "tile_size" : 1024, "disp_range_method" : "sift", "msk_erosion": 0, - "dsm_resolution": 0.5 + "dsm_resolution": 0.5, + "matching_algorithm": "hmsmnet", + "debug": true }