4

I have some double points (the same points in 2 different datums). I'd like to calculate the Helmert parameters but I don't understand how.

NOTE: I don't need software. I need an algorithm in order to write a small Delphi application.

9

The following relies on the Wikipedia article on seven-parameter Helmert transformations.

The data ("double points") consist of ordered pairs ((x,y,z), (x',y',z')) where (x,y,z) are earth-centered Cartesian coordinates in the source datum and (x',y',z') are the corresponding points in the target datum, all measured in meters. The latter are presumed derived from the former by a combination of a rotation, an isothety (uniform rescaling), and translation.

There are several ideas underlying the algorithm:

  1. Because a change in datums is so slight, all three kinds of transformations can be replaced by infinitesimal equivalents.

    • The infinitesimal equivalent of a translation is the translation itself, given by a vector of displacements in meters, (cx, cy, cz).

    • A small isothety has a scale factor f very close to 1. Its infinitesimal equivalent is the first-order term in the expansion of log(f) around 1, equal to f - 1 = s. Typically s is given in parts per million (ppm); that is, it is multiplied by 10^6.

    • A small rotation has an infinitesimal equivalent of the form

      x' = x - rz*y + ry*z
      y' = y - rx*z + rz*x
      z' = z - ry*x + rx*y
      

      where (rx, ry, rz) are given in radians. (The Wikipedia article is not explicit about this.)

  2. The relationship between (x,y,z) and (x',y',z') is a linear function of the related parameters (cx, cy, cz, 1+s, (1+s)*rx, (1+s)*ry, (1+s)*rz) = (b1, b2, ..., b7) = b. Estimates of the last four Helmert parameters can be recovered from estimates of b with simple arithmetic: s = b4-1, rx=b5/b4, ry=b6/b4, and rz=b7/b4. Care is needed to convert properly among the units of measurement--see the code example below.

  3. Although technically solving for b requires multivariate least squares, it can actually be solved using a standard (univariate) ordinary least squares algorithm. This is done by turning each "double point" ((x,y,z), (x',y',z')) into three separate rows of data having seven independent variables and just one dependent variable, as follows:

    1, 0, 0, x, 0, z, -y; x'
    0, 1, 0, y, -z, 0, x; y'
    0, 0, 1, z, y, -x, 0; z'
    

Thus if there were n "double points" to begin with, this derived dataset will have 3*n rows and 8 columns. The first seven columns correspond to b while the eighth evidently gives the components of all the target points. For example, we can read the first row as stating that

1*cx + 0*cy + 0*cz + x*(1+s) + 0*(1+s)*rx + z*(1+s)*ry - y*(1+s)*rz = x'.

Equivalently,

1*b1 + 0*b2 + 0*b3 + x*b4 + 0*b5 + z*b6 - y*b7 = x'.

It's straightforward to check that this is precisely what the infinitesimal rotation (equal to x + ryz - rzy), followed by the rescaling (which multiplies everything by 1+s), followed by a translation by cx, does in order to transform x into x'. The second and third rows work similarly to describe the transformations from y to y' and z to z'. By combining the data this way, we implicitly make the assumption that the errors in the three coordinates are independent and, at least approximately, have similar variances.

After having found the least-squares estimate of b using any standard solver, you can back-calculate the estimates of s, rx, ry, and rz. Note that this least-squares model does not have an intercept term. (It actually has three distinct intercepts, cx, cy, and cz, and so it doesn't need any more!)


Here is working R code to carry out the calculations involved in (1) applying a Helmert transformation and (2) estimating the parameters from a set of "double points." Of especial note are the conversions among the various units of measurement: meters, radians, arcseconds, and ppm. These are the kinds of details that are easy to mess up (and glossed over in the Wikipedia article, I'm afraid).

#
# Apply a 7-parameter Helmert transformation.
#
helmert <- function(x, c0=c(0,0,0), s=0, r=c(0,0,0)) {
  # x has 3D coordinates in *columns*
  # c is displacement in meters
  # s is infinitesimal scale in ppm
  # r is infinitesimal rotation in arcseconds
  rr <- r/(60*60*180) * pi # Convert to radians
  a <- matrix(1, 3, 3)
  a[upper.tri(a)] <- c(-rr[3], rr[2], -rr[1])
  a[lower.tri(a)] <- -a[upper.tri(a)]
  return ((1 + s/10^6) * a %*% x + c0)
}
#
# Estimate a Helmert transformation.
#
helmert.fit <- function(x, xp) {
  x.d <- data.frame(matrix(apply(x, 2, 
                                 function(y) matrix(c(1,0,0,y[1],0,y[3],-y[2],
                                                      0,1,0,y[2],-y[3],0,y[1],
                                                      0,0,1,y[3],y[2],-y[1],0),
                                                    nrow=3)), ncol=7, byrow=TRUE))
  x.d$y <- as.vector(xp)
  names(x.d) <- c("cx", "cy", "cz", "s", "rx", "ry", "rz", "y")
  rownames(x.d) <- outer(c("x","y","z"), 1:ncol(x), paste, sep="") # `ncol(x)` was `n`
  fit <- lm(y ~ . - 1, data=x.d)
  #
  # Convert the solution back to (cx, cy, cz, s, rx, ry, rz) in standard
  # units of measurement (ppm for `s` and arcseconds for `rx`, `ry`, `rz`).
  #
  parameters <- coef(fit) * c(1,1,1,1,60*60*180/pi*c(1,1,1))
  s1 <- parameters["s"]
  parameters <- parameters / c(1,1,1,1,s1,s1,s1)
  parameters["s"] <- (s1 - 1)*10^6
  return (list(parameters=parameters, fit=fit))
}

To illustrate its use, let's pick three points on the earth at random, apply a known Helmert transformation to them, add a little measurement error to the result, fit the parameters, and compare the fit to the known values.

#
# WGS84 --> Germany Bessel 1841:
#
c0 <- c(-582, -105, -414) # Meters
s <- -1.1                 # ppm
r <- c(-0.02, 0.26, 0.13) # Arc-seconds
actual <- c(c0, s, r)
names(actual) <-c("cx", "cy", "cz", "s", "rx", "ry", "rz")
#
# Generate sample data.
#
set.seed(17)
n <- 3
R <- 6371000
x <- apply(matrix(rnorm(3*n), nrow=3), 2, function(y) R * y / sqrt(sum(y*y)))
xp <- helmert(x, c0, s, r) +  rnorm(n*3, sd=0.1)
#
# Fit a seven-point Helmert transformation.
#
solution <- helmert.fit(x, xp)
rbind(actual, estimate=solution$parameters)

The output is

               cx       cy        cz         s          rx        ry       rz
actual   -582.000 -105.000 -414.0000 -1.100000 -0.02000000 0.2600000 0.130000
estimate -581.948 -104.893 -414.0355 -1.110138 -0.01565243 0.2635923 0.131022

You can see it's pretty good, even based on just three points. (It helps that they were spaced far apart on the earth's surface.) As a double-check, re-run this calculation with no error:

xp <- helmert(x, c0, s, r)
solution <- helmert.fit(x, xp)
rbind(actual, estimate=solution$parameters)

           cx   cy   cz    s    rx   ry   rz
actual   -582 -105 -414 -1.1 -0.02 0.26 0.13
estimate -582 -105 -414 -1.1 -0.02 0.26 0.13

The fit is perfect.

For those who understand least-squares regression, it is instructive to extract error estimates of the parameters from the fit. In R this can be done by post-processing solution$fit, which is an lm object. The residuals (differences between fitted and actual values of (x',y',z')) are especially interesting, so I will show them here (for the first fit):

> round(residuals(solution$fit), 2)

   x1    y1    z1    x2    y2    z2    x3    y3    z3 
-0.06  0.01 -0.01  0.07 -0.04  0.01 -0.01  0.03  0.00 

They are labeled by component and data point: x1 is the x-component of the first double point, for instance. These values are in meters and they are pretty small, given that the typical measurement error added in this simulation was 0.1 meters. This gives us a rough sense of the errors that might typically be incurred in applying this transformation to other points.

  • 1
    In the earlier answer by @whuber I would like to suggest a correction. I think it should be rownames(x.d) <- outer(c("x","y","z"), 1:ncol(x), paste, sep="") and not rownames(x.d) <- outer(c("x","y","z"), 1:n, paste, sep="") the function helmert.fit() does not know what n is. n is defined outside the function and that's why it works. If you work with real data and have not defined n, helmert.fit() would not work. – user2955884 Jun 23 '15 at 9:03
  • @user2955884 Thank you for noticing that! I believe you are correct, so for the benefit of future readers I will edit this post to implement your suggestion. – whuber Jun 23 '15 at 14:39
  • In the earlier answer by @whuber I would like to ask the following questions: 1. We can take rx, ry, rz to be angles and not rotation matrices because we assume the angles are very small (see i.e. www.maths.dundee.ac.uk/gawatson/helmertrev.pdf). But how small is small? 2. I understand "For example, we can read the first row as stating that 1*cx + 0*cy + 0*cz + x*(1+s) + 0*(1+s)*rx + z*(1+s)*ry - y*(1+s)*rz = x' " but not the way it is implemented in helmert.fit(): `x.d <- data.frame(matrix(apply(x, 2, function(y) matrix(c(1,0,0,y[1],0,y[3],-y[2], 0,1,0,y[2],-y[3],0,y[1], 0,0,1,y[3],y[2],-y[1 – user2955884 Jun 24 '15 at 9:57
  • @user2955884 (1) The infinitesimal rotation matrix r is a first-order approximation to the rotation exp(r) = 1 + r + r ^2/2! + ... . By (the multivariate version of) Taylor's Theorem, the error equals | s ^2/2| for another matrix s close to r. Thus, "small" means that ignoring the squared coefficients of r would give positions within an acceptable error range. Since the coordinates are less than 10^7 m, 60 arcsec=3e-4 radians still gives errors less than a meter. (2) In helmert.fit, the vector y[1:3] corresponds to x, y, and z in the exposition. – whuber Jun 24 '15 at 12:54
2

The guidelines I use for "small" angles are:

  1. The difference between the radian value of the angle and sin of the angle is close enough to zero
  2. The cos of the angle is close enough to one.

In the unlikely case you are not dealing with small angles here is a list of some papers I found helpful when working on 3D transformations (rigid body):

  1. "Deakin, R.E. (1998) 3D Coordinate Transformations..." http://www.researchgate.net/publication/228608056_3-D_Coordinate_Transformations
  2. "Dewitt, B.A. (1996) Initial Approximations for the Three-Dimensional Conformal Coordinate Transformation..." http://asprs.org/a/publications/pers/96journal/january/1996_jan_79-83.pdf
  3. "Arun, K.S, et.al. (1987) Least-Squares Fitting of Two 3-D Point Sets..." http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=4767965&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D4767965
  4. "Sjöberg, L.E. (2013) Closed-form and iterative weighted least squares solutions of Helmert transformation parameters..." http://www.degruyter.com/view/j/jogs.2013.3.issue-1/jogs-2013-0002/jogs-2013-0002.xml

protected by Aaron Jun 24 '15 at 11:15

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.