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

ABSTRACT

Minkowski functionals encompass standard geometric parameters such as volume, area, length and the Euler-Poincare 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-Poincare 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-Poincare 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-Poincare characteristic of the structure is approximated by the Euler-Poincare characteristic of its reconstruction. Different types of connectivity may be considered: 4-neighbours, 6-neighbours, 8-neighbours, leading to different approximations of the Euler-Poincare characteristic. In practice, using the Euler formula the Euler-Poincare characteristic of this reconstruction can be computed from 2 x 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-Poincare 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 1 OMITTED]

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-Poincare measure, a key argument is Zahle's definition of Gaussian curvatures for a large class of structures which includes structures with smooth boundaries, as well as polygons and polyhedra (Zahle, 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 [R.sup.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 (Zahle, 1984; 1986a). Such a set is associated with d + 1 Minkowski functionals and measures on [R.sup.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 [intersection] B and for any real-valued function [phi] defined on [R.sup.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 derivative]X of X. Hence U (B) is the length of [partial derivative]X [intersection] B and U([phi]) is the integral of [phi] along [partial derivative]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-Poincare measure. The total mass of [chi] is the Euler-Poincare characteristic of X. The measure [chi] is supported by the boundary [partial derivative]X. For a planar Borel subset B:

[chi] (B) 1 / 2[pi] [[integral].sub.[partial derivative]X[intersection]B] [kappa] (x) dx (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 derivative]X, e.g. when [partial derivative]X is twice continuously differentiable. The general definition of the Euler-Poincare measure for an arbitrary union of sets with positive reach is given in Zahle (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 derivative]X. The measure [bar.b] is also supported by [partial derivative]X. For a Borelian subset B of [R.sup.3], [bar.b](B) is the integral over [partial derivative]X [intersection] 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-Poincare measure:

[chi] (B) = 1 / 4[pi] [[integral].sub.[partial derivative] X [intersection] B] g(x) dx (2)

where g(x) is the Gaussian curvature of [partial derivative]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-Poincare 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-Poincare characteristic of X [intersection] L, the intersection of the line with the structure:

U = [pi] [integral] [chi] (L) dL. (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-Poincare characteristic of the one-dimensional subset X [chi square] 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 [intersection] 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-Poincare measure associated with X [intersection] L.

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

S = 4 [integral] [chi] (L) dL, (4)

and

[bar.b] = [integral] [chi] (P) dP, (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-Poincare measures.

BINARY IMAGES AND DISCRETE MEASURES

A 2D or 3D binary image can be considered as a subset of a rectangular grid [L.sup.d] a where d = 2,3. Such a grid can be written as

[L.sup.d] = [[DELTA].sub.1] Z x ... X [[DELTA].sub.d]Z,

where the d-uple ([[DELTA].sub.1], ..., [[DELTA].sub.d]) defines the pixel or voxel size. A digitization scheme represents how a continuous structure X is digitized into a binary image [X.sub.d]. According to the lattice scheme, [X.sub.d] = X [intersection] [L.sup.d]. According to the covering scheme a grid point x belongs to [X.sub.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 [L.sup.d] a is defined by a sequence of weights [mu] (x), x [member of] [L.sup.d] a and the relationship

[mu] (B) = [summation over (x[member of]B)] [mu] (x), B [subset] [L.sup.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.sub.d]. All approximations will be measures on [L.sup.d], i. e., images of weights. Hence the computations can be seen as image transforms: [X.sub.d] [??] [A.sub.d], [X.sub.d] [??] [U.sub.d] etc. where [A.sub.d], [U.sub.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.sub.d] (B) = [[DELTA].sub.1][[DELTA].sub.2]2 # {[X.sub.d] [intersection] B} (6)

[V.sub.d] (B) = [[DELTA].sub.1][[DELTA].sub.2][[DELTA].sub.3] # {[X.sub.d] [intersection] B} B [subset] [L.sup.d] (7)

where #{[X.sub.d] [intersection] B} is the number of pixels or voxels contained in B.

EULER-POINCARE MEASURES

The approximations [[chi].sub.d] are Euler-Poincare measures computed on polygonal and polyhedral reconstructions [X.sub.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.sub.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 2 OMITTED]

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 [X.sub.r] let us define its envelope [X.sub.r] as the union of all cells of [X.sub.r]. The approximation [[chi].sub.d] is defined as the Euler-Poincare measure associated with [X.sub.r].

Note that the Euler-Poincare measure associated with [X.sub.r] can not be defined by Eqs. 1 and 2 since the boundary of [X.sub.r] is not smooth. Instead one can use Zahle's extension. For polygonal and polyhedral sets the Euler-Poincare measure is supported by the vertices of [X.sub.r]. For a vertex x where [X.sub.r] is convex:

[[chi].sub.d] (x) = [alpha] ([X.sub.r],x] / [s.sub.d], (8)

where [alpha] ([X.sub.r],x) is the normal angle of [X.sub.r] at x, defined as the solid angle of the set of vectors opposite to the tangent cone, and [s.sub.d] is the surface of the unit ball in dimension d.

When [X.sub.r] is not convex at x, [[chi].sub.d](x) depends on the tangent cone of [X.sub.r] at x. It is a matter simple but tedious to show that [[chi].sub.d](x) can be expressed as:

[[chi].sub.d](x) = [summation over (C[member of] [X.sub.r] (x))] [(-1).sup.dimC [alpha] (C,x) / [S.sub.d] (9)

where [X.sub.r] (x) is the subset of cells containing x. The measure [[chi].sub.d] is obtained by:

[[chi].sub.d] (B) = [summation over (x [member of B])] [[chi].sub.d] (x). (10)

It is easy to check that [[chi].sub.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].sub.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].sub.d](x).

By summing [X.sub.d](x) on all vertices of [X.sub.r] one obtains the Euler-Poincare characteristic of [X.sub.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].sub.d] can be rewritten as the alternate sum of numbers of vertices (V), edges (E), faces (F) or solids (S) of the reconstruction. We obtain the classical Euler formula:

[[chi].sub.d] = #V -#E +#F -#S. (11)

The approximation of the Euler-Poincare 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 3 OMITTED]

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

[pi] / [summation over (k)] [chi] ([L.sub.k]) / [lambda] ([L.sub.k]), (12)

where k = 0, ..., n - 1 are indices of the n directions, and [lambda] ([L.sub.k]) is the density of the lines parallel to the line [L.sub.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.sub.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 Mucklich (2000) was used: each direction is represented by a point on the sphere, and weights [c.sub.k], k = 1, ..., 13 are chosen according to the relative area of the corresponding spherical Voronoi cell, normalized such that [SIGMA][c.sub.k] = 1 (Fig. 4). Note that line and plane densities vary with the direction.

[FIGURE 4 OMITTED]

The full discretization of Crofton formula requires the computation of the Euler-Poincare measure restricted to lines or planes. Let us consider the first approximation (Eq. 12) for the perimeter. For each direction k, let [X.sub.r,k] be the one-dimensional reconstruction of X [intersection] [L.sub.k]. This reconstruction consists of vertices of [X.sub.d] [intersection] [L.sub.k] and edges joining neighbour vertices. The measures [chi]([L.sub.k]) are approximated by the Euler-Poincare measures [[chi].sub.d,k] associated with the envelops [X.sub.r,k]. The final approximation of the perimeter measure can be written:

[U.sub.d](x) = [pi] / n [summation over (k)] 1 / [lambda]([L.sub.k]) [summation over (C[member of] [X.sub.r,k](x))] [(-1).sup.dimC] [alpha] (C,x) / 2[pi], (13)

where [lambda] ([L.sub.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.sub.d](x) = 4 [summation over (k)] [c.sub.k] / [lambda]([L.sub.k]) [summation over (C[member of] [X.sub.r,k](x))] [(-1).sup.dimC] [alpha] (C,x) / 4[pi], (14)

[[bar.b].sub.d](x) = [summation over (k)] [c.sub.k] / a([P.sub.k]) [summation over (C[member of] [X.sub.r,k](x))] [(-1).sup.dimC] [alpha] (C,x) / 4[pi], (15)

where [lambda] ([L.sub.k]) is the density of lines in the k-th direction, i.e., 1/[lambda] ([L.sub.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.sub.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.sub.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 Mucklich, 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.sub.r] or [X.sub.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.sub.r](x) in Eq. 9 or [X.sub.rk](x) in Eqs. 13-15 is based only on the 3 x 3 (resp. 3 x 3 x 3) neighbourhood of x.

Each 3 x 3 (resp. 3 x 3 x 3) neighbourhood is decomposed in 2 x 2 (resp. 2 x 2 x 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 5 OMITTED]

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.

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.

Results about area and volume measure approximations are already available (see, e.g., Kieu and Mora 2005; 2006). For a large class of planar structures, the approximation error for [A.sub.d] tends to 0 as fast as [([[DELTA].sub.1][[DELTA].sub.2]).sup.3/4]. Concerning spatial structures, the approximation error for [V.sub.d] tends to zero as fast as [([[DELTA].sub.1][[DELTA].sub.2][[DELTA].sub.3]).sup.2/3.

The convergence of the Euler-Poincare 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-Poincare measure in 2D, the special case where X is a disk has been considered. A function [psi] on the plane has been chosen, and the convergence of [[chi].sub.d]([psi]) towards [chi]([psi]) has been assessed. The chosen function [psi] 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].sub.d] ([psi]) converges as the grid resolution tends to 0 to a value which differs from [chi]([psi]) (Fig. 6).

Zahle (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-Poincare 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-Poincare measures. In view of the non-convergence of the latter measure, it is likely that the mean breadth measure approximation does not converge.

[FIGURE 6 OMITTED]

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:

[pi] / n [n.summation over (k=1)] [H.sub.k], (16)

where [H.sub.k] denotes the total projection length measure on a line parallel to the k-th direction. As shown by (Moran, 1966), this error [epsilon] is bounded by -0.2146 x U [less than or equal to] [epsilon] [less than or equal to] 0.1107 x U for 2 directions and -0.0519 x U [less than or equal to] [epsilon] [less than or equal to] 0.0262 x U for 4 directions.

The same types of errors are involved in the surface area measure approximation. By simple calculation the directional errors e can be bounded by -0.3333 x S [less than or equal to] [epsilon] [less than or equal to] 0.1547 x S for 3 directions and -0.0733 x S [less than or equal to] [epsilon] [less than or equal to] 0.0227 x S for 13 directions.

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 7 OMITTED]

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 mechnical 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 8 OMITTED]

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-Poincare 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 (1).

(Accepted Tune 1, 2007)

REFERENCES

Cheniclet C, Rong WY, Causse M, Frangne, N, Bolling L, Cade J, Renaudin JP (2005). Cell expansion and endoreduplication show a large genetic variablilty 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.

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

Kieu 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 Euter-Poincare characteristic of spatial images. J Microsc 198:54-62.

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

Ohser J, Nagel W, Schladitz K (2002). The Enter 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 Enter number on thin slabs. Adv Appl Prob 36:715-24.

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.

Zahle M (1984). Curvature measures and random sets, I. Math Nach 119:327-39.

Zahle M (1986a). Curvature measures and random sets, II. Probab Theory Re171:37-58.

Zahle M (1986b). Integral and current representation of Federer's curvature measures. Arch Math 46:557-71.

Zahle M (1987). Curvatures and currents for unions of sets of positive reach. Geom Ded 23:557-67.

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

(1) http://wwwjouy.inrafr/unites/miaj/article.php3?id_article=683

DAVID LEGLAND (1), KIEN KIEU (2) AND MARIE-FRANCOISE DEVAUX (1)

(1) INRA, UR1268 Research Unit on Biopolymers--Interactions and Assemblies, BP 71627, F-44316 Nantes, (2) INRA, UR341 Applied Mathematics and Computer Science Unit, F-78352 Jouy-en-Josas

e-mail: david.legland@nantes.inra.fr
```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.

shape U [U.sup.2.sub.d] [U.sup.4.sub.d]

disk 251.33 251.31 251.32
ring 376.99 376.99 377.02
square 0 120.00 94.24 112.66
square 10 -- 109.96 120.51
square 20 -- 120.95 122.68
square 30 -- 127.23 120.82
square 40 -- 130.38 115.73
square 45 -- 131.95 113.18

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 azimut [phi].

shape [phi] [theta] S [S.sup.3.sub.d] [S.sup.13.sub.d]

ball 0 0 5026.5 5024.5 5023.6
cube 0 0 2400.0 1600.0 2160.4
cube 45 0 -- 1989.3 2146.8
cube 45 45 -- 2489.3 2249.0
torus 0 0 5921.7 5713.3 5895.8
torus 45 0 -- 5893.3 5878.2
torus 45 45 -- 6028.0 5906.2

Table 3. Differences between exact and approximated values for the mean

[[bar.b].sup.3 [[bar.b].sup.13
shape [psi] [theta] [bar.b] .sub.d] .sub.d]

ball 0 0 40.0 39.9 39.9
cube 0 0 30.0 20.0 27.2
cube 45 0 -- 25.0 27.2
cube 45 45 -- 31.3 28.4
torus 0 0 47.1 40.0 46.2
torus 45 0 -- 48.0 46.2
torus 45 45 -- 48.7 47.0
```
COPYRIGHT 2007 Slovenian Society for Stereology and Quantitative Image Analysis
No portion of this article can be reproduced without the express written permission from the copyright holder.