Single Scan Granulometry Estimation from an Asymmetric Distance Map

Abstract

Granulometry, by characterizing the distribution of object sizes, is a powerful tool for the analysis of binary images. It may, for example, extract pertinent features in the context of texture classification. The granulometry is similar to a sieving process that filters image details of increasing sizes. Its computation classically relies on morphological openings, i.e. the sequence of an erosion followed by a dilation for each ball radius.

It is well known that a distance map can be described as an “erosion transform” that summarizes the erosions with the balls of all radii. Using a Steiner formula, we show how a vector of parameters measured on an eroded contour can be extrapolated to the measures of the dilation. Instead of completing the openings from the distance map by computing a dilation for each ball size, we can estimate their cardinality from the area, perimeter and Euler-Poincaré characteristic of each erosion. We extract these measures for all radii at once from an asymmetric distance map. The result is a fast streaming algorithm that provides an estimate of the granulometry distribution, in a single scan of the image and with very limited memory footprint.

1 Preliminaries

In this paper, we consider images defined on the discrete square grid $\mathbb {Z}^2$. A binary image A, or simply, a point set, is a subset of $\mathbb {Z}^2$ and its complement is the set of grid points not in A: $\overline{A}=\mathbb {Z}^2\setminus A$. Either A or $\overline{A}$ must be finite. A grayscale image is a mapping from a finite rectangular domain of $\mathbb {Z}^2$ to $\mathbb N$.

1.1 Morphology Operators on Binary Sets

Definition 1

(Erosion, dilation). The erosion and dilation of the set A by B, referred to as the structuring element, respectively denoted by $\varepsilon _{B}\left( A\right)$ and $\delta _{B}\left( A\right)$, are the sets:

$$\begin{aligned} \varepsilon _{B}\left( A\right)&=\{p:\left( {B}\right) _{p} \subseteq A\},&\delta _{B}\left( A\right)&=\bigcup _{p\in B}\left( {A}\right) _{p} \end{aligned}$$

(1)

where $\left( {B}\right) _{p}$ denotes the translation of B by p: $\left( {B}\right) _{p}=\left\{ p+q:q\in B\right\}$.

Definition 2

(Opening). The opening of A by B, denoted by $\gamma _{B}(A)$, is the erosion of A with B followed by the dilation with B:

$$\begin{aligned} \gamma _{B}(A)=\delta _{B}\left( \varepsilon _{B}\left( A\right) \right) . \end{aligned}$$

(2)

$\gamma _{B}(A)$ can be seen as the union of all translated images $\left( {B}\right) _{p}$ included in A, where the details of A smaller than B are deleted:

$$\gamma _{B}(A)=\bigcup _{p:\left( {B}\right) _{p}\subset A}\left( {B}\right) _{p}\subset A.$$

Definition 3

(Granulometry). An opening-based granulometry on A is a family of openings of A with a sequence of structuring elements $B_r$:

$$\begin{aligned} \left( \gamma _{B_r}(A),r\in \mathbb N_*\right) \end{aligned}$$

(3)

with $B_1=\{O\}$ and $\forall r\ge 0,B_{r+1}=\gamma _{B_r}(B_{r+1})$. Obviously, $\gamma _{B_1}(A)=A$.

A granulometry acts on an image like a sequence of sieves, the value r specifies the size of details to be suppressed. For the sake of simplicity, we write $\gamma _r$ instead of $\gamma _{B_r}$ when the sequence of $B_r$ is clear from the context. Figure 1 pictures a family of openings with the increasing sequence of disks shown in Fig. 2.

Fig. 1.

Openings with the sequence of increasing disks depicted in Fig. 2.

Fig. 2.

Increasing sequence of disks of the octagonal distance used as structuring elements for morphological openings.

Definition 4

(Granulometry function, pattern spectrum). Given a sequence of sets $B_r$, the granulometry function on A, $G_A$, maps a positive r with the cardinality of $\gamma _r(A)$ and the pattern spectrum maps r with the cardinality of $\gamma _r(A)\setminus \gamma _{r+1}(A)$:

