1

I have a sparse 3D array of values. I am trying to turn each "point" into a fuzzy "sphere", by applying a Gaussian filter to the array.

I would like the original value at the point (x,y,z) to remain the same. I just want to create falloff values around this point... But applying the Gaussian filter changes the original (x,y,z) value as well.

I am currently doing this:

dataCube = scipy.ndimage.filters.gaussian_filter(dataCube, 3, truncate=8)

Is there a way for me to normalize this, or do something so that my original values are still in this new dataCube? I am not necessarily tied to using a Gaussian filter, if that is not the best approach.

Kalina
  • 5,504
  • 16
  • 64
  • 101
  • So you know how to filter but the problem is that the values change? Can't you just write the original values back into the new array? – MB-F Aug 17 '17 at 05:52
  • @kazemakase I want a 3d gaussian kernel which has a value of 1 at its peak – Kalina Aug 17 '17 at 15:40

1 Answers1

3

You can do this using a convolution with a kernel that has 1 as its central value, and a width smaller than the spacing between your data points.

1-d example:

import numpy as np
import scipy.signal
data = np.array([0,0,0,0,0,5,0,0,0,0,0])
kernel = np.array([0.5,1,0.5])
scipy.signal.convolve(data, kernel, mode="same")

gives

array([ 0. ,  0. ,  0. ,  0. ,  2.5,  5. ,  2.5,  0. ,  0. ,  0. ,  0. ])

Note that fftconvolve might be much faster for large arrays. You also have to specify what should happen at the boundaries of your array.

Update: 3-d example

import numpy as np
from scipy import signal

# first build the smoothing kernel
sigma = 1.0     # width of kernel
x = np.arange(-3,4,1)   # coordinate arrays -- make sure they contain 0!
y = np.arange(-3,4,1)
z = np.arange(-3,4,1)
xx, yy, zz = np.meshgrid(x,y,z)
kernel = np.exp(-(xx**2 + yy**2 + zz**2)/(2*sigma**2))

# apply to sample data
data = np.zeros((11,11,11))
data[5,5,5] = 5.
filtered = signal.convolve(data, kernel, mode="same")

# check output
print filtered[:,5,5]

gives

[ 0.          0.          0.05554498  0.67667642  3.0326533   5.          3.0326533
  0.67667642  0.05554498  0.          0.        ]
user8153
  • 4,049
  • 1
  • 9
  • 18
  • Thank you! That's a step in the right direction, but I'm still stuck. Any idea how to programmatically make a 3D, non-linear kernel? – Kalina Aug 16 '17 at 22:54
  • 1
    The second example in the [fftconvolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html#scipy.signal.fftconvolve) documentation shows how to build a two-dimensional Gaussian kernel. Alternatively, you can also stick to `gaussian_filter`, but multiply the output by `np.sqrt(2*np.pi * sigma**2)**3`. This should preserve the data values pretty well unless `sigma` is large and/or `truncate` is small. – user8153 Aug 17 '17 at 00:34
  • I see the example, thanks, but I don't know how to go from 2D to 3D, either. The example says `np.outer(signal.gaussian(70, 8), signal.gaussian(70, 8))`. So what do I do with that? – Kalina Aug 17 '17 at 15:17