Matching Algorithms

global_regression(input_images, output_images, *, calculation_dtype='float32', output_dtype=None, vector_mask=None, debug_logs=False, custom_nodata_value=None, image_parallel_workers=None, window_parallel_workers=None, window_size=None, save_as_cog=False, specify_model_images=None, custom_mean_factor=1.0, custom_std_factor=1.0, save_adjustments=None, load_adjustments=None)

Performs global radiometric normalization across overlapping images using least squares regression.

Parameters:

Name	Type	Description	Default
input_images	(str | List[str], required)	Defines input files from a glob path, folder, or list of paths. Specify like: "/input/files/.tif", "/input/folder" (assumes .tif), ["/input/one.tif", "/input/two.tif"].	required
output_images	(str | List[str], required)	Defines output files from a template path, folder, or list of paths (with the same length as the input). Specify like: "/input/files/$.tif", "/input/folder" (assumes $_Global.tif), ["/input/one.tif", "/input/two.tif"].	required
calculation_dtype	str	Data type used for internal calculations. Defaults to "float32".	'float32'
output_dtype	str | None	Data type for output rasters. Defaults to input image dtype.	None
vector_mask	Tuple[Literal['include', 'exclude'], str, Optional[str]] | None	Mask to limit stats calculation to specific areas in the format of a tuple with two or three items: literal "include" or "exclude" the mask area, str path to the vector file, optional str of field name in vector file that includes (can be substring) input image name to filter geometry by. Loaded stats won't have this applied to them. The matching solution is still applied to these areas in the output. Defaults to None for no mask.	None
debug_logs	bool	If True, prints debug information and constraint matrices. Defaults to False.	False
custom_nodata_value	float | int | None	Overrides detected NoData value. Defaults to None.	None
image_parallel_workers	Tuple[Literal["process", "thread"], Literal["cpu"] | int] | None = None	Parallelization strategy at the image level. Provide a tuple like ("process", "cpu") to use multiprocessing with all available cores. Threads are supported too. Set to None to disable.	None
window_parallel_workers	Tuple[Literal["process"], Literal["cpu"] | int] | None = None	Parallelization strategy at the window level within each image. Same format as image_parallel_workers. Threads are not supported. Set to None to disable.	None
window_size	int | Tuple[int, int] | Literal['internal'] | None	Tile size for reading and writing: int for square tiles, tuple for (width, height), "internal" to use raster's native tiling, or None for full image. "internal" enables efficient streaming from COGs.	None
save_as_cog	bool	If True, saves output as a Cloud-Optimized GeoTIFF using proper band and block order.	False
specify_model_images	Tuple[Literal['exclude', 'include'], List[str]] | None	First item in tuples sets weather to 'include' or 'exclude' the listed images from model building statistics. Second item is the list of image names (without their extension) to apply criteria to. For example, if this param is only set to 'include' one image, all other images will be matched to that one image. Defaults to no exclusion.	None
custom_mean_factor	float	Weight for mean constraints in regression. Defaults to 1.0.	1.0
custom_std_factor	float	Weight for standard deviation constraints in regression. Defaults to 1.0.	1.0
save_adjustments	str | None	The output path of a .json file to save adjustments parameters. Defaults to not saving.	None
load_adjustments	str | None	If set, loads saved whole and overlapping statistics only for images that exist in the .json file. Other images will still have their statistics calculated. Defaults to None.	None

Returns:

Type	Description
list	List[str]: Paths to the globally adjusted output raster images.

Source code in spectralmatch/match/global_regression.py