$$\begin{aligned} \forall r>1,\mathrm{G}_A(r)&=\left| \gamma _r(X))\right| ,\end{aligned}$$

(4)

$$\begin{aligned} \forall r>1,\mathrm{PS}_A(r)&=\left| \gamma _r(X))\right| -\left| \gamma _{r+1}(X))\right| . \end{aligned}$$

(5)

Figure 3 depicts the granulometric spectrum corresponding to the sequence of morphological openings pictured in Fig. 1.

Fig. 3.

Granulometric spectrum of the sample image shown in Fig. 1a

In the following, we will use distance disks $B_<(p,r)=\big \{q:d(q,p)<r\big \}$ of some distance d as structuring elements for the opening. This definition and properties of distances imply that $B_<(p,0)=\emptyset$ and $B_<(p,1)=\{p\}$.

Definition 5

(Distance transform). The distance transform $\mathrm {DT}_A$ of the set A is the function that maps each point p to its distance from the complement of A and, equivalently, to the radius of the largest disk of center p included in A:

$$\begin{aligned} \mathrm {DT}_A(p)&=\min \big \{d(q,p):q\in \overline{A} \big \}\end{aligned}$$

(6)

$$\begin{aligned}&=\max \big \{r:B_<(p,r)\subset A\big \}. \end{aligned}$$

(7)

It is clear from Eq. (7) that $\mathrm {DT}_A(p)\ge r\iff B_<(p,r)\subset A\iff p\in \varepsilon _{B(O,r)}\left( A\right)$. Consequently, the erosion of A with $B_<(O,r)$ can be obtained by thresholding $DT_A$, i.e., $DT_A$ summarizes all erosions of A with $B_<(O,r),\forall r$.

A Neighborhood Sequence (NS) distance is a distance which disks are produced by dilations with a few neighborhoods, typically the 4-neighborhood $\mathcal {N}_4$ and the 8-neighborhood $\mathcal {N}_8$, arranged in a sequence N(r): $B_<(O,r+1)= \delta _{\mathcal {N}_{\mathrm{N}(r)}}\left( B_<(O,r)\right)$. The 2D NS distance produced by the strict alternation between $\mathcal {N}_4$ and $\mathcal {N}_8$ is called the octagonal distance (see Fig. 2). A so-called “translated” NS distance transform with asymmetric disks produced by translated neighborhoods $\mathcal {N}'_4$ and $\mathcal {N}'_8$, as shown in Fig. 4c, was described previously [6, 7]. It was proven that those distance transforms are equivalent in terms of included disks as in Eq. (7) [7]. An efficient implementation of the translated NS distance transform that requires a single scan of the image and whose memory requirements are limited to a single line of image was described in [6].

Definition 6

(Opening function). Consider d, a distance or pseudo-distance, e.g., asymmetric. The opening function, or opening transform, $\mathrm {OT}_A$ of the set X is a function that maps each point p to the radius of the largest opening disk of d included in X that contains p:

$$\begin{aligned} \begin{array}{rcl} \mathrm {OT}_A(p)= & {} \max \big \{r:\exists q,B(_)<(q,r)\subset A, p\in B(_)<(q,r)\}\,. \end{array} \end{aligned}$$

(8)

Equivalently, $\mathrm {OT}_A(p)$ is the smallest parameter r for which p is filtered out by the morphological opening of the image:

$$\begin{aligned} \mathrm {OT}_A(p)=\min \big \{r:p\not \in \gamma _r(X)\}.\end{aligned}$$

(9)

Fig. 4.

Opening transform $\mathrm {OT}_A$, distance transform $\mathrm {DT}_A$ and translated distance transform $\mathrm {DT}^\prime _A$ of a binary image A. All use the octagonal distance. Thresholding $\mathrm {DT}_A$ and $\mathrm {DT}^\prime _A$ results in equal sets up to a translation. Delimited shapes correspond to pixels whose value is 3 or more. Section 2 describes how the area and perimeter of these shapes are linked.

