Table of Contents
 See Also
 opals::IICP
Aim of module
Performs orientation (coregistration) of multiple point cloud datasets using the Iterative Closest Point algorithm.
Requirements for usage of opalsICP
This module requires the installation of the Matlab runtime.
General description
The Iterative Closest Point (ICP) algorithm (Chen & Medioni (1991), Besl & McKay (1992)) improves the alignment of two (or more) point clouds by minimizing iteratively the discrepancies within the overlap area of these point clouds. Nowadays the term ICP does not necessarily refer to the algorithm presented in the original publications, but rather to a group of surface matching algorithms which have in common the following aspects:

I
: correspondences are established iteratively (i.e. not only once) 
C
: as correspondence the closest point, or more generally, the corresponding point, is used 
P
: correspondences are established on a point basis (instead of using e.g. interpolated surfaces or rasters)
The following flow diagram illustrates the functionality of this module. The main iteration loop consists of 7 main steps, which are described below.
In the simplest case, the ICP algorithm is used to improve the orientation of a single loose point cloud with respect to a single fixed point cloud. The main steps of the algorithm are explained on the basis of this case. The initial orientation of the two point clouds is visualized in Figure 2. As can be seen, the point clouds are already roughly aligned, which is a main requirement of the ICP algorithm.
Step 1: Find overlap: The overlap area of the point clouds is determined by using socalled voxel hulls. Voxel hulls are a low resolution representation of the volume occupied by a point cloud. For the computation of the voxel hull the object space is subdivided into a voxel structure. The voxel hull of a point cloud consists of all voxels which contain at least one point of the point cloud. The parameter voxelSize defines the edge length of a single voxel.
Step 2: Selection: A subset of points is selected within the overlap area in one point cloud. For this, points are first selected uniformly in object space, which leads to a homogeneous distribution of the points within the overlap area (uniform sampling). The mean sampling distance in each coordinate direction is defined by the parameter samplingDist. The selected subset of points can be further refined with the strategies normal space sampling or maximum leverage sampling, which are described in section Advanced selection strategies.
Step 3: Matching: Find the closest points (= nearest neighbors) of the selected subset in the other point cloud. This step leads to a set of correspondences which possibly also includes some outliers.
Step 4: Rejection: Rejection of false correspondences (outliers) based on the compatibility of points. Correspondences are rejected if:
 the euclidean distance between corresponding points is larger than the value specified by the parameter maxDist.
 the roughness attribute of one of the corresponding points is larger than the value specified by the parameter maxSigma. (Note: as roughness attribute the standard deviation in normal direction of the points used for plane fitting is used.)
 the angle between the normals of corresponding points is larger than the value specified by the parameter maxAngleDev.
Since it is not guaranteed that with these strategies all outliers are rejected, a robust adjustment method is used in step 6 (Minimization) for the detection and deactivation of the remaining ones.
Step 5: Weighting: In this step a weight ( ) is assigned to each correspondence. Weights are based on:
 the roughness attribute of corresponding points
 the angle between the normals of corresponding points.
Step 6: Minimization: Estimation of the transformation parameters (for the loose point cloud) by a least squares adjustment. The adjustment minimizes the sum of squared pointtoplane distances. The pointtoplane distance of two corresponding points is defined as the orthogonal distance of one point to the fitted plane of the other point.
Step 7: Transformation: Transformation of the loose point cloud with the estimated parameters. The transformation model can be chosen with the parameter trafoType.
After step 7, a convergence criterion is tested. If it is not met, the process restarts with a new iteration from step 1 (Figure 1). However, if subsets are defined by the parameter subsetRadius, the iteration loop restarts from step 3 and the first two steps are only performed once. The total number of iterations heavily depends on the initial orientation of the points clouds, but typically, between 3 and 10 iterations should be sufficient to reach the global minimum.
Note that opalsICP is not restricted to two point clouds, but can handle an arbitrary number of point clouds simultaneously.
Advanced selection strategies
The selection of points (step 2) can heavily affect the final alignment accuracy. Therefore, opalsICP offers two advanced selection strategies, which are especially useful for point clouds where one normal direction is predominating, but the data still includes some valuable features for the alignment. First, points are always selected with the uniform sampling strategy described above, where the mean sampling distance is defined by the parameter samplingDist. To select a subset of this selection, one of the two following strategies can be applied:
 normal space sampling: The aim of this strategy is to select points such that the distribution of their normals in angular space is as uniform as possible. The percentage of points selected by this strategy can be defined with the parameter normalSpaceSampling.
 maximum leverage sampling: This strategy selects those points, which are best suited for the estimation of the transformation parameters. For this, the leverage of each point on the parameter estimation is considered. The percentage of points selected by this strategy can be defined with the parameter maxLeverageSampling.
