Surface-Based Computation of the Euler Characteristic in the BCC Grid

As opposed to the 3D cubic grid, the body-centered cubic (BCC) grid has some favorable topological properties: each set of voxels in the grid is a 3-manifold, with 2-manifold boundary. Thus, the Euler characteristic of an object O in this grid can be computed as half of the Euler characteristic of its boundary ∂O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial O$$\end{document}. We propose three new algorithms to compute the Euler characteristic in the BCC grid with this surface-based approach: one based on (critical point) Morse theory and two based on the discrete Gauss–Bonnet theorem. We provide a comparison between the three new algorithms and the classic approach based on counting the number of cells, either of the 3D object or of its 2D boundary surface.


Introduction
The Euler characteristic is a basic yet important, locally computable topological descriptor of both continuous and discrete shapes.It is equal to the alternating sum of (independent) i-dimensional holes.For shapes in 3D, it is equal to the number of connected components minus the number of tunnels plus the number of cavities.It can be computed as the alternating sum of the number of i-cells in a cell decomposition of the shape.
The commonly used 3D cubic grid has some unfavorable geometric and topological properties, mainly due to the existence of three different types of adjacencies between the cubes.This causes many image processing algorithms to become more involved, needing more computation [22].The body-centered cubic (BCC) grid is one feasible alternative to the traditional 3D cubic one, with many topological advantages.It has been effectively used in computer graphics [3,10,13], discrete geometry [9], tomographic reconstruction from projections [4], physically-based simulations [25], voxelization [15], distances by neighborhood sequences [36] and in various other fields including discretization, ray tracing and ray casting, volume rendering and repairing [7,10,13,16,17,32].A software system for processing and viewing 3D data on this grid has also been developed [29].
The BCC grid is a tessellation of the space R 3 through (space-filling) truncated octahedra.Two voxels in the BCC grid are either disjoint, or they share a whole face.This ensures that the boundary ∂ O of an object (a finite set of voxels) O in the BCC grid is always a 2-manifold (and O is a 3-manifold).Therefore, the Euler characteristic of O can be computed as half the Euler characteristic of its boundary ∂ O, which is likely to have significantly fewer cells than O in many practical cases.
Here, we address the computation of the Euler characteristic in the BCC grid with this surface-based approach.Apart from the known algorithm based on counting the boundary cells, we propose to use two frameworks for computing the Euler characteristic of a polyhedron (the boundary ∂ O of an object O) in this grid.
The first framework is the piecewise linear [1] critical point Morse theory [30,33], which provides relationships between the topology of a manifold M and the critical points of a scalar field f defined on M. The second is the discrete version of the Gauss-Bonnet theorem, which states that the sum of angular vertex deficits in a polyhedron is equal to 2π times its Euler characteristic.We propose two algorithms based on this framework.
We provide a detailed comparison of the proposed algorithms, and a comparison of them with the algorithm implementing the classic approach based on counting the number of cells, in a volume-based or surface-based version (i.e., on the 3D object or on its 2D boundary).
In particular, one of the two new algorithms based on the Gauss-Bonnet theorem is the fastest surface-based method, with execution times equal to about 40% of the naive 2D cell counting.The surface-based approach to the Euler characteristic computation is convenient over the volume-based one when the number of boundary faces is less than 25% (if we include the cost of boundary extraction) or less than 66% (if we consider only the Euler characteristic computation) than the number of voxels composing the object.

