pointfit: Least squares fit of point clouds, or the 'Procrustes'...

In cwhmisc: Miscellaneous Functions for Math, Plotting, Printing, Statistics, Strings, and Tools

Description

Find a transformation which consists of a translation tr and a rotation Q multiplied by a positive scalar f which maps a set of points x into the set of points xi: xi = f * Q * x + tr + error. The resulting error is minimized by least squares.

Usage

  pointfit(xi,x)

Arguments

x	Matrix of points to be mapped. Each row corresponds to one point.
xi	Matrix of target points. Each row corresponds to one point.

Details

The optimisation is least squares for the problem xi: xi = f * Q * x + tr. The expansion factor f is computed as the geometric mean of the quotients of corresponding coordinate pairs. See the program code.

Value

A list containing the following components:

Q	The rotation.
f	The expansion factor.
tr	The translation vector.
res	The residuals xi - f * Q * x + tr.

Author(s)

Walter Gander, [e-mail],
http://www.inf.ethz.ch/personal/gander/ adapted by Christian W. Hoffmann <[e-mail]>

References

"Least squares fit of point clouds" in: W. Gander and J. Hrebicek, ed., Solving Problems in Scientific Computing using Maple and Matlab, Springer Berlin Heidelberg New York, 1993, third edition 1997.

See Also

rotL to generate rotation matrices

Examples

  # nodes of a pyramid
  A <- matrix(c(1,0,0,0,2,0,0,0,3,0,0,0),4,3,byrow=TRUE)
  nr <- nrow(A)
  v <- c(1,2,3,4,1,3,4,2)  # edges to be plotted
#  plot
  # points on the pyramid
  x <-
matrix(c(0,0,0,0.5,0,1.5,0.5,1,0,0,1.5,0.75,0,0.5,
2.25,0,0,2,1,0,0),
    7,3,byrow=TRUE)
  # simulate measured points
  # theta <- runif(3)
  theta <- c(pi/4, pi/15, -pi/6)
  # orthogonal rotations to construct Qr
  Qr <- rotL(theta[3]) %*% rotL(theta[2],1,3) %*% rotL(theta[1],1)
  # translation vector
  # tr <- runif(3)*3
  tr <- c(1,3,2)
  # compute the transformed pyramid
  fr <- 1.3
  B <- fr * A %*% Qr + outer(rep(1,nr),tr)
  # distorted points
  # xi <- fr * x + outer(rep(1,nr),tr) + rnorm(length(x))/10
  xi <- matrix(c(0.8314,3.0358,1.9328,0.9821,4.5232,2.8703,1.0211,3.8075,1.0573,
0.1425,4.4826,1.5803,0.2572,5.0120,3.1471,0.5229,4.5364,3.5394,1.7713,
3.3987,1.9054),7,3,byrow=TRUE)
  (pf <- pointfit(xi,x))
  # the fitted pyramid
  (C <- A %*% pf$Q + outer(rep(1,nrow(A)),pf$tr))  ## !!!!!!  %*% instead of %*%
  # As a final check we generate the orthogonal matrix S from the computed angles
  # theta and compare it with the result pf$Q
  Ss <- rotL(theta[3]) 
  range(svd(Ss*pf$factor - pf$Q)$d) #  6.652662e-17 1.345509e-01

Example output

Loading required package: lattice
Loading required package: grid
$Q
           [,1]       [,2]       [,3]
[1,]  0.9391757 -0.3939876 -0.1804864
[2,]  0.3931310  0.5933228  0.7505123
[3,] -0.1823447 -0.7500630  0.6884828

$tr
[1] 0.5222852 4.6685166 1.3809854

$factor
[1] 1.034337

$res
             [,1]        [,2]        [,3]
[1,] -0.026496150 -0.01443422 -0.09182002
[2,] -0.074654401  0.15063180 -0.09587187
[3,]  0.087603596 -0.03262252 -0.12608467
[4,]  0.010950084 -0.02050265  0.16441235
[5,]  0.002392087 -0.02354832 -0.05157485
[6,]  0.025976675 -0.01485885  0.13781436
[7,] -0.025771890 -0.04466524  0.06312470

            [,1]     [,2]     [,3]
[1,]  1.46146090 4.274529 1.200499
[2,]  1.30854719 5.855162 2.882010
[3,] -0.02474899 2.418328 3.446434
[4,]  0.52228516 4.668517 1.380985
[1] 6.035174e-17 8.471175e-01

cwhmisc documentation built on May 1, 2019, 7:55 p.m.