1.2 Computation of the Granulometry Function

Morphological Openings. A simple algorithm consists in computing morphological openings for each structuring element $B_r$. Even with an efficient implementation of the opening operator, say linear with the size of the image domain N, the order of complexity is $\mathcal O(RN)$, also linear with the maximal radius R. If the structuring elements are chosen to be the disks of a distance, a speedup is achieved by computing the distance transform, that provides all the erosions at once. However, a dilation still has to be performed for each r and the order of complexity remains the same.

Opening Function. The granulometry function can be easily deduced from the opening function which it is simply the histogram.

$$\begin{aligned} \forall r>0,f_A(r) = {{\,\mathrm{Card}\,}}\left( \left\{ p:\mathrm {OT}_A(p)=r\right\} \right) \end{aligned}$$

(10)

However, despite the similarity between the distance and opening transforms $\mathrm {DT}_A$ and $\mathrm {OT}_A$, as in Eqs. (7) and (8), there is no linear time algorithm to compute the latter (except for special cases [11]).

In principle, it suffices to fill each disk $B(_)<(p,r)$ with its radius (given in the distance map) while keeping the maximal radius when several disks overlap. Because disks overlap, each pixel is visited more than once which increases the computational cost. In order to reduce the computation, we can observe that non maximal disks do not need to be filled since they are, by definition, included in larger disks. It is then enough to fill the maximal disks obtained in the medial axis $\mathrm {MAT}_A$ (or $\mathrm{CMB}_A$, center of maximal balls) [1]. The number of overlaps decreases without reaching 0, in general.

2 Estimation of the Pattern Spectrum from the Distance Transform

In this section we propose a method whose purpose is to estimate the pattern spectrum from the distance transform without actually computing the opening transform nor dilations of the distance transform. For a given r, the eroded image by the ball of radius r, $\varepsilon _r(A)$, is simply the result of binarizing the distance map with threshold r. Instead of actually dilating $\varepsilon _r(A)$, the measures on the opening $\gamma _r(A)$, i.e. the dilated of $\varepsilon _r(A)$, are extrapolated from those of $\varepsilon _r(A)$.

After first exposing the linear behavior of the Steiner’s formula in Sect. 2.1, we present a continuous interpolation of binary images based on $2\times 2$ cells where a version of Steiner’s formula can be adapted to our problem in Sect. 2.2. Section 2.3 shows how measures in the distance transform can be extracted from $2\times 2$ local configurations (as in [3, 4]), taking advantage of both the additive property and the continuous model. Finally, Sect. 2.4 presents the algorithm.

2.1 Rationale

Steiner’s formula for a convex K and a Euclidean ball of radius B(r) links the area and perimeter of the dilated convex $\delta _{B(r)}\left( K\right)$ with the area $\mathcal A$ and perimeter $\mathcal P$ of K and B(r) [10]. In the 2D case:

$$\begin{aligned} \mathcal A(\delta _{B(r)}\left( K\right) )&=\mathcal A(K)+r\mathcal P(K)+\mathcal A(B(r)),\end{aligned}$$

(11)

$$\begin{aligned} \mathcal P(\delta _{B(r)}\left( K\right) )&=\mathcal P(K)+P(B(r)). \end{aligned}$$

(12)

Federer later extended this result to non convex compacts and introduced the notion of reach, i.e., the maximal value of r for which Steiner’s formula holds [2]. A further generalization to erosions was stated by Matheron: Steiner’s formula holds for non convex compacts K and the erosion with B(r) if and only if K is the mathematical opening of a compact with B(r) [5]. This is equivalent to the notion a negative reach. If, equivalently to considering the erosion of K, we consider the dilation of its complement, $\overline{K}=\mathbb {R}^2\setminus K$, then, by duality of the erosion and dilation:

