Alignment

The alignment API is the primary Python interface for pairwise raster coregistration.

Import from the top-level package:

from coregix import align_image_pair, AlignmentResult

coregix

Pairwise raster coregistration for geospatial imagery.

AlignmentResult dataclass

Summary of pairwise alignment execution.

Attributes:

Name	Type	Description
output_image_path	str	Final aligned image path.
temp_dir	Optional[str]	Retained temporary directory path when keep_temp_dir=True; otherwise None.

Source code in coregix/pipelines/alignment.py

@dataclass
class AlignmentResult:
    """Summary of pairwise alignment execution.

    Attributes:
        output_image_path: Final aligned image path.
        temp_dir: Retained temporary directory path when ``keep_temp_dir=True``; otherwise ``None``.
    """

    output_image_path: str
    temp_dir: Optional[str]

align_image_pair(moving_image_path, fixed_image_path, output_image_path, *, band_index=0, moving_band_index=None, fixed_band_index=None, moving_nodata=None, fixed_nodata=None, output_nodata=None, min_valid_fraction=0.01, temp_dir=None, keep_temp_dir=False, log_to_console=False, clip_fixed_to_moving=True, output_on_moving_grid=True, trim_edge_invalid=False, edge_trim_depth=8, edge_trim_detection_band_index=0, edge_trim_invalid_below=None, edge_trim_invalid_above=None, enforce_mutual_valid_mask=True, use_edge_proxies=True, split_factor=0, solve_resolution=None)

Coregister a source raster to a reference raster.

Coregix estimates a translation followed by a rigid transform between a source raster and a reference raster. By default, registration is estimated from edge-proxy images over the mutual valid-data overlap. The resulting transform is applied to all source bands and written to a coregistered GeoTIFF.

Parameters:

Name	Type	Description	Default
moving_image_path	str	Path to the source raster that will be transformed.	required
fixed_image_path	str	Path to the reference raster used for alignment.	required
output_image_path	str	Path for the coregistered output raster.	required
band_index	int	0-based band index used for registration when the same band should be used from both rasters.	0
moving_band_index	Optional[int]	Optional 0-based source-raster band index used for registration. Overrides band_index for the source raster.	None
fixed_band_index	Optional[int]	Optional 0-based reference-raster band index used for registration. Overrides band_index for the reference raster.	None
moving_nodata	Optional[float]	Optional source-raster nodata override for mask generation.	None
fixed_nodata	Optional[float]	Optional reference-raster nodata override for mask generation.	None
output_nodata	Optional[float]	Optional output nodata override. Defaults to source nodata, then reference nodata, else 0.	None
min_valid_fraction	float	Minimum valid-mask fraction required in the registration ROI.	0.01
temp_dir	Optional[str]	Optional parent directory for temporary working artifacts.	None
keep_temp_dir	bool	If True, keep the temporary working directory for inspection.	False
log_to_console	bool	If True, emit registration backend logs to stdout.	False
clip_fixed_to_moving	bool	If True, restrict the reference domain to source-raster bounds.	True
output_on_moving_grid	bool	If True, write final output on the source-raster grid with the same transform, dimensions, and pixel size as the source.	True
trim_edge_invalid	bool	If True, post-process the final output by setting pixels adjacent to irregular invalid boundaries to nodata.	False
edge_trim_depth	int	Number of pixels to trim around each invalid boundary.	8
edge_trim_detection_band_index	int	0-based band used to detect edge artifacts.	0
edge_trim_invalid_below	Optional[float]	Optional lower threshold for edge artifact detection.	None
edge_trim_invalid_above	Optional[float]	Optional upper threshold for edge artifact detection.	None
enforce_mutual_valid_mask	bool	If True, constrain both registration masks to the mutual valid-data overlap of source and reference rasters.	True
use_edge_proxies	bool	If True, register on edge-proxy images rather than raw intensities.	True
split_factor	int	Split the registration solve and final resampling domains into 2**split_factor chunks. 0 disables chunking.	0
solve_resolution	Optional[float]	Optional target pixel size, in raster CRS units, for the registration solve. When omitted, the reference-raster ROI resolution is used.	None