def global_regression(
    input_images: Universal.SearchFolderOrListFiles,
    output_images: Universal.CreateInFolderOrListFiles,
    *,
    calculation_dtype: Universal.CalculationDtype = "float32",
    output_dtype: Universal.CustomOutputDtype = None,
    vector_mask: Universal.VectorMask = None,
    debug_logs: Universal.DebugLogs = False,
    custom_nodata_value: Universal.CustomNodataValue = None,
    image_parallel_workers: Universal.ImageParallelWorkers = None,
    window_parallel_workers: Universal.WindowParallelWorkers = None,
    window_size: Universal.WindowSize = None,
    save_as_cog: Universal.SaveAsCog = False,
    specify_model_images: Match.SpecifyModelImages = None,
    custom_mean_factor: float = 1.0,
    custom_std_factor: float = 1.0,
    save_adjustments: str | None = None,
    load_adjustments: str | None = None,
) -> list:
    """
    Performs global radiometric normalization across overlapping images using least squares regression.

    Args:
        input_images (str | List[str], required): Defines input files from a glob path, folder, or list of paths. Specify like: "/input/files/*.tif", "/input/folder" (assumes *.tif), ["/input/one.tif", "/input/two.tif"].
        output_images (str | List[str], required): Defines output files from a template path, folder, or list of paths (with the same length as the input). Specify like: "/input/files/$.tif", "/input/folder" (assumes $_Global.tif), ["/input/one.tif", "/input/two.tif"].
        calculation_dtype (str, optional): Data type used for internal calculations. Defaults to "float32".
        output_dtype (str | None, optional): Data type for output rasters. Defaults to input image dtype.
        vector_mask (Tuple[Literal["include", "exclude"], str, Optional[str]] | None): Mask to limit stats calculation to specific areas in the format of a tuple with two or three items: literal "include" or "exclude" the mask area, str path to the vector file, optional str of field name in vector file that *includes* (can be substring) input image name to filter geometry by. Loaded stats won't have this applied to them. The matching solution is still applied to these areas in the output. Defaults to None for no mask.
        debug_logs (bool, optional): If True, prints debug information and constraint matrices. Defaults to False.
        custom_nodata_value (float | int | None, optional): Overrides detected NoData value. Defaults to None.
        image_parallel_workers (Tuple[Literal["process", "thread"], Literal["cpu"] | int] | None = None): Parallelization strategy at the image level. Provide a tuple like ("process", "cpu") to use multiprocessing with all available cores. Threads are supported too. Set to None to disable.
        window_parallel_workers (Tuple[Literal["process"], Literal["cpu"] | int] | None = None): Parallelization strategy at the window level within each image. Same format as image_parallel_workers. Threads are not supported. Set to None to disable.
        window_size (int | Tuple[int, int] | Literal["internal"] | None): Tile size for reading and writing: int for square tiles, tuple for (width, height), "internal" to use raster's native tiling, or None for full image. "internal" enables efficient streaming from COGs.
        save_as_cog (bool): If True, saves output as a Cloud-Optimized GeoTIFF using proper band and block order.
        specify_model_images (Tuple[Literal["exclude", "include"], List[str]] | None ): First item in tuples sets weather to 'include' or 'exclude' the listed images from model building statistics. Second item is the list of image names (without their extension) to apply criteria to. For example, if this param is only set to 'include' one image, all other images will be matched to that one image. Defaults to no exclusion.
        custom_mean_factor (float, optional): Weight for mean constraints in regression. Defaults to 1.0.
        custom_std_factor (float, optional): Weight for standard deviation constraints in regression. Defaults to 1.0.
        save_adjustments (str | None, optional): The output path of a .json file to save adjustments parameters. Defaults to not saving.
        load_adjustments (str | None, optional): If set, loads saved whole and overlapping statistics only for images that exist in the .json file. Other images will still have their statistics calculated. Defaults to None.

    Returns:
        List[str]: Paths to the globally adjusted output raster images.
    """

    print("Start global regression")

    # Validate params
    Universal.validate(
        input_images=input_images,
        output_images=output_images,
        save_as_cog=save_as_cog,
        debug_logs=debug_logs,
        vector_mask=vector_mask,
        window_size=window_size,
        custom_nodata_value=custom_nodata_value,
        image_parallel_workers=image_parallel_workers,
        window_parallel_workers=window_parallel_workers,
        calculation_dtype=calculation_dtype,
        output_dtype=output_dtype,
    )

    Match.validate_match(
        specify_model_images=specify_model_images,
    )

    Match.validate_global_regression(
        custom_mean_factor=custom_mean_factor,
        custom_std_factor=custom_std_factor,
        save_adjustments=save_adjustments,
        load_adjustments=load_adjustments,
    )

    # Input and output paths
    input_image_paths = _resolve_paths(
        "search", input_images, kwargs={"default_file_pattern": "*.tif"}
    )
    output_image_paths = _resolve_paths(
        "create",
        output_images,
        kwargs={
            "paths_or_bases": input_image_paths,
            "default_file_pattern": "$_Global.tif",
        },
    )

    if debug_logs:
        print(f"Input images: {input_image_paths}")
    if debug_logs:
        print(f"Output images: {output_image_paths}")

    input_image_names = [
        os.path.splitext(os.path.basename(p))[0] for p in input_image_paths
    ]
    input_image_path_pairs = dict(zip(input_image_names, input_image_paths))
    output_image_path_pairs = dict(zip(input_image_names, output_image_paths))

    # Check raster requirements
    _check_raster_requirements(
        input_image_paths,
        debug_logs,
        check_geotransform=True,
        check_crs=True,
        check_bands=True,
        check_nodata=True,
    )

    nodata_val = _get_nodata_value(
        list(input_image_path_pairs.values()), custom_nodata_value
    )

    # Determine multiprocessing and worker count
    image_parallel, image_backend, image_max_workers = _resolve_parallel_config(
        image_parallel_workers
    )
    window_parallel, window_backend, window_max_workers = _resolve_parallel_config(
        window_parallel_workers
    )

    # Find loaded and input files if load adjustments
    loaded_model = {}
    if load_adjustments:
        with open(load_adjustments, "r") as f:
            loaded_model = json.load(f)
        _validate_adjustment_model_structure(loaded_model)
        loaded_names = set(loaded_model.keys())
        input_names = set(input_image_names)
    else:
        loaded_names = set([])
        input_names = set(input_image_names)

    matched = input_names & loaded_names
    only_loaded = loaded_names - input_names
    only_input = input_names - loaded_names
    if debug_logs:
        print(
            f"Total images: input images: {len(input_names)}, loaded images {len(loaded_names)}: "
        )
        print(
            f"    Matched adjustments (to override) ({len(matched)}):", sorted(matched)
        )
        print(
            f"    Only in loaded adjustments (to add) ({len(only_loaded)}):",
            sorted(only_loaded),
        )
        print(
            f"    Only in input (to calculate) ({len(only_input)}):", sorted(only_input)
        )

    # Find images to include in model
    included_names = list(matched | only_loaded | only_input)
    if specify_model_images:
        mode, names = specify_model_images
        name_set = set(names)
        if mode == "include":
            included_names = [n for n in input_image_names if n in name_set]
        elif mode == "exclude":
            included_names = [n for n in input_image_names if n not in name_set]
        excluded_names = [n for n in input_image_names if n not in included_names]
    if debug_logs:
        print("Images to influence the model:")
        print(
            f"    Included in model ({len(included_names)}): {sorted(included_names)}"
        )
        if specify_model_images:
            print(
                f"    Excluded from model ({len(excluded_names)}): {sorted(excluded_names)}"
            )
        else:
            print(f"    Excluded from model (0): []")

    if debug_logs:
        print("Calculating statistics")
    with rasterio.open(list(input_image_path_pairs.values())[0]) as src:
        num_bands = src.count

    # Get images bounds
    all_bounds = {}
    for name, path in input_image_path_pairs.items():
        with rasterio.open(path) as ds:
            all_bounds[name] = ds.bounds

    # Overlap stats
    overlapping_pairs = _find_overlaps(all_bounds)
    all_overlap_stats = {}

    # Load overlap stats
    if load_adjustments:
        for name_i, model_entry in loaded_model.items():
            if name_i not in input_image_path_pairs:
                continue

            for name_j, bands in model_entry.get("overlap_stats", {}).items():
                if name_j not in input_image_path_pairs:
                    continue

                all_overlap_stats.setdefault(name_i, {})[name_j] = {
                    int(k.split("_")[1]): {
                        "mean": bands[k]["mean"],
                        "std": bands[k]["std"],
                        "size": bands[k]["size"],
                    }
                    for k in bands
                }

    # Calculate overlap stats
    parallel_args = [
        (
            window_parallel,
            window_max_workers,
            window_backend,
            num_bands,
            input_image_path_pairs[name_i],
            input_image_path_pairs[name_j],
            name_i,
            name_j,
            all_bounds[name_i],
            all_bounds[name_j],
            nodata_val,
            nodata_val,
            vector_mask,
            window_size,
            debug_logs,
        )
        for name_i, name_j in overlapping_pairs
        if name_i not in loaded_model
        or name_j not in loaded_model.get(name_i, {}).get("overlap_stats", {})
    ]

    if image_parallel:
        with _get_executor(image_backend, image_max_workers) as executor:
            futures = [
                executor.submit(_overlap_stats_process_image, *args)
                for args in parallel_args
            ]
            for future in as_completed(futures):
                stats = future.result()
                for outer, inner in stats.items():
                    all_overlap_stats.setdefault(outer, {}).update(inner)
    else:
        for args in parallel_args:
            stats = _overlap_stats_process_image(*args)
            for outer, inner in stats.items():
                all_overlap_stats.setdefault(outer, {}).update(inner)

    # Load whole stats
    all_whole_stats = {
        name: {
            int(k.split("_")[1]): {
                "mean": loaded_model[name]["whole_stats"][k]["mean"],
                "std": loaded_model[name]["whole_stats"][k]["std"],
                "size": loaded_model[name]["whole_stats"][k]["size"],
            }
            for k in loaded_model[name]["whole_stats"]
        }
        for name in input_image_path_pairs
        if name in loaded_model
    }

    # Calculate whole stats
    parallel_args = [
        (
            window_parallel,
            window_max_workers,
            window_backend,
            image_path,
            nodata_val,
            num_bands,
            image_name,
            vector_mask,
            window_size,
            debug_logs,
        )
        for image_name, image_path in input_image_path_pairs.items()
        if image_name not in loaded_model
    ]

    # Compute whole stats
    if image_parallel:
        with _get_executor(image_backend, image_max_workers) as executor:
            futures = [
                executor.submit(_whole_stats_process_image, *args)
                for args in parallel_args
            ]
            for future in as_completed(futures):
                result = future.result()
                all_whole_stats.update(result)
    else:
        for args in parallel_args:
            result = _whole_stats_process_image(*args)
            all_whole_stats.update(result)

    # Get image names
    all_image_names = list(dict.fromkeys(input_image_names + list(loaded_model.keys())))
    num_total = len(all_image_names)

    # Print model sources
    if debug_logs:
        print(
            f"\nCreating model for {len(all_image_names)} total images from {len(included_names)} included:"
        )
        print(f"    {'ID':<4}\t{'Source':<6}\t{'Inclusion':<8}\tName")
        for i, name in enumerate(all_image_names):
            source = "load" if name in (matched | only_loaded) else "calc"
            included = "incl" if name in included_names else "excl"
            print(f"    {i:<4}\t{source:<6}\t{included:<8}\t{name}")

    # Build model
    all_params = _solve_global_model(
        num_bands,
        num_total,
        all_image_names,
        included_names,
        input_image_names,
        all_overlap_stats,
        all_whole_stats,
        custom_mean_factor,
        custom_std_factor,
        overlapping_pairs,
        debug_logs,
    )

    # Save adjustments
    if save_adjustments:
        _save_adjustments(
            save_path=save_adjustments,
            input_image_names=list(input_image_path_pairs.keys()),
            all_params=all_params,
            all_whole_stats=all_whole_stats,
            all_overlap_stats=all_overlap_stats,
            num_bands=num_bands,
            calculation_dtype=calculation_dtype,
        )

    # Apply corrections
    if debug_logs:
        print(f"Apply adjustments and saving results for:")
    parallel_args = [
        (
            name,
            img_path,
            output_image_path_pairs[name],
            np.array([all_params[b, 2 * idx, 0] for b in range(num_bands)]),
            np.array([all_params[b, 2 * idx + 1, 0] for b in range(num_bands)]),
            num_bands,
            nodata_val,
            window_size,
            calculation_dtype,
            output_dtype,
            window_parallel,
            window_backend,
            window_max_workers,
            save_as_cog,
            debug_logs,
        )
        for idx, (name, img_path) in enumerate(input_image_path_pairs.items())
    ]

    if image_parallel:
        with _get_executor(image_backend, image_max_workers) as executor:
            futures = [
                executor.submit(_apply_adjustments_process_image, *args)
                for args in parallel_args
            ]
            for future in as_completed(futures):
                future.result()
    else:
        for args in parallel_args:
            _apply_adjustments_process_image(*args)

    return output_image_paths

