16.49. pc_align

This tool can be used to align two point clouds. The algorithms employed are one of the several flavors of Iterative Closest Point (ICP), based on the libpointmatcher library [PCSM13]:

It also implements the Fast Global Registration algorithm from:

In addition, it supports feature-based alignment (terrains are hillshaded and interest point matches are found among them), and alignment using least squares. It can handle a scale change in addition to rotations and translations. For joint alignment of more than two clouds, the related tool n_align can be used (Section 16.40).

Usage:

pc_align --max-displacement <float> [other options]    \
  <reference cloud> <source cloud> -o <output prefix>}

The denser cloud must be the first one to be passed to this tool. This program is very sensitive to the value of --max-displacement (Section 16.49.2).

An example of using this tool is in Section 6.1.12. Evaluation of alignment is discussed in Section 16.49.9 and Section 16.49.10.

See the related tool image_align (Section 16.30) for performing alignment of images.

Several important things need to be kept in mind if pc_align is to be used successfully and give accurate results, as described below.

16.49.1. The input point clouds

Due to the nature of ICP, the first input point cloud, that is, the reference (fixed) cloud, should be denser than the second, source (movable) point cloud, to get the most accurate results. This is not a serious restriction, as one can perform the alignment this way and then simply invert the obtained transform if desired (pc_align outputs both the direct and inverse transform, and can output the reference point cloud transformed to match the source and vice-versa).

The user can choose how many points to pick from the reference and source point clouds to perform the alignment. The amount of memory and processing time used by pc_align is directly proportional to these numbers, ideally the more points the better. Pre-cropping to judiciously chosen regions may improve the accuracy and/or run-time.

16.49.2. The max displacement option

In many typical applications, the source and reference point clouds are already roughly aligned, but the source point cloud may cover a larger area than the reference. The user should provide to pc_align the expected maximum distance (displacement) source points may move by as result of alignment, using the option --max-displacement.

This number will help remove source points too far from the reference point cloud which may not match successfully and may degrade the accuracy. If in doubt, this value can be set to something large but still reasonable, as the tool is able to throw away a certain number of unmatched outliers.

At the end of alignment, pc_align will display the observed maximum displacement, a multiple of which can be used to seed the tool in a subsequent run. If an initial transform is applied to the source cloud (Section 16.49.6), the outliers are thrown out after this operation. The observed maximum displacement is also between the source points with this transform applied and the source points after alignment to the reference.

16.49.3. Alignment method

The default alignment method is Point-to-Plane ICP, which may be more robust to large translations than Point-to-Point ICP, though the latter can be good enough if the input point clouds have small alignment errors and it is faster and uses less memory as well. The tool also accepts an option named --highest-accuracy which will compute the normals for Point-to-Plane ICP at all points rather than about a tenth of them. This option is not necessary most of the time, but may result in better alignment at the expense of using more memory and processing time.

The default alignment transform is rigid, that is, a combination of rotation and translation. With Point-to-Point ICP, it is also possible to solve for a scale change (to obtain a so-called similarity transform). It is suggested this approach be used only when a scale change is expected. It can be turned on by setting --alignment-method to similarity-point-to-plane or similarity-point-to-point. (The first of these is better than the second one.)

For very large scale difference or translation among the two clouds, both of these algorithms may fail. If the clouds are DEMs, one may specify the option --initial-transform-from-hillshading string which will hillshade the two DEMs, find interest point matches among them, and use that to compute an initial transform between the clouds (Section 16.49.6), which may or may not contain scale, after which the earlier algorithms will be applied to refine the transform. See an example in Section 8.24.9.

This functionality is implemented with ASP’s hillshade, ipfind, and ipmatch tools, and pc_align has options to pass flags to these programs, such as to increase the number interest points being found, if the defaults are not sufficient. If the two clouds look too different for interest point matching to work, they perhaps can be re-gridded to use the same (coarser) grid, as described in Section 16.49.13, to obtain the initial transform which can then be applied to the original clouds.