Returns:

Type	Description
AlignmentResult	AlignmentResult summary with output path and retained temporary directory,
AlignmentResult	when requested.

Raises:

Type	Description
ValueError	If argument values are out of range or images are incompatible.

Source code in coregix/pipelines/alignment.py

def align_image_pair(
    moving_image_path: str,
    fixed_image_path: str,
    output_image_path: str,
    *,
    band_index: int = 0,
    moving_band_index: Optional[int] = None,
    fixed_band_index: Optional[int] = None,
    moving_nodata: Optional[float] = None,
    fixed_nodata: Optional[float] = None,
    output_nodata: Optional[float] = None,
    min_valid_fraction: float = 0.01,
    temp_dir: Optional[str] = None,
    keep_temp_dir: bool = False,
    log_to_console: bool = False,
    clip_fixed_to_moving: bool = True,
    output_on_moving_grid: bool = True,
    trim_edge_invalid: bool = False,
    edge_trim_depth: int = 8,
    edge_trim_detection_band_index: int = 0,
    edge_trim_invalid_below: Optional[float] = None,
    edge_trim_invalid_above: Optional[float] = None,
    enforce_mutual_valid_mask: bool = True,
    use_edge_proxies: bool = True,
    split_factor: int = 0,
    solve_resolution: Optional[float] = None,
) -> AlignmentResult:
    """Coregister a source raster to a reference raster.

    Coregix estimates a translation followed by a rigid transform between a
    source raster and a reference raster. By default, registration is estimated
    from edge-proxy images over the mutual valid-data overlap. The resulting
    transform is applied to all source bands and written to a coregistered
    GeoTIFF.

    Args:
        moving_image_path: Path to the source raster that will be transformed.
        fixed_image_path: Path to the reference raster used for alignment.
        output_image_path: Path for the coregistered output raster.
        band_index: 0-based band index used for registration when the same band
            should be used from both rasters.
        moving_band_index: Optional 0-based source-raster band index used for
            registration. Overrides ``band_index`` for the source raster.
        fixed_band_index: Optional 0-based reference-raster band index used for
            registration. Overrides ``band_index`` for the reference raster.
        moving_nodata: Optional source-raster nodata override for mask generation.
        fixed_nodata: Optional reference-raster nodata override for mask generation.
        output_nodata: Optional output nodata override. Defaults to source nodata,
            then reference nodata, else ``0``.
        min_valid_fraction: Minimum valid-mask fraction required in the registration ROI.
        temp_dir: Optional parent directory for temporary working artifacts.
        keep_temp_dir: If ``True``, keep the temporary working directory for inspection.
        log_to_console: If ``True``, emit registration backend logs to stdout.
        clip_fixed_to_moving: If ``True``, restrict the reference domain to
            source-raster bounds.
        output_on_moving_grid: If ``True``, write final output on the source-raster
            grid with the same transform, dimensions, and pixel size as the source.
        trim_edge_invalid: If ``True``, post-process the final output by setting pixels
            adjacent to irregular invalid boundaries to nodata.
        edge_trim_depth: Number of pixels to trim around each invalid boundary.
        edge_trim_detection_band_index: 0-based band used to detect edge artifacts.
        edge_trim_invalid_below: Optional lower threshold for edge artifact detection.
        edge_trim_invalid_above: Optional upper threshold for edge artifact detection.
        enforce_mutual_valid_mask: If ``True``, constrain both registration masks
            to the mutual valid-data overlap of source and reference rasters.
        use_edge_proxies: If ``True``, register on edge-proxy images rather than
            raw intensities.
        split_factor: Split the registration solve and final resampling domains
            into ``2**split_factor`` chunks. ``0`` disables chunking.
        solve_resolution: Optional target pixel size, in raster CRS units, for the
            registration solve. When omitted, the reference-raster ROI resolution is used.

    Returns:
        AlignmentResult summary with output path and retained temporary directory,
        when requested.

    Raises:
        ValueError: If argument values are out of range or images are incompatible.
    """
    if band_index < 0:
        raise ValueError("band_index must be >= 0 (0-based).")
    if moving_band_index is not None and moving_band_index < 0:
        raise ValueError("moving_band_index must be >= 0 (0-based).")
    if fixed_band_index is not None and fixed_band_index < 0:
        raise ValueError("fixed_band_index must be >= 0 (0-based).")
    if min_valid_fraction <= 0 or min_valid_fraction > 1:
        raise ValueError("min_valid_fraction must be in (0, 1].")
    if split_factor < 0:
        raise ValueError("split_factor must be >= 0.")
    if edge_trim_depth <= 0:
        raise ValueError("edge_trim_depth must be > 0.")
    if edge_trim_detection_band_index < 0:
        raise ValueError("edge_trim_detection_band_index must be >= 0.")
    if split_factor > 0:
        from coregix.pipelines.alignment_large_main import (
            align_image_pair as align_image_pair_large_main,
        )

        return align_image_pair_large_main(
            moving_image_path,
            fixed_image_path,
            output_image_path,
            band_index=band_index,
            moving_band_index=moving_band_index,
            fixed_band_index=fixed_band_index,
            moving_nodata=moving_nodata,
            fixed_nodata=fixed_nodata,
            output_nodata=output_nodata,
            min_valid_fraction=min_valid_fraction,
            temp_dir=temp_dir,
            keep_temp_dir=keep_temp_dir,
            log_to_console=log_to_console,
            clip_fixed_to_moving=clip_fixed_to_moving,
            output_on_moving_grid=output_on_moving_grid,
            trim_edge_invalid=trim_edge_invalid,
            edge_trim_depth=edge_trim_depth,
            edge_trim_detection_band_index=edge_trim_detection_band_index,
            edge_trim_invalid_below=edge_trim_invalid_below,
            edge_trim_invalid_above=edge_trim_invalid_above,
            enforce_mutual_valid_mask=enforce_mutual_valid_mask,
            use_edge_proxies=use_edge_proxies,
            split_factor=split_factor,
            solve_resolution=solve_resolution,
        )
    if solve_resolution is not None and solve_resolution <= 0:
        raise ValueError("solve_resolution must be > 0 when provided.")

    temp_ctx = None
    work_dir: str
    if keep_temp_dir:
        work_dir = tempfile.mkdtemp(prefix="vhr_align_", dir=temp_dir)
    else:
        temp_ctx = tempfile.TemporaryDirectory(prefix="vhr_align_", dir=temp_dir)
        work_dir = temp_ctx.name

    with rasterio.open(fixed_image_path) as fixed_src, rasterio.open(moving_image_path) as moving_src:
        moving_band_1based = (moving_band_index if moving_band_index is not None else band_index) + 1
        fixed_band_1based = (fixed_band_index if fixed_band_index is not None else band_index) + 1
        if fixed_band_1based > fixed_src.count:
            raise ValueError(
                f"Requested reference band index={fixed_band_1based - 1}, but reference raster has {fixed_src.count} band(s)."
            )
        if moving_band_1based > moving_src.count:
            raise ValueError(
                f"Requested source band index={moving_band_1based - 1}, but source raster has {moving_src.count} band(s)."
            )
        if fixed_src.crs is None or moving_src.crs is None:
            raise ValueError("Both reference and source rasters must have CRS.")
        if fixed_src.crs != moving_src.crs:
            raise ValueError(
                "Reference and source rasters must share the same CRS for tile-window extraction. "
                f"reference={fixed_src.crs}, source={moving_src.crs}"
            )

        moving_nodata_value = _resolve_nodata(moving_src, moving_nodata)
        fixed_nodata_value = _resolve_nodata(fixed_src, fixed_nodata)
        out_nodata = (
            output_nodata
            if output_nodata is not None
            else moving_nodata_value
            if moving_nodata_value is not None
            else fixed_nodata_value
            if fixed_nodata_value is not None
            else 0
        )

        if clip_fixed_to_moving:
            fixed_domain_window = _to_int_window(
                from_bounds(
                    left=moving_src.bounds.left,
                    bottom=moving_src.bounds.bottom,
                    right=moving_src.bounds.right,
                    top=moving_src.bounds.top,
                    transform=fixed_src.transform,
                ),
                max_width=fixed_src.width,
                max_height=fixed_src.height,
            )
            if fixed_domain_window.width <= 0 or fixed_domain_window.height <= 0:
                raise ValueError(
                    "No overlap between source-raster bounds and reference-raster grid."
                )
        else:
            fixed_domain_window = Window(
                col_off=0, row_off=0, width=fixed_src.width, height=fixed_src.height
            )

        os.makedirs(os.path.dirname(output_image_path) or ".", exist_ok=True)
        temp_output_image_path = os.path.join(work_dir, os.path.basename(output_image_path))
        if output_on_moving_grid:
            out_profile = _make_output_profile(
                moving_src.profile,
                count=moving_src.count,
                dtype=moving_src.dtypes[0],
                nodata=out_nodata,
            )
        else:
            out_profile = _make_output_profile(
                fixed_src.profile,
                count=moving_src.count,
                dtype=moving_src.dtypes[0],
                nodata=out_nodata,
                width=int(fixed_domain_window.width),
                height=int(fixed_domain_window.height),
                transform=fixed_src.window_transform(fixed_domain_window),
            )

        fixed_window = fixed_domain_window
        core_out_window = Window(
            col_off=0,
            row_off=0,
            width=int(fixed_domain_window.width),
            height=int(fixed_domain_window.height),
        )
        fixed_window_transform = fixed_src.window_transform(fixed_window)
        fixed_bounds = fixed_src.window_bounds(fixed_window)
        moving_window = _to_int_window(
            from_bounds(
                left=fixed_bounds[0],
                bottom=fixed_bounds[1],
                right=fixed_bounds[2],
                top=fixed_bounds[3],
                transform=moving_src.transform,
            ),
            max_width=moving_src.width,
            max_height=moving_src.height,
        )
        if moving_window.width <= 0 or moving_window.height <= 0:
            raise ValueError("No overlap between reference-raster ROI and source-raster grid.")
        moving_window_transform = moving_src.window_transform(moving_window)
        solve_width, solve_height, solve_transform = _resolve_solve_grid(
            base_transform=fixed_window_transform,
            base_width=int(fixed_window.width),
            base_height=int(fixed_window.height),
            solve_resolution=solve_resolution,
        )

        with rasterio.open(temp_output_image_path, "w+", **out_profile) as out_dst:
            # Preserve radiometric/band metadata from source raster and clear stale stats
            # that can cause misleading display stretches in GIS viewers.
            try:
                out_dst.colorinterp = moving_src.colorinterp
            except Exception:
                pass
            try:
                out_dst.scales = moving_src.scales
                out_dst.offsets = moving_src.offsets
            except Exception:
                pass
            for b in range(1, moving_src.count + 1):
                desc = moving_src.descriptions[b - 1]
                if desc:
                    out_dst.set_band_description(b, desc)
                band_tags = moving_src.tags(b).copy()
                # Remove stale stats tags from source and output; they are often invalid
                # after reprojection/warping and can distort visualization.
                for k in list(band_tags.keys()):
                    if k.upper().startswith("STATISTICS_"):
                        band_tags.pop(k, None)
                out_dst.update_tags(b, **band_tags)
                out_dst.update_tags(
                    b,
                    STATISTICS_MINIMUM="",
                    STATISTICS_MAXIMUM="",
                    STATISTICS_MEAN="",
                    STATISTICS_STDDEV="",
                )

            if output_on_moving_grid:
                # Preserve source-raster pixels outside the reference overlap domain.
                for b in range(1, moving_src.count + 1):
                    for _, block_window in out_dst.block_windows(b):
                        src_block = moving_src.read(b, window=block_window)
                        out_dst.write(src_block.astype(out_profile["dtype"]), b, window=block_window)
            else:
                for b in range(1, moving_src.count + 1):
                    for _, block_window in out_dst.block_windows(b):
                        fill = np.full(
                            (int(block_window.height), int(block_window.width)),
                            out_nodata,
                            dtype=out_profile["dtype"],
                        )
                        out_dst.write(fill, b, window=block_window)

            fixed_band = fixed_src.read(fixed_band_1based, window=fixed_window)
            moving_band = moving_src.read(moving_band_1based, window=moving_window)
            fixed_valid = fixed_src.read_masks(fixed_band_1based, window=fixed_window) > 0
            moving_valid = moving_src.read_masks(moving_band_1based, window=moving_window) > 0
            if fixed_nodata_value is not None:
                fixed_valid &= fixed_band != fixed_nodata_value
            if moving_nodata_value is not None:
                moving_valid &= moving_band != moving_nodata_value

            fixed_reg_data = np.full(
                (int(solve_height), int(solve_width)),
                fixed_nodata_value if fixed_nodata_value is not None else out_nodata,
                dtype=np.float32,
            )
            reproject(
                source=fixed_band.astype(np.float32),
                destination=fixed_reg_data,
                src_transform=fixed_window_transform,
                src_crs=fixed_src.crs,
                dst_transform=solve_transform,
                dst_crs=fixed_src.crs,
                src_nodata=fixed_nodata_value,
                dst_nodata=fixed_nodata_value if fixed_nodata_value is not None else out_nodata,
                resampling=Resampling.nearest,
            )
            fixed_valid_reprojected = np.zeros((int(solve_height), int(solve_width)), dtype=np.uint8)
            reproject(
                source=fixed_valid.astype(np.uint8),
                destination=fixed_valid_reprojected,
                src_transform=fixed_window_transform,
                src_crs=fixed_src.crs,
                dst_transform=solve_transform,
                dst_crs=fixed_src.crs,
                src_nodata=0,
                dst_nodata=0,
                resampling=Resampling.nearest,
            )
            moving_valid_reprojected = np.zeros((int(solve_height), int(solve_width)), dtype=np.uint8)
            reproject(
                source=moving_valid.astype(np.uint8),
                destination=moving_valid_reprojected,
                src_transform=moving_window_transform,
                src_crs=moving_src.crs,
                dst_transform=solve_transform,
                dst_crs=fixed_src.crs,
                src_nodata=0,
                dst_nodata=0,
                resampling=Resampling.nearest,
            )
            moving_mask_for_elastix: np.ndarray
            fixed_mask_for_elastix: np.ndarray
            moving_on_fixed = np.full(
                (int(solve_height), int(solve_width)),
                moving_nodata_value if moving_nodata_value is not None else out_nodata,
                dtype=np.float32,
            )
            reproject(
                source=moving_band.astype(np.float32),
                destination=moving_on_fixed,
                src_transform=moving_window_transform,
                src_crs=moving_src.crs,
                dst_transform=solve_transform,
                dst_crs=fixed_src.crs,
                src_nodata=moving_nodata_value,
                dst_nodata=moving_nodata_value if moving_nodata_value is not None else out_nodata,
                resampling=Resampling.nearest,
            )
            moving_on_fixed_valid = moving_valid_reprojected > 0
            fixed_valid_for_registration = fixed_valid_reprojected > 0
            fixed_mask_for_elastix = fixed_valid_for_registration.astype(np.uint8)
            moving_mask_for_elastix = moving_on_fixed_valid.astype(np.uint8)
            if use_edge_proxies:
                fixed_reg_data = _edge_proxy(fixed_reg_data, fixed_valid_for_registration)
                moving_reg_data = _edge_proxy(moving_on_fixed, moving_on_fixed_valid)
            else:
                moving_reg_data = moving_on_fixed
            if enforce_mutual_valid_mask:
                mutual = (fixed_mask_for_elastix > 0) & (moving_mask_for_elastix > 0)
                fixed_mask_for_elastix = mutual.astype(np.uint8)
                moving_mask_for_elastix = mutual.astype(np.uint8)

            min_valid_pixels = int(max(1, min_valid_fraction * (solve_width * solve_height)))
            if int(fixed_mask_for_elastix.sum()) < min_valid_pixels:
                raise ValueError("Insufficient valid reference-raster support in the registration ROI.")
            if int(moving_mask_for_elastix.sum()) < min_valid_pixels:
                raise ValueError("Insufficient valid source-raster support in the registration ROI.")

            fixed_reg_path = os.path.join(work_dir, "fixed_reg.tif")
            moving_reg_path = os.path.join(work_dir, "moving_reg.tif")
            fixed_mask_path = os.path.join(work_dir, "fixed_mask.tif")
            moving_mask_path = os.path.join(work_dir, "moving_mask.tif")

            _write_single_band_tif(
                fixed_reg_path,
                fixed_reg_data.astype("float32"),
                crs=fixed_src.crs,
                transform=solve_transform,
                dtype="float32",
                nodata=fixed_nodata_value,
            )
            _write_single_band_tif(
                moving_reg_path,
                moving_reg_data.astype("float32"),
                crs=fixed_src.crs,
                transform=solve_transform,
                dtype="float32",
                nodata=moving_nodata_value,
            )
            _write_single_band_tif(
                fixed_mask_path,
                fixed_mask_for_elastix.astype("uint8"),
                crs=fixed_src.crs,
                transform=solve_transform,
                dtype="uint8",
                nodata=0,
            )
            _write_single_band_tif(
                moving_mask_path,
                moving_mask_for_elastix.astype("uint8"),
                crs=fixed_src.crs,
                transform=solve_transform,
                dtype="uint8",
                nodata=0,
            )

            try:
                transform_parameter_object = estimate_elastix_transform(
                    fixed_image_path=fixed_reg_path,
                    moving_image_path=moving_reg_path,
                    parameter_map=DEFAULT_ALIGNMENT_PARAMETER_MAPS,
                    force_nearest_resample=True,
                    fixed_mask_path=fixed_mask_path,
                    moving_mask_path=moving_mask_path,
                    log_to_console=log_to_console,
                )
            except Exception as exc:
                raise RuntimeError(f"Elastix registration failed: {exc}") from exc

            if output_on_moving_grid:
                deformation_field = deformation_field_from_transform(
                    fixed_reg_path,
                    transform_parameter_object,
                    output_directory=work_dir,
                ).astype(np.float32)
                fixed_dx = deformation_field[..., 0]
                fixed_dy = deformation_field[..., 1]

            for b in range(1, moving_src.count + 1):
                moving_band_data = moving_src.read(b, window=moving_window).astype(np.float32)
                if output_on_moving_grid:
                    moving_band_valid = moving_src.read_masks(b, window=moving_window).astype(np.float32)
                    if moving_nodata_value is not None:
                        moving_band_valid *= (moving_band_data != moving_nodata_value).astype(np.float32)
                    block_width, block_height = out_dst.block_shapes[b - 1]
                    for row_off in range(0, int(moving_window.height), int(block_height)):
                        for col_off in range(0, int(moving_window.width), int(block_width)):
                            win_w = min(int(block_width), int(moving_window.width) - col_off)
                            win_h = min(int(block_height), int(moving_window.height) - row_off)
                            block_window = Window(col_off=col_off, row_off=row_off, width=win_w, height=win_h)
                            block_transform = rasterio.windows.transform(block_window, moving_window_transform)
                            x_world, y_world = _pixel_centers_world(block_transform, win_h, win_w)
                            solve_rows, solve_cols = _world_to_array_coords(
                                solve_transform,
                                x_world,
                                y_world,
                            )
                            dx_block, dx_valid = _sample_bilinear(
                                fixed_dx,
                                solve_rows,
                                solve_cols,
                                fill_value=0.0,
                            )
                            dy_block, dy_valid = _sample_bilinear(
                                fixed_dy,
                                solve_rows,
                                solve_cols,
                                fill_value=0.0,
                            )
                            field_valid = dx_valid & dy_valid
                            source_solve_rows = solve_rows + dy_block
                            source_solve_cols = solve_cols + dx_block
                            source_x_world, source_y_world = _array_to_world(
                                solve_transform,
                                source_solve_rows,
                                source_solve_cols,
                            )
                            source_moving_rows, source_moving_cols = _world_to_array_coords(
                                moving_window_transform,
                                source_x_world,
                                source_y_world,
                            )
                            remapped_block, moving_valid = _sample_bilinear(
                                moving_band_data,
                                source_moving_rows,
                                source_moving_cols,
                                fill_value=out_nodata,
                            )
                            sampled_mask, mask_valid = _sample_bilinear(
                                moving_band_valid,
                                source_moving_rows,
                                source_moving_cols,
                                fill_value=0.0,
                            )
                            source_block_window = Window(
                                col_off=int(moving_window.col_off + col_off),
                                row_off=int(moving_window.row_off + row_off),
                                width=win_w,
                                height=win_h,
                            )
                            valid = field_valid & moving_valid & mask_valid & (sampled_mask > 0.0)
                            combined = np.where(valid, remapped_block, out_nodata)
                            out_dst.write(
                                combined.astype(out_profile["dtype"]),
                                b,
                                window=source_block_window,
                            )
                else:
                    # Non-moving-grid output uses the solve grid, then remaps to fixed output grid.
                    moving_band_on_fixed = np.full(
                        (int(solve_height), int(solve_width)),
                        out_nodata,
                        dtype=np.float32,
                    )
                    reproject(
                        source=moving_band_data,
                        destination=moving_band_on_fixed,
                        src_transform=moving_window_transform,
                        src_crs=moving_src.crs,
                        dst_transform=solve_transform,
                        dst_crs=fixed_src.crs,
                        src_nodata=moving_nodata_value,
                        dst_nodata=out_nodata,
                        resampling=Resampling.nearest,
                    )
                    warped_full = apply_elastix_transform_array(
                        moving_image=moving_band_on_fixed,
                        transform_parameter_object=transform_parameter_object,
                        log_to_console=log_to_console,
                    )
                    if solve_width == int(fixed_window.width) and solve_height == int(fixed_window.height):
                        out_dst.write(warped_full.astype(out_profile["dtype"]), b, window=core_out_window)
                    else:
                        warped_on_fixed = np.full(
                            (int(fixed_window.height), int(fixed_window.width)),
                            out_nodata,
                            dtype=np.float32,
                        )
                        reproject(
                            source=warped_full,
                            destination=warped_on_fixed,
                            src_transform=solve_transform,
                            src_crs=fixed_src.crs,
                            src_nodata=out_nodata,
                            dst_transform=fixed_window_transform,
                            dst_crs=fixed_src.crs,
                            dst_nodata=out_nodata,
                            resampling=Resampling.bilinear,
                        )
                        out_dst.write(warped_on_fixed.astype(out_profile["dtype"]), b, window=core_out_window)

        if trim_edge_invalid:
            from coregix.postprocess import trim_edge_invalid_pixels

            trim_edge_invalid_pixels(
                input_image_path=temp_output_image_path,
                output_image_path=output_image_path,
                edge_depth=edge_trim_depth,
                detection_band_index=edge_trim_detection_band_index,
                invalid_below=edge_trim_invalid_below,
                invalid_above=edge_trim_invalid_above,
                nodata_value=out_nodata,
            )
        else:
            shutil.copyfile(temp_output_image_path, output_image_path)

    if temp_ctx is not None:
        temp_ctx.cleanup()
        kept_temp = None
    else:
        kept_temp = work_dir

    return AlignmentResult(
        output_image_path=output_image_path,
        temp_dir=kept_temp,
    )