Module DirectGeoref

Table of Contents

See Also
opals::IDirectGeoref

Aim of module

Transforms points from scanner coordinates into geo-referenced coordinates by applying the sensor model (flight trajectory, mounting parameters).

General description

For direct georeferencing of ALS data, the following input is required:

  • ALS data given in the laser scanner's coordinate system
  • mounting calibration
  • flight trajectory
  • a coordinate reference system (coord_ref_sys)

As output, we obtain the ALS ground points with respect to the defined spatial reference system (e.g. UTM33 North or geocentric WGS84), commonly termed as "world system".

Mounting calibration

We suppose that we have a strapdown inertial navigation system (INS), i.e. the inertial platform is fixed with respect to the aircraft's body. Consequently, its measurements are given in a body(-fixed) system $B$. Its axes are defined as follows: $x^B$: forward, $y^B$: right, $z^B$: down. Furthermore, we suppose that its origin corresponds to the current (interpolated) position of the flight trajectory.

Since the laser scanner data are usually not directly given in the body system $B$, the so-called mounting calibration is required, which contains the transformation chain required to describe the transition from the scanner's coordinate system $S$ to the body system $B$. Figure 1 shows these two coordinate systems and two other intermediate systems $S_0$ and $M$.

mountingCalibration_100dpi.png
Fig. 1: Mounting calibration

$S_0$ is the scanner's coordinate system in zero position in case if a scanner tilting device is used (if not, the systems $S_0$ and $S$ are identical). In general, the system $S$ is rotated (tilt rotation) and shifted (tilt shift) with respect to the system $S_0$.

The mounting system $M$ is defined as follows: Its axes result from cyclic permutation (scanner system) of those of system $S_0$ in such a way that they are approximately parallel to those of the body system $B$. The systems $M$ and $S_0$ have the same origin. The mount system is (usually slightly) rotated (mounting rotation) and shifted (mounting shift) with respect to the body system.

Using matrix notation, the relations between these coordinate systems may be described as follows:

$\mathbf{x}^B=\mathbf{R}_M^B\cdot\mathbf{x}^M + t_M^B$

$ \mathbf{x}^{M}=\mathbf{R}_{S_0}^M\cdot\mathbf{x}^{S_0}$

$ \mathbf{x}^{S_0}=\mathbf{R}_S^{S_0}\cdot\mathbf{x}^S + t_S^{S_0} $

The direct relation between the systems $S_0$ and $B$ is:

$ \mathbf{x}^B=\mathbf{R}_M^B\cdot\mathbf{R}_{S_0}^M\cdot\mathbf{x}^{S_0}+t_M^B=\mathbf{R}_{S_0}^B\cdot\mathbf{x}^{S_0}+t_{S_0}^B $

General notes on notation:

  • The columns of rotation matrix $\mathbf{R}_J^K$ correspond to the unit vectors of the axes of system J expressed in system K
  • The translation (shift) vector $\mathbf{t}_J^K$ describes the origin of system J with respect to system K

Apart from the above-mentioned geometric transformation elements, the mounting calibration contains the time offset $\Delta t$ (time lag) of the trajectory time system with respect to the laserscanner time system, too.

Specifying the parameters of the mounting calibration

The parameters of the mounting calibration may be specified using a string according to the following EBNF (Extended Backus-Naur Form) syntax:

<mounting pars> =
[ <time lag> ],
[ <scanner system> ],
[ <mounting rotation> ],
[ <mounting shift> ],
[ <tilt rotation> ],
[ <tilt shift> ];
(* may be specified in arbitrary order *)
<time lag> = "TIMELAG", "(", <real>, ")";
<scanner system> = "SCANNERSYS", "(", <dir sequence>, ")";
<mounting rotation> = "MOUNTROTATION", [<refsys>], "(", <rotation>, ")";
<mounting shift> = "MOUNTSHIFT", [<refsys>], "(", <shift>, ")";
<tilt rotation> = "TILTROTATION", [<refsys>], "(", <rotation>, ")";
<tilt shift> = "TILTSHIFT", [<refsys>], "(", <shift>, ")";
(* default values: *)
(* "TIMELAG(0.0)" *)
(* "SCANNERSYS(F-R-D)" *)
(* "MOUNTROTATION=GLOBAL(MATRIX(1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0))" *)
(* "MOUNTSHIFT=GLOBAL(0.0 0.0 0.0)" *)
(* "TILTROTATION=GLOBAL(MATRIX(1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0))" *)
(* "TILTSHIFT=GLOBAL(0.0 0.0 0.0)" *)
<dir sequence> = "F-R-D" | "F-D-L" | "F-L-U" | "F-U-R" |
"B-D-R" | "B-L-D" | "B-U-L" | "B-R-U" |
"L-F-D" | "L-D-B" | "L-B-U" | "L-U-F" |
"R-D-F" | "R-B-D" | "R-U-B" | "R-F-U" |
"U-F-L" | "U-L-B" | "U-B-R" | "U-R-F" |
"D-L-F" | "D-B-L" | "D-R-B" | "D-F-R";
(* F...FRONT; B...BACK; L...LEFT; R...RIGHT; U...UP; D...DOWN *)
(* all 24 possible right-handed systems are allowed *)
<rotation> = <vectors> | <matrix> |
( <angles>, [<axis hierarchy>], [<rotation sense>], [<units>] );
(* default values in case of <angles>: *)
(* "AXISHIERARCHY(X-Y-Z), SENSEOFROT(CCW), UNITS(DEG)" *)
<shift> = <dx>, <dy>, <dz>;
<refsys> = "=", "GLOBAL" | "LOCAL"; (* default: "=GLOBAL" *)
<vectors> = "VECTORS", "(", <vector sequence>, ")";
<matrix> = "MATRIX", "(", <matrix elements>, ")";
<angles> = "ANGLES", "(", <alpha>, <beta>, <gamma>, ")";
<axis hierarchy> = "AXISHIERARCHY", "(", <axis sequence>, ")";
<axis sequence> = "X-Y-Z" | "X-Z-Y" | "Y-Z-X" | "Y-X-Z" | "Z-X-Y" | "Z-Y-X";
<rotation sense> = "SENSEOFROT", "(", ( "CCW" | "CW"), ")";
<units> = "UNITS", "(", ( "DEG", "GRAD", "RAD" ), ")";
<vector sequence> = ( <x axis> | <y axis> | <z axis> );
(* at least two of the three possible elements have to be specified *)
(* the entered vectors need not to be normalised but they must be orthogonal to each other with a certain accuracy *)
<matrix elements> = 9 * <real>;
(* to be specified in column major order, i.e. r11 r21 r31 r12... (e.g. r21 is the element of row 2 and col 1) *)
(* the entered matrix must fullfill orthogonality conditions with a certain accuracy *)
<x axis> = "XAXIS", "(", <shift>, ")";
<y axis> = "YAXIS", "(", <shift>, ")";
<z axis> = "ZAXIS", "(", <shift>, ")";
<dx> = <real>;
<dy> = <real>;
<dz> = <real>;
<alpha> = <real>;
<beta> = <real>;
<gamma> = <real>;
<real> = ? literal convertible to IEEE double precision floating point number ?;

The following mounting parameters may be specified:

