<figure><img src="image_1.png" /></figure>

HAL
open science

# Computation of Minkowski measures on 2D and 3D binary images

David Legland, Kiên Kiêu, Marie Francoise M. F. Devaux

<figure><img src="image_2.png" /></figure>

## To cite this version:

David Legland, Kiên Kiêu, Marie Francoise M. F. Devaux. Computation of Minkowski measures on 2D and 3D binary images. Image Analysis & Stereology, 2007, 26, pp.83-92. 10.5566/ias.v26.p83-92. hal-02664450

HAL Id: hal-02664450

https://hal.inrae.fr/hal-02664450

Submitted on 31 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COMPUTATION OF MINKOWSKI MEASURES ON 2D AND 3D BINARY IMAGES

DAVID LEGLAND¹, KIÊN KIÊU² AND MARIE-FRANÇOISE DEVAUX¹

¹INRA, UR1268 Research Unit on Biopolymers – Interactions and Assemblies, BP 71627, F-44316 Nantes,

²INRA, UR341 Applied Mathematics and Computer Science Unit, F-78352 Jouy-en-Josas

e-mail: david.legland@nantes.inra.fr

*(Accepted June 1, 2007)*

## ABSTRACT

Minkowski functionals encompass standard geometric parameters such as volume, area, length and the Euler-Poincaré characteristic. Software tools for computing approximations of Minkowski functionals on binary 2D or 3D images are now available based on mathematical methods due to Serra, Lang and Ohser. Minkowski functionals can not be used to describe spatial heterogeneity of structures. This description can be performed by using Minkowski measures, which are local versions of Minkowski functionals. In this paper, we discuss how to extend the computation of Minkowski functionals to the computation of Minkowski measures. Approximations of Minkowski measures are computed using filtering and look-up table transformations. The final result is represented as a grey-level image. Approximation errors are investigated based on numerical examples. Convergence and non convergence of the measure approximations are discussed. The measure of surface area is used to describe spatial heterogeneity of a synthetic structure, and of an image of tomato pericarp.

Keywords: 3D images, Crofton formula, discrete measures, local estimation, Minkowski measures, polyhedral reconstruction.

## INTRODUCTION

Standard parameters used in quantitative analysis of spatial structures are the area and the perimeter in 2D, volume, surface area and mean breadth in 3D. Another parameter is the Euler-Poincaré characteristic, related to the topology of the structure. These parameters form the so-called Minkowski functionals.

Measuring Minkowski functionals in 2D or 3D digitized images is not a trivial task. For instance, measuring perimeter by counting the number of boundary pixels may be quite inaccurate, see the example shown in Fig. 1. Another common problem is the computation of the number of connected components in a binary image. Results depend heavily on the type of connectivity used. Topological parameters measured on digital images may not fulfil standard relationships and statements holding for continuous case.

The Euler-Poincaré characteristic is a standard connectivity parameter. For a planar structure, it is equal to the number of its connected components minus its number of holes. In 2D digital space, its computation requires the reconstruction of the structure by a set of non-overlapping vertices, edges and polygonal faces. The Euler-Poincaré characteristic of the structure is approximated by the Euler-Poincaré characteristic of its reconstruction. Different types of connectivity may be considered: 4-neighbours,

6-neighbours, 8-neighbours, leading to different approximations of the Euler-Poincaré characteristic. In practice, using the Euler formula the Euler-Poincaré characteristic of this reconstruction can be computed from $2 \times 2$ configuration counts (Serra, 1982; Rosenfeld and Kak, 1982).

The Cauchy formula which relates perimeter and number of intercepts provides a simple method for perimeter computation (Serra, 1982). Counting intercepts in 2D binary images is straightforward: it is again a configuration count.

Extensions to 3D images have been derived (see Nagel *et al.*, 2000, Ohser *et al.*, 2002 for the Euler-Poincaré characteristic and Lang *et al.*, 2001 for all other parameters). Again, all needed computations are configuration counts. Methods for computing any parameter belonging to the family of the Minkowski functionals both in 2D and 3D images are now available.

In order to characterize the spatial distribution of the 2D or 3D structure, one may compare Minkowski functional values measured in different regions. This approach can be further extended to consider Minkowski measures instead of Minkowski functionals. In particular, it is possible to build maps and profiles for each parameter.
<figure><img src="image_4.png" /><figcaption>Fig. 1. The digitized disk is the set of pixels contained in the initial disk (represented as black filled circles). A polygonal reconstruction is obtained by joining adjacent pixels of the digitized disk by edges. The perimeter of the disk may be approximated by the total length of the external edges (thick lines). For a disk, this perimeter approximation tends to the perimeter of the square with the same diameter as the digitization grid tends to 0. Hence the relative approximation error is about 25%.</figcaption></figure>

