Centered finite difference approximation for edge enhancement I implemented a centered finite difference approximation for edge enhancement in images. The difference approximation has a step size h of the partial derivatives in the x- and y-directions of a 2D-image. However, the function that calculates the difference approximation is slow on big images. 

Here is the code:

import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib
from skimage import io
import numpy as np    

# read image as a numpy array
sample_img = io.imread('samlpe.png')

def findiff(image, stepsize):
    # Padding to handle the boundaries when incrementing with step size
    padded_img = np.pad(image, pad_width=1, mode='constant', constant_values=0)
    img = padded_img.astype(float)
    M, N = img.shape # Get the dimensions of the image.
    pdfx = np.zeros(shape=(M, N))
    pdfy = np.zeros(shape=(M, N))
    img_lst = []
    for x in range(M-1):
        for y in range(N-1):
            pdfx[x][y] = (img[x][y + stepsize] - img[x][y - stepsize]) / (2 * stepsize) # partial with relation to x
            pdfy[x][y] = (img[x + stepsize][y] - img[x - stepsize][y]) / (2 * stepsize) # partial with relation to y

    img_lst.append(pdfx)
    img_lst.append(pdfy)
    return img_lst


fig, img_lst = plt.subplots(1, 2) 
imgs = findiff(sample_img, 1)
img_lst[0].imshow(imgs[0], cmap="gray")
img_lst[1].imshow(imgs[1], cmap="gray")
img_lst[0].set_title('pdfx')
img_lst[1].set_title('pdfy')
plt.show()


This solution works for step size=1. How can I improve this code? Can I avoid the nested loops? How can I generalize the function so that it can handle different step sizes?

I'm new to Python by the way.
I'm not yet very comfortable with Python style and so on, so I'm going to focus only on speeding up the code.

First of all, the indexing

img[x] extracts a line from the image. Then img[x][y] extracts a pixel from that line. In contrast, img[x,y] directly extracts a single pixel from the image. Thus, we can rewrite your indexing expressions reduces the amount of work to be done:

pdfx[x,y] = (img[x, y + stepsize] - img[x, y - stepsize]) / (2*stepsize) # partial with relation to x
pdfy[x,y] = (img[x + stepsize, y] - img[x - stepsize, y]) / (2*stepsize) # partial with relation to y


Let's time the difference:

import timeit
timeit.timeit("imgs = findiff(sample_img, 1)", number=30, globals=globals())


The original code gives 3.46 seconds, the code with changed indexing gives 2.37 seconds. That is a significant difference!

Next, let's remove the double loop

You should always leverage the vectorized operations in Numpy. In this case it is not difficult. The lines of code become longer, but with practice one immediately sees what is being computed. This is the modified function:

def findiff(image, stepsize):
    # Padding to handle the boundaries when incrementing with step size
    padded_img = np.pad(image, pad_width=stepsize, mode='constant', constant_values=0)
    img = padded_img.astype(float)
    pdfx = (img[stepsize:-stepsize, 2*stepsize:] - img[stepsize:-stepsize, 0:-2*stepsize]) / (2*stepsize) # partial with relation to x
    pdfy = (img[2*stepsize:, stepsize:-stepsize] - img[0:-2*stepsize, stepsize:-stepsize]) / (2*stepsize) # partial with relation to x
    img_lst = []
    img_lst.append(pdfx)
    img_lst.append(pdfy)
    return img_lst


Note the computation of pdfx and pdfy. Here we do exactly the same computations as before, but we extract an array containing all the elements that appear on the left hand side of the minus operator as you go through the original loop, and another array containing all the right hand side elements. Note that img[stepsize:-stepsize, stepsize:-stepsize] is the original image before padding (since we add stepsize pixels to each side). On the LHS we want to index pixels stepsize elements to the right, that is img[stepsize:-stepsize, (stepsize:-stepsize)+stepsize] == img[stepsize:-stepsize, 2*stepsize:]. On the RHS we want to index pixels stepsize elements to the left which we obtain in a similar way.

The modified function runs in 0.021 seconds. That is two orders of magnitude faster than the loop code!

Note that the output images now are of the same size as the input image, not of the size of the padded input image as before. This is more useful in general. But if it is desired to output images with padding then it is possible to instead do:

pdfx = np.zeros(shape=(M, N))
pdfx[stepsize:-stepsize, stepsize:-stepsize] = ...


Another change in the code above is the use of pad_width=stepsize in the call to np.pad, so that the indexing works correctly for step sizes larger than 1.