# On Centered and Compact Signal and Image Derivatives for Feature Extraction

Konstantinos K. Delibasis, Aristides Kechriniotis, Ilias Maglogiannis

► To cite this version:

Konstantinos K. Delibasis, Aristides Kechriniotis, Ilias Maglogiannis. On Centered and Compact Signal and Image Derivatives for Feature Extraction. 9th Artificial Intelligence Applications and Innovations (AIAI), Sep 2013, Paphos, Greece. pp.318-327, 10.1007/978-3-642-41142-7_33. hal-01459628

HAL Id: hal-01459628

https://inria.hal.science/hal-01459628

Submitted on 7 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On Centered and Compact Signal and Image derivatives for feature extraction

Konstantinos K. Delibasis¹, Aristides Kechriniotis², Ilias Maglogiannis³

¹ Univ. of Central Greece, Dept. of Computer Science and Biomedical Informatics, Lamia, Greece

² TEI of Lamia, Dept. of Electronics, Lamia, Greece

³University of Piraeus, Dept. of Digital Systems, Piraeus, Greece
kdelibasis@yahoo.com; kechrin@teilam.gr; imaglo@unipi.gr

**Abstract** A great number of Artificial Intelligence applications are based on features extracted from signals or images. Feature extraction often requires differentiation of discrete signals and/or images in one or more dimensions. In this work we provide two Theorems for the construction of finite length (finite impulse response -FIR) masks for signal and image differentiation of any order, using central differences of any required length. Moreover, we present a very efficient algorithm for implementing the compact (implicit) differentiation of discrete signals and images, as infinite impulse response (IIR) filters. The differentiator operators are assessed in terms of their spectral properties, as well as in terms of the performance of corner detection in gray scale images, achieving higher sensitivity than standard operators. These features are considered very important for computer vision systems. The computational complexity for the centered and the explicit derivatives is also provided.

**Keywords:** Signal analysis, computer vision, image feature extraction, signal discrete derivatives, centered image derivatives, compact (implicit) image derivatives.

## 1 Introduction

Analysis and feature extraction of discrete signals and images is often an integral part of any signal or image based artificial intelligence system. Local features in signals and images are very frequently extracted using information about the derivatives of the single or multidimensional signal. For instance, image edge enhancement is based on derivatives and allows better segmentation and pattern recognition results for computer vision systems. The Harris corner detection [1] is a simple and reliable image feature point detector that uses the image gradient information. A number of quite popular feature extraction methods [2-4] utilize the image Hessian matrix, which is computed for each image pixel using first and second order partial image derivatives. The Scale Invariant Feature Transform (SIFT) of an image [5] approximates the required image Laplacian operator with difference of Gaussians (DoG).
The derivatives of any given sequence of numbers may be approximated by convolution of one-dimensional masks using centered finite differences [6]. Most image processing textbooks provide a set of small size two-dimensional masks, such as the Robert's or Sobel masks which are convolved with the image to provide its two partial derivatives (with respect to rows or columns). The use of the Laplacian, as well as the Laplacian of Gaussian operator (LoG) operator is very common for second order derivatives. However, as the size of the convolving mask increases, the more accurate the numerical approximation becomes.

The contribution of this work concerns the provision of an algorithm for calculating the masks to be used for approximating the derivatives of any order applied on a discrete signal/image. The size of the convolving mask may be arbitrarily chosen. We also consider the use of compact or implicit derivatives, as described in the work of Lele [7]. The advantage of implicit derivatives is that they achieve better spectral properties the centered differences, while they require less support points. Although compact derivatives have been applied to areas of numerical analysis, their application in image processing has been only recently investigated [8]. In this work we propose an implementation scheme based on an IIR prefiltering step as described in Unser et al [9, section IIA]. The proposed scheme is computationally more efficient than the standard matrix formulation, as it will be discussed in detail.

## 2 Proposed Methodology

In this section, we present the theorems for generating the finite impulse response (FIR) filters for signal derivatives of odd and even order. These filters are essentially symmetric finite differences; therefore we will use the terms, central, centered of FIR interchangeably. We will also present the Equations for generating the 1st and 2nd order compact derivatives. A very fast numerical implementation of central derivatives is also presented.

### 2.1 Central Derivatives

Consider the nxn matrices $O_n, E_n$ defined by