In this paper the methods for computing Minkowski functionals on 2D and 3D images are extended to the computation of Minkowski measures. Concerning the Euler-Poincaré measure, a key argument is Zähle’s definition of Gaussian curvatures for a large class of structures which includes structures with smooth boundaries, as well as polygons and polyhedra (Zähle, 1984; 1986a; 1987). The computation of higher-dimensional Minkowski measures related e.g. to perimeter or surface area is based on the Cauchy and Crofton formulae as for the corresponding Minkowski functionals.

The article is organized as follows. We first recall mathematical definitions of the Minkowski measures. Then the reconstruction of the structure as a complex of convex cells is presented. We develop the computation of Minkowski measures approximations based on this reconstruction, then we present an algorithmic implementation based on configuration counts. The convergence of the Minkowski measures approximations is investigated. Finally, the measure of surface area is used to describe the spatial heterogeneity of a synthetic structure and of a tomato pericarp.

# MINKOWSKI MEASURES

## DEFINITIONS

Let $X$ be a subset of the $d$-dimensional Euclidean space $\mathbb{R}^d$ ($d = 1,2,3,...$) belonging to the class of

unions of sets with positive reach. Finite unions of convex sets, and sets whose boundary is a finite union of twice continuously differentiable manifolds, fall into this class (Zähle, 1984; 1986a). Such a set is associated with $d+1$ Minkowski functionals and measures on $\mathbb{R}^d$. Hence a 2D structure ($d=2$) defines 3 Minkowski functionals/measures, while a 3D structure ($d=3$) defines 4 Minkowski functionals/measures. In 2D and 3D, the Minkowski functionals coincide up to known multiplicative constants with standard geometric parameters. Intrinsic volumes are an alternative set of definitions for the same quantities.

For a planar structure $X$, the first Minkowski functional is the area $A$ of $X$. The corresponding Minkowski measure is just defined as the restriction to $X$ of the area (i.e., Lebesgue) measure. If this measure is also denoted by $A$, then for any planar Borel set $B$, $A(B)$ is the area of $X \cap B$ and for any real-valued function $\phi$ defined on $\mathbb{R}^2$, $A(\phi)$ is the integral of $\phi$ onto $X$. Also note that the total mass of the measure $A$ is equal to the area of $X$.

The second Minkowski measure is equal to $U/2$, where $U$ is the length measure restricted to the boundary $\partial X$ of $X$. Hence $U(B)$ is the length of $\partial X \cap B$ and $U(\phi)$ is the integral of $\phi$ along $\partial X$. The total mass of $U$ is the boundary length of $X$.

The third Minkowski measure is equal to $\pi\chi$, where $\chi$ is the Euler-Poincaré measure. The total mass of $\chi$ is the Euler-Poincaré characteristic of $X$. The measure $\chi$ is supported by the boundary $\partial X$. For a planar Borel subset $B$:

$$ \chi(B) = \frac{1}{2\pi} \int_{\partial X \cap B} \kappa(x) \, dx \quad (1) $$

where $\kappa(x)$ refers to the (signed) curvature. Note that the definition above holds if the curvature is defined for all point of the boundary $\partial X$, e.g. when $\partial X$ is twice continuously differentiable. The general definition of the Euler-Poincaré measure for an arbitrary union of sets with positive reach is given in Zähle (1984).

A spatial structure yields 4 Minkowski measures. Let $V$ be the volume (i.e., Lebesgue) measure restricted to $X$ and let $S$ be the surface area measure restricted to the boundary $\partial X$. The measure $\bar{b}$ is also supported by $\partial X$. For a Borelian subset $B$ of $\mathbb{R}^3$, $\bar{b}(B)$ is the integral over $\partial X \cap B$ of the mean curvature divided by $2\pi$. If $X$ is convex, the total mass is the so-called mean breadth of $X$. Up to an alternative renormalization the measure $\bar{b}$ can also be considered as a length measure for thick
fibres. The last measure $\chi$ is again the Euler-Poincaré measure:

$$ \chi(B) = \frac{1}{4\pi} \int_{\partial X \cap B} g(x) \, dx \quad (2) $$

where $g(x)$ is the Gaussian curvature of $\partial X$ at $x$, i.e., the product of the two main curvatures at $x$.

All measures introduced above are defined over the continuous plane or 3D space. In practice data on structures are most often available as discrete images. Information is available only for a finite subset of points usually organized as a regular grid. We consider here only binary images, a particular case in which information is either true or false. Binary images are usually obtained after a segmentation procedure, i.e., each pixel or voxel is either inside or outside the structure of interest. Our goal is the computation of Minkowski measures based on discrete binary images.

This is a non trivial problem for measures $U, S, \bar{b}$ and $\chi$, because their support is restricted to the boundary of the structure. Discrete approximation of curves and surfaces is still an active field of research (Klette and Rosenfeld, 2004). Below we present a method for computing the measure $\chi$ from binary 2D or 3D images. The computation of the other measures is derived both from Euler-Poincaré measure computations and from the Crofton formula.

## CROFTON FORMULA

The total perimeter $U$ of a structure can be expressed by integrating, for all (affine) lines $L$ in the plane, the Euler-Poincaré characteristic of $X \cap L$, the intersection of the line with the structure:

$$ U = \pi \int \chi(L) \, dL. \quad (3) $$

where $dL$ is the measure over spatial lines invariant with respect to rigid motion and normalized such that the mass of lines hitting the unit circle equals 2. Note that the Euler-Poincaré characteristic of the one-dimensional subset $X \cap L$ is just the intercept number of $X$ by $L$. Eq. 3 is known as Cauchy formula. This formula concerns the whole perimeter of $X$. However it is easy to see that the formula holds also locally: $U$ can be replaced by $U(B)$ in the left-hand side and $\chi(L)$ by $\chi(B \cap L)$ in the right-hand side. Hence Cauchy formula (Eq. 3) holds also when $U$ is the perimeter measure instead of the whole perimeter. Then, $\chi(L)$ must be considered as the Euler-Poincaré measure associated with $X \cap L$.

The 3D extension of Cauchy formula is called Crofton formula. Special cases of Crofton formula are

$$ S = 4 \int \chi(L) \, dL, \quad (4) $$

and

$$ \bar{b} = \int \chi(P) \, dP, \quad (5) $$

where $dL$ is the measure over spatial lines invariant with respect to rigid motion and normalized such that the mass of lines hitting the unit sphere is equal to $\pi$ and $dP$ is the measure over planes invariant by rigid motion and normalized such that the mass of planes hitting the unit sphere is equal to 2. Again $S$ and $\bar{b}$ can be considered either as total (scalar) parameters or as well as measures.

Hence the measures $U, S$ and $\bar{b}$ can be expressed as integrals with respect to lines or planes of Euler-Poincaré measures.

## BINARY IMAGES AND DISCRETE MEASURES

A 2D or 3D binary image can be considered as a subset of a rectangular grid $\mathbb{L}^d$ where $d = 2, 3$. Such a grid can be written as

$$ \mathbb{L}^d = \Delta_1 \mathbb{Z} \times \dots \times \Delta_d \mathbb{Z}, $$

where the $d$-uple $(\Delta_1, \dots, \Delta_d)$ defines the pixel or voxel size. A digitization scheme represents how a continuous structure $X$ is digitized into a binary image $X_d$. According to the lattice scheme, $X_d = X \cap \mathbb{L}^d$. According to the covering scheme a grid point $x$ belongs to $X_d$ if the grid cell centered on $x$ hits $X$. Further details about digitization schemes and so-called digitalizable maps can be found in Serra (1982).

Any measure $\mu$ on the grid $\mathbb{L}^d$ is defined by a sequence of weights $\mu(x), x \in \mathbb{L}^d$ and the relationship

$$ \mu(B) = \sum_{x \in B} \mu(x), \quad B \subset \mathbb{L}^d, $$

i.e., $\mu$ is a weighted sum of Dirac measures. Note that such a discrete measure can be also represented as an image of weights. Below we will see how to compute approximations of the Minkowski measures associated with a continuous structure $X$ from its digitization $X_d$. All approximations will be measures on $\mathbb{L}^d$, i.e., images of weights. Hence the computations can be seen as image transforms: $X_d \mapsto A_d$, $X_d \mapsto U_d$ etc. where $A_d, U_d$ denote the discrete approximations of $A$ and $U$.

Approximations of the planar area measure $A$ ($d = 2$) and the volume measure $V$ ($d = 3$) are straightforward:

$$ A_d(B) = \Delta_1 \Delta_2 \#\{X_d \cap B\} \quad (6) $$

$$ V_d(B) = \Delta_1 \Delta_2 \Delta_3 \#\{X_d \cap B\} \quad B \subset \mathbb{L}^d \quad (7) $$

where $\#\{X_d \cap B\}$ is the number of pixels or voxels contained in $B$.
EULER-POINCARÉ MEASURES

The approximations $\chi_d$ are Euler-Poincaré measures computed on polygonal and polyhedral reconstructions $\mathcal{X}_r$.

Reconstructions considered here are collections of convex polygonal or polyhedral sets of different dimensions, which touch only by their lower-dimensional faces. Let us call the reconstruction elements cells. For convenience, all lower-dimensional faces of each cell are added to the reconstruction. The cells are chosen such that their vertices belong to $X_d$, and such that their edges belong to a given set of adjacencies (Fig. 2).

The 4-adjacency in 2 dimensions and the 6-adjacency in 3 dimensions, which consider only orthogonal neighbours, lead to cubical reconstruction. Each cell is either a vertex, an edge or a face parallel to the directions of the grid, or a cube corresponding to a tile of the grid. This reconstruction was called cubical complex by Khachan et al. (2000) (Fig. 2).

For 8- and 26-adjacencies, all vertices of the grid with a difference of coordinate at most equal to one for each orthogonal direction are considered as neighbours. In 2D, cells may be triangles and squares. In 3D, cells may be polyhedra with 4 to 8 faces (Fig. 2).

