Merge pull request #1355 from pierotofy/extent

Feat: user-defined boundary
pull/1358/head
Piero Toffanin 2021-10-13 15:53:38 -04:00 zatwierdzone przez GitHub
commit bd3069ee2b
Nie znaleziono w bazie danych klucza dla tego podpisu
ID klucza GPG: 4AEE18F83AFDEB23
16 zmienionych plików z 246 dodań i 111 usunięć

Wyświetl plik

@ -1 +1 @@
2.6.4
2.6.5

118
opendm/boundary.py 100644
Wyświetl plik

@ -0,0 +1,118 @@
import fiona
import fiona.crs
import os
import io
import json
from opendm import system
from pyproj import CRS
from opendm.location import transformer
from opendm.utils import double_quote
from osgeo import ogr
from opendm.shots import get_origin
def compute_boundary_from_shots(reconstruction_json, buffer=0, reconstruction_offset=(0, 0)):
if not os.path.isfile(reconstruction_json):
raise IOError(reconstruction_json + " does not exist.")
with open(reconstruction_json) as f:
data = json.load(f)
reconstruction = data[0]
mp = ogr.Geometry(ogr.wkbMultiPoint)
for shot_image in reconstruction['shots']:
shot = reconstruction['shots'][shot_image]
if shot['gps_dop'] < 999999:
camera = reconstruction['cameras'][shot['camera']]
p = ogr.Geometry(ogr.wkbPoint)
origin = get_origin(shot)
p.AddPoint_2D(origin[0] + reconstruction_offset[0], origin[1] + reconstruction_offset[1])
mp.AddGeometry(p)
if mp.GetGeometryCount() < 3:
return None
convexhull = mp.ConvexHull()
boundary = convexhull.Buffer(buffer)
return load_boundary(boundary.ExportToJson())
def load_boundary(boundary_json, reproject_to_proj4=None):
if not isinstance(boundary_json, str):
boundary_json = json.dumps(boundary_json)
with fiona.open(io.BytesIO(boundary_json.encode('utf-8')), 'r') as src:
if len(src) != 1:
raise IOError("Boundary must have a single polygon (found: %s)" % len(src))
geom = src[0]['geometry']
if geom['type'] != 'Polygon':
raise IOError("Boundary must have a polygon feature (found: %s)" % geom['type'])
rings = geom['coordinates']
if len(rings) == 0:
raise IOError("Boundary geometry has no rings")
coords = rings[0]
if len(coords) == 0:
raise IOError("Boundary geometry has no coordinates")
dimensions = len(coords[0])
if reproject_to_proj4 is not None:
t = transformer(CRS.from_proj4(fiona.crs.to_string(src.crs)),
CRS.from_proj4(reproject_to_proj4))
coords = [t.TransformPoint(*c)[:dimensions] for c in coords]
return coords
def boundary_offset(boundary, reconstruction_offset):
if boundary is None or reconstruction_offset is None:
return boundary
res = []
dims = len(boundary[0])
for c in boundary:
if dims == 2:
res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1]))
else:
res.append((c[0] - reconstruction_offset[0], c[1] - reconstruction_offset[1], c[2]))
return res
def as_polygon(boundary):
if boundary is None:
return None
return "POLYGON((" + ", ".join([" ".join(map(str, c)) for c in boundary]) + "))"
def export_to_bounds_files(boundary, proj4, bounds_json_file, bounds_gpkg_file):
with open(bounds_json_file, "w") as f:
f.write(json.dumps({
"type": "FeatureCollection",
"name": "bounds",
"features": [{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [boundary]
}
}]
}))
if os.path.isfile(bounds_gpkg_file):
os.remove(bounds_gpkg_file)
kwargs = {
'proj4': proj4,
'input': double_quote(bounds_json_file),
'output': double_quote(bounds_gpkg_file)
}
system.run('ogr2ogr -overwrite -f GPKG -a_srs "{proj4}" {output} {input}'.format(**kwargs))

Wyświetl plik

