Nothing
landmarkreg <- function(fdobj, ximarks, x0marks=xmeanmarks,
WfdPar=NULL, monwrd=FALSE, ylambda=1e-10,
returnMatrix=FALSE)
{
# Arguments:
# FDOBJ ... functional data object for curves to be registered
# XIMARKS ... N by NL array of times of interior landmarks for
# each observed curve
# XOMARKS ... vector of length NL of times of interior landmarks for
# target curve
# WFDPAR ... a functional parameter object defining a warping function.
# If NULL, registration is done using linear interpolation
# of lamdmark times in XIMARKS plotted against corresponding
# target times in X0MARKS.
# MONWRD ... If TRUE, warping functions are estimated by monotone smoothing,
# otherwise by regular smoothing. The latter is faster, but
# not guaranteed to produce a strictly monotone warping
# function. If MONWRD is 0 and an error message results
# indicating nonmonotonicity, rerun with MONWRD = 1.
# Default: TRUE
# YLAMBDA ... smoothing parameter to be used in computing the registered
# functions. For high dimensional bases, local wiggles may be
# found in the registered functions or its derivatives that are
# not seen in the unregistered functions. In this event, this
# parameter should be increased.
# Returns:
# FDREG ... a functional data object for the registered curves
# WARPFD ... a functional data object for the warping functions
# WFD ... a functional data object for the W functions defining the
# warping functions
# RETURNMATRIX ... If False, a matrix in sparse storage model can be returned
# from a call to function BsplineS. See this function for
# enabling this option.
# Last modified 14 November 2012 by Jim Ramsay
# check FDOBJ
if (!(inherits(fdobj, "fd"))) stop(
"Argument fdobj not a functional data object.")
# extract information from curve functional data object and its basis
coef <- fdobj$coefs
coefd <- dim(coef)
ndim <- length(coefd)
ncurve <- coefd[2]
if (ndim > 2) {
nvar <- coefd[3]
} else {
nvar <- 1
}
basisobj <- fdobj$basis
nbasis <- basisobj$nbasis
rangeval <- basisobj$rangeval
fdParobj <- fdPar(basisobj, 2, ylambda)
# check landmarks
if (is.vector(ximarks) | is.data.frame(ximarks) ) ximarks = as.matrix(ximarks)
ximarksd <- dim(ximarks)
if (ximarksd[1] != ncurve) stop(
"Number of rows of XIMARKS is incorrect.")
if (any(ximarks <= rangeval[1]) || any(ximarks >= rangeval[2])) stop(
"Some landmark values are not within the range.")
nlandm <- dim(ximarks)[2]
xmeanmarks <- apply(ximarks,2,mean)
if (length(x0marks) != nlandm) stop(
"Number of target landmarks not equal to number of curve landmarks.")
# set up default WfdPar
if (is.null(WfdPar)) {
basisobj <- fdobj$basis
rangex <- basisobj$rangeval
wnbasis <- length(x0marks) + 2
wbasis <- create.bspline.basis(rangex, wnbasis)
WfdParobj <- fdPar(wbasis)
}
# check WFDPAR
WfdPar <- fdParcheck(WfdPar)
# set up WFD0 and WBASIS
Wfd0 <- WfdPar$fd
wLfd <- WfdPar$Lfd
wbasis <- Wfd0$basis
# set up WLAMBDA
wlambda <- WfdPar$lambda
# check landmark target values
wrange <- wbasis$rangeval
if (any(rangeval != wrange)) stop(
"Ranges for FD and WFDPAR do not match.")
# set up analysis
n <- min(c(101,10*nbasis))
x <- seq(rangeval[1],rangeval[2],length=n)
y <- eval.fd(x, fdobj, returnMatrix=returnMatrix)
yregmat <- y
hfunmat <- matrix(0,n,ncurve)
wlambda <- max(wlambda,1e-10)
xval <- c(rangeval[1],x0marks,rangeval[2])
nwbasis <- wbasis$nbasis
Wcoef <- matrix(0,nwbasis,ncurve)
nval <- length(xval)
wval <- rep(1,nval)
# --------------------------------------------------------------------
# Iterate through curves to register
# --------------------------------------------------------------------
cat("Progress: Each dot is a curve\n")
for (icurve in 1:ncurve) {
cat(".")
# set up landmark times for this curve
yval <- c(rangeval[1],ximarks[icurve,],rangeval[2])
# smooth relation between this curve"s values and target"s values
if (monwrd) {
# use monotone smoother
Wfds <- smooth.morph(xval, yval, WfdPar)
Wfd <- Wfds$Wfdobj
h <- monfn(x, Wfd, returnMatrix=returnMatrix)
b <- (rangeval[2]-rangeval[1])/(h[n]-h[1])
a <- rangeval[1] - b*h[1]
h <- a + b*h
h[c(1,n)] <- rangeval
wcoefi <- Wfd$coef
Wcoef[,icurve] <- wcoefi
} else {
# use unconstrained smoother
warpfd <- smooth.basis(xval, yval, WfdPar, wval)$fd
# set up warping function by evaluating at sampling values
h <- as.vector(eval.fd(x, warpfd, returnMatrix=returnMatrix))
b <- (rangeval[2]-rangeval[1])/(h[n]-h[1])
a <- rangeval[1] - b*h[1]
h <- a + b*h
h[c(1,n)] <- rangeval
# check for monotonicity
deltah <- diff(h)
if (any(deltah <= 0)) stop(
paste("Non-increasing warping function estimated for curve",icurve,
" Try setting MONWRD to TRUE."))
wcoefi <- warpfd$coef
Wcoef[,icurve] <- wcoefi
}
hfunmat[,icurve] <- h
# compute h-inverse in order to register curves
if (monwrd) {
wcoef <- Wfd$coefs
Wfdinv <- fd(-wcoef,wbasis)
WfdParinv <- fdPar(Wfdinv, wLfd, wlambda)
Wfdinv <- smooth.morph(h, x, WfdParinv)$Wfdobj
hinv <- monfn(x, Wfdinv, returnMatrix=returnMatrix)
b <- (rangeval[2]-rangeval[1])/(hinv[n]-hinv[1])
a <- rangeval[1] - b*hinv[1]
hinv <- a + b*hinv
hinv[c(1,n)] <- rangeval
} else {
hinvfd <- smooth.basis(h, x, WfdPar)$fd
hinv <- as.vector(eval.fd(x, hinvfd, returnMatrix=returnMatrix))
b <- (rangeval[2]-rangeval[1])/(hinv[n]-hinv[1])
a <- rangeval[1] - b*hinv[1]
hinv <- a + b*hinv
hinv[c(1,n)] <- rangeval
deltahinv <- diff(hinv)
if (any(deltahinv <= 0)) stop(
paste("Non-increasing warping function estimated for curve",icurve))
}
# compute registered curves
if (length(dim(coef)) == 2) {
# single variable case
yregfd <- smooth.basis(hinv, y[,icurve], fdParobj)$fd
yregmat[,icurve] <- eval.fd(x, yregfd, 0, returnMatrix=returnMatrix)
}
if (length(dim(coef)) == 3) {
# multiple variable case
for (ivar in 1:nvar) {
# evaluate curve as a function of h at sampling points
yregfd <- smooth.basis(hinv, y[,icurve,ivar], fdParobj)$fd
yregmat[,icurve,ivar] <- eval.fd(x, yregfd,
returnMatrix=returnMatrix)
}
}
}
cat("\n")
# create functional data objects for the registered curves
fdParobj <- fdPar(basisobj, 2, ylambda)
regfdobj <- smooth.basis(x, yregmat, fdParobj)$fd
regnames <- fdobj$fdnames
names(regnames)[3] <- paste("Registered",names(regnames)[3])
regfdobj$fdnames <- regnames
# create functional data objects for the warping functions
warpfdobj <- smooth.basis(x, hfunmat, fdParobj)$fd
warpfdnames <- fdobj$fdnames
names(warpfdnames)[3] <- paste("Warped",names(regnames)[1])
warpfdobj$fdnames <- warpfdnames
Wfd <- fd(Wcoef, wbasis)
return( list("regfd" = regfdobj, "warpfd" = warpfdobj, "Wfd" = Wfd) )
}
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.