<figure><img src="image_5.png" /><figcaption>Fig. 2. Polyhedral reconstructions of a structure using (a) the 6-adjacency and (b) the 26-adjacency. The latter reconstruction contains various convex polyhedra, triangular, square and rectangular faces, and diagonal edges.</figcaption></figure>

Other reconstructions may be considered, for example the 6-adjacency in 2 dimensions, the 18-adjacency, or the 14.1 and 14.2-adjacencies which have been proposed by Ohser et al. (2002).

For a given reconstruction $\mathcal{X}_r$, let us define its envelope $X_r$ as the union of all cells of $\mathcal{X}_r$. The approximation $\chi_d$ is defined as the Euler-Poincaré measure associated with $X_r$.

Note that the Euler-Poincaré measure associated with $X_r$ can not be defined by Eqs. 1 and 2 since the boundary of $X_r$ is not smooth. Instead one can use Zähle's extension. For polygonal and polyhedral sets the Euler-Poincaré measure is supported by the vertices of $X_r$. For a vertex $x$ where $X_r$ is convex:

$$ \chi_d(x) = \frac{\alpha(X_r, x)}{s_d}, \qquad (8) $$

where $\alpha(X_r, x)$ is the normal angle of $X_r$ at $x$, defined as the solid angle of the set of vectors opposite to the tangent cone, and $s_d$ is the surface of the unit ball in dimension $d$.

When $X_r$ is not convex at $x$, $\chi_d(x)$ depends on the tangent cone of $X_r$ at $x$. It is a matter simple but tedious to show that $\chi_d(x)$ can be expressed as:

$$ \chi_d(x) = \sum_{C \in \mathcal{X}_r(x)} (-1)^{\dim C} \frac{\alpha(C, x)}{s_d}, \qquad (9) $$

where $\mathcal{X}_r(x)$ is the subset of cells containing $x$. The measure $\chi_d$ is obtained by:

$$ \chi_d(B) = \sum_{x \in B} \chi_d(x). \qquad (10) $$

It is easy to check that $\chi_d(x)$ equals 1 for an isolated vertex, 1/2 for an edge extremity. In 2D, other possible values are 1/4 for square corners and 1/8 for triangle corners when using 8-adjacency on a square grid. In 3D, $\chi_d(x)$ equals 1/8 for the corner of a cube. Using 26-adjacency on a cubic grid leads to 15 different values for $\chi_d(x)$.

By summing $\chi_d(x)$ on all vertices of $X_r$, one obtains the Euler-Poincaré characteristic of $X_r$ expressed as a double sum on vertices and cells. The summation order can be inverted. As the sum of normalized normal angles of a convex cell equals 1, the total mass of $\chi_d$ can be rewritten as the alternate sum of numbers of vertices ($\mathcal{V}$), edges ($\mathcal{E}$), faces ($\mathcal{F}$) or solids ($\mathcal{S}$) of the reconstruction. We obtain the classical Euler formula:

$$ \chi_d = \#\mathcal{V} - \#\mathcal{E} + \#\mathcal{F} - \#\mathcal{S}. \qquad (11) $$

The approximation of the Euler-Poincaré measure strongly depends on the choice of the adjacency system. As shown by Ohser et al. (2002), different adjacencies can yield very different results. One possibility is to choose the adjacency depending on the geometry of the structure. For example, one can choose 4-adjacency or 6-adjacency for thick structures, and 8-adjacency or 26-adjacency for structures presenting thin edges or isthmuses.
## LENGTH AND SURFACE AREA MEASURES

Measures of perimeter, surface area and mean breadth require discretization of Crofton formula (Eqs. 3–5). The integrals over lines or planes are approximated by sums over discrete sets of lines or planes.

In the plane, the set of lines can be discretized into horizontal and vertical lines (Fig. 3). An alternative is to use also the 2 diagonal directions. As the density of the lines varies with their direction, they must be weighted accordingly.

<figure><img src="image_6.png" /><figcaption>Fig. 3. Discrete set of lines. Dots are grid vertices. Discretization with 2 directions considers only horizontal and vertical lines (in black). Discretization with 4 directions also considers diagonals (in orange). Line density is higher in the diagonal directions.</figcaption></figure>

For instance, the perimeter expressed as the integral (Eq. 3) is approximated by:

$$ \frac{\pi}{n} \sum_k \frac{\chi(L_k)}{\lambda(L_k)}, \quad (12) $$

where $k = 0, ..., n-1$ are indices of the $n$ directions, and $\lambda(L_k)$ is the density of the lines parallel to the line $L_k$. This density depends on the pixel sizes and may vary with the line orientation.

Similarly, the set of lines (resp. planes) in the 3D space can be discretized into lines parallel to (resp. planes normal to) the 3 main directions of the grid. Each discrete direction is associated with a weight $c_k = 1/3$, $k = 1, 2, 3$. A more precise solution is to use 13 directions, following the (undirected) directions between 2 vertices of the unit cube. Note that this discretization is not isotropic. The choice of Ohser and Mücklich (2000) was used: each direction is represented by a point on the sphere, and weights $c_k$, $k = 1, ..., 13$ are chosen according to the relative area of the corresponding spherical Voronoi cell, normalized such that $\sum c_k = 1$ (Fig. 4). Note that line and plane densities vary with the direction.

