kopia lustrzana https://github.com/OpenDroneMap/ODM
503 wiersze
20 KiB
Python
503 wiersze
20 KiB
Python
import os
|
|
import shutil
|
|
import warnings
|
|
import numpy as np
|
|
from opendm import get_image_size
|
|
from opendm import location
|
|
from opendm.gcp import GCPFile
|
|
from pyproj import CRS
|
|
import xmltodict as x2d
|
|
from six import string_types
|
|
|
|
from opendm import log
|
|
from opendm import io
|
|
from opendm import system
|
|
from opendm import context
|
|
from opendm import multispectral
|
|
|
|
from opendm.progress import progressbc
|
|
from opendm.photo import ODM_Photo
|
|
|
|
# Ignore warnings about proj information being lost
|
|
warnings.filterwarnings("ignore")
|
|
|
|
class ODM_Reconstruction(object):
|
|
def __init__(self, photos):
|
|
self.photos = photos
|
|
self.georef = None
|
|
self.gcp = None
|
|
self.multi_camera = self.detect_multi_camera()
|
|
self.filter_photos()
|
|
|
|
def detect_multi_camera(self):
|
|
"""
|
|
Looks at the reconstruction photos and determines if this
|
|
is a single or multi-camera setup.
|
|
"""
|
|
band_photos = {}
|
|
band_indexes = {}
|
|
|
|
for p in self.photos:
|
|
if not p.band_name in band_photos:
|
|
band_photos[p.band_name] = []
|
|
if not p.band_name in band_indexes:
|
|
band_indexes[p.band_name] = str(p.band_index)
|
|
|
|
band_photos[p.band_name].append(p)
|
|
|
|
bands_count = len(band_photos)
|
|
|
|
# Band name with the minimum number of photos
|
|
max_band_name = None
|
|
max_photos = -1
|
|
for band_name in band_photos:
|
|
if len(band_photos[band_name]) > max_photos:
|
|
max_band_name = band_name
|
|
max_photos = len(band_photos[band_name])
|
|
|
|
if bands_count >= 2 and bands_count <= 10:
|
|
# Validate that all bands have the same number of images,
|
|
# otherwise this is not a multi-camera setup
|
|
img_per_band = len(band_photos[max_band_name])
|
|
|
|
mc = []
|
|
for band_name in band_indexes:
|
|
mc.append({'name': band_name, 'photos': band_photos[band_name]})
|
|
|
|
filter_missing = False
|
|
for band in band_photos:
|
|
if len(band_photos[band]) < img_per_band:
|
|
log.ODM_WARNING("Multi-camera setup detected, but band \"%s\" (identified from \"%s\") has only %s images (instead of %s), perhaps images are missing or are corrupted." % (band, band_photos[band][0].filename, len(band_photos[band]), len(band_photos[max_band_name])))
|
|
filter_missing = True
|
|
|
|
if filter_missing:
|
|
# Calculate files to ignore
|
|
_, p2s = multispectral.compute_band_maps(mc, max_band_name)
|
|
|
|
max_files_per_band = 0
|
|
|
|
for filename in p2s:
|
|
max_files_per_band = max(max_files_per_band, len(p2s[filename]))
|
|
|
|
for filename in p2s:
|
|
if len(p2s[filename]) < max_files_per_band:
|
|
photos_to_remove = p2s[filename] + [p for p in self.photos if p.filename == filename]
|
|
for photo in photos_to_remove:
|
|
log.ODM_WARNING("Excluding %s" % photo.filename)
|
|
|
|
self.photos = [p for p in self.photos if p != photo]
|
|
for i in range(len(mc)):
|
|
mc[i]['photos'] = [p for p in mc[i]['photos'] if p != photo]
|
|
|
|
log.ODM_INFO("New image count: %s" % len(self.photos))
|
|
|
|
# We enforce a normalized band order for all bands that we can identify
|
|
# and rely on the manufacturer's band_indexes as a fallback for all others
|
|
normalized_band_order = {
|
|
'RGB': '0',
|
|
'REDGREENBLUE': '0',
|
|
|
|
'RED': '1',
|
|
'R': '1',
|
|
|
|
'GREEN': '2',
|
|
'G': '2',
|
|
|
|
'BLUE': '3',
|
|
'B': '3',
|
|
|
|
'NIR': '4',
|
|
'N': '4',
|
|
|
|
'REDEDGE': '5',
|
|
'RE': '5',
|
|
|
|
'PANCHRO': '6',
|
|
|
|
'LWIR': '7',
|
|
'L': '7',
|
|
}
|
|
|
|
for band_name in band_indexes:
|
|
if band_name.upper() not in normalized_band_order:
|
|
log.ODM_WARNING(f"Cannot identify order for {band_name} band, using manufacturer suggested index instead")
|
|
|
|
# Sort
|
|
mc.sort(key=lambda x: normalized_band_order.get(x['name'].upper(), '9' + band_indexes[x['name']]))
|
|
|
|
for c, d in enumerate(mc):
|
|
log.ODM_INFO(f"Band {c + 1}: {d['name']}")
|
|
|
|
return mc
|
|
|
|
return None
|
|
|
|
def filter_photos(self):
|
|
if not self.multi_camera:
|
|
return # Nothing to do, use all images
|
|
|
|
else:
|
|
# Sometimes people might try process both RGB + Blue/Red/Green bands
|
|
# because these are the contents of the SD card from a drone (e.g. DJI P4 Multispectral)
|
|
# But we don't want to process both, so we discard the RGB files in favor
|
|
bands = {}
|
|
for b in self.multi_camera:
|
|
bands[b['name'].lower()] = b['name']
|
|
|
|
bands_to_remove = []
|
|
|
|
if 'rgb' in bands or 'redgreenblue' in bands:
|
|
if 'red' in bands and 'green' in bands and 'blue' in bands:
|
|
bands_to_remove.append(bands['rgb'] if 'rgb' in bands else bands['redgreenblue'])
|
|
|
|
# Mavic 3M's RGB camera lens are too different than the multispectral ones
|
|
# so we drop the RGB channel instead
|
|
elif self.photos[0].is_make_model("DJI", "M3M") and 'red' in bands and 'green' in bands:
|
|
bands_to_remove.append(bands['rgb'] if 'rgb' in bands else bands['redgreenblue'])
|
|
|
|
else:
|
|
for b in ['red', 'green', 'blue']:
|
|
if b in bands:
|
|
bands_to_remove.append(bands[b])
|
|
|
|
if len(bands_to_remove) > 0:
|
|
log.ODM_WARNING("Redundant bands detected, probably because RGB images are mixed with single band images. We will trim some bands as needed")
|
|
|
|
for band_to_remove in bands_to_remove:
|
|
self.multi_camera = [b for b in self.multi_camera if b['name'] != band_to_remove]
|
|
photos_before = len(self.photos)
|
|
self.photos = [p for p in self.photos if p.band_name != band_to_remove]
|
|
photos_after = len(self.photos)
|
|
|
|
log.ODM_WARNING("Skipping %s band (%s images)" % (band_to_remove, photos_before - photos_after))
|
|
|
|
def is_georeferenced(self):
|
|
return self.georef is not None
|
|
|
|
def has_gcp(self):
|
|
return self.is_georeferenced() and self.gcp is not None and self.gcp.exists()
|
|
|
|
def has_geotagged_photos(self):
|
|
for photo in self.photos:
|
|
if photo.latitude is None and photo.longitude is None:
|
|
return False
|
|
|
|
return True
|
|
|
|
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:
|
|
gcp = GCPFile(gcp_file)
|
|
if gcp.exists():
|
|
if gcp.entries_count() == 0:
|
|
raise RuntimeError("This GCP file does not have any entries. Are the entries entered in the proper format?")
|
|
|
|
gcp.check_entries()
|
|
|
|
# Convert GCP file to a UTM projection since the rest of the pipeline
|
|
# does not handle other SRS well.
|
|
rejected_entries = []
|
|
utm_gcp = GCPFile(gcp.create_utm_copy(output_gcp_file, filenames=[p.filename for p in self.photos], rejected_entries=rejected_entries, include_extras=True))
|
|
|
|
if not utm_gcp.exists():
|
|
raise RuntimeError("Could not project GCP file to UTM. Please double check your GCP file for mistakes.")
|
|
|
|
for re in rejected_entries:
|
|
log.ODM_WARNING("GCP line ignored (image not found): %s" % str(re))
|
|
|
|
if utm_gcp.entries_count() > 0:
|
|
log.ODM_INFO("%s GCP points will be used for georeferencing" % utm_gcp.entries_count())
|
|
else:
|
|
raise RuntimeError("A GCP file was provided, but no valid GCP entries could be used. Note that the GCP file is case sensitive (\".JPG\" is not the same as \".jpg\").")
|
|
|
|
self.gcp = utm_gcp
|
|
|
|
# Compute RTC offsets from GCP points
|
|
x_pos = [p.x for p in utm_gcp.iter_entries()]
|
|
y_pos = [p.y for p in utm_gcp.iter_entries()]
|
|
x_off, y_off = int(np.round(np.mean(x_pos))), int(np.round(np.mean(y_pos)))
|
|
|
|
# Create coords file, we'll be using this later
|
|
# during georeferencing
|
|
with open(output_coords_file, 'w') as f:
|
|
coords_header = gcp.wgs84_utm_zone()
|
|
f.write(coords_header + "\n")
|
|
f.write("{} {}\n".format(x_off, y_off))
|
|
log.ODM_INFO("Generated coords file from GCP: %s" % coords_header)
|
|
|
|
# Deprecated: This is mostly for backward compatibility and should be
|
|
# be removed at some point
|
|
shutil.copyfile(output_coords_file, output_model_txt_geo)
|
|
log.ODM_INFO("Wrote %s" % output_model_txt_geo)
|
|
else:
|
|
log.ODM_WARNING("GCP file does not exist: %s" % gcp_file)
|
|
return
|
|
else:
|
|
log.ODM_INFO("Coordinates file already exist: %s" % output_coords_file)
|
|
log.ODM_INFO("GCP file already exist: %s" % output_gcp_file)
|
|
self.gcp = GCPFile(output_gcp_file)
|
|
|
|
self.georef = ODM_GeoRef.FromCoordsFile(output_coords_file)
|
|
return self.georef
|
|
|
|
def georeference_with_gps(self, images_path, output_coords_file, output_model_txt_geo, rerun=False):
|
|
try:
|
|
if not io.file_exists(output_coords_file) or rerun:
|
|
location.extract_utm_coords(self.photos, images_path, output_coords_file)
|
|
else:
|
|
log.ODM_INFO("Coordinates file already exist: %s" % output_coords_file)
|
|
|
|
# Deprecated: This is mostly for backward compatibility and should be
|
|
# be removed at some point
|
|
if not io.file_exists(output_model_txt_geo) or rerun:
|
|
with open(output_coords_file, 'r') as f:
|
|
with open(output_model_txt_geo, 'w+') as w:
|
|
w.write(f.readline()) # CRS
|
|
w.write(f.readline()) # Offset
|
|
else:
|
|
log.ODM_INFO("Model geo file already exist: %s" % output_model_txt_geo)
|
|
|
|
self.georef = ODM_GeoRef.FromCoordsFile(output_coords_file)
|
|
except:
|
|
log.ODM_WARNING('Could not generate coordinates file. The orthophoto will not be georeferenced.')
|
|
|
|
self.gcp = GCPFile(None)
|
|
return self.georef
|
|
|
|
def save_proj_srs(self, file):
|
|
# Save proj to file for future use (unless this
|
|
# dataset is not georeferenced)
|
|
if self.is_georeferenced():
|
|
with open(file, 'w') as f:
|
|
f.write(self.get_proj_srs())
|
|
|
|
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:
|
|
if p.filename == filename:
|
|
return p
|
|
|
|
|
|
class ODM_GeoRef(object):
|
|
@staticmethod
|
|
def FromCoordsFile(coords_file):
|
|
# check for coordinate file existence
|
|
if not io.file_exists(coords_file):
|
|
log.ODM_WARNING('Could not find file %s' % coords_file)
|
|
return
|
|
|
|
srs = None
|
|
utm_east_offset = None
|
|
utm_north_offset = None
|
|
|
|
with open(coords_file) as f:
|
|
# extract reference system and utm zone from first line.
|
|
# We will assume the following format:
|
|
# 'WGS84 UTM 17N' or 'WGS84 UTM 17N \n'
|
|
line = f.readline().rstrip()
|
|
srs = location.parse_srs_header(line)
|
|
|
|
# second line is a northing/easting offset
|
|
line = f.readline().rstrip()
|
|
utm_east_offset, utm_north_offset = map(float, line.split(" "))
|
|
|
|
return ODM_GeoRef(srs, utm_east_offset, utm_north_offset)
|
|
|
|
def __init__(self, srs, utm_east_offset, utm_north_offset):
|
|
self.srs = srs
|
|
self.utm_east_offset = utm_east_offset
|
|
self.utm_north_offset = utm_north_offset
|
|
self.transform = []
|
|
|
|
def proj4(self):
|
|
return self.srs.to_proj4()
|
|
|
|
def utm_offset(self):
|
|
return (self.utm_east_offset, self.utm_north_offset)
|
|
|
|
class ODM_Tree(object):
|
|
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')
|
|
|
|
# modules paths
|
|
|
|
# here are defined where all modules should be located in
|
|
# order to keep track all files al directories during the
|
|
# whole reconstruction process.
|
|
self.dataset_raw = os.path.join(self.root_path, 'images')
|
|
self.opensfm = os.path.join(self.root_path, 'opensfm')
|
|
self.openmvs = os.path.join(self.opensfm, 'undistorted', 'openmvs')
|
|
self.odm_meshing = os.path.join(self.root_path, 'odm_meshing')
|
|
self.odm_texturing = os.path.join(self.root_path, 'odm_texturing')
|
|
self.odm_25dtexturing = os.path.join(self.root_path, 'odm_texturing_25d')
|
|
self.odm_georeferencing = os.path.join(self.root_path, 'odm_georeferencing')
|
|
self.odm_filterpoints = os.path.join(self.root_path, 'odm_filterpoints')
|
|
self.odm_orthophoto = os.path.join(self.root_path, 'odm_orthophoto')
|
|
self.odm_report = os.path.join(self.root_path, 'odm_report')
|
|
|
|
# important files paths
|
|
|
|
# benchmarking
|
|
self.benchmarking = os.path.join(self.root_path, 'benchmark.txt')
|
|
self.dataset_list = os.path.join(self.root_path, 'img_list.txt')
|
|
|
|
# opensfm
|
|
self.opensfm_image_list = os.path.join(self.opensfm, 'image_list.txt')
|
|
self.opensfm_reconstruction = os.path.join(self.opensfm, 'reconstruction.json')
|
|
self.opensfm_reconstruction_nvm = os.path.join(self.opensfm, 'undistorted/reconstruction.nvm')
|
|
self.opensfm_geocoords_reconstruction = os.path.join(self.opensfm, 'reconstruction.geocoords.json')
|
|
self.opensfm_topocentric_reconstruction = os.path.join(self.opensfm, 'reconstruction.topocentric.json')
|
|
|
|
# OpenMVS
|
|
self.openmvs_model = os.path.join(self.openmvs, 'scene_dense_dense_filtered.ply')
|
|
|
|
# filter points
|
|
self.filtered_point_cloud = os.path.join(self.odm_filterpoints, "point_cloud.ply")
|
|
self.filtered_point_cloud_stats = os.path.join(self.odm_filterpoints, "point_cloud_stats.json")
|
|
|
|
# odm_meshing
|
|
self.odm_mesh = os.path.join(self.odm_meshing, 'odm_mesh.ply')
|
|
self.odm_meshing_log = os.path.join(self.odm_meshing, 'odm_meshing_log.txt')
|
|
self.odm_25dmesh = os.path.join(self.odm_meshing, 'odm_25dmesh.ply')
|
|
self.odm_25dmeshing_log = os.path.join(self.odm_meshing, 'odm_25dmeshing_log.txt')
|
|
|
|
# texturing
|
|
self.odm_textured_model_obj = 'odm_textured_model_geo.obj'
|
|
self.odm_textured_model_glb = 'odm_textured_model_geo.glb'
|
|
|
|
# odm_georeferencing
|
|
self.odm_georeferencing_coords = os.path.join(
|
|
self.odm_georeferencing, 'coords.txt')
|
|
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) or io.find('align.las', self.root_path) or io.find('align.tif', self.root_path)
|
|
|
|
self.odm_georeferencing_proj = 'proj.txt'
|
|
self.odm_georeferencing_model_txt_geo = os.path.join(
|
|
self.odm_georeferencing, 'odm_georeferencing_model_geo.txt')
|
|
self.odm_georeferencing_xyz_file = os.path.join(
|
|
self.odm_georeferencing, 'odm_georeferenced_model.csv')
|
|
self.odm_georeferencing_model_laz = os.path.join(
|
|
self.odm_georeferencing, 'odm_georeferenced_model.laz')
|
|
self.odm_georeferencing_model_las = os.path.join(
|
|
self.odm_georeferencing, 'odm_georeferenced_model.las')
|
|
self.odm_georeferencing_alignment_matrix = os.path.join(
|
|
self.odm_georeferencing, 'alignment_matrix.json'
|
|
)
|
|
|
|
# odm_orthophoto
|
|
self.odm_orthophoto_render = os.path.join(self.odm_orthophoto, 'odm_orthophoto_render.tif')
|
|
self.odm_orthophoto_tif = os.path.join(self.odm_orthophoto, 'odm_orthophoto.tif')
|
|
self.odm_orthophoto_corners = os.path.join(self.odm_orthophoto, 'odm_orthophoto_corners.txt')
|
|
self.odm_orthophoto_log = os.path.join(self.odm_orthophoto, 'odm_orthophoto_log.txt')
|
|
self.odm_orthophoto_tif_log = os.path.join(self.odm_orthophoto, 'gdal_translate_log.txt')
|
|
|
|
# tiles
|
|
self.orthophoto_tiles = os.path.join(self.root_path, "orthophoto_tiles")
|
|
|
|
# Split-merge
|
|
self.submodels_path = os.path.join(self.root_path, 'submodels')
|
|
|
|
# Tiles
|
|
self.entwine_pointcloud = self.path("entwine_pointcloud")
|
|
self.ogc_tiles = self.path("3d_tiles")
|
|
|
|
def path(self, *args):
|
|
return os.path.join(self.root_path, *args)
|
|
|
|
|
|
class ODM_Stage:
|
|
def __init__(self, name, args, progress=0.0, **params):
|
|
self.name = name
|
|
self.args = args
|
|
self.progress = progress
|
|
self.params = params
|
|
if self.params is None:
|
|
self.params = {}
|
|
self.next_stage = None
|
|
self.prev_stage = None
|
|
|
|
def connect(self, stage):
|
|
self.next_stage = stage
|
|
stage.prev_stage = self
|
|
return stage
|
|
|
|
def rerun(self):
|
|
"""
|
|
Does this stage need to be rerun?
|
|
"""
|
|
return (self.args.rerun is not None and self.args.rerun == self.name) or \
|
|
(self.args.rerun_all) or \
|
|
(self.args.rerun_from is not None and self.name in self.args.rerun_from)
|
|
|
|
def run(self, outputs = {}):
|
|
start_time = system.now_raw()
|
|
log.logger.log_json_stage_run(self.name, start_time)
|
|
|
|
log.ODM_INFO('Running %s stage' % self.name)
|
|
|
|
self.process(self.args, outputs)
|
|
|
|
# The tree variable should always be populated at this point
|
|
if outputs.get('tree') is None:
|
|
raise Exception("Assert violation: tree variable is missing from outputs dictionary.")
|
|
|
|
try:
|
|
system.benchmark(start_time, outputs['tree'].benchmarking, self.name)
|
|
except Exception as e:
|
|
log.ODM_WARNING("Cannot write benchmark file: %s" % str(e))
|
|
|
|
log.ODM_INFO('Finished %s stage' % self.name)
|
|
self.update_progress_end()
|
|
|
|
# Last stage?
|
|
if self.args.end_with == self.name or self.args.rerun == self.name:
|
|
log.ODM_INFO("No more stages to run")
|
|
return
|
|
|
|
# Run next stage?
|
|
elif self.next_stage is not None:
|
|
self.next_stage.run(outputs)
|
|
|
|
def delta_progress(self):
|
|
if self.prev_stage:
|
|
return max(0.0, self.progress - self.prev_stage.progress)
|
|
else:
|
|
return max(0.0, self.progress)
|
|
|
|
def previous_stages_progress(self):
|
|
if self.prev_stage:
|
|
return max(0.0, self.prev_stage.progress)
|
|
else:
|
|
return 0.0
|
|
|
|
def update_progress_end(self):
|
|
self.update_progress(100.0)
|
|
|
|
def update_progress(self, progress):
|
|
progress = max(0.0, min(100.0, progress))
|
|
progressbc.send_update(self.previous_stages_progress() +
|
|
(self.delta_progress() / 100.0) * float(progress))
|
|
|
|
def last_stage(self):
|
|
if self.next_stage:
|
|
return self.next_stage.last_stage()
|
|
else:
|
|
return self
|
|
|
|
|
|
def process(self, args, outputs):
|
|
raise NotImplementedError
|
|
|