layout: article.njk

Removing Seams in Multi-Scene Mosaics with Feathering

Removing seams in multi-scene mosaics with feathering requires generating a distance-weighted alpha mask for each overlapping tile, then blending pixel values proportionally to their distance from the nearest valid-data boundary. In Python, this is implemented by computing Euclidean distance transforms on valid-pixel masks, normalizing the resulting gradients, and applying a weighted average across overlap zones: V_final = Σ(V_i * α_i) / Σ(α_i). The technique eliminates hard radiometric discontinuities caused by atmospheric differences, sun-angle variations, and sensor calibration drift without requiring manual seamline digitization or global histogram equalization.

How Distance-Weighted Blending Works

Feathering operates on a purely geometric weighting principle that prioritizes pixel reliability based on spatial position:

  • Interior Priority: Pixels near a tile’s center are assigned higher weights because acquisition edges typically suffer from vignetting, cloud shadows, and atmospheric scattering.
  • Smooth Transitions: In overlap zones, the blending kernel gradually shifts dominance from one scene’s radiometry to another’s, preventing abrupt brightness or color jumps.
  • Deterministic & Lightweight: Unlike histogram matching or seamline optimization, feathering requires no iterative convergence or manual intervention. It executes in linear time relative to pixel count.
  • Precondition Dependency: The method assumes all inputs have already undergone radiometric normalization, orthorectification, and cloud masking. It composites aligned grids; it does not correct geometric misregistration.

Production-Ready Python Implementation

The following function reads spatially aligned rasters, computes per-tile distance weights, and writes a feathered mosaic. It handles arbitrary band counts, respects nodata values, and preserves original geospatial metadata.

import rasterio
import numpy as np
from scipy.ndimage import distance_transform_edt
from pathlib import Path

def feather_mosaic(input_paths: list[str], output_path: str, nodata: float = -9999.0) -> None:
    """
    Create a distance-weighted feathered mosaic from spatially aligned rasters.
    Inputs must share identical CRS, resolution, and grid alignment.
    """
    if len(input_paths) < 2:
        raise ValueError("At least two overlapping rasters are required.")

    mosaic_num = None
    mosaic_den = None
    profile = None

    for path in input_paths:
        with rasterio.open(path) as src:
            # Initialize accumulators on first pass
            if profile is None:
                profile = src.profile.copy()
                profile.update(dtype='float32', nodata=np.nan, count=src.count)
                mosaic_num = np.zeros((src.count, src.height, src.width), dtype=np.float32)
                mosaic_den = np.zeros((src.height, src.width), dtype=np.float32)

            data = src.read().astype(np.float32)
            
            # 2D valid mask: True where ALL bands contain valid data
            valid_mask = np.all(data != nodata, axis=0)
            if not np.any(valid_mask):
                continue

            # Euclidean distance from each valid pixel to the nearest invalid edge
            # ~valid_mask sets valid=0, invalid=1. EDT computes distance from 0 to nearest 1.
            dist = distance_transform_edt(~valid_mask)

            # Normalize weights to [0, 1] range
            alpha = dist / (dist.max() + 1e-6)
            alpha_3d = alpha[np.newaxis, :, :]  # Broadcast across bands

            # Accumulate weighted pixel values and weight sums
            mosaic_num += data * alpha_3d
            mosaic_den += alpha

    # Compute final weighted average, safely handling zero denominators
    with np.errstate(divide='ignore', invalid='ignore'):
        mosaic_final = np.where(mosaic_den > 0, mosaic_num / mosaic_den, np.nan)

    # Write output with updated profile
    with rasterio.open(output_path, 'w', **profile) as dst:
        dst.write(mosaic_final.astype(profile['dtype']))

Key Implementation Notes

  • Mask Logic: np.all(data != nodata, axis=0) ensures only pixels valid across every band participate in blending. This prevents spectral leakage from cloud-masked or saturated bands.
  • Distance Transform: The scipy.ndimage.distance_transform_edt routine (official documentation) computes exact Euclidean distances in O(N) time. It is significantly faster than iterative dilation approaches.
  • Numerical Stability: Adding 1e-6 to dist.max() prevents division-by-zero in edge cases where a tile contains only a single valid pixel. The output dtype is forced to float32 to preserve fractional weights during accumulation.

Memory Management & Large-Scale Processing

Full-scene distance transforms load complete arrays into RAM. For tiles exceeding 5000×5000 pixels or multi-terabyte coverage, sequential accumulation becomes a bottleneck. Implement these scaling strategies:

  1. Windowed I/O: Replace src.read() with src.read(window=...) to process tiles in 1024×1024 blocks. Accumulate numerator/denominator arrays per window, then flush to disk incrementally.
  2. Dask Integration: Swap numpy for dask.array to enable out-of-core computation. Dask’s lazy evaluation automatically chunks distance transforms and distributes blending across available cores.
  3. CRS & Grid Validation: Feathering assumes pixel-perfect alignment. Misaligned grids or differing resolutions will produce ghosting artifacts. Always verify alignment using rasterio.warp.reproject or gdalwarp before compositing.

Workflow Integration & Best Practices

Feathering is not a standalone correction; it is a terminal compositing operation. In production environments, it executes after atmospheric correction, BRDF normalization, and precise geometric registration. When integrated into broader Seamless Mosaicking and Edge Blending pipelines, it guarantees radiometric continuity across acquisition dates and sensor passes.

For teams managing petabyte-scale archives, this technique scales efficiently within automated Satellite Processing Workflows & Index Pipelines, where it runs as a deterministic, idempotent step. Always validate output against ground truth or reference orthomosaics to confirm that weight gradients haven’t suppressed legitimate high-contrast features (e.g., sharp coastlines or urban boundaries). If feature preservation is critical, consider hybrid approaches that combine distance weighting with gradient-domain blending or seamline optimization.