<figure><img src="image_7.png" /><figcaption>Fig. 4. Partition of the unit sphere into 6 areas (a), and 26 areas (b), corresponding to the directions of the 6 or 26 neighbours of the central voxel. The areas of the spherical sectors are the same in the case of 6-neighbourhood.</figcaption></figure>

The full discretization of Crofton formula requires the computation of the Euler-Poincaré measure restricted to lines or planes. Let us consider the first approximation (Eq. 12) for the perimeter. For each direction $k$, let $\mathcal{X}_{r,k}$ be the one-dimensional reconstruction of $X \cap L_k$. This reconstruction consists of vertices of $X_d \cap L_k$ and edges joining neighbour vertices. The measures $\chi(L_k)$ are approximated by the Euler-Poincaré measures $\chi_{d,k}$ associated with the envelops $X_{r,k}$. The final approximation of the perimeter measure can be written:

$$ U_d(x) = \frac{\pi}{n} \sum_k \frac{1}{\lambda(L_k)} \sum_{C \in \mathcal{X}_{r,k}(x)} (-1)^{\dim C} \frac{\alpha(C, x)}{2\pi}, \quad (13) $$

where $\lambda(L_k)$ is the distance between two pixels on a discrete line oriented in the $k$-th direction. The final approximations for the surface area measure and the mean breadth measure in 3D are given by:

$$ S_d(x) = 4 \sum_k \frac{c_k}{\lambda(L_k)} \sum_{C \in \mathcal{X}_{r,k}(x)} (-1)^{\dim C} \frac{\alpha(C, x)}{4\pi}, \quad (14) $$

$$ \bar{b}_d(x) = \sum_k \frac{c_k}{a(P_k)} \sum_{C \in \mathcal{X}_{r,k}(x)} (-1)^{\dim C} \frac{\alpha(C, x)}{4\pi}, \quad (15) $$

where $\lambda(L_k)$ is the density of lines in the $k$-th direction, i.e., $1/\lambda(L_k)$ is the area of a tile of the lattice formed by the intersection of the lines with a plane perpendicular to them. The term $a(P_k)$ corresponds to the density of planes with the $k$-th direction, i.e., $1/a(P,k)$ is the distance between 2 neighbour plates. The coefficients $c_k$ are the weights given to each direction, and are proportional to the area of Voronoi cells computed on
the unit sphere, with germs depending on the discrete directions (Ohser and Mücklich, 2000).

Note that for approximations of perimeter and surface area measures, there is a choice on the number of directions to use. For the mean breadth measure, one must choose the number of plane directions and the adjacency type for each plane.

## ALGORITHMIC IMPLEMENTATION

The Minkowski measure approximations can be computed from the global reconstructions $X_r$ or $X_{r,k}$. For large images, the reconstruction of the entire set, and the computation of normal angles for each cell, can lead to huge computation time. In practice, for each pixel (resp. voxel) $x$, the computation of $X_r(x)$ in Eq. 9 or $X_{r,k}(x)$ in Eqs. 13–15 is based only on the $3 \times 3$ (resp. $3 \times 3 \times 3$) neighbourhood of $x$.

Each $3 \times 3$ (resp. $3 \times 3 \times 3$) neighbourhood is decomposed in $2 \times 2$ (resp. $2 \times 2 \times 2$) tiles. A cell may belong to several tiles. In 2D, a pixel belongs to 4 tiles, an isothetic edge (parallel to one of the main directions of the grid) to 2 tiles, a diagonal edge to 1 tile, and a face to 1 tile (Fig. 5). In 3D, a voxel belongs to 8 tiles, an edge to 1, 2, or 4 tiles, depending on its direction, a face to 1 or 2 tiles, and a solid cell to 1 tile. The contribution of each cell is weighted according to the number of tiles it belongs to.

<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) $3 \times 3$ neighbourhood of a pixel (b) its decomposition in $4 \times 2 \times 2$ tiles. The central pixel belongs to the 4 tiles. Vertical edges belong to 2 tiles, the diagonal edge belongs to only one tile. The triangular face belongs to only 1 tile.</figcaption></figure>

Thus, approximations (Eq. 9, Eqs. 13–15) are decomposed as sums of 4 (resp. 8) tile contributions. The whole computation is achieved by linear filtering and look-up table transformations. The objective of the linear filtering is to identify each tile configuration as described in Lang *et al.* (2001). Tile contributions are computed from the filtered image using 4 (resp. 8) look-up table transformations.

## APPROXIMATION ACCURACY

First we provide numerical results about the accuracy of the functional approximations. Then we focus on the convergence of the measure approximations.

In two dimensions, three types of shapes were generated and digitized: disks with radius equal to 40 pixels, rings composed of the difference of two concentric disks with radii 40 and 20 pixels, and rotated squares with side length 30 pixels. The comparison of approximations with exact values is given in Table 1.

