From 97711f5d893dab8ac587da46fa1b53bcea3cf4f7 Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Wed, 17 Mar 2021 20:55:31 +0000 Subject: [PATCH] Add README for rectification module --- contrib/orthorectify/README.md | 60 ++++++++++++++++++++++++++++ contrib/orthorectify/orthorectify.py | 40 ++++++++++++++++--- 2 files changed, 95 insertions(+), 5 deletions(-) create mode 100644 contrib/orthorectify/README.md mode change 100644 => 100755 contrib/orthorectify/orthorectify.py diff --git a/contrib/orthorectify/README.md b/contrib/orthorectify/README.md new file mode 100644 index 00000000..548749bd --- /dev/null +++ b/contrib/orthorectify/README.md @@ -0,0 +1,60 @@ +# Orthorectification Tool + +This tool is capable of orthorectifing 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?). + +## Usage + +After running a reconstruction using ODM: + +``` +docker run -ti --rm -v /home/youruser/datasets:/datasets opendronemap/odm --project-path /datasets project +``` + +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 +``` + +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: + +``` +usage: orthorectify.py [-h] [--dem DEM] [--no-alpha NO_ALPHA] + [--interpolation {nearest,bilinear}] + [--outdir OUTDIR] [--image-list IMAGE_LIST] + [--images IMAGES] + dataset + +Orthorectification Tool + +positional arguments: + dataset Path to ODM dataset + +optional arguments: + -h, --help show this help message and exit + --dem DEM Absolute path to DEM to use to + orthorectify images. Default: + odm_dem/dsm.tif + --no-alpha NO_ALPHA Don't output an alpha channel + --interpolation {nearest,bilinear} + Type of interpolation to use to sample + pixel values.Default: bilinear + --outdir OUTDIR Output directory where to store results. + Default: orthorectified + --image-list IMAGE_LIST + Path to file that contains the list of + image filenames to orthorectify. By + default all images in a dataset are + processed. Default: img_list.txt + --images IMAGES Comma-separeted list of filenames to + rectify. Use as an alternative to --image- + list. Default: process all images. +``` + +## Roadmap + +Help us improve this module! We could add: + + - [ ] GPU support for faster processing + - [ ] Visibility checks + - [ ] Merging of multiple orthorectified images (blending, filtering, seam leveling) \ No newline at end of file diff --git a/contrib/orthorectify/orthorectify.py b/contrib/orthorectify/orthorectify.py old mode 100644 new mode 100755 index 6da8c2f1..fe2f41fd --- a/contrib/orthorectify/orthorectify.py +++ b/contrib/orthorectify/orthorectify.py @@ -31,7 +31,7 @@ parser.add_argument('--interpolation', type=str, choices=('nearest', 'bilinear'), default='bilinear', - help="Type of interpolation to use to sample pixel values. Can be one of: %(choice)s. Default: %(default)s") + help="Type of interpolation to use to sample pixel values.Default: %(default)s") parser.add_argument('--outdir', type=str, default=default_outdir, @@ -43,7 +43,7 @@ parser.add_argument('--image-list', parser.add_argument('--images', type=str, default="", - help="Comma-separeted list of filenames to rectify. Use as an alternative to --image-list. Default: %(default)s") + help="Comma-separeted list of filenames to rectify. Use as an alternative to --image-list. Default: process all images.") args = parser.parse_args() @@ -105,6 +105,27 @@ with rasterio.open(dem_path) as dem_raster: dem = dem_raster.read()[0] h, w = dem.shape + crs = dem_raster.profile.get('crs') + dem_offset_x, dem_offset_y = (0, 0) + + if crs: + print("DEM has a CRS: %s" % str(crs)) + + # Read coords.txt + coords_file = os.path.join(dataset_path, "odm_georeferencing", "coords.txt") + if not os.path.exists(coords_file): + print("Whoops! Cannot find %s (we need that!)" % coords_file) + exit(1) + + with open(coords_file) as f: + line = f.readline() # discard + + # second line is a northing/easting offset + line = f.readline().rstrip() + dem_offset_x, dem_offset_y = map(float, line.split(" ")) + + print("DEM offset: (%s, %s)" % (dem_offset_x, dem_offset_y)) + print("DEM dimensions: %sx%s pixels" % (w, h)) # Read reconstruction @@ -131,6 +152,7 @@ with rasterio.open(dem_path) as dem_raster: img_h, img_w, num_bands = shot_image.shape 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, h, w), np.nan) @@ -146,6 +168,14 @@ with rasterio.open(dem_path) as dem_raster: Xa, Ya = dem_raster.xy(j, i) Za = dem[j][i] + # Skip nodata + if has_nodata and Za == dem_raster.nodata: + continue + + # Remove offset (our cameras don't have the geographic offset) + Xa -= dem_offset_x + Ya -= dem_offset_y + # Colinearity function http://web.pdx.edu/~jduh/courses/geog493f14/Week03.pdf dx = (Xa - Xs) dy = (Ya - Ys) @@ -157,13 +187,12 @@ with rasterio.open(dem_path) as dem_raster: if x >= 0 and y >= 0 and x <= img_w - 1 and y <= img_h - 1: - # for b in range(num_bands): if interpolation == 'bilinear': xi = img_w - 1 - x yi = img_h - 1 - y values = bilinear_interpolate(shot_image, xi, yi) else: - # Linear + # nearest xi = img_w - 1 - int(round(x)) yi = img_h - 1 - int(round(y)) values = shot_image[yi][xi] @@ -231,7 +260,8 @@ with rasterio.open(dem_path) as dem_raster: 'dtype': imgout.dtype.name, 'transform': rasterio.transform.Affine(dem_transform[0], dem_transform[1], offset_x, dem_transform[3], dem_transform[4], offset_y), - 'nodata': None + 'nodata': None, + 'crs': crs } outfile = os.path.join(cwd_path, shot.id)