Cross modal satellite coregistration using cv2.optflow

I am experimenting with the library cv2.optflow to manage coregistration between two satellite images (optical vs. SAR). To date, I have tried calcOpticalFlowFarneback with modest results. Do you have any suggestions of more advanced functions available in outflow that fit my coregistration scenario? Thanks in advance.

last I checked, “DIS” optical flow was the most advanced in the module.

optical flow is probably too general for what you need.

look for findTransformECC instead. the sigma parameter affects how much initial misalignment the ECC algorithm can tolerate, but also what detail it can see.

if the misalignment is large enough, you’ll have to feature-match first.

for the “cross modal” situation, just run a laplacian on both pics before calculating alignment. the gradient maps probably look more similar than the pictures themselves.

Thanks for your contribution.
Just a specific clarification: sigma here refers to a Gaussian blur applied before alignment ? I believe it is not directly related to findTransformECC.

no, I mean gaussFiltSize of the function itself

when working with gradient maps (laplacian), expect ECC to have trouble converging. you’ll need several passes, starting from large sigmas and going down (recommend geometric progression/powers of two). the thing is basically a line drawing, and the difference of two such only gives worse gradient data (another gradient) for minimization.

you might want to provide a pair of pictures (both optical and SAR) for experimentation.

I am afraid that gaussFiltSize is not supported in OpenCV 4.11.0 (my current version)

how did you determine this?

the docs say it’s there. it’s always been there.

From what I can see, that parameter is not available with Python bindings. Anyway, I have some trouble with convergence in a single test case, where I get The algorithm stopped before its convergence. The correlation is going to be minimized. Images may be uncorrelated or non-overlapped in function 'findTransformECC'

I don’t believe that. it’s there. show me what you’re doing.

“test case”? MRE please. data and code.

Yes, I understand, but cannot post data depending on a signed NDA.
This is what I am doing with findTransformECC

def run_ecc_with_initial_transform(optical, sar, warp_matrix, number_of_levels=3, warp_mode=cv2.MOTION_AFFINE):
    warp_matrix = sanitize_affine_matrix(warp_matrix, max_translation=500, max_scale=1.5)  # Lower thresholds

    if not is_valid_affine(warp_matrix):
        print("[!] Feature matrix suspicious, using identity for ECC.")
        warp_matrix = np.eye(2, 3, dtype=np.float32)

    criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 500, 1e-6)

    for level in reversed(range(number_of_levels)):
        scale = 1 / (2 ** level)
        optical_resized = cv2.resize(optical, (0, 0), fx=scale, fy=scale, interpolation=cv2.INTER_LINEAR)
        sar_resized = cv2.resize(sar, (0, 0), fx=scale, fy=scale, interpolation=cv2.INTER_LINEAR)
        print(f"[~] ECC level {level} input warp matrix:\n{warp_matrix}")

        # Reset warp_matrix at the highest resolution (level 0) if necessary
        if level == number_of_levels - 1:
            warp_matrix = np.eye(2, 3, dtype=np.float32) 

        try:
            warp_matrix = warp_matrix.astype(np.float32)
            cc, warp_matrix = cv2.findTransformECC(
                templateImage=optical_resized,
                inputImage=sar_resized,
                warpMatrix=warp_matrix,
                motionType=warp_mode,
                criteria=criteria
            )
            print(f"[✓] ECC converged at level {level} (scale={scale:.2f}) with cc={cc:.5f}")

            # Additional diagnostics for the warp matrix
            print(f"[~] ECC level {level} warp matrix after transformation:\n{warp_matrix}")
            print(f"[~] Warp matrix translation at level {level}: {warp_matrix[0, 2]}, {warp_matrix[1, 2]}")
            print(f"[~] Warp matrix scale at level {level}: {np.linalg.norm(warp_matrix[:2, :2])}")
            
        except cv2.error as e:
            print(f"[!] ECC failed at level {level}: {e}")
            continue

        # Scale up warp matrix for next level
        if level > 0:
            warp_matrix[0, 2] *= 2
            warp_matrix[1, 2] *= 2

    return warp_matrix

As said the parameter gaussFiltSize is not recognized as valid, even if I can see it in C version.
Later I will add a working snippet, just to make easier and testable a complete flow.

Thanks

This is a working script that I use for testing:

import cv2
import numpy as np
import rasterio

def read_image_grayscale(path):
    with rasterio.open(path) as src:
        image = src.read(1).astype(np.float32)
        image = np.log1p(image) if 'sar' in path.lower() else image  # log-enhance SAR
        image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX)
        return image.astype(np.uint8)

def preprocess(image):
    clahe = cv2.createCLAHE(clipLimit=3.0, tileGridSize=(8,8))
    return clahe.apply(image)

def resize_to_match(base, target):
    return cv2.resize(target, (base.shape[1], base.shape[0]))

def build_pyramid(image, levels):
    pyramid = [image]
    for _ in range(1, levels):
        image = cv2.pyrDown(image)
        pyramid.append(image)
    return pyramid[::-1]  # coarse to fine

def run_multiscale_ecc(optical, sar, levels=3, warp_mode=cv2.MOTION_TRANSLATION):
    # Ensure same size
    sar = resize_to_match(optical, sar)

    # Preprocess both
    optical = preprocess(optical)
    sar = preprocess(sar)

    # Build pyramids
    optical_pyr = build_pyramid(optical, levels)
    sar_pyr = build_pyramid(sar, levels)

    # Initialize warp matrix
    warp_matrix = np.eye(2, 3, dtype=np.float32)

    for i in range(levels):
        scale = 2 ** (levels - i - 1)
        optical_level = optical_pyr[i]
        sar_level = sar_pyr[i]

        try:
            cc, warp_matrix = cv2.findTransformECC(
                templateImage=optical_level,
                inputImage=sar_level,
                warpMatrix=warp_matrix,
                motionType=warp_mode,
                criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 200, 1e-6),
                inputMask=None,
                gaussFiltSize=5  # Gaussian smoothing for robustness
            )
            print(f"[✓] Level {i+1}/{levels}, Scale 1/{scale}, Corr={cc:.6f}")
        except cv2.error as e:
            print(f"[✗] ECC failed at level {i+1}: {e}")
            break

        # Scale warp for next level if needed
        if i < levels - 1:
            warp_matrix[:2, 2] *= 2.0

    print("Final warp matrix:\n", warp_matrix)
    return warp_matrix

if __name__ == "__main__":
    import sys
    if len(sys.argv) != 3:
        print("Usage: python ecc_pyramid.py <optical.tif> <sar.tif>")
        sys.exit(1)
optical_path = sys.argv[1]
sar_path = sys.argv[2]

optical = read_image_grayscale(optical_path)
sar = read_image_grayscale(sar_path)

run_multiscale_ecc(optical, sar, levels=3, warp_mode=cv2.MOTION_TRANSLATION)

I am testing two pair of images, on the first one I get decent results:

[✓] Level 1/3, Scale 1/4, Corr=0.695040
[✓] Level 2/3, Scale 1/2, Corr=0.647860
[✓] Level 3/3, Scale 1/1, Corr=0.600530
Final warp matrix:
 [[1.         0.         0.07731839]
 [0.         1.         0.4849937 ]]

but on the second, even I tried several tricks, I always get this:

[✗] ECC failed at level 1: OpenCV(4.11.0) /io/opencv/modules/video/src/ecc.cpp:589: error: (-7:Iterations do not converge) The algorithm stopped before its convergence. The correlation is going to be minimized. Images may be uncorrelated or non-overlapped in function 'findTransformECC'

Final warp matrix:
 [[1. 0. 0.]
 [0. 1. 0.]]

Thank you again.

p.s. you were right regarding gaussFiltSize.