[SciPy-user] filling array without loop...

[Anne Archibald]

On 22/04/07, fred <fredmfp at gmail.com> wrote:

> Each array cell is a convolution that I wrote as a scalar product
> between KW, a weights filter vector, and input_data,
> the "raw" data (3D); output_data is the filtered response given by the
> filter.

First: it is almost certainly possible to write this without a loop.
But the algorithm will still be O(n**3*m**3). You can rewrite
convolutions using the Fourier transform to get an algorithm that is
O((n+m)**3*log(n+m)**3), which will be much faster. Doing it for a
three-dimensional matrix will be a bit messy.

One problem with convolutions is you need to specify what happens at
the boundaries. You can use periodic boundaries, you can pad the input
with zeros, or you can discard everything in the output that's too
close to the edge.

For your particular problem, I'd say, use the FFT approach. So:
* Pad your arrays to the same size.
** The FFT method will give you periodic boundary conditions, that is,
if your convolution blurs your image, the FFT will blur the left side
into the right side. If the blur is only 2n pixels wide, then trimming
off the n pixels on either side (as you do now) will get rid of this.
** For comprehensibility you may want to center your "blur" array at
index zero (negative indices go at the end of the array). Not doing
this will shift the result array.
* Take numpy.fft.rfftn of both arrays.
* Multiply the resulting (complex) arrays.
* Take numpy.fft.irfftn of the result.

A good test case is using as your "blur" array something like
[1,0,0,0] (should leave the array unchanged), [0,1,0,0] (should shift
the array over by one), or [1,0,1,0] (should add a shifted copy to the
original), or rather, their three-dimensional equivalents.

For more on what you can do with the FFT in terms of convolution (for
example, if your convolution kernel is small you can make this run
faster), look at Numerical Recipes in C (which is available for free
on the Web).

Good luck,
Anne