Adding alternate matching methods (#11)

- Added alternate matching methods to s2p (LoFTR, SuperGlue and SIFT rectification is now more stable!)

---------

Co-authored-by: twonzo72@gmail.com <timiles@tims-air-5.mynet>

Author: tonzowonzo

s2p/__init__.py

@@ -127,7 +127,8 @@ def pointing_correction(tile, i):
     method = 'relative' if cfg['relative_sift_match_thresh'] is True else 'absolute'
     A, m = pointing_accuracy.compute_correction(
         img1, img2, rpc1, rpc2, x, y, w, h, method,
-        cfg['sift_match_thresh'], cfg['max_pointing_error']
+        cfg['sift_match_thresh'], cfg['max_pointing_error'], cfg['matching_method'],
+        cfg['min_value'], cfg['max_value'], cfg['confidence_threshold']
     )
 
     if A is not None:  # A is the correction matrix
@@ -645,17 +646,18 @@ def main(user_cfg, start_from=0, merge_matches=False):
         print('2) correcting pointing globally...')
         global_pointing_correction(tiles)
         common.print_elapsed_time()
-    
+        
     # Create matches GeoJSON.
-    print("Creating matches GeoJSON")
     if merge_matches:
+        print("Creating matches GeoJSON")
         merge_all_match_files()
         matches_to_geojson(f"{cfg['out_dir']}/merged_sift_matches.txt",
                         cfg['images'][0]['rpcm'],
                         10,
                         [0, 1],
                         f"{cfg['out_dir']}/matches.geojson"
         )
+
     # rectification step:
     if start_from <= 3:
         print('3) rectifying tiles...')

s2p/config.py

@@ -177,3 +177,13 @@
 # If the out_crs is not set this parameter determines if the output CRS uses the EGM2008 geoid vertical datum (if True)
 # or the WGS84 ellipsoid vertical datum (if False). If out_crs is set, this parameter is ignored.
 cfg['out_geoid'] = False
+
+# LoFTR parameters, points whose confidence is below this value are ignored.
+cfg["confidence_threshold"] = 0.5
+
+# Matching parameters
+# Matching method to use, choose from ["sift", "loftr", "superglue", "all"]
+cfg["matching_method"] = "sift"
+# The min and max valeus to scale the image between 0 and 1 with, only used for non-SIFT methods.
+cfg["max_value"] = 3000
+cfg["min_value"] = 200

s2p/image_coordinates_to_coordinates.py

@@ -1,10 +1,16 @@
+import glob
+import os
+
+import rasterio
 import numpy as np
 import geopandas as gpd
 import shapely
-import srtm4
-from rpcm import RPCModel, rpc_from_rpc_file
+from shapely.geometry import Point
 import geopy
 import geopy.distance
+import srtm4
+
+from rpcm import RPCModel, rpc_from_rpc_file
 
 
 def localize_row_col_geometry(

s2p/pointing_accuracy.py

@@ -100,7 +100,8 @@ def local_translation(r1, r2, x, y, w, h, m):
 
 
 def compute_correction(img1, img2, rpc1, rpc2, x, y, w, h,
-                       method, sift_thresh, epipolar_threshold):
+                       method, sift_thresh, epipolar_threshold,
+                       matching_method, min_value, max_value, confidence_threshold):
     """
     Computes pointing correction matrix for specific ROI
 
@@ -125,7 +126,8 @@ def compute_correction(img1, img2, rpc1, rpc2, x, y, w, h,
         raise ValueError(f"width or height <= 0 for:\n{img1}\n{img2}\nx={x}, y={y}, w={w}, h={h}. Try a different"
                          f"tilesize or different ROI.")
     m = sift.matches_on_rpc_roi(img1, img2, rpc1, rpc2, x, y, w, h,
-                                method, sift_thresh, epipolar_threshold)
+                                method, sift_thresh, epipolar_threshold,
+                                matching_method, min_value, max_value, confidence_threshold)
 
     if m is not None:
         A = local_translation(rpc1, rpc2, x, y, w, h, m)

s2p/rectification.py

@@ -320,11 +320,19 @@ def rectify_pair(im1, im2, rpc1, rpc2, x, y, w, h, out1, out2, A=None, sift_matc
 
     else:
         raise Exception("Unknown value {} for argument 'method'".format(method))
+        
+       
+            
+    # If there are still no matches, raise the no matches error.
+    if matches is None or len(matches) < 4 or len(matches.shape) != 2:
+        # find virtual matches from RPC camera models
+        matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h,
+                                             cfg['n_gcp_per_axis'])
 
-    if matches is None or len(matches) < 4:
-        raise NoRectificationMatchesError(
-            "No or not enough matches found to rectify image pair"
-        )
+        # correct second image coordinates with the pointing correction matrix
+        if A is not None:
+            matches[:, 2:] = common.points_apply_homography(np.linalg.inv(A),
+                                                            matches[:, 2:])
 
     # compute rectifying homographies
     H1, H2, F = rectification_homographies(matches, x, y, w, h)

s2p/sift.py

@@ -10,11 +10,25 @@
 import numpy as np
 import rasterio as rio
 from numpy.ctypeslib import ndpointer
+import cv2
 import ransac
 
 from s2p import rpc_utils
 from s2p import estimation
 
+try:
+    from matching.matcher import get_keypoints_loftr
+except ImportError as e:
+    print("LoFTR matching library could not be imported.")
+    print(e)
+
+try:
+    from matching.matcher import get_keypoints_superglue
+except ImportError as e:
+    print("SuperGlue matching library could not be imported.")
+    print(e)
+
+
 # Locate sift4ctypes library and raise an ImportError if it can not be
 # found This call will raise an exception if library can not be found,
 # at import time
@@ -240,13 +254,12 @@ def keypoints_match_from_nparray(k1, k2, method, sift_threshold,
 
 
 def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h,
-                       method, sift_thresh, epipolar_threshold):
+                       method, sift_thresh, epipolar_threshold, matching_method="sift",
+                       min_value=200, max_value=3000, confidence_threshold=0.5):
     """
     Compute a list of SIFT matches between two images on a given roi.
-
     The corresponding roi in the second image is determined using the rpc
     functions.
-
     Args:
         im1, im2: paths to two large tif images
         rpc1, rpc2: two instances of the rpcm.RPCModel class
@@ -255,7 +268,6 @@ def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h,
             of the rectangle.
         method, sift_thresh, epipolar_threshold: see docstring of
             s2p.sift.keypoints_match()
-
     Returns:
         matches: 2D numpy array containing a list of matches. Each line
             contains one pair of points, ordered as x1 y1 x2 y2.
@@ -267,18 +279,100 @@ def matches_on_rpc_roi(im1, im2, rpc1, rpc2, x, y, w, h,
     rpc_matches = rpc_utils.matches_from_rpc(rpc1, rpc2, x, y, w, h, 5)
     F = estimation.affine_fundamental_matrix(rpc_matches)
 
-    # if less than 10 matches, lower thresh_dog. An alternative would be ASIFT
-    thresh_dog = 0.0133
-    for _ in range(2):
-        p1 = image_keypoints(im1, x, y, w, h, thresh_dog=thresh_dog)
-        p2 = image_keypoints(im2, x2, y2, w2, h2, thresh_dog=thresh_dog)
-        matches = keypoints_match(p1, p2, method, sift_thresh, F,
-                                  epipolar_threshold=epipolar_threshold,
-                                  model='fundamental')
-        if matches is not None and matches.ndim == 2 and matches.shape[0] > 10:
-            break
-        thresh_dog /= 2.0
-    else:
-        print("WARNING: sift.matches_on_rpc_roi: found no matches.")
-        return None
+    if matching_method == "loftr":
+        p1, p2 = get_keypoints_loftr(im1, im2, min_value, max_value, confidence_threshold, 
+                                     x, x2, y, y2, w, w2, h, h2, rpc_match=True)
+        if len(p1) < 10 or len(p2) < 10:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            return None
+        p1[:, 0] += x
+        p1[:, 1] += y
+        p2[:, 0] += x2
+        p2[:, 1] += y2
+        matches = np.hstack((p1, p2))
+        
+        if len(matches) == 0 or matches.ndim != 2:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            return None
+        elif len(matches) > 10:
+            inliers = ransac.find_fundamental_matrix(matches, ntrials=1000,
+                                                 max_err=0.3)[0]
+            matches = matches[inliers]
+        else:
+            print("WARNING: sift.matches_on_rpc_roi: not enough matches.")
+            return None
+        
+        
+    elif matching_method == "superglue":
+        p1, p2 = get_keypoints_superglue(im1, im2, min_value, max_value, x, x2, y, y2, w, w2, h, h2, rpc_match=True)
+        if len(p1) < 10 or len(p2) < 10:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            return None
+        
+        p1[:, 0] += x
+        p1[:, 1] += y
+        p2[:, 0] += x2
+        p2[:, 1] += y2
+        matches = np.hstack((p1, p2))
+        
+        if matches.ndim != 2:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            return None
+
+    elif matching_method == "sift":
+        # if less than 10 matches, lower thresh_dog. An alternative would be ASIFT
+        thresh_dog = 0.0133
+        for _ in range(2):
+            p1 = image_keypoints(im1, x, y, w, h, thresh_dog=thresh_dog)
+            p2 = image_keypoints(im2, x2, y2, w2, h2, thresh_dog=thresh_dog)
+            matches = keypoints_match(p1, p2, method, sift_thresh, F,
+                                    epipolar_threshold=epipolar_threshold,
+                                    model='fundamental')
+            if matches is not None and matches.ndim == 2 and matches.shape[0] > 10:
+                break
+            thresh_dog /= 2.0
+        else:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            return None
+
+        
+    elif matching_method == "all":
+        p1, p2 = get_keypoints_loftr(im1, im2, min_value, max_value, confidence_threshold, 
+                                     x, x2, y, y2, w, w2, h, h2, rpc_match=True)
+        p1[:, 0] += x
+        p1[:, 1] += y
+        p2[:, 0] += x2
+        p2[:, 1] += y2
+        matches = np.hstack((p1, p2))
+        if len(matches) == 0 or matches.ndim != 2:
+            print("WARNING: sift.matches_on_rpc_roi: found no matches.")
+            matches = None
+        elif len(matches) > 10:
+            inliers = ransac.find_fundamental_matrix(matches, ntrials=1000,
+                                                 max_err=0.3)[0]
+            matches = matches[inliers]
+        
+        p1_sg, p2_sg = get_keypoints_superglue(im1, im2, min_value, max_value, x, x2, y, y2, w, w2, h, h2, rpc_match=True)
+        p1_sg[:, 0] += x
+        p1_sg[:, 1] += y
+        p1_sg[:, 0] += x2
+        p1_sg[:, 1] += y2
+        matches_sg = np.hstack((p1_sg, p2_sg))
+
+        thresh_dog = 0.0133
+        for _ in range(2):
+            p1_sift = image_keypoints(im1, x, y, w, h, thresh_dog=thresh_dog)
+            p2_sift = image_keypoints(im2, x2, y2, w2, h2, thresh_dog=thresh_dog)
+            matches_sift = keypoints_match(p1_sift, p2_sift, method, sift_thresh, F,
+                                    epipolar_threshold=epipolar_threshold,
+                                    model='fundamental')
+            if matches_sift is not None and matches_sift.ndim == 2 and matches_sift.shape[0] > 10:
+                break
+            thresh_dog /= 2.0
+        
+        matches = np.vstack((matches, matches_sg, matches_sift))
+        
+        if len(matches) < 10:
+            return None
+
     return matches

setup.py

@@ -54,6 +54,7 @@ def finalize_options(self):
                 'rpcm @ https://github.com/20treeAI/rpcm/archive/refs/tags/v1.4.8.tar.gz',
                 'srtm4 @ https://github.com/20treeAI/srtm4/archive/refs/tags/1.2.4.tar.gz',
                 'requests',
+                'opencv-python',
                 'geopandas',
                 'geopy']
 
@@ -62,7 +63,7 @@ def finalize_options(self):
 }
 
 setup(name="s2p",
-      version="1.3.9",
+      version="1.4.0",
       description="Satellite Stereo Pipeline.",
       long_description=readme(),
       long_description_content_type='text/markdown',

tests/rectification_test.py

@@ -39,50 +39,6 @@ def test_rectification_homographies(matches):
         np.testing.assert_allclose(variable, expected, rtol=0.01, atol=1e-6, verbose=True)
 
 
-def test_rectify_pair_no_matches(tmp_path, images):
-    """
-    Test running rectification.rectify_pair() where no matches are found.
-    """
-    im1, rpc1, im2, rpc2 = images
-    with pytest.raises(s2p.rectification.NoRectificationMatchesError):
-        s2p.rectification.rectify_pair(
-            im1=im1,
-            im2=im2,
-            rpc1=rpc1,
-            rpc2=rpc2,
-            x=100,
-            y=100,
-            w=200,
-            h=200,
-            out1=str(tmp_path / 'out1.tiff'),
-            out2=str(tmp_path / 'out2.tiff'),
-            sift_matches=None,
-            method='sift',
-        )
-
-
-def test_rectify_pair_few_matches(tmp_path, matches, images):
-    """
-    Test running rectification.rectify_pair() where less than 4 matches are found.
-    """
-    im1, rpc1, im2, rpc2 = images
-    with pytest.raises(s2p.rectification.NoRectificationMatchesError):
-        s2p.rectification.rectify_pair(
-            im1=im1,
-            im2=im2,
-            rpc1=rpc1,
-            rpc2=rpc2,
-            x=100,
-            y=100,
-            w=200,
-            h=200,
-            out1=str(tmp_path / 'out1.tiff'),
-            out2=str(tmp_path / 'out2.tiff'),
-            sift_matches=matches[:3],
-            method='sift',
-        )
-
-
 def test_rectify_pair_with_matches(tmp_path, matches, images):
     """
     Test running rectification.rectify_pair() with some matches.