$$\begin{aligned} \mathcal A(\delta _{B(r)}\left( \overline{K}\right) )&=\mathcal A(\overline{K})+r\mathcal P(\overline{K})-\mathcal A(B(r)),\end{aligned}$$

(13)

$$\begin{aligned} \mathcal P(\delta _{B(r)}\left( \overline{K}\right) )&=\mathcal P(\overline{K})-P(B(r)) \end{aligned}$$

(14)

where, by convention, $\forall K,\mathcal A(\overline{K})=-\mathcal A(K)$.

For any compact, or complement of a compact, equal to an opening with B(r), we build a vector of characteristics $\begin{pmatrix}\mathcal A,\mathcal P,\chi \end{pmatrix}^T$ with positive area $\mathcal A$, $\chi =1$ for compacts and negative area and $\chi =-1$ for complements. Then

$$\begin{aligned} \begin{pmatrix}\mathcal A'\\ \mathcal P'\\ \chi \end{pmatrix}= \begin{pmatrix} 1&{}\,r&{}\,\pi r^2\\ 0&{}\,1&{}\,2\pi r\\ 0&{}\,0&{}\,1\end{pmatrix} \begin{pmatrix}\mathcal A\\ \mathcal P\\ \chi \end{pmatrix} \end{aligned}$$

(15)

where $\begin{pmatrix}\mathcal A',\mathcal P',\chi \end{pmatrix}^T$ is the corresponding vector of characteristics of the set dilated by B(r). Consider now a composite shape K, an opening with B(r). We treat K as a mixture of its outer and inner contours and build a vector of characteristics $\begin{pmatrix}\mathcal A,\mathcal P,\chi \end{pmatrix}^T$ where $\mathcal A$ is the area of K, i.e., the sum of the areas of the interiors of outer contours minus the areas of the interiors of the inner contours, $\mathcal P$ the total perimeter, and $\chi$ is the count of outer minus inner contours. It is clear that $\begin{pmatrix}\mathcal A,\mathcal P,\chi \end{pmatrix}^T$ is the sum of the vectors of characteristics of the contours. If the contours do not merge by dilation with K, then, by linearity, Eq. (15) holds.

2.2 Continuous Model

Let A be a binary image. Every $2\times 2$ configuration of points of A is mapped to its convex hull $\mathfrak C_A(k,l)=\mathrm {conv}(\{(k,l),(k+1,l),(k,l+1),(k+1,l+1)\}\cap A)$. We define $\mathfrak A$, the continuous interpolation of A, as the union of all $\mathfrak C_A(k,l)$:

$$\begin{aligned} \mathfrak A=\bigcup _{k,l}\mathfrak C_A(k,l)=\bigcup _{k,l}\mathrm {conv}(\{(k,l),(k+1,l),(k,l+1),(k+1,l+1)\}\cap A). \end{aligned}$$

(16)

It is clear that A is the Gauss discretization of $\mathfrak A$: $A=\mathfrak A\cap \mathbb {Z}^2$. Figure 5 pictures two binary images, their continuous interpolation, inner and outer contours.

Consider P, the polygon with edges $\varvec{e_i},1\le i\le n$, obtained by an opening with the line segment $L=[O,O+\varvec{v}]$. Then

$$\begin{aligned} \mathcal A(\delta _{L}\left( P\right) )=\mathcal A(P)+\frac{1}{2}\sum _{i=1}^n|\det (\varvec{e_i},\varvec{v})| \end{aligned}$$

(17)

whether P is an outer or inner contour. In the former case, $\delta _{L}\left( P\right)$ contains 2 more edges of directions $\varvec{v}$ and $-\varvec{v}$ than P. In the latter case, areas are negative and the dilated contour contains two less edges of directions $\varvec{v}$ and $-\varvec{v}$ than P.