<table><caption>Table 1. Differences between exact and approximated values for perimeter of planar shapes: a disk, a ring, a square rotated with 0, 10...45 degrees. Approximations are computed using 2 and 4 directions.</caption><thead><tr><th>shape</th><th>U</th><th>U<sup>2</sup><sub>d</sub></th><th>U<sup>4</sup><sub>d</sub></th></tr></thead><tbody><tr><td>disk</td><td>251.33</td><td>251.31</td><td>251.32</td></tr><tr><td>ring</td><td>376.99</td><td>376.99</td><td>377.02</td></tr><tr><td>square 0</td><td>120.00</td><td>94.24</td><td>112.66</td></tr><tr><td>square 10</td><td>—</td><td>109.96</td><td>120.51</td></tr><tr><td>square 20</td><td>—</td><td>120.95</td><td>122.68</td></tr><tr><td>square 30</td><td>—</td><td>127.23</td><td>120.82</td></tr><tr><td>square 40</td><td>—</td><td>130.38</td><td>115.73</td></tr><tr><td>square 45</td><td>—</td><td>131.95</td><td>113.18</td></tr></tbody></table>

For three-dimensional shapes, spheres with constant radius of $R = 30$ voxels, and cubes with side length $L = 30$ were used. The discretization of a torus obtained by rotating a 10-radius disk along a perpendicular 15-radius circle was also considered. Comparison of true values and their approximations is given in Tables 2 and 3. If differences are small in the case of the ball, it can be bigger for cubes or torus. In these cases, the effect of the directions discretization is important. Note that for the considered structures, it is reduced when using 13 directions instead of 3.

<table><caption>Table 2. Differences between exact and approximated values for surface area of 3D shapes: a ball, a cube and a torus. The orientation of the shapes is given by the colatitude $\theta$ and the azimuth $\phi$.</caption><thead><tr><th>shape</th><th>φ</th><th>θ</th><th>S</th><th>S<sup>3</sup><sub>d</sub></th><th>S<sup>13</sup><sub>d</sub></th></tr></thead><tbody><tr><td>ball</td><td>0</td><td>0</td><td>5026.5</td><td>5024.5</td><td>5023.6</td></tr><tr><td>cube</td><td>0</td><td>0</td><td>2400.0</td><td>1600.0</td><td>2160.4</td></tr><tr><td>cube</td><td>45</td><td>0</td><td>—</td><td>1989.3</td><td>2146.8</td></tr><tr><td>cube</td><td>45</td><td>45</td><td>—</td><td>2489.3</td><td>2249.0</td></tr><tr><td>torus</td><td>0</td><td>0</td><td>5921.7</td><td>5713.3</td><td>5895.8</td></tr><tr><td>torus</td><td>45</td><td>0</td><td>—</td><td>5893.3</td><td>5878.2</td></tr><tr><td>torus</td><td>45</td><td>45</td><td>—</td><td>6028.0</td><td>5906.2</td></tr></tbody></table>
<table><caption>Table 3. Differences between exact and approximated values for the mean breadth of 3D shapes.</caption><thead><tr><th>shape</th><th>φ</th><th>θ</th><th>b&#773;</th><th>b&#773;<sup>3</sup><sub>d</sub></th><th>b&#773;<sup>13</sup><sub>d</sub></th></tr></thead><tbody><tr><td>ball</td><td>0</td><td>0</td><td>40.0</td><td>39.9</td><td>39.9</td></tr><tr><td>cube</td><td>0</td><td>0</td><td>30.0</td><td>20.0</td><td>27.2</td></tr><tr><td>cube</td><td>45</td><td>0</td><td>—</td><td>25.0</td><td>27.2</td></tr><tr><td>cube</td><td>45</td><td>45</td><td>—</td><td>31.3</td><td>28.4</td></tr><tr><td>torus</td><td>0</td><td>0</td><td>47.1</td><td>40.0</td><td>46.2</td></tr><tr><td>torus</td><td>45</td><td>0</td><td>—</td><td>48.0</td><td>46.2</td></tr><tr><td>torus</td><td>45</td><td>45</td><td>—</td><td>48.7</td><td>47.0</td></tr></tbody></table>

Results about area and volume measure approximations are already available (see, e.g., Kiêu and Mora 2005; 2006). For a large class of planar structures, the approximation error for $A_d$ tends to 0 as fast as $(\Delta_1\Delta_2)^{3/4}$. Concerning spatial structures, the approximation error for $V_d$ tends to zero as fast as $(\Delta_1\Delta_2\Delta_3)^{2/3}$.

The convergence of the Euler-Poincaré characteristic was established for the 2D case by Serra (1982), and extended to the 3D for different types of adjacencies (Nagel et al., 2000; Ohser et al., 2002; Rataj, 2004).

