From 673a25b79ca63e5868cce798c28314e192910ba8 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Fri, 21 Jan 2022 19:39:52 +0000 Subject: [PATCH] Added visibility test to orthorectification tool --- contrib/orthorectify/README.md | 13 ++++-- contrib/orthorectify/orthorectify.py | 68 ++++++++++++++++++++++------ 2 files changed, 63 insertions(+), 18 deletions(-) diff --git a/contrib/orthorectify/README.md b/contrib/orthorectify/README.md index 3a74c69c..f239793b 100644 --- a/contrib/orthorectify/README.md +++ b/contrib/orthorectify/README.md @@ -2,7 +2,7 @@ ![image](https://user-images.githubusercontent.com/1951843/111536715-fc91c380-8740-11eb-844c-5b7960186391.png) -This tool is capable of orthorectifying individual images (or all images) from an ODM reconstruction. It does not account for visibility occlusion, so you will get artifacts near buildings (help us improve this?). +This tool is capable of orthorectifying individual images (or all images) from an existing ODM reconstruction. ![image](https://user-images.githubusercontent.com/1951843/111529183-3ad6b500-8738-11eb-9960-b1aa676f863b.png) @@ -17,7 +17,7 @@ docker run -ti --rm -v /home/youruser/datasets:/datasets opendronemap/odm --proj You can run the orthorectification module by running: ``` -docker run -ti --rm -v /home/youruser/datasets:/datasets --entrypoint /code/contrib/orthorectify/orthorectify.py opendronemap/odm /datasets/project +docker run -ti --rm -v /home/youruser/datasets:/datasets --entrypoint /code/contrib/orthorectify/run.sh opendronemap/odm /datasets/project ``` This will start the orthorectification process for all images in the dataset. See additional flags you can pass at the end of the command above: @@ -26,7 +26,8 @@ This will start the orthorectification process for all images in the dataset. Se usage: orthorectify.py [-h] [--dem DEM] [--no-alpha NO_ALPHA] [--interpolation {nearest,bilinear}] [--outdir OUTDIR] [--image-list IMAGE_LIST] - [--images IMAGES] + [--images IMAGES] [--threads THREADS] + [--skip-visibility-test SKIP_VISIBILITY_TEST] dataset Orthorectification Tool @@ -53,6 +54,9 @@ optional arguments: --images IMAGES Comma-separated list of filenames to rectify. Use as an alternative to --image- list. Default: process all images. + --skip-visibility-test SKIP_VISIBILITY_TEST + Skip visibility testing (faster but leaves + artifacts due to relief displacement) ``` ## Roadmap @@ -60,5 +64,6 @@ optional arguments: Help us improve this module! We could add: - [ ] GPU support for faster processing - - [ ] Visibility checks - [ ] Merging of multiple orthorectified images (blending, filtering, seam leveling) + - [ ] Faster visibility test + - [ ] Different methods for orthorectification (direct) \ No newline at end of file diff --git a/contrib/orthorectify/orthorectify.py b/contrib/orthorectify/orthorectify.py index 7b488dd8..f09a4483 100755 --- a/contrib/orthorectify/orthorectify.py +++ b/contrib/orthorectify/orthorectify.py @@ -6,12 +6,14 @@ import os import sys sys.path.insert(0, os.path.join("..", "..", os.path.dirname(__file__))) +from math import sqrt import rasterio import numpy as np import numpy.ma as ma import multiprocessing import argparse import functools +from skimage.draw import line from opensfm import dataset default_dem_path = "odm_dem/dsm.tif" @@ -50,7 +52,9 @@ parser.add_argument('--threads', type=int, default=multiprocessing.cpu_count(), help="Number of CPU processes to use. Default: %(default)s") - +parser.add_argument('--skip-visibility-test', + type=bool, + help="Skip visibility testing (faster but leaves artifacts due to relief displacement)") args = parser.parse_args() dataset_path = args.dataset @@ -112,11 +116,16 @@ with rasterio.open(dem_path) as dem_raster: dem_has_nodata = dem_raster.profile.get('nodata') is not None if dem_has_nodata: - dem_min_value = ma.array(dem, mask=dem==dem_raster.nodata).min() + m = ma.array(dem, mask=dem==dem_raster.nodata) + dem_min_value = m.min() + dem_max_value = m.max() else: dem_min_value = dem.min() - + dem_max_value = dem.max() + print("DEM Minimum: %s" % dem_min_value) + print("DEM Maximum: %s" % dem_max_value) + h, w = dem.shape crs = dem_raster.profile.get('crs') @@ -132,18 +141,18 @@ with rasterio.open(dem_path) as dem_raster: exit(1) with open(coords_file) as f: - line = f.readline() # discard + l = f.readline() # discard # second line is a northing/easting offset - line = f.readline().rstrip() - dem_offset_x, dem_offset_y = map(float, line.split(" ")) + l = f.readline().rstrip() + dem_offset_x, dem_offset_y = map(float, l.split(" ")) print("DEM offset: (%s, %s)" % (dem_offset_x, dem_offset_y)) print("DEM dimensions: %sx%s pixels" % (w, h)) - + # Read reconstruction - udata = dataset.UndistortedDataSet(dataset.DataSet(os.path.join(dataset_path, "opensfm"))) + udata = dataset.UndistortedDataSet(dataset.DataSet(os.path.join(dataset_path, "opensfm")), undistorted_data_path=os.path.join(dataset_path, "opensfm", "undistorted")) reconstructions = udata.load_undistorted_reconstruction() if len(reconstructions) == 0: raise Exception("No reconstructions available") @@ -160,6 +169,8 @@ with rasterio.open(dem_path) as dem_raster: r = shot.pose.get_rotation_matrix() Xs, Ys, Zs = shot.pose.get_origin() + cam_grid_y, cam_grid_x = dem_raster.index(Xs + dem_offset_x, Ys + dem_offset_y) + a1 = r[0][0] b1 = r[0][1] c1 = r[0][2] @@ -170,15 +181,26 @@ with rasterio.open(dem_path) as dem_raster: b3 = r[2][1] c3 = r[2][2] + if not args.skip_visibility_test: + distance_map = np.full((h, w), np.nan) + + for j in range(0, h): + for i in range(0, w): + distance_map[j][i] = sqrt((cam_grid_x - i) ** 2 + (cam_grid_y - j) ** 2) + distance_map[distance_map==0] = 1e-7 + print("Camera pose: (%f, %f, %f)" % (Xs, Ys, Zs)) img_h, img_w, num_bands = shot_image.shape + half_img_w = (img_w - 1) / 2.0 + half_img_h = (img_h - 1) / 2.0 print("Image dimensions: %sx%s pixels" % (img_w, img_h)) f = shot.camera.focal * max(img_h, img_w) has_nodata = dem_raster.profile.get('nodata') is not None def process_pixels(step): imgout = np.full((num_bands, dem_bbox_h, dem_bbox_w), np.nan) + minx = dem_bbox_w miny = dem_bbox_h maxx = 0 @@ -192,13 +214,14 @@ with rasterio.open(dem_path) as dem_raster: im_i = i - dem_bbox_minx # World coordinates - Xa, Ya = dem_raster.xy(j, i) Za = dem[j][i] # Skip nodata if has_nodata and Za == dem_raster.nodata: continue + Xa, Ya = dem_raster.xy(j, i) + # Remove offset (our cameras don't have the geographic offset) Xa -= dem_offset_x Ya -= dem_offset_y @@ -209,11 +232,27 @@ with rasterio.open(dem_path) as dem_raster: dz = (Za - Zs) den = a3 * dx + b3 * dy + c3 * dz - x = (img_w - 1) / 2.0 - (f * (a1 * dx + b1 * dy + c1 * dz) / den) - y = (img_h - 1) / 2.0 - (f * (a2 * dx + b2 * dy + c2 * dz) / den) + x = half_img_w - (f * (a1 * dx + b1 * dy + c1 * dz) / den) + y = half_img_h - (f * (a2 * dx + b2 * dy + c2 * dz) / den) if x >= 0 and y >= 0 and x <= img_w - 1 and y <= img_h - 1: - + # Visibility test + if not args.skip_visibility_test: + check_dem_points = np.column_stack(line(i, j, cam_grid_x, cam_grid_y)) + check_dem_points = check_dem_points[np.all(np.logical_and(np.array([0, 0]) <= check_dem_points, check_dem_points < [w, h]), axis=1)] + + visible = True + for p in check_dem_points: + ray_z = Zs + (distance_map[p[1]][p[0]] / distance_map[j][i]) * dz + if ray_z > dem_max_value: + break + + if dem[p[1]][p[0]] > ray_z: + visible = False + break + if not visible: + continue + if interpolation == 'bilinear': xi = img_w - 1 - x yi = img_h - 1 - y @@ -291,6 +330,7 @@ with rasterio.open(dem_path) as dem_raster: # Merge image imgout, _ = results[0] + for j in range(dem_bbox_miny, dem_bbox_maxy + 1): im_j = j - dem_bbox_miny resimg, _ = results[j % max_workers] @@ -308,10 +348,10 @@ with rasterio.open(dem_path) as dem_raster: miny = min(bounds[1], miny) maxx = max(bounds[2], maxx) maxy = max(bounds[3], maxy) - + print("Output bounds: (%s, %s), (%s, %s) pixels" % (minx, miny, maxx, maxy)) if minx <= maxx and miny <= maxy: - imgout = imgout[:,miny:maxy,minx:maxx] + imgout = imgout[:,miny:maxy+1,minx:maxx+1] if with_alpha: alpha = np.zeros((imgout.shape[1], imgout.shape[2]), dtype=np.uint8)