A non-ICP algorithm supported by ASP is Fast Global Registration, accessible with --alignment-method fgr, and customizable using the --fgr-options field (see the table below for more details). This approach can perform better than ICP when the clouds are close enough to each other but there is a large number of outliers, since it does a cross-check, so it can function with very large --max-displacement. It does worse if the clouds need a big shift to align.

This one is being advertised as less sensitive to outliers, hence it should give good results with a larger value of the maximum displacement.

Another option is to use least squares (with outlier handling using a robust cost function) to find the transform, if the reference cloud is a DEM. For this, one should specify the alignment method as least-squares or similarity-least-squares (the latter also solves for scale). It is suggested that the input clouds be very close or otherwise the --initial-transform option be used, for the method to converge, and use perhaps on the order of 10-20 iterations and a smaller value for --max-num-source-points (perhaps a few thousand) for this approach to converge reasonably fast.

16.49.4. File formats

The input point clouds can be in one of several formats: ASP’s point cloud format (the output of stereo), DEMs as GeoTIFF or ISIS cub files, LAS files, or plain-text CSV files (with .csv or .txt extension).

By default, CSV files are expected to have on each line the latitude and longitude (in degrees), and the height above the datum (in meters), separated by commas or spaces. Alternatively, the user can specify the format of the CSV file via the --csv-format option. Entries in the CSV file can then be (in any order) (a) longitude, latitude (in degrees), height above datum (in meters), (b) longitude, latitude, distance from planet center (in meters or km), (c) easting, northing and height above datum (in meters), in this case a PROJ.4 string must be set via --csv-proj4, (d) Cartesian coordinates \((x, y, z)\) measured from planet center (in meters). The precise syntax is described in the table below. The tool can also auto-detect the LOLA RDR PointPerRow format.

