#' align two 3D-pointclouds/meshes by their principal axes
#'
#' align two 3D-pointclouds/meshes by their principal axes
#' @param x matrix or mesh3d
#' @param y matrix or mesh3d, if missing, x will be centered by its centroid and aligned by its princial axis.
#' @param optim logical if TRUE, the RMSE between reference and target will be minimized testing all possible axes alignments and (if iterations > 0) followed by a rigid ICP procedure.
#' @param subsample integer: use subsampled points to decrease computation time of optimization.
#' @param iterations integer: number of iterations for optimization (the higher the more accurate but also more time consuming).
#' @param mc.cores use parallel processing to find best alignment to original shape.
#' @return rotated and translated version of x to the center and principal axes of y.
#' @details \code{x} and \code{y} will first be centered and aligned by their PC-axes. If \code{optim=TRUE},all possible 8 ordinations of PC-axes will be tested and the one with the smallest RMSE between the transformed version of \code{x} and the closest points on \code{y} will be used. Then the rotated version of \code{x} is translated to the original center of mass of \code{y}.
#' @examples
#' data(boneData)
#' blm1 <- pcAlign(boneLM[,,1],boneLM[,,2])
#' \dontrun{
#' require(rgl)
#' spheres3d(boneLM[,,1])#original position
#' spheres3d(blm1,col=2)#aligned configuration
#' spheres3d(boneLM[,,2],col=3)#target
#' }
#' @rdname pcAlign
#' @importFrom Rvcg vcgKDtree
#' @export
pcAlign <- function(x,y,optim=TRUE,subsample=NULL,iterations=10,mc.cores=2) UseMethod("pcAlign")
#' @rdname pcAlign
#' @export
pcAlign.matrix <- function(x, y,optim=TRUE,subsample=NULL,iterations=10, mc.cores=2) {
if (!missing(y)) {
if (inherits(y,"mesh3d"))
y <- vert2points(y)
if (.Platform$OS.type == "windows")
mc.cores <- 1
pca1 <- prcompfast(x, retx=FALSE)
pca2 <- prcompfast(y, retx=FALSE)
x <- scale(x, scale=F)
y <- scale(y, scale=F)
rotx <- pca1$rotation
roty <- pca2$rotation
chk <- det(rotx%*%roty)
if (chk < 0)
rotx[,3] <- -rotx[,3]
x <- x%*%rotx
y <- y%*%roty
rotlist <- list(
arot=getTrafoRotaxis(pt1=c(1,0,0),pt2=c(0,0,0),theta=pi),
brot=getTrafoRotaxis(pt1=c(0,1,0),pt2=c(0,0,0),theta=pi),
crot=getTrafoRotaxis(pt1=c(0,0,1),pt2=c(0,0,0),theta=pi))
tests <- as.matrix(expand.grid(c(1,0),c(1,0),c(1,0)))
tmpfun <- function(x,rotlist){
for (i in 1:3) {
if (x[i] == 0)
rotlist[[i]] <- diag(4)
}
return(rotlist)
}
subs <- 1:nrow(x)
if (!is.null(subsample)) {
subsample <- min(nrow(x)-1,subsample)
subs <- fastKmeans(x,k=subsample,iter.max = 100,threads=1)$selected
}
dists <- 1e10
fintrafo <- diag(4)
optifun <- function(i) {
#for (i in 1:8) {
rottmp <- tmpfun(tests[i,],rotlist)
trafotmp <- rottmp[[1]]%*%rottmp[[2]]%*%rottmp[[3]]
xtmp <- applyTransform(x,trafotmp)
if (iterations > 0) {
xtmp1 <- icpmat(xtmp,y,iterations=iterations,subsample=subsample)
trafoicp <- computeTransform(xtmp1,xtmp)
trafotmp <- trafoicp%*%trafotmp
disttmp <- mean(vcgKDtree(y,xtmp1[subs,],k=1)$dist^2)
} else {
disttmp <- mean(vcgKDtree(y,xtmp[subs,],k=1)$dist^2)
}
out <- list(dist=disttmp,trafo=trafotmp)
}
if (optim) {
checkit <- mclapply(1:8,optifun,mc.cores=mc.cores)
testdists <- sapply(checkit,function(x) x <- x$dist)
mindist <- which(testdists == min(testdists))[1]
fintrafo <- checkit[[mindist]]$trafo
}
x <- applyTransform(x,fintrafo)
x <- x%*%t(pca2$rotation)
x <- t(t(x)+pca2$center)
return(x)
} else {
x <- scale(x,scale=F)
rotms <- eigen(crossprod(x))$vectors
if (det(rotms) < 0)
rotms[,1] <- rotms[,1]*-1
x <- x%*%rotms
return(x)
}
}
#' @rdname pcAlign
#' @export
pcAlign.mesh3d <- function(x,y,optim=TRUE,subsample=NULL,iterations=10,mc.cores=2) {
xorig <- x
x <- vert2points(x)
tmpverts <- pcAlign(x,y,optim=optim,subsample=subsample,iterations=iterations,mc.cores=mc.cores)
xorig$vb[1:3,] <- t(tmpverts)
xorig <- vcgUpdateNormals(xorig)
return(xorig)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.