Practically, for neighborhood sequence distances, each structuring element, namely the 4- and 8-neighborhood, combine two dilations with orthogonal vectors, respectively $(-1,1)$ and (1, 1) on one hand and (0, 2) and (2, 0) on the other hand. Let $\nu _0(\mathfrak A)$ (resp. $\nu _1(\mathfrak A)$) be the number of horizontal and vertical (resp. diagonal) vectors. For each continuous analogue of a contour (outer or inner), we build a vector of characteristics $V(\mathfrak A)=\big (\mathcal A(\mathfrak A), \nu _0(\mathfrak A), \nu _1(\mathfrak A), \chi (\mathfrak A)\big )^T$ where $A(\mathfrak A)>0$ and $\chi S(\mathfrak A)=1$ for outer contours, $A(\mathfrak A)<0$ and $\chi S(\mathfrak A)=-1$ for inner contours. Then, a dilation of $\mathfrak A$ with $\mathcal {N}_4$ (resp. $\mathcal {N}_8$) corresponds to the matrix product of $M_{\mathcal {N}_4}$ (resp. $M_{\mathcal {N}_8}$) with $V(\mathfrak A)$ where:

(18)

Let B be a disk produced by m 4-neighborhoods and n 8-neighborhoods, then for each contour K of $\varepsilon _{B}\left( A\right)$:

$$\begin{aligned} V(\delta _{B}\left( K\right) )=(M_{\mathcal {N}_4})^m(M_{\mathcal {N}_8})^nV(K) \end{aligned}$$

(19)

because the dilated contour $\delta _{B}\left( K\right)$ is an opening with B. Like in the previous section, Eq. (19) holds for a mixture of contours if the contours do not merge during the dilation, otherwise some quantities may be counted multiple times. For example, parameter vectors measured on the opening transform (Fig. 4a and Table 2 left) are correctly extrapolated from parameter vectors measured on the distance transform (Fig. 4c and Table 2 right). For instance, for radius 3:

$$ M_{\mathcal {N}_4}\cdot M_{\mathcal {N}_8} \cdot \begin{pmatrix}5/2&2&3&1\end{pmatrix}^T =\begin{pmatrix}59/2&10&7&1\end{pmatrix}^T. $$

Fig. 5.

Each line presents a binary image (left), its continuous interpolation (second column), its outer (third column) and inner (right) contours. Their vectors of characteristics are, from (b) to (d) and (f) to (h): $V_a=(22,10,18,0)^T$, $V_b=(34,8,10,1)^T$, $V_c=(-12,2,8,-1)^T$, $V_d=(50,10,18,0)^T$, $V_e=(54,8,14,1)^T$ and $V_f=(-4,2,4,-1)^T$ with $V_a=V_b+V_c$ and $V_d=V_e+V_f$. Each shape on the first row is the eroded of the shape underneath with the 4-neighborhood or its continuous interpolation, and conversely, each shape on the second row is the dilated of the shape above with the 4-neighborhood or its continuous interpolation. Thus, $V_d=M_{\mathcal {N}_4}V_a$, $V_e=M_{\mathcal {N}_4}V_b$ and $V_f=M_{\mathcal {N}_4}V_c$, with $M_{\mathcal {N}_4}$ from Eq. (18).

2.3 Local Configurations

For an arbitrary image A, the following quantities are measured on it continuous analogue $\mathfrak A$: the area $\mathcal A(\mathfrak A)$, the number of horizontal and vertical elementary edges $\nu _0(\mathfrak A)$, the count of diagonal elementary edges $\nu _1(\mathfrak A)$ and the difference between the number of outer and inner contours, i.e., the Euler-Poincaré characteristic $\chi (\mathfrak A)$. For these four measures, the additivity property hold:

$$\begin{aligned} \forall \mu \in \big \{\mathcal A,\nu _0,\nu _1,\chi \big \},\forall \mathfrak A,\mathfrak B,\mu (\mathfrak A\cup \mathfrak B)=\mu (\mathfrak A)+\mu (\mathfrak B)-\mu (\mathfrak A\cap \mathfrak B). \end{aligned}$$

