16.32. image_align

The program image_align aligns a second image to a first image. In the produced aligned second image, each feature has the same row and column coordinates as in the first image. It can return a transform in pixel space and one in the planet’s ECEF coordinate system.

16.32.1. Overview

Several alignment transforms are supported, including rigid, translation, similarity, etc. The alignment transform is determined with subpixel precision and is applied using bilinear interpolation.

Features are matched among the images using either interest points or a disparity produced with ASP’s correlation algorithms.

If the first image is georeferenced, the second aligned image will use the same georeference as the first one. The first image and second aligned image can then be blended with dem_mosaic (Section 16.20).

Since the first image is kept fixed, if portions of the second aligned image move higher or more to the left than the upper-left corner of the first image, those extra portions will be cut. In that case it is suggested to reverse the order of images when invoking this tool.

The alignment transform can be saved, and a custom alignment transform can be applied instead of the one found automatically. The interest point matches which determine the alignment transform can be saved as well.

This tool extends the co-registration functionality of CASP-GO (Section 15.1).

16.32.2. Assumptions

The images should be at a similar enough resolution and projection (if georeferenced). Otherwise, regrid the second image onto the grid of the first (the reference), for example with gdalwarp and cubic spline interpolation.

If the level of detail or illumination in the two images is too different, the alignment will fail.

When the images are notably different, an approach based on dense pyramid correlation is recommended (Section 16.32.4).

Use the option --individually-normalize if the range of pixel values in the input images differs.

The images are expected to have a single band and have float or integer values. If the images have more than one band, only the first one will be read. The processing is done in double precision. The default output pixel value is float32, as casting to integer may result in precision loss.

16.32.3. Interest point-based alignment

In this mode, features found based on interest points are matched.

image_align                           \
  --alignment-transform rigid         \
  --ip-detect-method 0                \
  --ip-per-image 20000                \
  --individually-normalize            \
  image1.tif image2.tif               \
  --output-prefix out_image_align/run \
  -o image2_align.tif

ASP supports 3 interest point detection methods: OBALoG (0), SIFT (1), and ORB (2) (Section 17.1).

The directory out_image_align will contain the interest point matches (that are cached for future runs), the computed transform, and other auxiliary data.

This program can read externally produced match files in binary or plain-text format (Section 19.10), if they already exist and respect the naming convention (Section 19.11).

16.32.4. Disparity-based alignment

Alternatively, instead of using interest points for alignment, use a (dense) disparity produced from correlation (Section 16.17). This method can be more robust to differences in illumination.

For this path the two images must be on the same grid (same projection, extent, and ground sample distance). Otherwise parallel_stereo --correlator-mode stops with an error. Regrid them first if they differ (see the Assumptions above).

Example:

parallel_stereo              \
  --correlator-mode          \
  --stereo-algorithm asp_mgm \
  --cost-mode 4              \
  --subpixel-mode 9          \
  image1.tif image2.tif      \
  run/run-corr

The disparity will be computed from the first to second image, but the alignment transform is from the second to first image, so the disparity and this transform will show opposite trends.

Run image alignment with the produced filtered disparity run/run-corr-F.tif (Section 19.5):

image_align                                       \
  image1.tif image2.tif                           \
  --disparity-params "run/run-corr-F.tif 1000000" \
  --output-prefix run/run                         \
  --output-image image2_align.tif

For very precise subpixel alignment, use --subpixel-mode 2 above, but this is very slow. See Section 6.1 for the choices when it comes to stereo algorithms and subpixel methods, and Section 16.17 for the image correlator functionality.

When the images are very different, the choice of --cost-mode to use with the asp_mgm algorithm can be quite important. The ternary census cost (--cost-mode 4) has performed better than normalized cross-correlation in our experiments.

For noisy images the asp_bm algorithm should also be considered. It has a larger correlation window size.

See Section 16.32.8 for attaching a per-match uncertainty to the dense matches.

16.32.5. Application for alignment of DEMs

Given a DEM, it can be treated as an image with float values, or an image can be produced from it with the hillshade command (Section 16.28), or, if the DEM is obtained from a stereo point cloud with point2dem, an orthoimage in one-to-one correspondence with this DEM can be created with the --orthoimage option of this tool (Section 16.56).

In either case, given two DEMs, dem1.tif and dem2.tif, their corresponding images can be aligned to each other as:

image_align image1.tif image2.tif --output-prefix run \
  --alignment-transform rigid -o image2_align.tif

Then, the alignment transform can be used to align the second DEM to the first, as:

image_align dem1.tif dem2.tif             \
  --input-transform run/run-transform.txt \
  --output-prefix run/run                 \
  -o dem2_align.tif

It appears that applying this tool on the DEMs themselves may result in more accurate results than if applied on their hillshaded images. (Consider also using for hillshading the tool gdaldem hillshade, Section 16.25.)

If the DEMs have very different grids and projections, regridding them with gdalwarp may make them more similar and easier to align (invoke this tool with cubic spline interpolation).

Note that the alignment transform is a 3x3 matrix and can be examined and edited. Its inputs and outputs are 2D pixels in homogeneous coordinates, that is, of the form (x, y, 1). It is able to model affine and homography transforms in the pixel plane.