Basic Notions
A k-cell in R n (with 0 ≤ k ≤ n) is a homeomorphic image of the closed unit k-ball { p ∈ R k : ||p|| ≤ 1}.A cell k-complex Γ is a collection of i-cells, with max{i} = k, that fit nicely together: the boundary of each cell and the intersection of any two cells (if non-empty) is composed of cells of lower dimension.We are interested in the cases n = 3 and k = 2, 3 that are used to represent solid objects and surfaces in R 3 .
A subset I ⊆ R n is a k-manifold (with boundary) if each of its points has a neighborhood homeomorphic to R k (or to In R 3 , we are interested in 3-manifolds with boundary, representing a solid object, and 2-manifolds, representing surfaces.The boundary surface of a 3-manifold object is a (closed) 2-manifold without boundary.
The set of points in a cell complex Γ in R 3 is denoted as |Γ |.If |Γ | is a manifold, we say that Γ is a manifold cell complex.

The Cubic Grid
The 3D cubic grid is the regular tessellation of R 3 into closed unit cubes centered at points in Z 3 , with faces parallel to the coordinate planes [20].Three types of adjacency relation are defined between the cubes in the grid, depending on their intersection.Two cubes are face-, edge-, or vertex-adjacent if they share at least one face, edge or ver-  [26] if it has no critical vertices (incident with exactly two strictly vertex-adjacent black or exactly two strictly vertex-adjacent white cubes) and no critical edges (incident with exactly two strictly edge-adjacent black and exactly two strictly edgeadjacent white cubes).
For an object O, the associated cubical complex Q is a cell 3-complex that consists of the cubes in O together with all of their square faces, edges and vertices.An object O is well-composed if and only if |Q| is a 3-manifold.The two definitions of well-composed 3D objects (discrete, through the critical configurations and continuous, through the notion of a manifold) are equivalent [2].

The BCC Grid
The body-centered cubic (BCC) [19] grid is a Voronoi tessellation of R 3 associated with points in (2Z) 3 ∪ (2Z + 1) 3 , which can be seen as the centers of cubes of two interlaced cubic grids, rescaled by factor 2. The cubes of the first (second) grid are centered at points with even (odd) Cartesian coordinates (see Fig. 1a).The voxel of the BCC grid is a truncated octahedron, with eight regular hexagonal faces and six square faces, 36 edges and 24 vertices (see Fig. 1b).
An object O in the BCC grid is defined similarly to an object in the cubic grid: it is a finite set of (black) voxels in the grid.Its associated cell 3-complex Q is composed of all voxels in O and all their faces, edges and vertices.
The BCC grid has only one type of adjacency relation: if two voxels share a vertex or an edge, they also share a Fig. 2 Two BCC voxels which are not disjoint share an entire a square or b hexagonal face face (see Fig. 2).Due to this property, there are no critical configurations in the BCC grid, and every object in this grid is well-composed.The shared face of two adjacent voxels may be hexagonal or square, and the two voxels are called hex-adjacent and quad-adjacent, respectively.

The Euler Characteristic
The Euler characteristic is a basic topological (homological) invariant, extensively used in many application domains, such as topological data analysis, image processing and pattern recognition and classification, to name just a few.It can be defined in several equivalent ways, and on different spaces, e.g., topological manifolds, cell complexes or polyhedra.
For a cell 3-complex Γ with c 0 vertices, c 1 edges, c 2 faces and c 3 3-cells (voxels), the Euler characteristic χ(Γ ) is equal to For a general 3-complex Γ , this definition assumes that two cells are adjacent if they share at least one vertex.If applied to a cubical complex Q, it assumes vertex-adjacency for Q.
Alternatively, if Γ has β 0 connected components, β 1 tunnels and β 2 cavities, the Euler characteristic χ(Γ ) is equal to β 0 − β 1 + β 2 [21].A cavity is formally a bounded connected component of the background, and intuitively a hole inside the object, completely surrounded by it.A tunnel is intuitively a hole traversing the object from side to side, while a formal definition requires notions from homology theory.Intuitively, the number of tunnels is equal to the maximal number of non-separating cuts that can be made in the object.
In 3D, there is a connection between the Euler characteristic of a manifold 3-complex Γ and the Euler characteristic of its boundary ∂Γ [5,14,28,31], namely and the Euler characteristic of ∂Γ can be computed as where, here, c 0 , c 1 and c 2 denote the number of vertices, edges and faces of ∂Γ .

Morse Theory
Morse theory studies the relationship between the topological shape of a manifold and (the critical points of) a function defined on the manifold, both in the continuous [30,33] and in the piecewise linear [1,11] setting.For a C 2 -differentiable real-valued function f defined over a closed compact manifold domain M ⊆ R d , a point p ∈ M is a critical point of f if all first order partial derivatives of f (in a suitable coordinate system around p) vanish at p.The function f is a Morse function if the Hessian matrix Hess p f of the second derivatives of f at p is non-singular.Usually, it is supposed that there are no two critical points with the same function value.This property allows analyzing the topology changes at each critical point separately.
The number of negative eigenvalues of Hess p f is called the index of a critical point p.The corresponding eigenvectors point in the directions in which f is decreasing.A critical point p is a minimum or a maximum if it has index 0 or d, respectively.Otherwise, if the index of p is λ, 0 < λ < d, p is a λ-saddle.
Let us consider a closed orientable 2-dimensional manifold M in R 3 , which is the boundary of a 3-manifold Γ , i.e., M = ∂Γ .On M, there are three types of critical points: minima, maxima and saddles.The topology (homotopy type) of the lower level sets M ≤α = {x ∈ M : f (x) ≤ α} of f changes only at the critical points.
For example, let M be a (hollow) torus with horizontal axis of symmetry and let f be the vertical elevation.This is a 2-manifold without boundary and its Euler characteristic is 0. Figure 3 presents four stages of the evolution of the lower level set (shown in shaded color) while increasing the elevation.At the lowest point of M, the lower level set is created, as it initially a 2-manifold with boundary, homotopic to a disc (χ = 1).At the lower saddle, it becomes homotopic to a cylinder (χ = 0).At the higher saddle, two boundary circles of the cylinder touch to create a torus with a hole (χ = −1).Until here the lower level set has been a 2-manifold with boundary.The highest point closes up the surface creating the internal cavity (χ = 0 again).
In summary, minima create new connected components and maxima close off (create) cavities.Thus, each extremal critical point increases the Euler characteristic by 1. Saddles either merge two different connected components, or create a hole (a tunnel) by merging two parts of the same connected component.Thus, each saddle decreases the Euler characteristic by 1.
If m λ denotes the number of critical points of f with index λ, then the Euler characteristic χ(M) of M is given by i.e., the number of minima, minus the number of saddles, plus the number of maxima.

The Discrete Gauss-Bonnet Theorem
The Gauss-Bonnet theorem gives a connection between the Gaussian curvature on a manifold M and its topology expressed through its Euler characteristic.It states that where M is a closed 2-manifold without boundary, K G is its Gaussian curvature and the integral is a surface integral over M [35].
The Descartes theorem gives a connection between the angles of the faces of a polyhedron P and its topology.It states that, if P is homeomorphic to the unit 2-sphere S 2 , then where v is a vertex of P and the angular deficit δ(v) is the amount by which the sum of the face angles at v differs from 2π , i.e., where f are the faces incident to the vertex v and α(v, f ) is the internal angle of f at v.
The discrete version of the Gauss-Bonnet theorem extends the Descartes theorem to arbitrary (manifold) polyhedra.It where P is the given polyhedron [23] and v are the vertices of P. Intuitively, the Gaussian curvature in the interior of the faces and edges of P is equal to zero.The curvature is concentrated at the vertices of P and is equal to the vertex deficit.It follows that

Relation Between BCC and Cubic Grids
As common 3D acquisition devices provide data in a cubic grid, it is necessary to convert them to the BCC grid.This can be done in different ways.
A naive conversion can be performed in two ways: either we keep only cubes with all even or all odd coordinates and expand them to BCC voxels, or we insert a new datum at each cube vertex (with interpolated value) and place a BCC voxel at each vertex and each cube center.In the first case, we have half the resolution and we loose data, while in the second case we double the resolution.However there is no warranty that the Euler characteristic is preserved.
A method has been proposed in [7] which performs the conversion in the second way (i.e., by doubling the resolution), but choosing the color assigned to the cube vertices in such a way to preserve the Euler characteristic of the original cubic object according to vertex-adjacency or faceadjacency.Alternatively, Edelsbrunner and Kerber [12] have proposed a method transforming each cubic voxel into a (combinatorial) BCC voxel, by slightly shrinking it from one diagonal direction.The resolution is the same as the origi-nal cubic grid, but in general the Euler characteristic is not preserved.

Related Work
A lot of algorithms have been proposed for the computation of the Euler characteristic of 3D objects in the cubic grid.The main drawback of this grid is that the value of the Euler characteristic depends on the adjacency model considered (classically vertex-adjacency or face-adjacency, but other adjacency models have been proposed as well [34,37]), and it is unique only for 3-manifold objects.Here, we briefly review only the algorithms for computing the Euler characteristic in the cubic grid that are based either on the discrete version of the Gauss-Bonnet theorem, or on Morse theory.
The algorithm by Chen and Rong [6] works on a wellcomposed object O, with manifold boundary Thus, the vertices incident to four boundary faces do not affect the Euler characteristic, and the discrete Gauss-Bonnet theorem implies that The algorithm by Imiya and Eckhardt [18] makes a finer classification of the vertices in the boundary of a wellcomposed object O to obtain the same formula.
The algorithm by Lee at al. [27] works for objects with face-adjacency (or with vertex-adjacency, by considering the complement of the object), not necessarily well-composed.It is based on smoothing the black cubes (slightly inflating them and rounding the corners and edges) and applying the continuous Gauss-Bonnet theorem.It reduces to using a lookup table with vertex contributions to χ(O) for each possible configuration of 2 × 2 × 2 cubes.
The algorithm by Čomić and Magillo [8] is based on counting only the boundary faces and vertices.The vertex count of the vertices where ∂ O is non-manifold is adjusted depending on the chosen adjacency relation (face-or vertexadjacency).

Algorithms for the Computation of the Euler Characteristic in the BCC Grid
We describe five algorithms to compute the Euler characteristic of an object in the BCC grid.The first two algorithms are the classic ones based on cell counting.The other three algorithms are new: one is based on Morse theory and critical points, and two algorithms are based on summing angles.
Our implementation of all algorithms uses the 4-valued coordinate system proposed in [7].In it, every voxel, face, edge and vertex in the BCC grid is represented by four dependent coordinates, whose sum is zero.Of the four coordinates of a BCC voxel centered in the first interlaced cubic grid, two belong to Z 8 and two to Z 8 + 4; for a BCC voxel centered in the second interlaced cubic grid, two coordinates belong to Z 8 + 2 and two to Z 8 + 6.The coordinates of a face are the average of the coordinates of its two incident BCC voxels.The coordinates of an edge (shared by two hex-and one quadface) are the average of the coordinates of its two incident hex-faces.The coordinates of a vertex are the average of the coordinates of its four incident BCC voxels.All adjacency and incidence relations are retrieved by arithmetic operations on the four coordinates of the involved cells.
The input object O is given as a list of (black) BCC voxels.In addition, each voxel of (the portion of) the BCC grid containing O is marked as black or white.In this setting, in order to find, for example, the black voxels adjacent to a given voxel x, we access the 14 neighbors of x and test their color in constant time.The surface ∂ O, needed by the boundary-based algorithms, is represented as a list of faces.

Algorithm Based on Counting Cells
The straightforward classic approach simply computes the alternating sum of the number of i-cells.It admits a volumebased version, which computes the Euler characteristic of a three-dimensional object O in the BCC grid with Formula (1), and a surface-based version, that does it by computing the Euler characteristic of the boundary surface ∂ O of O with Formula (3).

Volume-Based Version (VOL)
This algorithm that we denote as VOL counts each black voxel and its faces, edges and vertices, paying attention to count a face, edge or vertex in just one of the voxels containing it.
Our implementation scans the black voxels and processes each one.Processing of a voxel x consists of the following actions: -count x and mark it; -access the 14 adjacent voxels of x, and check the configuration.Based on the following quantities: -N 4 and N 6 , the number of quad-and hex-faces, respectively, of x, whose other incident voxel is black and already marked; we call them black faces; -N 6,6 , the number of edges of x, shared by two hexfaces, such that their two incident voxels, different from x, are black and already marked; -N 4,6 , the number edges of x, shared by a quad-and a hex-face, such that their two incident voxels, different from x, are black and already marked; -N 4,6,6 , the number of vertices of x (each shared by one quad-and two hex-faces of x) such that their three incident voxels, different from x, are black and already marked; -count 14 − N 4 − N 6 faces for x (from the 14 faces of x, subtract the ones already counted in the adjacent black voxels); -count 36−4N 4 −6N 6 + N 6,6 + N 4,6 edges for x (from the 36 edges of x, subtract the ones belonging to black faces to avoid counting them twice, but those edges shared by two black faces must not be subtracted twice), -count 24 − 4N 4 − 6N 6 + 2(N 6,6 + N 4,6 ) − N 4,6,6 vertices for x (from the 24 vertices of x, subtract the ones belonging to black faces, with similar considerations as for edges, but a vertex may be shared by up to three black faces).

Surface-Based Version (SUR)
This version, that we denote as SUR assumes that the boundary surface of the object (consisting of hexagonal and square faces) is available.The algorithm performs a loop on input boundary faces.For each face f , it counts the face f itself and potentially counts each edge and vertex of f .Suitable attention is payed to count each edge or vertex only at one of the faces it belongs to.For such purpose, in our implementation we consider the lexicographic (total) order of the face coordinates.An edge e is counted only if the current boundary face precedes the other boundary face incident in e, and a vertex v is counted only if the current boundary face is minimum among the boundary faces incident in v.
Thanks to the fact that ∂ O is 2-manifold (each edge belongs to two faces), we could count 1/2 for each edge of a face.Nevertheless, a mechanism such as the total order is still necessary to count the vertices correctly.

Algorithm Based on Morse Theory (MOR)
The new algorithm based on Morse theory, that we denote as MOR, computes the Euler characteristic of the boundary ∂ O by recording the changes in the lower level sets with a suitable elevation function f , defined in such a way that no two vertices have the same function value (as unicity is a usual requirement for Morse functions, see Sect.2.5).This is done by considering the lexicographic order of the (z, y, x) triplets obtained from the coordinates of all the vertices, and function f maps each vertex onto its "elevation" defined as its position in the sorted list of all vertices.Intuitively, this is equivalent to sweeping ∂ O through a set of parallel almost horizontal planes moving in the direction of the positive z axis.Such planes are tilted slightly so that the plane encounters one vertex at a time among the ones with a given z.The change of the Euler characteristic when sweeping over a vertex v depends on the relative position (elevation) of the neighbors of v.Each vertex v in ∂ O is adjacent to either three or four other vertices in ∂ O, and it is processed as follows: -If all neighbors come after v (they are higher than v w.r.The justification for such increments is given in Sect.2.5. The implementation first loops on the faces of the given object boundary ∂ O, extracts the vertices of each hex-face (this guarantees to get each vertex at least once) and puts them into an array.Then, the array is sorted lexicographically and duplicates are removed.
Before starting the main loop on the sorted vertices, χ = 0 is set.Then, a sweep over the vertices of ∂ O, contained in the sorted array, is performed.These are the actions performed while sweeping over a vertex v: -Get the (at most four) boundary faces incident in v; -For each incident boundary face f , get the vertices w 1 , w 2 preceding and following v along the boundary of f .If w 1 < v < w 2 or w 1 > v > w 2 in the lexicographic order, then record one switch; -Let s be the total number of recorded switches at all incident boundary faces of v (s is even and 0 ≤ s ≤ 4).If s = 0 then increment χ ; if s > 2 then decrement χ ; otherwise, leave χ unchanged.
The same formula could be obtained from the known relation 2 by noting that each edge is shared by two faces, because ∂ O is manifold.Therefore, 2c 1 = 4c The algorithm visits all faces in ∂ O and counts the number c q 2 of quad-faces and the number c h 2 of hex-faces.It also counts the number c 0 of vertices in ∂ O, by counting, for each face, all its vertices which are not yet marked, and then marking them.Then, the Euler characteristic is computed as For example, in Fig. 5b we have 44 vertices, 10 quad-faces and 16 hex-faces; therefore, χ(∂ O) = 44 − 10 − 32 = 2.

Computational Complexity
We evaluate the time complexity of the algorithms in the worst case.We denote with n O the number of voxels in the object O and with n ∂ O the number of faces in its boundary surface ∂ O.
Algorithm VOL is in O(n O ) and SUR is in O(n ∂ O ), because they perform a constant amount of work for each voxel in O, and for each face in ∂ O, respectively.SUR does not need to compute the sorted sequence of faces in the lexicographic order as it just tests whether one face comes before or after another (adjacent) one.

Experiments and Results
We implemented all the algorithms for computing the Euler characteristic of an object in the BCC grid, described in  We also implemented an algorithm BND to extract the boundary surface of a given BCC object, which is needed to produce the input for surface-based algorithms to compute the Euler characteristic (i.e., all the above ones, but VOL).
The algorithms have been coded in C language and run on a PC equipped with an Intel CPU i7-2600K CPU at 3.4 Gigahertz with 32 Gigabytes of RAM.
Our experiments are aimed at: 1. Comparing the four surface-based algorithms among them, identifying the most efficient one.2. Comparing the surface-based and the volume-based approach to the Euler characteristic computation, identifying the conditions under which the former or the latter is more efficient.

Comparing the Four Boundary-Based Algorithms
Our test set is obtained by discretizing some shapes from the Digital Shape Workbench, 1 shown in Fig. 6.Such shapes are given as tetrahedral meshes filling an object, or as triangle meshes bounding it.All the shapes have been discretized by first rescaling them so that their bounding box fits a cubic grid of side N (i.e., N 3 cubes) and then considering the BCC grid, where a BCC voxel is centered at each cube center and at each cube vertex.A BCC voxel has been set to black if and only if its center lies inside the object.The used grid sides are N = 50, 100, 150, 200, 250.
Table 1 summarizes the information about the resulting BCC images: the object name, the grid side N , the number n O of black voxels, and the ratio r = n ∂ O /n O (where n ∂ O is the number of boundary faces).As the volume of an object grows roughly with N 3 and the area of its boundary surface with N 2 , the ratio r decreases with resolution and is inversely dependent on the grid side N .The last column in Table 1 shows the value of r × N , which is roughly constant for each shape.
The running times, in milliseconds, are shown in Fig. 7.As expected, the running times are linearly dependent on the size n ∂ O of the boundary surfaces, with a constant that is roughly 0.0062 for SUR, 0.0056 for MOR, 0.0059 for GB1and 0.0024 for GB2 (with variations involving the fourth decimal digit).Therefore, each of the three new boundarybased algorithms represents an improvement over the classic one SUR based on cell counting.
The fastest method is GB2, i.e., the second algorithm using the Gauss-Bonnet theorem.The running time of GB2 is less than 39% of SUR, less than 44% of MOR, and less than 42% of GB1.Algorithms MOR and GB1 are not much faster than SUR, with the former slightly faster than the latter (95% and 96%, respectively).

Comparison with the Volume-Based Approach
As expected, the running times of VOL are linearly dependent on the number n O of black voxels, with a constant roughly equal to 0.0016 (with variations involving the fifth decimal digit).
We compare the fastest boundary-based algorithm, i.e., GB2, with the volume-based one VOL, with the aim of understanding under which conditions the former is convenient over the latter.If the boundary surface is either given, or has to be extracted for other purposes, it is sufficient that the running time of GB2 is smaller than the running time of VOL.Otherwise, it is necessary that the sum of execution times of GB2 and of boundary extraction BND is smaller than VOL. Figure 10 shows the above running times, obtained on the same test set used in Sect.5.1.Referring to Fig. 10, we compare the bottom yellow rectangle (VOL) with either the red rectangle (GB2) or the union of the red and orange rectangles (GB2+BND), i.e., the left part of the top rectangle, or the entire top rectangle.
For an easier comparison, Fig. 8 plots the running times of VOL on the x-axis and of GB2 and GB2+BND on the yaxis.A line connects the five resolutions of each test object.
The cases in which GB2 (or GB2+BND) is faster than VOL are the dots lying below the line y = x.
The boundary-based approach to Euler computation is more convenient than the volume-based one when the boundary surface is small compared to its enclosed volume (in terms of the number of cells), i.e., when the ratio r = n ∂ O /n O is small.The plots in Fig. 9 show the ratio of running times GB2/VOL and (GB2+BND) /VOL as a function of r .Indeed, Fig. 9 shows a linear dependency on r , with factor 1.5 for GB2/VOL and 4.1 for (GB2+BND)/VOL.GB2 is faster than VOL when r < 0.66.The total running time GB2+BND, including boundary extraction, is smaller than that of VOL when r < 0.25.
Small values of r are more likely to happen when the shape is "fat" and when the resolution of the discretization (the

Analysis of the Impact of the Object Shape
In Sect.5.2 we noted that the boundary-based algorithm GB2 outperforms the volume-based one VOL when the ratio r = n ∂ O /n O is small enough.To further analyze the relation between the shape of an object, the resolution of the discretization, and the ratio r , we built three artificial test sets: ), but using another discretization scheme at higher resolutions: the width (difference between outer and inner radius) remains constant in terms of the number of voxels, at all resolution.The cavity gets larger at higher resolutions.

512K voxels
).This probably depends on the hidden constants, and on the larger amount of fixed work in VOL, that becomes relevant for small n O .The ratio of running times (GB2+BND)/VOL ranges from 55 to 60.5, with similar, but more noisy, variations when increasing N .These behaviors are plotted in Fig. 14.Table 2 summarizes the information for the (hollow) spheres with the two discretization schemes, showing the same information as Table 1.
For Ball and BubbleA, BubbleB, BubbleC, the behavior of r is similar to the shapes in   The running times of VOL confirm to be linear in the number n O of black voxels, with a multiplicative factor 0.0016.The running times of GB2 are roughly linear in the number n ∂ O of boundary faces, with a factor that is slightly smaller for Ball, BubbleA, BubbleB, BubbleC (0.0021-0.0022) than for BubbleA', BubbleB', BubbleC' (0.0023-0.0025).
Confirming the results in Sect.5.2, GB2 is faster than VOL when r < 0.7 and GB2+BND is faster than VOL when r < 0.25.The plots in Figs. 15 and 16 compare the running times of GB2 and GB2+BND, with the running times of VOL.
In Ball, r = 0.77 at the lowest resolution, and at higher resolutions it ranges from 0.39 to 0.15; therefore, GB2 is almost always the best algorithm.In the hollow spheres Bub-bleA, BubbleB, BubbleC, r becomes < 0.7 at a resolution that is coarser if the inner cavity is smaller.This confirms that the boundary-based approach is preferable for fat objects.
In BubbleA', BubbleB', BubbleC', the ratio r is almost constant w.r.t. the resolution, and so is the ratio of running times of GB2 and (GB2+BND) over that of VOL.For Bub-bleA' and BubbleB', the running time of GB2 is about 9 times and: 2.2-2.3 times that of VOL, respectively.For BubbleC' the two running times are almost equal (GB2 being slightly slower).
The value of the ratio r = n ∂ O /n O , i.e., the size of the boundary surface of an object O over its inner volume, is the Fig. 16 Running times for the three hollow spheres BubbleA',B',C' with the second discretization scheme.The ratios of running times of GB2 and of GB2+BND as a function of the running time of VOL key point to determine whether the boundary-based approach or the volume-based one will be more efficient.
If the surface of the object is already available, we can compute r and, based on its value, select the appropriate algorithm for computing the Euler characteristic.
As the value of the ratio r is inversely dependent on the grid resolution N , this property suggests that the availability of 3D images with higher and higher resolution will make the surface-based approach more convenient in the future.

Concluding Remarks
We explored the computation of the Euler characteristic of a 3D object with the surface-based approach, i.e., by processing only the boundary cells of the cell 3-complex associated with the object, rather than all cells.This requires that the input object is well-composed, and therefore, we considered objects in the BCC grid.We described the classic method based on cell counting and three new ones: one based on an alternative definition of the Euler characteristic provided by Morse theory, and two based on the discrete Gauss-Bonnet theorem.
Among the proposed methods, the one based on discrete Gauss-Bonnet theorem and cycling on faces, outperforms the classic 2D cell counting algorithm by a factor 0.4 for the tested data sets.Experiments have also shown that the boundary-based approach is faster than the classic volumetric one when the number of boundary cells is less than 66% the number of interior cells (or less than 25%, if the boundary has to be extracted just for this purpose), which is likely to happen if the represented 3D shape has no thin parts and the discretization has a sufficient resolution.
The algorithms proposed in this paper are not limited to objects in the BCC grid.The same approaches can be applied to well-composed objects in grids of any type (including well-composed cubic ones), and in general to 3-manifold polyhedral objects.For example, the Gauss-Bonnet theorem based on faces for a well-composed cubic object would compute χ(∂ O) as c 0 −c 2 (as already noted in [8]).We expect that similar results will be obtained experimentally, perhaps with different ratios for discriminating the case when a surfacebased computation of the Euler characteristic performs better than a volume-based one.

Fig. 1 a
Fig. 1 a The centers of cubes of two interlaced cubic grids (rescaled by factor 2), one composed of points with even and the other composed of points with odd Cartesian coordinates (in different colors).b BCC voxels centered at those points tex, respectively.They are strictly edge-or vertex-adjacent if they are edge-or vertex-adjacent but are not face-or edgeadjacent, respectively.Connectivity relation is the transitive closure of adjacency.A (binary) object O is a finite collection of cubes in the cubic grid.The cubes in O are called black.The cubes in the complement of O are called white.Connected components of O are the maximal connected subsets of O with respect to the chosen adjacency.An object O is called well-composed[26] if it has no critical vertices (incident with exactly two strictly vertex-adjacent black or exactly two strictly vertex-adjacent white cubes) and no critical edges (incident with exactly two strictly edge-adjacent black and exactly two strictly edgeadjacent white cubes).For an object O, the associated cubical complex Q is a cell 3-complex that consists of the cubes in O together with all of their square faces, edges and vertices.An object O is well-composed if and only if |Q| is a 3-manifold.The two definitions of well-composed 3D objects (discrete, through the critical configurations and continuous, through the notion of a manifold) are equivalent[2].

Fig. 3
Fig. 3 Evolution of the lower level set of a torus surface.At each stage, the current lower level set is shaded.The lower level set is homeomorphic a to a disc, b to a cylinder, c to a torus with a hole, d to a torus four, five or six boundary faces.The sets of such vertices are denoted M 3 , M 4 , M 5 and M 6 , and the numbers of such vertices are |M 3 |, |M 4 |, |M 5 | and |M 6 |, respectively.Since each boundary face is a square, and each face angle is π/2, the deficit of each boundary vertex v incident to k boundary faces, 3 ≤ k ≤ 6, is

Fig. 4 Fig. 6
Fig.4 The object composed of one voxel is swept by moving a plane which stops at the vertices

Fig. 7 Fig. 8
Fig. 7 Running times of the boundary-based algorithms SUR, MOR, GB1 and GB2.Times are in milliseconds

Fig. 9
Fig.9 The ratios of running times of GB2 and of GB2+BND over the running time of VOL as functions of r = n ∂ O /n O .The right image is a detail of the left one, focusing on r < 1

Fig. 10
Fig. 10 Running times of the fastest boundary-based algorithm GB2 (red), the algorithm for extracting the boundary surface (orange), and the volume-based algorithm for Euler computation VOL (yellow).Times are in milliseconds

Fig. 13
Fig.13 The two different discretization schemes for hollow spheres BubbleA and BubbleA'.At the lower resolution N = 25 the width is one voxel and 8% of the outer radius.At double resolution N = 50, the width is two voxels and still 8% of the outer radius in the first scheme.In the second scheme, the width remains one voxel and becomes 4% of the outer radius

Fig. 14
Fig. 14 Isolated voxels.The ratios of running times of GB2 and of GB2+BND as a function of the running time of VOL

Fig. 15
Fig.15 Running times for the solid sphere Ball and the three hollow spheres BubbleA,B,C with the first discretization scheme.a The ratios of running times of GB2 and of GB2+BND as a function of the running t. the elevation function f ), then v creates a new component of the lower level set, and the Euler characteristic is increased by 1.-If all neighbors come before v (they are lower than v w.r.t.f ), then v closes off a cavity in the lower level set, and the Euler characteristic is also increased by 1. -If some neighbors come before v and some after v (some are higher and some are lower than v w.r.t.f ), and the ones coming before (after) v are consecutive around v, then v is a regular vertex inducing no change in the Euler characteristic of the lower level set.
-If v has four neighbors in ∂ O, which are alternating before-after-before-after v w.r.t.f , then v either merges two different connected components, or it creates a tunnel by connecting two pieces of the same connected component.The Euler characteristic of the lower level set is decreased by 1.

Table 1
Information about the discretized shapes composing the first test set used in the experiments.N denotes the side of the reference grid, n O the number of black voxels, n ∂ O the number of boundary faces, r = n ∂ O /n O , and r × N

Table 2
Information about the BCC objects representing solid and hollow spheres.N denotes the side of the reference grid, n O the number of black voxels, n ∂ O the number of boundary faces, r = n ∂ O /n O , and (for the first discretization scheme) r × N