It is mentioned in some literature that we should always use central difference when computing the derivatives of an image instead of forward or backward difference. Does anyone knows why is that?

Central difference = $\frac{df(x)}{dx} = \frac{f(x+h) - f(x-h)}{2h}$

Forward difference = $\frac{df(x)}{dx} = \frac{f(x+h) - f(x)}{h}$

Backward difference = $\frac{df(x)}{dx} = \frac{f(x) - f(x-h)}{h}$

  • $\begingroup$ There are much more ways to approximate the derivative : search for 'derivative filter design'. And you can use the one you'd like, if you understand what are the difference (phase, frequency response, effect on the probabilistic model for the signal). So the advantage of the central difference is that it has 0 phase : its Fourier transform is real. The 2 others (being the central difference shifted and strechted !) have a linear phase, it is not bad too. In general, we tend to avoid filters with chaotic phase (varying too much with the frequency) $\endgroup$ – reuns Jan 30 '17 at 15:39

A few things come to mind:

  • For smooth $f$, the central difference scheme is second order in $h$, whereas the other two you mentioned are first order in $h$. In other words, if $f$ is smooth, the (real space) error for the centered difference scheme is $O(h^2)$ whereas for the forward/backward schemes it is $O(h)$. Thus at fixed sufficiently small $h$, the central difference will have a better (real space) error than the other two. (Specifically, if $f \leq c_1 h$ and $g \leq c_2 h^2$ then $g<f$ once $h<\frac{c_2}{c_1}$. In this context, $c_1$ is proportional to a bound on $f''$ while $c_2$ is proportional to a bound on $f'''$.)
  • Central differencing uses the same number of points as the other two you mentioned, so there is no loss in efficiency compared to those.
  • There are higher order methods even than central differencing, and actually some of these may be even better. But very high order methods are not practical in floating point arithmetic, because they get their high order from subtraction of large nearly equal numbers, which causes loss of precision. "Middle order" methods (of order 2 through 6) are generally preferred when using double precision, because they do a good job of balancing accuracy in exact arithmetic with roundoff error when the calculations are done in double precision.

There are some other aspects of this that are more specific to signal processing. If this doesn't answer your question by itself then I can mention a little bit about that.

Upon request, here is a figure comparing the (signed) error in the central, forward, and backward derivative estimates for $e^x$ at $x=1$. The centered difference is in black, the forward difference is in blue, and the backward difference is in red.

enter image description here

  • 1
    $\begingroup$ @user1952009 No, this is just blatantly wrong. $f(x+h)=f(x)+hf'(x)+1/2 h^2 f''(x)+O(h^3)$. $f(x-h)=f(x)-hf'(x)+1/2 h^2 f''(x)+O(h^3)$. Therefore $f(x+h)-f(x-h)=2hf'(x)+O(h^3)$, so the centered difference has an error of $O(h^2)$ for three times continuously differentiable functions. Doing essentially the same calculation reveals that the forward difference has an error of $O(h)$ for twice continuously differentiable functions. So for reasonably small $h$, the centered difference scheme will perform better on the real space side. $\endgroup$ – Ian Jan 30 '17 at 16:53
  • 1
    $\begingroup$ @user1952009 No, $f(x+2h)=f(x)+f'(x)2h+O(h^2)$. Thus taking the difference and dividing by $2h$ gives you an error of $O(h)$. With the centered difference (again for three times continuously differentiable functions) the error is merely $O(h^3)$, because the second derivative terms cancel one another. $\endgroup$ – Ian Jan 30 '17 at 17:02
  • 1
    $\begingroup$ Again the second derivative terms cancel one another in the centered case: $(f(x)+hf'(x)+1/2 h^2 f''(x)+O(h^3))-(f(x)-hf'(x)+1/2 h^2 f''(x)+O(h^3))=2hf'(x)+O(h^3)$. This is exceptionally basic numerical analysis, I'm not sure why I have to explain it 5 times to someone who clearly knows some heavy math. $\endgroup$ – Ian Jan 30 '17 at 17:03
  • 1
    $\begingroup$ @user1952009 The fact that $h$ is fixed does not matter. Once $h$ is small enough, centered difference will perform better (in real space error metric). I see what you mean about one being a shift of the other, so that on the Fourier side one can convert between the two through a phase factor...but this can be viewed as a phase discrepancy, i.e. a flaw in the forward/backward schemes which happens to be easy to correct on the Fourier side. I also disagree that there is enough context here to really be certain of what exactly the OP wants. $\endgroup$ – Ian Jan 30 '17 at 18:10
  • 1
    $\begingroup$ @Ray.R.Chua Is this better? $\endgroup$ – Ian Jan 31 '17 at 15:34

You already have got a couple of good relevant points, so I'm just gonna add one I haven't seen so far among the answers.

  • The result of an operator with a well defined center pixel is on the same grid where you could argue that forward or backward difference are off by a fraction of 1/2 samples in either dimension (compared to the in-grid), this could be impractical for many reasons and in many circumstances.

Here we can see how the grid moves from before and to after forward x-difference and forward y-difference respectively: None of the 3 grids are the same! enter image description here

And for the central difference scheme all of them align nicely on top of each other: enter image description here

And as many things in engineering and science estimating a differential is one of many operations which need to be compatible or align in some sense for them to be useful.


As proposed by @Ian in the comment to avoid the grids not being aligned, one could use (expressed in the functions Z-transform):

$$\mathcal{Z}\{d_k(x_1,x_2)\} = 1-{z_k}^2$$ But we see that it will simply be a lazy filtering by one grid step of the central difference: $$\mathcal{Z}\{d_k(x_1,x_2)\} = \underset{\text{lazy}}{\underbrace{z_k}}\underset{\text{central diff.}}{\underbrace{({z_k}^{-1}-{z_k}^1)}}$$

So in practice we would be calculating the exact same thing, just moving the result one step to the left (or down).

  • $\begingroup$ I don't see what you mean here: why not just have $h$ be the same as the grid size? $\endgroup$ – Ian Jan 30 '17 at 18:11
  • $\begingroup$ Yes then we would at least get from one set to the same set of middle points in the convolution results. So in some sense that is better. But in practice it really just becomes a permutation of the central difference filter in the respective dimensions. $\endgroup$ – mathreadler Jan 30 '17 at 18:18
  • $\begingroup$ @mathreadler, what is $z^{-1}_{k}$ here? $\endgroup$ – Ray.R.Chua Jan 31 '17 at 14:51
  • $\begingroup$ It's coefficient is the filter tap at coordinate -1 for dimension k and 0 for all other dimensions (their exponents are all 0 so they vanish). $\endgroup$ – mathreadler Jan 31 '17 at 14:59

One of the reasons, not unique to image processing, is that the central difference produces a 2nd order error, while the others - 1st order, so usually the central difference schemes have better convergence and/or stability properties.

  • $\begingroup$ Not always, you may need more regularization for higher frequencies. $\endgroup$ – mathreadler Feb 2 '17 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.