Any line in a CSV file starting with the pound character (#) is ignored.

If none of the input files have a geoheader with datum information, and the input files are not in Cartesian coordinates, the datum needs to be specified via the --datum option, or by setting --semi-major-axis and --semi-minor-axis.

16.49.5. The alignment transform

The transform obtained by pc_align is output to a text file as a 4 × 4 matrix with the upper-left 3 × 3 submatrix being the rotation (and potentially also a scale, per Section 16.49.3) and the top three elements of the right-most column being the translation. It is named <output prefix>-transform.txt.

This transform, if applied to the source point cloud, will bring it in alignment with the reference point cloud. The transform assumes the 3D Cartesian coordinate system with the origin at the planet center (known as ECEF). This matrix can be supplied back to the tool as an initial guess (Section 16.49.6).

The inverse transform, from the reference cloud to the source cloud is saved as well, as <output prefix>-inverse-transform.txt.

These two transforms can be used to move cameras from one cloud’s coordinate system to another one’s, as shown in Section 16.49.14.

16.49.6. Applying an initial transform

The transform output by pc_align can be supplied back to the tool as an initial guess via the --initial-transform option, with the same clouds as earlier, or some supersets or subsets of them. If it is desired to simply apply this transform without further work, one can specify --num-iterations 0.

This may be useful, for example, in first finding the alignment transform over a smaller, more reliable region (e.g., over rock, excluding moving ice), then applying it over the entire available dataset. To illustrate this, consider a DEM, named dem.tif, obtained with ASP, from whom just a portion, dem_crop.tif is known to have reliable measurements, which are stored, for example, in a file called meas.csv. Hence, pc_align is first used on the smaller DEM, as:

pc_align <other options> dem_crop.tif meas.csv -o run/run

Then, the command:

pc_align --max-displacement -1 --num-iterations 0 \
  --save-transformed-source-points                \
  --save-inv-transformed-reference-points         \
  --initial-transform run/run-transform.txt       \
  --csv-format <csv format string>                \
  dem.tif meas.csv -o run_full/run

will transform the full dem.tif into the coordinate system of meas.csv, and meas.csv into the coordinate system of ref.tif with no further iterations. See also Section 16.49.14 for how to use such transforms with cameras.

If an initial transform is used, with zero or more iterations, the output transform produced by such an invocation will be from the source points before the initial transform, hence the output alignment transform will incorporate the initial transform.

Using --max-displacement -1 should be avoided, as that will do no outlier filtering in the source cloud. Here that is not necessary, as this invocation simply moves the DEM according to the specified transform.

If a good initial alignment is found, it is suggested to use a smaller value for --max-displacement to refine the alignment, as the clouds will already be mostly on top of each other after the initial transform is applied.

16.49.7. Applying an initial specified translation or rotation

One can apply to the source cloud an initial shift, expressed in the North-East-Down coordinate system at the centroid of the source points, before the alignment algorithm is invoked. Hence, if it is desired to first move the source cloud North by 5 m, East by 10 m, and down by 15 m relative to the point on planet surface which is the centroid of the source points, the continue with alignment, one can invoke pc_align with:

--initial-ned-translation "5 10 15"

(Notice the quotes.)

The option --initial-rotation-angle can be used analogously.

As in Section 16.49.6, one can simply stop after such an operation, if using zero iterations. In either case, such initial transform will be incorporated into the transform file output by pc_align, hence that one will go from the source cloud before user’s initial transform to the reference cloud.

16.49.8. Interpreting the transform

The alignment transform, with its origin at the center of the planet, can result in large movements on the planet surface even for small angles of rotation. Because of this it may be difficult to interpret both its rotation and translation components.

The pc_align program outputs the translation component of this transform, defined as the vector from the centroid of the original source points (before any initial transform applied to them) to the centroid of the source points with the computed alignment transform applied to them. This translation component is displayed in three ways (a) Cartesian coordinates with the origin at the planet center, (b) Local North-East-Down coordinates at the centroid of the source points (before any initial transform), and (c) Latitude-Longitude-Height differences between the two centroids. If the effect of the transform is small (e.g., the points moved by at most several hundred meters) then the representation in the form (b) above is most amenable to interpretation as it is in respect to cardinal directions and height above ground if standing at a point on the planet surface.

This program prints to screen the Euler angles of the rotation transform, and also the axis of rotation and the angle measured against that axis. It can be convenient to interpret the rotation as being around the center of gravity of the reference cloud, even though it was computed as a rotation around the planet center, since changing the point around which a rigid transform is applied will only affect its translation component, which is relative to that point, but not the rotation matrix.

16.49.9. Error metrics and outliers

The tool outputs to CSV files the lists of errors together with their locations in the source point cloud, before the alignment of the source points (but after applying any initial transform), and also after the alignment computed by the tool. They are named <output prefix>-beg_errors.csv and <output prefix>-end_errors.csv. An error is defined as the distance from a source point used in alignment to the closest reference point (measured in meters).

The format of output CSV files is the same as of input CSV files, or as given by --csv-format, although any columns of extraneous data in the input files are not saved on output. The first line in these files shows the names of the columns.

See Section 16.64.6 for how to visualize these files. By default, this tool shows the 4th column in these files, which is the absolute error difference. Run, for example:

stereo_gui --colorbar run/run-end_errors.csv

The program prints to screen and saves to a log file the 16th, 50th, and 84th error percentiles as well as the means of the smallest 25%, 50%, 75%, and 100% of the errors.

When the reference point cloud is a DEM, a more accurate computation of the errors from source points to the reference cloud is used. A source point is projected onto the datum of the reference DEM, its longitude and latitude are found, then the DEM height at that position is interpolated. That way we determine the closest point on the reference DEM that interprets the DEM not just as a collection of points but rather as a polyhedral surface going through those points. These errors are what is printed in the statistics. To instead compute errors as done for other type of point clouds, use the option --no-dem-distances.

By default, when pc_align discards outliers during the computation of the alignment transform, it keeps the 75% of the points with the smallest errors. As such, a way of judging the effectiveness of the tool is to look at the mean of the smallest 75% of the errors before and after alignment.

16.49.10. Evaluation of aligned clouds

The pc_align program can save the source cloud after being aligned to the reference cloud and vice-versa, via --save-transformed-source-points and --save-inv-transformed-reference-points.

To validate that the aligned source cloud is very close to the reference cloud, DEMs can be made out of them with point2dem (Section 16.52), and those can be overlaid as georeferenced images in stereo_gui (Section 16.64) for inspection. A GIS tool can be used as well.

Alternatively, the geodiff program (Section 16.23) can be used to compute the (absolute) difference between aligned DEMs, which can be colorized with colormap (Section 16.14), or colorized on-the-fly and displayed with a colorbar in stereo_gui (Section 16.64.5).

The geodiff tool can take the difference between a DEM and a CSV file as well. The obtained error differences can be visualized in stereo_gui (Section 16.64.6).

16.49.11. Output point clouds and convergence history

The transformed input point clouds (the source transformed to match the reference, and the reference transformed to match the source) can also be saved to disk if desired. If an input point cloud is in CSV, ASP point cloud format, or LAS format, the output transformed cloud will be in the same format. If the input is a DEM, the output will be an ASP point cloud, since a gridded point cloud may not stay so after a 3D transform. The point2dem program can be used to re-grid the obtained point cloud back to a DEM.

As an example, assume that pc_align is run as:

pc_align --max-displacement 100              \
  --csv-format '1:x 2:y 3:z'                 \
  --save-transformed-source-points           \
  --save-inv-transformed-reference-points    \
  ref_dem.tif source.csv                     \
  -o run/run

This will save run/run-trans_reference.tif which is a point cloud in the coordinate system of the source dataset, and run/run-trans_source.csv which is in reference coordinate system of the reference dataset.

Care is needed, as before, with setting --max-displacement.

The convergence history for pc_align (the translation and rotation change at each iteration) is saved to disk with a name like:

<output prefix>-iterationInfo.csv

and can be used to fine-tune the stopping criteria.

16.49.12. Manual alignment

If automatic alignment fails, for example, if the clouds are too different, or they differ by a scale factor, a manual alignment can be computed as an initial guess transform (and one can stop there if pc_align is invoked with 0 iterations).

For that, the input point clouds should be first converted to DEMs using point2dem, unless in that format already. Then, stereo_gui can be called to create manual point correspondences (interest point matches) from the reference to the source DEM (hence they should be displayed in the GUI in this order, from left to right, and one can hillshade them to see features better).

Once the match file is saved to disk, it can be passed to pc_align via the --match-file option, which will compute an initial transform before continuing with alignment. This transform can also be used for non-DEM clouds once it is found using DEMs obtained from those clouds.

16.49.13. Creating a point cloud from a DEM

Given a DEM, if one invokes pc_align as follows:

pc_align dem.tif dem.tif --max-displacement -1 --num-iterations 0 \
   --save-transformed-source-points -o run/run

this will create a point cloud out of the DEM. This cloud can then be re-gridded using point2dem at a lower resolution or with a different projection.

16.49.14. Applying the pc_align transform to cameras

If pc_align is used to align a DEM obtained with ASP to a preexisting reference DEM, the obtained alignment transform can be applied to the cameras used to create the ASP DEM, so the cameras then become aligned with the pre-existing DEM. That is accomplished by running bundle adjustment with the options --initial-transform and --apply-initial-transform-only.

Please note that the way this transform is applied depends on the order of DEMs in pc_align and on whether the cameras have been bundle-adjusted or not. Precise commands are given below.

First, assume, for example, that the reference DEM is ref.tif, and the ASP DEM is created without bundle adjustment, as:

parallel_stereo left.tif right.tif left.xml right.xml output/run
point2dem output/run-PC.tif

(See further down for when the cameras have been bundle-adjusted.)

It is very important to distinguish the cases when the obtained DEM is the first or second argument of pc_align.

If the ASP DEM output/run-DEM.tif is aligned to the reference DEM as:

pc_align --max-displacement 1000 ref.tif output/run-DEM.tif \
  -o align/run

then, the alignment is applied to cameras the following way:

bundle_adjust left.tif right.tif left.xml right.xml \
  --initial-transform align/run-transform.txt       \
  --apply-initial-transform-only -o ba_align/run

This should create the adjusted cameras incorporating the alignment transform:

ba_align/run-left.adjust, ba_align/run-right.adjust

(see Section 16.5.9 for discussion of .adjust files).

If pc_align was invoked with the two DEMs in reverse order, the transform to use is:

align/run-inverse-transform.txt

The idea here is that run-transform.txt goes from the second DEM passed to pc_align to the first, hence, bundle_adjust invoked with this transform would move cameras from second DEM’s coordinate system’s to first. And vice-versa, if run-inverse-transform.txt is used, cameras from first DEM’s coordinate system would be moved to second’s.

After applying a transform this way, cameras which are now aligned with the reference DEM can be used to mapproject onto it, hopefully with no registration error as:

mapproject ref.tif left.tif left_map.tif \
  --bundle-adjust-prefix ba_align/run

and in the same way for the right image.

If, the initial stereo was done with cameras that already were bundle-adjusted, with output prefix initial_ba/run, so the stereo command had the option:

--bundle-adjust-prefix initial_ba/run

we need to integrate those initial adjustments with this alignment transform. To do that, again need to consider two cases, as before.

If the just-created stereo DEM is the second argument to pc_align, then run the slightly modified command:

bundle_adjust left.tif right.tif left.xml right.xml \
  --initial-transform align/run-transform.txt       \
  --input-adjustments-prefix initial_ba/run         \
  --apply-initial-transform-only -o ba_align/run

Otherwise, if the stereo DEM is the first argument to pc_align, use instead align/run-inverse-transform.txt as input to --initial-transform.

Note that this way bundle adjustment will not do any further camera refinements after the initial transform is applied.

A stereo run can be reused after the cameras have been modified as above, with the option --prev-run-prefix. Only triangulation will then be redone. Ensure the option --bundle-adjust-prefix ba_align/run is used to point to the new cameras. See Section 8.27.11 and Section 6.1.7.8.

16.49.15. Troubleshooting

Remember that filtering is applied only to the source point cloud. If you have an input cloud with a lot of noise, make sure it is being used as the source cloud.

If you are not getting good results with pc_align, something that you can try is to convert an input point cloud into a smoothed DEM. Use point2dem to do this and set --search-radius-factor if needed to fill in holes in the DEM. For some input data this can significantly improve alignment accuracy.

16.49.16. Command-line options for pc_align

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

Maximum number of iterations.

--max-displacement <float>

Maximum expected displacement (horizonal + vertical) of source points as result of alignment, in meters (after the initial guess transform is applied to the source points). Used for removing gross outliers in the source (movable) point cloud.

-o, --output-prefix <filename>

Specify the output file prefix.

--outlier-ratio <float (default: 0.75)>

Fraction of source (movable) points considered inliers (after gross outliers further than max-displacement from reference points are removed).

--max-num-reference-points <integer (default: 10^8)>

Maximum number of (randomly picked) reference points to use.

--max-num-source-points <integer (default: 10^5)>

Maximum number of (randomly picked) source points to use (after discarding gross outliers).

--alignment-method <string (default: point-to-plane)>

The type of iterative closest point method to use. Choices: point-to-plane, point-to-point, similarity-point-to-plane, similarity-point-to-point, fgr, least-squares, similarity-least-squares.

--highest-accuracy

Compute with highest accuracy for point-to-plane (can be much slower).

--datum <string>

Sets the datum for CSV files. Options:

  • WGS_1984

  • D_MOON (1,737,400 meters)

  • D_MARS (3,396,190 meters)

  • MOLA (3,396,000 meters)

  • NAD83

  • WGS72

  • NAD27

  • Earth (alias for WGS_1984)

  • Mars (alias for D_MARS)

  • Moon (alias for D_MOON)

--semi-major-axis <float>

Explicitly set the datum semi-major axis in meters.

--semi-minor-axis <float>

Explicitly set the datum semi-minor axis in meters.

--csv-format <string>

Specify the format of input CSV files as a list of entries column_index:column_type (indices start from 1). Examples: 1:x 2:y 3:z (a Cartesian coordinate system with origin at planet center is assumed, with the units being in meters), 5:lon 6:lat 7:radius_m (longitude and latitude are in degrees, the radius is measured in meters from planet center), 3:lat 2:lon 1:height_above_datum, 1:easting 2:northing 3:height_above_datum (need to set --csv-proj4; the height above datum is in meters). Can also use radius_km for column_type, when it is again measured from planet center.

--csv-proj4 <string>

The PROJ.4 string to use to interpret the entries in input CSV files, if those files contain Easting and Northing fields.

--compute-translation-only

Compute the transform from source to reference point cloud as a translation only (no rotation).

--save-transformed-source-points

Apply the obtained transform to the source points so they match the reference points and save them.

--save-inv-transformed-reference-points

Apply the inverse of the obtained transform to the reference points so they match the source points and save them.

--initial-transform <string>

The file containing the transform to be used as an initial guess. It can come from a previous run of the tool.

--initial-ned-translation <string>

Initialize the alignment transform based on a translation with this vector in the North-East-Down coordinate system around the centroid of the reference points. Specify it in quotes, separated by spaces or commas.

--initial-rotation-angle <double (default: 0.0)>

Initialize the alignment transform as the rotation with this angle (in degrees) around the axis going from the planet center to the centroid of the point cloud. If --initial-ned-translation is also specified, the translation gets applied after the rotation.

--initial-transform-from-hillshading <string>

If both input clouds are DEMs, find interest point matches among their hillshaded versions, and use them to compute an initial transform to apply to the source cloud before proceeding with alignment. Specify here the type of transform, as one of: ‘similarity’ (rotation + translation + scale), ‘rigid’ (rotation + translation) or ‘translation’. See the options further down for tuning this.

--hillshade-options

Options to pass to the hillshade program when computing the transform from hillshading. Default: --azimuth 300 --elevation 20 --align-to-georef.

--ipfind-options

Options to pass to the ipfind program when computing the transform from hillshading. Default: --ip-per-image 1000000 --interest-operator sift --descriptor-generator sift.

--ipmatch-options

Options to pass to the ipmatch program when computing the transform from hillshading. Default: --inlier-threshold 100 --ransac-iterations 10000 --ransac-constraint similarity.

--initial-transform-ransac-params <num_iter factor (default: 10000 1.0)>

When computing an initial transform based on hillshading, use this number of RANSAC iterations and outlier factor. A smaller factor will reject more outliers.

--match-file

Compute an initial transform from the source to the reference point cloud using manually selected point correspondences (obtained for example using stereo_gui). The type of transform can be set via --initial-transform-from-hillshading string. It may be desired to change --initial-transform-ransac-params if it rejects as outliers some manual matches.

--fgr-options

Options to pass to the Fast Global Registration algorithm, if used. Default: div_factor: 1.4 use_absolute_scale: 0 max_corr_dist: 0.025 iteration_number: 100 tuple_scale: 0.95 tuple_max_cnt: 10000.

--diff-rotation-error <float (default: 1e-8)>

Change in rotation amount below which the algorithm will stop (if translation error is also below bound), in degrees.

--diff-translation-error <float (default: 1e-3)>

Change in translation amount below which the algorithm will stop (if rotation error is also below bound), in meters.

--no-dem-distances

For reference point clouds that are DEMs, don’t take advantage of the fact that it is possible to interpolate into this DEM when finding the closest distance to it from a point in the source cloud (the text above has more detailed information).

--skip-shared-box-estimation

Do not estimate the shared bounding box of the two clouds. This estimation can be costly for large clouds but helps with eliminating outliers.

--config-file <file.yaml>

This is an advanced option. Read the alignment parameters from a configuration file, in the format expected by libpointmatcher, over-riding the command-line options.

--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 bigtiffs.

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

TIFF compression method.

-v, --version

Display the version of software.

-h, --help

Display this help message.