From d204ab77916b111386cb5676a0c478af9a54bbbf Mon Sep 17 00:00:00 2001 From: Piero Toffanin Date: Mon, 18 Jun 2018 09:54:09 -0400 Subject: [PATCH] Skirts improvements, DEM default improvements --- opendm/dem/commands.py | 4 ++++ opendm/mesh.py | 27 ++++++++++++++------------- scripts/odm_dem.py | 4 ++-- 3 files changed, 20 insertions(+), 15 deletions(-) diff --git a/opendm/dem/commands.py b/opendm/dem/commands.py index 071188d8..a53ce55c 100644 --- a/opendm/dem/commands.py +++ b/opendm/dem/commands.py @@ -126,6 +126,10 @@ def gap_fill(filenames, fout): return_distances=False, return_indices=True) arr = arr[tuple(indices)] + + # Median filter + from scipy import signal + arr = signal.medfilt(arr, 5) # write output imgout = gippy.GeoImage.create_from(imgs[0], fout) diff --git a/opendm/mesh.py b/opendm/mesh.py index 230b81de..8239109f 100644 --- a/opendm/mesh.py +++ b/opendm/mesh.py @@ -1,11 +1,12 @@ from __future__ import absolute_import -import os, shutil, sys, struct +import os, shutil, sys, struct, random from io import BytesIO from gippy import GeoImage from opendm.dem import commands from opendm import system from opendm import log -from scipy import signal +from scipy import signal, ndimage +import numpy as np def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples=1, verbose=False): # Create DSM from point cloud @@ -18,7 +19,7 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples= os.mkdir(tmp_directory) log.ODM_INFO('Created temporary directory: %s' % tmp_directory) - radius_steps = [dsm_resolution / 4.0, dsm_resolution / 2.0, dsm_resolution, dsm_resolution * 2] + radius_steps = [dsm_resolution / 4.0, dsm_resolution / 2.0, dsm_resolution] log.ODM_INFO('Creating DSM for 2.5D mesh') @@ -29,6 +30,7 @@ def create_25dmesh(inPointCloud, outMesh, dsm_resolution=0.05, depth=8, samples= gapfill=True, outdir=tmp_directory, resolution=dsm_resolution, + products=['idw'], verbose=verbose ) @@ -48,11 +50,6 @@ def dem_to_points(inGeotiff, outPointCloud): arr = image['z'].read_raw() resolution = max(abs(image.resolution().x()), abs(image.resolution().y())) - # Median filter - log.ODM_INFO('Applying median filter...') - - arr = signal.medfilt2d(arr, 5) - log.ODM_INFO('Writing points...') mem_file = BytesIO() @@ -64,7 +61,7 @@ def dem_to_points(inGeotiff, outPointCloud): skirt_points = 0 skirt_height_threshold = 1 # meter - skirt_increments = resolution + skirt_increments = 0.1 for x in range(2, arr_width - 2): for y in range(2, arr_height - 2): @@ -74,16 +71,20 @@ def dem_to_points(inGeotiff, outPointCloud): mem_file.write(struct.pack('ffffff', tx, ty, z, 0, 0, 1)) # Skirting - for (nx, ny) in ((y + 1, x), (y - 1, x), (y, x + 1), (y, x - 1)): + for (nx, ny) in ((x, y + 1), (x, y - 1), (x + 1, y), (x - 1, y)): current_z = z - neighbor_z = arr[nx][ny] + neighbor_z = arr[ny][nx] + if current_z - neighbor_z > skirt_height_threshold: - stop_at_z = neighbor_z + skirt_increments - while current_z > stop_at_z: + while current_z > neighbor_z: current_z -= skirt_increments mem_file.write(struct.pack('ffffff', tx, ty, current_z, 0, 0, 1)) skirt_points += 1 + mem_file.write(struct.pack('ffffff', tx, ty, neighbor_z, 0, 0, 1)) + skirt_points += 1 + + with open(outPointCloud, "wb") as f: f.write("ply\n") f.write("format binary_%s_endian 1.0\n" % sys.byteorder) diff --git a/scripts/odm_dem.py b/scripts/odm_dem.py index ec039eeb..08864b6a 100644 --- a/scripts/odm_dem.py +++ b/scripts/odm_dem.py @@ -90,9 +90,9 @@ class ODMDEMCell(ecto.Cell): if args.dsm: products.append('dsm') if args.dtm: products.append('dtm') - radius_steps = [args.dem_resolution] + radius_steps = [args.dem_resolution / 4.0] for _ in range(args.dem_gapfill_steps - 1): - radius_steps.append(radius_steps[-1] * 3) # 3 is arbitrary, maybe there's a better value? + radius_steps.append(radius_steps[-1] * 2) # 2 is arbitrary, maybe there's a better value? for product in products: commands.create_dems(