time lag: This is the time offset $\Delta t$ of the trajectory time system with respect to the laser scanner system (i.e., this value has to be subtracted from a trajectory time stamp in order to obtain the same point of time in the scanner's time system)

scanner system: A sequence of 3 directions out of the set {F(ront),B(ack),L(eft),R(ight),U(p),D(own)} is used to specify the approximate attitude of the zero-tilted scanner system $S_0$ - note: NOT that of the scanner system $S$ (!) - with respect to the body system $B$ (see figure 1). The direction sequence is restricted to the 24 possible right-handed systems. Based on this input, the rotation matrix $\mathbf{R}_{S_0}^M$ (mentioned in the previous section) is determined. For example, in case of figure 1, the 3 axes of system $S_0$ point down, back, and left (here, it is not specified which of the axes is X,Y,Z). Since we support only right-handed systems, the set-up of this example corresponds to one of the 3 possible sequences "D-B-L", "B-L-D", "L-D-B".

For example, "SCANNERSYS(D-B-L)" yields $ \mathbf{R}_{S_0}^M = \begin{pmatrix} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 1 & 0 & 0 \end{pmatrix} $.

mounting rotation: The rotation from the system $M$ into the system $B$

mounting shift: The translation (shift) from the system $S_0$ into the system $B$

tilt rotation: The rotation from the system $S$ into the system $S_0$

tilt shift: The translation (shift) from the system $S$ into the system $S_0$

Hence, from the point of view of these transformation elements, the above-mentioned systems take role as a global or a local system.

The next section describes how to specify a rotation and/or a shift using the provided syntax.

How to specify a rotation and/or a shift

Using the above-introduced matrix notation, a congruent (6-parameter) transformation from a local system $L$ into a global system $G$ may be described as:

$ \mathbf{x}^G=\mathbf{R}_L^G\cdot\mathbf{x}^L + t_L^G \Longleftrightarrow \mathbf{x}^L={(\mathbf{R}_L^G)}^T\cdot(\mathbf{x}^G-t_L^G) $

In the same way, the inverse transformation may be described as:

$ \mathbf{x}^L=\mathbf{R}_G^L\cdot\mathbf{x}^G + t_G^L \Longleftrightarrow \mathbf{x}^G={(\mathbf{R}_G^L)}^T\cdot(\mathbf{x}^L-t_G^L) $

Consequently, the following relations hold:

$ \mathbf{R}_L^G={(\mathbf{R}_G^L)}^T $ and $ t_L^G=-\mathbf{R}_L^G\cdot t_G^L $

The above-decribed syntax allows the user to specify a rotation and/or a shift in different ways:

Specifying a shift:

Depending on the refsys parameter, a shift vector may be specified either as $t_L^G$ ("GLOBAL": the shift vector corresponds to the position of the local system's origin expressed in the global system) or as $t_G^L$ ("LOCAL": the shift vector corresponds to the position of the global system's origin expressed in the local system).

Note: Since the system $M$ is just an auxiliary system with hardly practical relevance, the mounting shift expressed in the local system is NOT $t_B^M$ BUT $t_B^{S_0}$ (!). Hence, the mounting shift expressed in the global system $B$ may be determined by $ t_{S_0}^B=t_M^B=-\mathbf{R}_{S_0}^B\cdot t_B^{S_0}=-\mathbf{R}_M^B\cdot\mathbf{R}_{S_0}^M\cdot t_B^{S_0} $

Specifying a rotation:

Basically, there are three possibilities to specify a rotation:

  • by 2 or 3 vectors defining the direction(s) of 2 or 3 axes of one system expressed in the other system (using the refsys parameter as described above) where the vectors need not be normalized but they must be orthogonal to each other with a certain accuracy (1.e-6)
  • by a matrix given in column major order (i.e. r11 r21 r31 r12... where r11, r21, r31 are the elements of the first column) which must fulfill all orthogonality conditions with a certain accuracy (1.e-6). Depending on the refsys parameter, the matrix is interpreted either as $\mathbf{R}_L^G$ ("GLOBAL") or $\mathbf{R}_G^L$ ("LOCAL")
  • by 3 rotation angles (as described right below)

Specifying a rotation by 3 angles:

Note: when specifying a rotation by 3 angles $\alpha$, $\beta$, $\gamma$, the resulting matrix will be anyway (!) $\mathbf{R}_L^G$ (i.e. rotation from the local into the global system). The refsys parameter specifies which system the angles are related to (i.e. around the axes of which system the rotations are performed). Depending on the refsys parameter, the entire rotation matrix is created either by right or by left multiplication of the sub-matrices $\mathbf{R}_1(\alpha)$, $\mathbf{R}_2(\beta)$, $\mathbf{R}_3(\gamma)$:

"GLOBAL": $ \mathbf{R}_L^G=\mathbf{R}_1(\alpha)\cdot\mathbf{R}_2(\beta)\cdot\mathbf{R}_3(\gamma) $

"LOCAL": $ \mathbf{R}_L^G=\mathbf{R}_3(\gamma)\cdot\mathbf{R}_2(\beta)\cdot\mathbf{R}_1(\alpha) $

Depending on the parameter axis hierarchy, the submatrices $\mathbf{R}_k(\varphi_k)$ are replaced by $\mathbf{R}_x(\varphi_x)$, $\mathbf{R}_y(\varphi_y)$, $\mathbf{R}_z(\varphi_z)$, where

$ \mathbf{R}_x(\varphi)=\begin{pmatrix}1 & 0 & 0 \\ 0 & \cos\varphi_x & -\sin\varphi_x \\ 0 & \sin\varphi_x & \cos\varphi_x \end{pmatrix} $ and $ \mathbf{R}_y(\varphi)=\begin{pmatrix}\cos\varphi_y & 0 & \sin\varphi_y \\ 0 & 1 & 0 \\ -\sin\varphi_y & 0 & \cos\varphi_y \end{pmatrix} $ and $ \mathbf{R}_z(\varphi)=\begin{pmatrix}\cos\varphi_z & -\sin\varphi_z & 0 \\ \sin\varphi_z & \cos\varphi_z & 0 \\ 0 & 0 & 1 \end{pmatrix} $.

The default rotation sense of the angles $\alpha$, $\beta$, $\gamma$ is counter-clockwise ("CCW") (mathematically positive when looking against the respective rotation axis). It may be as well specified as clockwise ("CW") using the parameter rotation sense: in this case, all the 3 angles are multiplied by (-1) before being evaluated in the submatrices.

The parameter units allows the user to select the angle units, where "DEG", "GRAD", or "RAD" may be specified (the default is "DEG").

For example, the specification of the mounting rotation matrix by three angles $\alpha$ (roll), $\beta$ (pitch), $\gamma$ (yaw) according to the aviation ARINC 705 (see section Alignment of the aircraft) corresponds to refsys="=LOCAL", axis hierarchy="AXISHIERARCHY(X-Y-Z)" and rotation sense="SENSEOFROT(CCW)". $ \Longrightarrow \mathbf{R}_L^G=\mathbf{R}_z(\gamma)\cdot\mathbf{R}_y(\beta)\cdot\mathbf{R}_x(\alpha) $

Examples for specifying the mounting calibration

For the following examples, we define $ \rho=\frac{180}{\pi} $.

Example 1:

User input:

"SCANNERSYS(D-F-R), MOUNTROTATION=LOCAL(ANGLES(0.07346 0.2479 -0.37684)), MOUNTSHIFT(-0.7834 0.193422 0.07165)"

This is equivalent to the complete input:

"TIMELAG(0.0),
SCANNERSYS(D-F-R),
MOUNTROTATION=LOCAL(ANGLES(0.07346 0.2479 -0.37684), AXISHIERARCHY(X-Y-Z), SENSEOFROT(CCW), UNITS(DEG)),
MOUNTSHIFT=GLOBAL(-0.7834 0.193422 0.07165),
TILTROTATION=GLOBAL(MATRIX(1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0)),
TILTSHIFT=GLOBAL(0.0 0.0 0.0)"

It yields:

$\Delta t=0.0$, $ \mathbf{R}_{S_0}^M = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} $, $ \mathbf{R}_M^B = \mathbf{R}_z(\frac{-0.37684}{\rho})\cdot\mathbf{R}_y(\frac{0.2479}{\rho})\cdot\mathbf{R}_x(\frac{0.07346}{\rho})= \begin{pmatrix} 0.9999690 & 0.0065826 & 0.0043181 \\ -0.0065770 & 0.9999775 & -0.0013105 \\ -0.0043267 & 0.0012821 & 0.9999898 \end{pmatrix} $, $ \mathbf{t}_M^B = \mathbf{t}_{S_0}^B = \begin{pmatrix}-0.7834 \\ 0.193422 \\ 0.07165\end{pmatrix} $, $ \mathbf{R}_S^{S_0} = \mathbf{I} $, $ \mathbf{t}_S^{S_0} = \mathbf{o} $

Example 2:

User input:

"TIMELAG(-0.007), SCANNERSYS(R-F-U), MOUNTROTATION=GLOBAL(ANGLES(0.07346 0.2479 -0.37684), AXISHIERARCHY(Y-X-Z)),
MOUNTSHIFT(-0.927 0.014 0.053),TILTROTATION(VECTORS(XAXIS(1 0 2),ZAXIS(-2 0 1))),TILTSHIFT=LOCAL(0.1204 0.0564 -0.0134)"

This is equivalent to the complete input:

"TIMELAG(-0.007),
SCANNERSYS(R-F-U),
MOUNTROTATION=GLOBAL(ANGLES(0.07346 0.2479 -0.37684), AXISHIERARCHY(Y-X-Z), SENSEOFROT(CCW), UNITS(DEG)),
MOUNTSHIFT=GLOBAL(-0.927 0.014 0.053),
TILTROTATION=GLOBAL(VECTORS(XAXIS(1.0 0.0 2.0), YAXIS(0.0 1.0 0.0), ZAXIS(-2.0 0.0 1.0))),
TILTSHIFT=LOCAL(0.1204 0.0564 -0.0134)"

It yields:

$\Delta t=-0.007$, $ \mathbf{R}_{S_0}^M = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix} $, $ \mathbf{R}_M^B = \mathbf{R}_y(\frac{0.07346}{\rho})\cdot\mathbf{R}_x(\frac{0.2479}{\rho})\cdot\mathbf{R}_z(\frac{-0.37684}{\rho})= \begin{pmatrix} 0.9999775 & 0.0065826 & 0.0012821 \\ -0.0065770 & 0.9999690 & -0.0043267 \\ -0.0013105 & 0.0043181 & 0.9999898 \end{pmatrix} $, $ \mathbf{t}_M^B = \mathbf{t}_{S_0}^B = \begin{pmatrix}-0.927 \\ 0.014 \\ 0.053\end{pmatrix} $,

$ \mathbf{R}_S^{S_0} = \begin{pmatrix} 0.4472136 & 0.0000000 & -0.8944272 \\ 0.0000000 & 1.0000000 & 0.0000000 \\ 0.8944272 & 0.0000000 & 0.4472136 \end{pmatrix} $, $ \mathbf{t}_{S_0}^S = \begin{pmatrix} 0.1204 \\ 0.0564 \\ -0.0134 \end{pmatrix} \Longrightarrow \mathbf{t}_S^{S_0} = -\mathbf{R}_S^{S_0}\cdot\mathbf{t}_{S_0}^S = \begin{pmatrix}-0.0658 \\ -0.0564 \\ -0.1017 \end{pmatrix} $

Example 3:

User input:

"MOUNTROTATION(MATRIX(0.9999690 -0.0065770 -0.0043267 0.0065826 0.9999775 0.0012821 0.0043181 -0.0013105 0.9999898)),
SCANNERSYS(D-F-R), MOUNTSHIFT=LOCAL(-0.068013 0.784958 -0.188353)"

This is equivalent to the complete input:

"TIMELAG(0.0),
SCANNERSYS(D-F-R),
MOUNTROTATION=GLOBAL(MATRIX(0.999969 -0.006577 -0.0043267 0.0065826 0.9999775 0.0012821 0.0043181 -0.0013105 0.9999898)),
MOUNTSHIFT=LOCAL(-0.068013 0.784958 -0.188353),
TILTROTATION=GLOBAL(MATRIX(1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0)),
TILTSHIFT=GLOBAL(0.0 0.0 0.0)"

It yields the same mounting calibration as example 1:

$\Delta t=0.0$, $ \mathbf{R}_{S_0}^M = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix} $, $ \mathbf{R}_M^B =\begin{pmatrix} 0.9999690 & 0.0065826 & 0.0043181 \\ -0.0065770 & 0.9999775 & -0.0013105 \\ -0.0043267 & 0.0012821 & 0.9999898 \end{pmatrix} $, $ \mathbf{t}_B^{S_0}= \begin{pmatrix} -0.068013 \\ 0.784958 \\ -0.188353 \end{pmatrix} \Longrightarrow \mathbf{t}_M^B =\mathbf{t}_{S_0}^B =-\mathbf{R}_M^B\cdot\mathbf{R}_{S_0}^M\cdot\mathbf{t}_B^{S_0}=\begin{pmatrix}-0.7834 \\ 0.193422 \\ 0.07165\end{pmatrix} $,

$ \mathbf{R}_S^{S_0} = \mathbf{I} $, $ \mathbf{t}_S^{S_0} = \mathbf{o} $

Example 4:

User input:

"MOUNTROTATION=LOCAL(MATRIX(0.999969 -0.006577 -0.0043267 0.0065826 0.9999775 0.0012821 0.0043181 -0.0013105 0.9999898)),
SCANNERSYS(D-B-L),MOUNTSHIFT=LOCAL(-0.068013 0.784958 -0.188353),TILTROTATION=LOCAL(VECTORS(XAXIS(5 0 1),YAXIS(0 1 0)))"

This is equivalent to the complete input:

"TIMELAG(0.0),
SCANNERSYS(D-B-L),
MOUNTROTATION=LOCAL(MATRIX(0.999969 -0.006577 -0.0043267 0.0065826 0.9999775 0.0012821 0.0043181 -0.0013105 0.9999898)),
MOUNTSHIFT=LOCAL(-0.068013 0.784958 -0.188353),
TILTROTATION=LOCAL(VECTORS(XAXIS(5.0 0.0 1.0),YAXIS(0.0 1.0 0.0), ZAXIS(-1.0 0.0 5.0))),
TILTSHIFT=GLOBAL(0.0 0.0 0.0)"

It yields:

$\Delta t=0.0$, $ \mathbf{R}_{S_0}^M = \begin{pmatrix} 0 & -1 & 0 \\ 0 & 0 & -1 \\ 1 & 0 & 0 \end{pmatrix} $, $ \mathbf{R}_B^M =\begin{pmatrix} 0.9999690 & 0.0065826 & 0.0043181 \\ -0.0065770 & 0.9999775 & -0.0013105 \\ -0.0043267 & 0.0012821 & 0.9999898 \end{pmatrix}$ $ \Longrightarrow \mathbf{R}_M^B={(\mathbf{R}_B^M)}^T=\begin{pmatrix} 0.9999690 & -0.0065770 & -0.0043267 \\ 0.0065826 & 0.9999775 & 0.0012821 \\ 0.0043181 & -0.0013105 & 0.9999898 \end{pmatrix} $,

$ \mathbf{t}_B^{S_0}= \begin{pmatrix} -0.068013 \\ 0.784958 \\ -0.188353 \end{pmatrix} \Longrightarrow \mathbf{t}_M^B =\mathbf{t}_{S_0}^B =-\mathbf{R}_M^B\cdot\mathbf{R}_{S_0}^M\cdot\mathbf{t}_B^{S_0}=\begin{pmatrix}-0.7859 \\ -0.183094 \\ 0.07165\end{pmatrix} $,

$ \mathbf{R}_{S_0}^S = \begin{pmatrix} 0.9805807 & 0.0000000 & -0.1961161 \\ 0.0000000 & 1.0000000 & 0.0000000 \\ 0.1961161 & 0.0000000 & 0.4472136 \end{pmatrix} \Longrightarrow \mathbf{R}_S^{S_0} = {(\mathbf{R}_{S_0}^S)}^T = \begin{pmatrix} 0.9805807 & 0.0000000 & 0.1961161 \\ 0.0000000 & 1.0000000 & 0.0000000 \\ -0.1961161 & 0.0000000 & 0.4472136 \end{pmatrix}$, $ \mathbf{t}_S^{S_0} = \mathbf{o} $

Alignment of the aircraft

According to the aviation norm ARINC 705 (Bäumker and Heimes, 2001), we use the following definitions (Figure 2):

  • Primary rotation: roll angle $\phi$ around the aircraft's longitudinal axis (x), where $\phi$ counts positive if the right wing points down (w.r.t. the local horizon)
  • Secondary rotation: pitch angle $\theta$ around the aircraft's lateral axis (y), where $\theta$ counts positive if the nose points up (w.r.t. the local horizon)
  • Tertiary rotation: yaw angle $\psi$ around the aircraft's vertical axis (z), where $\psi$ is the azimuth of the flight direction (counted clockwise starting from north)
alignmentOfAircraft_100dpi.png
Fig. 2: Alignment of aircraft described by aviation norm ARINC 705

The entire rotation matrix $\mathbf{R}_B^H$ transforms a point from the body system $B$ into the system of the local horizon $H$:

$\mathbf{x}^H=\mathbf{R}_B^H\cdot\mathbf{x}^B$

where

$\mathbf{x}^B={\begin{pmatrix}x^B & y^B & z^B \end{pmatrix}}^T$, $\mathbf{x}^H={\begin{pmatrix}x^H & y^H & z^H \end{pmatrix}}^T$

and

$\mathbf{R}_B^H=\mathbf{R}_z(\psi)\cdot\mathbf{R}_y(\theta)\cdot\mathbf{R}_x(\phi)=$

$=\begin{pmatrix}\cos\psi & -\sin\psi & 0 \\ \sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{pmatrix}\cdot \begin{pmatrix}\cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{pmatrix}\cdot \begin{pmatrix}1 & 0 & 0 \\ 0 & \cos\phi & -\sin\phi \\ 0 & \sin\phi & \cos\phi \end{pmatrix}=$ $\begin{pmatrix}\cos\psi\cdot\cos\theta & \cos\psi\cdot\sin\theta\cdot\sin\phi-\sin\psi\cdot\cos\phi & \cos\psi\cdot\sin\theta\cdot\cos\phi+\sin\psi\cdot\sin\phi \\ \sin\psi\cdot\cos\theta & \sin\psi\cdot\sin\theta\cdot\sin\phi+\cos\psi\cdot\cos\phi & \sin\psi\cdot\sin\theta\cdot\cos\phi-\cos\psi\cdot\sin\phi \\ -\sin\theta & \cos\theta\cdot\sin\phi & \cos\theta\cdot\cos\phi \end{pmatrix}$

Note: In compliance with the notation of the previous section, the columns of the rotation matrix $\mathbf{R}_B^H$ correspond to the unit vectors of the axes of system B expressed in system H.

Note: If a map projection is used as global coordinate system (see next section), the yaw angle $\psi$ will be automatically corrected for the meridian convergence. In this case, the system $H$ is not oriented towards True North but towards Grid North.

World coordinate system

The world coordinate system may be specified by the global parameter coord_ref_sys.

Basically, two types of world coordinate systems are supported by this module:

  • Geocentric, ECEF (Earth-Centered, Earth-Fixed) Cartesian coordinate systems (e.g. WGS84-XYZ)
  • Conformal Transverse Mercator projections (e.g. UTM)

Note: Geographic coordinates (latitude, longitude, ellipsoidal height) have to be converted into one of the above-mentioned system types first.

The specification of the world coordinate system should contain the parameters $a$ (semi major axis) and $f$ (flattening) of the reference ellipsoid (if not, the parameters of the WGS84 ellipsoid are used instead). In case of a transverse mercator projection, the parameters false_easting (default: 0.0), false_northing (default: 0.0), and scale_factor (default: 1.0) should be contained, too.

For example, the system "UTM zone 33 - northern hemisphere" may be specified either by the WKT string

WKT:PROJCS["UTM Zone 33, Northern Hemisphere",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",15],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]

or by its EPSG code

EPSG:32633

In case of a geocentric system, a ground point given in the local horizon system $H$ is transformed into the world system $W$ by:

$\mathbf{x}^W=\mathbf{R}_H^W\cdot\mathbf{x}^H+\mathbf{t}_H^W$

where $\mathbf{R}_H^W$ is a rotation matrix, which describes the attitude of the local horizon system $H$ with respect to the world system $W$, and where $\mathbf{t}_H^W=\mathbf{t}_B^W$ is the current (interpolated) trajectory point given in the world system $W$.

In case of a Transverse Mercator map projection, beside the meridian convergence (already considered by correction of the yaw angle), the projection's length distortion and earth curvature are taken into account to transform a ground point into the "projected" horizon system $H_P$ (having the same origin as $H$), from where it is finally translated into the world system $W$:

$ \mathbf{x}^W = \mathbf{x}^{H_P}+\mathbf{t}_H^W = F(\mathbf{x}^H, \mathbf{t}_H^W, a, f)+\mathbf{t}_H^W $

where $ F(\mathbf{x}^H, \mathbf{t}_H^W, a, f) $ contains the correction terms necessary due to the (space-dependent) length distortion and earth curvature

Trajectory file

At least one trajectory file is required for this module.

A trajectory file must contain a set of at least 2 records containing the following data:

  • 3 coordinates (x,y,z) of the position given in a supported world coordinate system (see previous section)
  • the gps time stamp (t)
  • 3 rotation angles [radiant] which must comply with the above-mentioned aviation norm ARINC 705 (roll, pitch, yaw)

The module tries to recognise the trajectory format based on its content if not specified. The standard input format is the "trajectory file format (trj)" which supports ASCII and binary mode. The 7 columns can be given in two different orders:

  • x, y, z, t, roll, pitch, yaw (angle unit: degree)
  • t, x, y, z, roll, pitch, yaw (angle unit: degree)

Since the gps time is excepted to be in ascending order, the two formats can be easily differentiated without any specification. Whereas all angles columns are stored in radiant within the ODM (see ODM as a database table), the trajectory format assumes the angle to be given in degrees. Hence, while reading the angles are automatically converted. Please note, that in binary mode the the first 4 values are read as doubles (8 bytes per values) and the 3 angles as floats (4 bytes per value).

In case your trajecotry file doesn't match the format described above, it is possible to create an appropriate OPALS format description (see OPALS Generic Format) and specify those as tFormat parameter. Please note that angle values are expected in radians which can be simply realised specifying an appropriate transformation in the OPALS format description file.

If two ore more trajectory files are given, their time ranges (intervals between first and last entry) must not overlap. Furthermore, the time range of a scanner file must be completely within the time range of one trajectory.

Parameter description

-inFileinput file(s)
Type: opals::Vector<opals::Path>
Remarks: mandatory
input file(s) containing laser scanner data given in device coordinates
-iFormatinput file format [auto, sdc, odm, las, fwf]
Type: opals::Vector<opals::String>
Remarks: estimable
input file format [auto, sdc, odm, las, fwf]
Estimation rule: from input file extension
-mountingmounting calibration parameters
Type: opals::MountingPars
Remarks: default=TIMELAG(0), SCANNERSYS(F-R-D), MOUNTSHIFT=GLOBAL (+0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00), MOUNTROTATION=GLOBAL (+1.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +1.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +1.000000000000e+00), TILTSHIFT=GLOBAL (+0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00), TILTROTATION=GLOBAL (+1.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +1.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +0.000000000000e+00 +1.000000000000e+00),
mounting calibration parameters
-trjFiletrajectory file(s)
Type: opals::Vector<opals::Path>
Remarks: mandatory
ODM file(s) containing trajectory data
-tFormattrajectory file format [auto, trj, odm, <opals format def. xml file>]
Type: opals::Vector<opals::String>
Remarks: estimable
In case of multiple trajectory files with present tFormat parameter, the number of tFormat and trjFiles entries have to match or only a single tFormat entry is specified (=all trajectory files have the same format).
Estimation rule: auto. The file content is used to recognise the file format.
-outFileoutput file(s)
Type: opals::Vector<opals::Path>
Remarks: estimable
output file(s) containing laser scanner data given in world system (given by trajectory)
Estimation rule: inFile_+'dirGeoref'+.ext (.ext derived from oFormat)
-oFormatoutput file format [sdw, odm, las, fwf]
Type: opals::String
Remarks: mandatory
output file format [sdw, odm, las, fwf]
Estimation rule: from output file extension (default='odm' if outFile is not given)

Examples

The data used in the following examples are located in the demo directory.

Example 1

We start from the given demo file 060622_173547.sdf. First, Module Fullwave has to be applied to this file in order to generate the file 060622_173547.sdc:

opalsFullwave -infile 060622_173547.sdf -detectThrLow 10 -outfile 060622_173547.sdc

The associated trajectory file F002_20060622_UTM33_part.txt is given in UTM 33, the global parameter coord_ref_sys may be set as described in the example above.

In this example, we have a Down-Front-Right scanner system (without tilting). The time lag is -3.2 milliseconds.

opalsDirectGeoref -inFile 060622_173547.sdc -mounting "SCANNERSYS(D-F-R), MOUNTROTATION=LOCAL(ANGLES(-0.060 0.072 -0.325)), MOUNTSHIFT(-0.052 0.151 0.083), TIMELAG(-0.0032)" -trjFile F002_20060622_UTM33_part.txt -oformat las -coord_ref_sys EPSG:32633

As result, we obtain the LAS file 060622_173547_dirGeoref.las containing the georeferenced point cloud.

Example 2

We are starting from the demo file stripPart.sdc. Since the trajectory file trjPart_geoc.txt is given in geocentric WGS 84 coordinates, we set the global parameter coord_ref_sys to

WKT:GEOCCS["WGS 84 (geocentric)",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],
UNIT["m",1.0],
AXIS["Geocentric X",OTHER],
AXIS["Geocentric Y",EAST],
AXIS["Geocentric Z",NORTH],
AUTHORITY["EPSG","4328"]]

or to

EPSG:4328

In this example, we have a Down-Back-Left scanner system (without tilting). There is no time lag.

opalsDirectGeoref -inFile stripPart.sdc -mounting "SCANNERSYS(D-B-L), MOUNTROTATION=LOCAL(ANGLES(0.09685 0.4549 -0.21325)), MOUNTSHIFT(-0.6633 0.146125 0.05575)" -trjFile trjPart_geoc.txt -oformat LAS_1.4_ptRecFmt_6.xml -coord_ref_sys EPSG:4328

As result, we obtain the LAS file stripPart_dirGeoref.las containing the georeferenced point cloud.

References

Bäumker, M., F.-J. Heimes (2001): New Calibration and Computing Method for Direct Georeferencing of Image and Scanner Data Using the Position and Angular Data of an Hybrid Navigation System. In: Proceedings OEEPE-Workshop Integrated Sensor Orientation, Hannover, 17.-18. Sept. 2001.