See the related tool pc_align (Section 16.53) for alignment of point clouds. That one is likely to perform better than image_align, as it makes use of the 3D nature of point clouds, the inputs need not be gridded, and one of the clouds can be sparse.

16.32.6. Determination of ECEF transform

If the images are georeferenced, this program can find the approximate 3D transform around the planet that brings the second image in alignment with the first one. It is assumed that there exist DEMs associated with these images, from which the 3D coordinates of the locations of interest point matches are determined, and the best-fit 3D transform is computed.

Example:

image_align img1.tif img2.tif \
  -o img2_align.tif           \
  --alignment-transform rigid \
  --ecef-transform-type rigid \
  --dem1 dem1.tif             \
  --dem2 dem2.tif             \
  --output-prefix run/run

This will save run/run-ecef-transform.txt in the pc_align format (rotation + translation + scale, Section 16.53.5). This transform can be passed to pc_align in order to transform a point cloud (Section 16.53.6), and to bundle_adjust if desired to transform cameras (Section 16.53.14). Use zero iterations with these tools to apply the transform without further refinement.

It is important to keep in mind that the ECEF transform is from the second cloud to the first, hence pc_align should have the clouds in the same order as for image_align in order to use this transform.

The inverse of this transform is saved as well, if desired to transform the clouds or cameras from the coordinate system of the first image to the one of the second image.

If no DEMs exist, the images themselves can be used in their place. The grayscale values will be interpreted as heights above the datum in meters. The image_calc program (Section 16.33) can modify these values before the DEMs are passed to image_align.

If only DEMs exist, their hillshaded versions (Section 16.28) can be used as images. As earlier, the more similar visually the images are, the better the results.

It is suggested to use --alignment-transform rigid and --ecef-transform-type rigid if it is thought that a rotational component exists, and the value translation for these options if no rotation is expected.

Note that this will produce a rotation + translation around planet center, rather than a local “in-plane” transform, so it can be hard to interpret. A similarity transform can be used when there is a difference in scale.

Note that this transform is an approximation. It is not possible to precisely convert a 2D transform between images to a 3D transform in ECEF unless the underlying terrain is perfectly flat.

16.32.7. Evaluating matches between orthoimages

This program can evaluate the geolocation agreement between two georeferenced orthoimages, such as a reference basemap and a source image. In this mode the actual alignment is skipped. The interest point matches are found and saved to a GeoPackage, recording where each matched feature lands in each image, in the units of the georeference in the input images (meters or degrees), and the offset between them.

image_align reference.tif source.tif         \
  --output-prefix run/run                    \
  --match-points-geopackage run/matches.gpkg \
  --individually-normalize                   \
  --ip-per-image 0 --ip-per-tile 200         \
  --matches-per-tile 100

Here, the output image (-o) is not set, so no alignment transform is applied. If -o is also set, the aligned image is produced as well.

The matching uses the same interest point engine as for alignment, so the options --individually-normalize (helpful for images of differing modality or wavelength), --ip-per-tile, --matches-per-tile, and the RANSAC filtering controls (--num-ransac-iterations, --inlier-threshold) all apply. Both images must be georeferenced.

The matches are filtered with RANSAC. The option --inlier-threshold (default: 50 pixels) sets how far a match may deviate from the best-fit transform and still be kept as an inlier. When the two images are already approximately co-registered, a smaller value (a few pixels) is recommended, as it rejects spurious matches caused by illumination or content differences between the images. A larger value should be used when a substantial misregistration between the images is expected.

For images of differing modality or wavelength, finding matches from a dense disparity (Section 16.32.4) rather than sparse interest points is often more robust, and it can attach a per-match uncertainty to each GeoPackage feature (Section 16.32.8).

16.32.8. Per-match uncertainty

To create dense matches with an uncertainty, while producing a GeoPackage report as in Section 16.32.7, run the parallel_stereo disparity step from Section 16.32.4 with the extra option:

--save-left-right-disparity-difference

This writes run/run-corr-L-R-disp-diff.tif, an integer per-pixel reliability measure (Section 14.4.3). Then pass that file as the third entry of --disparity-params:

image_align                          \
  reference.tif source.tif           \
  --disparity-params                 \
    "run/run-corr-F.tif 1000000
     run/run-corr-L-R-disp-diff.tif" \
  --output-prefix run/run            \
  --match-points-geopackage run/matches.gpkg

The sigma for each match in the GeoPackage (Section 16.32.8.1) is set to the corresponding left-right disparity difference, with a lower bound of 0.5 pixels (so it is never 0). A larger value means a less self-consistent (so less reliable) match.

16.32.8.1. GeoPackage format

The output is a GeoPackage (.gpkg) with one point feature per inlier match. The point geometry is the source-image location, expressed in the projection of the reference image. When the two images are in different projections, the source coordinates are reprojected into the reference one, so the dx, dy offset is always a valid difference within a single coordinate system. All coordinates and offsets are in the units of the reference projection (meters or degrees). Each feature has the following fields:

• ref_x, ref_y: the coordinates of the match in the reference image projection.

• src_x, src_y: the source-image match location, in the reference projection as well.

• dx, dy: the offset, as source minus reference. This is the local geolocation discrepancy between the two images.

• ref_col, ref_row, src_col, src_row: the pixel locations of the match in each image.

• sigma: the per-match sigma, in pixels. For disparity-based matches given a left-right disparity difference (Section 16.32.8), this is that difference, with a minimum of 0.5 pixels: a per-match uncertainty, larger where the match is less self-consistent. Otherwise it is the interest point detector scale, a feature size rather than a co-registration uncertainty.

The GeoPackage can be inspected or converted with ogrinfo and ogr2ogr (Section 16.25), or opened in QGIS.

The same interest point matches are also written as a match file under the output prefix (binary, or plain text with --matches-as-txt); see Section 19.10 for that format.

16.32.9. Usage

image_align [options] <reference image> <source image> \
  --output-prefix <prefix> -o <aligned source image>

16.32.10. Command-line options for image_align

--output-image, -o <string (default: “”)>

Specify the output image.

--output-prefix <string (default: “out_image_align/run”)>

Save the interest point matches, computed transform, and other auxiliary data at this prefix. These are cached for future runs.

--alignment-transform <string (default: “rigid”)>

Specify the transform to use to align the second image to the first. Options: translation, rigid (translation + rotation), similarity (translation + rotation + scale), affine, homography.

--output-data-type, -d <string (default: “float32”)>

The data type of the output file. Options: uint8, uint16, uint32, int16, int32, float32, float64. The values are clamped (and also rounded for integer types) to avoid overflow.

--ip-detect-method <integer (default: 0)>

Interest point detection algorithm (0: Integral OBALoG (default), 1: OpenCV SIFT, 2: OpenCV ORB).

--ip-per-image <integer (default: 20000)>

How many interest points to detect in each image (the resulting number of matches will be much less).

--ip-per-tile <integer (default: 0)>

How many interest points to detect in each 1024^2 image tile (default: automatic determination). This is before matching. Not all interest points will have a match. See also --matches-per-tile.

--matches-per-tile <integer (default: 0)>

How many interest point matches to compute in each image tile (of size normally 1024^2 pixels). Use a value of --ip-per-tile a few times larger than this. See also --matches-per-tile-params.

--matches-per-tile-params <string (default: “1024 1280”)>

To be used with --matches-per-tile. The first value is the image tile size for both images. A larger second value allows each right tile to further expand to this size, resulting in the tiles overlapping. This may be needed if the homography alignment between these images is not great, as this transform is used to pair up left and right image tiles.

--individually-normalize

Individually normalize the input images instead of using common values.

--matches-as-txt

Read and write match files as plain text instead of binary. See Section 19.10.

--num-ransac-iterations <integer (default: 1000)>

How many iterations to perform in RANSAC when finding interest point matches.

--inlier-threshold <double (default: 50.0)>

The inlier threshold (in pixels) to separate inliers from outliers when computing interest point matches. A smaller threshold will result in fewer inliers.

--min-matches <integer (default: 10)>

Set the minimum number of inlier matches between images for successful matching.

--disparity-params <string (default: “”)>

Find the alignment transform by using, instead of interest points, a disparity, such as produced by parallel_stereo --correlator-mode. Specify as a string in quotes, in the format: “disparity.tif num_samples”. An optional third entry, the left-right disparity difference (from parallel_stereo --save-left-right-disparity-difference), sets each match’s sigma to that value, with a lower bound of 0.5 pixels (Section 16.32.8).

--input-transform <string (default: “”)>

Instead of computing an alignment transform, read and apply the one from this file. Must be stored as a 3x3 matrix.

--ecef-transform-type <string (default: “”)>

Save the ECEF transform corresponding to the image alignment transform to <output prefix>-ecef-transform.txt. The type can be: ‘translation’, ‘rigid’ (rotation + translation), or ‘similarity’ (rotation + translation + scale). See Section 16.32.6 for an example.

--dem1 <string (default: “”)>

The DEM associated with the first image. To be used with --ecef-transform-type.

--dem2 <string (default: “”)>

The DEM associated with the second image. To be used with --ecef-transform-type.

--nodata-value <float (default: NaN)>

Pixels with values less than or equal to this number are treated as no-data. This overrides the no-data values from input images.

--match-points-geopackage <string (default: “”)>

Write the inlier interest point matches to this GeoPackage (.gpkg) file, in the units of the georeference (meters or degrees). In this mode the actual alignment is skipped, unless -o is set. See Section 16.32.7 for details.

--threads <integer (default: 0)>

Select the number of threads to use for each process. If 0, use the value in ~/.vwrc.

--cache-size-mb <integer (default = 1024)>

Set the system cache size, in MB.

--tile-size <integer (default: 256 256)>

Image tile size used for multi-threaded processing.

--no-bigtiff

Tell GDAL to not create BigTiff files.

--tif-compress <None|LZW|Deflate|Packbits (default: LZW)>

TIFF compression method.

-v, --version

Display the version of software.

-h, --help

Display this help message.