local_block_adjustment(input_images, output_images, *, calculation_dtype='float32', output_dtype=None, vector_mask=None, debug_logs=False, custom_nodata_value=None, image_parallel_workers=None, window_parallel_workers=None, window_size=None, save_as_cog=False, number_of_blocks=100, alpha=1.0, correction_method='gamma', save_block_maps=None, load_block_maps=None, override_bounds_canvas_coords=None, block_valid_pixel_threshold=0.001)

Performs local radiometric adjustment on a set of raster images using block-based statistics.

Parameters:

Name	Type	Description	Default
input_images	(str | List[str], required)	Defines input files from a glob path, folder, or list of paths. Specify like: "/input/files/.tif", "/input/folder" (assumes .tif), ["/input/one.tif", "/input/two.tif"].	required
output_images	(str | List[str], required)	Defines output files from a template path, folder, or list of paths (with the same length as the input). Specify like: "/input/files/$.tif", "/input/folder" (assumes $_Local.tif), ["/input/one.tif", "/input/two.tif"].	required
calculation_dtype	str	Precision for internal calculations. Defaults to "float32".	'float32'
output_dtype	str | None	Data type for output rasters. Defaults to input image dtype.	None
vector_mask	Tuple[Literal['include', 'exclude'], str, Optional[str]] | None	A mask limiting pixels to include when calculating stats for each block in the format of a tuple with two or three items: literal "include" or "exclude" the mask area, str path to the vector file, optional str of field name in vector file that includes (can be substring) input image name to filter geometry by. It is only applied when calculating local blocks, as the reference map is calculated as the mean of all local blocks. Loaded block maps won't have this applied unless it was used when calculating them. The matching solution is still applied to these areas in the output. Defaults to None for no mask.	None
debug_logs	bool	If True, prints progress. Defaults to False.	False
custom_nodata_value	float | int | None	Overrides detected NoData value. Defaults to None.	None
image_parallel_workers	Tuple[Literal["process", "thread"], Literal["cpu"] | int] | None = None	Parallelization strategy at the image level. Provide a tuple like ("process", "cpu") to use multiprocessing with all available cores. Threads are supported too. Set to None to disable.	None
window_parallel_workers	Tuple[Literal["process"], Literal["cpu"] | int] | None = None	Parallelization strategy at the window level within each image. Same format as image_parallel_workers. Threads are not supported. Set to None to disable.	None
window_size	int | Tuple[int, int] | Literal['block'] | None	Tile size for processing: int for square tiles, (width, height) for custom size, or "block" to set as the size of the block map, None for full image. Defaults to None.	None
save_as_cog	bool	If True, saves as COG. Defaults to False.	False
number_of_blocks	int | tuple | Literal['coefficient_of_variation']	int as a target of blocks per image, tuple to set manually set total blocks width and height, coefficient_of_variation to find the number of blocks based on this metric.	100
alpha	float	Blending factor between reference and local means. Defaults to 1.0.	1.0
correction_method	Literal['gamma', 'linear']	Local correction method. Defaults to "gamma".	'gamma'
save_block_maps	tuple(str, str) | None	If enabled, saves block maps for review, to resume processing later, or to add additional images to the reference map. - First str is the path to save the global block map. - Second str is the path to save the local block maps, which must include "$" which will be replaced my the image name (because there are multiple local maps).	None
load_block_maps	Tuple[str, List[str]] | Tuple[str, None] | Tuple[None, List[str]] | None	Controls loading of precomputed block maps. Can be one of: - Tuple[str, List[str]]: Load both reference and local block maps. - Tuple[str, None]: Load only the reference block map. - Tuple[None, List[str]]: Load only the local block maps. - None: Do not load any block maps. This supports partial or full reuse of precomputed block maps: - Local block maps will still be computed for each input image that is not linked to a local block map by the images name being included in the local block maps name (file name). - The reference block map will only be calculated (mean of all local blocks) if not set. - The reference map defines the reference block statistics and the local maps define per-image local block statistics. - Both reference and local maps must have the same canvas extent and dimensions which will be used to set those values.	None
override_bounds_canvas_coords	Tuple[float, float, float, float] | None	Manually set (min_x, min_y, max_x, max_y) bounds to override the computed/loaded canvas extent. If you wish to have a larger extent than the current images, you can manually set this, along with setting a fixed number of blocks, to anticipate images will expand beyond the current extent.	None
block_valid_pixel_threshold	float	Minimum fraction of valid pixels required to include a block (0–1).	0.001

Returns:

Type	Description
list	List[str]: Paths to the locally adjusted output raster images.

Source code in spectralmatch/match/local_block_adjustment.py

def local_block_adjustment(
    input_images: Universal.SearchFolderOrListFiles,
    output_images: Universal.CreateInFolderOrListFiles,
    *,
    calculation_dtype: Universal.CalculationDtype = "float32",
    output_dtype: Universal.CustomOutputDtype = None,
    vector_mask: Universal.VectorMask = None,
    debug_logs: Universal.DebugLogs = False,
    custom_nodata_value: Universal.CustomNodataValue = None,
    image_parallel_workers: Universal.ImageParallelWorkers = None,
    window_parallel_workers: Universal.WindowParallelWorkers = None,
    window_size: Universal.WindowSizeWithBlock = None,
    save_as_cog: Universal.SaveAsCog = False,
    number_of_blocks: int | Tuple[int, int] | Literal["coefficient_of_variation"] = 100,
    alpha: float = 1.0,
    correction_method: Literal["gamma", "linear"] = "gamma",
    save_block_maps: Tuple[str, str] | None = None,
    load_block_maps: (
        Tuple[str, List[str]] | Tuple[str, None] | Tuple[None, List[str]] | None
    ) = None,
    override_bounds_canvas_coords: Tuple[float, float, float, float] | None = None,
    block_valid_pixel_threshold: float = 0.001,
) -> list:
    """
    Performs local radiometric adjustment on a set of raster images using block-based statistics.

    Args:
        input_images (str | List[str], required): Defines input files from a glob path, folder, or list of paths. Specify like: "/input/files/*.tif", "/input/folder" (assumes *.tif), ["/input/one.tif", "/input/two.tif"].
        output_images (str | List[str], required): Defines output files from a template path, folder, or list of paths (with the same length as the input). Specify like: "/input/files/$.tif", "/input/folder" (assumes $_Local.tif), ["/input/one.tif", "/input/two.tif"].
        calculation_dtype (str, optional): Precision for internal calculations. Defaults to "float32".
        output_dtype (str | None, optional): Data type for output rasters. Defaults to input image dtype.
        vector_mask (Tuple[Literal["include", "exclude"], str, Optional[str]] | None): A mask limiting pixels to include when calculating stats for each block in the format of a tuple with two or three items: literal "include" or "exclude" the mask area, str path to the vector file, optional str of field name in vector file that *includes* (can be substring) input image name to filter geometry by. It is only applied when calculating local blocks, as the reference map is calculated as the mean of all local blocks. Loaded block maps won't have this applied unless it was used when calculating them. The matching solution is still applied to these areas in the output. Defaults to None for no mask.
        debug_logs (bool, optional): If True, prints progress. Defaults to False.
        custom_nodata_value (float | int | None, optional): Overrides detected NoData value. Defaults to None.
        image_parallel_workers (Tuple[Literal["process", "thread"], Literal["cpu"] | int] | None = None): Parallelization strategy at the image level. Provide a tuple like ("process", "cpu") to use multiprocessing with all available cores. Threads are supported too. Set to None to disable.
        window_parallel_workers (Tuple[Literal["process"], Literal["cpu"] | int] | None = None): Parallelization strategy at the window level within each image. Same format as image_parallel_workers. Threads are not supported. Set to None to disable.
        window_size (int | Tuple[int, int] | Literal["block"] | None): Tile size for processing: int for square tiles, (width, height) for custom size, or "block" to set as the size of the block map, None for full image. Defaults to None.
        save_as_cog (bool, optional): If True, saves as COG. Defaults to False.
        number_of_blocks (int | tuple | Literal["coefficient_of_variation"]): int as a target of blocks per image, tuple to set manually set total blocks width and height, coefficient_of_variation to find the number of blocks based on this metric.
        alpha (float, optional): Blending factor between reference and local means. Defaults to 1.0.
        correction_method (Literal["gamma", "linear"], optional): Local correction method. Defaults to "gamma".
        save_block_maps (tuple(str, str) | None): If enabled, saves block maps for review, to resume processing later, or to add additional images to the reference map.
            - First str is the path to save the global block map.
            - Second str is the path to save the local block maps, which must include "$" which will be replaced my the image name (because there are multiple local maps).
        load_block_maps (Tuple[str, List[str]] | Tuple[str, None] | Tuple[None, List[str]] | None, optional):
            Controls loading of precomputed block maps. Can be one of:
                - Tuple[str, List[str]]: Load both reference and local block maps.
                - Tuple[str, None]: Load only the reference block map.
                - Tuple[None, List[str]]: Load only the local block maps.
                - None: Do not load any block maps.
            This supports partial or full reuse of precomputed block maps:
                - Local block maps will still be computed for each input image that is not linked to a local block map by the images name being *included* in the local block maps name (file name).
                - The reference block map will only be calculated (mean of all local blocks) if not set.
                - The reference map defines the reference block statistics and the local maps define per-image local block statistics.
                - Both reference and local maps must have the same canvas extent and dimensions which will be used to set those values.
        override_bounds_canvas_coords (Tuple[float, float, float, float] | None): Manually set (min_x, min_y, max_x, max_y) bounds to override the computed/loaded canvas extent. If you wish to have a larger extent than the current images, you can manually set this, along with setting a fixed number of blocks, to anticipate images will expand beyond the current extent.
        block_valid_pixel_threshold (float): Minimum fraction of valid pixels required to include a block (0–1).

    Returns:
        List[str]: Paths to the locally adjusted output raster images.
    """

    print("Start local block adjustment")

    # Validate params
    Universal.validate(
        input_images=input_images,
        output_images=output_images,
        save_as_cog=save_as_cog,
        debug_logs=debug_logs,
        vector_mask=vector_mask,
        window_size=window_size,
        custom_nodata_value=custom_nodata_value,
        image_parallel_workers=image_parallel_workers,
        window_parallel_workers=window_parallel_workers,
        calculation_dtype=calculation_dtype,
        output_dtype=output_dtype,
    )

    Match.validate_local_block_adjustment(
        number_of_blocks=number_of_blocks,
        alpha=alpha,
        correction_method=correction_method,
        save_block_maps=save_block_maps,
        load_block_maps=load_block_maps,
        override_bounds_canvas_coords=override_bounds_canvas_coords,
        block_valid_pixel_threshold=block_valid_pixel_threshold,
    )

    # Determine multiprocessing and worker count
    image_parallel, image_backend, image_max_workers = _resolve_parallel_config(
        image_parallel_workers
    )
    window_parallel, window_backend, window_max_workers = _resolve_parallel_config(
        window_parallel_workers
    )

    input_image_paths = _resolve_paths(
        "search", input_images, kwargs={"default_file_pattern": "*.tif"}
    )
    output_image_paths = _resolve_paths(
        "create",
        output_images,
        kwargs={
            "paths_or_bases": input_image_paths,
            "default_file_pattern": "$_Local.tif",
        },
    )

    if debug_logs:
        print(f"Input images: {input_image_paths}")
    if debug_logs:
        print(f"Output images: {output_image_paths}")

    input_image_names = [
        os.path.splitext(os.path.basename(p))[0] for p in input_image_paths
    ]
    input_image_path_pairs = dict(zip(input_image_names, input_image_paths))
    output_image_path_pairs = dict(zip(input_image_names, output_image_paths))

    _check_raster_requirements(
        input_image_paths,
        debug_logs,
        check_geotransform=True,
        check_crs=True,
        check_bands=True,
        check_nodata=True,
    )

    if isinstance(window_size, int):
        window_size = (window_size, window_size)
    nodata_val = _get_nodata_value(input_image_paths, custom_nodata_value)
    projection = rasterio.open(input_image_paths[0]).crs
    if debug_logs:
        print(f"Global nodata value: {nodata_val}")
    with rasterio.open(input_image_paths[0]) as ds:
        num_bands = ds.count

    # Load data from precomputed block maps if set
    if load_block_maps:
        (
            loaded_block_local_means,
            loaded_block_reference_mean,
            loaded_num_row,
            loaded_num_col,
            loaded_bounds_canvas_coords,
        ) = _get_pre_computed_block_maps(load_block_maps, calculation_dtype, debug_logs)
        loaded_names = list(loaded_block_local_means.keys())
        block_reference_mean = loaded_block_reference_mean

        matched = list(
            (
                soft_matches := {
                    input_name: loaded_name
                    for input_name in input_image_names
                    for loaded_name in loaded_names
                    if input_name in loaded_name
                }
            ).keys()
        )
        only_loaded = [
            l for l in loaded_names if not any(n in l for n in input_image_names)
        ]
        only_input = [
            n for n in input_image_names if not any(n in l for l in loaded_names)
        ]

    else:
        only_input = input_image_names
        matched = []
        only_loaded = []
        block_reference_mean = None

    if debug_logs:
        print(
            f"Total images: input images: {len(input_image_names)}, loaded local block maps: {len(loaded_names) if load_block_maps else 0}:"
        )
        print(
            f"    Matched local block maps (to override) ({len(matched)}):",
            sorted(matched),
        )
        print(
            f"    Only in loaded local block maps (to use) ({len(only_loaded)}):",
            sorted(only_loaded),
        )
        print(
            f"    Only in input (to compute) ({len(only_input)}):", sorted(only_input)
        )

    # Unpack path to save block maps
    if save_block_maps:
        reference_map_path, local_map_path = save_block_maps

    # Create image bounds dict
    bounds_images_coords = {
        name: rasterio.open(path).bounds
        for name, path in input_image_path_pairs.items()
    }

    # Get bounds canvas coords
    if not override_bounds_canvas_coords:
        if not load_block_maps:
            bounds_canvas_coords = _get_bounding_rectangle(input_image_paths)
        else:
            bounds_canvas_coords = loaded_bounds_canvas_coords
    else:
        bounds_canvas_coords = override_bounds_canvas_coords
        if load_block_maps:
            if bounds_canvas_coords != loaded_bounds_canvas_coords:
                raise ValueError(
                    "Override bounds canvas coordinates do not match loaded block maps bounds"
                )

    # Calculate the number of blocks
    if not load_block_maps:
        if isinstance(number_of_blocks, int):
            num_row, num_col = _compute_block_size(
                input_image_paths, number_of_blocks, bounds_canvas_coords
            )
        elif isinstance(number_of_blocks, tuple):
            num_row, num_col = number_of_blocks
        elif isinstance(number_of_blocks, str):
            num_row, num_col = _compute_mosaic_coefficient_of_variation(
                input_image_paths, nodata_val, debug_logs
            )  # This is the approach from the paper to compute bock size
    else:
        num_row, num_col = loaded_num_row, loaded_num_col

    if debug_logs:
        print("Computing local block maps:")

    # Compute local blocks
    local_blocks_to_calculate = {
        k: v for k, v in input_image_path_pairs.items() if k in only_input
    }
    local_blocks_to_load = {
        **{k: loaded_block_local_means[soft_matches[k]] for k in matched},
        **{k: loaded_block_local_means[k] for k in only_loaded},
    }

    if local_blocks_to_calculate:
        args = [
            (
                name,
                path,
                bounds_canvas_coords,
                num_row,
                num_col,
                num_bands,
                window_size,
                debug_logs,
                nodata_val,
                calculation_dtype,
                vector_mask,
                block_valid_pixel_threshold,
                window_parallel,
                window_backend,
                window_max_workers,
            )
            for name, path in local_blocks_to_calculate.items()
        ]

        if image_parallel:
            with _get_executor(image_backend, image_max_workers) as executor:
                futures = [
                    executor.submit(_calculate_block_process_image, *arg)
                    for arg in args
                ]
                results = [f.result() for f in futures]
        else:
            results = [_calculate_block_process_image(*arg) for arg in args]

        block_local_means = {name: mean for name, mean, _ in results}

        overlap = set(block_local_means) & set(local_blocks_to_load)
        if overlap:
            raise ValueError(
                f"Duplicate keys when merging loaded and computed blocks: {overlap}"
            )

        block_local_means = {**block_local_means, **local_blocks_to_load}
    else:
        block_local_means = local_blocks_to_load

    bounds_images_block_space = _get_bounding_rect_images_block_space(block_local_means)

    # Compute reference block
    if debug_logs:
        print("Computing reference block map")
    if block_reference_mean is None:
        block_reference_mean = _compute_reference_blocks(
            block_local_means,
            calculation_dtype,
        )

    if save_block_maps:
        _download_block_map(
            (
                np.nan_to_num(block_reference_mean, nan=nodata_val)
                if nodata_val is not None
                else block_reference_mean
            ),
            bounds_canvas_coords,
            reference_map_path,
            projection,
            calculation_dtype,
            nodata_val,
            num_col,
            num_row,
        )
        for name, block_local_mean in block_local_means.items():
            _download_block_map(
                (
                    np.nan_to_num(block_local_mean, nan=nodata_val)
                    if nodata_val is not None
                    else block_local_mean
                ),
                bounds_canvas_coords,
                local_map_path.replace("$", name),
                projection,
                calculation_dtype,
                nodata_val,
                num_col,
                num_row,
            )
            # _download_block_map(
            #     np.nan_to_num(block_local_count, nan=nodata_val),
            #     bounds_canvas_coords,
            #     os.path.join(output_image_folder, "BlockLocalCount", f"{input_image_name}_BlockLocalCount.tif"),
            #     projection,
            #     calculation_dtype,
            #     nodata_val,
            #     num_col,
            #     num_row,
            # )

    # block_local_mean = _smooth_array(block_local_mean, nodata_value=global_nodata_value)

    # Apply adjustments to images
    if debug_logs:
        print(f"Computing local correction, applying, and saving:")
    args = [
        (
            name,
            input_image_path_pairs[name],
            output_image_path_pairs[name],
            num_bands,
            block_reference_mean,
            block_local_means[name],
            bounds_images_block_space[name],
            bounds_canvas_coords,
            window_size,
            num_row,
            num_col,
            nodata_val,
            alpha,
            correction_method,
            calculation_dtype,
            output_dtype,
            debug_logs,
            window_parallel,
            window_backend,
            window_max_workers,
            save_as_cog,
        )
        for name in input_image_path_pairs
    ]

    if image_parallel:
        with _get_executor(image_backend, image_max_workers) as executor:
            futures = [
                executor.submit(_apply_adjustment_process_image, *arg) for arg in args
            ]
            for future in as_completed(futures):
                future.result()
    else:
        for arg in args:
            _apply_adjustment_process_image(*arg)

    return output_image_paths