In order to assess the asymptotic behaviour of the Euler-Poincaré measure in 2D, the special case where $X$ is a disk has been considered. A function $\varphi$ on the plane has been chosen, and the convergence of $\chi_d(\varphi)$ towards $\chi(\varphi)$ has been assessed. The chosen function $\varphi$ is a triangle function of $\theta$, defined as the angle between the line joining the origin to $x$ and the horizontal axis. The support of this function is a circular sector, whose axis of symmetry is chosen equal to $\pi/8$ and angular extent to $\pi/4$. This numerical experiments shows that $\chi_d(\varphi)$ converges as the grid resolution tends to 0 to a value which differs from $\chi(\varphi)$ (Fig. 6).

Zähle (1990) provides a criterion which guarantees that Minkowski measures of a reconstruction converge towards the true value. This criterion involves the normals of the reconstruction. For the reconstructions considered in this paper, normals have a limited number of directions, which does not increase with the resolution. The same counter-example can be considered in 3D space. Hence the 3D approximation of the Euler-Poincaré measure does not converge either. More general results on the approximation errors would be of interest.

The mean breadth measure approximation involves approximations of 2D Euler-Poincaré measures. In view of the non-convergence of the latter measure, it is likely that the mean breadth measure approximation does not converge.

<figure><img src="image_9.png" /><figcaption>Fig. 6. Non-convergence of the 2D Euler-Poincaré measure approximation. The stars represent $\chi_d(\varphi)$ at different resolutions computed using the 8-adjacency, averaged over 100 realizations. See the main text for description of $\varphi$. The horizontal line shows the true value of $\chi(\varphi)$. The dashed lines show the interval containing 95% of the measurements.</figcaption></figure>

The perimeter measure is concerned with 3 levels of discretization error. At the highest level, the set of line directions is discretized into 2 or 4 directions. Next, the set of lines in given direction is discretized by a series of parallel lines. The line spacing is determined by the grid resolution. At the lowest level, the number of intercepts is approximated by configuration counts. Approximations of the number of intercepts converge as the grid resolution increases. At the intermediate level, the approximation error also converges to zero as the resolution increases. However, the error involved at the highest level does not depend on the spatial resolution. The latter error is expressed as follows:

$$ \frac{\pi}{n} \sum_{k=1}^{n} H_k, \qquad (16) $$

where $H_k$ denotes the total projection length measure on a line parallel to the $k$-th direction. As shown by (Moran, 1966), this error $\varepsilon$ is bounded by $-0.2146 \times U \le \varepsilon \le 0.1107 \times U$ for 2 directions and $-0.0519 \times U \le \varepsilon \le 0.0262 \times U$ for 4 directions.

The same types of errors are involved in the surface area measure approximation. By simple calculation the directional errors $\varepsilon$ can be bounded by $-0.3333 \times S \le \varepsilon \le 0.1547 \times S$ for 3 directions and $-0.0733 \times S \le \varepsilon \le 0.0227 \times S$ for 13 directions.
GRADIENT PROFILE OF SURFACE AREA

Many observed structures are not stationary and exhibit heterogeneity. This heterogeneity can be described by the variations of a geometric parameter with the position. Practical solutions often consist in computing the parameter in a subset of the image, defined by a window, and to make the position of the restricted window vary. The use of measures generalizes this principle and overcomes problems such as the choice of the window size, by producing a map of contributions for each pixel or voxel of the image.

The interpretation of Minkowski measures for heterogeneity analysis is not straightforward. Spatial variations at a scale smaller than the size of the objects of interest are not meaningful. Small scale variations can be removed by smoothing. When investigating spatial heterogeneity along a given direction, it is proposed to plot projections of Minkowski measures onto the axis direction. Below, such plots are refered to as gradient profiles.

The interest of gradient profiles in one direction is shown through a synthetic example of spheres with heterogeneous distribution and through the application to a real image of tomato pericarp.

SYNTHETIC STRUCTURE

Fig. 7 shows a synthetic ball process presenting a strong heterogeneity. The balls are generated with a density linearly increasing along the x-axis, and by keeping a minimal distance between centroids such that balls do not overlap. The radius of the balls are chosen according to the ball density, such that the volume occupied by the balls remains theoretically constant. The structure is discretized such that the diameter of smaller balls corresponds to 8 pixels.

The approximation of the corresponding surface area measure is computed from the discrete image using 3 and 13 discrete directions. Contributions of voxels with the same x-coordinate are summed in order to produce a profile of surface area. The profiles for the true value and the approximations nicely superimpose. The approximation errors are small and do not disturb the interpretation of the spatial heterogeneity.

<figure><img src="image_10.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b)</figcaption></figure>

Fig. 7. (a) An inhomogeneous ball process. The density decreases from left to right, as the typical ball volume increases. A minimal distance between centroids is kept such that the balls do not overlap. (b) The volume and the surface area profiles are estimated along the x-axis and compared with the true measures.

TOMATO PERICARP

