pointfit {cwhmisc}R Documentation

Least squares fit of point clouds, or the Procrustes problem.

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 \dot Qx + tr. 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 \dot Qx + 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, gander@inf.ethz.ch,
http://www.inf.ethz.ch/personal/gander/ adapted by Christian W. Hoffmann <c-w.hoffmann@sunrise.ch>
http://www.wsl.ch/personal_homepages/hoffmann/index_EN

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

rotm to generate rotation matrices (mathlib), rotangle to determine them from orthogonal matrix.

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
  # thetar <- runif(3)
  thetar <- c(pi/4, pi/15, -pi/6)
  # orthogonal rotations to construct Qr
  Qr <- rotm(3,1,2,thetar[3]) %*% rotm(3,1,3,thetar[2]) %*% rotm(3,2,3,thetar[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))
  (theta <- rotangle(pf$Q))
  thetar
  # as a final check we generate the orthogonal matrix S from the computed angles
  # theta and compare it with the result pf$Q
  Ss <- rotm(3,1,2,theta[3]) %*% rotm(3,1,3,theta[2]) %*% rotm(3,2,3,theta[1])
  max(svd(Ss*pf$factor-pf$Q)$d) #>  4.84082e-16  but why pf$factor ??

[Package cwhmisc version 2.0.1 Index]