Artifacts in convolution

Question

I am using a direct convolution algorithm to compute the convolution between this image:

and this kernel:

I am using the implementation in astropy for the direct convolution.

This results in the following convolution, leaving all settings (including boundary handling) to the defaults, i.e. astropy.convolution.convolve(image,kernel):

This convolution has some puzzling artifacts. In particular, there is a 'square' pattern at an offset of about 50 pixels from the edge. It seems to me that this is due to the extent of the kernel; even though the kernel is formally 249x249 in size, most information is clearly contained within a radius of about 100 pixels -- which means that we presumably run into trouble when the kernel is applied to the edges.

Which brings me to my questions:

1. Is this hypothesis correct - that it is indeed an edge problem?
2. How would I go about solving this? I don't know how to justify using different edge-handling (zero padding, interpolating, wrapping, ...) I'm sure different cases require different solutions, but I'm not sure how to decide on this...
3. Just... trying to understand the difference between using a direct algorithm and a FFT convolution. If the kernel and image are equally big, no zero padding is necessary for the FT convolution, no edge-effects will show up. For the direct method, you will inadvertently have some edge handling... so are there results even equivalent? Because in principle only their performance should differ, right?

Answer 1

Yes, this is an edge-effect issue which comes about because you have negative values in your kernel. As soon as the kernel is partially off the edge, the mean value of the kernel starts changing.

One solution would be to use boundary='fill' and fill_value=(mean of your image) or something along those lines. It might not completely remove these artifacts, but it ought to reduce them.

For the FFT convolution part of your question - the FFT convolution will do the same thing. However, edge padding is necessary for the FFT convolution, because otherwise the boundary will wrap. Not padding (e.g., convolve_fft(..., boundary='wrap')) will actually get rid of your artifacts, but it will do it in a way that may surprise you, since it will be averaging pixels from the right side of the image with the left side.

astropy's convolve and convolve_fft both will do the same thing given the same boundary conditions, but the naive fft convolution (i.e., conv = ifft(fft(im) * fft(kernel))) is equivalent to using boundary='wrap'.