SHAPE-BASED IMAGE MATCHING USING HEAT KERNELS AND DIFFUSION MAPS

: 2D image matching problem is often stated as an image-to-shape or shape-to-shape matching problem. Such shape-based matching techniques should provide the matching of scene image fragments registered in various lighting, weather and season conditions or in different spectral bands. Most popular shape-to-shape matching technique is based on mutual information approach. Another well-known approach is a morphological image-to-shape matching proposed by Pytiev. In this paper we propose the new image-to-shape matching technique based on heat kernels and diffusion maps. The corresponding Diffusion Morphology is proposed as a new generalization of Pytiev morphological scheme. The fast implementation of morphological diffusion filtering is described. Experimental comparison of new and aforementioned shape-based matching techniques is reported applying to the TV and IR image matching problem.


INTRODUCTION
This paper addresses the problem of shape-based image matching with no dependence on the concrete intensity or radiometric pixel values.For example, one can compare images of one scene captured in different seasons, times of day, weather and lighting conditions, spectral ranges and so on.In such cases the traditional image matching techniques like least square matching (LSM) or zero mean normalized cross correlation (ZNCC) are usually unstable due to the variability of intensity.The most popular technique for image shape comparison utilizes the mutual information (MI) measure based on probabilistic reasoning and information theory (Maes et al, 1997).The other approach was proposed by Pytiev (so called "Pytiev morphology") based on geometrical and algebraic reasoning (Pyt'ev, 1993).This approach expresses the geometrical idea of image shape in the evident form and provides the image-to-shape matching technique based on morphological correlation coefficient (MCC).In our paper (Vizilter, Zheltov, 2012) this Pytiev morphological approach was applied for obtaining the shape-to-shape matching technique based on mean square effective morphological correlation coefficient (MSEMCC).Unfortunately, the experiments with TV and IR image matching demonstrate that all these shape-based matching techniques are not robust enough relative to noise and high-frequency damage of images.We suppose that such insufficient robustness is a result of poor quality of image segmentation, because the traditional morphological 2D shape description as a set of flat zones of frame tessellation is too sensitive relative to noise, segmentation scheme and its parameters.From the other hand, in the area of manifold learning and data shape matching some stable shape description and matching techniques are known, based on heat kernels and diffusion maps (Belkin, Niyogi, 2001), (Lafon, 2004), (Coifman et al, 2007), (Memoli, 2011).These techniques do not require preliminary data segmentation.So, there is a chance that modified morphological tools based on such "diffusion" ideas will be more robust.
In this paper we propose a new generalization of Pytiev morphological scheme based on heat kernels and diffusion maps.In the framework of this Diffusion Morphology the new image-to-shape matching technique is implemented and tested relative to the aforementioned matching techniques applying to the TV and IR images of the same scene.

RELATED WORKS
This section contains the brief overview of both shape-based 2D image comparison techniques and data description and matching techniques based on heat kernels and diffusion maps.

Mutual information
The approach for 2D shape matching based on mutual information (MI) notion from information theory was proposed in (Maes et al, 1997).Mutual information I(A,B) estimates the dependence of two random variables A and B by measuring the distance between the joint distribution p AB (a,b) and the distribution of complete independence p A (a)p B (b): where H(A) is an entropy of A, H(B) is an entropy of B, and H(A,B) is their joint entropy.For two intensity values a and b of corresponding pixels in the two images, the required empirical estimations of joint and marginal distributions can be obtained by normalization of the joint (2D) and marginal (1D) histograms of compared image fragments.Maximal I(A,B) value over a set of relative fragment positions corresponds to the best geometrical matching of image fragments.This approach does not describe the image shape in the evident form, but different successful shape matching applications were created based on this MI approach in recent years, for example (Goebel, 2005), (Ji, 2005).

Morphological image analysis and geometrical correlation
In the framework of Morphological Image Analysis (MIA) proposed by Pytiev (Pyt'ev, 1993) images are considered as piecewise-constant 2D functions where nnumber of non-intersecting connected regions of tessellation F of the frame , F={F 1 ,…,F n }; f=(f 1 ,…,f n )corresponding vector of real-valued region intensities; This tessellation of image is supposed to be obtained by some image segmentation procedure.Set of images with the same tessellation F is a convex and close subspace FL 2 () called shape-tessellation, mosaic shape or simply shape: For any image g(x,y)L 2 () the projection onto the shape F is determined as Pytiev morphological comparison of images f(x,y) and g(x,y) is performed using the normalized morphological correlation coefficients (MCC) of the following form The first formula estimates the closeness of image g to the "shape" of image f.Second formula measures the closeness of image f to the "shape" of image g.For elimination of constant non-informative part of image brightness following image normalization is usually performed: where P O fprojection of image f onto the "empty" shape O with one flat zone.This projection is a constant-valued image filled by mean value of projected image.
Later some advances and modifications of this approach were proposed in (Pyt'ev et al, 2006), (Falomkin, Pyt'ev, 2007), (Vizilter, Zheltov, 2009), but their basic ideas are still based on projection of one image to the shape of other image or image class.
In (Vizilter, Zheltov, 2012) where square of normalized morphological correlation for pair of regions.

Manifold learning, diffusion maps and shape matching
The application of heat kernels and diffusion maps to shape analysis is initially inspired by the manifold learning technique developed for nonlinear dimensionality reduction (NLDR).Most interesting NLDR techniques are the following: Isomap (Tenenbaum et al, 2000), Locally Linear Embedding (LLE) (Roweis et al 2000), Kernel Principle Component Analysis (Scholkopf, 1999), Laplacian Eigenmaps (Belkin, Niyogi, 2001), Hessian LLE (Donoho, 2003), Manifold Sculpting (Gashler, 2008) and some other.The terms "heat kernel" and "heat dissipation" were introduced in (Belkin, Niyogi, 2001) in the context of Laplacian Eigenmap.In this concept they play the role of some manifold shape characteristics.Based on this, authors of (Lafon, 2004), (Coifman, Lafon 2006), (Coifman et al, 2007) introduced and developed the theory of diffusion maps (DM).Let the manifold be described by some set of points X={x i } in a high-dimensional space.The solution of NLDR problem in DM approach has a following form:  Generate a neighborhood graph G.  Form a heat kernel (matrix of pairwise similarity weights) H = || h ij || using the rule.If i-th and j-th points are connected in G, then Normalize the heat kernel, and obtain the diffusion kernel Select the scale parameter t, and form the t-degree diffusion matrix P t , t1. Compute the spectral decomposition of P t with eigenvalues and corresponding eigenfunctions {ψ i }  Map the data to the low-dimensional vector space via selection of l eigenvalues and forming new coordinates based on corresponding eigenfunctions:


The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-3, 2014ISPRS Technical Commission III Symposium, 5 -7 September 2014, Zurich, Switzerland This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XL-3-357-2014 The normalization of heat kernel provides here both the linear smoothing diffusion operator P t , and its interpretation in terms of Markov chain with transition matrix P t for t steps of random walking.Thus, the Euclidean distance in this new space can be interpreted as a probability of point-to-point transition in t steps of this random walking.Such distance is called a diffusion distance, and it is very popular now in the area of high dimensional data analysis and machine learning.Being inspired by NLDR task, the DM approach was later successfully applied for other types of data analysis problems.In particular, the following image restoration technique was outlined in (Coifman et al, 2007).Let I(p) be a 2D image with p=(x p ,y p ).And let pixels be described by some feature vector v(p).Then for given >0 the diffusion kernel is defined as Let's note that idea of linear smoothing with adaptive kernel is utilized in different image restoration techniques, for example in (Takeda et al, 2007), (Milanfar, 2013).But the DM approach provides a unified way for shape description, shape comparison and shape-based image restoration.Let's look at the data shape matching techniques based on heat kernels, diffusion maps and their spectral features.A number of such approaches including heat kernel signature (HKS), heat kernel spectrum (set of eigenvalues), heat kernel signature distribution (HKSD) and heat trace (HT) are proposed and developed in (Sun et al, 2009), (de Goes et al, 2008), (Lieu, Saito, 2008), (Reuter et al, 2006) and other papers.In the brilliant paper (Memoli, 2011) the overview, classification and unified mathematical description of these techniques are given.Moreover, the special Gromov-Wasserstein distances are proposed for shape matching based on such spectral characteristics, and the stability (robustness) of such distances and matching procedures are theoretically proved (Memoli, 2011).These theoretical results are supported by impressive experiments with 3D models and real data collections.These results allow supposing, that the combination of DM and MIA approaches could create some effective 2D image and shape matching tools.

DIFFUSION MORPHOLOGY AND IMAGE MATCHING
In the first part of this section we propose a new generalized morphological framework based on diffusion shape models.In the second part the corresponding new morphological image and shape matching technique is described.

Generalized diffusion morphology
Let the image be a 2D function where Rset of real numbers, R 2image plane, rectangular frame region of image plane.Images are elements of Hilbert space L 2 () with scalar product (f,g) and norm || f || = (f,f) 1/2 .Let's introduce the generalized definitions of basic MIA notions by substitution of "mosaic shape" to "diffusion shape" in order to obtain the relaxed version of MIA.The relational model of diffusion shape F for image f is a heat kernel h F (x,y,u,v): , and the unique basic similarity measurement function  exists providing The operator model of diffusion shape F is a diffusion operator P F with normalized diffusion kernel p F (x,y,u,v), such that So, any relational model h F corresponds to the unique operator model with normalized kernel p F :


The spatial model of diffusion shape F for image f(x,y) with precision n is an eigenspace of diffusion operator where { 1 ,…, n }n first eigenvalues, а { 1 (x,y),…, n (x,y)}n first eigenfunctions of morphological diffusion operator P F : It is easy to see, that in particular case of heat kernel of the form , , , ( * it expresses the binary relation "points of equal values", and the diffusion-based morphological definitions stated above degrade to the MIA definitions given in section 2.2.Really, in this case the similarity relation  becomes an equivalent relation on the image frame points splitting them to a set of non-overlapping regions F={F 1 ,…,F n }, where n is a number of regions of frame tessellation F. Hence, the morphological diffusion operator becomes a morphological projector P F = P F P F with n 1-valued eigenvalues and characteristic (support) functions of n regions as eigenfunctions.So, the diffusion shape in this particular case becomes a mosaic shape F of the form (3).And for any g(x,y)L 2 () the diffusion filtering P F becomes a morphological projection onto the mosaic shape F of the classic form (4). Thus, the Pytiev MIA approach is a particular case of generalized diffusion morphology described in this section.So, all shape analysis schemes and tools of MIA can be recovered on this wider basis just using the diffusion operator instead of Pytiev morphological projector.
Let's note that the formal definition of diffusion morphology can start directly from a heat kernel of the classic form where d ijbasic distance between i-th and j-th points of discrete digital image.In this formulation of diffusion morphology the MIA case corresponds to a special selection of basic distance as a discrete distance by pixel value: Then the heat kernel takes a form h ij = h F *(x i ,y i ,u j ,v j ) of binary relation "points of equal values" that forces the transformation of diffusion operator to the mosaic projector.The final note here is that different choice of descriptors and metrics for pixel or/and region comparison provides the design of different diffusion morphologies with different semantic properties.Additionally, one can use the diffusion filter P t with scale parameter t for morphological scale space analysis.Thus, this generalized morphological approach provides more information about the image shape than just the information about shape of frame tessellation exploited by MI or original MIA.And if this new information is robust relative to noise, it will support the higher quality of matching.

Image and shape matching based on diffusion morphology
In Pytiev morphology the comparison of image g(x,y) and shape of image f(x,y) is performed using the normalized morphological correlation coefficient of the following form and this coefficient satisfies the property K M (f,F)=1 due to the fact that f = P F f.The diffusion morphological operator is not a projector, but it is a smoothing filter with || P F f ||  || f ||.Nevertheless, it is natural to suppose that the smoothing power of P F will be essentially less for images with similar shapes than for images with different shapes.So, the morphological diffusion correlation coefficient (MDCC) is defined as a ratio of Pytiev morphological coefficients , where K M (f,F) describes the power of self-smoothing of f by F. Let's note that MDCC is a correct generalization of MCC, because in case of projective morphology || f || = || P F f || and K MD (g,F)=K M (g,F).As in MCC, for elimination of non-informative part of image brightness images should be normalized before comparison: where P O fmorphological filtering of image f by the "empty" diffusion shape O.But in diffusion morphology such normalization is not the trivial subtraction of global mean value (like in Pytiev MIA).It is a subtraction of mean value in a sliding window determined by support neighborhood (effective size) of heat kernel.This subtraction preserves the local informative features only (in the corresponding scale of analysis).These informative elements of image g will be passed (if the shape G is similar to shape F) or extremely smoothed (if the shape G is essentially different) by the diffusion filter P F .This trick is called morphological image normalization (Fig. 2).).
If the effective size of heat kernel is small, such diffusion image-to-shape matching technique uses the local features only like the points-based and contour-based matching techniques.Fortunately, as we stated above, the theory of diffusion maps has a natural instrument for multi-scale data analysisparameter t (number of Markov random walking steps).The description of the image shape by the set of different scale diffusion operators {P t } allows performing the morphological scale-space analysis.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XL-3, 2014ISPRS Technical Commission III Symposium, 5 -7 September 2014, Zurich, Switzerland This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XL-3-357-2014 Finally, as well as for Pytiev correlation coefficient, we prefer to use the square of MDCC instead of MDCC, because both K M 2 (g,F) and K MD 2 (g,F) can be interpreted as statistical coefficient of determination between the model (shape F) and observed data (image g).So, the generalized morphological technique for image-to-shape matching is formed.

IMPLEMENTATION AND EXPERIMENTAL RESULTS
In the first part of this section we describe the original fast implementation of diffusion filters and morphological diffusion operators.This implementation provides a possibility for realtime image processing.It was utilized in our comparative experiments with different shape matching techniques for TV and IR image matching task.These experiments and corresponding results are outlined in the second part of this section.

Fast implementation of morphological diffusion operators
The computation of diffusion filtering with heat kernel of the form where v(p) is a neighborhood of image point p, is an extremely time-consuming procedure even for reasonable size of v(p).In this paper we propose to substitute such computationally unpleasant descriptors by the new type of point feature descriptors -iLBP (intensity + LBP): where m(p)mean value of v(p); LBP(p)threshold LBP (Ahonen, 2004) for v(p).The local binary pattern (LBP) is calculated here as an 8-bit vector for each pixel p based on a comparison of its value and values of its 8 nearest neighbors.If the value of neighbor pixel is less than the value of central pixel and the difference between them is greater than threshold, then the corresponding bit is set to 1, otherwiseto 0. Correspondingly we substitute the original neighborhood matching metrics by the special iLBP matching metrics: where d ham -Hamming distance, βimportance weight for intensity part of iLBP.Local binary patterns are stored as bit fields, and the computation of Hamming distance is performed via bitwise XOR operation.The exponent is calculated using table values.Mean value in a sliding window is computed by a fast algorithm with sliding sum recalculation.Due to this, the usage of iLBP allows both increasing the computational speed and obtaining heat kernels very similar to original.The software implementation of this idea provides the processing time about 300 ms for image size 640×480 and PC with Core i7 -860 /8Gb/ GeForce GTX 680 in one thread.The application of CUDA technology (NVIDIA, 2014) additionally improves the performance of diffusion filtering.Our current CUDA-based implementation provides the processing time about 40 ms for image size 640×480 and the same PC configuration.Thus, it allows performing the diffusion filtering and correspondingly the proposed shape-based matching procedure in real-time.

TV and IR matching experiments and results
The scheme of experimental comparison of matching techniques over a set of real images is following.Some clear fragments of TV image ("etalon" fragments) are compared with all equalsized fragments of corresponding IR image without segmentation.TV fragments are segmented with n =4 levels using least square optimal segmentation procedure.For segmented image fragments f (TV) and non-segmented g (IR) at each fragment g position (x,y) following similarity measures are calculated: mutual information MI(F,G) (1); square of centered Pytiev morphological coefficient ).The quality of these measures is estimated by following statistics of 2D correlation function С(x,y)signal-to-noise ratio (SNR) and exceeding of first maximum to second one (E):    In some experiments MDCC provides SNR in 2-3 times better.
Examples 2 and 3 demonstrate the higher stability of proposed MDCC relative to noise in comparison with mutual information (1), Pytiev MCC (4) and MSEMCC (5).

CONCLUSION
This paper presents the new image-to-shape matching technique based on heat kernels and diffusion maps.
The theoretical contribution of the paper is a formulation of Diffusion Morphology.It is designed as a generalization of Pytiev morphological image analysis (MIA) and contains MIA as a particular case.In our experiments with airborne TV and IR imagery the proposed MDCC totally outperforms all compared matching scores by both statistics, especially, in SNR.
The shape of correlation field is essentially better (sharper) than for all concurrent coefficients.In particular, relative to the stateof-art mutual information approach this new morphological technique provides SNR in 2-3 times better.Thus, we can conclude that proposed MDCC is a best correlation score for TV and IR matching problem.
Figure 1.Examples of diffusion filtering for denoising of TV and IR imagesFig.1demonstratesthe examples of TV and IR images diffusion filtering.In this case, like in paper(Coifman et al, 2007), v(p) is just a 5×5 vector of grayscale image values in a 5×5 pixel neighborhood.Both for TV and for IR image the shape was preserved and the noise was essentially removed.The reason of such success in noise suppression is an adaptive smoothing with high kernel weights for similar neighbors and low or zero weights for dissimilar.Let's note that idea of linear smoothing with adaptive kernel is utilized in different image restoration techniques, for example in(Takeda et al, 2007),(Milanfar, 2013).But the DM approach provides a unified way for shape description, shape comparison and shape-based image restoration.Let's look at the data shape matching techniques based on heat kernels, diffusion maps and their spectral features.A number of such approaches including heat kernel signature (HKS), heat kernel spectrum (set of eigenvalues), heat kernel signature distribution (HKSD) and heat trace (HT) are proposed and developed in(Sun et al, 2009), (de Goes et al, 2008),(Lieu, Saito, 2008),(Reuter et al, 2006) and other papers.In the brilliant paper(Memoli, 2011) the overview, classification and unified mathematical description of these techniques are given.Moreover, the special Gromov-Wasserstein distances are

Figure 2 .
Figure 2. Example of morphological image normalization: a) image f; b) image g; c) image f normalized by self-shape F ( f P f P O F 

Table 3
Numeric data for TV-IR matching (Example 3).