KLT-based Algorithm for Registration of Images from Observing Systems (KARIOS)
In general, quality assessment process are fundamental to appreciate how well data fit for Earth Observation application purposes. Those assessment dedicated to geometric accuracy are rarely opened to community. As consequence, it is always difficult to inter compare data based on the same process and it is always difficult to compare results based on harmonized mapping accuracy metrics.
To overcome this situation, thanks to funding of ESA / EDAP EDAP project, the KARIOS initiative has been promoted and a user tool is now available.
The KARIOS tool has been designed to analyse geometric deformations within optical / radar images. For this purpose, the tool performs image matching and generate several key graphical presentations and compute accuracy statistics.
Image matching process does not follow traditional approach because it is based on feature point matching (corner). A KLT implementation available in OpenCV library is used in KARIOS. Also, the Candidate point selection is done with GoodFeaturesToTrack function and matching is done with calcOpticalFlowPyrLK function.
As show in the following picture, KARIOS makes KLT algorithm compatible with remote sensing images embedding suitable pre-processing (image filtering) / post-processing (outlier filtering).
As an optional and experimental feature, KARIOS have the capability to detect large shift between images. If a large shift is detected, the monitored image is shifted according to the offsets it found, and then apply the KLT matching.
To enable large shift detection, use
--enable-large-shift-detection
program argument.NOTICE that it could use lots of memory in case of large image such as Sentinel2 B04 (11GB)
Furthermore, KARIOS analyses displacements between the two input image grids both in line and pixel direction outputing, providing user with the three following items:
The geometric accuracy report includes the following accuracy metrics, in both directions when relevant:
- Root Mean Square Error
- Minimum / Maximum Error
- Mean Error
- Standard deviation Error
- Circular Error @90 percentile
The Circular Error (CE) at the 90% level confidence graphic is used for horizontal accuracy in image products.
This representation is relevant for image expressed within cartographic system grid.
Because, with the CE representation, it is straightforward to evaluate mapping accuracy, considering reference data with known accuracy.Rather, in case of images with no cartographic system grid, the CE graphic representation becomes less informative.
The CE graphic is still generated, and equally spaced sample data is assumed.
This hypothesis is not obvious, when details on image grids are unknown.
This tool is a Python application, to run it you should have a dedicated conda environnement.
To do so, install conda or miniconda, then run the following command to create the karios conda env
conda env create -f environment.yml
Then activate the env:
conda activate karios
KARIOS takes as mandatory inputs :
- monitored sensor image file
- reference sensor image file
Input files shall contain only one layer of data, and the format shall recognized by gdal library.
NOTICE: Inputs images grids should be comparable. The user should take care of its data preparation. That means geo coded images must have the same footprint, same geo transform information (same EPSG code) and same resolution. Image pixel resolution should also be square (same X,Y) and unit meter.
This is also applicable to DEM for
--dem-file-path
option
Examples:
- Without DEM:
python karios/karios.py \
/data/12SYH/LS9_OLIL2F_20220824T175017_N0403_R035_T12SYH_20230214T164934.SAFE/GRANULE/L2F_T12SYH_A000000_20220824T175017_LS9_R035/IMG_DATA/L2F_T12SYH_20220824T175017_LS9_R035_B04_10m.TIF \
/data/References/GRI/T12SYH_20220514T175909_B04.jp2
- With DEM:
python karios/karios.py \
--dem-file-path /data/12SYH/dem10.tiff \
--dem-description "Copernicus DSM_30 resample from 90m to 10m" \
/data/12SYH/LS9_OLIL2F_20220824T175017_N0403_R035_T12SYH_20230214T164934.SAFE/GRANULE/L2F_T12SYH_A000000_20220824T175017_LS9_R035/IMG_DATA/L2F_T12SYH_20220824T175017_LS9_R035_B04_10m.TIF \
/data/References/GRI/T12SYH_20220514T175909_B04.jp2
- csv file: list of key points and associated dx/dy deviations
- png file: visualisation of the deviations
- geojson of KP found by co registration process
usage: karios.py [-h] [--conf CONF] [--resume] [--mask MASK_FILE_PATH] [--input-pixel-size PIXEL_SIZE] [--out OUT] [--generate-key-points-mask] [--generate-intermediate-product] [--title-prefix TITLE_PREFIX]
[--dem-file-path DEM_FILE_PATH] [--dem-description DEM_DESCRIPTION] [--enable-large-shift-detection] [--no-log-file] [--debug] [--log-file-path LOG_FILE_PATH]
MONITORED_IMAGE_PATH REFERENCE_IMAGE_PATH
options:
-h, --help show this help message and exit
Mandatory arguments:
MONITORED_IMAGE_PATH Path to the monitored sensor product
REFERENCE_IMAGE_PATH Path to the reference sensor product
Processing options:
--conf CONF Configuration file path (default: /opt/karios/karios/configuration/processing_configuration.json)
--resume Do not run KLT matcher, only accuracy analysis and report generation (default: False)
--mask MASK_FILE_PATH
Path to the mask to apply to the reference image (default: None)
--input-pixel-size PIXEL_SIZE, -pxs PIXEL_SIZE
Input image pixel size in meter. Ignored if image resolution can be read from input image (default: None)
Output options:
--out OUT Output results folder path (default: /opt/karios/results)
--generate-key-points-mask, -kpm
Generate a tiff mask based on KP from KTL (default: False)
--generate-intermediate-product, -gip
Generate a two band tiff based on KP with band 1 dx and band 2 dy (default: False)
--title-prefix TITLE_PREFIX, -tp TITLE_PREFIX
Add prefix to title of generated output charts (limited to 26 characters) (default: None)
DEM arguments (optional):
--dem-file-path DEM_FILE_PATH
DEM file path. If given, "shift mean by altitude group plot" is generated. (default: None)
--dem-description DEM_DESCRIPTION
DEM source name. It is added in "shift mean by altitude group plot" DEM source (example: COPERNICUS DEM resample to 10m). Ignored if --dem-file-path is not given (default: None)
Experimental (optional):
--enable-large-shift-detection
If enabled, KARIOS looks for large pixel shift between reference and monitored image. When a significant shift is detected, KARIOS shifts the monitored image according to the offsets it computes
and then process to the matching (default: False)
Logging arguments (optional):
--no-log-file Do not log in file (default: False)
--debug, -d Enable Debug mode (default: False)
--log-file-path LOG_FILE_PATH
Log file path (default: /opt/karios/karios.log)
Notice that default folder and file path are adapted to your KARIOS install folder.
The default configuration is located in karios/configuration/processing_configuration.json
processing_configuration.shift_image_processing
parameters (Large Shift Matching processing parameters)
bias_correction_min_threshold
: number of pixel threshold from which large shift is applied.
xStart
: image X margin to apply (margin is skipped by the matcher)tile_size
: tile size to process by KTL in the input imagelaplacian_kernel_size
: Aperture size used to compute the second-derivative filters of Laplacian process
The following parameter allows to control how to find the most prominent corners in the reference image, as described by the OpenCV documentation goodFeaturesToTrack, after applying Laplacian.
minDistance
: Minimum possible Euclidean distance between the returned corners.blocksize
: Size of an average block for computing a derivative covariation matrix over each pixel neighbourhood.maxCorners
: Maximum number of corners to extract. If there are more corners than are found, the strongest of them is returned.maxCorners = 0
implies that no limit on the maximum is set and all detected corners are returned.qualityLevel
: Parameter characterizing the minimal accepted quality of image corners. The parameter value is multiplied by the best corner quality measure. The corners with the quality measure less than the product are rejected. For example, if the best corner has the quality measure = 1500, and the qualityLevel=0.01, then all the corners with the quality measure less than 15 are rejected.matching_winsize
: size of the search window during matching corners in the reference and the monitored Laplacian images.outliers_filtering
: whether to filter or not the outliers points found during the matching.
Refer to section KLT param leverage for details
confidence_threshold
: max score for points found by the matcher to use to compute statistics written in correl_res.txt. IfNone
, not applied.
fig_size
: Size of the generated figure in inchesshift_colormap
: matplotlib color map name for the KP shift error scatter plotshift_auto_axes_limit
: auto compute KP shift error colorbar scaleshift_axes_limit
: KP shift error colorbar maximum limit, N/A ifshift_auto_axes_limit
istrue
theta_colormap
: matplotlib color map name for the KP theta error scatter plot
fig_size
: Size of the generated figure in inchesscatter_colormap
: matplotlib color map name for the KP shift scatter plotscatter_auto_limit
: auto compute KP shift scatter plot limitscatter_min_limit
: KP shift scatter plot minimum limit, N/A ifscatter_auto_limit
istrue
scatter_max_limit
: KP shift scatter plot maximum limit, N/A ifscatter_auto_limit
istrue
histo_mean_bin_size
: KP shift histogram bin size (number of image row/col for the histogram bin)
fig_size
: Size of the generated figure in inchesshow_fliers
: draw fliers of box plothisto_mean_bin_size
: KP altitude histogram bin size (altitude ranges size)
fig_size
: Height size of the generated figure in inches, width is 5/3 of the heightce_scatter_colormap
: matplotlib color map name for the KP shift density scatter plot
In order to have a lower memory usage during KLT process, it is possible to define a tile size to process for KLT.
For example, a tile_size of 10000 for an image having a size of 20000x20000 pixels will result of 4 tiles to process.
In this context, the KLT process will look in each tiles for maxCorners
.
While an image of 20000x20000 pixels results of 4 equals tiles, an image of 20000x15000 pixels also result of 4 tiles, but with different size, two of 10000x10000 pixels and two of 10000x5000 pixels.
The consequence is that the density for matching point will not be the same each tiles, the bigger tiles will have a lower matching point than the smallest.
You may also consider that the image can content empty parts where KLT will not find any matching point. So tiles having a large empty parts will also results to a bigger matching point density.
In order to avoid density difference in the final result, you can define a tile_size
largest than the image with an hight maxCorners
, or a small tile_size
and maxCorners
in order to have tiles with almost same size.
For example, for image of 20000x15000 pixels, you should consider a tile_size
of 20000 (1 tile), or 5000 (12 equal tiles)
This output use box plot to show statistics of KP on altitudes groups.
The box extends from the first quartile (Q1) to the third quartile (Q3) of the data, with a line at the median. The whiskers extend from the box to the farthest data point lying within 1.5x the inter-quartile range (IQR) from the box. Flier points are those past the end of the whiskers. See https://en.wikipedia.org/wiki/Box_plot for reference.
Q1-1.5IQR Q1 median Q3 Q3+1.5IQR
|-----:-----|
o |--------| : |--------| o o
|-----:-----|
flier <-----------> fliers
IQR
credits https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.boxplot.html