SpatialVx (version 1.0-2)

rigider: Rigid Transformation

Description

Find the optimal rigid transformation for a spatial field (e.g. an image).

Usage

rigider(x1, x0, p0, init = c(0, 0, 0), type = c("regular", "fast"),
    translate = TRUE, rotate = FALSE, loss, loss.args = NULL,
    interp = "bicubic", stages = TRUE,
    verbose = FALSE, ...)
# S3 method for rigided
plot(x, ...)
# S3 method for rigided
print(x, ...)
# S3 method for rigided
summary(object, ...)
rigidTransform(theta, p0, N, cen)

Value

A list object of class “rigided” is returned with components:

call

the function call.

translation.only

If stages argument is true, this part is the optimal translation before rotation.

rotate

optimal translation and rotation together, if stages argument is true.

initial

initial values used.

interp.method

same as input argument interp.

optim.args

optional arguments passed to nlminb.

loss, loss.args

same as input arguments.

par

optimal parameter values found.

value

value of loss function at optimal parameters.

x0, x1, p0

same as input arguments.

p1

transformed p0 coordinates.

x1.transformed

The field F(W(s)).

Arguments

x1, x0

matrices of same dimensions giving the forecast (or 1-energy) and observation (or 0-energy) fields, resp.

x, object

list object of class “rigided” as output by rigider.

N

(optional) the dimension of the fields (i.e., if x1 and x0 are n by m, then N is the product m * n).

cen

N by 2 matrix whoes rows are all the same giving the center of the field (used to subtract before determining rotations, etc.).

p0

N by 2 matrix giving the coordinates for the 0-energy (observed) field.

init

(optional) numeric vector of length equal to the number of parameters (e.g., 2 for translation only, 3 for both, and 1 for rotation only). If missing, then these will be estimated by taking the difference in centroids (translation) and the difference in orientation angles (rotation) as determined using image moments by way of imomenter.

theta

numeric vector of length 1, 2 or 3 (depending on whether you want to translate only (2), rotate only (1) or both (3)) giving the rigid transformation parameters.

type

character stating whether to optimize a loss function or just find the centroid (and possibly orientation angle) difference(s).

translate, rotate

logical, should the optimal translation/rotation be found?

loss

character naming a loss function (see details) to use in optimizing the rigid transformation (defaults to square error loss.

loss.args

named list giving any optional arguments to loss.

interp

character naming the 2-d interpolation method to use in calls to Fint2d. Must be one of “round” (default), “bilinear” or “bicubic”.

stages

logical. Should the optimal translation be found before finding both the optimal tranlsation and rotation?

verbose

logical. Should progress information be printed to the screen?

...

optional arguments to nlminb.

Author

Eric Gilleland

Details

A rigid transformation translates coordinates of values in a matrix and/or rotates them. That is, if (r, s) are coordinates in a field with center (c1, c2), then the rigid transformation with parameters (x, y) and theta is given by:

(r, s) + (x, y) + Phi ((r, s) - (c1, c2)),

where Phi is the matrix with first column given by (cos( theta ), - sin( theta)) and second column given by (sin( theta), cos( theta )).

The optimal transformation is found by way of numerical optimization using the nlminb function on the loss function given by loss. If no value is given for loss, then square error loss is assumed. In this case, the loss function is based on an assumption of Gaussian errors, but this assumption is only important if you try to make inferences based on this model, in which case you should probably think much harder about what you are doing. In particular, the default objective function, Q, is given by:

Q = - sum( ( F(W(s)) - O(s) )^2 / (2 * sigma^2) - (N / 2) * log( sigma^2 ),

where s are the coordinates, W(s) are the rigidly transformed coordinates, F(W(s)) is the value of the 1-enegy field (forecast) evaluated at W(s) (which is interpolated as the translations typically do not give integer translations), O(s) is the 0-energy (observed) field evaluated at coordinate s, and sigma^2 is the estimated variance of the error field. A good alternative is to use “QcorrRigid”, which calculates the correlation between F and O instead, and has been found by some to give better performance.

The function rigidTransform performs a rigid transform for given parameter values. It is intended as an internal function, but may be of use to some users.

See Also

nlminb, Fint2d

Examples

# Simple uninteresting example for the R robots.
x <- y <- matrix(0, 20, 40)

x[ 12:18, 2:3 ] <- 1

y[ 13:19, 5:6 ] <- 1

xycoords <- cbind(rep(1:20, 40), rep(1:40, each = 20))

tmp <- rigider(x1 = x, x0 = y, p0 = xycoords)
tmp
plot(tmp)

# Rotate a coordinate system.
data( "geom000" )

loc <- cbind(rep(1:601, 501), rep(1:501, each = 601))

# Rotate the coordinates by pi / 4.
th <- c(0, 0, pi / 4)
names(th) <- c("x", "y", "rotation")
cen <- colMeans(loc[ geom000 > 0, ])
loc2 <- rigidTransform(theta = th, p0 = loc, cen = cen)

geom101 <- Fint2d(X = geom000, Ws = loc2, s = loc, method = "round")

if (FALSE) {

image.plot(geom101)

# Try to find the optimal rigid transformation.
# First, allow a translation as well as rotation.

tmp <- rigider(x1 = geom101, x0 = geom000, p0 = loc,
    rotate = TRUE, verbose = TRUE)
tmp
plot(tmp)

# Now, only allow rotation, which does not work as
# well as one would hope.
tmp <- rigider(x1 = geom101, x0 = geom000, p0 = loc,
    translate = FALSE, rotate = TRUE, verbose = TRUE)
tmp
plot(tmp)

# Using correlation.
tmp <- rigider(x1 = geom101, x0 = geom000, p0 = loc,
    rotate = TRUE, loss = "QcorrRigid", verbose = TRUE)
tmp
summary(tmp)
plot(tmp)

##
## Examples from ICP phase 1.
##
## Geometric cases.
##

data( "geom001" )
data( "geom002" )
data( "geom003" )
data( "geom004" )
data( "geom005" )


tmp <- rigider(x1 = geom001, x0 = geom000, p0 = loc, verbose = TRUE)
tmp
plot(tmp)

tmp <- rigider(x1 = geom002, x0 = geom000, p0 = loc, verbose = TRUE)
tmp
plot(tmp)

tmp <- rigider(x1 = geom003, x0 = geom000, p0 = loc, verbose = TRUE)
tmp
plot(tmp)

tmp <- rigider(x1 = geom004, x0 = geom000, p0 = loc, verbose = TRUE)
tmp
plot(tmp)

# Note: Above is a scale error rather than a rotation, but can we
# approximate it with a rotation?
tmp <- rigider(x1 = geom004, x0 = geom000, p0 = loc, rotate = TRUE,
    verbose = TRUE)
tmp
plot(tmp)
# Ok, maybe need to give it better starting values?  Or, run it again
# with just the translation.

tmp <- rigider(x1 = geom005, x0 = geom000, p0 = loc, verbose = TRUE)
tmp
plot(tmp)


}