$$O_n = \begin{bmatrix} 1 & 2 & \cdots & n \\ 1 & 2^3 & \cdots & n^3 \\ \vdots \\ 1 & 2^{2n-1} & \cdots & n^{2n-1} \end{bmatrix}, \quad E_n = \begin{bmatrix} 1 & 2^2 & \cdots & n^2 \\ 1 & 2^4 & \cdots & n^4 \\ \vdots \\ 1 & 2^{2n} & \cdots & n^{2n} \end{bmatrix} \quad (1)$$

For $i = 1,...,n; k = 1,...,n-1$ we denote by $DO_{i,k}$ and $DE_{i,k}$ the determinant of order $n-1$, which follows by deleting the $i^{th}$ row, by replacing the $k^{th}$ column with the nth column and then deleting the nth column of $O_n$ and $E_n$ respectively. Further, we denote by $DO_{i,n}$ and $DE_{i,n}$ the negative of the determinant of order $n-1$, obtained by deleting the ith row and nth column of $O_n$ and $E_n$ respectively.
**Example.** For $n=5$, we present the 5x5 matrix $O_5$ according to (1) and the exemplar matrices $DO_{2,4}$ and $DO_{2,5}$:

$$O_5 = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 \\ 1 & 2^3 & 3^3 & 4^3 & 5^3 \\ 1 & 2^5 & 3^5 & 4^5 & 5^5 \\ 1 & 2^7 & 3^7 & 4^7 & 5^7 \\ 1 & 2^9 & 3^9 & 4^9 & 5^9 \end{bmatrix}, \quad DO_{2,4} = \begin{bmatrix} 1 & 2 & 3 & 5 \\ 1 & 2^5 & 3^5 & 5^5 \\ 1 & 2^7 & 3^7 & 5^7 \\ 1 & 2^9 & 3^9 & 5^9 \end{bmatrix}, \quad DO_{2,5} = - \begin{bmatrix} 1 & 2 & 3 & 4 \\ 1 & 2^5 & 3^5 & 4^9 \\ 1 & 2^7 & 3^7 & 4^9 \\ 1 & 2^9 & 3^9 & 4^9 \end{bmatrix}$$

We can now present the theorems for approximating the odd and even-order derivative of any given sequence.

**Theorem 1.** Given a real sequence $f(x)$, its $n$-point derivative $f^{(2i-1)}(x)$ of any odd order $(2i-1)$ with $i=1,...,n-1$, is given by

$$f^{(2i-1)}(x) = \frac{(2i-1)!}{2h^{2i-1} \sum_{k=1}^{n} k^{2i-1} DO_{i,k}} \sum_{k=1}^{n} DO_{i,k} (f(x+kh) - f(x-kh)) + O(h^{2n+1}) \quad (2)$$

**Theorem 2.** Given a real sequence $f(x)$, its $n$-point derivative $f^{(2i)}(x)$ of any even order $(2i)$ with $i=1,...,n-1$, is given by:

$$f^{(2i)}(x) = \frac{(2i)!}{2h^{2i} \sum_{k=1}^{n} k^{2i} DE_{i,k}} \left[ \frac{1}{2} \sum_{k=1}^{n} DE_{i,k} (f(x+kh) - f(x-kh)) + f(x) \sum_{k=1}^{n} DE_{i,k} \right] + O(h^{2n+1}) \quad (3)$$

Assuming that $f(x)$ is sampled with a constant sampling frequency, we may set $h=1$. Now, equations (2) and (5) can be written as linear convolution:

$$f^{(2i-1)}(x) = (f * M_{\text{odd}})(x) \quad (4)$$

$$f^{(2i)}(x) = (f * M_{\text{even}})(x) \quad (5)$$

where $M_{\text{even}}(k)$ and $M_{\text{odd}}(k)$ for $k=-n,...,n$ are defined as

$$M_{\text{even}}(k) = \frac{(2i)! DE_{i,k}}{2 \sum_{k=1}^{n} k^{2i} DE_{i,k}} \left[ -1, -1, \dots, -1, \frac{1}{2}, 1, \dots, 1, 1 \right] \quad (6)$$

$$M_{\text{odd}}(k) = \frac{(2i-1)! DO_{i,k}}{2 \sum_{k=1}^{n} k^{2i-1} DO_{i,k}} \left[ -1, -1, \dots, -1, 0, 1, \dots, 1, 1 \right] \quad (7)$$

The one-dimensional masks for the first and second order derivatives, created by Equations (6) and (7) using up to 11 points ($n=5$), are provided in Table 1, in rational form. The 3x3 Sobel operator can be generated by the following linear convolution, by exploiting its separability property, using the first mask of Table 1, as
$[1/2 \ 0 \ -1/2]*[1 \ 2 \ 1]^T$. The well-known Laplacian operator can be generated by the following linear convolution: $[1 \ -2 \ 1]*[1 \ -2 \ 1]^T$, using the 3-point 2nd order derivative mask from Table 1.

<table><caption>**Table 1.** The resulting masks $M_{even}$ and $M_{odd}$ for first and second order differentiation.</caption><thead><tr><th>order of deriv.</th><th>n</th><th>Num. of points in Mask</th><th>Mask</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>3</td><td>[1/2, 0, -1/2]</td></tr><tr><td>1</td><td>2</td><td>5</td><td>[-1, 8, 0, -8, 1]/12</td></tr><tr><td>1</td><td>3</td><td>7</td><td>[1, -9, 45, 0, -45, 9, -1]/60</td></tr><tr><td>1</td><td>4</td><td>9</td><td>[-1/280, 4/105, -1/4, 0, 1/4, -4/105, 1/280]</td></tr><tr><td>1</td><td>5</td><td>11</td><td>[4, -5, +30, -120, +430, 0, -430, +120, -30, 5, -4]/504</td></tr><tr><td>2</td><td>1</td><td>3</td><td>[1, -2, 1]</td></tr><tr><td>2</td><td>2</td><td>5</td><td>[-1/12, 4/3, -2.5, 4/3, -1/12]</td></tr><tr><td>2</td><td>3</td><td>7</td><td>[1/90, -3/20, 1.5, -49/18, 1.5, -3/20, 1/90]</td></tr><tr><td>2</td><td>4</td><td>9</td><td>[-1/560, 8/315, -1/5, 8/5, 205/144, 8/5, -1/5, 8/315, -1/560]</td></tr></tbody></table>

## 2.2 Compact (Implicit) derivatives

Let us denote by $f_n=f(n)$ the values of a discrete signal with $n$ integer and by $f'$ the values of the signal's derivative of the first order. According to Eq.(2.1) in the work of Lele [7], the following holds:

$$q_2 f'_{i-2} + q_1 f'_{i-1} + f'_{i} + q_1 f'_{i+1} + q_2 f'_{i+2} = \left( \frac{c_3}{6}, \frac{c_2}{4}, \frac{c_1}{2}, 0, -\frac{c_1}{2}, -\frac{c_2}{4}, -\frac{c_3}{6} \right) * f \quad (8)$$

where $q_1$, $q_2$ and $c_1$, $c_2$, $c_3$ are coefficients, whose value defines the order of accuracy of the 1st derivative approximation and the symbol $*$ stands for linear convolution. The order of accuracy for derivative approximation is not to be confused with the order of the derivative itself. The second order derivative $f_i''$ can also be described in a manner similar to (8), according to (2.2) in [7], for the one-dimensional case:

$$q_2 f''_{i-2} + q_1 f''_{i-1} + f''_{i} + q_1 f''_{i+1} + q_2 f''_{i+2} = \left( \frac{c_3}{9}, \frac{c_2}{4}, \frac{c_1}{2}, -2 \left( \frac{c_3}{9} + \frac{c_2}{4} + \frac{c_1}{2} \right), \frac{c_1}{2}, \frac{c_2}{4}, \frac{c_3}{9} \right) * f \quad (9)$$

According to Equations (8), (9), the values of $f'$ and $f''_i$ are not obtained explicitly with a simple convolution (as in the case of centered difference (4), (5)), thus the term *implicit*. More specifically, the right part of (8) and (9) is a linear convolution with a symmetric kernel of up to 7 points, whereas the left hand side involves inverse convolution with a symmetric kernel $Mask_1$ of up to 5 points, where $Mask_1[k] = (q_2, q_1, 1, q_1, q_2), k=-2, -1, \dots, 2$.
The 6th, 8th and 10th order of accuracy of the 1st derivative can be obtained by setting the $q_1$, $q_2$ and $c_1$, $c_2$, $c_3$ parameters in (8) according to (2.1.12) and (2.1.14) in Lele's work [7]. Similarly, the 4th, 6th, 8th and 10th order of accuracy of the 2nd derivative can be obtained by setting the $q_1$, $q_2$ and $c_1$, $c_2$, $c_3$ parameters in (9) according to (2.2.7) - (2.1.11) in [7]. All the necessary parameters in (8) and (9) are provided in Table 2.

<table><caption>**Table 2.** The values of the parameters of Eq.(8) and Eq.(9) for the 1st and 2nd order derivative.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="3">1st order Derivative (8)<br>with order of accuracy</th><th colspan="3">2nd order Derivative (9)<br>with order of accuracy</th></tr><tr><th>6</th><th>8</th><th>10</th><th>4</th><th>6</th><th>8</th><th>10</th></tr></thead><tbody><tr><td>$q_1$</td><td>1/3</td><td>4/9</td><td>1/2</td><td>0</td><td>2/11</td><td>344/1179</td><td>334/899</td></tr><tr><td>$q_2$</td><td>0</td><td>1/36</td><td>1/20</td><td>0</td><td>0</td><td>23/2358</td><td>43/1798</td></tr><tr><td>$c_1$</td><td>14/9</td><td>40/27</td><td>17/12</td><td>4/3</td><td>12/11</td><td>320/393</td><td>1065/1798</td></tr><tr><td>$c_2$</td><td>1/9</td><td>25/54</td><td>101/105</td><td>-1/3</td><td>3/11</td><td>310/393</td><td>1038/899</td></tr><tr><td>$c_3$</td><td>0</td><td>0</td><td>1/100</td><td>0</td><td>0</td><td>1/213</td><td>79/1798</td></tr></tbody></table>

In matrix form, the problem obtaining the required values $f'$ and $f''$ from (8) and (9) respectively, is reformulated as solving a linear system of simultaneous equations, whose matrix is symmetric pentadiagonal. Thus, it can be solved efficiently using Gauss elimination with complexity O(10N) [10]. The contribution of this work concerning implicit derivatives, is the implementation based on the Infinite impulse response (IIR) prefiltering algorithm described in the work of Unser et al [9, section IIA]. This algorithm was proposed as a very efficient implementation of the b-spline interpolation in signals and images and as it will be shown, it achieves complexity of O(8N). According to this approach, the left hand side of (8) and (9) may be written as filtering with the following transfer function,

$$H(z) = \frac{1}{q_2 z^2 + q_1 z + 1 + q_1 z^{-1} + q_2 z^{-2}} \quad (10)$$

The roots of the denominator in (11) are real pairs of reciprocals: $\rho_1, \rho_1^{-1}, \rho_2, \rho_2^{-1}$ with $|\rho_1|, |\rho_2| < 1$. In this case, the algorithm (2.5) in [9] may be applied directly, according to which, the calculation of the unknown $f'$ and $f''$ is performed as a cascade of a first order causal and an anti-causal linear filter with infinite impulse response (IIR). Let us denote by $x(k)$ the signal/image values (which is considered as input) and by $y^+(k)$ and $y(k)$ the output of the causal and the anti-causal part of (10), respectively (see Fig. 1). The linear difference equations that implements the $H(z)$ is shown in two steps (forward and inverse):

$$\begin{aligned} y^+(k) &= x(k) + \rho_i y^+(k-1), k = 2, 3, \dots, K \\ y(k) &= \rho_i (y(k+1) - y^+(k)), k = K-1, \dots, 1 \end{aligned} \quad (11)$$

The necessary initial conditions for (11) are:
$$y^{+}(1) = \sum_{k=1}^{K_0} z_i^{|k|-1} x(k), y(K) = -\frac{z_i}{1-z_i^2} (2y^{+}(K) - x(K)) \quad (12)$$

where $K_0$ a small positive integer [9]. Equations (11) and (12) are repeated for $i=1,2$, which correspond to the two roots $|\rho_1|, |\rho_2|$ of the denominator of (10).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** The implementation of the implicit derivative as an IIR filter, decomposed into a cascade of a causal (forward) and an anticausal (inverse) filters</figcaption></figure>

The application of image partial derivatives takes place in a straightforward manner, applying the one-dimensional derivative approximation along image lines and subsequently along columns. For instance, in order to calculate the quantity $\frac{\partial^3 I}{\partial x \partial y^2}$ (image 2nd partial derivative along rows and 1st partial derivative along columns), we apply the one-dimensional calculation of first order derivative along each image line and subsequently the result is used as input column-wise into the 1-D approximation of the second order derivative.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Magnitude of Frequency response of the 1st and 2nd order central derivatives and compact derivatives (the order of accuracy for the later is also provided).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** The gradient along columns and lines of the well-known cameraman image, using the implicit (order of accuracy=8) and the centered (7-point) filters.</figcaption></figure>

