PoC ECC and homography scoring band alignment working

opendm/config.py

@@ -157,15 +157,6 @@ def config(argv=None, parser=None):
                             'Can be one of: %(choices)s. Default: '
                             '%(default)s'))
     
-    parser.add_argument('--matcher-type',
-                        metavar='<string>',
-                        action=StoreValue,
-                        default='flann',
-                        choices=['flann', 'bow'],
-                        help=('Set feature matcher algorithm. FLANN is more robust, but slower. BOW is much faster, but can miss some valid matches. '
-                            'Can be one of: %(choices)s. Default: '
-                            '%(default)s'))
-
     parser.add_argument('--matcher-neighbors',
                         metavar='<integer>',
                         action=StoreValue,

opendm/multispectral.py

@@ -1,11 +1,16 @@
 import math
 import re
 import cv2
+import os
 from opendm import dls
 import numpy as np
 from opendm import log
 from opensfm.io import imread
 
+from skimage import exposure
+from skimage.morphology import disk
+from skimage.filters import rank, gaussian
+
 # Loosely based on https://github.com/micasense/imageprocessing/blob/master/micasense/utils.py
 
 def dn_to_radiance(photo, image):
@@ -181,9 +186,11 @@ def get_primary_band_name(multi_camera, user_band_name):
     return band_name_fallback
 
 
-def compute_primary_band_map(multi_camera, primary_band):
+def compute_band_maps(multi_camera, primary_band):
     """
-    Computes a map of { photo filename --> associated primary band photo }
+    Computes maps of: 
+     - { photo filename --> associated primary band photo } (s2p)
+     - { primary band filename --> list of associated secondary band photos } (p2s)
     by looking at capture time or filenames as a fallback
     """
     band_name = get_primary_band_name(multi_camera, primary_band)
@@ -196,7 +203,8 @@ def compute_primary_band_map(multi_camera, primary_band):
     # Try using capture time as the grouping factor
     try:
         capture_time_map = {}
-        result = {}
+        s2p = {}
+        p2s = {}
 
         for p in primary_band_photos:
             t = p.get_utc_time()
@@ -221,15 +229,19 @@ def compute_primary_band_map(multi_camera, primary_band):
                 if capture_time_map.get(t) is None:
                     raise Exception("Unreliable capture time detected (no primary band match)")
 
-                result[p.filename] = capture_time_map[t]
+                s2p[p.filename] = capture_time_map[t]
 
-        return result
+                if band['name'] != band_name:
+                    p2s.setdefault(capture_time_map[t].filename, []).append(p)
+
+        return s2p, p2s
     except Exception as e:
         # Fallback on filename conventions
         log.ODM_WARNING("%s, will use filenames instead" % str(e))
 
         filename_map = {}
-        result = {}
+        s2p = {}
+        p2s = {}
         file_regex = re.compile(r"^(.+)[-_]\w+(\.[A-Za-z]{3,4})$")
 
         for p in primary_band_photos:
@@ -251,15 +263,65 @@ def compute_primary_band_map(multi_camera, primary_band):
                 if filename_without_band == p.filename:
                     raise Exception("Cannot match bands by filename on %s, make sure to name your files [filename]_band[.ext] uniformly." % p.filename)
 
-                result[p.filename] = filename_map[filename_without_band]
-        
-        return result
+                s2p[p.filename] = filename_map[filename_without_band]
 
-def compute_aligned_image(image, align_image_filename, feature_retention=0.15):
-    if len(image.shape) != 3:
-        raise ValueError("Image should have shape length of 3 (got: %s)" % len(image.shape))
+                if band['name'] != band_name:
+                    p2s.setdefault(filename_map[filename_without_band].filename, []).append(p)
 
+        return s2p, p2s
+
+def compute_alignment_matrices(multi_camera, primary_band_name, images_path, s2p, p2s, max_samples=9999):
+    log.ODM_INFO("Computing band alignment")
+
+    alignment_info = {}
+
+    # For each secondary band
+    for band in multi_camera:
+        if band['name'] != primary_band_name:
+            matrices = []
+
+            # if band['name'] != "NIR":
+            #     continue # TODO REMOVE
+
+            # 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
+                
+                warp_matrix, score, dimension = compute_homography(os.path.join(images_path, p.filename),
+                                                            os.path.join(images_path, primary_band_photo.filename))
+
+                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))
+                    
+                if len(matrices) >= max_samples:
+                    log.ODM_INFO("Got enough samples for %s (%s)" % (band['name'], max_samples))
+                    break
+            
+            # Sort
+            matrices.sort(key=lambda x: x['score'], reverse=False)
+            
+            if len(matrices) > 0:
+                alignment_info[band['name']] = matrices[0]
+                print(matrices[0])
+            else:
+                log.ODM_WARNING("Cannot find alignment matrix for band %s, The band will likely be misaligned!" % band['name'])
+
+    return alignment_info, p2s
+
+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:
@@ -270,7 +332,60 @@ def compute_aligned_image(image, align_image_filename, feature_retention=0.15):
         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)
