OF: Optical Flow Verification

In SpatialVx: Spatial Forecast Verification

Optical Flow Verification

Description

Perform verification using optical flow as described in Marzban and Sandgathe (2010).

Usage

OF(x, ...)

## Default S3 method:
OF(x, ..., xhat, W = 5, grads.diff = 1, center = TRUE, 
    cutoffpar = 4, verbose = FALSE)

## S3 method for class 'SpatialVx'
OF(x, ..., time.point = 1, obs = 1, model = 1, W = 5, grads.diff = 1,
    center = TRUE, cutoffpar = 4, verbose = FALSE)

## S3 method for class 'OF'
plot(x, ...)

## S3 method for class 'OF'
print(x, ...)

## S3 method for class 'OF'
hist(x, ...)

## S3 method for class 'OF'
summary(object, ...)

Arguments

x, xhat	Default: m by n matrices describing the verification and forecast fields, resp. The forecast field is considered the initial field that is morphed into the final (verification) field. OF.SpatialVx: list object of class “SpatialVx”. plot, hist and print methods: list object as returned by OF.
object	list object as returned by OF.
W	numeric/integer giving the window size (should be no smaller than 5).
grads.diff	1 or 2 describing whether to use first or second differences in finding the first derivative.
center	logical, should the fields be centered before performing the optical flow?
cutoffpar	numeric, set to NaN everything exceeding median +/- cutoffpar*sd.
verbose	logical, should progress information be printed to the screen?
time.point	numeric or character indicating which time point from the “SpatialVx” verification set to select for analysis.
obs, model	numeric indicating which observation/forecast model to select for the analysis.
...	For OF: optional arguments to the optim function (cannot be par, fn, gr or method). See details section for plot and hist method functions. Not used by the summary method function.

Details

Estimates the optical flow of the forecast field into the verification field. Letting I_o(x,y) and I_f(x,y) represent the intensities of each field at coordinate (x,y), the collection of pairs (dx, dy) is the optical flow field, where:

I_o(x,y) ~ I_f(x,y) + [partial(I_f) wrt x]*dx + [partial(I_f) wrt y]*dy.

The procedure follows that proposed by Lucas and Kanade (1981) whereby for some window, W, it is assumed that all dx (dy) are assumed constant, and least squares estimation is used to estimate dx and dy (see Marzban and Sandgathe, 2010 for more on this implementation). This function iteratively calls optflow for each window in the field.

The above formulation is linear in the parameters. Marzban and Sandgathe (2010) also introduce an additive error component, which leads to a nonlinear version of the above. Namely,

I_o(x,y) ~ I_f(x,y) + [partial(I_f) wrt x]*dx + [partial(I_f) wrt y]*dy + A(x,y).

See Marzban and Sandgathe for more details.

The plot method function can produce a figure like that of Fig. 1, 5, and 6 in Marzban and Sandgathe (2010) or with option full=TRUE, even more plots. Optional arguments that may be passed in via the ellipses include: full (logical, produce a figure analogous to Fig. 1, 5 and 6 from Marzban and Sandgathe (2010) (FALSE/default) or make more plots (TRUE)), scale (default is 1 or no scaling, any numeric value by which the fields are divided/scaled before plotting), of.scale (default is 1, factor by which display vectors can be magnified), of.step (plot OF vectors every of.step, default is 4), prop (default is 2, value for prop argument in the call to rose.diag from package CircStats), nbins (default is 40, number of bins to use in the call to rose.diag).

The hist method function produces a two-dimensional histogram like that of Fig. 3 and 7 in Marzban and Sandgathe (2010). It can also take various arguments passed via the ellipses. They include: xmin, xmax, ymin, ymax (lower and upper bounds for the histogram breaks in the x- (angle) and y- (magnitude/displacement error) directions, resp. Defaults to (0,360) and (0,4)), nbreaks (default is 100, the number of breaks to use).

The summary method mostly uses the stats function from package fields to summarize results of the errors, but also uses circ.summary from package CircStats for the angular errors.

Value

OF returns a list object of class “OF” with components:

data	list with components x and xhat containing the data.
data.name	character vector giving the names of the verification and forecast fields.
call	object of class “call” giving the original function call.
rows, cols	numeric vector giving the rows and columns used for finding the centers of windows. Needed by the plot and hist method functions.
err.add.lin	m by n matrix giving the linear additive errors (intensities).
err.mag.lin	m by n matrix giving the linear magnitude (displacement) errors.
err.ang.lin	m by n matrix giving the linear angular errors.
err.add.nlin, err.mag.nlin, err.ang.nlin	same as above but for nonlinear errors.
err.vc.lin, err.vr.lin, err.vc.nlin, err.vr.nlin	m by n matrices giving the x- and y- direction movements for the linear and nonlinear cases, resp.

The hist method function invisibly returns a list object of class “OF” that contains the same object that was passed in along with new components:

breaks	a list with components x and y giving the breaks in each direction
hist.vals	itself a list with components xb, yb (the number of breaks -1 used for each direction), and nb (the histogram values for each break)

The plot and summary mehtod functions do not return anything.

Author(s)

Caren Marzban, marzban “at” u.washington.edu, with modifications by Eric Gilleland

References

Lucas, B D. and Kanade, T. (1981) An iterative image registration technique with an application to stereo vision. Proc. Imaging Understanding Workshop, DARPA, 121–130.

Marzban, C. and Sandgathe, S. (2010) Optical flow for verification. Wea. Forecasting, 25, 1479–1494, doi:10.1175/2010WAF2222351.1.

See Also

optflow, optim, circ.summary

Examples

## Not run: 
data(hump)
initial <- hump$initial
final <- hump$final
look <- OF(final, xhat=initial, W=9, verbose=TRUE)
plot(look) # Compare with Fig. 1 in Marzban and Sandgathe (2010).
par(mfrow=c(1,1))
hist(look) # 2-d histogram.
plot(look, full=TRUE) # More plots.
summary(look)

# Another way to skin the cat.
hold <- make.SpatialVx( final, initial, field.type = "Bi-variate Gaussian",
    obs.name = "final", model.name = "initial" )

look2 <- OF(hold, W=9, verbose=TRUE)
plot(look2)
par(mfrow=c(1,1))
hist(look2)
plot(look2, full=TRUE)
summary(look2)

## End(Not run)

SpatialVx documentation built on May 29, 2024, 9:31 a.m.