mirror
/
OpenDroneMap-ODM
kopia lustrzana https://github.com/OpenDroneMap/ODM
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

Pseudo-georeferencing for datasets missing GPS/GCP information

pull/1068/head
Piero Toffanin 2020-02-04 17:56:20 +00:00
rodzic b4116da492
commit 4d8012026a
6 zmienionych plików z 78 dodań i 35 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

39
opendm/gsd.py
Wyświetl plik

@ -21,7 +21,7 @@ def rounded_gsd(reconstruction_json, default_value=None, ndigits=0, ignore_gsd=F
return default_value
def image_max_size(photos, target_resolution, reconstruction_json, gsd_error_estimate = 0.5, ignore_gsd=False):
def image_max_size(photos, target_resolution, reconstruction_json, gsd_error_estimate = 0.5, ignore_gsd=False, has_gcp=False):
"""
:param photos images database
:param target_resolution resolution the user wants have in cm / pixel
@ -36,7 +36,7 @@ def image_max_size(photos, target_resolution, reconstruction_json, gsd_error_est
if ignore_gsd:
isf = 1.0
else:
isf = image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimate)
isf = image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimate, has_gcp=has_gcp)
for p in photos:
max_width = max(p.width, max_width)
@ -44,7 +44,7 @@ def image_max_size(photos, target_resolution, reconstruction_json, gsd_error_est
return int(math.ceil(max(max_width, max_height) * isf))
def image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimate = 0.5):
def image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimate = 0.5, has_gcp=False):
"""
:param target_resolution resolution the user wants have in cm / pixel
:param reconstruction_json path to OpenSfM's reconstruction.json
@ -52,7 +52,7 @@ def image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimat
:return A down-scale (<= 1) value to apply to images to achieve the target resolution by comparing the current GSD of the reconstruction.
If a GSD cannot be computed, it just returns 1. Returned scale values are never higher than 1.
"""
gsd = opensfm_reconstruction_average_gsd(reconstruction_json)
gsd = opensfm_reconstruction_average_gsd(reconstruction_json, use_all_shots=has_gcp)
if gsd is not None and target_resolution > 0:
gsd = gsd * (1 + gsd_error_estimate)
@ -61,7 +61,7 @@ def image_scale_factor(target_resolution, reconstruction_json, gsd_error_estimat
return 1
def cap_resolution(resolution, reconstruction_json, gsd_error_estimate = 0.1, ignore_gsd=False):
def cap_resolution(resolution, reconstruction_json, gsd_error_estimate = 0.1, ignore_gsd=False, ignore_resolution=False, has_gcp=False):
"""
:param resolution resolution in cm / pixel
:param reconstruction_json path to OpenSfM's reconstruction.json
@ -73,11 +73,11 @@ def cap_resolution(resolution, reconstruction_json, gsd_error_estimate = 0.1, ig
if ignore_gsd:
return resolution
gsd = opensfm_reconstruction_average_gsd(reconstruction_json)
gsd = opensfm_reconstruction_average_gsd(reconstruction_json, use_all_shots=has_gcp or ignore_resolution)
if gsd is not None:
gsd = gsd * (1 - gsd_error_estimate)
if gsd > resolution:
if gsd > resolution or ignore_resolution:
log.ODM_WARNING('Maximum resolution set to GSD - {}% ({} cm / pixel, requested resolution was {} cm / pixel)'.format(gsd_error_estimate * 100, round(gsd, 2), round(resolution, 2)))
return gsd
else:
@ -88,7 +88,7 @@ def cap_resolution(resolution, reconstruction_json, gsd_error_estimate = 0.1, ig
@lru_cache(maxsize=None)
def opensfm_reconstruction_average_gsd(reconstruction_json):
def opensfm_reconstruction_average_gsd(reconstruction_json, use_all_shots=False):
"""
Computes the average Ground Sampling Distance of an OpenSfM reconstruction.
:param reconstruction_json path to OpenSfM's reconstruction.json
@ -114,18 +114,19 @@ def opensfm_reconstruction_average_gsd(reconstruction_json):
gsds = []
for shotImage in reconstruction['shots']:
shot = reconstruction['shots'][shotImage]
if shot['gps_dop'] < 999999:
camera = reconstruction['cameras'][shot['camera']]
if not use_all_shots and shot['gps_dop'] >= 999999:
continue
shot_height = shot['translation'][2]
focal_ratio = camera.get('focal', camera.get('focal_x'))
if not focal_ratio:
log.ODM_WARNING("Cannot parse focal values from %s. This is likely an unsupported camera model." % reconstruction_json)
return None
gsds.append(calculate_gsd_from_focal_ratio(focal_ratio,
shot_height - ground_height,
camera['width']))
camera = reconstruction['cameras'][shot['camera']]
shot_height = shot['translation'][2]
focal_ratio = camera.get('focal', camera.get('focal_x'))
if not focal_ratio:
log.ODM_WARNING("Cannot parse focal values from %s. This is likely an unsupported camera model." % reconstruction_json)
return None
gsds.append(calculate_gsd_from_focal_ratio(focal_ratio,
shot_height - ground_height,
camera['width']))
if len(gsds) > 0:
mean = np.mean(gsds)

3
opendm/types.py
Wyświetl plik

@ -190,6 +190,9 @@ class ODM_Reconstruction(object):
def is_georeferenced(self):
return self.georef is not None
def has_gcp(self):
return self.is_georeferenced() and self.gcp is not None
def georeference_with_gcp(self, gcp_file, output_coords_file, output_gcp_file, rerun=False):
if not io.file_exists(output_coords_file) or not io.file_exists(output_gcp_file) or rerun:
gcp = GCPFile(gcp_file)

41
stages/odm_dem.py
Wyświetl plik

@ -9,28 +9,48 @@ from opendm import types
from opendm import gsd
from opendm.dem import commands, utils
from opendm.cropper import Cropper
from opendm import pseudogeo
class ODMDEMStage(types.ODM_Stage):
def process(self, args, outputs):
tree = outputs['tree']
las_model_found = io.file_exists(tree.odm_georeferencing_model_laz)
reconstruction = outputs['reconstruction']
dem_input = tree.odm_georeferencing_model_laz
pc_model_found = io.file_exists(dem_input)
ignore_resolution = False
pseudo_georeference = False
if not pc_model_found:
log.ODM_WARNING("Georeferenced point cloud not found, trying ungeoreferenced...")
dem_input = tree.path("odm_filterpoints", "point_cloud.ply")
pc_model_found = io.file_exists(dem_input)
ignore_resolution = True
pseudo_georeference = True
resolution = gsd.cap_resolution(args.dem_resolution, tree.opensfm_reconstruction,
gsd_error_estimate=-3,
ignore_gsd=args.ignore_gsd,
ignore_resolution=ignore_resolution,
has_gcp=reconstruction.has_gcp())
log.ODM_INFO('Classify: ' + str(args.pc_classify))
log.ODM_INFO('Create DSM: ' + str(args.dsm))
log.ODM_INFO('Create DTM: ' + str(args.dtm))
log.ODM_INFO('DEM input file {0} found: {1}'.format(tree.odm_georeferencing_model_laz, str(las_model_found)))
log.ODM_INFO('DEM input file {0} found: {1}'.format(dem_input, str(pc_model_found)))
# define paths and create working directories
odm_dem_root = tree.path('odm_dem')
if not io.dir_exists(odm_dem_root):
system.mkdir_p(odm_dem_root)
if args.pc_classify and las_model_found:
if args.pc_classify and pc_model_found:
pc_classify_marker = os.path.join(odm_dem_root, 'pc_classify_done.txt')
if not io.file_exists(pc_classify_marker) or self.rerun():
log.ODM_INFO("Classifying {} using Simple Morphological Filter".format(tree.odm_georeferencing_model_laz))
commands.classify(tree.odm_georeferencing_model_laz,
log.ODM_INFO("Classifying {} using Simple Morphological Filter".format(dem_input))
commands.classify(dem_input,
args.smrf_scalar,
args.smrf_slope,
args.smrf_threshold,
@ -49,7 +69,7 @@ class ODMDEMStage(types.ODM_Stage):
self.update_progress(progress)
# Do we need to process anything here?
if (args.dsm or args.dtm) and las_model_found:
if (args.dsm or args.dtm) and pc_model_found:
dsm_output_filename = os.path.join(odm_dem_root, 'dsm.tif')
dtm_output_filename = os.path.join(odm_dem_root, 'dtm.tif')
@ -61,15 +81,14 @@ class ODMDEMStage(types.ODM_Stage):
if args.dsm or (args.dtm and args.dem_euclidean_map): products.append('dsm')
if args.dtm: products.append('dtm')
resolution = gsd.cap_resolution(args.dem_resolution, tree.opensfm_reconstruction, gsd_error_estimate=-3, ignore_gsd=args.ignore_gsd)
radius_steps = [(resolution / 100.0) / 2.0]
for _ in range(args.dem_gapfill_steps - 1):
radius_steps.append(radius_steps[-1] * 2) # 2 is arbitrary, maybe there's a better value?
for product in products:
commands.create_dem(
tree.odm_georeferencing_model_laz,
dem_input,
product,
output_type='idw' if product == 'dtm' else 'max',
radiuses=map(str, radius_steps),
@ -99,6 +118,10 @@ class ODMDEMStage(types.ODM_Stage):
commands.compute_euclidean_map(unfilled_dem_path,
io.related_file_path(dem_geotiff_path, postfix=".euclideand"),
overwrite=True)
if pseudo_georeference:
# 0.1 is arbitrary
pseudogeo.add_pseudo_georeferencing(dem_geotiff_path, 0.1)
progress += 30
self.update_progress(progress)

5
stages/odm_meshing.py
Wyświetl plik

@ -42,7 +42,10 @@ class ODMeshingStage(types.ODM_Stage):
if not io.file_exists(tree.odm_25dmesh) or self.rerun():
log.ODM_INFO('Writing ODM 2.5D Mesh file in: %s' % tree.odm_25dmesh)
ortho_resolution = gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd) / 100.0
ortho_resolution = gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction,
ignore_gsd=args.ignore_gsd,
ignore_resolution=not reconstruction.is_georeferenced(),
has_gcp=reconstruction.has_gcp()) / 100.0
dsm_multiplier = max(1.0, gsd.rounded_gsd(tree.opensfm_reconstruction, default_value=4, ndigits=3, ignore_gsd=args.ignore_gsd))

23
stages/odm_orthophoto.py
Wyświetl plik

@ -10,6 +10,7 @@ from opendm import orthophoto
from opendm.concurrency import get_max_memory
from opendm.cutline import compute_cutline
from pipes import quote
from opendm import pseudogeo
class ODMOrthoPhotoStage(types.ODM_Stage):
def process(self, args, outputs):
@ -21,6 +22,16 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
system.mkdir_p(tree.odm_orthophoto)
if not io.file_exists(tree.odm_orthophoto_tif) or self.rerun():
gsd_error_estimate = 0.1
if not reconstruction.is_georeferenced():
# Match DEMs
gsd_error_estimate = -3
resolution = 1.0 / (gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction,
gsd_error_estimate=gsd_error_estimate,
ignore_gsd=args.ignore_gsd,
ignore_resolution=not reconstruction.is_georeferenced(),
has_gcp=reconstruction.has_gcp()) / 100.0)
# odm_orthophoto definitions
kwargs = {
@ -28,7 +39,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
'log': tree.odm_orthophoto_log,
'ortho': tree.odm_orthophoto_render,
'corners': tree.odm_orthophoto_corners,
'res': 1.0 / (gsd.cap_resolution(args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd) / 100.0),
'res': resolution,
'bands': '',
'verbose': verbose
}
@ -43,7 +54,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
if io.file_exists(odm_georeferencing_model_txt_geo_file):
reconstruction.georef.extract_offsets(odm_georeferencing_model_txt_geo_file)
else:
log.ODM_WARNING('Cannot read UTM offset from {}. An orthophoto will not be generated.'.format(odm_georeferencing_model_txt_geo_file))
log.ODM_WARNING('Cannot read UTM offset from {}.'.format(odm_georeferencing_model_txt_geo_file))
models = []
@ -146,8 +157,10 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
geotiffcreated = True
if not geotiffcreated:
log.ODM_WARNING('No geo-referenced orthophoto created due '
'to missing geo-referencing or corner coordinates.')
if io.file_exists(tree.odm_orthophoto_render):
# 0.1 is arbitrary
pseudogeo.add_pseudo_georeferencing(tree.odm_orthophoto_render, 0.1)
log.ODM_INFO("Renaming %s --> %s" % (tree.odm_orthophoto_render, tree.odm_orthophoto_tif))
os.rename(tree.odm_orthophoto_render, tree.odm_orthophoto_tif)
else:
log.ODM_WARNING('Found a valid orthophoto in: %s' % tree.odm_orthophoto_tif)

2
stages/run_opensfm.py
Wyświetl plik

@ -48,7 +48,7 @@ class ODMOpenSfMStage(types.ODM_Stage):
# Make sure it's capped by the depthmap-resolution arg,
# since the undistorted images are used for MVS
outputs['undist_image_max_size'] = max(
gsd.image_max_size(photos, args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd),
gsd.image_max_size(photos, args.orthophoto_resolution, tree.opensfm_reconstruction, ignore_gsd=args.ignore_gsd, has_gcp=reconstruction.has_gcp()),
args.depthmap_resolution
)