1 Introduction
Sparse representation (SR), the linear decomposition of a signal using a few atoms of a dictionary, has led to notable results for various image processing and computer vision tasks [1, 2]. While significant steps have been taken towards expanding the theory of SR, such representations in nonEuclidean spaces have received comparatively little attention. This paper tackles the problem of sparse coding within the space of symmetric positive definite (SPD) matrices.
SPD matrices are fundamental building blocks in computer vision and machine learning. A notable example is the covariance descriptor
[3], which offer a compact way of describing regions/cuboids in images/videos and fusion of multiple features. Covariance descriptors have been exploited in several applications, such as diffusion tensor imaging [4], action recognition [5, 6, 7], pedestrian detection [3], face recognition [8, 9], texture classification [9, 10], and tracking [11].SPD matrices form a cone of zero curvature and can be analysed using the geometry of Euclidean space. However, several studies have shown that a Riemannian structure of negative curvature is more suitable for analysing SPD matrices [4, 12]. More specifically, Pennec et al. [4] introduced the Affine Invariant Riemannian Metric (AIRM) and showed that the induced Riemannian structure is invariant to inversion and similarity transforms. The AIRM is perhaps the most widely used similarity measure for SPD matrices. Nevertheless, efficiently and accurately handling the Riemannian structure is nontrivial as basic computations on Riemannian manifolds (such as similarities and distances) involve nonlinear operators. This not only hinders the development of optimisation algorithms but also incurs a significant numerical burden.
To address the above drawbacks, in this paper we propose to perform the sparse coding of SPD matrices by embedding Riemannian manifolds into reproducing kernel Hilbert spaces (RKHS) [13]. This is in contrast to directly embedding into Euclidean spaces [7, 6, 14].
Related Work. Sra et al. [14] used the cone of SPD matrices and the Frobenius norm as a measure of similarity between SPD matrices. While this results in a regularised nonnegative leastsquares approach, it does not consider the Riemannian geometry induced by AIRM.
Guo et al. [6] and Yuan et al. [7] separately proposed to solve sparse representation by a logEuclidean approach, where a Riemannian problem is converted to an Euclidean one by embedding manifolds into tangent spaces. While logEuclidean approaches benefit from simplicity, the true geometry of the manifold is not taken into account. More specifically, on a tangent space only distances to the pole of space are true geodesic distances. As such, the pairwise distances between arbitrary points on the tangent space do not represent the structure of the manifold.
Sivalingam et al. [9] used Burg divergence [15] as a metric and reformulated the Riemannian^{1}^{1}1We loosely use ‘Riemannian’ to refer to the Riemannian manifold formed by SPD matrices. SR problem as a determinant maximisation problem. This has the advantage of avoiding the explicit manifold embedding, as well as resulting in a convex MAXDET problem [16] that can be solved by interior point methods. However, there are two downsides: the solution is computationally very expensive, and the relations between Burg divergence and the geometry of Riemannian manifolds were not well established.
Contributions. With the aid of the recently introduced Stein kernel [17], which is related to AIRM via a tight bound, we propose a Riemannian sparse solver by embedding Riemannian manifolds into RKHS. We show that the embedding leads to a convex and kernelised version of the Lasso problem [1], which can be solved efficiently. We furthermore propose a sparsitymaximising algorithm for dictionary learning within the space of SPD matrices, closely tied to the Stein kernel. Lastly, we show that the proposed sparse coding approach obtains superior performance on several visual classification tasks (face recognition, texture classification, person reidentification), in comparison to several stateoftheart methods: tensor sparse coding [9], logEuclidean sparse representation [6, 7], Gabor feature based sparse representation [18], and Riemannian locality preserving projection [10].
We continue the paper as follows. Section 2 begins with an overview of Bregman divergence and the Stein kernel. Section 3 describes the proposed kernel solution of Riemannian sparse coding, followed by Section 4, which covers the problem of dictionary learning on Riemannian manifolds. In Section 5 we compare the performance of the proposed method with previous approaches on several visual classification tasks. The main findings and possible future directions are summarised in Section 6.
2 Background
In this section we first overview the properties of Bregman matrix divergences, including a special case known as the symmetric Stein divergence. This leads to the Stein kernel, which can be used to embed Riemannian manifolds into RKHS.
2.1 Bregman Matrix Divergences
The Bregman matrix divergence for two symmetric matrices and is defined as [15]:
(1) 
where and is a real valued, strictly convex function on symmetric matrices. Bregman divergences are nonnegative, definite, and in general asymmetric. Among the several ways to symmetrise them, the JensenShannon symmetrisation is often used [15]:
(2) 
If , then the symmetric Stein divergence is obtained from (2):
(3) 
The space induced by AIRM on symmetric positive definite matrices of dimension is a Riemannian manifold of negative curvature. For two points , the AIRM is defined as , where . The symmetric Stein divergence and Riemannian metric over manifolds are related in several aspects. Two important properties are summarised below.
Property 1
Let , and be the Thompson metric [19] with
representing the vector of eigenvalues of
. The following sandwiching inequality between the symmetric Stein divergence and Riemannian metric exists [17]:(4) 
Property 2
The curve parameterises the unique geodesic between the SPD matrices and . On this curve the Riemannian geodesic distance satisfies [12]. The symmetric Stein divergence satisfies a similar but slightly weaker result, .
The first property establishes a bound between the geodesic distance and Stein divergence, providing motivation for addressing Riemannian problems via the divergence. The second property reinforces the motivation by explaining that the behaviour of Stein divergences along geodesic curves is similar to true Riemannian geometry.
2.2 Stein Kernel
Definition 1
Let be a nonempty set on Riemannian manifold . A function is a Riemannian kernel if is symmetric for all , ie., , and the following inequality is satisfied for all :
Under a mild condition (explained afterwards), the following function forms a Riemannian kernel [17]:
(5) 
We shall refer to this kernel as the Stein kernel from here on. The following theorem states the condition under which Stein kernel is positive definite.
Theorem 2.1
Let be a set of Riemannian points. The matrix , with , defined in (5), is positive definite iff:
(6) 
Interested readers can follow the proof in [17]. For values of outside of the above set, it is possible to convert a pseudo kernel into a true kernel, as discussed for example in [20]. The determinant of an SPD matrix can be efficiently computed by Cholesky decomposition in . As such, the complexity of computing Stein kernel is .
3 Kernel Sparse Coding
Sparse coding on Riemannian manifolds in general means that a given query point on a manifold can be expressed as a sparse “combination” of dictionary elements. Our idea here is to embed the manifold into RKHS and replace the idea of “combination” on manifolds with the general concept of linear combination in Hilbert spaces. More specifically, given a Riemannian dictionary , and an embedding function , for a Riemannian point we seek for a sparse vector such that admits the sparse representation over . In other words, we are interested in solving the following kernelised version of the Lasso problem [1]:
(7) 
The first term in (7) can be expanded as:
(8) 
where and . This reveals that the optimisation problem in (7) is convex and similar to its counterpart in Euclidean space, except for the definition of and . Consequently, greedy or relaxation solutions can be adapted to obtain the sparse codes [1]. To solve problem (7) we used CVX [21], a package for specifying and solving convex programs^{2}^{2}2The SPAMS package can also be used: http://spamsdevel.gforge.inria.fr/.
3.1 Classification Using Sparse Codes
There are two main approaches for classification based on the obtained sparse codes (vectors) for a given query sample: (i) directly, and (ii)
indirectly, with the aid of an Euclideanbased classifier. We elucidate the two approaches below.
(i) If the atoms in sparse dictionary have associated class labels (ie. each atom in the dictionary is a training sample), the sparse codes can be directly used for classification. This approach is applicable only to closedset identification tasks. Let be the classspecific sparse codes, where is the class label of atom and is the discrete Dirac function [22]. An efficient way of using classspecific sparse codes is through computing residual errors [2]. In this case, the residual error of query sample for class is defined as:
(9) 
which can be computed via the use of a Riemannian kernel in a similar manner to (8). The class with the minimum residual error is deemed to represent the query. Alternatively, the similarity between query sample to class can be defined as . The function can be linear like or even nonlinear like .
(ii) If the atoms in the sparse dictionary are not labelled (eg. is a generic dictionary not tied to any particular class [23]
), the generated sparse codes (vectors) for both training and query data can be fed to Euclideanbased classifiers, such as support vector machines
[22]. The sparse code is hence interpreted as a feature vector, which in essence means that a classification problem on a Riemannian manifold is converted to an Euclidean classification problem. This approach is applicable to both closedset and openset classification tasks.4 Learning Riemannian Dictionaries
If the indirect classification of sparse codes is required (as elucidated in the preceding section) a Riemannian dictionary is first required. Given a set of Riemannian points , learning a dictionary can be formulated as jointly minimising the energy function
(10) 
over the dictionary and the sparse codes , ie., .
Among the various solutions to the problem of dictionary learning in Euclidean spaces, iterative methods like KSVD have received much attention [1]. Borrowing the idea from Euclidean spaces, we propose to minimise the energy in (10) iteratively. After initialising the dictionary , for example by Riemannian clustering using the Karcher mean [4], we iterate between a sparse coding step and a dictionary update step. In the sparse coding step, is fixed and is computed. In the dictionary update step, is fixed while is updated, with each dictionary atom updated independently.
The derivative of (10) with respect to , while and other atoms are fixed, is:
(11) 
As , (11) can be further simplified to:
(12) 
Since (12) contains linear and nonlinear terms of (eg. inverse and kernel terms), a closedform solution for computing its root, ie., , cannot be sought. As such, we propose an alternative solution by exploiting previous values of and in the updating step. More specifically, rearranging (12
) and estimating
as well as by their previous values, atom at iteration is updated using:(13) 
where
(14)  
(15) 
To avoid the degenerative case (due to numerical inconsistency), atoms are normalised by their second norms at the end of each iteration. Algorithm 1 assembles all the above details into pseudocode for dictionary learning.
5 Experiments
Two sets of experiments^{3}^{3}3Matlab/Octave source code is available at http://itee.uq.edu.au/~uqmhara1 are presented in this section. In the first set, we evaluate the performance of the proposed Riemannian SR (RSR) method (as described in Section 3) without dictionary learning. Each atom in the dictionary is a training sample. This is to contrast RSR to previous stateoftheart methods on several popular closedset classification tasks. We use the residual error approach for classification, as described in Eqn. (9).
In the second set, the performance of the RSR method is evaluated in conjunction with dictionaries learned via three methods: random, Riemannian kmeans, and the proposed dictionary learning technique (as described in Section 4).
5.1 Riemannian Sparse Representation
5.1.1 Synthetic Data.
We first consider a multiclass classification problem over using synthetic data. We compared the proposed RSR against Tensor Sparse Coding (TSC) [9] and logEuclidean Sparse Representation (logESR) [6, 7]. The data used in the experiments constitutes 512 random samples from 4 classes. Half of the samples were used for training and the rest were used as test data.
To create a Riemannian manifold, samples were generated over a particular tangent space and then mapped back to the manifold using the exponential map [12]
. The positions of tangent spaces were chosen randomly and samples in each class obeyed a normal distribution. By fixing the mean of each class and increasing the class variance we created two classification problems: ‘easy’ and ‘hard’. To draw useful statistics, the data creation process was repeated 100 times.
Table 1
shows the average recognition accuracy and the total running time (in seconds). All algorithms were implemented in Matlab and executed on a 3 GHz Intel CPU. In terms of recognition accuracy, RSR obtains superior performance when compared with previous stateoftheart approaches. We note that by increasing the class variance, samples from the four classes are intertwined, leading to a decrease in recognition accuracy. The performance of logESR is higher than TSC, which might be due the to fact that the generated data can be modelled by Gaussian distribution over tangent space, hence favouring the tangentbased solution.
Focusing on run time, Table 1 suggests that logESR has the lowest complexity while TSC has the highest. The proposed RSR method is substantially faster than TSC, while delivering the highest recognition accuracy.
logESR [6, 7]  TSC [9]  RSR (proposed)  