+
+        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, 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, None)
+        
+        ratio = svd[0] / svd[-1]
+        if ratio > 100000:
+            return None, None, (None, None)
+
+        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)
+
+def find_ecc_homography(image_gray, align_image_gray, number_of_iterations=5000, termination_eps=1e-8):
+    image_gray = to_8bit(gradient(gaussian(image_gray)))
+    align_image_gray = to_8bit(gradient(gaussian(align_image_gray)))
+
+    # Define the motion model
+    warp_matrix = np.eye(3, 3, dtype=np.float32)
+
+    # Define termination criteria
+    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)
+
+    return warp_matrix
+    
+
+def find_features_homography(image_gray, align_image_gray, feature_retention=0.25):
     # Detect SIFT features and compute descriptors.
     detector = cv2.SIFT_create(edgeThreshold=10, contrastThreshold=0.1)
     kp_image, desc_image = detector.detectAndCompute(image_gray, None)
@@ -301,12 +416,62 @@ def compute_aligned_image(image, align_image_filename, feature_retention=0.15):
 
     # Find homography
     h, _ = cv2.findHomography(points_image, points_align_image, cv2.RANSAC)
+    return h
 
-    # Use homography
-    height, width = align_image_gray.shape
-    aligned_image = cv2.warpPerspective(image, h, (width, height))
+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)
+    grad_y = cv2.Sobel(im,cv2.CV_32F,0,1,ksize=ksize)
+    grad = cv2.addWeighted(np.absolute(grad_x), 0.5, np.absolute(grad_y), 0.5, 0)
+    return grad
+
+def local_normalize(im):
+    width, _ = im.shape
+    disksize = int(width/5)
+    if disksize % 2 == 0:
+        disksize = disksize + 1
+    selem = disk(disksize)
+    im = rank.equalize(im, selem=selem)
+    return im
+
+
+def align_image(image, warp_matrix, dimension):
+    if warp_matrix.shape == (3, 3):
+        return cv2.warpPerspective(image, warp_matrix, dimension)
+    else:
+        return cv2.warpAffine(image, warp_matrix, dimension)
 
-    return aligned_image
 
 def to_8bit(image):
     if image.dtype == np.uint8:
@@ -323,6 +488,10 @@ def to_8bit(image):
     image = image.astype(np.float32)
     image *= 255.0 / value_range
     np.around(image, out=image)
+    image[image > 255] = 255
+    image[image < 0] = 0
     image = image.astype(np.uint8)
 
-    return image
+    return image
+
+

opendm/osfm.py

@@ -105,7 +105,7 @@ class OSFMContext:
                 except Exception as e:
                     log.ODM_WARNING("Cannot set camera_models_overrides.json: %s" % str(e))
 
-            use_bow = args.matcher_type == "bow"
+            use_bow = bool(reconstruction.multi_camera)
             feature_type = "SIFT"
 
             # GPSDOP override if we have GPS accuracy information (such as RTK)
@@ -315,16 +315,67 @@ class OSFMContext:
         else:
             log.ODM_INFO("Already extracted cameras")
     
-    def convert_and_undistort(self, rerun=False, imageFilter=None):
+    def convert_and_undistort(self, rerun=False, imageFilter=None, image_list=None):
         log.ODM_INFO("Undistorting %s ..." % self.opensfm_project_path)
         undistorted_images_path = self.path("undistorted", "images")
 
         if not io.dir_exists(undistorted_images_path) or rerun:
-            undistort.run_dataset(DataSet(self.opensfm_project_path), "reconstruction.json", 
+            ds = DataSet(self.opensfm_project_path)
+
+            if image_list is not None:
+                ds._set_image_list(image_list)
+
+            undistort.run_dataset(ds, "reconstruction.json", 
                                   0, None, "undistorted", imageFilter)
         else:
             log.ODM_WARNING("Found an undistorted directory in %s" % undistorted_images_path)
 
+    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:
+            reconstruction = json.loads(f.read())
+
+        # Augment reconstruction.json
+        for recon in reconstruction:
+            shots = recon['shots']
+            sids = list(shots)
+            
+            for shot_id in sids:
+                secondary_photos = p2s.get(shot_id)
+                if secondary_photos is None:
+                    log.ODM_WARNING("Cannot find secondary photos for %s" % shot_id)
+                    continue
+
+                for p in secondary_photos:
+                    shots[p.filename] = shots[shot_id]
+
+        with open(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()

stages/mvstex.py

@@ -5,6 +5,7 @@ from opendm import io
 from opendm import system
 from opendm import context
 from opendm import types
+from opendm.multispectral import get_primary_band_name
 
 class ODMMvsTexStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -36,8 +37,9 @@ class ODMMvsTexStage(types.ODM_Stage):
                 }]
 
         if reconstruction.multi_camera:
+
             for band in reconstruction.multi_camera:
-                primary = band == reconstruction.multi_camera[0]
+                primary = band['name'] == get_primary_band_name(reconstruction.multi_camera, args.primary_band)
                 nvm_file = os.path.join(tree.opensfm, "undistorted", "reconstruction_%s.nvm" % band['name'].lower())
                 add_run(nvm_file, primary, band['name'].lower())
         else:

stages/odm_georeferencing.py

@@ -9,6 +9,7 @@ from opendm import system
 from opendm import context
 from opendm.cropper import Cropper
 from opendm import point_cloud
+from opendm.multispectral import get_primary_band_name
 
 class ODMGeoreferencingStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -45,7 +46,7 @@ class ODMGeoreferencingStage(types.ODM_Stage):
         
         if reconstruction.multi_camera:
             for band in reconstruction.multi_camera:
-                primary = band == reconstruction.multi_camera[0]
+                primary = band['name'] == get_primary_band_name(reconstruction.multi_camera, args.primary_band)
                 add_run(primary, band['name'].lower())
         else:
             add_run()

stages/odm_orthophoto.py

@@ -11,6 +11,7 @@ from opendm.concurrency import get_max_memory
 from opendm.cutline import compute_cutline
 from pipes import quote
 from opendm import pseudogeo
+from opendm.multispectral import get_primary_band_name
 
 class ODMOrthoPhotoStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -72,7 +73,7 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
 
             if reconstruction.multi_camera:
                 for band in reconstruction.multi_camera:
-                    primary = band == reconstruction.multi_camera[0]
+                    primary = band['name'] == get_primary_band_name(reconstruction.multi_camera, args.primary_band)
                     subdir = ""
                     if not primary:
                         subdir = band['name'].lower()

stages/openmvs.py

@@ -8,6 +8,7 @@ from opendm import point_cloud
 from opendm import types
 from opendm.utils import get_depthmap_resolution
 from opendm.osfm import OSFMContext
+from opendm.multispectral import get_primary_band_name
 
 class ODMOpenMVSStage(types.ODM_Stage):
     def process(self, args, outputs):
@@ -30,8 +31,9 @@ class ODMOpenMVSStage(types.ODM_Stage):
             cmd = 'export_openmvs'
             if reconstruction.multi_camera:
                 # Export only the primary band
-                primary = reconstruction.multi_camera[0]
-                image_list = os.path.join(tree.opensfm, "image_list_%s.txt" % primary['name'].lower())
+                primary_band = get_primary_band_name(reconstruction.multi_camera, args.primary_band)
+
+                image_list = os.path.join(tree.opensfm, "image_list_%s.txt" % primary_band.lower())
                 cmd += ' --image_list "%s"' % image_list
             octx.run(cmd)
             

stages/run_opensfm.py

@@ -63,7 +63,8 @@ class ODMOpenSfMStage(types.ODM_Stage):
 
         # Undistorted images will be used for texturing / MVS
 
-        primary_band_map = None
+        alignment_info = None
+        primary_band_name = None
         undistort_pipeline = []
 
         def undistort_callback(shot_id, image):
@@ -77,38 +78,52 @@ class ODMOpenSfMStage(types.ODM_Stage):
 
 
         def align_to_primary_band(shot_id, image):
-            primary_band_photo = primary_band_map.get(shot_id)
+            photo = reconstruction.get_photo(shot_id)
+
+            # No need to align primary
+            if photo.band_name == primary_band_name:
+                return image
+
+            ainfo = alignment_info.get(photo.band_name)
+            if ainfo is not None:
+                return multispectral.align_image(image, ainfo['warp_matrix'], ainfo['dimension'])
+            else:
+                log.ODM_WARNING("Cannot align %s, no alignment matrix could be computed. Band alignment quality might be affected." % (shot_id))
+                return image
 
-            try:
-                if primary_band_photo is None:
-                    raise Exception("Cannot find primary band file for %s" % shot_id)
-                
-                if primary_band_photo.filename == shot_id:
-                    # This is a photo from the primary band, skip
-                    return image
-                
-                align_image_filename = os.path.join(tree.dataset_raw, primary_band_map[shot_id].filename)
-                image = multispectral.compute_aligned_image(image, align_image_filename)
-            except Exception as e:
-                log.ODM_WARNING("Cannot align %s: %s. Output quality might be affected." % (shot_id, str(e)))
-            
-            return image
 
         if args.radiometric_calibration != "none":
             undistort_pipeline.append(radiometric_calibrate)
         
-        if reconstruction.multi_camera:
-            primary_band_map = multispectral.compute_primary_band_map(reconstruction.multi_camera, args.primary_band)
-            
-            # TODO:
-            # if (not flag_file exists or rerun)
-            # Change reconstruction.json and tracks.csv by adding copies of the primary
-            # band camera to the reconstruction
+        image_list_override = None
 
+        if reconstruction.multi_camera:
             
+            # Undistort only secondary bands
+            image_list_override = [os.path.join(tree.dataset_raw, p.filename) for p in photos] # if p.band_name.lower() != primary_band_name.lower()
+
+            # We backup the original reconstruction.json, tracks.csv
+            # then we augment them by duplicating the primary band
+            # camera shots with each band, so that exports, undistortion,
+            # etc. include all bands
+            # 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)
+                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)
 
-        octx.convert_and_undistort(self.rerun(), undistort_callback)
+        octx.convert_and_undistort(self.rerun(), undistort_callback, image_list_override)
 
         self.update_progress(80)
 
@@ -126,9 +141,6 @@ class ODMOpenSfMStage(types.ODM_Stage):
                 else:
                     log.ODM_WARNING("Found a valid image list in %s for %s band" % (image_list_file, band['name']))
                 
-                # TODO!! CHANGE
-                print("TODO!")
-                exit(1)
                 nvm_file = octx.path("undistorted", "reconstruction_%s.nvm" % band['name'].lower())
                 if not io.file_exists(nvm_file) or self.rerun():
                     octx.run('export_visualsfm --points --image_list "%s"' % image_list_file)
@@ -136,6 +148,17 @@ 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)
+
         if not io.file_exists(tree.opensfm_reconstruction_nvm) or self.rerun():
             octx.run('export_visualsfm --points')
         else: