Nothing
Q <- function( p, p0, s, Im0, Im1, imethod, B, beta, Cmat = NULL, do.gr = FALSE, dF = NULL, ... ) {
# 'p' are the 1-energy (p1) column-stacked locations: nc * d
# numeric vector (although, for now, d = 2). Last value is
# the nuisance parameter, sigma_epsilon.
# 'p0' are the 0-energy control locations: nc X d matrix.
# 's' are the N X d full set of locations (N >= nc).
# 'B' the Tps "B" matrix (calculated by 'warpTpsMatrices').
# 'beta' single numeric giving the penalty term
# (chosen a priori by user).
# 'Im0' k X m (where k * m = N) matrix giving the 0-energy image.
# 'Im1' k X m matrix giving the 1-energy image.
# 'imethod' character string naming the interpolation method.
# 'Cmat' nc X nc precision matrix for penalty function.
# 'do.gr' logical, should the gradients also be returned.
# '...' not used.
np <- length( p )
nc <- nrow( p0 )
# p1 <- cbind( p[ 1:nc ], p[ (nc + 1):( np - 1 ) ] )
p1 <- cbind( p[ 1:nc ], p[ (nc + 1):np ] )
# Technically, there is a nuisance parameter, sigma, but it is not really
# needed if not using a true likelihood (cf. 'Qgauss'). Leaving it here
# for posterity in case it can be useful later.
sigma <- sigma2 <- 1
# Intensity part of the loss function.
wpts <- warpTps( p1 = p1, B = B )
Im1.def <- Fint2d( Im1, Ws = wpts, s = s, method = imethod, derivs = do.gr )
if( do.gr ) {
# Calculate components of the gradient function.
# Gradient associated with the interpolation part.
dX <- Im1.def$dx
dY <- Im1.def$dy
# Convert Im1.def back to its form for calculating the
# objective function.
Im1.def <- Im1.def$xy
# Gradient associated with the pixels themselves.
# This part is done before warping!
if( is.null( dF ) ) dF <- dF( Im1 )
} # end of if 'do.gr' stmt.
ImD <- Im1.def - Im0
ImD2 <- ImD^2
if( !do.gr ) {
ll <- sum( colSums( ImD2 ) ) / sigma2
if( is.na( ll ) ) ll <- sum( colSums( ImD2, na.rm = TRUE ), na.rm = TRUE ) / sigma2
} # end of if calculate the objective function stmt.
# Penalty part of the loss function.
xdelta <- p1[,1, drop = FALSE] - p0[,1, drop = FALSE]
ydelta <- p1[,2, drop = FALSE] - p0[,2, drop = FALSE]
if( do.gr ) {
dLx = dF$dFx * ( ImD / ( sigma2 ) ) * dX
dLy = dF$dFy * ( ImD / ( sigma2 ) ) * dY
dL = c( t( B ) %*% cbind( c( dLx ), c( dLy ) ) )
} # end of if 'do.gr' stmt.
if( beta > 0 ) {
if( is.null( Cmat ) ) stop( "Q: must specify Cmat argument when beta > 0." )
if( !do.gr ) {
pen <- beta * ( t(xdelta) %*% Cmat %*% xdelta +
t(ydelta) %*% Cmat %*% ydelta ) / 2 # + nc * log( beta )
} else dpen <- beta * c( t(xdelta) %*% Cmat, t(ydelta) %*% Cmat )
} else if( beta == 0 ) pen <- dpen <- 0
else stop("Q: beta must be non-negative.")
if( !do.gr ) res <- ll + pen
else res <- dL + dpen
return( res )
} # end of 'Q' function.
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.