Dear Nick,
Thanks a lot for your help. Your advise works very good also for my data.
Your sincerely,
Catalina
Nick.Ellis@csiro.au wrote:
>
> 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
> >
|