diff --git a/SuperBuild/cmake/External-OpenMVS.cmake b/SuperBuild/cmake/External-OpenMVS.cmake index b8cc89ac..2a08de76 100644 --- a/SuperBuild/cmake/External-OpenMVS.cmake +++ b/SuperBuild/cmake/External-OpenMVS.cmake @@ -20,7 +20,7 @@ ExternalProject_Add(${_proj_name} #--Download step-------------- DOWNLOAD_DIR ${SB_DOWNLOAD_DIR} GIT_REPOSITORY https://github.com/OpenDroneMap/openMVS - GIT_TAG 210 + GIT_TAG 230 #--Update/Patch step---------- UPDATE_COMMAND "" #--Configure step------------- diff --git a/opendm/concurrency.py b/opendm/concurrency.py index 910cf6be..d09555fa 100644 --- a/opendm/concurrency.py +++ b/opendm/concurrency.py @@ -24,7 +24,7 @@ def get_max_memory_mb(minimum = 100, use_at_most = 0.5): """ return max(minimum, (virtual_memory().available / 1024 / 1024) * use_at_most) -def parallel_map(func, items, max_workers=1): +def parallel_map(func, items, max_workers=1, single_thread_fallback=True): """ Our own implementation for parallel processing which handles gracefully CTRL+C and reverts to @@ -85,7 +85,7 @@ def parallel_map(func, items, max_workers=1): stop_workers() - if error is not None: + if error is not None and single_thread_fallback: # Try to reprocess using a single thread # in case this was a memory error log.ODM_WARNING("Failed to run process in parallel, retrying with a single thread...") diff --git a/opendm/multispectral.py b/opendm/multispectral.py index 9f3aca3e..9ccbec83 100644 --- a/opendm/multispectral.py +++ b/opendm/multispectral.py @@ -5,6 +5,7 @@ import os from opendm import dls import numpy as np from opendm import log +from opendm.concurrency import parallel_map from opensfm.io import imread from skimage import exposure @@ -270,7 +271,7 @@ def compute_band_maps(multi_camera, primary_band): return s2p, p2s -def compute_alignment_matrices(multi_camera, primary_band_name, images_path, s2p, p2s, max_samples=9999): +def compute_alignment_matrices(multi_camera, primary_band_name, images_path, s2p, p2s, max_concurrency=1, max_samples=30): log.ODM_INFO("Computing band alignment") alignment_info = {} @@ -280,93 +281,129 @@ def compute_alignment_matrices(multi_camera, primary_band_name, images_path, s2p if band['name'] != primary_band_name: matrices = [] - # if band['name'] != "NIR": - # continue # TODO REMOVE + def parallel_compute_homography(p): + try: + if len(matrices) >= max_samples: + log.ODM_INFO("Got enough samples for %s (%s)" % (band['name'], max_samples)) + return - # Find good matrix candidates for alignment - for p in band['photos']: - primary_band_photo = s2p.get(p.filename) - if primary_band_photo is None: - log.ODM_WARNING("Cannot find primary band photo for %s" % p.filename) - continue + # Find good matrix candidates for alignment - warp_matrix, score, dimension = compute_homography(os.path.join(images_path, p.filename), - os.path.join(images_path, primary_band_photo.filename)) + primary_band_photo = s2p.get(p['filename']) + if primary_band_photo is None: + log.ODM_WARNING("Cannot find primary band photo for %s" % p['filename']) + return - if warp_matrix is not None: - log.ODM_INFO("%s --> %s good match (score: %s)" % (p.filename, primary_band_photo.filename, score)) - matrices.append({ - 'warp_matrix': warp_matrix, - 'score': score, - 'dimension': dimension - }) - else: - log.ODM_INFO("%s --> %s cannot be matched" % (p.filename, primary_band_photo.filename)) + warp_matrix, dimension, algo = compute_homography(os.path.join(images_path, p['filename']), + os.path.join(images_path, primary_band_photo.filename)) - if len(matrices) >= max_samples: - log.ODM_INFO("Got enough samples for %s (%s)" % (band['name'], max_samples)) - break + if warp_matrix is not None: + log.ODM_INFO("%s --> %s good match" % (p['filename'], primary_band_photo.filename)) + + matrices.append({ + 'warp_matrix': warp_matrix, + 'eigvals': np.linalg.eigvals(warp_matrix), + 'dimension': dimension, + 'algo': algo + }) + else: + log.ODM_INFO("%s --> %s cannot be matched" % (p['filename'], primary_band_photo.filename)) + except Exception as e: + log.ODM_WARNING("Failed to compute homography for %s: %s" % (p['filename'], str(e))) + + parallel_map(parallel_compute_homography, [{'filename': p.filename} for p in band['photos']], max_concurrency, single_thread_fallback=False) + + # Choose winning algorithm (doesn't seem to yield improvements) + # feat_count = 0 + # ecc_count = 0 + # for m in matrices: + # if m['algo'] == 'feat': + # feat_count += 1 + # if m['algo'] == 'ecc': + # ecc_count += 1 + + # algo = 'feat' if feat_count >= ecc_count else 'ecc' + + # log.ODM_INFO("Feat: %s | ECC: %s | Winner: %s" % (feat_count, ecc_count, algo)) + # matrices = [m for m in matrices if m['algo'] == algo] + + # Find the matrix that has the most common eigvals + # among all matrices. That should be the "best" alignment. + for m1 in matrices: + acc = np.array([0.0,0.0,0.0]) + e = m1['eigvals'] + + for m2 in matrices: + acc += abs(e - m2['eigvals']) + + m1['score'] = acc.sum() # Sort matrices.sort(key=lambda x: x['score'], reverse=False) if len(matrices) > 0: alignment_info[band['name']] = matrices[0] - print(matrices[0]) + log.ODM_INFO("%s band will be aligned using warp matrix %s (score: %s)" % (band['name'], matrices[0]['warp_matrix'], matrices[0]['score'])) else: log.ODM_WARNING("Cannot find alignment matrix for band %s, The band will likely be misaligned!" % band['name']) return alignment_info def compute_homography(image_filename, align_image_filename): - # try: - # Convert images to grayscale if needed - image = imread(image_filename, unchanged=True, anydepth=True) - if image.shape[2] == 3: - image_gray = to_8bit(cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)) - else: - image_gray = to_8bit(image[:,:,0]) + try: + # Convert images to grayscale if needed + image = imread(image_filename, unchanged=True, anydepth=True) + if image.shape[2] == 3: + image_gray = to_8bit(cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)) + else: + image_gray = to_8bit(image[:,:,0]) - align_image = imread(align_image_filename, unchanged=True, anydepth=True) - if align_image.shape[2] == 3: - align_image_gray = to_8bit(cv2.cvtColor(align_image, cv2.COLOR_BGR2GRAY)) - else: - align_image_gray = to_8bit(align_image[:,:,0]) + align_image = imread(align_image_filename, unchanged=True, anydepth=True) + if align_image.shape[2] == 3: + align_image_gray = to_8bit(cv2.cvtColor(align_image, cv2.COLOR_BGR2GRAY)) + else: + align_image_gray = to_8bit(align_image[:,:,0]) - def compute_using(algorithm): - h = algorithm(image_gray, align_image_gray) - if h is None: - return None, None, (None, None) + def compute_using(algorithm): + h = algorithm(image_gray, align_image_gray) + if h is None: + return None, (None, None) - det = np.linalg.det(h) + det = np.linalg.det(h) + + # Check #1 homography's determinant will not be close to zero + if abs(det) < 0.25: + return None, (None, None) + + # Check #2 the ratio of the first-to-last singular value is sane (not too high) + svd = np.linalg.svd(h, compute_uv=False) + if svd[-1] == 0: + return None, (None, None) + + ratio = svd[0] / svd[-1] + if ratio > 100000: + return None, (None, None) + + return h, (align_image_gray.shape[1], align_image_gray.shape[0]) - # Check #1 homography's determinant will not be close to zero - if abs(det) < 0.25: - return None, None, (None, None) + algo = 'feat' + result = compute_using(find_features_homography) - # Check #2 the ratio of the first-to-last singular value is sane (not too high) - svd = np.linalg.svd(h, compute_uv=False) - if svd[-1] == 0: - return None, None, (None, None) + if result[0] is None: + algo = 'ecc' + log.ODM_INFO("Can't use features matching, will use ECC (this might take a bit)") + result = compute_using(find_ecc_homography) + if result[0] is None: + algo = None - ratio = svd[0] / svd[-1] - if ratio > 100000: - return None, None, (None, None) + warp_matrix, dimension = result + return warp_matrix, dimension, algo - return h, compute_alignment_score(h, image_gray, align_image_gray), (align_image_gray.shape[1], align_image_gray.shape[0]) - - result = compute_using(find_features_homography) - if result[0] is None: - log.ODM_INFO("Can't use features matching, will use ECC") - result = compute_using(find_ecc_homography) - - return result + except Exception as e: + log.ODM_WARNING("Compute homography: %s" % str(e)) + return None, None, (None, None) - # except Exception as e: - # log.ODM_WARNING("Compute homography: %s" % str(e)) - # return None, None, (None, None) - -def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=5000, termination_eps=1e-8): +def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=2500, termination_eps=1e-9): image_gray = to_8bit(gradient(gaussian(image_gray))) align_image_gray = to_8bit(gradient(gaussian(align_image_gray))) @@ -377,10 +414,7 @@ def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=5000, criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps) - try: - (cc, warp_matrix) = cv2.findTransformECC (image_gray,align_image_gray,warp_matrix, cv2.MOTION_HOMOGRAPHY, criteria, inputMask=None, gaussFiltSize=1) - except: - (cc, warp_matrix) = cv2.findTransformECC (image_gray,align_image_gray,warp_matrix, cv2.MOTION_HOMOGRAPHY, criteria) + _, warp_matrix = cv2.findTransformECC (image_gray,align_image_gray,warp_matrix, cv2.MOTION_HOMOGRAPHY, criteria, inputMask=None, gaussFiltSize=9) return warp_matrix @@ -418,37 +452,6 @@ def find_features_homography(image_gray, align_image_gray, feature_retention=0.2 h, _ = cv2.findHomography(points_image, points_align_image, cv2.RANSAC) return h -def compute_alignment_score(warp_matrix, image_gray, align_image_gray, apply_gradient=True): - projected = align_image(image_gray, warp_matrix, (align_image_gray.shape[1], align_image_gray.shape[0])) - borders = projected==0 - - if apply_gradient: - image_gray = to_8bit(gradient(gaussian(image_gray))) - align_image_gray = to_8bit(gradient(gaussian(align_image_gray))) - - # cv2.imwrite("/datasets/micasense/opensfm/undistorted/align_image_gray.jpg", align_image_gray) - # cv2.imwrite("/datasets/micasense/opensfm/undistorted/projected.jpg", projected) - - # Threshold - align_image_gray[align_image_gray > 128] = 255 - projected[projected > 128] = 255 - align_image_gray[align_image_gray <= 128] = 0 - projected[projected <= 128] = 0 - - # Mark borders - align_image_gray[borders] = 0 - projected[borders] = 255 - - - # cv2.imwrite("/datasets/micasense/opensfm/undistorted/threshold_align_image_gray.jpg", align_image_gray) - # cv2.imwrite("/datasets/micasense/opensfm/undistorted/threshold_projected.jpg", projected) - - # cv2.imwrite("/datasets/micasense/opensfm/undistorted/delta.jpg", projected - align_image_gray) - - # Compute delta --> the more the images overlap perfectly, the lower the score - return (projected - align_image_gray).sum() - - def gradient(im, ksize=5): im = local_normalize(im) grad_x = cv2.Sobel(im,cv2.CV_32F,1,0,ksize=ksize) diff --git a/opendm/osfm.py b/opendm/osfm.py index 94db2397..891e38e0 100644 --- a/opendm/osfm.py +++ b/opendm/osfm.py @@ -315,11 +315,11 @@ class OSFMContext: else: log.ODM_INFO("Already extracted cameras") - def convert_and_undistort(self, rerun=False, imageFilter=None, image_list=None): + def convert_and_undistort(self, rerun=False, imageFilter=None, image_list=None, runId="nominal"): log.ODM_INFO("Undistorting %s ..." % self.opensfm_project_path) - undistorted_images_path = self.path("undistorted", "images") + done_flag_file = self.path("undistorted", "%s_done.txt" % runId) - if not io.dir_exists(undistorted_images_path) or rerun: + if not io.file_exists(done_flag_file) or rerun: ds = DataSet(self.opensfm_project_path) if image_list is not None: @@ -327,32 +327,35 @@ class OSFMContext: undistort.run_dataset(ds, "reconstruction.json", 0, None, "undistorted", imageFilter) + + self.touch(done_flag_file) else: - log.ODM_WARNING("Found an undistorted directory in %s" % undistorted_images_path) + log.ODM_WARNING("Already undistorted (%s)" % runId) + + def restore_reconstruction_backup(self): + if os.path.exists(self.recon_backup_file()): + # This time export the actual reconstruction.json + # (containing only the primary band) + if os.path.exists(self.recon_file()): + os.remove(self.recon_file()) + os.rename(self.recon_backup_file(), self.recon_file()) + log.ODM_INFO("Restored reconstruction.json") + + def backup_reconstruction(self): + if os.path.exists(self.recon_backup_file()): + os.remove(self.recon_backup_file()) + + log.ODM_INFO("Backing up reconstruction") + shutil.copyfile(self.recon_file(), self.recon_backup_file()) + + def recon_backup_file(self): + return self.path("reconstruction.backup.json") + + def recon_file(self): + return self.path("reconstruction.json") def add_shots_to_reconstruction(self, p2s): - recon_file = self.path("reconstruction.json") - recon_backup_file = self.path("reconstruction.backup.json") - # tracks_file = self.path("tracks.csv") - # tracks_backup_file = self.path("tracks.backup.csv") - # image_list_file = self.path("image_list.txt") - # image_list_backup_file = self.path("image_list.backup.txt") - - log.ODM_INFO("Adding shots to reconstruction") - if os.path.exists(recon_backup_file): - os.remove(recon_backup_file) - # if os.path.exists(tracks_backup_file): - # os.remove(tracks_backup_file) - # if os.path.exists(image_list_backup_file): - # os.remove(image_list_backup_file) - - log.ODM_INFO("Backing up reconstruction, tracks, image list") - shutil.copyfile(recon_file, recon_backup_file) - - # shutil.copyfile(tracks_file, tracks_backup_file) - # shutil.copyfile(image_list_file, image_list_backup_file) - - with open(recon_file) as f: + with open(self.recon_file()) as f: reconstruction = json.loads(f.read()) # Augment reconstruction.json @@ -369,13 +372,9 @@ class OSFMContext: for p in secondary_photos: shots[p.filename] = shots[shot_id] - with open(recon_file, 'w') as f: + with open(self.recon_file(), 'w') as f: f.write(json.dumps(reconstruction)) - return True #(recon_file, recon_backup_file) - # (tracks_file, tracks_backup_file), - # (image_list_file, image_list_backup_file)] - def update_config(self, cfg_dict): cfg_file = self.get_config_file_path() diff --git a/requirements.txt b/requirements.txt index 15bd2a9d..a082dafd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -23,6 +23,7 @@ PyYAML==5.1 rasterio==1.1.8 repoze.lru==0.7 scikit-learn==0.23.2 +scikit-image==0.17.2 scipy==1.5.4 xmltodict==0.12.0 diff --git a/stages/run_opensfm.py b/stages/run_opensfm.py index 3a1635a6..e1021453 100644 --- a/stages/run_opensfm.py +++ b/stages/run_opensfm.py @@ -91,7 +91,6 @@ class ODMOpenSfMStage(types.ODM_Stage): log.ODM_WARNING("Cannot align %s, no alignment matrix could be computed. Band alignment quality might be affected." % (shot_id)) return image - if args.radiometric_calibration != "none": undistort_pipeline.append(radiometric_calibrate) @@ -109,16 +108,16 @@ class ODMOpenSfMStage(types.ODM_Stage): # We finally restore the original files later added_shots_file = octx.path('added_shots_done.txt') - + if not io.file_exists(added_shots_file) or self.rerun(): primary_band_name = multispectral.get_primary_band_name(reconstruction.multi_camera, args.primary_band) s2p, p2s = multispectral.compute_band_maps(reconstruction.multi_camera, primary_band_name) - alignment_info = multispectral.compute_alignment_matrices(reconstruction.multi_camera, primary_band_name, tree.dataset_raw, s2p, p2s) + alignment_info = multispectral.compute_alignment_matrices(reconstruction.multi_camera, primary_band_name, tree.dataset_raw, s2p, p2s, max_concurrency=args.max_concurrency) + + log.ODM_INFO("Adding shots to reconstruction") + + octx.backup_reconstruction() octx.add_shots_to_reconstruction(p2s) - - # TODO: what happens to reconstruction.backup.json - # if the process fails here? - octx.touch(added_shots_file) undistort_pipeline.append(align_to_primary_band) @@ -148,16 +147,8 @@ class ODMOpenSfMStage(types.ODM_Stage): else: log.ODM_WARNING("Found a valid NVM file in %s for %s band" % (nvm_file, band['name'])) - recon_file = octx.path("reconstruction.json") - recon_backup_file = octx.path("reconstruction.backup.json") - - if os.path.exists(recon_backup_file): - # This time export the actual reconstruction.json - # (containing only the primary band) - if os.path.exists(recon_backup_file): - os.remove(recon_file) - os.rename(recon_backup_file, recon_file) - octx.convert_and_undistort(self.rerun(), undistort_callback) + octx.restore_reconstruction_backup() + octx.convert_and_undistort(self.rerun(), undistort_callback, runId='primary') if not io.file_exists(tree.opensfm_reconstruction_nvm) or self.rerun(): octx.run('export_visualsfm --points') @@ -194,6 +185,9 @@ class ODMOpenSfMStage(types.ODM_Stage): if args.optimize_disk_space: os.remove(octx.path("tracks.csv")) + if io.file_exists(octx.recon_backup_file()): + os.remove(octx.recon_backup_file()) + if io.dir_exists(octx.path("undistorted", "depthmaps")): files = glob.glob(octx.path("undistorted", "depthmaps", "*.npz")) for f in files: