Image Anal Stereol 2007;26:83-92 Original Research Paper COMPUTATION OF MINKOWSKI MEASURES ON 2D AND 3D BINARY IMAGES David Legland1, Kien Kieu2 and Marie-Françoise Devaux1 1INRA, UR1268 Research Unit on Biopolymers – Interactions and Assemblies, BP 71627, F-44316 Nantes, 2INRA, 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 × 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. 83 Legland D et al: Computation of Minkowski measures on 2D and 3D binary images 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%. 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 Rd (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 Rd. 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 n B and for any real-valued function 0 defined on R2, A(0) is the integral of 0 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 dX of X. Hence U(B) is the length of dXf]B and U( along dX. The total mass of U is the boundary length of X. The third Minkowski measure is equal to n%, where % is the Euler-Poincaré measure. The total mass of x is the Euler-Poincaré characteristic of X. The measure % is supported by the boundary dX. For a planar Borel subset B: X(B) = — K(x)dx (1) 2n JdxnB where k(x) refers to the (signed) curvature. Note that the definition above holds if the curvature is defined for all point of the boundary dX, e.g. when dX 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 dX. The measure b is also supported by dX. For a Borelian subset B of R3, b(B) is the integral over dX n B of the mean curvature divided by 2%. If X is convex, the total mass is the so-called mean breadth of X. Up to an alternative renormalization the measure b can also be considered as a length measure for thick 84 Image Anal Stereol 2007;26:83-92 fibres. The last measure % is again the Euler-Poincaré measure: X(B) = — f g(x)dx (2) 47T JdxnB where g(x) is the Gaussian curvature of dX 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, b — and %, 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 % 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 n L, the intersection of the line with the structure: U = n I %{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-Poincaré characteristic of the one-dimensional subset X n 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 %{L) by X(BnL) in the right-hand side. Hence Cauchy formula (Eq. 3) holds also when U is the perimeter measure instead of the whole perimeter. Then, x(L) must be considered as the Euler-Poincaré measure associated withXHL. The 3D extension of Cauchy formula is called Crofton formula. Special cases of Crofton formula are S = 4 fx(L) dL, (4) and b —=[x(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 n 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 b can be considered either as total (scalar) parameters or as well as measures. Hence the measures U, S and 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 hd where d = 2,3. Such a grid can be written as l/=A1Zx...xAdZ, where the d-uple (A1,...,Ad) defines the pixel or voxel size. A digitization scheme represents how a continuous structure X is digitized into a binary image Xd. According to the lattice scheme, Xd = X n I/. According to the covering scheme a grid point x belongs to Xd 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 ß on the grid hd is defined by a sequence of weights p(x),x G hd and the relationship /i(B) = L/i(x), Bel/, xeB i.e., m 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 Xd. All approximations will be measures on hd, i.e., images of weights. Hence the computations can be seen as image transforms: Xd ^ Ad, Xd ^ Ud etc. where Ad, Ud denote the discrete approximations of A andU. Approximations of the planar area measure A (d = 2) and the volume measure V (d = 3) are straightforward: Ad(B) = A1A2#{Xdr\B} (6) Vd{B) = A1A2A3#lXdHB} Bel/ (7) where #{Xd n B} is the number of pixels or voxels contained in B. 85 Legland D et al: Computation of Minkowski measures on 2D and 3D binary images EULER-POINCARÉ MEASURES The approximations ?d are Euler-Poincaré measures computed on polygonal and polyhedral reconstructions Xr. 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 Xd, 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). (a) (b) 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. 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 Xr, let us define its envelope Xr as the union of all cells of Xr. The approximation ?d is defined as the Euler-Poincaré measure associated with Xr. Note that the Euler-Poincaré measure associated with Xr can not be defined by Eqs. 1 and 2 since the boundary of Xr 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 Xr. For a vertex x where Xr is convex: ?(Xr,x) ?d(x)= , (8) sd where ?(Xr,x) is the normal angle of Xr at x, defined as the solid angle of the set of vectors opposite to the tangent cone, and sd is the surface of the unit ball in dimension d. When Xr is not convex at x, ?d(x) depends on the tangent cone of Xr at x. It is a matter simple but tedious to show that ?d(x) can be expressed as: dimC ?(C,x) ?d (x)= C??Xr(x)(-1) sd , (9) where Xr(x) is the subset of cells containing x. The measure ?d is obtained by: ?d(B)= ??d(x). (10) x?B It is easy to check that ?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, ?d(x) equals 1/8 for the corner of a cube. Using 26-adjacency on a cubic grid leads to 15 different values for ?d(x). By summing ?d(x) on all vertices of Xr, one obtains the Euler-Poincaré characteristic of Xr 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 ?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: ?d =#V -#E+#F-#S . (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. 86 Image Anal Stereol 2007;26:83-92 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. 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. For instance, the perimeter expressed as the integral (Eq. 3) is approximated by: Wk) (12) where k = 0,...,n - 1 are indices of the n directions, and ?(Lk) is the density of the lines parallel to the line Lk. 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 ck = 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 ck, k = 1,...,13 are chosen according to the relative area of the corresponding spherical Voronoi cell, normalized such that ?ck = 1 (Fig. 4). Note that line and plane densities vary with the direction. (a) (b) 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. 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 JTr fe be the one-dimensional reconstruction of XHLk. This reconstruction consists of vertices of Xd n Lk and edges joining neighbour vertices. The measures %(Lk) are approximated by the Euler-Poincaré measures Xd,k associated with the envelops X1,k. The final approximation of the perimeter measure can be written: (13) where X(Lk) 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: Sd(x) %Wk) ? (-1)dimC C?Xr,k(x) a(C,x) bA ( x) = / ck a(Pk) ? (-1)dimC C?Xr,k(x) a(C,x) (14) (15) where ?(Lk) is the density of lines in the k-th direction, i.e., 1/?(Lk) 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(Pk) 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 ck are the weights given to each direction, and are proportional to the area of Voronoi cells computed on 87 Legland D et al: Computation of Minkowski measures on 2D and 3D binary images 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 Xr or Xr,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 Xr(x) in Eq. 9 or Xr,k(x) in Eqs. 13–15 is based only on the 3×3 (resp. 3×3×3) neighbourhood of x. Each 3 × 3 (resp. 3 × 3 × 3) neighbourhood is decomposed in 2×2 (resp. 2×2×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. (b) Fig. 5. (a) 3 × 3 neighbourhood of a pixel (b) its decomposition in 4 2 × 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. 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 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 Ud2 Ud4 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 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 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 ? and the azimut ?. shape