On Sep 24, 1:19 pm, Luds <lud...@uvic.ca> wrote:
> I've been trying for a couple days now to write a Gaussian-smoothing
> algorithm to smooth a cube of (scalar) data with periodic boundary
> conditions (this is needed for my task since "structure" in the data
> that straddles an edge of the cube appears on two+ sides of the box).
> I've made it so far, but now can't seem to get around excessive For-
> loop's...
>
> For example, say the box of scalars values runs from (0,1) in x,y, and
> z, and has N^3 points. To smooth at point (x,y,z) in the box I
> generate a 3-D Gaussian with its centroid (mean) at point x,y,z:
>
> Gauss_field = rebin(periodic_gauss_func(X,[[sig],[x]]),N,N,N) * $
> rebin(reform(periodic_gauss_func(X,[[sig],[y]]),
> 1,N),N,N,N) * $
> rebin(reform(periodic_gauss_func(X,[[sig],[z]]),
> 1,1,N),N,N,N)
>
> where periodic_gauss_func is a 1-D Gaussian kernel function that wraps
> around the box edge, X=(0,1,...N-1)... sig=sigma. (i.e. this just does
> separate Gaussian smoothing along each direction and combines the
> result).
>
> Then the smoothed field at point (x,y,z) is something like
>
> Smoothed(x,y,z) = TOTAL(TOTAL(scalar_field*Gauss_field,1))
>
> What I can't figure out is an efficient way to do this for all (x,y,z)
> - for a N=1024^3 grid it takes a couple seconds to generate
> Gauss_field. Realistically, I'll have N=1024^3, so For-loops are
> pretty much useless(???), and memory is a bit of an issue too.
>
> Does anyone know of any "canned" routines to do this type of Gaussian
> smoothing? Or of an efficient way to convolve my 3D Gaussian field
> with my scalar field for all (x,y,z)? (I must stress that the Gaussian
> kernel must not be affected by, or truncated at, the box edge)
>
> Many thanks!!
>
> Aaron
Wouldn't the Fourier convolution theorem approach work here? FFT your
data cube, FFT your 3D Gaussian kernel, multiply them, and reverse FFT
them back out? You may need to judiciously use TEMPORARY and/or the /
OVERWRITE keyword if memory is an issue.
-Jeremy.
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