(20)

As a consequence, $\mathcal A$, $\nu _0$, $\nu _1$ and $\chi$ can be evaluated on $\mathfrak A$ from the examination of subsets of $\mathfrak A$ and their intersections. We divide $\mathfrak A$ along elementary grid cells, according to its definition in Eq. (16). Then

$$\begin{aligned} \mu (\mathfrak A)=&+\sum _{k,l}\mu (\mathfrak C_A(k,l))\\&-\sum _{k,l}\mu (\mathfrak C_A(k,l)\cap \mathfrak C_A(k+1,l)) +\mu (\mathfrak C_A(k,l)\cap \mathfrak C_A(k,l+1))\\&+\sum _{k,l}\mu (\mathfrak C_A(k,l)\cap \mathfrak C_A(k+1,l)\cap \mathfrak C_A(k,l+1)\cap \mathfrak C_A(k+1,l+1)). \end{aligned}$$

Note that each $\mathfrak C_A(k,l)$ determines the intersections with its neighbors, e.g., $\mathfrak C_A(k,l)\cap \mathfrak C_A(k+1,l)=\mathfrak C_A(k,l)\cap \mathrm {conv}(\{(k+1,l),(k+1,l+1)\})$. We chose to include in the measures associated to each $\mathfrak C_A(k,l)$, the shares of the intersections with its neighbors, equally distributed (see Fig. 6). This avoids a separate count for measures on $\mathfrak C_A(k,l)$ and on their intersections.

By convention, $\mathfrak C_A(k,l)$ is associated with the configuration index $2^0\mathbf 1_A(k,l)+2^1\mathbf 1_A(k+1,l)+2^2\mathbf 1_A(k,l+1)+2^3\mathbf 1_A(k+1,l+1)$, where $\mathbf 1_A$ is the indicator function of the set A. Table 1 presents the 16 possible configurations and their contributions to measures.

Fig. 6.

Total contribution to the Euler-Poincaré characteristic in a grid cell including the shares of intersections with neighbor cells. Left. the cell contribution to the Euler-Poincaré characteristic for itself is 1. The intersection is non-empty in the four edges. The contribution of each one is $-1$, shared by two cells. The intersection is non-empty in three vertices. The contribution of each one is 1, shared by four cells. The global measure for the cell is then $\chi =1 - 4\times \frac{1}{2} + 3\times \frac{1}{4}=-\frac{1}{4}$. Middle. $\chi =1 - 3\times \frac{1}{2} + 2\times \frac{1}{4}=0$. Right. $\chi =1 - 2\times \frac{1}{2} + 1\times \frac{1}{4}=\frac{1}{4}$.

Table 1. Contributions of the 16 local configurations to image measures.

2.4 Algorithm

The measures $\varvec{\mu }_\varepsilon =(\mathcal A,\nu _0,\nu _1,\chi )$ have now to be computed for each erosion of A with $B_<(O,r)$, i.e., in each binarization of $\mathrm {DT}^\prime _A$ (or $\mathrm {DT}_A$) with threshold r.

Each local configuration $\mathcal C_A(k,l)$ is determined by a $2\times 2$ patch in $\mathrm {DT}^\prime _A$ and the threshold r. The multiple choices of r generate at most 5 configurations according to how r compares to the values of the 4 pixels. To avoid iterating over all values of r, Algorithm 1 enumerates the 5 possible values of the configuration index c(r) and update the first differences of $\varvec{\mu }_\varepsilon$, $d\varvec{\mu }_\varepsilon :r\mapsto \varvec{\mu }_\varepsilon (r)-\varvec{\mu }_\varepsilon (r-1)$. Let v contains the 4 values of the patch in increasing order of weights of the configuration index and $\sigma$ be the reciprocal permutation of the ranks of elements of v in increasing order i.e., $\forall i\in [0,2], v(\sigma _i)\le v(\sigma _{i+1})$. Starting with $r=0$ and $c(r)=15$, c(r) is decreased by $2^{\sigma _i}$ each time r exceeds some $v(\sigma _i)$. For each interval $[r_1,r_2)$ of r, $d\varvec{\mu }_\varepsilon (r1)$ (resp. $d\varvec{\mu }_\varepsilon (r2)$) has to be increased (resp. decreased) with the value of $\mu _c$. These computation are performed in the first loop of Algorithm 1. An algorithm similar in principle was presented by Snidaro and Foresti [9], although limited to measuring the Euler-Poincaré characteristic.

