Scale-Space Theories for Scalar and Vector Images
Scale-Space Theories for Scalar and Vector Images
Luc M.J. Florack Utrecht University, PO Box 80010 NL-3584 CD Utrecht, The Netherlands [email protected] http://www.math.uu.nl/people/florack/index.html
Abstract. We define mutually consistent scale-space theories for scalar and vector images. Consistency pertains to the connection between the already established scalar theory and that for a suitably defined scalar field induced by the proposed vector scale-space. We show that one is compelled to reject the Gaussian scale-space paradigm in certain cases when scalar and vector fields are mutually dependent. Subsequently we investigate the behaviour of critical points of a vectorvalued scale-space image—i.e. points at which the vector field vanishes— as well as their singularities and unfoldings in linear scale-space.
1
Introduction
It has been argued by Koenderink [8] that a C21 (IRn×IR+ ) function1 u(x; s)—with scale parameter s ∈ IR+ —is a reasonable choice for a multiscale representation of a scalar image f (x) if, at the location of spatial extrema, us ∆u > 0 . Koenderink proposed to take the simplest instance of such a representation, which led him to consider the linear heat equation2 , (
∂s u = ∆u , lim u = f ,
(1)
s→0
1
2
A function u : IRn×IR+ → IR is in C21 (IRn×IR+ ) if it is twice continuously differentiable with respect to x ∈ IRn and once with respect to s ∈ IR+ . If u 6∈ C21 (IRn ×IR+ ) then there are other possibilities, cf. Van den Boomgaard et al. [1, 2]. Note that ∆u 6= 0 at a spatial extremum.
M. Kerckhove (Ed.): Scale-Space 2001, LNCS 2106, pp. 193–204, 2001. c Springer-Verlag and IEEE/CS 2001
194
Luc M.J. Florack
as the generating equation for a scale-space representation. Thus the argument pertains to the behaviour of extrema, and, in the case of the heat equation, basically boils down to the (strong) maximum principle 3 [4]. In this article we wish to generalise the scale-space paradigm to vector-valued images. A nontrivial demand to be reckoned with is mutual consistency: If the vector field implies the existence of a scalar field, as is for instance the case in a vector space endowed with a scalar product, then a multiscale representation of the latter must be consistent with the one induced by a multiscale representation of the former. We show that this implies that one is compelled to reject the Gaussian scale-space paradigm in certain cases. Subsequently we concentrate on linear vector scale-space representations and investigate the behaviour of its critical points—i.e. points at which the vector field vanishes—as well as their singularities and unfoldings in scale-space.
2
Theory
2.1
Vector Scale-Space
Once one appreciates that the addition of first order terms on the right hand side of Eq. (1) does not violate Koenderink’s causality principle, it becomes relatively straightforward to generalise the causality argument to vector images in such a way that consistency is manifest. That is to say, the scale-space of an arbitrary scalar field generated by the vector field (e.g. its magnitude) should be compatible with the already established theory for scalar images. In order to appreciate why first order terms are sometimes necessary, consider the following example. Example 1. Suppose we would require the magnitude of a suitably defined multiscale vector-valued image to satisfy the linear Eq. (1), i.e. without a first order nonlinearity, then we are confronted with a dilemma: As the norm of the highresolution vector image (at scale s = 0, say) is a positive function4 with isolated zeros only, its finite-scale representation must be everywhere strictly positive as a consequence of the maximum principle. Thus wherever defined, the scale-space representation of the vector field itself can have no critical points—spatial loci at which the vector field’s components vanish for fixed scale s > 0—not even at 3
4
It should be stressed that this does not imply that extrema cannot be created as scale increases, as is sometimes claimed in the literature, cf. Simmons et al. [10]. Damon [3] has shown that such creations do in fact occur generically, although they are typically outnumbered by annihilations. R A function f is called positive if f (x) ≥ 0 for almost all x ∈ IRn and f (x) dx 6= 0.
Scale-Space Theories for Scalar and Vector Images
195
“infinitesimal” scales. Elsewhere [6] it is explained that if the components v and w satisfy the following coupled system of p.d.e.’s 1 ∂s v = ∆v − (w∇v − v∇w) · ∇w v 2 + w2 1 ∂s w = ∆w − (v∇w − w∇v) · ∇v . v 2 + w2 √ then u = v 2 +w2 indeed satisfies the linear Eq. (1). The problem alluded to above is reflected in the denominator occuring on the right hand sides. In fact, these equations do not admit an initial condition with critical points. The scaling behaviour of the vector field’s critical points outlined in the example is a highly undesirable situation, and so we must reject Eq. (1) as a viable scalespace representation for the norm of a vector field (or any other scalar induced by the vector field for that matter). Critical points of a vector field should be able to survive a range of physical scales. The Poincar´e-Hopf theorem [11] bears witness to the fact that in spite of their zero measure they are not something “infinitesimal”. The way out of the dilemma is of course to relax the constraint on scalars induced by the vector field. None of these should be subjected to the linear heat equation for reasons amply discussed. Indeed, insisting on Koenderink’s scale causality demand does not compel us to assume the existence of any scalar satisfying Eq. (1), notably linearity, as this is only a special case of a more general class. Only if the scalar image is not some quantity derived from another physical observable, i.c. a vector field, the simplifying assumption of linearity may be justified, as there is no additional consistency demand to be met in that case. Vice versa, if we are given a scalar image and (thus) adopt Eq. (1) as our paradigm, the equations for the induced vector field v = ∇u are obviously linear, too: ( ∂s v = ∆v , (2) lim v = v0 , s→0
in which v0 is the high resolution gradient image. We may subsequently adopt Eq. (2) as the paradigm for vector fields in general, i.e. beyond the class of gradient fields. Apart from being the most straightforward thing to do in the first place, it is easy to verify that this is indeed consistent with the scalar theory, provided we refrain from the linearity condition √ for vector-induced scalar fields. Indeed, any function of the magnitude ρ = v · v satisfies an admissible semilinear diffusion equation. Let us consider some examples. Example 2. In two dimensions we may introduce a complex scalar field, u = v+iw ∈ C (respectively f = v0 +iw0 ), where (v, w) are the Cartesian components
196
Luc M.J. Florack
of v. The form of Eq. (2) then formally reduces to that Eq. (1). If we write u = % eiφ then we have ∂s % = ∆% − %∇φ · ∇φ , 2 ∂s φ = ∆φ + ∇% · ∇φ , % lim (%, φ) = (%0 , φ0 ) . s→0
The first order terms in the equation for % prevent the magnitude from becoming instantly everywhere positive, as opposed to the case in which they are absent. Thus critical points in the initial condition (%0 (xc ) = 0) are able to survive a finite amount of blur, as it should. Note also that the weighted angle ψ = %φ satisfies an equation of the same form as that of %. Using this function instead of φ one may thus circumvent the singularity % = 0 (implying φ is ill-defined but ψ = 0) on the right hand side.
We have seen that linear equations for the vector image imply nonlinear equations for induced scalar images. The example of a scalar image and its gradient illustrates the possibility of linearity in both domains. The next example, finally, shows that it may also be necessary to add nonlinear terms to the defining equations for a multiscale vector field derived from a linear scalar image.
Example 3. Elsewhere [7, 9] a theory has been proposed for multiscale motion extraction consistent with the scale-space paradigm for the underlying scalar image u(x; s). It has been noted that the operationally defined multiscale motion field v(x; s) does not satisfy Eq. (2), but that the flux field j(x; s) = u(x; s)v(x; s) by construction does. From the fact that u and j satisfy the linear diffusion equation it follows that the motion field itself satisfies ∂ v = ∆v + 2 ∇u · ∇v , s u lim v = v0 . s→0
Here the occurence of u in the denominator poses no problem since it is manifestly positive at scales s > 0, cf. Eq. (1).
All the examples given indicate that one must be cautious about which scalar or vector field one should subject to linear equations, Eq. (1) respectively Eq. (2), at least in cases where mutual dependencies exist (vector field versus magnitude field, scalar field versus motion field, et cetera). Otherwise a linear equation is the natural choice a priori.
Scale-Space Theories for Scalar and Vector Images
2.2
197
Behaviour of Critical Points
Next we study the behaviour of critical points given Eq. (2). Recall that a critical point is defined as the spatial locus xc at which v(xc ; s) = 0 for any fixed s. In the generic case, assuming v(x; s) is a Morse function for some given scale s, the critical points are all isolated. As scale increases each such point will move along a critical path. The situation is similar to that of critical points in a scalar image—in that context defined as spatial loci at which the image gradient vanishes—and can be analysed in a similar fashion [3, 5]. The situation for non-gradient fields is slightly more complicated, however. Theorem 1. Let J be the Jacobian matrix with components viµ = ∂i v µ , with row and column indices i, µ = 1, . . . , n, respectively. Furthermore, let J be its associated cofactor matrix, with components v jν —in which upper indices are now row indices—and det J the Jacobian determinant. A tangent to a parametrised critical path Γ : IR → IRn × IR+ : λ 7→ (x(λ); s(λ)) is then given by x˙ −J ∆v = s˙ det J in which a dot indicates differentiation with respect to λ. It is understood that the right hand side is evaluated at the location of the critical point. Note that ∆v = ∂s v. The definition of the cofactor matrix is reviewed in Appendix A. The essential property is J J = det J I. Proof. The derivation is essentially the same as presented elsewhere for critical paths in a scalar image [5]. The theorem is easily verified by noticing that the tangent must satisfy the n equations Jx˙ + ∆vs˙ = 0 , which follows by inspection of the first order Taylor terms of v at the location of a spatial critical point. Insertion of the right hand side in Theorem 1, using J J = det J I, readily shows that it indeed satisfies this constraint. It follows from Theorem 1 that as long as the Jacobian does not degenerate the critical path intersects the planes of constant scale transversally. However, as soon as det J = 0, the critical path becomes horizontal, and the critical point changes its character as it “reverses in scale”. This indicates either an annihilation or a creation—depending on the sign of curvature of the critical path at the singularity—of a pair of critical points with opposite Poincar´e indices. Theorem 1 applies to arbitrary dimensions. Let us scrutinise the 2-dimensional situation for simplicity. In that case the explicit form of the scale-space tangent
198
Luc M.J. Florack
to the critical path becomes ∂y v∆w − ∆v∂y w x˙ y˙ = ∆v∂x w − ∂x v∆w . ∂x v∂y w − ∂y v∂x w s˙ Note that the right hand side is just the outer product x˙ = ∇v ×∇w, in which ∇ is the scale-space gradient with components (∂x , ∂y , ∂s ), and in which v and w are the Cartesian components of the vector field. Let us now consider the classification of possible critical points. The eigenvalue equation for J is λ2 − tr J λ + det J = 0 , with tr J = div v. Analysis then shows that we may distinguish the following cases (it is understood that v = 0 in the points of interest): 1. If tr 2 J < 4 det J then λ1 , λ2 ∈ C\IR, λ1 = λ∗2 , say λ1,2 ≡ λ±iµ with λ, µ ∈ IR and µ 6= 0. 2. If tr 2 J = 4 det J then λ1 , λ2 ∈ IR, λ1 = λ2 ≡ λ 6= 0, so we have a symmetric nodal point: 2a. If λ < 0 we have a symmetric sink. 2b. If λ > 0 we have a symmetric source. 3. If tr 2 J > 4 det J then λ1 , λ2 ∈ IR, λ1 6= λ2 . The first and last case can be further subdivided as follows. 1a. If λ < 0 we have a spiral sink. 1b. If λ = 0 we have a central point. 1c. If λ > 0 we have a spiral source. 3a. If det J < 0 then sign λ1 = −sign λ2 , so that we have a saddle point. 3b. If det J = 0 then λ1 = 0, λ2 = div v or vice versa, which corresponds to a generic degeneracy. 3c. If, finally, det J > 0 then sign λ1 = sign λ2 , so that we have a nodal point. Still a further subdivision can be made (note that sign λ1,2 = sign div v): 3c1. If λ1 < 0 we have an asymmetric sink. 3c2. If λ1 > 0 we have an asymmetric source. Figure 1 shows all types of nondegenerate critical points. (All figures may be arbitrarily rotated.) The nature of the generic degeneracy (case 3b) is particularly interesting, for it is precisely at such a point in scale-space that the vector field changes topologically. In two dimensions the vanishing of one eigenvalue of the Jacobian matrix along
Scale-Space Theories for Scalar and Vector Images
199
Fig. 1. Nondegenerate Critical Points. Upper row, from left to right: spiral sink, centre point, spiral source. Lower row: saddle, (asymmetric) sink, (asymmetric) source. a critical path implies that the second eigenvalue must be real at the singularity, for if λ1 ∈ C\IR then λ2 = λ∗1 . Moreover, in the generic case under consideration, λ2 does not vanish at the same time as λ1 does. This follows from the fact that if λ1 = 0 then λ2 = λ1 + λ2 = tr J = div v. In the generic case the common points of the scale-space surface det J = 0 and the curve v = 0 are isolated. The probability that such a point lies within the surface div v = 0 is zero. This admits only one possible type of event, viz. a saddle point colliding with (i.e. annihilating with or emerging from) a nodal point, i.e. a sink or source node. Whether we have an annihilation or a creation depends on the (sign of the) curvature of the critical path at the singularity (v.i.). Figure 2 shows the vector field in the neighbourhood of a degenerate critical point.
Fig. 2. Generic Degeneracy (Unperturbed).
The above analysis implies that a spiral node (case 1a or 1c.) is stable in the sense that it can only cease to exist by first transforming into a critical point of
200
Luc M.J. Florack
type 3c1, respectively 3c2, via a corresponding symmetric nodal point, i.e. case 2, after which an annihilation with a saddle point may occur. This sequence of events corresponds to a pair of complex conjugate eigenvalues approaching the real axis and scattering off in horizontal direction after collision. As soon as one of the eigenvalues reaches 0 we get an event of type 3b. Vice versa, spiral nodes cannot be created spontaneously; the abovesketched sequence must be traversed in opposite direction, starting out from a creation of a saddle and a nodal point, whereby the latter turns into a spiral node via an intermediate symmetric nodal point. See Figure 3.
c i 0
1
Fig. 3. Mutation diagram of (the eigenvalues associated with) a critical point in the C-plane, showing a spiral node (nonreal complex conjugate points approaching the IR-axis) transforming into an asymmetric nodal point (pair of distinct real points moving along the IR-axis) via a symmetric nodal point (scatterpoint on the IR-axis). As soon as one of the scattered points reaches the origin we obtain a degeneracy. The velocities of corresponding points are mirror-symmetric relative to the real axis.
Theorem 1 implies that at the scale-space location of a catastrophe at which a saddle and a nodal point collide the tangent to the critical path is horizontal (s˙ = 0). In order to distinguish between creation and annihilation events involving a saddle and a nodal point we need a higher order local analysis. Example 4. Consider the vector field germs, together with their perturbations, 2 2 x + 2s x − 2y 2 − 2s 0 0 + + and vc (x, y, s) = . va (x, y, s) = 0 −4xy y y The second term on each right hand side is the canonical form of a typical perturbation on a full scale-space neighbourhood of the fiducial origin, (x, y, s) = a creation event (0, 0, 0). The field va (x, y, s) captures an annihilation, vc (x, y, s) √ : (x, y, s) = (± −2s, 0, s), s ≤ at the origin. The critical paths are given by Γ a √ 0, and Γc : (x, y, s) = (± 2s, 0, s), s ≥ 0, respectively. Figure 4 illustrates the behaviour of these fields near the origin. These germs and their perturbations are analogous to those describing the unfoldings of corresponding annihilation and creation events in scalar images, cf. Damon [3].
Scale-Space Theories for Scalar and Vector Images
201
Fig. 4. Unfolding of Generic Singularities. Upper row: annihilation event. Lower row: creation event. Scale (resolution) increases (decreases) from left to right.
Critical points of a vector field all coincide with global minima of the vector field’s magnitude, but the latter generally has additional critical points (maxima, saddles and local minima), viz. those points where v · ∇v = 2 JT v = 0 and v 6= 0, i.e. where the vector field itself happens to be nonvanishing and orthogonal to its gradient (note that this implies degeneracy of the Jacobian). Inspection of the gradients and Hessian matrices of the scalar images ua,c (x, y, s) = kva,c (x, y, s)k2 as a function of scale reveals that at a singularity of the vector field there is an intersection of critical paths with the shape of a pitchfork pointing downward in the annihilation event and upward in the creation event. In the first case a saddle collides with two global minima (the critical points of va ) at (x, y, s) = (0, 0, 0), after which a single local minimum emerges. In the second case a local minimum spawns two global minima (the critical points of vc ) at (x, y, s) = (0, 0, 0), together with a saddle. In both cases the Poincar´e index of a spatial region containing (only) the critical points under consideration remains invariant, as it should: Figure 5. (In this argument “local minima” may be replaced by “local maxima” depending on the case of interest.) The critical points of the vector field’s magnitude that are not critical points of the vector field itself define critical paths in scale-space implicitly defined by def
q(x; s) = v(x; s) · ∇v(x; s) = 0 .
202
Luc M.J. Florack saddle
local min/max global min
global min global min
global min local min/max
saddle
Fig. 5. Unfolding of the generic singularities of vector field’s magnitude that correspond to those of the vector field itself. Left: annihilation event. Right: creation event. Scale (resolution) increases (decreases) upward. The singularity always involves two global minima as well as a saddle/local minimum/maximum pair. In the two-dimensional case at hand we have two equations, one for each component of v = (v, w), and three unknowns, (x; s) = (x, y, s). Consequently we again expect to find critical paths in scale-space, which can be analysed as previously. If we expand the defining equation to second order we obtain q(x; s) = q + ∇q · x +
1 T x · ∇∇T q · x + ∂s q s + h.o.t. 2
in which the coefficients on the right hand side are partial derivatives evaluated at the origin. Setting the left hand side equal to zero, and defining the matrix def
Q = ∇q ,
(3)
with components qij = ∂i (v · ∂j v), it follows that a scale-space tangent to the critical path is given by the following theorem, the derivation of which is completely analogous to that of Theorem 1. Theorem 2. A tangent to a parametrised critical path Γ : IR → IRn ×IR+ : λ 7→ (x(λ); s(λ)) through a local critical point of the magnitude field u(x; s) = kv(x; s)k is given by x˙ −Q ∂s q = s˙ det Q in which a dot indicates differentiation with respect to the curve parameter λ, and in which Q is the matrix as defined in Eq. (3). It is understood that the right hand side is evaluated at the location of the critical point. Note that ∂s q = ∆v · ∇v+v · ∇∆v. The critical path is transversal to s-planes if and only if det Q 6= 0, in which case we may use λ = s itself as a valid curve parameter. If det Q = 0 the curve becomes horizontal. Assuming that v 6= 0 at that singularity (the other case has been discussed), we have either an annihilation or a creation of local critical points, one of which is a saddle, the other a maximum or local minimum. The
Scale-Space Theories for Scalar and Vector Images
203
criterion is again the sign of the critical path’s curvature at the singularity. This case falls within the scope of Damon’s analysis [3]. See Figure 6.
local min/max
local min/max
saddle
saddle
Fig. 6. Unfolding of the generic singularities of vector field’s magnitude that do not involve critical points of the vector field itself. Left: annihilation event. Right: creation event. Scale (resolution) increases (decreases) upward. The singularity always involves a saddle/local minimum/maximum pair. The situation is similar to that encountered in scalar images.
A
Cofactor Matrix
Let A be a square n × n matrix with components aµν . Then we define the transposed cofactor matrix A as follows. In order to obtain the matrix entry aµν skip the µ-th column and ν-th row of A, evaluate the determinant of the resulting submatrix, and multiply by (−1)µ+ν (“checkerboard pattern”). Using tensor index notation, µν def
A
=
1 εµµ1 ...µn−1 ενν1 ...νn−1 Aµ1 ν1 . . . Aµn−1 νn−1 , (n − 1)!
in which εµ1 ...µn is the spatial L´evi-Civita tensor, defined as the normalised, completely antisymmetric symbol with ε1...n = 1. An important property of the cofactor matrix is A A = A A = det A I , in which I is the n×n identity matrix. Thus if a matrix is invertible, its cofactor matrix equals its inverse times its determinant. Unlike the inverse matrix, however, the cofactor matrix is always well-defined, because its coefficients are homogeneous polynomials of degree n−1 relative to the coefficients of the original matrix. From the definition it follows that I = I. Moreover, r = r for any number r ∈ IR, and if B = A then A = B. The following lemma has been used in the foregoing text.
204
Luc M.J. Florack
Lemma 1. For any matrix A we have ∇ det A = tr A∇A . Proof. For invertible matrices the proof goes as follows. Write det A = exp tr ln A. Taking the gradient yields ∇ det A = ∇(tr ln A) exp(tr ln A) = tr (∇ ln A) det A = tr (Ainv ∇A) det A = tr (A ∇A). If A is not invertible, consider the regularised, invertible matrix Aε = A+ε I instead and apply the theorem. Left and right hand sides are polynomials in ε; taking the limit ε → 0 establishes the result.
References [1] R. van den Boomgaard. The morphological equivalent of the Gauss convolution. Nieuw Archief voor Wiskunde, 10(3):219–236, November 1992. [2] R. van den Boomgaard and A. W. M. Smeulders. The morphological structure of images, the differential equations of morphological scale-space. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(11):1101–1113, November 1994. [3] J. Damon. Local Morse theory for solutions to the heat equation and Gaussian blurring. Journal of Differential Equations, 115(2):368–401, January 1995. [4] L. C. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics. American Mathematical Society, Providence, Rhode Island, 1998. [5] L. Florack and A. Kuijper. The topological structure of scale-space images. Journal of Mathematical Imaging and Vision, 12(1):65–79, February 2000. [6] L. M. J. Florack. Non-linear scale-spaces isomorphic to the linear case. In B. K. Ersbøll and P. Johansen, editors, Proceedings of the 11th Scandinavian Conference on Image Analysis (Kangerlussuaq, Greenland, June 7–11 1999), volume 1, pages 229–234, Lyngby, Denmark, 1999. [7] L. M. J. Florack, W. J. Niessen, and M. Nielsen. The intrinsic structure of optic flow incorporating measurement duality. International Journal of Computer Vision, 27(3):263–286, May 1998. [8] J. J. Koenderink. The structure of images. Biological Cybernetics, 50:363–370, 1984. [9] M. Nielsen and O. F. Olsen. The structure of the optic flow field. In H. Burkhardt and B. Neumann, editors, Proceedings of the Fifth European Conference on Computer Vision (Freiburg, Germany, June 1998), volume 1407 of Lecture Notes in Computer Science, pages 281–287. Springer-Verlag, Berlin, 1998. [10] A. Simmons, S. R. Arridge, P. S. Tofts, and G. J. Barker. Application of the extremum stack to neurological MRI. IEEE Transactions on Medical Imaging, 17(3):371–382, June 1998. [11] M. Spivak. Differential Geometry, volume 1. Publish or Perish, Berkeley, 1975.
Luc M.J. Florack Utrecht University, PO Box 80010 NL-3584 CD Utrecht, The Netherlands [email protected] http://www.math.uu.nl/people/florack/index.html
Abstract. We define mutually consistent scale-space theories for scalar and vector images. Consistency pertains to the connection between the already established scalar theory and that for a suitably defined scalar field induced by the proposed vector scale-space. We show that one is compelled to reject the Gaussian scale-space paradigm in certain cases when scalar and vector fields are mutually dependent. Subsequently we investigate the behaviour of critical points of a vectorvalued scale-space image—i.e. points at which the vector field vanishes— as well as their singularities and unfoldings in linear scale-space.
1
Introduction
It has been argued by Koenderink [8] that a C21 (IRn×IR+ ) function1 u(x; s)—with scale parameter s ∈ IR+ —is a reasonable choice for a multiscale representation of a scalar image f (x) if, at the location of spatial extrema, us ∆u > 0 . Koenderink proposed to take the simplest instance of such a representation, which led him to consider the linear heat equation2 , (
∂s u = ∆u , lim u = f ,
(1)
s→0
1
2
A function u : IRn×IR+ → IR is in C21 (IRn×IR+ ) if it is twice continuously differentiable with respect to x ∈ IRn and once with respect to s ∈ IR+ . If u 6∈ C21 (IRn ×IR+ ) then there are other possibilities, cf. Van den Boomgaard et al. [1, 2]. Note that ∆u 6= 0 at a spatial extremum.
M. Kerckhove (Ed.): Scale-Space 2001, LNCS 2106, pp. 193–204, 2001. c Springer-Verlag and IEEE/CS 2001
194
Luc M.J. Florack
as the generating equation for a scale-space representation. Thus the argument pertains to the behaviour of extrema, and, in the case of the heat equation, basically boils down to the (strong) maximum principle 3 [4]. In this article we wish to generalise the scale-space paradigm to vector-valued images. A nontrivial demand to be reckoned with is mutual consistency: If the vector field implies the existence of a scalar field, as is for instance the case in a vector space endowed with a scalar product, then a multiscale representation of the latter must be consistent with the one induced by a multiscale representation of the former. We show that this implies that one is compelled to reject the Gaussian scale-space paradigm in certain cases. Subsequently we concentrate on linear vector scale-space representations and investigate the behaviour of its critical points—i.e. points at which the vector field vanishes—as well as their singularities and unfoldings in scale-space.
2
Theory
2.1
Vector Scale-Space
Once one appreciates that the addition of first order terms on the right hand side of Eq. (1) does not violate Koenderink’s causality principle, it becomes relatively straightforward to generalise the causality argument to vector images in such a way that consistency is manifest. That is to say, the scale-space of an arbitrary scalar field generated by the vector field (e.g. its magnitude) should be compatible with the already established theory for scalar images. In order to appreciate why first order terms are sometimes necessary, consider the following example. Example 1. Suppose we would require the magnitude of a suitably defined multiscale vector-valued image to satisfy the linear Eq. (1), i.e. without a first order nonlinearity, then we are confronted with a dilemma: As the norm of the highresolution vector image (at scale s = 0, say) is a positive function4 with isolated zeros only, its finite-scale representation must be everywhere strictly positive as a consequence of the maximum principle. Thus wherever defined, the scale-space representation of the vector field itself can have no critical points—spatial loci at which the vector field’s components vanish for fixed scale s > 0—not even at 3
4
It should be stressed that this does not imply that extrema cannot be created as scale increases, as is sometimes claimed in the literature, cf. Simmons et al. [10]. Damon [3] has shown that such creations do in fact occur generically, although they are typically outnumbered by annihilations. R A function f is called positive if f (x) ≥ 0 for almost all x ∈ IRn and f (x) dx 6= 0.
Scale-Space Theories for Scalar and Vector Images
195
“infinitesimal” scales. Elsewhere [6] it is explained that if the components v and w satisfy the following coupled system of p.d.e.’s 1 ∂s v = ∆v − (w∇v − v∇w) · ∇w v 2 + w2 1 ∂s w = ∆w − (v∇w − w∇v) · ∇v . v 2 + w2 √ then u = v 2 +w2 indeed satisfies the linear Eq. (1). The problem alluded to above is reflected in the denominator occuring on the right hand sides. In fact, these equations do not admit an initial condition with critical points. The scaling behaviour of the vector field’s critical points outlined in the example is a highly undesirable situation, and so we must reject Eq. (1) as a viable scalespace representation for the norm of a vector field (or any other scalar induced by the vector field for that matter). Critical points of a vector field should be able to survive a range of physical scales. The Poincar´e-Hopf theorem [11] bears witness to the fact that in spite of their zero measure they are not something “infinitesimal”. The way out of the dilemma is of course to relax the constraint on scalars induced by the vector field. None of these should be subjected to the linear heat equation for reasons amply discussed. Indeed, insisting on Koenderink’s scale causality demand does not compel us to assume the existence of any scalar satisfying Eq. (1), notably linearity, as this is only a special case of a more general class. Only if the scalar image is not some quantity derived from another physical observable, i.c. a vector field, the simplifying assumption of linearity may be justified, as there is no additional consistency demand to be met in that case. Vice versa, if we are given a scalar image and (thus) adopt Eq. (1) as our paradigm, the equations for the induced vector field v = ∇u are obviously linear, too: ( ∂s v = ∆v , (2) lim v = v0 , s→0
in which v0 is the high resolution gradient image. We may subsequently adopt Eq. (2) as the paradigm for vector fields in general, i.e. beyond the class of gradient fields. Apart from being the most straightforward thing to do in the first place, it is easy to verify that this is indeed consistent with the scalar theory, provided we refrain from the linearity condition √ for vector-induced scalar fields. Indeed, any function of the magnitude ρ = v · v satisfies an admissible semilinear diffusion equation. Let us consider some examples. Example 2. In two dimensions we may introduce a complex scalar field, u = v+iw ∈ C (respectively f = v0 +iw0 ), where (v, w) are the Cartesian components
196
Luc M.J. Florack
of v. The form of Eq. (2) then formally reduces to that Eq. (1). If we write u = % eiφ then we have ∂s % = ∆% − %∇φ · ∇φ , 2 ∂s φ = ∆φ + ∇% · ∇φ , % lim (%, φ) = (%0 , φ0 ) . s→0
The first order terms in the equation for % prevent the magnitude from becoming instantly everywhere positive, as opposed to the case in which they are absent. Thus critical points in the initial condition (%0 (xc ) = 0) are able to survive a finite amount of blur, as it should. Note also that the weighted angle ψ = %φ satisfies an equation of the same form as that of %. Using this function instead of φ one may thus circumvent the singularity % = 0 (implying φ is ill-defined but ψ = 0) on the right hand side.
We have seen that linear equations for the vector image imply nonlinear equations for induced scalar images. The example of a scalar image and its gradient illustrates the possibility of linearity in both domains. The next example, finally, shows that it may also be necessary to add nonlinear terms to the defining equations for a multiscale vector field derived from a linear scalar image.
Example 3. Elsewhere [7, 9] a theory has been proposed for multiscale motion extraction consistent with the scale-space paradigm for the underlying scalar image u(x; s). It has been noted that the operationally defined multiscale motion field v(x; s) does not satisfy Eq. (2), but that the flux field j(x; s) = u(x; s)v(x; s) by construction does. From the fact that u and j satisfy the linear diffusion equation it follows that the motion field itself satisfies ∂ v = ∆v + 2 ∇u · ∇v , s u lim v = v0 . s→0
Here the occurence of u in the denominator poses no problem since it is manifestly positive at scales s > 0, cf. Eq. (1).
All the examples given indicate that one must be cautious about which scalar or vector field one should subject to linear equations, Eq. (1) respectively Eq. (2), at least in cases where mutual dependencies exist (vector field versus magnitude field, scalar field versus motion field, et cetera). Otherwise a linear equation is the natural choice a priori.
Scale-Space Theories for Scalar and Vector Images
2.2
197
Behaviour of Critical Points
Next we study the behaviour of critical points given Eq. (2). Recall that a critical point is defined as the spatial locus xc at which v(xc ; s) = 0 for any fixed s. In the generic case, assuming v(x; s) is a Morse function for some given scale s, the critical points are all isolated. As scale increases each such point will move along a critical path. The situation is similar to that of critical points in a scalar image—in that context defined as spatial loci at which the image gradient vanishes—and can be analysed in a similar fashion [3, 5]. The situation for non-gradient fields is slightly more complicated, however. Theorem 1. Let J be the Jacobian matrix with components viµ = ∂i v µ , with row and column indices i, µ = 1, . . . , n, respectively. Furthermore, let J be its associated cofactor matrix, with components v jν —in which upper indices are now row indices—and det J the Jacobian determinant. A tangent to a parametrised critical path Γ : IR → IRn × IR+ : λ 7→ (x(λ); s(λ)) is then given by x˙ −J ∆v = s˙ det J in which a dot indicates differentiation with respect to λ. It is understood that the right hand side is evaluated at the location of the critical point. Note that ∆v = ∂s v. The definition of the cofactor matrix is reviewed in Appendix A. The essential property is J J = det J I. Proof. The derivation is essentially the same as presented elsewhere for critical paths in a scalar image [5]. The theorem is easily verified by noticing that the tangent must satisfy the n equations Jx˙ + ∆vs˙ = 0 , which follows by inspection of the first order Taylor terms of v at the location of a spatial critical point. Insertion of the right hand side in Theorem 1, using J J = det J I, readily shows that it indeed satisfies this constraint. It follows from Theorem 1 that as long as the Jacobian does not degenerate the critical path intersects the planes of constant scale transversally. However, as soon as det J = 0, the critical path becomes horizontal, and the critical point changes its character as it “reverses in scale”. This indicates either an annihilation or a creation—depending on the sign of curvature of the critical path at the singularity—of a pair of critical points with opposite Poincar´e indices. Theorem 1 applies to arbitrary dimensions. Let us scrutinise the 2-dimensional situation for simplicity. In that case the explicit form of the scale-space tangent
198
Luc M.J. Florack
to the critical path becomes ∂y v∆w − ∆v∂y w x˙ y˙ = ∆v∂x w − ∂x v∆w . ∂x v∂y w − ∂y v∂x w s˙ Note that the right hand side is just the outer product x˙ = ∇v ×∇w, in which ∇ is the scale-space gradient with components (∂x , ∂y , ∂s ), and in which v and w are the Cartesian components of the vector field. Let us now consider the classification of possible critical points. The eigenvalue equation for J is λ2 − tr J λ + det J = 0 , with tr J = div v. Analysis then shows that we may distinguish the following cases (it is understood that v = 0 in the points of interest): 1. If tr 2 J < 4 det J then λ1 , λ2 ∈ C\IR, λ1 = λ∗2 , say λ1,2 ≡ λ±iµ with λ, µ ∈ IR and µ 6= 0. 2. If tr 2 J = 4 det J then λ1 , λ2 ∈ IR, λ1 = λ2 ≡ λ 6= 0, so we have a symmetric nodal point: 2a. If λ < 0 we have a symmetric sink. 2b. If λ > 0 we have a symmetric source. 3. If tr 2 J > 4 det J then λ1 , λ2 ∈ IR, λ1 6= λ2 . The first and last case can be further subdivided as follows. 1a. If λ < 0 we have a spiral sink. 1b. If λ = 0 we have a central point. 1c. If λ > 0 we have a spiral source. 3a. If det J < 0 then sign λ1 = −sign λ2 , so that we have a saddle point. 3b. If det J = 0 then λ1 = 0, λ2 = div v or vice versa, which corresponds to a generic degeneracy. 3c. If, finally, det J > 0 then sign λ1 = sign λ2 , so that we have a nodal point. Still a further subdivision can be made (note that sign λ1,2 = sign div v): 3c1. If λ1 < 0 we have an asymmetric sink. 3c2. If λ1 > 0 we have an asymmetric source. Figure 1 shows all types of nondegenerate critical points. (All figures may be arbitrarily rotated.) The nature of the generic degeneracy (case 3b) is particularly interesting, for it is precisely at such a point in scale-space that the vector field changes topologically. In two dimensions the vanishing of one eigenvalue of the Jacobian matrix along
Scale-Space Theories for Scalar and Vector Images
199
Fig. 1. Nondegenerate Critical Points. Upper row, from left to right: spiral sink, centre point, spiral source. Lower row: saddle, (asymmetric) sink, (asymmetric) source. a critical path implies that the second eigenvalue must be real at the singularity, for if λ1 ∈ C\IR then λ2 = λ∗1 . Moreover, in the generic case under consideration, λ2 does not vanish at the same time as λ1 does. This follows from the fact that if λ1 = 0 then λ2 = λ1 + λ2 = tr J = div v. In the generic case the common points of the scale-space surface det J = 0 and the curve v = 0 are isolated. The probability that such a point lies within the surface div v = 0 is zero. This admits only one possible type of event, viz. a saddle point colliding with (i.e. annihilating with or emerging from) a nodal point, i.e. a sink or source node. Whether we have an annihilation or a creation depends on the (sign of the) curvature of the critical path at the singularity (v.i.). Figure 2 shows the vector field in the neighbourhood of a degenerate critical point.
Fig. 2. Generic Degeneracy (Unperturbed).
The above analysis implies that a spiral node (case 1a or 1c.) is stable in the sense that it can only cease to exist by first transforming into a critical point of
200
Luc M.J. Florack
type 3c1, respectively 3c2, via a corresponding symmetric nodal point, i.e. case 2, after which an annihilation with a saddle point may occur. This sequence of events corresponds to a pair of complex conjugate eigenvalues approaching the real axis and scattering off in horizontal direction after collision. As soon as one of the eigenvalues reaches 0 we get an event of type 3b. Vice versa, spiral nodes cannot be created spontaneously; the abovesketched sequence must be traversed in opposite direction, starting out from a creation of a saddle and a nodal point, whereby the latter turns into a spiral node via an intermediate symmetric nodal point. See Figure 3.
c i 0
1
Fig. 3. Mutation diagram of (the eigenvalues associated with) a critical point in the C-plane, showing a spiral node (nonreal complex conjugate points approaching the IR-axis) transforming into an asymmetric nodal point (pair of distinct real points moving along the IR-axis) via a symmetric nodal point (scatterpoint on the IR-axis). As soon as one of the scattered points reaches the origin we obtain a degeneracy. The velocities of corresponding points are mirror-symmetric relative to the real axis.
Theorem 1 implies that at the scale-space location of a catastrophe at which a saddle and a nodal point collide the tangent to the critical path is horizontal (s˙ = 0). In order to distinguish between creation and annihilation events involving a saddle and a nodal point we need a higher order local analysis. Example 4. Consider the vector field germs, together with their perturbations, 2 2 x + 2s x − 2y 2 − 2s 0 0 + + and vc (x, y, s) = . va (x, y, s) = 0 −4xy y y The second term on each right hand side is the canonical form of a typical perturbation on a full scale-space neighbourhood of the fiducial origin, (x, y, s) = a creation event (0, 0, 0). The field va (x, y, s) captures an annihilation, vc (x, y, s) √ : (x, y, s) = (± −2s, 0, s), s ≤ at the origin. The critical paths are given by Γ a √ 0, and Γc : (x, y, s) = (± 2s, 0, s), s ≥ 0, respectively. Figure 4 illustrates the behaviour of these fields near the origin. These germs and their perturbations are analogous to those describing the unfoldings of corresponding annihilation and creation events in scalar images, cf. Damon [3].
Scale-Space Theories for Scalar and Vector Images
201
Fig. 4. Unfolding of Generic Singularities. Upper row: annihilation event. Lower row: creation event. Scale (resolution) increases (decreases) from left to right.
Critical points of a vector field all coincide with global minima of the vector field’s magnitude, but the latter generally has additional critical points (maxima, saddles and local minima), viz. those points where v · ∇v = 2 JT v = 0 and v 6= 0, i.e. where the vector field itself happens to be nonvanishing and orthogonal to its gradient (note that this implies degeneracy of the Jacobian). Inspection of the gradients and Hessian matrices of the scalar images ua,c (x, y, s) = kva,c (x, y, s)k2 as a function of scale reveals that at a singularity of the vector field there is an intersection of critical paths with the shape of a pitchfork pointing downward in the annihilation event and upward in the creation event. In the first case a saddle collides with two global minima (the critical points of va ) at (x, y, s) = (0, 0, 0), after which a single local minimum emerges. In the second case a local minimum spawns two global minima (the critical points of vc ) at (x, y, s) = (0, 0, 0), together with a saddle. In both cases the Poincar´e index of a spatial region containing (only) the critical points under consideration remains invariant, as it should: Figure 5. (In this argument “local minima” may be replaced by “local maxima” depending on the case of interest.) The critical points of the vector field’s magnitude that are not critical points of the vector field itself define critical paths in scale-space implicitly defined by def
q(x; s) = v(x; s) · ∇v(x; s) = 0 .
202
Luc M.J. Florack saddle
local min/max global min
global min global min
global min local min/max
saddle
Fig. 5. Unfolding of the generic singularities of vector field’s magnitude that correspond to those of the vector field itself. Left: annihilation event. Right: creation event. Scale (resolution) increases (decreases) upward. The singularity always involves two global minima as well as a saddle/local minimum/maximum pair. In the two-dimensional case at hand we have two equations, one for each component of v = (v, w), and three unknowns, (x; s) = (x, y, s). Consequently we again expect to find critical paths in scale-space, which can be analysed as previously. If we expand the defining equation to second order we obtain q(x; s) = q + ∇q · x +
1 T x · ∇∇T q · x + ∂s q s + h.o.t. 2
in which the coefficients on the right hand side are partial derivatives evaluated at the origin. Setting the left hand side equal to zero, and defining the matrix def
Q = ∇q ,
(3)
with components qij = ∂i (v · ∂j v), it follows that a scale-space tangent to the critical path is given by the following theorem, the derivation of which is completely analogous to that of Theorem 1. Theorem 2. A tangent to a parametrised critical path Γ : IR → IRn ×IR+ : λ 7→ (x(λ); s(λ)) through a local critical point of the magnitude field u(x; s) = kv(x; s)k is given by x˙ −Q ∂s q = s˙ det Q in which a dot indicates differentiation with respect to the curve parameter λ, and in which Q is the matrix as defined in Eq. (3). It is understood that the right hand side is evaluated at the location of the critical point. Note that ∂s q = ∆v · ∇v+v · ∇∆v. The critical path is transversal to s-planes if and only if det Q 6= 0, in which case we may use λ = s itself as a valid curve parameter. If det Q = 0 the curve becomes horizontal. Assuming that v 6= 0 at that singularity (the other case has been discussed), we have either an annihilation or a creation of local critical points, one of which is a saddle, the other a maximum or local minimum. The
Scale-Space Theories for Scalar and Vector Images
203
criterion is again the sign of the critical path’s curvature at the singularity. This case falls within the scope of Damon’s analysis [3]. See Figure 6.
local min/max
local min/max
saddle
saddle
Fig. 6. Unfolding of the generic singularities of vector field’s magnitude that do not involve critical points of the vector field itself. Left: annihilation event. Right: creation event. Scale (resolution) increases (decreases) upward. The singularity always involves a saddle/local minimum/maximum pair. The situation is similar to that encountered in scalar images.
A
Cofactor Matrix
Let A be a square n × n matrix with components aµν . Then we define the transposed cofactor matrix A as follows. In order to obtain the matrix entry aµν skip the µ-th column and ν-th row of A, evaluate the determinant of the resulting submatrix, and multiply by (−1)µ+ν (“checkerboard pattern”). Using tensor index notation, µν def
A
=
1 εµµ1 ...µn−1 ενν1 ...νn−1 Aµ1 ν1 . . . Aµn−1 νn−1 , (n − 1)!
in which εµ1 ...µn is the spatial L´evi-Civita tensor, defined as the normalised, completely antisymmetric symbol with ε1...n = 1. An important property of the cofactor matrix is A A = A A = det A I , in which I is the n×n identity matrix. Thus if a matrix is invertible, its cofactor matrix equals its inverse times its determinant. Unlike the inverse matrix, however, the cofactor matrix is always well-defined, because its coefficients are homogeneous polynomials of degree n−1 relative to the coefficients of the original matrix. From the definition it follows that I = I. Moreover, r = r for any number r ∈ IR, and if B = A then A = B. The following lemma has been used in the foregoing text.
204
Luc M.J. Florack
Lemma 1. For any matrix A we have ∇ det A = tr A∇A . Proof. For invertible matrices the proof goes as follows. Write det A = exp tr ln A. Taking the gradient yields ∇ det A = ∇(tr ln A) exp(tr ln A) = tr (∇ ln A) det A = tr (Ainv ∇A) det A = tr (A ∇A). If A is not invertible, consider the regularised, invertible matrix Aε = A+ε I instead and apply the theorem. Left and right hand sides are polynomials in ε; taking the limit ε → 0 establishes the result.
References [1] R. van den Boomgaard. The morphological equivalent of the Gauss convolution. Nieuw Archief voor Wiskunde, 10(3):219–236, November 1992. [2] R. van den Boomgaard and A. W. M. Smeulders. The morphological structure of images, the differential equations of morphological scale-space. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(11):1101–1113, November 1994. [3] J. Damon. Local Morse theory for solutions to the heat equation and Gaussian blurring. Journal of Differential Equations, 115(2):368–401, January 1995. [4] L. C. Evans. Partial Differential Equations, volume 19 of Graduate Studies in Mathematics. American Mathematical Society, Providence, Rhode Island, 1998. [5] L. Florack and A. Kuijper. The topological structure of scale-space images. Journal of Mathematical Imaging and Vision, 12(1):65–79, February 2000. [6] L. M. J. Florack. Non-linear scale-spaces isomorphic to the linear case. In B. K. Ersbøll and P. Johansen, editors, Proceedings of the 11th Scandinavian Conference on Image Analysis (Kangerlussuaq, Greenland, June 7–11 1999), volume 1, pages 229–234, Lyngby, Denmark, 1999. [7] L. M. J. Florack, W. J. Niessen, and M. Nielsen. The intrinsic structure of optic flow incorporating measurement duality. International Journal of Computer Vision, 27(3):263–286, May 1998. [8] J. J. Koenderink. The structure of images. Biological Cybernetics, 50:363–370, 1984. [9] M. Nielsen and O. F. Olsen. The structure of the optic flow field. In H. Burkhardt and B. Neumann, editors, Proceedings of the Fifth European Conference on Computer Vision (Freiburg, Germany, June 1998), volume 1407 of Lecture Notes in Computer Science, pages 281–287. Springer-Verlag, Berlin, 1998. [10] A. Simmons, S. R. Arridge, P. S. Tofts, and G. J. Barker. Application of the extremum stack to neurological MRI. IEEE Transactions on Medical Imaging, 17(3):371–382, June 1998. [11] M. Spivak. Differential Geometry, volume 1. Publish or Perish, Berkeley, 1975.
