Smooth discrete 2D array

Question

I need to smooth a 2D numpy array containing elevations at discrete steps. Here I have the z-value = 1 (elevation) in position (0, 0) however, it can be any value between 0 and 1. My array looks like this:

['1.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']

I tried applying a cosine kernel with scipy, and then using convolve2d. Like this:

import numpy as np
from scipy import signal

peak_array = np.array([
   [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

def elevation_grid():
    """Creates a smooth elevation-grid, from an array of peaks."""
    def gkern():
        """Returns a 2D Triangle kernel array."""
        # Decay is a cosine.
        # It is centered, so multiply by two, and add one
        # This way, we cover the entire grid, and assert
        # equal smoothing (due to an odd width)
        gkern1d_c = signal.cosine(peak_array.shape[0]*2 + 1)
        gkern1d_r = signal.cosine(peak_array.shape[1]*2 + 1)

        gkern2d = np.outer(gkern1d_c, gkern1d_r)
        return gkern2d

    kernel = gkern()
    grad = signal.convolve2d(peak_array, kernel, mode='same')

    # Normalize the grid
    grad -= np.amin(grad)
    grad /= np.amax(grad)

    return grad

def print_readable(array):
    """Prints the map to a human-readable format."""
    for row in range(0, array.shape[0]):
        # Round to two decimals
        r = ["%.2f" % array[col][row] for col in range(0, array.shape[1])]
        print(r)

smooth_array = elevation_grid()

print_readable(smooth_array)

Which results in an array looking like this:

['1.00', '0.99', '0.95', '0.90', '0.82', '0.72', '0.60', '0.47', '0.33', '0.18']
['0.99', '0.98', '0.94', '0.89', '0.81', '0.71', '0.60', '0.47', '0.33', '0.18']
['0.95', '0.94', '0.91', '0.85', '0.78', '0.68', '0.57', '0.45', '0.32', '0.17']
['0.90', '0.89', '0.85', '0.80', '0.73', '0.64', '0.54', '0.42', '0.29', '0.16']
['0.82', '0.81', '0.78', '0.73', '0.67', '0.59', '0.49', '0.38', '0.27', '0.14']
['0.72', '0.71', '0.68', '0.64', '0.59', '0.51', '0.43', '0.33', '0.23', '0.12']
['0.60', '0.60', '0.57', '0.54', '0.49', '0.43', '0.36', '0.28', '0.19', '0.09']
['0.47', '0.47', '0.45', '0.42', '0.38', '0.33', '0.28', '0.21', '0.14', '0.06']
['0.33', '0.33', '0.32', '0.29', '0.27', '0.23', '0.19', '0.14', '0.09', '0.03']
['0.18', '0.18', '0.17', '0.16', '0.14', '0.12', '0.09', '0.06', '0.03', '0.00']

Which is the expected result. However if i place peaks at each corner:

['1.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '1.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00']
['1.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '0.00', '1.00']

They add up in the middle resulting in a larger peak in the center, and i would expect the lowest value in the middle...

['0.00', '0.17', '0.30', '0.39', '0.44', '0.44', '0.39', '0.30', '0.17', '0.00']
['0.17', '0.36', '0.51', '0.61', '0.66', '0.66', '0.61', '0.51', '0.36', '0.17']
['0.30', '0.51', '0.67', '0.77', '0.83', '0.83', '0.77', '0.67', '0.51', '0.30']
['0.39', '0.61', '0.77', '0.89', '0.94', '0.94', '0.89', '0.77', '0.61', '0.39']
['0.44', '0.66', '0.83', '0.94', '1.00', '1.00', '0.94', '0.83', '0.66', '0.44']
['0.44', '0.66', '0.83', '0.94', '1.00', '1.00', '0.94', '0.83', '0.66', '0.44']
['0.39', '0.61', '0.77', '0.89', '0.94', '0.94', '0.89', '0.77', '0.61', '0.39']
['0.30', '0.51', '0.67', '0.77', '0.83', '0.83', '0.77', '0.67', '0.51', '0.30']
['0.17', '0.36', '0.51', '0.61', '0.66', '0.66', '0.61', '0.51', '0.36', '0.17']
['0.00', '0.17', '0.30', '0.39', '0.44', '0.44', '0.39', '0.30', '0.17', '0.00']

How can I smooth out my peaks? It needs to be efficient, since I have large arrays (up to 1000x1000).

Answer 1

Therefore, if k is your kernel, and * is convolution, and v1,..v4 input matrices like the first one you posted, you can sum them and obtain another matrix as w

w=v1+v2+v3+v4

You are applying convolution, which is a linear operation.

output=k * w = k * (v1+v2+v3+v4)= k * v1 + k2 * v2 ...

If you notice, you input with four 1s at the corner (Inp2) is basically the same as the sum of four of the initial matrix with only one 1 (Inp1).

You could build Inp2 by simply rotating by 90 degrees Inp1 and summing ( in numpy, transposes and fliprl). Therefore, you could just rotate the results of the first computation and add them up to find the correct final results.

It is very unlikely that the sum of the four rotated version of Inp1 is going to be different from the convolution apply to the Inp2.

If that is the case, it could be a numerical problem (if you are using those numbers it should not happen). Check well how you rescale the data with max and min. You are dividing by max, so make sure that is nonnegative, by more than an eps