Nothing
#' relax one specific landmark configuration against a reference
#'
#' relax one specific landmark configuration against a reference (e.g. a
#' sample mean)
#'
#'
#' @param lm k x 3 or k x 2 matrix containing landmark data to be slidden - or a triangular mesh of class "mesh3d". See details
#' @param reference k x 3 or k x 2 matrix containing landmark of the reference, or a mesh with the same amount of vertices as there are landmarks in \code{lm}.
#' @param SMvector A vector containing the row indices of (semi-) landmarks on the curve(s) that are
#' allowed to slide
#' @param outlines A vector (or if threre are several curves) a list of vectors
#' (containing the rowindices) of the (Semi-)landmarks forming the curve(s) in
#' the successive position on the curve - including the beginning and end
#' points, that are not allowed to slide.
#' @param surp integer vector containing the row indices of semi-landmarks positioned on surfaces.
#' @param sur.name character: containing the filename of the corresponding
#' surface.When specified, mesh has to be NULL. If \code{sur.name=NULL} and \code{mesh=NULL}, the tangent planes will be estimated from the point cloud.
#' @param mesh triangular mesh of class "mesh3d" loaded into the R workspace,
#' when specified, "sur.name" has to be NULL.
#' @param tol numeric: Threshold for convergence in the sliding proces. Full
#' Procrustes distance between actual result and previous iteration.
#' @param deselect Logical: if TRUE, the SMvector is interpreted as those
#' landmarks, that are not allowed to slide.
#' @param inc.check Logical: if TRUE, the program stops when convergence
#' criterion starts increasing and reports result from last iteration.
#' @param iterations integer: maximum amounts the algorithm runs - even when
#' 'tol' is not reached. When iterations=0, the algorithm runs until
#' convergence.
#' @param fixRepro logical: if \code{TRUE}, fix landmarks will also be
#' projected onto the surface. If you have landmarks not on the surface, select
#' \code{fixRepro=FALSE}
#' @param missing vector of integers, specifying row indices of missing (semi-)landmarks. They will be relaxed freely in 3D and not projected onto the target (works only for 2D data).
#' @param bending if TRUE, bending energy will be minimized, Procrustes distance otherwise (not suggested with large shape differences)
#' @param stepsize integer: dampening factor for the amount of sliding.
#' Useful to keep semi-landmarks from sliding too far off the surface.
#' The displacement is calculated as \eqn{\Upsilon = \Upsilon^0 + stepsize * UT}{Y = Y0 + stepsize * UT}.
#' Default is set to 1 for bending=TRUE and 0.5 for bending=FALSE.
#' @param use.lm indices specifying a subset of (semi-)landmarks to be used in the rotation step - only used if \code{bending=FALSE}.
#' @param silent logical: if TRUE, console output is suppressed.
#' @param ... additonal arguments - currently unused
#' @return returns kx3 matrix of slidden landmarks
#' @author Stefan Schlager
#' @seealso \code{\link{slider3d}}
#' @references Gunz, P., P. Mitteroecker, and F. L. Bookstein. 2005.
#' Semilandmarks in Three Dimensions, in Modern Morphometrics in Physical
#' Anthropology. Edited by D. E. Slice, pp. 73-98. New York: Kluwer
#' Academic/Plenum Publishers.
#' @details if \code{lm} is a surface mesh, all vertices will be treated as semilandmarks and a allowed to freely slide along the surface.
#' @examples
#'
#' require(rgl)
#' data(nose)
#' ### relax shornose against longnose
#'
#' # define fix landmarks
#' fix <- c(1:5,20:21)
#' # define surface patch by specifying row indices of matrices
#' # all except those defined as fix
#' surp <- c(1:dim(shortnose.lm)[1])[-fix]
#'
#' relax <- relaxLM(shortnose.lm,
#' longnose.lm, mesh=shortnose.mesh, iterations=1,
#' SMvector=fix, deselect=TRUE, surp=surp)
#'
#' ## example minimizing Procrustes distance when displacement is not
#' ## dampened by stepsize
#' relaxProcD <- relaxLM(shortnose.lm,
#' longnose.lm, mesh=shortnose.mesh, iterations=1,
#' SMvector=fix, deselect=TRUE, surp=c(1:623)[-fix],bending=FALSE,stepsize=1)
#'
#' \dontrun{
#' # visualize differences red=before and green=after sliding
#' deformGrid3d(shortnose.lm, relax, ngrid=0)
#'
#'
#' # visualize differences minimizing Procrusted distances red=before and green=after sliding
#'
#' deformGrid3d(shortnose.lm, relaxProcD, ngrid=0)
#' ## no smooth displacement, now let's check the distances:
#' rot2ref <- rotonto(relaxProcD,longnose.lm)
#' angle.calc(rot2ref$X,rot2ref$Y)
#' # 0.2492027 Procrustes distance between reference and slided shape
#' # (minimizing Procrustes distance)
#' rot2refBend <- rotonto(relax,longnose.lm)
#' angle.calc(rot2refBend$X,rot2refBend$Y)
#' # 0.2861322 Procrustes distance between reference and slided shape
#' # (minimizing bending energy)
#'
#' rot2refOrig <- rotonto(shortnose.lm,longnose.lm)
#' angle.calc(rot2refOrig$X,rot2refOrig$Y)
#' # 0.3014957 Procrustes distance between reference and original shape
#' ##result: while minimizing Procrustes distance, displacement is not
#' ##guaranteed to be smooth
#'
#' # add surface
#' wire3d(shortnose.mesh, col="white")
#'
#'
#' ## finally relax two meshes with corresponding vertices:
#'
#' mediumnose.mesh <- tps3d(shortnose.mesh,shortnose.lm, (shortnose.lm+longnose.lm)/2,threads=1)
#' ## we use Procrustes distance as criterion as bending energy is pretty slow because
#' ## of too many coordinates (more than 3000 is very unreasonable).
#' relaxMesh <- relaxLM(shortnose.mesh,mediumnose.mesh,iterations=2,bending=FALSE,stepsize=0.05)
#' }
#' @rdname relaxLM
#' @export
relaxLM <- function(lm,...)UseMethod("relaxLM")
#' @rdname relaxLM
#' @export
relaxLM.matrix <- function(lm,reference,SMvector,outlines=NULL,surp=NULL,sur.name=NULL,mesh=NULL,tol=1e-05,deselect=FALSE,inc.check=TRUE,iterations=0, fixRepro=TRUE, missing=NULL, bending=TRUE,stepsize=ifelse(bending,1,0.5),use.lm=NULL,silent=FALSE,...) {
nomesh <- FALSE
if (is.null(mesh) && is.null(sur.name))
nomesh <- TRUE
if(inherits(reference,"mesh3d"))
reference <- vert2points(reference)
k <- dim(lm)[1]
m <- dim(lm)[2]
weights <- NULL
if (!is.null(use.lm)) {
weights <- rep(0,nrow(lm))
weights[use.lm] <- 1
}
free <- NULL
p1 <- 10^12
lm.orig <- lm
reference <- scale(reference, scale=FALSE)
if (bending)
L <- CreateL(reference,output="Lsubk3")
if (deselect)
fixLM <- SMvector
else if (length(SMvector) < k)
fixLM <- c(1:k)[-SMvector]
else
fixRepro <- TRUE
if (iterations == 0)
iterations <- 1e10
if (m == 3) {
if (!silent)
if (!silent && !nomesh)
cat(paste("Points will be initially projected onto surfaces","\n","-------------------------------------------","\n"))
if (is.null(mesh) && !nomesh) {
tmp <- projRead(lm, sur.name)
} else if (!nomesh) {
tmp <- projRead(lm,mesh)
} else {
tmp <- vcgUpdateNormals(lm,silent=TRUE)
message("no surfaces specified - surface is approximated from point cloud")
}
vs <- vert2points(tmp)
vn <- t(tmp$normals[1:3,])
if (!fixRepro)# use original positions for fix landmarks
vs[fixLM,] <- lm.orig[fixLM,]
if (length(missing)) {
free <- missing
vs[missing,] <- lm.orig[missing,]
}
} else {
vs <- lm
}
count <- 1
while (p1 > tol && count <= iterations) {
lm_old <- vs
if (!silent)
cat(paste("Iteration",count,sep=" "),"..\n") # reports which Iteration is calculated
if (!bending) {
rot <- rotonto(reference,vs,reflection=FALSE,scale=TRUE,weights=weights,centerweight = TRUE)
vs <- rot$yrot
if (m == 3)
vn <- vn%*%rot$gamm
}
if (m == 3)
U <- .calcTang_U_s(vs,vn,SMvector=SMvector,outlines=outlines,surface=surp,deselect=deselect,free=free)
else
U <- .calcTang_U(vs,SMvector=SMvector,outlines=outlines,deselect=deselect)
if (bending)
dataslido <- calcGamma(U$Gamma0,L$Lsubk3,U$U,dims=m,stepsize=stepsize)
else {
dataslido <- calcProcDGamma(U$U,U$Gamma0,reference,dims=m,stepsize=stepsize)
dataslido <- rotreverse(dataslido,rot)
}
if (m == 3) {
if (is.null(mesh) && !nomesh) {
tmp <- projRead(dataslido, sur.name)
} else if (!nomesh){
tmp <- projRead(dataslido,mesh)
} else {
tmp <- vcgUpdateNormals(dataslido,silent=TRUE)
}
vs <- vert2points(tmp)
vn <- t(tmp$normals[1:3,])
if (!fixRepro)# use original positions for fix landmarks
vs[fixLM,] <- lm.orig[fixLM,]
if (length(missing))
vs[missing,] <- dataslido[missing,]
} else {
vs <- dataslido
}
p1_old <- p1
testproc <- rotonto(lm_old,vs)
p1 <- sum(diag(crossprod((testproc$X/cSize(testproc$X))-(testproc$Y/cSize(testproc$Y)))))### check for increasing convergence criterion ###
if (inc.check) {
if (p1 > p1_old) {
vs <- lm_old
if (!silent)
cat(paste("Distance between means starts increasing: value is ",p1, ".\n Result from last iteration step will be used. \n"))
p1 <- 0
count <- count+1
} else {
if (!silent)
cat(paste("squared distance between iterations:",p1,sep=" "),"\n","-------------------------------------------","\n")
count <- count+1
}
} else {
if (!silent)
cat(paste("squared distance between iterations:",p1,sep=" "),"\n","-------------------------------------------","\n")
count <- count+1
}
}
gc()
return(vs)
}
#' @rdname relaxLM
#' @export
relaxLM.mesh3d <- function(lm,reference,tol=1e-05,deselect=FALSE,inc.check=TRUE,iterations=0, fixRepro=TRUE, missing=NULL, bending=FALSE,stepsize=ifelse(bending,1,0.5),use.lm=NULL,silent=FALSE, ...){
lmtmp <- vert2points(lm)
mesh <- lm
SMvector <- surp <- 1:ncol(lm$vb)
sur.name <- NULL
outlines <- NULL
reltmp <- relaxLM(lmtmp,reference=reference,SMvector=SMvector,surp=surp,mesh=mesh,tol=tol,deselect=deselect,inc.check=inc.check,iterations=iterations, fixRepro=fixRepro, missing=missing, bending=bending,stepsize=stepsize,use.lm=use.lm,silent=silent)
lm$vb[1:3,] <- t(reltmp)
lm <- vcgUpdateNormals(lm)
return(lm)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.