1 Introduction

Mathematical morphology combines nonlinear image processing and filtering techniques based on ideas from (random) set theory, topology, and stochastic geometry. It is widely applied for the analysis of spatial structures, e.g. from geology, biology, or materials science. The foundations of mathematical morphology are laid in the books by Matheron [1], Serra [2, 3], and Soille [4]. Sternberg [5] generalised mathematical morphology to numerical functions via the umbra method, i.e. the application of set morphology to the graph of the function. Ronse [6] and Goutsias et al. [7] derived a further generalisation to complete lattices, that is, partially ordered sets with the property that all subsets have a supremum and an infimum.

The progress of imaging methods from binary via grey-scale images to vector-valued images, such as colour or ultra-spectral images, requires modified morphological operators. In particular, diffusion tensor imaging (DTI) in neuroimaging [8] results in vector-valued images showing the diffusion of water molecules in the brain. These images can be used to better understand neural diseases as the diffusion in diseased brains is disturbed or altered compared to healthy brains. For instance, Zhang et al. [9] proposed a novel Alzheimer’s disease multi-class classification framework based on multi-modal neuroimaging.

In mathematical morphology, the step from univariate to multivariate pixel values involves the challenging task of ordering vectors to define the notions of minimum and maximum of a set of vectors. An ordering structure for vectors can be derived from the concept of depth [10]. For a set of vectors, depth functions assign to each vector a value that measures its “centrality” within the set. The “centre” is the vector maximising the depth function. Thus, a centre-outward ordering based on the depth values of each vector yields a sound definition of minimum and maximum.

Velasco-Forero et al. [11] defined multivariate mathematical morphological operators via random projection depth. They illustrated their approach on colour and hyperspectral images. Concepts of depth for directional data, i.e. unit vectors in \(\mathbb {R}^d\), have been studied, for instance, by Liu et al. [12] and Pandolfo et al. [13]. Directional statistics [14, 15], which deals with the analysis of directional data, provides the theoretical framework for such studies. Ley et al. [16] introduced quantiles for directional data and the angular Mahalanobis depth.

Based on an ordering on the unit sphere \({\mathcal {S}^{d-1}}\), mathematical morphology was extended to directional images by several authors. Roerdink [17] introduced mathematical morphology on the sphere via generalised Minkowski operations. His motivation is morphological processing of earth pictures under consideration of the curvature. Morphological operators for angle-valued images were introduced by Peters [18] and Hanbury et al. [19]. However, a generalisation of their results to the unit sphere \({\mathcal {S}^{d-1}}\) with \(d> 2\) is not straightforward. Frontera-Pons and Angulo [20] defined morphological operators via a local partial ordering. They used the Fréchet–Karcher barycenter as local origin \(\mu \) on the sphere. The maximum and minimum of a set of vectors are then found via projecting the vectors into the tangent space at \(\mu \). A drawback of their approach is that the minimum and maximum derived this way are not necessarily elements of the given vector set, which seems unnatural. Angulo [21] extended the concept of multi-scale morphological operators, i.e. a scale-dependent morphology and filters, from \(\mathbb {R}\)-valued to \(\mathbb {R}^d\)-valued images. He defined morphological scale space operators on metric Maslov measurable spaces for images supported on point clouds. An application of his operators to RGB-valued point clouds, e.g. for object extraction, showed the usefulness of applying the theory of multi-scale morphological operators to \(\mathbb {R}^d\)-valued images.

In this work, we introduce morphological operators on directional images. Section 2 contains basics about mathematical morphology for grey-scale images and statistical depth functions on \({\mathcal {S}^{d-1}}\). In Sect. 3, we introduce the concept of mathematical morphology on directional data using the directional projection depth. These morphological operators are extended to multi-scale operators in Sect. 4. We illustrate and interpret our findings on synthetic and real-world \(\mathcal {S}^2\)-valued images in Sect. 5. A short conclusion is given in Sect. 6.

2 Basics and Notation

2.1 Images

We denote by \({\mathcal {S}^{d-1}}= \{x\in \mathbb {R}^d: \Vert x \Vert _2 = 1\}\) the \((d-1)\)-dimensional unit sphere and by \(\Vert x \Vert _2 = \sqrt{x^T x}\) the Euclidean norm. An image is a mapping \(I:E \rightarrow S\) that is defined on a q-dimensional domain E with \(E\subset \mathbb {R}^q\) or \(E\subset \mathbb {Z}^q\). The set S is the set of pixel values. We call I a binary image if \(S=\{0,1\}\), a grey-scale image if \(S\subset \mathbb {R}\), a vector-valued image if \(S\subset \mathbb {R}^d\) and a directional (\({\mathcal {S}^{d-1}}\)-valued) image if \(S \subset {\mathcal {S}^{d-1}}\). Usually, \(q,d\in \{2,3\}\). In this paper, we consider examples of q-dimensional \({\mathcal {S}^{d-1}}\)-valued images with \(q=2,3\) and \(d=3\). We denote by

  • \(\max I\) the maximal possible pixel value of an image I. For instance, \(\max I=1\) if I is a binary image or \(\max I=255\) if I is an 8-bit grey-scale image.

  • \(\overline{\mathbb {R}}= \mathbb {R}\cup \{-\infty , \infty \}\) the extended real numbers.

  • SO(d) the rotation group on \(\mathbb {R}^d\).

  • \(x\times y\) the cross product between two vectors \(x,y\in \mathbb {R}^d\).

  • \(\mathcal{G}\mathcal{C}_{\mu ,x}= \{y\in \mathbb {R}^d: (x\times \mu )^Ty=0\} \cap {\mathcal {S}^{d-1}}\) the great circle containing the vectors \(\mu \) and \(x\in {\mathcal {S}^{d-1}}\). For \(d=3\), the great circle \(\mathcal{G}\mathcal{C}_{\mu ,x}\) corresponds to a closed curve on the surface of \(\mathcal {S}^2\) created by the intersection of \(\mathcal {S}^2\) and a 2-dimensional hyperplane H passing through the origin \(0_3\) [22].

  • \(\check{b}\) the reflection of a function \(b:\mathbb {R}^q\rightarrow \overline{\mathbb {R}}\), i.e. \(\check{b}(i) = b(-i)\) for all \(i\in \mathbb {R}^q\), and by \(\check{B}\) the reflection of a set \(B\subset \mathbb {R}^q\), i.e. \(\check{B} = -B\). A function b is called symmetric if \(b=\check{b}\) and a set B is symmetric if \(B=\check{B}\).

