kopia lustrzana https://github.com/OpenDroneMap/ODM
Merge pull request #1068 from pierotofy/pseudogeo
Pseudo-georeferencing for datasets missing GPS/GCP information, AKA putting bananas on the mappull/1073/head
Piero Toffanin
GitHub
commit
ed551d6218
|
|
@ -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,9 +114,8 @@ def opensfm_reconstruction_average_gsd(reconstruction_json):
|
|||
gsds = []
|
||||
for shotImage in reconstruction['shots']:
|
||||
shot = reconstruction['shots'][shotImage]
|
||||
if shot['gps_dop'] < 999999:
|
||||
if use_all_shots or shot['gps_dop'] < 999999:
|
||||
camera = reconstruction['cameras'][shot['camera']]
|
||||
|
||||
shot_height = shot['translation'][2]
|
||||
focal_ratio = camera.get('focal', camera.get('focal_x'))
|
||||
if not focal_ratio:
|
||||
|
|
|
|||
|
|
@ -0,0 +1,23 @@
|
|||
import osr
|
||||
import gdal
|
||||
from gdalconst import GA_Update
|
||||
from opendm import io
|
||||
from opendm import log
|
||||
|
||||
def add_pseudo_georeferencing(geotiff, scale=1.0):
|
||||
if not io.file_exists(geotiff):
|
||||
log.ODM_WARNING("Cannot add pseudo georeferencing, %s does not exist" % geotiff)
|
||||
return
|
||||
|
||||
try:
|
||||
log.ODM_INFO("Adding pseudo georeferencing (raster should show up at the equator) to %s" % geotiff)
|
||||
|
||||
dst_ds = gdal.Open(geotiff, GA_Update)
|
||||
srs = osr.SpatialReference()
|
||||
srs.ImportFromProj4('+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
|
||||
dst_ds.SetProjection( srs.ExportToWkt() )
|
||||
dst_ds.SetGeoTransform( [ 0.0, scale, 0.0, 0.0, 0.0, -scale ] )
|
||||
dst_ds = None
|
||||
|
||||
except Exception as e:
|
||||
log.ODM_WARNING("Cannot add psuedo georeferencing to %s (%s), skipping..." % (geotiff, str(e)))
|
||||
|
|
@ -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)
|
||||
|
|
@ -238,7 +241,7 @@ class ODM_Reconstruction(object):
|
|||
|
||||
self.georef = ODM_GeoRef.FromCoordsFile(output_coords_file)
|
||||
except:
|
||||
log.ODM_WARNING('Could not generate coordinates file. An orthophoto will not be generated.')
|
||||
log.ODM_WARNING('Could not generate coordinates file. The orthophoto will not be georeferenced.')
|
||||
|
||||
self.gcp = GCPFile(None)
|
||||
return self.georef
|
||||
|
|
|
|||
|
|
@ -9,28 +9,55 @@ 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 reconstruction.is_georeferenced():
|
||||
# Special case to clear previous run point cloud
|
||||
# (NodeODM will generate a fake georeferenced laz during postprocessing
|
||||
# with non-georeferenced datasets). odm_georeferencing_model_laz should
|
||||
# not be here! Perhaps we should improve this.
|
||||
if io.file_exists(tree.odm_georeferencing_model_laz) and self.rerun():
|
||||
os.remove(tree.odm_georeferencing_model_laz)
|
||||
|
||||
log.ODM_WARNING("Not georeferenced, using ungeoreferenced point cloud...")
|
||||
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 +76,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 +88,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 +125,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)
|
||||
|
|
|
|||
|
|
@ -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))
|
||||
|
||||
|
|
|
|||
|
|
@ -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,18 @@ 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
|
||||
ignore_resolution = False
|
||||
if not reconstruction.is_georeferenced():
|
||||
# Match DEMs
|
||||
gsd_error_estimate = -3
|
||||
ignore_resolution = True
|
||||
|
||||
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=ignore_resolution,
|
||||
has_gcp=reconstruction.has_gcp()) / 100.0)
|
||||
|
||||
# odm_orthophoto definitions
|
||||
kwargs = {
|
||||
|
|
@ -28,7 +41,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 +56,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 +159,12 @@ 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("Could not generate an orthophoto (it did not render)")
|
||||
else:
|
||||
log.ODM_WARNING('Found a valid orthophoto in: %s' % tree.odm_orthophoto_tif)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue