📖 Overview

This repository contains a Python-based pipeline designed to align (coregister) multiple satellite images with high precision. Unlike traditional global homography approaches, this tool uses a hierarchical, grid-based strategy to account for local non-linear distortions and topographic parallax. It employs a combination of Semantic Segmentation, DIScrete Keypoints (DISK) feature extraction and LightGlue feature matching to optimize resources and achieve state-of-the-art performance.

Key Features

1. Dual-Purpose Semantic Masking: Uses Deep Learning to both localize the Region of Interest and mask out "hyper-stable" elevated terrain that causes parallax errors.

2. Hierarchical Pipeline: Configurable multi-stage alignment (e.g. Coarse Grid $\to$ Fine Grid $\to$ Final Warp).

3. Robust Matching: Uses GPU-enabled implemenatations of DISK and LightGlue for quick scale-agnostic feature matching.

4. Warping: Achieves pixel-perfect alignment in 3D imaging geometry by utilizing the Thin Plate Spline (TPS) transformation.

5. Geospatial Support: Built on top of GDAL to handle georeferenced JP2 and GeoTIFF inputs.

🧠 Method Details

Figure 1: Visualization of a displacement vector field used in image co-registration. Red dots and green arrows denote local tie points and their corresponding shift trajectories. Anomalous mappings have been removed using a Vector Median Filter (VMF), yielding a coherent vector field to be used as input for a final Thin Plate Spline (TPS) warping transform.

Key Algorithms

1. Initial Feature Matching with DISK & LightGlue: We utilize DISK (DIScrete Keypoints) to extract a dense and highly discriminative set of local features from both the target and reference images. This reinforcement learning-driven extractor ensures robust keypoint detection even across challenging terrain textures or varying viewpoints. To pair these features, we employ LightGlue, a state-of-the-art deep neural network matcher that dynamically prunes unmatchable points to save memory and computation. Together, this combination establishes a robust, high-quality initial set of tie points required for the downstream geometric filtering stages.

2. Vector Median Filter (VMF): To refine the initial keypoints, we process the matched coordinates through a two-stage Vector Median Filter that evaluates the geometric consistency of the displacement vectors. First, a global hard threshold eliminates obvious outliers by discarding any matches where the translation magnitude exceeds a predefined maximum shift. Next, a spatial KD-Tree identifies the nearest neighbors for each remaining point to calculate the localized median flow vector. The algorithm compares each point's actual vector against its local median using a dynamic tolerance threshold—allowing basic wiggle room for minor camera jitter while scaling proportionally to accommodate larger regional movements. Any points deviating beyond this adaptive boundary are removed, leaving only a highly coherent vector field free of anomalous shifts.

3. Thin Plate Spline (TPS) warping: To finalize the co-registration, we apply a Thin Plate Spline (TPS) transformation to warp the target image into precise alignment with the reference image based on our filtered vector field. Unlike simple affine or homography transforms that assume a flat, two-dimensional scene, TPS computes a non-linear, locally adaptive spatial deformation modeled after the physical bending of a thin metal sheet. This flexibility makes it particularly powerful for handling imagery captured with complex 3D camera movements over uneven terrain, which inherently introduces severe parallax distortions.

Hierarchical Matching

Figure 2: The ROI (left) is subdivided into individual tiles (right) to facilitate local alignment. This structure supports a multi-pass approach with shrinking grid sizes, where local offsets derived from coarser grids serve as reference points for finer tiles in subsequent iterations.

The pipeline fits a grid over the identified region of interest (ROI) by following these steps:

• (1/5) Crop: Extract corresponding patches from Reference and Source images. On the first iteration, crops are identical; otherwise, query the KDTree with the patch center to retrieve a local offset from the Translation Matrix, shifting the Source crop to a pre-aligned starting position.
• (2/6) Match & filter: Compute matches with DISK & LightGlue, filter them with a custom Vector Median Filter.
• (3) Solve (translation): Calculate the local offsets.
• (4) Interpolate: Translation for hard to allign tiles are smoothly interpolated based on their neighbours.
• (7) Solve (affine): Calculate the local warp using Thin Plate Spline transform (TPS).
• (8) Warp & save: Apply the final transform to cropped tiles and save the resulting images.

   YES ------------ Repeat with finer grid?
    |                                     |
    V                                     |
1 -----+   2 ------+   3 ------+   4 ------------+   5 -----+   6 ------+   7 ------+   8 -----+
| Crop | > | Match | > | Solve | > | Interpolate | > | Crop | > | Match | > | Solve | > | Warp |
+------+   +-------+   +-------+   +-------------+   +------+   +-------+   +-------+   +------+

📂 Repository structure

config/
│   ├─__init__.py --------------- # Parse the config.yaml
│   ├─cli.py -------------------- # Command Line Interface
│   ├─config.py ----------------- # BaseConfig() to store the defaults
│   └─config.yaml --------------- # Default parameters
|
core/
│   ├─convert_to_tif.py ---------- # Paralellized .JP2 to .tif conversion 
│   ├─coreg_tools.py ------------- # Main coregistration algorithm
│   ├─gdal_utils.py -------------- # GDAL integration 
│   └─uitls.py ------------------- # Misc helper functions
|
models/
│   ├─dataloader.py -------------- # Inference dataloader 
│   ├─MobaNet.py ----------------- # PyTorch model architechture 
│   └─segmnet.py ----------------- # Main segmentation routine
│
weights/
│   └─model.pth ------------------ # Pre-trained MobaNet weights
|
coregister.py -------------------- # Pipeline entry-point
README.md

