kopia lustrzana https://github.com/OpenDroneMap/ODM
PoC codem integration
Piero Toffanin
rodzic
d3832e1cf5
commit
9ddb4554c4
|
|
@ -0,0 +1,34 @@
|
|||
import os
|
||||
import shutil
|
||||
import codem
|
||||
import dataclasses
|
||||
from opendm import log
|
||||
|
||||
def compute_alignment_matrix(input_laz, align_file, tmp_dir):
|
||||
if os.path.exists(tmp_dir):
|
||||
shutil.rmtree(tmp_dir)
|
||||
os.mkdir(tmp_dir)
|
||||
|
||||
conf = dataclasses.asdict(codem.CodemRunConfig(align_file, input_laz, MIN_RESOLUTION=0.2, OUTPUT_DIR=tmp_dir)) # TODO: how to compute this
|
||||
fnd_obj, aoi_obj = codem.preprocess(conf)
|
||||
fnd_obj.prep()
|
||||
aoi_obj.prep()
|
||||
log.ODM_INFO("Aligning reconstruction to %s" % align_file)
|
||||
log.ODM_INFO("Coarse registration...")
|
||||
dsm_reg = codem.coarse_registration(fnd_obj, aoi_obj, conf)
|
||||
log.ODM_INFO("Fine registration...")
|
||||
icp_reg = codem.fine_registration(fnd_obj, aoi_obj, dsm_reg, conf)
|
||||
|
||||
app_reg = codem.registration.ApplyRegistration(
|
||||
fnd_obj,
|
||||
aoi_obj,
|
||||
icp_reg.registration_parameters,
|
||||
icp_reg.residual_vectors,
|
||||
icp_reg.residual_origins,
|
||||
conf,
|
||||
None,
|
||||
)
|
||||
|
||||
return app_reg.get_registration_transformation()
|
||||
# if os.path.exists(tmp_dir):
|
||||
# shutil.rmtree(tmp_dir)
|
||||
|
|
@ -470,6 +470,14 @@ def config(argv=None, parser=None):
|
|||
'EPSG:<code> or <+proj definition>\n'
|
||||
'image_name geo_x geo_y geo_z [omega (degrees)] [phi (degrees)] [kappa (degrees)] [horz accuracy (meters)] [vert accuracy (meters)]\n'
|
||||
'Default: %(default)s'))
|
||||
|
||||
parser.add_argument('--align-to',
|
||||
metavar='<path string>',
|
||||
action=StoreValue,
|
||||
default=None,
|
||||
help=('Path to a GeoTIFF DEM or a LAS/LAZ point cloud '
|
||||
'that the reconstruction outputs should be automatically aligned to. '
|
||||
'Default: %(default)s'))
|
||||
|
||||
parser.add_argument('--use-exif',
|
||||
action=StoreTrue,
|
||||
|
|
|
|||
|
|
@ -241,7 +241,7 @@ class ODM_GeoRef(object):
|
|||
return (self.utm_east_offset, self.utm_north_offset)
|
||||
|
||||
class ODM_Tree(object):
|
||||
def __init__(self, root_path, gcp_file = None, geo_file = None):
|
||||
def __init__(self, root_path, gcp_file = None, geo_file = None, align_file = None):
|
||||
# root path to the project
|
||||
self.root_path = io.absolute_path_file(root_path)
|
||||
self.input_images = os.path.join(self.root_path, 'images')
|
||||
|
|
@ -296,6 +296,7 @@ class ODM_Tree(object):
|
|||
self.odm_georeferencing_gcp = gcp_file or io.find('gcp_list.txt', self.root_path)
|
||||
self.odm_georeferencing_gcp_utm = os.path.join(self.odm_georeferencing, 'gcp_list_utm.txt')
|
||||
self.odm_geo_file = geo_file or io.find('geo.txt', self.root_path)
|
||||
self.odm_align_file = align_file or io.find('align.laz', self.root_path)
|
||||
|
||||
self.odm_georeferencing_proj = 'proj.txt'
|
||||
self.odm_georeferencing_model_txt_geo = os.path.join(
|
||||
|
|
|
|||
|
|
@ -46,7 +46,7 @@ def load_images_database(database_file):
|
|||
class ODMLoadDatasetStage(types.ODM_Stage):
|
||||
def process(self, args, outputs):
|
||||
outputs['start_time'] = system.now_raw()
|
||||
tree = types.ODM_Tree(args.project_path, args.gcp, args.geo)
|
||||
tree = types.ODM_Tree(args.project_path, args.gcp, args.geo, args.align_to)
|
||||
outputs['tree'] = tree
|
||||
|
||||
if io.file_exists(tree.benchmarking):
|
||||
|
|
|
|||
|
|
@ -18,6 +18,7 @@ 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
|
||||
from opendm.align import compute_alignment_matrix
|
||||
|
||||
class ODMGeoreferencingStage(types.ODM_Stage):
|
||||
def process(self, args, outputs):
|
||||
|
|
@ -30,143 +31,158 @@ class ODMGeoreferencingStage(types.ODM_Stage):
|
|||
gcp_gml_export_file = tree.path("odm_georeferencing", "ground_control_points.gml")
|
||||
gcp_geojson_export_file = tree.path("odm_georeferencing", "ground_control_points.geojson")
|
||||
|
||||
if reconstruction.has_gcp() and (not io.file_exists(gcp_export_file) or self.rerun()):
|
||||
octx = OSFMContext(tree.opensfm)
|
||||
gcps = octx.ground_control_points(reconstruction.georef.proj4())
|
||||
# if reconstruction.has_gcp() and (not io.file_exists(gcp_export_file) or self.rerun()):
|
||||
# octx = OSFMContext(tree.opensfm)
|
||||
# gcps = octx.ground_control_points(reconstruction.georef.proj4())
|
||||
|
||||
if len(gcps):
|
||||
gcp_schema = {
|
||||
'geometry': 'Point',
|
||||
'properties': OrderedDict([
|
||||
('id', 'str'),
|
||||
('observations_count', 'int'),
|
||||
('observations_list', 'str'),
|
||||
('error_x', 'float'),
|
||||
('error_y', 'float'),
|
||||
('error_z', 'float'),
|
||||
])
|
||||
}
|
||||
# if len(gcps):
|
||||
# gcp_schema = {
|
||||
# 'geometry': 'Point',
|
||||
# 'properties': OrderedDict([
|
||||
# ('id', 'str'),
|
||||
# ('observations_count', 'int'),
|
||||
# ('observations_list', 'str'),
|
||||
# ('error_x', 'float'),
|
||||
# ('error_y', 'float'),
|
||||
# ('error_z', 'float'),
|
||||
# ])
|
||||
# }
|
||||
|
||||
# Write GeoPackage
|
||||
with fiona.open(gcp_export_file, 'w', driver="GPKG",
|
||||
crs=fiona.crs.from_string(reconstruction.georef.proj4()),
|
||||
schema=gcp_schema) as f:
|
||||
for gcp in gcps:
|
||||
f.write({
|
||||
'geometry': {
|
||||
'type': 'Point',
|
||||
'coordinates': gcp['coordinates'],
|
||||
},
|
||||
'properties': OrderedDict([
|
||||
('id', gcp['id']),
|
||||
('observations_count', len(gcp['observations'])),
|
||||
('observations_list', ",".join([obs['shot_id'] for obs in gcp['observations']])),
|
||||
('error_x', gcp['error'][0]),
|
||||
('error_y', gcp['error'][1]),
|
||||
('error_z', gcp['error'][2]),
|
||||
])
|
||||
})
|
||||
# # Write GeoPackage
|
||||
# with fiona.open(gcp_export_file, 'w', driver="GPKG",
|
||||
# crs=fiona.crs.from_string(reconstruction.georef.proj4()),
|
||||
# schema=gcp_schema) as f:
|
||||
# for gcp in gcps:
|
||||
# f.write({
|
||||
# 'geometry': {
|
||||
# 'type': 'Point',
|
||||
# 'coordinates': gcp['coordinates'],
|
||||
# },
|
||||
# 'properties': OrderedDict([
|
||||
# ('id', gcp['id']),
|
||||
# ('observations_count', len(gcp['observations'])),
|
||||
# ('observations_list', ",".join([obs['shot_id'] for obs in gcp['observations']])),
|
||||
# ('error_x', gcp['error'][0]),
|
||||
# ('error_y', gcp['error'][1]),
|
||||
# ('error_z', gcp['error'][2]),
|
||||
# ])
|
||||
# })
|
||||
|
||||
# Write GML
|
||||
try:
|
||||
system.run('ogr2ogr -of GML "{}" "{}"'.format(gcp_gml_export_file, gcp_export_file))
|
||||
except Exception as e:
|
||||
log.ODM_WARNING("Cannot generate ground control points GML file: %s" % str(e))
|
||||
# # Write GML
|
||||
# try:
|
||||
# system.run('ogr2ogr -of GML "{}" "{}"'.format(gcp_gml_export_file, gcp_export_file))
|
||||
# except Exception as e:
|
||||
# log.ODM_WARNING("Cannot generate ground control points GML file: %s" % str(e))
|
||||
|
||||
# Write GeoJSON
|
||||
geojson = {
|
||||
'type': 'FeatureCollection',
|
||||
'features': []
|
||||
}
|
||||
# # Write GeoJSON
|
||||
# geojson = {
|
||||
# 'type': 'FeatureCollection',
|
||||
# 'features': []
|
||||
# }
|
||||
|
||||
from_srs = CRS.from_proj4(reconstruction.georef.proj4())
|
||||
to_srs = CRS.from_epsg(4326)
|
||||
transformer = location.transformer(from_srs, to_srs)
|
||||
# from_srs = CRS.from_proj4(reconstruction.georef.proj4())
|
||||
# to_srs = CRS.from_epsg(4326)
|
||||
# transformer = location.transformer(from_srs, to_srs)
|
||||
|
||||
for gcp in gcps:
|
||||
properties = gcp.copy()
|
||||
del properties['coordinates']
|
||||
# for gcp in gcps:
|
||||
# properties = gcp.copy()
|
||||
# del properties['coordinates']
|
||||
|
||||
geojson['features'].append({
|
||||
'type': 'Feature',
|
||||
'geometry': {
|
||||
'type': 'Point',
|
||||
'coordinates': transformer.TransformPoint(*gcp['coordinates']),
|
||||
},
|
||||
'properties': properties
|
||||
})
|
||||
# geojson['features'].append({
|
||||
# 'type': 'Feature',
|
||||
# 'geometry': {
|
||||
# 'type': 'Point',
|
||||
# 'coordinates': transformer.TransformPoint(*gcp['coordinates']),
|
||||
# },
|
||||
# 'properties': properties
|
||||
# })
|
||||
|
||||
with open(gcp_geojson_export_file, 'w') as f:
|
||||
f.write(json.dumps(geojson, indent=4))
|
||||
# with open(gcp_geojson_export_file, 'w') as f:
|
||||
# f.write(json.dumps(geojson, indent=4))
|
||||
|
||||
else:
|
||||
log.ODM_WARNING("GCPs could not be loaded for writing to %s" % gcp_export_file)
|
||||
# else:
|
||||
# log.ODM_WARNING("GCPs could not be loaded for writing to %s" % gcp_export_file)
|
||||
|
||||
if not io.file_exists(tree.odm_georeferencing_model_laz) or self.rerun():
|
||||
cmd = ('pdal translate -i "%s" -o \"%s\"' % (tree.filtered_point_cloud, tree.odm_georeferencing_model_laz))
|
||||
stages = ["ferry"]
|
||||
params = [
|
||||
'--filters.ferry.dimensions="views => UserData"',
|
||||
'--writers.las.compression="lazip"',
|
||||
]
|
||||
# cmd = ('pdal translate -i "%s" -o \"%s\"' % (tree.filtered_point_cloud, tree.odm_georeferencing_model_laz))
|
||||
# stages = ["ferry"]
|
||||
# params = [
|
||||
# '--filters.ferry.dimensions="views => UserData"',
|
||||
# '--writers.las.compression="lazip"',
|
||||
# ]
|
||||
|
||||
if reconstruction.is_georeferenced():
|
||||
log.ODM_INFO("Georeferencing point cloud")
|
||||
# if reconstruction.is_georeferenced():
|
||||
# log.ODM_INFO("Georeferencing point cloud")
|
||||
|
||||
stages.append("transformation")
|
||||
params += [
|
||||
'--filters.transformation.matrix="1 0 0 %s 0 1 0 %s 0 0 1 0 0 0 0 1"' % reconstruction.georef.utm_offset(),
|
||||
'--writers.las.offset_x=%s' % reconstruction.georef.utm_east_offset,
|
||||
'--writers.las.offset_y=%s' % reconstruction.georef.utm_north_offset,
|
||||
'--writers.las.scale_x=0.001',
|
||||
'--writers.las.scale_y=0.001',
|
||||
'--writers.las.scale_z=0.001',
|
||||
'--writers.las.offset_z=0',
|
||||
'--writers.las.a_srs="%s"' % reconstruction.georef.proj4()
|
||||
]
|
||||
# stages.append("transformation")
|
||||
# params += [
|
||||
# '--filters.transformation.matrix="1 0 0 %s 0 1 0 %s 0 0 1 0 0 0 0 1"' % reconstruction.georef.utm_offset(),
|
||||
# '--writers.las.offset_x=%s' % reconstruction.georef.utm_east_offset,
|
||||
# '--writers.las.offset_y=%s' % reconstruction.georef.utm_north_offset,
|
||||
# '--writers.las.scale_x=0.001',
|
||||
# '--writers.las.scale_y=0.001',
|
||||
# '--writers.las.scale_z=0.001',
|
||||
# '--writers.las.offset_z=0',
|
||||
# '--writers.las.a_srs="%s"' % reconstruction.georef.proj4()
|
||||
# ]
|
||||
|
||||
if reconstruction.has_gcp() and io.file_exists(gcp_gml_export_file):
|
||||
log.ODM_INFO("Embedding GCP info in point cloud")
|
||||
params += [
|
||||
'--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM_GCP\\\", \\\"description\\\": \\\"Ground Control Points (GML)\\\"}"' % gcp_gml_export_file.replace(os.sep, "/")
|
||||
]
|
||||
# if reconstruction.has_gcp() and io.file_exists(gcp_gml_export_file):
|
||||
# log.ODM_INFO("Embedding GCP info in point cloud")
|
||||
# params += [
|
||||
# '--writers.las.vlrs="{\\\"filename\\\": \\\"%s\\\", \\\"user_id\\\": \\\"ODM_GCP\\\", \\\"description\\\": \\\"Ground Control Points (GML)\\\"}"' % gcp_gml_export_file.replace(os.sep, "/")
|
||||
# ]
|
||||
|
||||
system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
|
||||
# system.run(cmd + ' ' + ' '.join(stages) + ' ' + ' '.join(params))
|
||||
|
||||
self.update_progress(50)
|
||||
# self.update_progress(50)
|
||||
|
||||
if args.crop > 0:
|
||||
log.ODM_INFO("Calculating cropping area and generating bounds shapefile from point cloud")
|
||||
cropper = Cropper(tree.odm_georeferencing, 'odm_georeferenced_model')
|
||||
# if args.crop > 0:
|
||||
# log.ODM_INFO("Calculating cropping area and generating bounds shapefile from point cloud")
|
||||
# cropper = Cropper(tree.odm_georeferencing, 'odm_georeferenced_model')
|
||||
|
||||
if args.fast_orthophoto:
|
||||
decimation_step = 4
|
||||
else:
|
||||
decimation_step = 40
|
||||
# if args.fast_orthophoto:
|
||||
# decimation_step = 4
|
||||
# else:
|
||||
# decimation_step = 40
|
||||
|
||||
# More aggressive decimation for large datasets
|
||||
if not args.fast_orthophoto:
|
||||
decimation_step *= int(len(reconstruction.photos) / 1000) + 1
|
||||
decimation_step = min(decimation_step, 95)
|
||||
# # More aggressive decimation for large datasets
|
||||
# if not args.fast_orthophoto:
|
||||
# decimation_step *= int(len(reconstruction.photos) / 1000) + 1
|
||||
# decimation_step = min(decimation_step, 95)
|
||||
|
||||
try:
|
||||
cropper.create_bounds_gpkg(tree.odm_georeferencing_model_laz, args.crop,
|
||||
decimation_step=decimation_step)
|
||||
except:
|
||||
log.ODM_WARNING("Cannot calculate crop bounds! We will skip cropping")
|
||||
args.crop = 0
|
||||
# try:
|
||||
# cropper.create_bounds_gpkg(tree.odm_georeferencing_model_laz, args.crop,
|
||||
# decimation_step=decimation_step)
|
||||
# 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")
|
||||
# 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))
|
||||
|
||||
# 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))
|
||||
|
||||
|
||||
if tree.odm_align_file is not None:
|
||||
tmp_dir = tree.path("odm_georeferencing", "codem")
|
||||
a_matrix = compute_alignment_matrix(tree.odm_georeferencing_model_laz, tree.odm_align_file, tmp_dir)
|
||||
if a_matrix is not None:
|
||||
print(a_matrix)
|
||||
exit(1)
|
||||
else:
|
||||
log.ODM_WARNING("Alignment to %s will be skipped." % tree.odm_align_file)
|
||||
# Align
|
||||
# - Compute alignment
|
||||
# - Transform point cloud
|
||||
# - Transform textured model(s)
|
||||
#
|
||||
|
||||
point_cloud.post_point_cloud_steps(args, tree, self.rerun())
|
||||
else:
|
||||
log.ODM_WARNING('Found a valid georeferenced model in: %s'
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue