Convolution by a gaussian is equivalent to multiplication by a gaussian in
the Fourier domain. Since Fourier transforms can be done very efficiently
via the Fast Fourier Transform, gaussian smoothing can be achieved
efficiently by an fft, followed by a multiplication, followed by an inverse
fft. The examples below use 64=2^6 points for which fft is fast, but Splus
fft() is not restricted to power-of-2 grids.
#
# 1-d example of gaussian smoothing via FFT
#
one <- rep(0,64)
one[32] <- 1
ONE <- fft(one)
sigma <- .1
gauwin <- exp(-0.5*(-31:32)^2/32^2/sigma^2)/sqrt(2*pi)/sigma
gauwin <- gauwin[c(32:64,1:31)] # see DETAILS in ?fft for explanation of
this reordering
plot(gauwin)
ones.smooth <- fft(ONE*gauwin,inverse=T)/64
range(Im(ones.smooth))
plot(Re(ones.smooth))
# some experiments to check where is Nyquist frequency stored
one <- matrix(0,8,8)
one[] <- rep(1:2,c(8,8)) # Nyquist frequency data in x, constant in y
one
fft(one)/64
one[] <- 1:2 # Nyquist frequency data in y, constant in x
one
fft(one)/64
one[] <- 1
one[(row(one)+col(one))%%2==1] <- 2
one
fft(one)/64
#
# impulse response example
#
one <- matrix(0,64,64)
one[32,32] <- 1
ONE <- fft(one)
sigma <- .1
gauwin <- exp(-0.5*(-31:32)^2/32^2/sigma^2)/sqrt(2*pi)/sigma
gauwin <- gauwin[c(32:64,1:31)]
ONE <- sweep(ONE,1,gauwin,"*")
ONE <- sweep(ONE,2,gauwin,"*") # 2 sweeps are equivalent to multiplying by
2d gaussian
ones.smooth <- fft(ONE*gauwin,inverse=T)/64
image(one)
range(Im(ones.smooth))
image(Re(ones.smooth))
#
# noisy data example
#
one <- matrix(0,64,64)
one[32:33,32:33] <- 1
one <- one+rnorm(64*64)/10
ONE <- fft(one)
sigma <- .2
gauwin <- exp(-0.5*(-31:32)^2/32^2/sigma^2)/sqrt(2*pi)/sigma
gauwin <- gauwin[c(32:64,1:31)]
ONE <- sweep(ONE,1,gauwin,"*")
ONE <- sweep(ONE,2,gauwin,"*")
ones.smooth <- fft(ONE*gauwin,inverse=T)/64
image(one)
range(Im(ones.smooth))
image(Re(ones.smooth))
Nick Ellis
CSIRO Marine Research mailto:Nick.Ellis@csiro.au
PO Box 120 ph +61 (07) 3826 7260
Cleveland QLD 4163 fax +61 (07) 3826 7222
Australia http://www.marine.csiro.au
> -----Original Message-----
> From: Catalina Mesina [mailto:cmesina@cs.vu.nl]
> Sent: Tuesday, November 26, 2002 1:54 AM
> To: S-news list
> Subject: [S] Gaussian smoothing
>
>
>
> Dear all,
>
> I have tried a 3D smoothing of my data, using "loess", or
> smooth but I did not get the right smoothed image.
> Suppose that my data is a matrix "myim" 50x50 and at each
> point (x,y) in the grid I have an intensity value.
>
> myim<-matrix(rnorm(2500),50,50)
>
> x<-seq(1,50,1)
> y<-seq(1,50,1)
> grid<-expand.grid(x,y)
> myd<-data.frame(grid[,1],grid[,2],as.vector(myim))
> sm.myd<-fitted(loess(myd$X3~myd$X1*myd$X2,data=myd))
> image(matrix(myd$X3,50,50))
>
> Then what I would like to have is a gaussian smoothed image
> (matrix) of myim.
>
> I will appreciate every help answer from you.
>
> Catalina
> --------------------------------------------------------------------
> This message was distributed by s-news@lists.biostat.wustl.edu. To
> unsubscribe send e-mail to s-news-request@lists.biostat.wustl.edu with
> the BODY of the message: unsubscribe s-news
>
|