recognition 
accuracy 

easy  
hard  
runtime 


easy+hard  sec  sec  sec 
5.1.2 Face Recognition.
We used the ‘b’ subset of the FERET dataset [25], which includes 1400 images from 198 subjects. The images were closely cropped around the face and downsampled to . Examples are shown in Figure 1.
We performed four tests with various pose angles. Training data was composed of images marked ‘ba’, ‘bj’ and ‘bk’ (ie., frontal faces with expression and illumination variations). Images with ‘bd’, ‘be’, ‘bf’ and ‘bg’ labels (ie., nonfrontal faces) were used as test data. For Riemannianbased methods, a covariance descriptor described a face image, using the following features:
where is the intensity value at position and is the response of a 2D Gabor wavelet [26] centered at with orientation and scale :
with and .
Table 2 shows a comparison of RSR against logESR [6, 7], TSC [9], and two purely Euclidean sparse representations, PCASRC [2] and Gabor SR (GSR) [18]. In all cases the proposed RSR method obtains the highest accuracy. Furthermore, the proposed approach significantly outperforms stateoftheart Euclidean solutions, especially for test images with label ‘bg’.
PCASR [2]  GSR [18]  logESR [6, 7]  TSC [9]  RSR (proposed)  

bg  
bf  
be  
bd  
average 
5.1.3 Texture Classification.
We performed a classification task using the Brodatz texture dataset [27]. Examples are shown in Fig. 3. We followed the test protocol devised in [9] and generated nine test scenarios with various number of classes. This includes 5texture (‘5c’, ‘5m’, ‘5v’, ‘5v2’, ‘5v3’), 10texture (‘10’, ‘10v’) and 16texture (‘16c’, ‘16v’) mosaics. To create a Riemannian manifold, each image was first downsampled to and then split into 64 regions of size . The feature vector for any pixel is . Each region is described by a covariance descriptor of these features. For each test scenario, five covariance matrices per class were randomly selected as training data and the rest was used for testing. The random selection of training/testing data was repeated 20 times.
5.1.4 Person Reidentification.
We used the modified ETHZ dataset [31]. The original ETHZ dataset was captured using a moving camera [28], providing a range of variations in the appearance of people. The dataset is structured into 3 sequences. Sequence 1 contains 83 pedestrians (4,857 images), Sequence 2 contains 35 pedestrians (1,936 images), and Sequence 3 contains 28 pedestrians (1,762 images). See Fig. 4 for examples.
We downsampled all images to pixels. For each subject we randomly selected 10 images for training and used the rest for testing. Random selection of training and testing data was repeated 20 times to obtain reliable statistics. To describe each image, the covariance descriptor was computed using the following features:
where is the position of a pixel, while , and represent the corresponding colour information. The gradient and Laplacian for colour are represented by and , respectively.
We compared the proposed RSR method with several techniques previously used for pedestrian detection: Histogram Plus Epitome (HPE) [29], SymmetryDriven Accumulation of Local Features (SDALF) [30], and Riemannian Locality Preserving Projection (RLPP) [10]. The performance of logESR was below HPE method and is not shown. The results for TSC could not be generated in a timely manner, due to the heavy computational load of the algorithm.
Results for Sequence 1 and 2 are shown in Fig. 5, in terms of cumulative matching characteristic (CMC) curves. The CMC curve represents the expectation of finding the correct match in the top matches. The proposed method obtains the highest accuracy. For Sequence 3 (not shown), very similar performance is obtained by SDALF, RLPP and the proposed RSR, with HPE having the lowest performance.
5.2 Dictionary Learning
Here we compare the performance of the proposed Riemannian dictionary learning technique (as described in Section 4), with the performances of dictionaries obtained by random sampling and Riemannian kmeans. We first use synthetic data to show that the proposed method obtains a lower representation error in RKHS, followed by classification experiments on texture data.
5.2.1 Synthetic Data.
We synthesised 512 Riemannian samples from a set of 32 source points in . The source points can be considered as a form of groundtruth. The synthesised samples were then used for dictionary creation by Riemannian kmeans [24] and the proposed algorithm.
To generate each source point, an SPD matrix was created by computing the covariance of 100 random samples of a 5 dimensional normal distribution. The mean and variance of the distribution are different for each source point. To synthesise each of the 512 Riemannian samples, we uniformly selected source points and combined them with random positive weights, where the weights obeyed a normal distribution with zero mean and unit variance.
5.2.2 Texture Classification.
Here we consider a multiclass classification problem, using 111 texture images of the Brodatz texture dataset [27]. From each image we randomly extracted 50 blocks of size . To train the dictionary, 20 blocks from each image were randomly selected, resulting in a dictionary learning problem with 2200 samples. From the remaining blocks, 20 per image were used as probe data and 10 as gallery samples. The process of random block creation and dictionary generation was repeated twenty times. The average recognition accuracies over probe data are reported here. In the same manner as in Section 5.1, we used the feature vector to create the covariance, where the first dimension is the grayscale intensity, and the remaining dimensions capture first and second order gradients.
We used the proposed RSR approach to obtain the sparse codes, coupled with a dictionary generated via three separate methods: random dictionary generation, Riemannian kmeans algorithm [24], and the proposed learning algorithm (Section 4). The sparse codes were then classified using a nearestneighbour classifier.
For the randomly generated dictionary case, the classification rates are averaged over 10 runs, with each run using a different random dictionary. For all methods, dictionaries of size were trained. The best results for each approach (ie., the results for the dictionary size that obtained the highest recognition accuracy) are reported in Table 3. For the random dictionary, ; for the kmeans dictionary, ; for the proposed dictionary learning algorithm, . The results show that the proposed algorithm leads to a considerable gain in accuracy.
random  kmeans  learning 