🚀 Getting Started

Prerequisites

This project relies heavily on GDAL and PyTorch. It is recommended to use conda for environment management to handle GDAL dependencies gracefully.

# 1. Create environment
(base) x@y:~$ conda create -n corg python=3.10
(base) x@y:~$ conda activate corg
(corg) x@y:~$ 

# 2. Install GDAL (3.8.5 is stable)
(corg) x@y:~$ conda install -c conda-forge gdal=3.8.5

# 3. Install PyTorch (Adjust for your CUDA version)
(corg) x@y:~$ pip install torch torchvision --index-url https://download.pytorch.org/whl/cu118

Install Dependencies

(corg) x@y:~$ pip install -r requirements.txt

🛠️ Usage

1. Data Preparation

Data is controlled via a CSV file that defines the batch of images to be coregistered. The pipeline reads this file to locate images and understand their spatial footprint.

Directory Structure

root/
├─ data_dir/
│   ├─ csv/ ------------------------- # Input CSVs (e.g., cid_3.csv)
│   └─ rdr/ 
|       ├─ cid_1/
|       ├─ cid_2/
|       ├─ ...
|       └─ cid_n/ ------------------- # Unprocessed JP2 files
|
└─ out_dir/
    └─ cid_1/
        ├─ crop_1.png --------------- # Coregistered .png tile
        ├─ crop_2.png 
        ├─ ... 
        ├─ crop_n.png 
        └─ _debug/
            └─ test_grid/
                ├─ grid.png --------- # Visualization of coregistration grid
                ├─ log.json --------- # Debug info for a given grid
                ├─ offsets.npz ------ # Saved translation matrix
                └─ id.png ----------- # Debug snapshots of failed tiles

Note

Debug and result folders are auto-created in out_dir/ at runtime.

CSV Specifications

The input CSV (e.g., cid_3.csv) must contain the following columns:

Column Name	Description
PRODUCT_ID	The unique identifier for the image. Must match the filename in the rdr director.
MY	The Mars Year (or temporal index) to sort images in chronological order..
C1_X, C1_Y	Map coordinates for Corner 1 (Top-Left).
C2_X, C2_Y	Map coordinates for Corner 2 (Top-Right).
C3_X, C3_Y	Map coordinates for Corner 3 (Bottom-Right).
C4_X, C4_Y	Map coordinates for Corner 4 (Bottom-Left).

2. Configuration (config.yaml)

Edit config/config.yaml to define multiple passes with different resolutions.

# Example Configuration
coregistration_pipeline:
  - name: "coarse_grid"
    grid_res: 1280      # Pixel size of grid cells
    grid_overlap: 640   # Overlap between cells
    min_pts: 10         # Min matches to estimate offset
    max_pts: 2048       # Max K-pts/img to extract
    bucket_size: 160    # Spatial bucketing
    vmf_k: 10           # Compare to 10 neighbors
    max_shift: 640      # Max allowed dx/dy

  - name: "final_warp"
    grid_res: 640     
    grid_overlap: 0
    min_pts: 20         # Min matches to attempt a warp
    max_pts: 6400
    bucket_size: 40
    vmf_k: 20
    max_shift: 50
    warp: True          # Perfrom a warp

3. Running the Pipeline

For quick adjustment, you can perform single-pass debug runs by manually setting the coregistration parameters:

(corg) x@y:~$ python main.py --grid_res 640 --bucket_size 80 --max_shift 100

To execute the entire coregistration pipeline run the script with --run_pipeline flag:

(corg) x@y:~$ python main.py --run_pipeline

🛰️ Example output

The pipeline generates co-registered tiles, typically achieving pixel-perfect precision. However, performance can vary depending on dataset resolution and feature richness. In our specific use case (0.25 m/pixel resolution with variable terrain, shifting backgrounds, and changing lighting conditions), a 320-pixel grid size proved to be the cutoff point for reliable performance. Beyond this limit, the lack of persistent, unique features causes degradation. A useful rule of thumb is to determine the maximum zoom level at which an area retains distinct features suitable for matching; surpassing this often leads to alignment failures.

Figure 3: Final output showing co-registered image patches across four distinct Mars Years (MY 29–33)

Note

For "feature-starved" regions sandwiched between feature-rich areas, a larger grid resolution (coarser grid) can be more effective. Because the transform maintains a degree of global rigidity, aligning the well-defined "flanks" will naturally interpolate and align the featureless middle section.

📝 License

Distributed under the Apache 2.0 License. See LICENSE for more information.

✉️ Contact

Martynchuk Oleksii - martyn.chuckie@gmail.com

🤝 Acknowledgements

This project was made possible thanks to the support and resources provided by:

• Technische Universität Berlin (TU Berlin)
• German Aerospace Center (DLR) Berlin
• HiRISE (High Resolution Imaging Science Experiment) team at the University of Arizona
• HEIBRIDS School for Data Science

Additional thanks to the open‑source community and all contributors who help improve this project.