image.smooth
Kernel smoother for irregular 2-d data
Description
Takes an image matrix and applies a kernel smoother to it. Missing values are handled using the Nadaraya/Watson normalization of the kernel.
Usage
## S3 method for class 'smooth'
image(x, wght = NULL, dx = 1, dy = 1,
kernel.function = double.exp,
theta = 1, grid = NULL, tol = 1e-08, xwidth = NULL, ywidth = NULL,
weights = NULL,...)
setup.image.smooth(nrow = 64, ncol = 64, dx = 1, dy = 1,
kernel.function = double.exp,
theta = 1, xwidth = nrow * dx, ywidth = ncol * dx, lambda=NULL, ...)
Arguments
x |
A matrix image. Missing values can be indicated by NAs. |
wght |
FFT of smoothing kernel. If this is NULL the default is to compute this object. |
grid |
A list with x and y components. Each are equally spaced and define the rectangular. ( see grid.list) |
dx |
Grid spacing in x direction |
dy |
Grid spacing in x direction |
kernel.function |
An R function that takes as its argument the squared distance between two points divided by the bandwidth. The default is exp( -abs(x)) yielding a normal kernel |
theta |
the bandwidth or scale parameter. |
xwidth |
Amount of zero padding in horizontal dimension in units of the grid spacing. If NULL the default value is equal to the width of the image the most conservative value but possibly inefficient for computation. Set this equal to zero to get periodic wrapping of the smoother. This is useful to smooth a Mercator map projection. |
ywidth |
Same as xwidth but for the vertical dimension. |
weights |
Weights to apply when smoothing. |
tol |
Tolerance for the weights of the N-W kernel. This avoids kernel estimates that are “far” away from data. Grid points with weights less than tol are set to NA. |
nrow |
X dimension of image in setting up smoother weights |
ncol |
Y dimension of image |
lambda |
Smoothing parameter if smoother is interpreted in a spline-like way. |
... |
Other arguments to be passed to the kernel function |
Details
The function works by taking convolutions using an FFT. The missing pixels are taken into account and the kernel smoothing is correctly normalized for the edge effects following the classical Nadaraya-Watson estimator. For this reason the kernel doe snot have to be a desity as it is automatically normalized when the kernel weight function is found for the data. If the kernel has limited support then the width arguments can be set to reduce the amount of computation. (See example below.) For multiple smoothing compute the fft of the kernel just once using setup.image.smooth and pass this as the wght argument to image.smooth. this will save an FFT in computations.
Value
The smoothed image in R image format. ( A list with components x, y and z.) setup.image.smooth returns a list with components W a matrix being the FFT of the kernel, dx, dy, xwidth and ywidth.
See Also
as.image, sim.rf, image.plot
Examples
# first convert precip data to the 128X128 discretized image format ( with
# missing values to indicate where data is not observed)
#
out<- as.image( RMprecip$y, x= RMprecip$x, nrow=128, ncol=128)
# out$z is the image matrix
dx<- out$x[2]- out$x[1]
dy<- out$y[2] - out$y[1]
#
# grid scale in degrees and choose kernel bandwidth to be .25 degrees.
look<- image.smooth( out, theta= .25)
image.plot(look)
points( RMprecip$x)
US( add=TRUE, col="grey", lwd=2)
# to save on computation, decrease the padding with zeroes
# only pad 32 grid points around the margins ofthe image.
look<- image.smooth(out$z, dx=dx, dy=dy, theta= .25, xwidth=32*dx,ywidth=32*dy)
# the range of these data is ~ 10 degrees and so
# with a padding of 32 grid points 32*( 10/128) = 2.5
# about 10 standard deviations of the normal kernel so there is still
# lots of room for padding
# a minimal choice might be xwidth = 4*(.25)= 1 4 SD for the normal kernel
# creating weighting object outside the call
# this is useful when one wants to smooth different data sets but on the
# same grid with the same kernel function
#
#
# random fields from smoothing white noise with this filter.
#
set.seed(123)
test.image<- matrix( rnorm(128**2),128,128)
wght<- setup.image.smooth( nrow=128, ncol=128, dx=dx, dy=dy,
theta=.25, xwidth=2.5, ywidth=2.5)
#
look<- image.smooth( test.image, dx=dx, dy=dy, wght)
# NOTE: this is the same as using
#
# image.smooth( test.image , 128,128), xwidth=2.5,
# ywidth=2.5, dx=dx,dy=dy, theta=.25)
#
# but the call to image.smooth is faster because fft of kernel
# has been precomputed.
# periodic smoothing in the horizontal dimension
look<- image.smooth( test.image , xwidth=1.5,
ywidth=2.5, dx=dx,dy=dy, theta=1.5)
look2<- image.smooth( test.image , xwidth=0,
ywidth=2.5, dx=dx,dy=dy, theta=1.5)
# compare these two
set.panel( 1,2)
image.plot( look)
title("free boundaries")
image.plot( look2) # look for periodic continuity at edges!
title("periodic boundary in horizontal")
set.panel(1,1)
Results
R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(fields)
Loading required package: spam
Spam version 0.29-3 (2013-04-23) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
Attaching package: 'spam'
The following objects are masked from 'package:base':
backsolve, forwardsolve
Loading required package: maps
> png(filename="/home/oogasawa/snapshot/RGM3/R_CC/result/fields/image.smooth.Rd_%03d_medium.png", width=480, height=480)
> ### Name: image.smooth
> ### Title: Kernel smoother for irregular 2-d data
> ### Aliases: image.smooth setup.image.smooth
> ### Keywords: smooth
>
> ### ** Examples
>
> # first convert precip data to the 128X128 discretized image format ( with
> # missing values to indicate where data is not observed)
> #
> out<- as.image( RMprecip$y, x= RMprecip$x, nrow=128, ncol=128)
> # out$z is the image matrix
>
> dx<- out$x[2]- out$x[1]
> dy<- out$y[2] - out$y[1]
>
> #
> # grid scale in degrees and choose kernel bandwidth to be .25 degrees.
>
> look<- image.smooth( out, theta= .25)
>
> image.plot(look)
> points( RMprecip$x)
> US( add=TRUE, col="grey", lwd=2)
>
> # to save on computation, decrease the padding with zeroes
> # only pad 32 grid points around the margins ofthe image.
>
> look<- image.smooth(out$z, dx=dx, dy=dy, theta= .25, xwidth=32*dx,ywidth=32*dy)
>
> # the range of these data is ~ 10 degrees and so
> # with a padding of 32 grid points 32*( 10/128) = 2.5
> # about 10 standard deviations of the normal kernel so there is still
> # lots of room for padding
> # a minimal choice might be xwidth = 4*(.25)= 1 4 SD for the normal kernel
> # creating weighting object outside the call
> # this is useful when one wants to smooth different data sets but on the
> # same grid with the same kernel function
> #
>
> #
> # random fields from smoothing white noise with this filter.
> #
> set.seed(123)
> test.image<- matrix( rnorm(128**2),128,128)
>
>
> wght<- setup.image.smooth( nrow=128, ncol=128, dx=dx, dy=dy,
+ theta=.25, xwidth=2.5, ywidth=2.5)
> #
> look<- image.smooth( test.image, dx=dx, dy=dy, wght)
>
> # NOTE: this is the same as using
> #
> # image.smooth( test.image , 128,128), xwidth=2.5,
> # ywidth=2.5, dx=dx,dy=dy, theta=.25)
> #
> # but the call to image.smooth is faster because fft of kernel
> # has been precomputed.
>
>
>
> # periodic smoothing in the horizontal dimension
>
> look<- image.smooth( test.image , xwidth=1.5,
+ ywidth=2.5, dx=dx,dy=dy, theta=1.5)
> look2<- image.smooth( test.image , xwidth=0,
+ ywidth=2.5, dx=dx,dy=dy, theta=1.5)
> # compare these two
> set.panel( 1,2)
plot window will lay out plots in a 1 by 2 matrix
> image.plot( look)
> title("free boundaries")
> image.plot( look2) # look for periodic continuity at edges!
> title("periodic boundary in horizontal")
> set.panel(1,1)
plot window will lay out plots in a 1 by 1 matrix
>
>
>
>
>
>
> dev.off()
null device
1
>
  |