6 Main Findings and Future Directions
With the aim of addressing sparse representation on Riemannian manifolds, proposed to seek the solution through embedding the manifolds into RKHS, with the aid of the recently introduced Stein kernel. This led to a relaxed and extended version of the Lasso problem [1] on Riemannian manifolds.
Experiments on several classification tasks (face recognition, texture classification, person reidentification) show that the proposed approach achieves notable improvements in discrimination accuracy, in comparison to stateoftheart methods such as tensor sparse coding, Riemannian locality preserving projection, and symmetrydriven accumulation of local features. We conjuncture that this stems from better exploitation of Riemannian geometry, as the Stein kernel is related to geodesic distances via a tight bound. The proposed sparse coding method is also considerably faster than the stateoftheart MAXDET reformulation used by Tensor Sparse Coding [9].
We have furthermore proposed an algorithm for learning a Riemannian dictionary, closely tied to the Stein kernel. In comparison to Riemannian kmeans [24], the proposed algorithm obtains a lower representation error in RKHS and leads to improved classification accuracies.
Future directions include using the Stein kernel for solving large margin classification problems on Riemannian manifolds. This translates to designing a machinery that maximises a margin on SPD matrices based on Stein divergence, which can be considered as an extension of support vector machines [13] to tensor spaces.
Acknowledgements. NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications and the Digital Economy, as well as by the Australian Research Council through the ICT Centre of Excellence program.
References
 [1] Elad, M.: Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing. Springer (2010)
 [2] Wright, J., Yang, A., Ganesh, A., Sastry, S., Ma, Y.: Robust face recognition via sparse representation. IEEE Trans. Pattern Analysis and Machine Intelligence (PAMI) 31(2) (2009) 210–227
 [3] Tuzel, O., Porikli, F., Meer, P.: Pedestrian detection via classification on Riemannian manifolds. IEEE Trans. Pattern Analysis and Machine Intelligence 30(10) (2008) 1713–1727
 [4] Pennec, X.: Intrinsic statistics on Riemannian manifolds: Basic tools for geometric measurements. Journal of Mathematical Imaging and Vision 25(1) (2006) 127–154
 [5] Lui, Y.M.: Advances in matrix manifolds for computer vision. Image and Vision Computing 30(6–7) (2012) 380–388
 [6] Guo, K., Ishwar, P., Konrad, J.: Action recognition using sparse representation on covariance manifolds of optical flow. In: IEEE Conf. Advanced Video and Signal Based Surveillance. (2010) 188–195
 [7] Yuan, C., Hu, W., Li, X., Maybank, S., Luo, G.: Human action recognition under logEuclidean Riemannian metric. In: Asian Conference on Computer Vision. Volume 5994 of Lecture Notes in Computer Science. Springer Berlin / Heidelberg (2010) 343–353
 [8] Pang, Y., Yuan, Y., Li, X.: Gaborbased region covariance matrices for face recognition. IEEE Transactions on Circuits and Systems for Video Technology 18(7) (2008) 989 –993
 [9] Sivalingam, R., Boley, D., Morellas, V., Papanikolopoulos, N.: Tensor sparse coding for region covariances. In: European Conference on Computer Vision (ECCV). Volume 6314 of Lecture Notes in Computer Science. (2010) 722–735
 [10] Harandi, M.T., Sanderson, C., Wiliem, A., Lovell, B.C.: Kernel analysis over Riemannian manifolds for visual recognition of actions, pedestrians and textures. In: IEEE Workshop on the Applications of Computer Vision (WACV). (2012) 433–439
 [11] Hu, W., Li, X., Luo, W., Zhang, X., Maybank, S., Zhang, Z.: Single and multiple object tracking using logEuclidean Riemannian subspace and blockdivision appearance model. IEEE Trans. Pattern Analysis and Machine Intelligence (in press) DOI: 10.1109/TPAMI.2012.42.
 [12] Bhatia, R.: Positive Definite Matrices. Princeton University Press (2007)
 [13] ShaweTaylor, J., Cristianini, N.: Kernel Methods for Pattern Analysis. Cambridge University Press (2004)
 [14] Sra, S., Cherian, A.: Generalized dictionary learning for symmetric positive definite matrices with application to nearest neighbor retrieval. In: Machine Learning and Knowledge Discovery in Databases. Volume 6913 of Lecture Notes in Computer Science. (2011) 318–332
 [15] Kulis, B., Sustik, M.A., Dhillon, I.S.: Lowrank kernel learning with Bregman matrix divergences. Journal of Machine Learning Reseach 10 (2009) 341–376
 [16] Boyd, S., Vandenberghe, L.: Convex Optimization. Cambridge University Press, New York, NY, USA (2004)
 [17] Sra, S.: Positive definite matrices and the symmetric Stein divergence. Preprint: [arXiv:1110.1773] (2012)
 [18] Yang, M., Zhang, L.: Gabor feature based sparse representation for face recognition with gabor occlusion dictionary. In: European Conference on Computer Vision (ECCV). Volume 6316 of Lecture Notes in Computer Science. Springer Berlin / Heidelberg (2010) 448–461
 [19] Thompson, A.C.: On certain contraction mappings in a partially ordered vector space. Proceedings of the American Mathematical Society 14 (1963) 438–443
 [20] Chen, Y., Garcia, E.K., Gupta, M.R., Rahimi, A., Cazzanti, L.: Similaritybased classification: Concepts and algorithms. Journal of Machine Learning Research 10 (2009) 747–776
 [21] Grant, M., Boyd, S.: CVX: Matlab software for disciplined convex programming, version 1.21. http://cvxr.com/cvx/ (April 2011)
 [22] Bishop, C.M.: Pattern Recognition and Machine Learning. Springer (2006)