## 3 Experimental Results

The frequency response of the central derivatives with 3 (Sobel), 5, 7, 9 and 11 points (equivalently *n*=2, 3, 4 and 5 in (4)) is shown in Fig. 2(a). In the same Figure, the frequency response of the compact derivative with order of accuracy equal to 6, 8 and 10 (equivalently 5, 5 and 7-point stencil) is also shown. It is clear that the central difference (FIR) derivatives follow closer the ideal differentiator as the number of points in the finite length mask increases. The compact derivatives exhibit better spectral properties than the central ones. The same results are shown for the 2nd order derivative in Fig. 2b. It can also be observed that the compact (IIR) derivative filter follows more accurately the frequency response of the ideal differentiator of the first and second order, compared to the centered (FIR) filters. The FIR convolution masks as well as the IIR derivative filters were evaluated for the tasks of edge detection and corner points extraction using the Harris detector [1]. Initial results, assessed visually, were obtained using various images from [http://www.imageprocessingplace.com/DIP-3E/dip3e_book_images_downloads.htm] and from routine clinical computer tomography (CT) studies.

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The results of corner detection in a typical gray scale photograph and in a transverse CT of the abdomen, using the implicit derivative mask with order of accuracy=10 (1st raw), the 7-point FIR derivative mask generated by the proposed algorithm (2nd raw) and the 7-point derivative proposed in [12] (3rd raw).</figcaption></figure>

In Fig. 3 the first order differentiator filters (first row: implicit with 8th order of accuracy, second row: 7-point centered) are applied for gradient calculation for the well-known cameraman image. Figure 4 shows the result of the application of the Harris corner detection [1], on two different gray scale images: a photograph of a building (routinely used for algorithm evaluation) and a typical clinical transverse slice of a CT abdomen study. The Harris corner detection is a simple and popular method for discovering salient points in images. Salient points are very useful for performing image-based robot localization, image spatial registration, scene recognition etc. The Harris corner detection is based on the image Hessian matrix, whose 1st and 2nd order derivatives are calculated using the implicit derivative mask with order of accuracy=8 (1st raw), the 7-point FIR derivative mask generated by the proposed algorithm (2nd raw) and the 7-point derivative proposed in [12] (3rd raw). The source code provided in [11] is used for the Harris detector. As it can be observed, the detector using the IIR-based implicit derivative is significantly more sensitive to detecting corners in images than the 7-point FIR differentiator generated by (7), whereas the detector based on the 7-point differentiator [12] is the least sensitive.
# 4 Discussion and Conclusions

As shown in the previous section (Fig. 4), the sensitivity of the corner detector is increased when using the IIR differentiation filter, compared with the centered differences (FIR) filter. The number of arithmetic operations per sample for the centered derivative filters of even and odd order is $2(n+1)$ and $2n+1$ respectively, where $n$ is defined as in Theorem 1 and 2 (the symmetry of the masks has been taken into consideration). The computational complexity of the implicit derivatives is decomposed into the cost of the left hand side (implicit part) and the right hand side (explicit linear convolution) of (8) and (9). The former is 2 arithmetic operations for the forward and 2 for the inverse filtering for each root of the denominator (pole) in (10) with magnitude less than 1, for each image pixel (for each image dimension). By observing the number of poles of (10) as a function of the accuracy order, Table 3 can be constructed. Thus, the arithmetic complexity of $O(8N)$ achieved by the proposed IIR algorithm is superior to the $O(10N)$ achieved by the Gaussian elimination of the penta-diagonal matrix, with $N$ denoting the number of pixels.

Furthermore, the 1st order compact (implicit) derivative filter with 6th order of accuracy with computational complexity of $O(9N)$, exhibits significantly better spectral properties than the 11-point FIR differentiator with $O(11N)$, as depicted in Fig. 2(a). The 2nd order implicit differentiator with order of accuracy equal to 6 with $O(10N)$ exhibits marginally better spectral properties, compared to 11-point centered filter with $O(12N)$. Finally, it has to be reminded that the IIR differentiation filters cannot be used in real-time signal processing, due to their non-causal component.

<table><caption>**Table 3.** Number of arithmetic operations for implicit derivatives, as a function of $N$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">1st order Derivative<br>with order of accuracy</th><th colspan="4">2nd order Derivative with<br>order of accuracy</th></tr><tr><th>6</th><th>8</th><th>10</th><th>4</th><th>6</th><th>8</th><th>10</th></tr></thead><tbody><tr><th>Left side of (8) and (9)</th><td>4N</td><td>8N</td><td>8N</td><td>0</td><td>4N</td><td>8N</td><td>8N</td></tr><tr><th>Right side of (8)</th><td>5N</td><td>5N</td><td>7N</td><td colspan="4">N/A</td></tr><tr><th>Right side of (9)</th><td colspan="3">N/A</td><td>6N</td><td>6N</td><td>8N</td><td>8N</td></tr><tr><th>Total</th><td>9N</td><td>13N</td><td>15N</td><td>6N</td><td>10N</td><td>16N</td><td>16N</td></tr></tbody></table>

# Appendix

We will present the proof for the odd-order derivatives. The proof for the even-order is in total analogy.

## Proof of Theorem 1

Using the Taylor's Theorem with an integral form of remainder we get

$$f(x+t) = \sum_{j=1}^{2n} \frac{t^j f^{(j)}(x)}{j!} + R_{2n}(t)$$
$$ \text{where } R_{2n}(t) = \frac{1}{(2n+1)!} \int_{x}^{x+t} f^{(2n+1)}(u)(x+t-u)^{2n} du $$

We apply the previous relation 2n times for $t=kh$, $k=\pm 1, \pm 2, \dots, \pm n$. Then for any real constants $c_1, \dots, c_n$ we get

$$ \frac{1}{2} \sum_{k=1}^{n} c_k (f(x+kh) - f(x-kh)) = \sum_{j=1}^{n} \frac{h^{2j-1} f^{(2j-1)}(x)}{(2j-1)!} \sum_{k=1}^{n} c_k k^{2j-1} + \frac{1}{2} \sum_{k=1}^{n} c_k (R_{2n}(kh) - R_{2n}(-kh)) $$

If we choose $c_n=-1$ and find $c_k$ for $k=1,2,\dots,n-1$ by solving the following system of linear equations, $\sum_{k=1}^{n} c_k k^{2i-1} = n^{2j-1}$, $j=1,2,\dots,n$, $i \neq j$, we obtain the required result.

## References

1. Harris C., Stephens M.: A combined corner and edge detector, in ALVEV Vision Conference, pp. 147–151 (1988)

2. Lindeberg T., Garding J.: Shape-adapted smoothing in estimation of 3-D shape cues from affine deformations of local 2-D brightness structure, Image and Vision Computing, 15(6), 415–434 (1997)

3. Mikolajczyk K., Schmid C.: Scale and affine invariant interest point detectors, International Journal of Computer Vision, 1(60), 63–86 (2004).

4. Bay H., Ess A., Tuytelaars T., Gool van L.: Speeded-up robust features (SURF), International Journal on Computer Vision and Image Understanding, 110(3), 346–359 (2008)

5. Lowe D., Distinctive Image Features from Scale-Invariant Keypoints, International Journal of Computer Vision, 60(2), 91–110 (2004)

6. Keller H. B., Pereyra V.: Symbolic Generation of Finite Difference Formulas, Mathematics of Computation, 32(144), 955–971 (1978)

7. Lele S.K.: Compact difference Schemes with Spectral-like Resolution, Journal of Computational Physics, 103, 16–42 (1992)

8. Belyaev A.: On implicit image derivatives and their applications, in Hoey J., McKenna S., Trucco E. (eds.) BMVC 2011, Dundee, Scotland, UK (2011).

9. Unser M., Aldroubi A., Eden M.: B-spline signal processing: Part II - Efficient Design and Applications, IEEE Trans. Signal Process., 41(2), 834–848 (1993).

10. Benkert K., Fischer R.: An Efficient Implementation of the Thomas-Algorithm for Block Penta-diagonal in vector Computers, in 7th International Conference on Computer Science, vol. 4487, pp. 144–151, Beijing (2007).

11. MATLAB and Octave Functions for Computer Vision and Image Processing http://www.csse.uwa.edu.au/~pk/research/matlabfnss/

12. Farid H., Simoncelli E.: Differentiation of Discrete Multi-Dimensional Signals, IEEE Trans. Image Processing 13(4): 496–508 (2004).