 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.
Figure 1: Functionality of opalsICP.
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.
Figure 2: Initial state of the point clouds.
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.
Figure 3: Overlap area of the point clouds (grey filled voxels) as intersection of the two voxel hulls (blue resp. red voxels).
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.
Figure 4: Selection of points in first point cloud.
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.
Figure 5: Matching of the selected points (blue / from previous step) to the closest points (red) in the second point cloud.
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 ( \(0 \le w \le 1\)) 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.
Figure 6: Minimization of the pointtoplane distances (green) between corresponding points.
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.
Figure 7: Final state of the point clouds.
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.
Figure 8: The advanced methods 'normal space sampling' and 'maximum leverage sampling' for the selection of points (step 2).
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.
inFilelist of input files
Type:
opals::Vector<
opals::Path>
Remarks: mandatory
The input files containing the unregistered point clouds (i.e. individual scan positions, flight strips...) can be provided in any OPALS supported Format. If not OPALS Data Manager (ODM) format is used, the reduction point currently has to be stated explicitely using redPoint. Note that wildcard characters (*,?) can be used for the selection of multiple files at once (e.g. scan*.odm).
redPoint3d coordinate reduction point
Type:
opals::Array<double, 3>
Remarks: estimable
The reduction point defines the origin of a local coordinate system in which the point clouds are stored internally  for different input formats than ODM, this parameter is mandatory. 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.
filtermodify the input using a (tree of) filter(s)
Type:
opals::String
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.
outDirectoryoutput directory
Type:
opals::Path
Remarks: mandatory
To specify the folder where the final transformed datset files and additional log and MATLAB files are stored. Note that the transformed datasets are only generated if transformData is activated and that they are furthermore created in the same data format as the input files.
tempDirectorytemporary data directory
Type:
opals::Path
Remarks: default=tempICP
To specify the folder where temporary data are stored (e.g. intermediate MATLAB files).
trafoTypetransformation type
Type:
opals::LSMTrafoType
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.
fixedFilefile name(s) of fixed point cloud(s)
Type:
opals::Vector<
opals::Path>
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.
voxelSizevoxel size of voxel hulls
Type: double
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.
samplingDistmean distance between corresponding points
Type:
opals::Vector<double>
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.
searchRadiussearch radius for plane fitting
Type: double
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.
maxItermaximum number of iterations
Type: int32
Remarks: default=5
Maximum number of iterations within the ICP algorithm.
subsetRadiusradius of selected subset areas
Type: double
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.
maxDistmaximum distance between corresponding points
Type: double
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.
maxSigmamaximum plane roughness
Type: double
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.
maxAngleDevmaximum angle (in degrees) between the normals of corresponding points
Type: double
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.
normalSpaceSamplingpercentage value for normal space sampling
Type: double
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.
maxLeverageSamplingpercentage value for maximum leverage sampling
Type: double
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.
plotvisualization of ICP results
Type: bool
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.
transformDataperform input data transformation
Type: bool
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).
outTrafParsoutput transformation parameters
Type:
opals::Vector<
opals::TrafPars3dAffine>
Remarks: output
The results of the ICP computation are the transformation parameters for each input dataset.
TrafPars3dAffine: Describes a relative or absolute 3Daffine transformation. This class contains the transformation parameters of an affine 3D (i.e. 12parameter) transformation, given as 12array
Examples
The data used in the following examples can be found in the $OPALS_ROOT/demo/
directory.
Example 1: ALS strip adjustment
In this example a small strip adjustment (3 strips) is calculated with opalsICP. Besides the mandatory parameters inFile and outDirectory, we also set the samplingDist to 5m and the reduction point (redPoint) since we are directly using LAS data:
Figure 9: Results obtained with opalsICP (set the parameter plot to true for this visualization).
Please note, that all values displayed in this visualization are also included in the command line output:
B  S : GLOBALICP > RUNICP > ICP ITERATION 5 OF 5 > RESULTS
B  T : 
B  T : GENERAL STATISTICS:
B  T : iteration corresp. std(dp) mean(dp) norm(dx)
B  T : 1 2422 0.04827 0.00189 6.42463
B  T : 2 2369 0.03621 0.00070 5.78050
B  T : 3 2363 0.03632 0.00073 0.94349
B  T : 4 2365 0.03605 0.00074 0.29409
B  T : new: 5 2362 0.03611 0.00075 0.11907
B  T : where dp = vector of all (signed) pointtoplane distances
B  T : dx = vector of parameter increments
B  T : 
B  T : POINT CLOUD DEPENDANT STATISTICS:
B  T : point_cloud corresp. std(dp) mean(dp) file
B  T : [1] 1121 0.03689 0.00034 G111
B  T : [2] 2059 0.03478 0.00087 G112
B  T : [3] 1544 0.03726 0.00090 G113
B  T : where dp = vector of all (signed) pointtoplane distances associated to a single point cloud
B  T : 
B  E : GLOBALICP > RUNICP > ICP ITERATION 5 OF 5 > RESULTS
The following files are created in the folder specified by the parameter outDirectory:
The following files are created in the temporary folder (not explicitly specified in the above command):
Please note that the transformed output files (icpALS\*.las
) are exported in the same data format (LAS version, point data record format) as the input LAS files. Please use Module Import to convert the LAS files to OPALS Datamanager format for further use within OPALS.
Example 2: 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). Please note: 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:
Figure 10: The orientation before and after running opalsICP.
The new OPALS Datamanager files are placed in the folder icpCloseRange
.
Example 3: TLS dataset
The orientation of three TLS scans of a room is improved in this example. As a first step, the data are imported to ODM format:
opalsImport inFile roomScan1.laz outFile roomScan1.odm
opalsImport inFile roomScan2.laz outFile roomScan2.odm
opalsImport inFile roomScan3.laz outFile roomScan3.odm
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
icpTLS
. Instead, the transformation parameters are written to the file icp_output.xml
. This file can afterwards be used to directly transform the input data into any format supported by OPALS.
Figure 11: A crosssection through the scanned interior before and after running opalsICP.
Now, we can directly create the transformed LAS files with Module Export using Parameter Mapping to set the transformation parameters trafo:
opalsExport inFile roomScan1.odm outFile icpTLS\roomScan1.las 
inParamFiles icp_output.xml 
paramMapping "trafo = opalsICP.outTrafPars me.inFile[0].stem().startswith( Path(you.getIdGridMov()).stem() )"
opalsExport inFile roomScan2.odm outFile icpTLS\roomScan2.las 
inParamFiles icp_output.xml 
paramMapping "trafo = opalsICP.outTrafPars me.inFile[0].stem().startswith( Path(you.getIdGridMov()).stem() )"
opalsExport inFile roomScan3.odm outFile icpTLS\roomScan3.las 
inParamFiles icp_output.xml 
paramMapping "trafo = opalsICP.outTrafPars me.inFile[0].stem().startswith( Path(you.getIdGridMov()).stem() )"
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.
 Author
 pg, gm, jo, mw
 Date
 11.11.2018