#' rotates, translates and scales one matrix onto an other using Procrustes
#' fitting
#'
#' rotate a matrix onto an other without loosing information about the location
#' of the targetmatrix and reverse this transformations using rotreverse
#'
#'
#' @title rotates, translates and scales one matrix onto an other using Procrustes
#' fitting
#' @param x k x m matrix to be rotated onto (targetmatrix)
#' @param y k x m matrix which will be rotated (reference matrix)
#' @param scale logical: scale matrix to minimize sums of squares
#' @param signref logical: report if reflections were involved in the rotation
#' @param mat matrix on which the reverse transformations have to be applied
#' @param rot an object resulting from the former application of rotonto
#' @param reflection allow reflections.
#' @param weights vector of length k, containing weights for each landmark.
#' @param centerweight logical or vector of weights: if weights are defined and centerweigths=TRUE,
#' the matrix will be centered according to these weights instead of the
#' barycenter. If centerweight is a vector of length \code{nrow(x)}, the barycenter will be weighted accordingly.
#' @param ... currently not used
#' @return
#' \item{yrot }{rotated and translated matrix}
#' \item{Y }{centred and rotated reference matrix}
#' \item{X }{centred target matrix}
#' \item{trans }{vector between original position of target and centered
#' reference (during rotation process)}
#' \item{transy }{vector between original position of reference and
#' centered reference (during rotation process)}
#' \item{gamm }{rotation matrix}
#' \item{bet }{scaling factor applied}
#' \item{reflect }{if \code{reflect = 1}, reflections are involved in the
#' superimposition. Else, reflect = 0}
#' @author Stefan Schlager
#' @seealso \code{\link{rotmesh.onto}}
#' @references Lissitz, R. W., Schoenemann, P. H., & Lingoes, J. C. (1976). A
#' solution to the weighted Procrustes problem in which the transformation is
#' in agreement with the loss function. Psychometrika, 41,547-550.
#' @note all lines containing NA, or NaN are ignored in computing the transformation.
#' @examples
#'
#' if (require(shapes)) {
#' lims <- c(min(gorf.dat[,,1:2]),max(gorf.dat[,,1:2]))
#' rot <- rotonto(gorf.dat[,,1],gorf.dat[,,2]) ### rotate the second onto the first config
#' plot(rot$yrot,pch=19,xlim=lims,ylim=lims) ## view result
#' points(gorf.dat [,,2],pch=19,col=2) ## view original config
#' rev1 <- rotreverse(rot$yrot,rot)
#' points(rev1,cex=2) ### show inversion by larger circles around original configuration
#' }
#'
#' @export
rotonto <- function(x,y,scale=FALSE,signref=TRUE,reflection=TRUE,weights=NULL,centerweight=FALSE,...) {
xorig <- x
yorig <- y
reflect=0
k <- nrow(x)
m <- ncol(x)
##check for missing entries
xrows <- rowSums(x)
yrows <- rowSums(y)
xbad <- which(as.logical(is.na(xrows) + is.nan(xrows)))
ybad <- which(as.logical(is.na(yrows) + is.nan(yrows)))
bad <- sort(unique(c(xbad,ybad)))
docenter <- FALSE
if (length(centerweight) == 1) {
if (centerweight)
docenter <- TRUE
centerweight <- weights
} else {
docenter <- TRUE
}
if (length(bad)) {
message("some landmarks are missing and ignored for calculating the transform")
x <- x[-bad,,drop=FALSE]
y <- y[-bad,,drop=FALSE]
if (!is.null(weights))
weights <- weights[-bad]
if (!is.logical(centerweight[1]))
centerweight <- centerweight[-bad]
}
if (!is.null(weights))
weights <- weights/sum(weights)
if (!is.null(centerweight))
centerweight <- centerweight/sum(centerweight)
if (nrow(x) > 1) {
X <- scale(x, scale=FALSE)
Y <- scale(y, scale=FALSE)
} else {
X <- x
Y <- y
}
if (docenter && !is.null(centerweight)) {
xcent <- apply(X,2,weighted.mean,w=centerweight)
ycent <- apply(Y,2,weighted.mean,w=centerweight)
X <- scale(X,scale=F,center=xcent)
Y <- scale(Y,scale=F,center=ycent)
}
if (!is.null(weights)) {
Dn <- diag(weights)
X1 <- Dn%*%X
Y1 <- Dn%*%Y
XY <- crossprod(X1,Y1)
} else {
XY <- crossprod(X,Y)
}
sv1 <- svd(XY)
gamm <- tcrossprod(sv1$v,sv1$u)
if(sign(det(gamm))<1)
{ reflect <- 1
if (signref && reflection)
cat("reflection involved\n")
if (!reflection) {
u <- sv1$u
v <- sv1$v
chk1 <- Re(prod(eigen(v)$values))
chk2 <- Re(prod(eigen(u)$values))
if ((chk1 < 0) && (chk2 > 0)) {
v[, dim(v)[2]] <- v[, dim(v)[2]] * (-1)
gamm <- v %*% t(u)
}
if ((chk2 < 0) && (chk1 > 0)) {
u[, dim(u)[2]] <- u[, dim(u)[2]] * (-1)
gamm <- v %*% t(u)
}
}
}
trans <- x[1,]-X[1,]
transy <- y[1,]-Y[1,]
del <- sv1$d
## translate original y config to origin of centroid of potentially partial config and
## apply transforms
Yorig <- sweep(yorig,2,transy)
ctrace <- function(MAT) sum(diag(crossprod(MAT)))
if (scale) {
if (!is.null(weights))
bet <- sum(del)/ctrace(Y1)
else
bet <- sum(del)/ctrace(Y)
yrot <- bet*Yorig%*%gamm
} else {
bet <- 1
yrot <- Yorig%*%gamm
}
## scaled and rotated Y
Y <- yrot
## also translated to x
yrot <- t(t(yrot)+trans)
matlist <- list(yrot=yrot,Y=Y,X=X)
## move original x to (partial) centroid of x
matlist$X <- sweep(xorig,2,trans)
out <- list(yrot=matlist$yrot,Y=matlist$Y,X=matlist$X,trans=trans,transy=transy,gamm=gamm,bet=bet,reflect=reflect)
class(out) <- "rotonto"
return(out)
}
#' @rdname rotonto
#' @export
rotreverse <- function(mat,rot)UseMethod("rotreverse")
#' @rdname rotonto
#' @export
rotreverse.matrix <- function(mat,rot){
hmat <- solve(getTrafo4x4(rot))
out <-homg2mat(hmat%*%mat2homg(mat))
return(out)
}
#' @rdname rotonto
#' @export
rotreverse.mesh3d <- function(mat,rot)
{
x <- rotreverse(vert2points(mat),rot)
mat$vb[1:3,] <- t(x)
if (!is.null(mat$normals))
mat <- vcgUpdateNormals(mat)
return(mat)
}
#' get 4x4 Transformation matrix
#'
#' get 4x4 Transformation matrix
#' @param x object of class "rotonto"
#'
#' @return returns a 4x4 transformation matrix
#' @examples
#' data(boneData)
#' rot <- rotonto(boneLM[,,1],boneLM[,,2])
#' trafo <- getTrafo4x4(rot)
#' @rdname getTrafo4x4
#' @export
getTrafo4x4 <- function(x)UseMethod("getTrafo4x4")
#' @rdname getTrafo4x4
#' @export
getTrafo4x4.rotonto <- function(x) {
m <- ncol(x$gamm)
hgamm <- rbind(cbind(x$gamm,0),0);hgamm[m+1,m+1] <- 1
htrans <- diag(m+1);htrans[1:m,m+1] <- c(-x$transy)
htrans2 <- diag(m+1);htrans2[1:m,m+1] <- c(x$trans)
scale <- diag(m+1);diag(scale)[1:m] <- x$bet
hall <- htrans2%*%scale%*%t(hgamm)%*%htrans
return(hall)
}
mat2homg <- function(x) {
x <- rbind(t(x),1)
return(x)
}
homg2mat <- function(x) {
m <- nrow(x)
x <- t(x[1:(m-1),])
return(x)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.