A comparison of the selection strategies offered by opalsICP is given in Figure 8. For most ICP variants, this ALS scene is rather difficult because only one feature  the ditch  can constrain the transformation at the finest level. Since normal space sampling and maximum leverage sampling consider the usefulness of points for the alignment process, the number of points can be dramatically reduced, without affecting the solubility of the adjustment. Further information can be found in Glira et al. 2015.
Parameter description
Note that the default parameter values were chosen to be well suited for ALS point clouds. For point clouds at other scales, the Examples section may serve as a reference.
Remarks: mandatory
The input files containing the unregistered point clouds (i.e. individual scan positions, flight strips...) are expected either OPALS Data Manager (ODM) or in binary XYZ (coordinates: 64bit double precision) format. Note that wildcard characters (*,?) can be used for the selection of multiple files at once (e.g. scan*.odm).
Remarks: estimable
The reduction point defines the origin of a local coordinate system in which the point clouds are stored internally. Using a reduction point is crucial for large coordinates, since otherwise transformation parameters are highly correlated.
Estimation rule: If not explicitly specified, the center of gravity of the bounding boxes of all involved input ODM datasets are used.
Remarks: optional
A filter string in EBNF syntax as described in section 'Filter syntax' can be passed to restrict the set of input points used within the ICP algorithm (e.g. only last echoes or only smooth areas). Note that this parameter does not affect the final ODM files in which (regardless of the filter settings) always all points are included.
Remarks: mandatory
To specify the folder where the final transformed ODM files and intermediate files are stored. Note that the transformed ODM files are only generated if transformData is activated.
Remarks: default=tempICP
To specify the folder where temporary data are stored (e.g. intermediate MATLAB files).
Remarks: default=full
Possible values: zShift .... compute zshift only shifts .... compute 3Dshifts only rigid ..... compute rigid 3D transformation (3 rotations + 3 shifts = 6 parameters) helmert ... compute 3D Helmert transformation (3 rotations + 3 shifts + 1 scale = 7 parameters) full ...... compute full 3D affine transformation (12 parameters)
Allows to specify the transformation model applied to the point clouds within the ICP algorithm.
Remarks: estimable
In order to define the datum of the adjustment, at least one point cloud has to be fixed in object space. With this parameter the file names of one or more fixed point clouds can be defined.
Estimation rule: If not specified, the first inFile is fixed.
Remarks: default=10
The edge length of a single voxel used for the voxel hull computation to detect overlaps. 100 must be a whole multiple of it: fmod(100,x) == 0.
Remarks: default=10
The selection of correspondences is based on a 3d uniform sampling strategy. This parameter defines the mean sampling distance in each coordinate direction. In the case that more than one value is passed, specific sampling distances are used for the n individual input datasets (inFile). If less than n values are specified the last value is used for the remaining datasets.
Remarks: default=2
The implemented ICP algorithm minimizes pointtoplane distances i.e. points in one point cloud are compared to planes in all overlapping point clouds. All points within a sphere of the specified radius are considered for the plane fitting. A good choice of the search radius is based on the point cloud density and the geometry of the scanned object. Note: to ensure a certain redundancy, planes are only considered if the search area contains at least 8 points.
Remarks: default=5
Maximum number of iterations within the ICP algorithm.
Remarks: default=0
For large point clouds (e.g. ALS flight block), it is usually not possible to load all point clouds in memory simultaneously. Thus, only a small subset of points is selected for each point cloud. For this points are selected around the established correspondences within a specific radius defined by this parameter. If subsets are used, the main iterations loop restarts from step 3 instead of step 1 (Figure 1), i.e. step 1 and step 2 are only performed once. This parameter can be set to zero if no subsets should be used, i.e. all point clouds can be kept in memory simultaneously.
Remarks: estimable
To reject false correspondences, a maximum allowed (pointtopoint) euclidean distance between two corresponding points can be defined with this parameter.
Estimation rule: 2*searchRadius.
Remarks: default=0.1
To reject false correspondences a maximum allowed roughness value (i.e. standard deviation of plane fitting) can be defined with this parameter. This way the ICP algorithm only considers smooth surfaces for estimating the transformation parameters. Whenever possible, a small threshold should be used to ensure highquality correspondences. Higher values, in turn, increase the number of correspondences.
Remarks: default=10
To reject false correspondences, a maximum deviation between the angles of normals of two corresponding points can be defined with this parameter.
Remarks: optional
Normal space sampling is an advanced selection strategy which can be used to reduce the points previously selected by the parameter samplingDist. This parameter defines the percentage of the remaining points.
Remarks: optional
Maximum leverage sampling is an advanced selection strategy which can be used to reduce the points previously selected by the parameter samplingDist. This parameter defines the percentage of the remaining points.
Remarks: default=0
If specified, after each ICP iteration a 3d point cloud visualization and the most relevant ICP results are displayed in a window.
Remarks: default=1
If activated, the estimated transformations are applied to each individual input data set. The resulting ODM files are stored in the folder specified by the outDirectory parameter (subfolder post).
Remarks: output
The results of the ICP computation are the transformation parameters for each input dataset.
Examples
The data used in the following examples can be found in the $OPALS_ROOT/demo/
directory.
ALS strip adjustment
In this example a small strip adjustment (3 strips) is calculated with opalsICP. As a first step, we have to import the strips, which are given in las format, using Module Import :
Now we can run opalsICP. Besides the mandatory parameters inFile and outDirectory, we also set the samplingDist to 5m:
Please note, that all values displayed in this visualization are also included in the command line output:
The following files are created in the folder specified by the parameter outDirectory:

Final (not filtered) transformed point clouds in odm format: version01\G111.odmversion01\G112.odmversion01\G113.odm

Command line output (also contained in
opalsLog.xml
):version01\ICPLog.txt 
ICP result file for further analysis in Matlab: version01\ICP.mat
The following files are created in the temporary folder (not explicitely specified in the above command):

(Filtered) input point clouds in mat format: These point cloud files correspond to the bxyz files and can be opened in Matlab (as objects of class pointCloud).tempICP\pre\G111.mattempICP\pre\G112.mattempICP\pre\G113.mat
Finally, we can export the new odm files to the las format using Module Export or Module Translate :
Close range TLS dataset
In this example the orientation of six overlapping close range TLS scans is improved by opalsICP. The initial orientation of the scans is shown on the lefthand side of Figure 10 (coordinate units = mm). First, we have to import the las files with Module Import :
Now we call opalsICP. But notice: since the default parameter values of this module were chosen to be well suited for ALS point clouds, compared to the previous example, a few more parameters have to be specified:
The new odm files are placed in the folder version01
.
TLS dataset
The orientation of three TLS scans of a room is improved in this example. First, let's import the data:
The following call of opalsICP differs from the previous examples in two points:
 A wildcard character is used to define the list of input files.

The input data is not transformed within opalsICP, which means that no odm files are created in the folder
version01\post
. Instead, we save the transformation parameters into the fileicp_output.xml
. This file can afterwards be used to directly transform the input data into any format supported by opals (e.g. the las format), without creating intermediate odm files.
Now, we can directly create the transformed las files with Module Export , using Parameter Mapping to set its parameter trafo:
References
Besl P., McKay N., 1992. A Method for Registration of 3D Shapes. In: IEEE Trans. PAMI, Vol. 14, No. 2: 239256.
Chen Y., Medioni, G., 1991. Object Modeling by Registration of Multiple Range Images. In: Proc. IEEE Conf. on Robotics and Automation, Sacramento, California: 27242729.
Glira P., Pfeifer N., Briese C., Ressl C., 2015: A correspondence framework for ALS strip adjustments based on variants of the ICP algorithm.
 Date
 23.2.2015