2.2 Mathematical Morphology

We first define morphological operators on grey-scale images I. See [1, 2, 4] for a detailed introduction. The two fundamental operations of mathematical morphology are erosion \(\tilde{\varepsilon }_b\) and dilation \(\tilde{\delta }_b\). They depend on a structuring function \(b:\mathbb {R}^q\rightarrow \overline{\mathbb {R}}\) and are defined as

$$\begin{aligned} \tilde{\varepsilon }_b(I)(i)&= \inf _{j\in E} \{I(j)-b(j-i)\}, \quad i\in E, \end{aligned}$$
(1)
$$\begin{aligned} \tilde{\delta }_b(I)(i)&= \sup _{j\in E} \{I(j)+b(j-i)\}, \quad i\in E. \end{aligned}$$
(2)
Fig. 1
figure 1

Binary image \(I =\{0 \text { (= black) },1 \text { (= white)}\}\) (a), flat erosion \(\tilde{\delta }_B\) (b), flat dilation \(\tilde{\varepsilon }_B\) (c), flat opening \(\tilde{\gamma }_B\) (d), and flat closing \(\tilde{\varphi }_B\) (e) with B a square (illustrated in blue in (b, c) as a flat structuring element. To illustrate the changes of the rectangles under the morphological operations, the border of the original rectangles is displayed in grey in b–d (Color figure online)

Full size image
Fig. 2
figure 2

Grey-scale image \(I_1\) (a), flat erosion \(\tilde{\delta }_B\) (b), flat dilation \(\tilde{\varepsilon }_B\) (c), flat opening \(\tilde{\gamma }_B\) (d), \(\tilde{\gamma }_{B_2}\) (e), flat closing \(\tilde{\varphi }_{B_2}\) (f) with B a square of edge length 3 and \(B_2\) a square of edge length 5

Full size image

In the discrete case, supremum and infimum can be replaced by maximum and minimum. We will not discuss edge effects here (see [23]). Throughout the paper, we assume that b is symmetric.

A well-known and widely used example for b is the flat structuring function. Given a structuring element \(B\subset \mathbb {R}^q\), it is defined as

$$\begin{aligned} b(i)={\left\{ \begin{array}{ll} 0 & \quad i\in B \\ -\infty & \quad i\not \in B. \\ \end{array}\right. } \end{aligned}$$

We assume that B is centred at the origin and symmetric. Morphological operators with flat structuring functions are called flat operators. Symbols for flat operators will be written with index B rather than b, e.g., a flat erosion with structuring element B is denoted by \(\tilde{\varepsilon }_B\). On the binary image I given in Fig. 1a, we illustrate the flat erosion \(\tilde{\varepsilon }_B\) and dilation \(\tilde{\delta }_B\), where B is a square, in Fig. 1b, c, respectively. Illustrations for a grey-scale image are shown in Fig. 2.

Non-flat or volumic structuring functions assign non-constant weights to the pixel values [2]. For instance,

$$\begin{aligned} b(i)&= -\left( \frac{\Vert i \Vert _2}{2}\right) ^2 \end{aligned}$$

is a non-flat structuring function. Note that the grey-scale ranges of images dilated or eroded by non-flat structuring functions are not bounded [4]. For instance, an erosion with a non-flat structuring function of an image with non-negative pixel values can result in negative pixel values. Flat structuring functions do not suffer from that problem. The output of a flat erosion or dilation is bounded by the grey-scale range of the input image.

The composition of erosion and dilation yields the morphological operators opening \(\tilde{\gamma }_b\) and closing \(\tilde{\varphi }_b\)

$$\begin{aligned} \tilde{\gamma }_b(I)(i)&= \tilde{\delta }_{\check{b}}(\tilde{\varepsilon }_b(I)) (i), \\ \tilde{\varphi }_b(I)(i)&= \tilde{\varepsilon }_{\check{b}}(\tilde{\delta }_b(I)) (i). \end{aligned}$$

For the binary image I given in Fig. 1a, we illustrate the flat opening \(\tilde{\gamma }_B\) and closing \(\tilde{\varphi }_B\) with B a square in Fig. 1b, c, respectively. Illustrations for a grey-scale image are shown in Fig. 2. In applications, the opening is used to remove bright noise and the closing to remove dark noise.

Some properties of the morphological operators are summarised in the following, see also [1, 2, 23]. Let \(i\in E\), \(I,I'\) grey-scale images, and \(\{I_l\}_{l\in \mathbb {N}}\) a family of grey-scale images. We write \(I \le I'\) if \(I(i) \le I'(i)\) for all \(i\in E\). Furthermore, \(\bigvee \) denotes the point-wise maximum operator and \(\bigwedge \) the point-wise minimum operator.

  1. 1.

    All the morphological operators are nonlinear, i.e. in general

    $$\begin{aligned} \tilde{\Psi }(aI+bI') \not = a\tilde{\Psi }(I)+b\tilde{\Psi }(I'), \end{aligned}$$

    where \(a,b\in \mathbb {R}\) and \(\tilde{\Psi }=\tilde{\delta }_b, \tilde{\varepsilon }_b, \tilde{\gamma }_b\), or \(\tilde{\varphi }_b\).

  2. 2.

    Dilation and erosion are dual w.r.t. complementation \(\tilde{{\,\textrm{C}\,}}\) with \(\tilde{{\,\textrm{C}\,}}I(i)=\max I -I(i)\), i.e.

    $$\begin{aligned} \tilde{\delta }_b (I) = \tilde{{\,\textrm{C}\,}}\tilde{\varepsilon }_b( \tilde{{\,\textrm{C}\,}}I) \end{aligned}$$

    with \(b(i) = b(-i)\) due to symmetry.

  3. 3.

    Opening and closing are dual w.r.t. complementation \(\tilde{{\,\textrm{C}\,}}\), i.e.

    $$\begin{aligned} \tilde{\gamma }_b (I) = \tilde{{\,\textrm{C}\,}}\tilde{\varphi }_b( \tilde{{\,\textrm{C}\,}}I) \end{aligned}$$
  4. 4.

    Opening and closing are idempotent, i.e.

    $$\begin{aligned} \tilde{\gamma }_b (\tilde{\gamma }_b (I))&= \tilde{\gamma }_b (I), \\ \tilde{\varphi }_b (\tilde{\varphi }_b (I))&= \tilde{\varphi }_b (I). \end{aligned}$$
  5. 5.

    The following distribution laws hold

    $$\begin{aligned} \tilde{\delta }_b \left( \bigvee _l I_l\right)&= \bigvee _l \tilde{\delta }_b (I_l), \\ \tilde{\varepsilon }_b \left( \bigwedge _l I_l\right)&= \bigwedge _l \tilde{\varepsilon }_b (I_l). \end{aligned}$$
  6. 6.

    Dilation is associative and erosion fulfils the chain rule, i.e.

    $$\begin{aligned} \tilde{\delta }_b (\tilde{\delta }_{b'} (I))&= \tilde{\delta }_{\tilde{\delta }_{\check{b}}(b')} (I) , \\ \tilde{\varepsilon }_b (\tilde{\varepsilon }_{b'} (I))&= \tilde{\varepsilon }_{\tilde{\delta }_{\check{b}}(b')} (I) . \end{aligned}$$
  7. 7.

    All morphological operators are increasing, i.e.

    $$\begin{aligned} I \le I' \Rightarrow \tilde{\Psi }(I) \le \tilde{\Psi }(I') , \end{aligned}$$

    where \(\tilde{\Psi }=\tilde{\delta }_b, \tilde{\varepsilon }_b, \tilde{\gamma }_b\), or \(\tilde{\varphi }_b\).

  8. 8.

    If b is defined at the origin and \(b(0)\ge 0\), dilation is extensive and erosion is anti-extensive, i.e.

    $$\begin{aligned} \tilde{\varepsilon }_b (I) \le I \le \tilde{\delta }_b (I). \end{aligned}$$
  9. 9.

    Closing is extensive and opening is anti-extensive, i.e.

    $$\begin{aligned} \tilde{\gamma }_b (I) \le I \le \tilde{\varphi }_b (I). \end{aligned}$$

The latter two properties can be combined to give the ordering relation

$$\begin{aligned} \tilde{\varepsilon }_b(I)&\le \tilde{\gamma }_b(I) \le I \le \tilde{\varphi }_b(I) \le \tilde{\delta }_b(I). \end{aligned}$$
(3)

2.3 Morphological Scale Space

In the context of mathematical morphology for grey-scale images, erosion and dilation can be written as scale-space operators [2, 24]

$$\begin{aligned} \tilde{\varepsilon }_{b_{t}}(I)(i)&= \inf _{j\in E} \{I(j)-{b_{t}}(j-i)\} \end{aligned}$$
(4)
$$\begin{aligned} \tilde{\delta }_{b_{t}}(I)(i)&= \sup _{j\in E} \{I(j)+{b_{t}}(j-i)\}, \end{aligned}$$
(5)

with a scaled structuring function

$$\begin{aligned} {b_{t}}:\mathbb {R}^q\rightarrow \overline{\mathbb {R}}, \quad t\ge 0 \end{aligned}$$

with scale parameter t. Following [23, p.15], we assume \({b_{t}}\) to have the following properties:

  1. 1.

    \(\{{b_{t}}\}_t\) is a one-parametric family of convex, continuous, symmetric functions and fulfils the semi-group property, i.e.

    $$\begin{aligned} {b_{t}}\dot{+}{b_{s}}&={b_{t+s}}\qquad t,s \ge 0, \end{aligned}$$
    (6)

    where \(\dot{+}\), \(+\) are group operations [24]. For instance, when choosing \(b_t\) to be a flat structuring function with structuring element tB, the group operations are the Minkowski addition and the standard addition in \(\mathbb {R}\):

    $$\begin{aligned} tB\oplus sB&=(t+s)B \quad t,s \ge 0. \end{aligned}$$
  2. 2.

    \({b_{t}}\) is non-positive and monotonically decreasing with a global maximum at the origin of value zero, i.e.

    $$\begin{aligned} {b_{t}}(i)&\le 0 \qquad ~~\forall ~t \ge 0, i\in \mathbb {R}^q, \end{aligned}$$
    (7)
    $$\begin{aligned} {b_{t}}(i)&\ge {b_{t}}(j) \quad \Vert i \Vert _2 < \Vert j \Vert _2, \end{aligned}$$
    (8)
    $$\begin{aligned} {b_{t}}(0)&=0. \end{aligned}$$
    (9)

An example of a non-flat scaled structuring function is the Poweroid structuring function [23, p.14, Def. 5]

$$\begin{aligned} {b_{t}}(i)&= - t\left( \frac{\Vert i \Vert _2}{t}\right) ^a \quad a\ge 0, t>0. \end{aligned}$$
(10)

2.4 Directional Data

Directional data appear across a variety of fields, including astronomy, biology, and engineering. Instances range from the arrival direction of cosmic rays [22] to the flight direction of migratory birds [14], or to the alignment of steel fibres within fibre-reinforced concrete [25].

Directional data are given as a set of vectors on the unit sphere \({\mathcal {S}^{d-1}}\) for some \(d\ge 2\). The vector \(x=(x_1,\dots ,x_d)^T\in {\mathcal {S}^{d-1}}\) can be represented as a point on the surface of \({\mathcal {S}^{d-1}}\). For x in Cartesian coordinates, a representation in angular coordinates for arbitrary d is given in [26].

For \(d=3\), we use the usual notation for spherical coordinates

$$\begin{aligned} x=x(\phi ,\theta )&=\begin{pmatrix} \sin {(\theta )}\cos {(\phi )}\\ \sin {(\theta )}\sin {(\phi )}\\ \cos {(\theta )} \end{pmatrix} \end{aligned}$$
(11)

with co-latitude \(\theta \in [0,\pi ]\) and longitude \(\phi \in [0,2\pi )\). A random vector \(X\in \mathcal {S}^2\) can be represented in spherical coordinates by random variables \(\Theta \) and \(\Phi \) with realisations \(\theta \) and \(\phi \).

2.5 Statistical Depth Functions

Depth functions are applied in multidimensional non-parametric robust data analysis and establish an ordering relation between vectors. See [10] for a detailed introduction.

The literature gives several depth functions, for instance, halfspace depth [27], simplicial depth [28], projection depth [29], spatial depth [30], or the Mahalanobis depth [10]. For defining a depth for directional data, we initially restrict attention to the class \(\mathcal {F}_\mu \) of distributions on \({\mathcal {S}^{d-1}}\) with a bounded density that admits a unique modal direction \(\mu \). Examples of distributions in \(\mathcal {F}_\mu \) are the von Mises–Fisher distribution, the Kent distribution, or the Bingham distribution [22]. We further assume that \(\mu \) coincides with the Fisher spherical median [31], that is,

$$\begin{aligned} \mu = \arg \min _{\gamma \in {\mathcal {S}^{d-1}}}E(\arccos (X^T\gamma )). \end{aligned}$$

Ley et al. [16] adapted the definition of a depth function given in [10, Definition 2.1] for directional data as follows.

Definition 1

A statistical depth function on \({\mathcal {S}^{d-1}}\) is a bounded, non-negative mapping \(\mathcal {D}_\cdot (\cdot ): \mathcal {F}_\mu \times {\mathcal {S}^{d-1}}\rightarrow \mathbb {R}\) satisfying

  1. 1.

    \(\mathcal {D}_{F_{AX}}(AX)=\mathcal {D}_{F_{X}}(X)\) holds for any random vector \(X\in {\mathcal {S}^{d-1}}\) and any rotation matrix \(A\in \mathbb {R}^{d\times d}\).

  2. 2.

    \(\mathcal {D}_F(\mu )=\sup _{y\in {\mathcal {S}^{d-1}}}\mathcal {D}_F(y)\) holds for any \(F\in \mathcal {F}_\mu \).

  3. 3.

    For any \(F\in \mathcal {F}_\mu \), \(\mathcal {D}_F(x)\le \mathcal {D}_F(c(\alpha ))\) for the unique geodesic \(c:[0,1]\rightarrow {\mathcal {S}^{d-1}}\) with \(c(0)=\mu \), \(c(1)=x\) and \(\alpha \in [0,1]\).

  4. 4.

    \(\mathcal {D}_F(-\mu ) = 0\) for each \(F\in \mathcal {F}_\mu \) where \(-\mu \) is the antipodal point of the centre \(\mu \).

Property 1 implies that the depth of a vector should be independent of the underlying coordinate system. Property 2 means that for any F with unique centre \(\mu \), the depth function is maximal at \(\mu \). Property 3 indicates that for any point X moving away from the centre \(\mu \) along any geodesic starting at \(\mu \), the depth of X decreases monotonically. Property 4 implies that the depth value of X is zero at \(-\mu \) since \(-\mu \) is at maximal distance from \(\mu \) on \({\mathcal {S}^{d-1}}\).

Depth functions for directional data are, for instance, directional distance-based depths [13] or the angular Mahalanobis depth [16]. The ordering property (3) for morphological operators can be violated if \(\mu \) is not fixed (see Appendix A). We use a re-scaled version of the cosine distance depth from [13] as defined in the following section and select \(\mu \) as a global parameter as described in Sect. 3.

2.6 The Angular Projection Depth

Let \(X\in {\mathcal {S}^{d-1}}\) be a random vector with \(X\sim F\in \mathcal {F}_\mu \). We define a depth for directional data by

$$\begin{aligned} D^{proj}_{F}(X) = \frac{1 + X^T\mu }{2}, \quad X\in {\mathcal {S}^{d-1}}. \end{aligned}$$
(12)

\(D^{proj}_F(X)\) provides a centre-outward ordering with \(D^{proj}_F(\mu ) = 1\), \(D^{proj}_F(-\mu ) = 0\) and is decreasing on a geodesic from \(\mu \) to \(-\mu \). Let \(A\in \mathbb {R}^{d\times d}\) be any rotation matrix. Then the distribution of the transformed vector AX has centre \(A\mu \). It follows that

$$\begin{aligned} D^{proj}_{F_{AX}}(AX)&=\frac{1 + (AX)^T A\mu }{2} =\frac{1 + X^TA^TA\mu }{2} \\&=\frac{1 + X^T\mu }{2} =D^{proj}_{F_{X}}(X). \end{aligned}$$

Thus, the properties of a directional depth function given in Definition 1 are fulfilled.

The empirical angular projection depth reads

$$\begin{aligned} D^{proj}(x)&= \frac{1 + x^T\hat{\mu }}{2}, \end{aligned}$$

where \(\hat{\mu }\) is the empirical Fisher spherical median [31]

$$\begin{aligned} \hat{\mu }= \arg \min _{\gamma \in \mathcal {S}^2}\sum _{i=1}^N\arccos (X_i^T\gamma ) \end{aligned}$$
(13)

obtained from an i.i.d. sample \(X_1,\dots ,X_n\) with distribution F.

3 Mathematical Morphology on Directional Images Using Directional Projection Depth

The definition of the standard operators of mathematical morphology requires an ordering relation between pixel values, see (1) and (2). Multivariate extensions of ordering and mathematical morphology were subject to several studies [7, 32,33,34,35]. In particular, orderings can be derived from depth functions.

Fig. 3
figure 3

Vector ordering by directional projection depth

Full size image

For \({\mathcal {S}^{d-1}}\) valued images, we use the projection depth \(D^{proj}\) to define an ordering relation between unit vectors. See Fig. 3 for an illustration. This allows for a definition of erosion, dilation, and the morphological operators derived from them as described below. For simplicity, we concentrate on flat morphological operators with structuring element B.

In a given application, the assumption \(I(i)\sim F \in \mathcal {F}_\mu \) may not be fulfilled. Nevertheless, a suitable central direction \(\mu \) may be derived from the experimental setup. Hence, \(\mu \) will be interpreted as a parameter in the following. We will write \(D_\mu ^{proj}\) for the projection depth using \(\mu \in {\mathcal {S}^{d-1}}\) as the centre. The parameter \(\mu \) can be selected globally for the whole image or in a locally adaptive manner for parts of the image. The latter is not considered here.

The erosion \(\varepsilon _B\) of a directional image \(I\) at pixel position \(i\in E\) is (implicitly) defined by

$$\begin{aligned} D_\mu ^{proj}(\varepsilon _B(I)) (i)&= \tilde{\varepsilon }_B(D_\mu ^{proj}(I))(i), \end{aligned}$$

where \(\tilde{\varepsilon }\) is the flat erosion of a grey-scale image. Analogously, the dilation \(\delta _B\) of a directional image \(I\) at pixel position \(i\in E\) is (implicitly) defined by

$$\begin{aligned} D_\mu ^{proj}(\delta _B(I))(i)&= \tilde{\delta }_B(D_\mu ^{proj}(I))(i), \end{aligned}$$

where \(\tilde{\delta }\) is the flat dilation of a grey-scale image. Hence, a dilation will select the most central direction covered by the structuring element while an erosion selects the most outlying direction.

The explicit definition reads

$$\begin{aligned} \varepsilon _B(I)(i)&= {D_\mu ^{proj}}^{-1}\left( \tilde{\varepsilon }_B(D_\mu ^{proj}(I))\right) (i),\\ \delta _B(I)(i)&= {D_\mu ^{proj}}^{-1}\left( \tilde{\delta }_B(D_\mu ^{proj}(I))\right) (i), \end{aligned}$$

where \({D_\mu ^{proj}}^{-1}\) refers to the preimage under \({D_\mu ^{proj}}\).

As \(D_\mu ^{proj}\) is not injective, the preimage may consist of more than one element. Consider vectors \(x_1,x_2\in {\mathcal {S}^{d-1}}\) with \(x_1^T\mu =x_2^T\mu \). Then, \(D_\mu ^{proj}(x_1)=D_\mu ^{proj}(x_2)\) but \(x_1=x_2\) is not necessarily the case. To resolve this issue, we use a lexicographic order for vectors of the same depth value which yields a total ordering [11]. Rotate the vectors such that \(\mu =(0,0,1)^T\). Among the (pairwise different) vectors \(\{x_i\}_{i=1,\dots ,n}\) with \(D_\mu ^{proj}(x_i)=D_\mu ^{proj}(x_j)\), we choose the vector \(x^*\) with the smallest longitude angle, i.e. \(x^*=\arg \min \limits _{x_i}\phi _i\) with \((\theta _i,\phi _i)\) spherical coordinates of \(x_i\), \(i=1,\dots ,n\). Then rotate \(x^*\) back.

Opening and closing are defined by

$$\begin{aligned} \gamma _B(I)(i)&= \delta _{\check{B}}(\varepsilon _B(I))(i),\\ \varphi _B(I)(i)&= \varepsilon _{\check{B}}(\delta _B(I))(i). \end{aligned}$$

The theory of h-orderings [7] immediately implies that the operators defined above fulfil the properties of morphological operators stated in Sect. 2.2, see [32]. In this theory, h is a map from \({\mathcal {S}^{d-1}}\) into a space with a partial order. Here we choose \(h=D^{proj}\).

Of course, further morphological filters have their depth analogy. For instance, the scalar difference between the depth of dilation and erosion defines the morphological gradient

$$\begin{aligned} {{\,\textrm{g}\,}}_{D_\mu ^{proj},B}(I)(i)&\quad = D_\mu ^{proj}\left( \delta _B(I)(i)\right) - D_\mu ^{proj}\left( \varepsilon _B(I)(i)\right) . \end{aligned}$$

The morphological Laplacian is defined by

$$\begin{aligned} \Delta _{D_\mu ^{proj},B}(I)(i)&= \Delta _\delta (i) - \Delta _\varepsilon (i) \end{aligned}$$

with \(\Delta _\delta (i) = D_\mu ^{proj}\left( \delta _B(I)(i)\right) - D_\mu ^{proj}(I)(i)\) and \(\Delta _\varepsilon (i) = D_\mu ^{proj}(I)(i) - D_\mu ^{proj}\left( \varepsilon _B(I)(i)\right) \). The shock filter (also known as Toggle Mapping) is defined in [3]

$$\begin{aligned} {{\,\textrm{sf}\,}}_{D_\mu ^{proj},B}(I)(i) ={\left\{ \begin{array}{ll} \varepsilon _B(I)(i) & \Delta _{D_\mu ^{proj},B}(I)(i) < 0, \\ \delta _B(I)(i) & \Delta _{D_\mu ^{proj},B}(I)(i) > 0, \\ I(i) & \text { otherwise.} \\ \end{array}\right. } \end{aligned}$$

Shock filtering is used to enhance edges. For grey-scale images, the idea is to dilate near local maxima and erode near local minima. This way, the contrast between regions of pixel values with high and low depth is enhanced. If we apply this to our depth approach, we can contrast between regions with different orientations (depending on \(\mu \)).

4 Pseudo-Morphological Multi-scale Operators for Directional Images

In the next step, we want to extend morphological operators to a multi-scale setting, leading to a morphological scale space (see [21]). A very simple approach to scaling is scaling of a flat structuring element B by a factor t. Here, we will consider an alternative approach which yields a structuring function which has a meaningful interpretation, fulfils the semi-group property (at least partially), and results in bounded non-flat operators. We call them pseudo-morphological operators since their behaviour reminds of erosion and dilation as demonstrated in Sect. 5 while some properties of morphological operators are lost (see Remark 1).

4.1 Structuring Function for Directional Images

In general, structuring functions \({b_{t}}:E\rightarrow \overline{\mathbb {R}}\) are unbounded. Thus, multi-scale dilation and erosion inherit this unboundedness if we use \(D_\mu ^{proj}(I(j))+{b_{t}}(j-i)\) and \(D_\mu ^{proj}(I(j))-{b_{t}}(j-i)\), respectively, \(i,j\in E\). However, we cannot associate a vector with a depth value outside [0, 1] since no unit vector could be found as preimage. In our approach, we therefore use structuring functions \({b_{t}}:E\rightarrow SO(d)\). We will define the operation \(I(j)+{b_{t}}(j-i)\) and \(I(j)-{b_{t}}(j-i)\) to be a vector rotation of \(I(j)\) where the rotation angle depends on the scale t. Therefore, the depths of \(I(j)+{b_{t}}(j-i)\) and \(I(j)-{b_{t}}(j-i)\) are both bounded.

The main idea is to rotate a pixel value \(x=I(i)\), \(i\in E\), about \(x\times \mu \) towards \(\mu \) or away from \(\mu \) on a great circle \(\mathcal{G}\mathcal{C}_{\mu ,x}\). This increases or decreases the depth value of the rotated x which gives our approach also a meaningful interpretation. Note that \(I(j)+{b_{t}}(j-i)\) and \(I(j)-{b_{t}}(j-i)\), \(i\in E\), result again in directional images which seem natural.

Fig. 4
figure 4

Illustration of the operations \(x - {b_{t}}\) given in (18) and \(x + {b_{t}}\) given in (19) with \(d=3\) and \(x=x(\phi ,\theta )\) in spherical coordinates (11). \(x + {b_{t}}\) increases the geodesic distance between x and \(\mu \) by \(\alpha _+\). \(x - {b_{t}}\) decreases the geodesic distance between x and \(\mu \) by \(\alpha _-\). x moves on the great circle \(\mathcal{G}\mathcal{C}_{\mu ,x}\)

Full size image

For a formal definition, let \(\mu \in {\mathcal {S}^{d-1}}\), \(i\in E\), \(x=I(i)\in {\mathcal {S}^{d-1}}\) with \(x^T\mu = \cos (\theta )\) and \(R_{\mu ,x}(\alpha )\in SO(d)\) be a rotation matrix which rotates x about the cross product \(\mu \times x\) and \(R_{\mu ,x}(\theta )x=\mu \). We define a structuring function \({b_{t}}\) for directional images by a mapping

$$\begin{aligned}&{b_{t}}: E \rightarrow SO(d) \nonumber \\&i \mapsto {b_{t}}(i) := R_{\mu ,x}({\alpha _t(i)}), \end{aligned}$$
(14)

with rotation angle \({\alpha _t(i)}\), \(t\ge 0\).

Since \(R_{\mu ,x}\) depends on x, the structuring function is locally adaptive. We choose

$$\begin{aligned} {\alpha _t(i)}&= \min \left( \frac{\Vert i \Vert _2^2}{t}, \pi \right) , \end{aligned}$$
(15)

where the truncation at \(\pi \) is discussed in more detail in Remark 1. Note the similarity of \(\alpha _t\) with quadratic scaling as discussed in [36].

Let \(f\Box g\) be the infimal convolution of two functions \(f,g: \mathbb {R}^q\rightarrow \mathbb {R}\) [37], i.e.

$$\begin{aligned} (f \Box g) (i) = \inf _{j\in \mathbb {R}^q} \left\{ f(i-j) + g(j) \right\} . \end{aligned}$$
(16)

Define

$$\begin{aligned} ({b_{t}}\dot{+} {b_{s}})(i) := R_{\mu ,x}( (\alpha _t \Box \alpha _s)(i)). \end{aligned}$$
(17)

The structuring element \({b_{t}}\) fulfils the semi-group property (6), i.e. \({b_{t}}(i) \dot{+} {b_{s}}(i) = {b_{t+s}}(i)\), for \(\Vert i \Vert _2^2 / (t+s) < \pi \) as shown in Appendix B.

We illustrate the construction for \(d=3\) and \(x=x(\phi ,\theta )\) in spherical coordinates as given in (11): The rotation matrix \(R_{\mu ,x}\) rotates x on the great circle \(\mathcal{G}\mathcal{C}_{\mu ,x}\). We use \(R_{\mu ,x}\) to move a pixel value \(x\in \mathcal {S}^2\) towards \(\mu \) or away from \(\mu \) on \(\mathcal{G}\mathcal{C}_{\mu ,x}\). Thus, the rotation matrix \(R_{\mu ,x}\) causes a change of the angle between x and \(\mu \), here the co-latitudinal angle \(\theta \). The longitude angle \(\phi \) is not changed. Therefore, we define a rotation of the vector x depending on the structuring function \({b_{t}}\) to this change of \(\theta \) as follows (see Fig. 4):

$$\begin{aligned} x(\phi ,\theta ) - {b_{t}}(i)&:= x(\phi , \theta -\alpha _-) \end{aligned}$$
(18)

with

$$\begin{aligned} \alpha _-&={\left\{ \begin{array}{ll} {\alpha _t(i)},& \quad \text {for } {\alpha _t(i)}\in [0,\theta ]\\ \theta , & \quad \text {for } {\alpha _t(i)}> \theta . \end{array}\right. } \nonumber \\ x(\phi ,\theta ) + {b_{t}}(i)&:= x(\phi , \theta +\alpha _+) \end{aligned}$$
(19)

with

$$\begin{aligned} \alpha _+&={\left\{ \begin{array}{ll} {\alpha _t(i)}, & \quad \text {for } {\alpha _t(i)}\in [0,\pi -\theta ]\\ \pi - \theta , & \quad \text {for } {\alpha _t(i)}> \pi -\theta . \end{array}\right. } \end{aligned}$$

\(x - {b_{t}}\) has a smaller geodesic distance to \(\mu \) than x and \(x + {b_{t}}\) has a larger geodesic distance to \(\mu \) than x.

Remark 1

(Truncation of rotation) The restriction of \(\alpha _-\) to \(\theta \) and of \(\alpha _+\) to \(\pi -\theta \) is necessary to achieve an analogous ordering relation as in (3). Otherwise, x can be rotated beyond \(\mu \) which would decrease \(D_\mu ^{proj}(x - {b_{t}}(i))\). Analogously, x can be rotated beyond \(-\mu \) which would increase \(D_\mu ^{proj}(x + {b_{t}}(i))\). To avoid this, we fix the vector at \(\mu \) or \(-\mu \), respectively, such that no rotation beyond is possible.

With truncation, the semi-group property is no longer valid. However, truncation ensures that

$$\begin{aligned} D_\mu ^{proj}(x + {b_{t}})&\le D_\mu ^{proj}(x) \le D_\mu ^{proj}(x - {b_{t}}). \end{aligned}$$
(20)

We obtain (in some sense) analogue properties to (7)–(9):

  1. 1.

    \({b_{t}}\) is non-positive w.r.t. \(D_\mu ^{proj}\) in the sense that for rotations of x away from \(\mu \) its depth decreases, i.e.

    $$\begin{aligned} D_\mu ^{proj}(x + {b_{t}}(i)) - D_\mu ^{proj}(x)&\le 0. \end{aligned}$$
  2. 2.

    \({b_{t}}\) is monotonically decreasing w.r.t. \(D_\mu ^{proj}\) in the sense that for increasing distance between pixel positions the depth decreases, i.e.

    $$\begin{aligned} D_\mu ^{proj}(x + {b_{t}}(i))&\ge D_\mu ^{proj}(x + {b_{t}}(j)) \end{aligned}$$

    for \(\Vert i \Vert _2 < \Vert j \Vert _2.\)

  3. 3.

    \({b_{t}}\) has a global maximum at the origin \(o\in \mathbb {R}^q\) w.r.t. \(D_\mu ^{proj}\) with

    $$\begin{aligned} D_\mu ^{proj}(x + {b_{t}}(o)) =D_\mu ^{proj}(x) \end{aligned}$$

4.2 Pseudo-Morphological Multi-scale Operators for Directional Images

Let \({b_{t}}\) be as in Eq. (14). The multi-scale erosion \(\varepsilon _{b_{t}}\) of a directional image \(I\) at pixel position \(i\in E\) is (implicitly) defined by

$$\begin{aligned} D_\mu ^{proj}(\varepsilon _{b_{t}}(I)) (i)&= \inf _{j\in E} \{ D_\mu ^{proj}(I(j) -{b_{t}}(j-i)) \}, \end{aligned}$$
(21)

where \(I(j) -{b_{t}}(j-i)\) is defined in Eq. (18). Analogously, the multi-scale dilation \(\delta _{b_{t}}\) of a directional image \(I\) at pixel position \(i\in E\) is (implicitly) defined by

$$\begin{aligned} D_\mu ^{proj}(\delta _{b_{t}}(I)) (i)&= \sup _{j\in E} \{ D_\mu ^{proj}(I(j) +{b_{t}}(j-i)) \}, \end{aligned}$$
(22)

where \(I(j) +{b_{t}}(j-i)\) is defined in Eq. (19). The interpretation of the multi-scale operators is as follows: We want to rotate the vector \(I(j)\) depending on \(||i-j||_2\) and the scale t towards \(\mu \) (\(D_\mu ^{proj}\) increases) or towards \(-\mu \) (\(D_\mu ^{proj}\) decreases).

Multi-scale opening, closing, morphological gradient, and shock filter can be defined as analogues to their flat counterparts.

5 Examples

In this section, we investigate the effect of the morphological operators on synthetic \(\mathcal {S}^2\)-valued images and compare them to their standard grey-scale counterparts. As real application examples, we detect misaligned regions in a glass fibre reinforced composite and enhance changes in the displacement field of a compressed glass foam.

5.1 Flat Morphological Operators

To investigate the newly defined flat morphological operators, we generate a 2-dimensional \(\mathcal {S}^2\)-valued image \({I}\) mimicking direction vectors obtained from two fibres on a homogeneous background. That is, we assume that the image contains two objects consisting of vectors with a small angular deviation from \(\mu \), see Fig. 5. In the following, we will interpret those as foreground. The background is formed by vectors directed along a plane perpendicular to \(\mu \). Flat dilation and flat erosion (Fig. 6) show similar behaviour as their standard grey-scale counterparts. Dilation expands objects directed along \(\mu \). That is, the directions of background pixels that are close to the edge of these objects are rotated towards \(\mu \). An erosion shrinks objects, i.e. object vectors at the edge are rotated so that they are assigned to the image background.

In a similar manner, the flat opening removes foreground objects that are smaller than the structuring element. Figure 6 reveals that object vectors in an object of size smaller than the structuring element B are rotated such that we assign them to the image background. The flat closing removes small holes in the foreground. Vectors with large angular deviation from \(\mu \) within a background object of size smaller than B are rotated to become part of the image foreground.

The results of shock filtering of \({I}\) are shown in Fig. 7. The edges between the two objects and the background pixels are enhanced.

Fig. 5
figure 5

Original \(\mathcal {S}^2\)-valued image \({I}\). The direction vectors within the fibres are rotationally symmetric around \(\mu \) (large red vector). Their angular deviation from \(\mu \) is uniformly drawn from \((0,\pi /8)\). Background vectors are also rotationally symmetric around \(\mu \). Their angular deviation from \(\mu \) is uniformly drawn from \([3/8\pi ,5/8\pi ]\) (Color figure online)

Full size image
Fig. 6
figure 6

Flat dilation \(\delta _B\) (a, b), flat erosion \(\varepsilon _B\) (c, d), flat opening \(\gamma _B\) (e, f) and flat closing \(\varphi _B\) (g, h) with B a square as a structuring element. The large red vector corresponds to \(\mu \) (Color figure online)

Full size image
Fig. 7
figure 7

Flat shock filter \({{\,\textrm{sf}\,}}_B({I})\) with B a square as structuring element

Full size image

5.2 Pseudo-Morphological Multi-scale Operators

We illustrate the multi-scale dilation and erosion on a directional image \({I}\). \(\delta _{b_{t}}({I})\) and \(\varepsilon _{b_{t}}({I})\) behave like their grey-scale counterparts: Fig. 8 shows that \(\delta _{b_{t}}({I})\) enlarges objects in the foreground and Fig. 9 shows that \(\varepsilon _{b_{t}}({I})\) shrinks objects in the foreground. Of course, the enlargement and shrinkage depend on the continuous scale parameter t.

Fig. 8
figure 8

\(\delta _{b_{t}}({I})\) at scales \(t=0.1,0.5,0.7,0.9,1.1,2\). The large red vector corresponds to \(\mu \) (Color figure online)

Full size image
Fig. 9
figure 9

\(\varepsilon _{b_{t}}({I})\) at scales \(t=0.1,0.5,0.7,0.9,1.1,2\). The large red vector corresponds to \(\mu \) (Color figure online)

Full size image

5.3 Segmentation of Aligned and Misaligned Fibre Regions of Glass Fibre Reinforced Polymers

Glass fibre reinforced polymers (GFRP) are frequently used in lightweight construction of, e.g., cars, planes, and wind power plants. In civil engineering, they have emerged as an alternative to steel reinforcement as they do not corrode in aggressive conditions such as marine environments [38]. GFRP materials are anisotropic and characterised by high tensile strength in the direction of the reinforcing fibres. When manufactured by injection moulding, the fibres follow the injection direction. However, the fibre orientation in a central layer differs depending on the production parameters. Analysing fibre orientation and quantifying misalignment, for instance by directional filtering [39], may help to optimise the production process and to improve GFRP.

We examine a \(\mu \)CT image of GFRP produced and imaged by the Leibniz Institute for Composite Materials (IVW) in Kaiserslautern (Germany) and discussed in [40]. See Fig. 10a for an illustration. Pixel-wise fibre orientations can be determined by using partial second derivatives as described in [41] resulting in an \(\mathcal {S}^2\)-valued image.

Our goal is to segment the region aligned along the injection direction, as desired, and the misaligned fibre region. Figure 10b shows a grey value encoding of the scalar product between the pixelwise fibre directions and \(\mu \). The correctly aligned regions are at the left and right borders while fibres in the central part are misaligned. The higher the value (brighter the image pixel), the more aligned the vectors are along \(\mu \). Simply thresholding this image results in a very noisy segmentation. Therefore, the directional image is smoothed prior to the segmentation by a morphological closing along \(\mu =(0,1,0)\) with B a \(7\times 7\times 7\) cube.

We choose a threshold value of 1.5 to segment areas aligned along \(\mu \), see Fig. 10b.

Fig. 10
figure 10

Sectional \(\mu \)CT image of GFRP produced and imaged by the Leibniz Institute for Composite Materials (IVW) in Kaiserslautern (Germany), the corresponding original and closed directional image

Full size image

5.4 Enhancement of Fault Zones in Compressed Glass Foam

To characterise the structural behaviour of complex materials during loading, mechanical tests and simultaneous \(\mu \)CT imaging (in situ CT) can be used to estimate local displacements of the material on the micro-scale. We analyse a glass foam which consists of very thin struts.

The \(\mu \)CT images were recorded while the glass foam was compressed in the z-direction. The foam is expected to fail very suddenly due to its thin struts. Nogatz et al. [42] computed the displacement field for the transition from strain level 1–3.8% which is shown in Fig. 11a. During compression, a fault zone forms which is now post-processed by directional morphology. To this end, our operators are applied to the directional component of the displacement vectors. That is, the vectors were normalised for the calculation of the morphological operators and were subsequently scaled back to their original length. As we mostly expect compression along the loading direction, we choose \(\mu =(0,0,1)\).

We enhance the fault zone with the morphological gradient, see Fig. 11b. If motion estimation algorithms fail to reconstruct these sharp edges, we can enhance them by erosion, as shown in Fig. 11c. Some other materials, however, are known to show creep before failure. Here, a smoother transition seems more reasonable, which is obtained by a dilation, see Fig. 11d.

Closing (Fig. 11e) and opening (Fig. 11f) remove misalignment in the transition between the two regions.

Fig. 11
figure 11

XZ-slices of displacement field based on [42]. Yellow colour indicates movement along the z-direction, and blue refers to the opposite direction. The computed displacement field reflects the influence on the microstructure during compression (a). The morphological gradient (with B a \(3\times 3\) square) enhances the fault zone (b). The results of other morphological operators are given in (c–f). We chose B as a \(15\times 15\) square (Color figure online)

Full size image

6 Conclusion

We have defined morphological operators for \({\mathcal {S}^{d-1}}\)-valued images. The required ordering for unit vectors is derived from the angular projection depth which enables a sound definition of \({\mathcal {S}^{d-1}}\)-valued morphological operators and filters. Relations of these morphological operators to their grey-scale counterparts are emphasised. Additionally, we segmented misaligned fibre regions of glass fibre reinforced polymers and enhanced the fault region in a compressed glass foam.

Future work could address the removal of disturbances or the highlighting of directional changes in DTI images by morphological operations [8, 9].