[23]
Wong, Y., Harandi, M.T., Sanderson, C., Lovell, B.C.:
On robust biometric identity verification via sparse encoding
of faces: holistic vs local approaches.
In: IEEE International Joint Conference on Neural Networks. (2012) 1762–1769
 [24] Goh, A., Vidal, R.: Clustering and dimensionality reduction on Riemannian manifolds. In: IEEE Conf. Computer Vision and Pattern Recognition. (2008) 1–7
 [25] Phillips, P., Moon, H., Rizvi, S., Rauss, P.: The FERET evaluation methodology for facerecognition algorithms. IEEE Trans. Pattern Analysis and Machine Intelligence 22(10) (2000) 1090–1104
 [26] Lee, T.S.: Image representation using 2d Gabor wavelets. IEEE Trans. Pattern Analysis and Machine Intelligence 18 (1996) 959–971
 [27] Randen, T., Husøy, J.H.: Filtering for texture classification: A comparative study. IEEE Trans. Pattern Analysis and Machine Intelligence 21(4) (1999) 291–310
 [28] Ess, A., Leibe, B., Gool, L.V.: Depth and appearance for mobile scene analysis. Int. Conf. Computer Vision (ICCV) (2007) 1–8
 [29] Bazzani, L., Cristani, M., Perina, A., Farenzena, M., Murino, V.: Multipleshot person reidentification by HPE signature. In: Int. Conf. Pattern Recognition. (2010) 1413–1416
 [30] Farenzena, M., Bazzani, L., Perina, A., Murino, V., Cristani, M.: Person reidentification by symmetrydriven accumulation of local features. IEEE Conf. Computer Vision and Pattern Recognition (2010) 2360–2367
 [31] Schwartz, W.R., Davis, L.S.: Learning discriminative appearancebased models using partial least squares. In: Brazilian Symposium on Computer Graphics and Image Processing. (2009) 322–329