Median of derivative in X axis of an image

Question

I computed derivatives using different methods such as :

1. convolution with an array [[-1, 1]].
2. Using the fourier theorem by computing DFT of the image and the array mentioned above, multiplying them and performing IDFT.
3. Directly through the derivative formula (Computing Fourier, multiplying by index and a constant and computing the inverse).

All methods seem to work almost identically, but have slight differences.

An explanation why they end up with slightly different results would be appreciated.

After computing those I started playing with the result to learn about it, and I found out something that confused me:

The main thing that baffles me is that when I try computing the median of this derivative, its ALWAYS 0.0.

Why is that?

I added the code I used to compute this (the first method at least) because maybe I'm doing something wrong.

from scipy.signal import convolve2d

im = sl.read_image(r'C:\Users\ahhal\Desktop\Essentials\Uni\year3\SemesterA\ImageProcessing\Exercises\Ex2\external\monkey.jpg', 1)


b = [[-1, 1]]

print(np.median(convolve2d(im, b)))

output: 0.0

The read_image function is my own and this is the implementation:

from imageio import imread
from skimage.color import rgb2gray
import numpy as np
def read_image(filename, representation):
    """
    Receives an image file and converts it into one of two given representations.
    :param filename: The file name of an image on disk (could be grayscale or RGB).
    :param representation: representation code, either 1 or 2 defining wether the output
    should be a grayscale image (1) or an RGB image (2). If the input image is grayscale,
    we won't call it with representation = 2.
    :return: An image, represented by a matrix of type (np.float64) with intensities
    normalized to the range [0,1].
    """
    assert representation in [1, 2]

    # reads the image
    im = imread(filename)
    if representation == 1:  # If the user specified they need grayscale image,
        if len(im.shape) == 3:  # AND the image is not grayscale yet
            im = rgb2gray(im)  # convert to grayscale (**Assuming its RGB and not a different format**)

    im_float = im.astype(np.float64)  # Convert the image type to one we can work with.

    if im_float.max() > 1:  # If image values are out of bound, normalize them.
        im_float = im_float / 255

    return im_float

Edit 2: I tried it on several different images, and got 0.0 at all of them. The image I'm using in the example is:

Answer 1

I computed derivatives using different methods such as :

1. convolution with an array [[-1, 1]].
2. Using the fourier theorem by computing DFT of the image and the array mentioned above, multiplying them and performing IDFT.
3. Directly through the derivative formula (Computing Fourier, multiplying by index and a constant and computing the inverse).

These derivative methods are all approximate and make different assumptions:

1. Convolution by [[-1, 1]] computes differences between adjacent elements,

   derivative ~= data[n+1] − data[n]
   

   You can interpret this like interpolating the data with a line segment, then taking the derivative of that interpolant:

   I(x) = data[n] + (data[n+1] − data[n]) * (x − n)
   

   So the approximation assumes the underlying function is locally linear. You can analyze the error by Taylor expansion to find that the error comes from the ignored higher-order terms. In other words, the approximation is accurate provided the function doesn't have strong nonlinear terms. This is a simple case of finite differences.

2. This is the same as 1, except with different boundary handling to handle convolution of samples near the edges of the image. By default, scipy.signal.convolve2d does zero padding (though you can use the boundary option to choose some other methods). However when computing the convolution through the DFT, then implicitly the boundary handling is periodic, wrapping around at the image edges. So the results of 1 and 2 differ for a margin of pixels near the edge because of the different boundary handling.

3. Computing the derivative through multiplying iω under the DFT representation can be interpreted like evaluating the derivative of the sinc interpolation the data. Sinc interpolation assumes the data is band limited. The error comes from spectra beyond the Nyquist frequency. Particularly, if there is a hard jump discontinuity from an object boundary, then the image is not bandlimited and the DFT-based derivative will have substantial error in the vicinity of the jump, appearing as ringing artifacts.

The main thing that baffles me is that when I try computing the median of this derivative, its ALWAYS 0.0.

I don't know why this happened here, but it shouldn't always be the case. For instance if each image row is the unit ramp data[n] = n, then the convolution by [[-1, 1]] is equal to 1 everywhere, except depending on boundary handling possibly not at the edges, so the median is 1.

Answer 2

Pascal already gave a wonderful explanation of the differences between the various approximations to the derivative. So I'll focus here on the "why always 0.0?" question.

The median of the derivative is 0.0 only by approximation. When I compute it, based on the finite difference approximation (method #1), I get -5.15e-5 as the median. Close to zero, but not exactly zero.

The derivative is 0 in uniform (flat) regions of the image such as the out-of-focus background. Other features in the image tend to have both a positive and a negative edge, making the histogram of the derivative image very symmetric:

This symmetry causes the median (as well as the mean) to be close to zero for such an image. However, this is not always the case. For example, if the image is brighter on the left edge than the right edge (or the other way around), then there must be a net gradient across the image, causing the mean or median to be different from zero.