@ -303,6 +303,24 @@ def config(argv=None, parser=None):
'Use 0 to disable cropping. '
'Default: %(default)s'))
parser.add_argument('--boundary',
default='',
metavar='<json>',
action=StoreValue,
type=path_or_json_string,
help='GeoJSON polygon limiting the area of the reconstruction. '
'Can be specified either as path to a GeoJSON file or as a '
'JSON string representing the contents of a '
'GeoJSON file. Default: %(default)s')
parser.add_argument('--auto-boundary',
action=StoreTrue,
nargs=0,
default=False,
help='Automatically set a boundary using camera shot locations to limit the area of the reconstruction. '
'This can help remove far away background artifacts (sky, background landscapes, etc.). See also --boundary. '
'Default: %(default)s')
parser.add_argument('--pc-quality',
metavar='<string>',
action=StoreValue,

Wyświetl plik

@ -5,6 +5,7 @@ from opendm.point_cloud import export_summary_json
from osgeo import ogr
import json, os
from opendm.concurrency import get_max_memory
from opendm.utils import double_quote
class Cropper:
def __init__(self, storage_dir, files_prefix = "crop"):
@ -41,9 +42,9 @@ class Cropper:
try:
kwargs = {
'gpkg_path': gpkg_path,
'geotiffInput': original_geotiff,
'geotiffOutput': geotiff_path,
'gpkg_path': double_quote(gpkg_path),
'geotiffInput': double_quote(original_geotiff),
'geotiffOutput': double_quote(geotiff_path),
'options': ' '.join(map(lambda k: '-co {}={}'.format(k, gdal_options[k]), gdal_options)),
'warpOptions': ' '.join(warp_options),
'max_memory': get_max_memory()
@ -252,10 +253,13 @@ class Cropper:
bounds_gpkg_path = os.path.join(self.storage_dir, '{}.bounds.gpkg'.format(self.files_prefix))
if os.path.isfile(bounds_gpkg_path):
os.remove(bounds_gpkg_path)
# Convert bounds to GPKG
kwargs = {
'input': bounds_geojson_path,
'output': bounds_gpkg_path,
'input': double_quote(bounds_geojson_path),
'output': double_quote(bounds_gpkg_path),
'proj4': pc_proj4
}

Wyświetl plik

@ -155,7 +155,7 @@ def run_pipeline(json, verbose=False):
cmd = [
'pdal',
'pipeline',
'-i %s' % jsonfile
'-i %s' % double_quote(jsonfile)
]
if verbose or sys.platform == 'win32':
system.run(' '.join(cmd))

Wyświetl plik

@ -126,7 +126,7 @@ def dem_to_mesh_gridded(inGeotiff, outMesh, maxVertexCount, verbose=False, maxCo
system.run('"{reconstructmesh}" -i "{infile}" '
'-o "{outfile}" '
'--remove-spikes 0 --remove-spurious 20 --smooth 0 '
'--remove-spikes 0 --remove-spurious 0 --smooth 0 '
'--target-face-num {max_faces} '.format(**cleanupArgs))
# Delete intermediate results

Wyświetl plik

@ -70,7 +70,7 @@ def generate_kmz(orthophoto_file, output_file=None, outsize=None):
'--config GDAL_CACHEMAX %s%% ' % (orthophoto_file, output_file, bandparam, get_max_memory()))
def post_orthophoto_steps(args, bounds_file_path, orthophoto_file, orthophoto_tiles_dir):
if args.crop > 0:
if args.crop > 0 or args.boundary:
Cropper.crop(bounds_file_path, orthophoto_file, get_orthophoto_vars(args), keep_original=not args.optimize_disk_space, warp_options=['-dstalpha'])
if args.build_overviews and not args.cog:

Wyświetl plik

@ -510,10 +510,11 @@ def get_submodel_argv(args, submodels_path = None, submodel_name = None):
tweaking --crop if necessary (DEM merging makes assumption about the area of DEMs and their euclidean maps that require cropping. If cropping is skipped, this leads to errors.)
removing --gcp (the GCP path if specified is always "gcp_list.txt")
reading the contents of --cameras
reading the contents of --boundary
"""
assure_always = ['orthophoto_cutline', 'dem_euclidean_map', 'skip_3dmodel', 'skip_report']
remove_always = ['split', 'split_overlap', 'rerun_from', 'rerun', 'gcp', 'end_with', 'sm_cluster', 'rerun_all', 'pc_csv', 'pc_las', 'pc_ept', 'tiles', 'copy-to', 'cog']
read_json_always = ['cameras']
read_json_always = ['cameras', 'boundary']
argv = sys.argv

Wyświetl plik

@ -1,4 +1,4 @@
import os, sys, shutil, tempfile, json, math
import os, sys, shutil, tempfile, math, json
from opendm import system
from opendm import log
from opendm import context
@ -7,6 +7,8 @@ from opendm import entwine
from opendm import io
from opendm.concurrency import parallel_map
from opendm.utils import double_quote
from opendm.boundary import as_polygon
from opendm.dem.pdal import run_pipeline
def ply_info(input_ply):
if not os.path.exists(input_ply):
@ -38,7 +40,8 @@ def ply_info(input_ply):
return {
'has_normals': has_normals,
'vertex_count': vertex_count,
'has_views': has_views
'has_views': has_views,
'header_lines': i + 1
}
@ -68,7 +71,7 @@ def split(input_point_cloud, outdir, filename_template, capacity, dims=None):
return [os.path.join(outdir, f) for f in os.listdir(outdir)]
def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=16, sample_radius=0, verbose=False, max_concurrency=1):
def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=16, sample_radius=0, boundary=None, verbose=False, max_concurrency=1):
"""
Filters a point cloud
"""
@ -86,6 +89,10 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
log.ODM_INFO("Filtering {} (statistical, meanK {}, standard deviation {})".format(input_point_cloud, meank, standard_deviation))
filters.append('outlier')
filters.append('range')
if boundary is not None:
log.ODM_INFO("Boundary {}".format(boundary))
filters.append('crop')
info = ply_info(input_point_cloud)
dims = "x=float,y=float,z=float,"
@ -116,7 +123,8 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
filter(pcs['path'], io.related_file_path(pcs['path'], postfix="_filtered"),
standard_deviation=standard_deviation,
meank=meank,
sample_radius=sample_radius,
sample_radius=sample_radius,
boundary=boundary,
verbose=verbose,
max_concurrency=1)
# Filter
@ -132,34 +140,40 @@ def filter(input_point_cloud, output_point_cloud, standard_deviation=2.5, meank=
shutil.rmtree(partsdir)
else:
# Process point cloud (or a point cloud submodel) in a single step
filterArgs = {
'inputFile': input_point_cloud,
'outputFile': output_point_cloud,
'stages': " ".join(filters),
'dims': dims
}
pipeline = []
cmd = ("pdal translate -i \"{inputFile}\" "
"-o \"{outputFile}\" "
"{stages} "
"--writers.ply.sized_types=false "
"--writers.ply.storage_mode=\"little endian\" "
"--writers.ply.dims=\"{dims}\" "
"").format(**filterArgs)
# Input
pipeline.append(input_point_cloud)
if 'sample' in filters:
cmd += "--filters.sample.radius={} ".format(sample_radius)
if 'outlier' in filters:
cmd += ("--filters.outlier.method=\"statistical\" "
"--filters.outlier.mean_k={} "
"--filters.outlier.multiplier={} ").format(meank, standard_deviation)
if 'range' in filters:
# Remove outliers
cmd += "--filters.range.limits=\"Classification![7:7]\" "
# Filters
for f in filters:
params = {}
if f == 'sample':
params = {'radius': sample_radius}
elif f == 'outlier':
params = {'method': 'statistical', 'mean_k': meank, 'multiplier': standard_deviation}
elif f == 'range':
params = {'limits': 'Classification![7:7]'}
elif f == 'crop':
params = {'polygon': as_polygon(boundary)}
else:
raise RuntimeError("Invalid filter in PDAL pipeline (this should not have happened, please report it: https://github.com/OpenDroneMap/ODM/issues")
pipeline.append(dict({
'type': "filters.%s" % f,
}, **params))
system.run(cmd)
# Output
pipeline.append({
'type': 'writers.ply',
'sized_types': False,
'storage_mode': 'little endian',
'dims': dims,
'filename': output_point_cloud
})
run_pipeline(pipeline, verbose=verbose)
if not os.path.exists(output_point_cloud):
log.ODM_WARNING("{} not found, filtering has failed.".format(output_point_cloud))

Wyświetl plik

@ -69,7 +69,7 @@ class ODM_Reconstruction(object):
return self.georef is not None
def has_gcp(self):
return self.is_georeferenced() and self.gcp is not None
return self.is_georeferenced() and self.gcp is not None and self.gcp.exists()
def georeference_with_gcp(self, gcp_file, output_coords_file, output_gcp_file, output_model_txt_geo, rerun=False):
if not io.file_exists(output_coords_file) or not io.file_exists(output_gcp_file) or rerun:
@ -158,6 +158,12 @@ class ODM_Reconstruction(object):
def get_proj_srs(self):
if self.is_georeferenced():
return self.georef.proj4()
def get_proj_offset(self):
if self.is_georeferenced():
return (self.georef.utm_east_offset, self.georef.utm_north_offset)
else:
return (None, None)
def get_photo(self, filename):
for p in self.photos:

Wyświetl plik

@ -9,7 +9,7 @@ from opendm import system
from opendm.geo import GeoFile
from shutil import copyfile
from opendm import progress
from opendm import boundary
def save_images_database(photos, database_file):
with open(database_file, 'w') as f:
@ -154,3 +154,11 @@ class ODMLoadDatasetStage(types.ODM_Stage):
reconstruction.save_proj_srs(os.path.join(tree.odm_georeferencing, tree.odm_georeferencing_proj))
outputs['reconstruction'] = reconstruction
# Try to load boundaries
if args.boundary:
if reconstruction.is_georeferenced():
outputs['boundary'] = boundary.load_boundary(args.boundary, reconstruction.get_proj_srs())
else:
args.boundary = None
log.ODM_WARNING("Reconstruction is not georeferenced, but boundary file provided (will ignore boundary file)")

Wyświetl plik

@ -111,14 +111,14 @@ class ODMDEMStage(types.ODM_Stage):
dem_geotiff_path = os.path.join(odm_dem_root, "{}.tif".format(product))
bounds_file_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.gpkg')
if args.crop > 0:
if args.crop > 0 or args.boundary:
# Crop DEM
Cropper.crop(bounds_file_path, dem_geotiff_path, utils.get_dem_vars(args), keep_original=not args.optimize_disk_space)
if args.dem_euclidean_map:
unfilled_dem_path = io.related_file_path(dem_geotiff_path, postfix=".unfilled")
if args.crop > 0:
if args.crop > 0 or args.boundary:
# Crop unfilled DEM
Cropper.crop(bounds_file_path, unfilled_dem_path, utils.get_dem_vars(args), keep_original=not args.optimize_disk_space)

Wyświetl plik

@ -6,6 +6,8 @@ from opendm import system
from opendm import context
from opendm import point_cloud
from opendm import types
from opendm import gsd
from opendm.boundary import boundary_offset, compute_boundary_from_shots
class ODMFilterPoints(types.ODM_Stage):
def process(self, args, outputs):
@ -21,12 +23,33 @@ class ODMFilterPoints(types.ODM_Stage):
else:
inputPointCloud = tree.openmvs_model
# Check if we need to compute boundary
if args.auto_boundary:
if reconstruction.is_georeferenced():
if not 'boundary' in outputs:
avg_gsd = gsd.opensfm_reconstruction_average_gsd(tree.opensfm_reconstruction)
outputs['boundary'] = compute_boundary_from_shots(tree.opensfm_reconstruction, avg_gsd * 20, reconstruction.get_proj_offset()) # 20 is arbitrary
if outputs['boundary'] is None:
log.ODM_WARNING("Cannot compute boundary from camera shots")
else:
log.ODM_WARNING("--auto-boundary set but so is --boundary, will use --boundary")
else:
log.ODM_WARNING("Not a georeferenced reconstruction, will ignore --auto-boundary")
point_cloud.filter(inputPointCloud, tree.filtered_point_cloud,
standard_deviation=args.pc_filter,
sample_radius=args.pc_sample,
boundary=boundary_offset(outputs.get('boundary'), reconstruction.get_proj_offset()),
verbose=args.verbose,
max_concurrency=args.max_concurrency)
# Quick check
info = point_cloud.ply_info(tree.filtered_point_cloud)
if info["vertex_count"] == 0:
extra_msg = ''
if 'boundary' in outputs:
extra_msg = '. Also, since you used a boundary setting, make sure that the boundary polygon you specified covers the reconstruction area correctly.'
raise system.ExitException("Uh oh! We ended up with an empty point cloud. This means that the reconstruction did not succeed. Have you followed best practices for data acquisition? See https://docs.opendronemap.org/flying/%s" % extra_msg)
else:
log.ODM_WARNING('Found a valid point cloud file in: %s' %
tree.filtered_point_cloud)

Wyświetl plik

@ -17,7 +17,7 @@ from opendm.cropper import Cropper
from opendm import point_cloud
from opendm.multispectral import get_primary_band_name
from opendm.osfm import OSFMContext
from opendm.boundary import as_polygon, export_to_bounds_files
class ODMGeoreferencingStage(types.ODM_Stage):
def process(self, args, outputs):
@ -152,6 +152,14 @@ class ODMGeoreferencingStage(types.ODM_Stage):
except:
log.ODM_WARNING("Cannot calculate crop bounds! We will skip cropping")
args.crop = 0
if 'boundary' in outputs and args.crop == 0:
log.ODM_INFO("Using boundary JSON as cropping area")
bounds_base, _ = os.path.splitext(tree.odm_georeferencing_model_laz)
bounds_json = bounds_base + ".bounds.geojson"
bounds_gpkg = bounds_base + ".bounds.gpkg"
export_to_bounds_files(outputs['boundary'], reconstruction.get_proj_srs(), bounds_json, bounds_gpkg)
else:
log.ODM_INFO("Converting point cloud (non-georeferenced)")
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))

Wyświetl plik

@ -152,7 +152,7 @@ class ODMReport(types.ODM_Stage):
overlap_color_map = os.path.join(report_assets, "overlap_color_map.txt")
bounds_file_path = os.path.join(tree.odm_georeferencing, 'odm_georeferenced_model.bounds.gpkg')
if args.crop > 0 and os.path.isfile(bounds_file_path):
if (args.crop > 0 or args.boundary) and os.path.isfile(bounds_file_path):
Cropper.crop(bounds_file_path, diagram_tiff, get_orthophoto_vars(args), keep_original=False)
system.run("gdaldem color-relief \"{}\" \"{}\" \"{}\" -of PNG -alpha".format(diagram_tiff, overlap_color_map, diagram_png))

Wyświetl plik

@ -105,71 +105,6 @@ class ODMSplitStage(types.ODM_Stage):
self.update_progress(50)
# TODO: this is currently not working and needs a champion to fix it
# https://community.opendronemap.org/t/filenotfound-error-cameras-json/6047/2
# resplit_done_file = octx.path('resplit_done.txt')
# if not io.file_exists(resplit_done_file) and bool(args.split_multitracks):
# submodels = mds.get_submodel_paths()
# i = 0
# for s in submodels:
# template = octx.path("../aligned_submodels/submodel_%04d")
# with open(s+"/reconstruction.json", "r") as f:
# j = json.load(f)
# for k in range(0, len(j)):
# v = j[k]
# path = template % i
# #Create the submodel path up to opensfm
# os.makedirs(path+"/opensfm")
# os.makedirs(path+"/images")
# #symlinks for common data
# images = os.listdir(octx.path("../images"))
# for image in images:
# os.symlink("../../../images/"+image, path+"/images/"+image)
# os.symlink("../../../opensfm/exif", path+"/opensfm/exif")
# os.symlink("../../../opensfm/features", path+"/opensfm/features")
# os.symlink("../../../opensfm/matches", path+"/opensfm/matches")
# os.symlink("../../../opensfm/reference_lla.json", path+"/opensfm/reference_lla.json")
# os.symlink("../../../opensfm/camera_models.json", path+"/opensfm/camera_models.json")
# shutil.copy(s+"/../cameras.json", path+"/cameras.json")
# shutil.copy(s+"/../images.json", path+"/images.json")
# with open(octx.path("config.yaml")) as f:
# doc = yaml.safe_load(f)
# dmcv = "depthmap_min_consistent_views"
# if dmcv in doc:
# if len(v["shots"]) < doc[dmcv]:
# doc[dmcv] = len(v["shots"])
# print("WARNING: Reduced "+dmcv+" to accommodate short track")
# with open(path+"/opensfm/config.yaml", "w") as f:
# yaml.dump(doc, f)
# #We need the original tracks file for the visualsfm export, since
# #there may still be point matches between the tracks
# shutil.copy(s+"/tracks.csv", path+"/opensfm/tracks.csv")
# #Create our new reconstruction file with only the relevant track
# with open(path+"/opensfm/reconstruction.json", "w") as o:
# json.dump([v], o)
# #Create image lists
# with open(path+"/opensfm/image_list.txt", "w") as o:
# o.writelines(list(map(lambda x: "../images/"+x+'\n', v["shots"].keys())))
# with open(path+"/img_list.txt", "w") as o:
# o.writelines(list(map(lambda x: x+'\n', v["shots"].keys())))
# i+=1
# os.rename(octx.path("../submodels"), octx.path("../unaligned_submodels"))
# os.rename(octx.path("../aligned_submodels"), octx.path("../submodels"))
# octx.touch(resplit_done_file)
mds = metadataset.MetaDataSet(tree.opensfm)
submodel_paths = [os.path.abspath(p) for p in mds.get_submodel_paths()]
@ -332,7 +267,7 @@ class ODMMergeStage(types.ODM_Stage):
if io.file_exists(dem_file):
# Crop
if args.crop > 0:
if args.crop > 0 or args.boundary:
Cropper.crop(merged_bounds_file, dem_file, dem_vars, keep_original=not args.optimize_disk_space)
log.ODM_INFO("Created %s" % dem_file)