The construction of gradient profiles is adapted for describing the changes of morphology in biological structures organized with respect to a reference surface, for example epidermis, organs, or fruits. Tomato pericarp is the fleshy part found around tomato fruit and is made of cells whose morphology varies with depth (Cheniclet, 2005). Quantification of variations of cell morphology is required to study relations with mechanical or sensorial texture properties.

Fig. 8 shows a 3D image of tomato pericarp acquired by confocal microscopy. The image was obtained after merging several adjacent 3D stacks
acquired from the external epidermis to the tomato
center. After noise filtering, cells were segmented
by a 3D watershed operation, resulting into two
components: cell interiors and cell walls. Small cells
were observed near the epidermises and large cells
more or less elongated in the pericarp centre.

<figure><img src="image_12.png" /><figcaption>Fig. 8. A 3D isosurface representation of a binary image obtained after 3D segmentation of a confocal set of tomato pericarp images. The surface area measure is plotted in the radial direction of the pericarp.</figcaption></figure>

The surface measure of cell walls, and the volume
measure of pericarp were computed from the image,
and the gradient profile of surface area density was
extracted according to the distance to the external
epidermis (Fig. 8).

Surface area density is much higher in the region
close to the external epidermis, has a lower value in
the middle of the pericarp, and increases again near
the internal epidermis. As the volume occupied by the
cells is nearly constant, and cells are nearly convex, the
increase of surface area density can be interpreted as a
decrease of cell mean size, and hence as an increase of
cell number.

CONCLUSION

In this paper, methods for approximating
Minkowski measures of 2D and 3D structures have
been developed. These methods are direct extensions
of previous works on Minkowski functionals. By using
measures, it is possible to compute profiles describing
spatial heterogeneity in a given direction.

Convergence results hold for area measure in 2D
and volume measure in 3D. Asymptotic error bounds
are provided for the perimeter measure in 2D and
the surface area measure in 3D. For most practical
cases their accuracy seems sufficient when 4 directions

in 2D or 13 directions in 3D are used. Based on
numerical counter-examples, it is shown that Euler-
Poincaré measure approximations in 2D and 3D, and
mean-breadth measure approximations in 3D do not
converge.

Further results about the accuracy of these measure
approximations are still lacking, making the use
of these approximations hazardous. An alternative
approach is to consider other types of reconstructions
(Klette and Rosenfeld, 2004), whose advantage is that
the number of normal directions increases with the
spatial resolution.

An implementation of the method for the Matlab
software is proposed. The complete functions, as well
as some examples, are provided as a library, available
on the Internet <fnref n="1" />.

REFERENCES

Cheniclet C, Rong WY, Causse M, Frangne, N, Bolling
L, Cade J, Renaudin JP (2005). Cell expansion and
endoreduplication show a large genetic variability in
pericarp and contribute strongly to tomato fruit growth.
Plant Physiol 139:1984-94.

Klette R, Rosenfeld A (2004). Digital Geometry: Geometric methods for digital picture analysis. Morgan Kaufman Publishers.

Khachan M, Chenin P, Deddi H (2000). Polyhedral representation and adjacency graph in *n*-dimensional digital images. Comp Vision Image Understanding 79(3):428–41.

Kiêu K, Mora M (2005). Stereological estimation of mean volume: precision of three simple sampling designs. Technical report 2005-1, INRA MIA. Available online.

Kiêu K, Mora M (2006). Precision of stereological planar area predictors. J Microsc 222(3):201–11.

Lang C, Ohser J, Hilfer R (2001). On the analysis of spatial binary images. J Microsc 203(3):303–13.

Moran PAP (1966). Measuring the length of a curve. Biometrika 53:359-64.

Nagel W, Ohser J, Pischang K (2000). An integral geometric approach for the Euler-Poincaré characteristic of spatial images. J Microsc 198:54–62.

Ohser J, Mücklich F (2000). Statistical Analysis of
Microstructures in Material Sciences. Chichester: J.
Wiley & Sons.

Ohser J, Nagel W, Schladitz K (2002). The Euler number of discretized sets—on the choice of adjacency in homogeneous lattices. Lect Notes Phys 600:275–98.

Rataj J (2004). On estimation of the Euler number on thin slabs. Adv Appl Prob 36:715-24.

<footnote n="1">¹http://www.jouy.inra.fr/unites/miaj/article.php3?id_article=683</footnote>
Rosenfeld A, Kak AC (1982). Digital picture processing, Vol. 2. Academic Press, second edition.

Serra J (1982). Image Analysis and Mathematical Morphology, Volume 1. Reading: Academic Press.

Zähle M (1984). Curvature measures and random sets, I. Math Nach 119:327–39.

Zähle M (1986a). Curvature measures and random sets, II. Probab Theory Rel 71:37–58.

Zähle M (1986b). Integral and current representation of Federer’s curvature measures. Arch Math 46:557–71.

Zähle M (1987). Curvatures and currents for unions of sets of positive reach. Geom Ded 23:557–67.

Zähle M (1990). Approximation and characterization of generalized Lipschitz-Killing curvatures. Ann Global Anal Geom 8:249–60.