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
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
