Module ICP
See Also
opals::IICP

Aim of module

Performs orientation (co-registration) 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.

steps_diagram_small.png
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.

initial_state_small.png
Figure 2: Initial state of the point clouds.

Step 1: Find overlap: The overlap area of the point clouds is determined by using so-called 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.

voxel_hulls_small.png
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.

selection_small.png
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.

matching_small.png
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 point-to-plane distances. The point-to-plane distance of two corresponding points is defined as the orthogonal distance of one point to the fitted plane of the other point.

minimization_small.png
Figure 6: Minimization of the point-to-plane 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.

final_state_small.png
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.

selection_advanced_small.png
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 (odm|bxyz)
Type: opals::Vector<opals::Path>
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: 64-bit double precision) format. 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. 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 ODM files and intermediate files are stored. Note that the transformed ODM files are only generated if transformData is activated.
-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 z-shift only
  shifts .... compute 3D-shifts 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 point-to-plane 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: int
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 (point-to-point) 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 high-quality 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.

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 :

opalsImport -inFile G111.las -outFile G111.odm
opalsImport -inFile G112.las -outFile G112.odm
opalsImport -inFile G113.las -outFile G113.odm

Now we can run opalsICP. Besides the mandatory parameters inFile and outDirectory, we also set the samplingDist to 5m:

opalsICP -inFile G111.odm G112.odm G113.odm -outDirectory version01 -samplingDist 5
exampleALS_small.png
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:

S : GLOBALICP > RUNICP > ICP ITERATION 5 OF 5 > RESULTS
T : GENERAL STATISTICS:
T : iteration corresp. std(dp) mean(dp) norm(dx)
T : 1 1215 0.03805 -0.00025 0.49176
T : 2 1198 0.03108 0.00021 0.07041
T : 3 1191 0.03057 0.00016 0.03387
T : 4 1192 0.03056 0.00022 0.01185
T : 5 1193 0.03062 0.00028 0.00697
T : where ... dp is the vector of all (signed) point-to-plane distances
T : ... dx is the vector of parameter increments
T : POINT CLOUD DEPENDANT STATISTICS:
T : point cloud corresp. std(dp) mean(dp) file
T : [1] 530 0.03004 0.00034 G111.mat
T : [2] 1063 0.03039 0.00061 G112.mat
T : [3] 793 0.03133 -0.00021 G113.mat
T : where ... dp is the vector of all (signed) point-to-plane distances associated to a single point cloud
E : GLOBALICP > RUNICP > ICP ITERATION 5 OF 5 > RESULTS

The following files are created in the folder specified by the parameter outDirectory:

  • Final (not filtered) transformed point clouds in odm format:
    version01\G111.odm
    version01\G112.odm
    version01\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:
    tempICP\pre\G111.mat
    tempICP\pre\G112.mat
    tempICP\pre\G113.mat
    These point cloud files correspond to the bxyz files and can be opened in Matlab (as objects of class pointCloud).

Finally, we can export the new odm files to the las format using Module Export or Module Translate :

opalsExport -inFile version01\G111.odm -outFile version01\G111.las
opalsExport -inFile version01\G112.odm -outFile version01\G112.las
opalsExport -inFile version01\G113.odm -outFile version01\G113.las

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 left-hand side of Figure 10 (coordinate units = mm). First, we have to import the las files with Module Import :

opalsImport -inFile lionScan1.las -outFile lionScan1.odm
opalsImport -inFile lionScan2.las -outFile lionScan2.odm
opalsImport -inFile lionScan3.las -outFile lionScan3.odm
opalsImport -inFile lionScan4.las -outFile lionScan4.odm
opalsImport -inFile lionScan5.las -outFile lionScan5.odm
opalsImport -inFile lionScan6.las -outFile lionScan6.odm

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:

opalsICP -inFile lionScan1.odm lionScan2.odm lionScan3.odm lionScan4.odm lionScan5.odm lionScan6.odm -outDirectory version01 -samplingDist 4 -searchRadius 2 -maxSigma 0.5 -trafoType rigid
exampleTLSCR_small.png
Figure 10: The orientation before and after running opalsICP.

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:

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 version01\post. Instead, we save the transformation parameters into the file icp_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.
opalsICP -inFile roomScan?.odm -trafoType rigid -voxelSize 0.2 -samplingDist 0.5 -searchRadius 0.1 -maxSigma 0.02 -outDirectory version01 -transformData false -outParamFile icp_output.xml
exampleTLS_small.png
Figure 11: A cross-section 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 its parameter trafo:

opalsExport -inFile roomScan1.odm -outFile version01\roomScan1.las -inParamFiles icp_output.xml -paramMapping "trafo = opalsICP.outTrafPars me.inFile[0].stem().startswith( Path(you.getIdGridMov()).stem() )"
opalsExport -inFile roomScan2.odm -outFile version01\roomScan2.las -inParamFiles icp_output.xml -paramMapping "trafo = opalsICP.outTrafPars me.inFile[0].stem().startswith( Path(you.getIdGridMov()).stem() )"
opalsExport -inFile roomScan3.odm -outFile version01\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 3-D Shapes. In: IEEE Trans. PAMI, Vol. 14, No. 2: 239-256.

Chen Y., Medioni, G., 1991. Object Modeling by Registration of Multiple Range Images. In: Proc. IEEE Conf. on Robotics and Automation, Sacramento, California: 2724-2729.

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
Date
23.2.2015