Table 2. Parameter vectors measured (left) on the opening transform depicted in Fig. 4a, (right) on the distance transforms shown in Figs. 4b and c. The number of pixels $\#$ is deduced from the other four values and Pick’s theorem [8]: $\#=\mathcal A+\frac{\nu _1+\nu _2}{2}+\chi$.

The second loop of Algorithm 1 simultaneously accumulates $d\varvec{\mu }_\varepsilon (r)$ in $\varvec{\mu }_\varepsilon (r)$, computes the matrix for the dilation with the disk $B_<(O,r)$ using either $M_{\mathcal {N}_4}$ or $M_{\mathcal {N}_8}$ according to which neighborhood is used for a specific radius and, finally, extrapolates the result vector of measures $\hat{\varvec{\mu }}_\gamma (r)$ from $\varvec{\mu }_\varepsilon (r)$.

The overall time complexity of the method, including the asymmetric distance transform [6], is of order $\mathcal O(KL+R)$, where KL is the number of pixels and R the maximal radius. The space complexity is of order $\mathcal O(K+R)$, as its only needs one line of K pixels for the distance transform, two lines of pixels and an array of size R for the extraction of vectors of characteristics. Moreover, the examination of $2\times 2$ image patches can be merged with the single scan of the distance transform, leading to a streaming algorithm where computations are performed continuously while the data is input.

3 Results and Discussion

The method was tested on sections of tomographic reconstructions of composite materials and compared to the results of a morphological granulometry. Figure 8 shows an example of such an image. Numerical results are presented in Table 3 and pictured in Fig. 7. On this class of images, discrepancies on the pixel count are low (at most 2.25 %). They attain at most 13.33 % on the count of vertical and horizontal segments. Higher disparities appear on the Euler-Poincaré characteristic. This quantity measures the difference between the number of components and holes. It is highly sensitive to merging contours, which occur frequently in images with a high count of touching components like the one displayed in Fig. 8. Note that the pixel count is the only relevant measure here for granulometry. The other quantities are only displayed for completeness.

The source code of the complete method is available at https://github.com/nnormand/DGtalTools-contrib.git.

Table 3. Values $\#$, $\nu _0$, $\nu _1$ and $\chi$ measured on the actual openings of the image in Fig. 8 and the corresponding estimated values, $\hat{\#}$, $\hat{\nu }_0$, $\hat{\nu }_1$ and $\hat{\chi }$, obtained from a translated distance map.

Fig. 7.

Measured and estimated granulometry parameters

Fig. 8.

Tomographic section of a composite material (orthogonal to fibers) and its openings ($r=2$ to 6).

4 Conclusion and Future Works

This article describes a method for computing an estimation of the granulometric spectrum of a binary image by extrapolating parameters measured on the distance transform with a neighborhood-sequence distance. The algorithmic structure makes the method fit for very large amounts of data, or continuously acquired data. The time complexity of the algorithm is linear with the number of pixels for the distance transform, linear with the maximal disk radius for the measure extrapolation. Its space complexity is very low: not more that the result size (the vector of estimated parameters) and two lines of image.

The collision of contours during dilation is the sole reason why the number of pixels can be overestimated. A preliminary analysis shows that such collisions can be identified by saddle points in the distance transform. Further investigation is needed both to upper bound the pixel count error and conclude if and how the effects of contour merging can be neutralized or at least mitigated.