Rachid Namane
email: rdnamane@gmail.com

Serge Miguet
email: serge.miguet@univ-lyon2.fr

Fatima Boumghar Oulebsir

F Boumghar
email: fboumghar@usthb.dz

A Fast Voxelization Algorithm For Trilinearly Interpolated Isosurfaces

Keywords: Voxelized Isosurface, Dividing Cubes, Trilinear Interpolant, Incremental Algorithm

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The accurate isosurface representation and its fast extraction have attracted the attention of researchers for several decades. The Marching Cubes (MC) algorithm [START_REF] Lorensen | Marching Cubes: A high resolution 3D surface construction algorithm[END_REF] is arguably the most popular isosurface extraction technique available. In a nutshell, an isosurface is discretized using a triangular mesh by combining all triangular patches produced from the grid cells, with the intention to make it both correct and topologically consistent. An isosurface is correct if it accurately matches the behavior of a known function (or some assumed interpolant) that describes the phenomenon sampled in the data set [START_REF] Newman | A survey of the marching cubes algorithm[END_REF]. If each component of an isosurface is continuous (i.e, there are no holes), then it is topologically consistent. It is possible for an isosurface to have a consistent topology but to not be correct [START_REF] Newman | A survey of the marching cubes algorithm[END_REF].

When dealing with isosurface extraction from 3D data grid, applied algorithms only operate on the discrete data at the vertices of the grid because no surface function is known a priori. Several works were developed to improve the topological correctness and accuracy of the isosurface representation and many authors assume that the ideal isosurface is the one homeomorphic to the surface induced by the trilinear interpolant [START_REF] Nielson | On Marching Cubes[END_REF][START_REF] Sreevalsan-Nair | Using Ray Intersection for Dual Isosurfacing[END_REF][START_REF] Lopes | Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing[END_REF][START_REF] Carr | Subdivision Analysis of the Trilinear Interpolant[END_REF]. This motivates alternatives to the original MC algorithm that produce triangular meshes that are geometrically closer and having the same topology as the ideal surface.

An alternative to the isosurface polygonal representation is the use points and normals to represent small isosurface patches. The use of points in iso-surface visualization was first proposed in the Dividing Cubes (DC) algorithm [START_REF] Cline | Two algorithms for three-dimensional reconstruction of tomograms[END_REF]. Unlike in MC, where the surface within an intersected 3D cell is approximated by a polygonal mesh guided by a local configuration, the DC algorithm approximates the surface by a set of sufficiently dense discrete points. The advantage of using points as compared to polygons is their structural simplicity which is easier to manage, parallelize and no pre-built lookup tables required for local topological guidance. The main drawbacks of point-based approaches, as in the DC algorithm, is the high complexity of the exhaustive scan performed in the extraction process with the large number of points generated.

The main contribution of our work is therefore to develop a fast and efficient incremental point-based method for the discrete approximation and generation of isosurfaces described by trilinear interpolation. In addition, we attempt to get a voxelized isosurfaces which is thin, correct and topologically consistent with respect to the trilinear interpolant without look-up tables. Our incremental generation approach is based on a Bresenhamlike algorithm [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF]. Using normal information, the surface is divided into different zones. For not missing any portion of the surface, the traversal direction in generating the points of each zone is chosen appropriately. The extraction process is repeated recursively until all surface components are correctly retrieved.

The structure of our paper is as follows: after a short presentation of related work in section 2, we introduce in section 3 the necessary terms traditionally used in digital topology. Then in section 4, both the exhaustive scan search and the graph-based search strategies are reviewed before giving a detailed explanation about our incremental approach. The theoretical analysis about the correctness and consistency of the obtained voxelized isosurfaces is discussed in section 5. A further optimization in the running time of the proposed incremental approach is presented in section 6. The efficiency of our approach will be illustrated in section 7 with the obtained experimental results in terms of both computational complexity and surface accuracy.

Related Work

Since the original Marching Cubes algorithm proposed in [START_REF] Lorensen | Marching Cubes: A high resolution 3D surface construction algorithm[END_REF], several variants of that algorithm were developed to resolve topological and geometrical problems related to the triangle mesh representation. Most of these works tend to resolve all possible face and internal ambiguities of every grid cell, and at the same time trying to achieve the best possible polygonal approximation of the trilinear interpolant [START_REF] Nielson | On Marching Cubes[END_REF][START_REF] Sreevalsan-Nair | Using Ray Intersection for Dual Isosurfacing[END_REF][START_REF] Lopes | Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing[END_REF].

The polygonal mesh refinement relies mostly on supplementary points introduced on the faces and in the interior of the cube. Based on this principle, A. Lopes and K. Brodie proposed in [START_REF] Lopes | Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing[END_REF] an extension version of the MC algorithm. A similar work has been proposed by M. G . Nielson in [START_REF] Nielson | On Marching Cubes[END_REF] but with less triangles since extra points are added just in the interior of the cuboid and not on the faces. Both works lead to topological correct surfaces and their supplementary points are determined based on the analytical behavioural analysis of the trilinear interpolant implicit function. In [START_REF] Carr | Subdivision Analysis of the Trilinear Interpolant[END_REF], H.Carr and M. G. Nielson have demonstrated and proved the correctness and completeness of the full possible configurations, classified in three hierarchical levels, already proposed in [START_REF] Nielson | On Marching Cubes[END_REF]. A more recent similar work of [START_REF] Custodio | Practical considerations on Marching Cubes 33 topological correctness[END_REF] presented by L.Custodio and all, in which the authors address some issues related to case ambiguity in the MC 33 algorithm [START_REF] Chernyaev | Marching Cubes 33: construction of topologically correct isosurfaces[END_REF] and proposed a disambiguation procedure to solve the raised issues and obtain topological correct surfaces.

We need to mention that in polygon-based approaches there is always a compromise between a smoother representation of the surface and a minimization of supplementary triangles. This means that a good approximation can be achieved with a finer mesh having a large number of triangles. However, rendering such meshes might lead to triangles whose projected area is less than one pixel, which results in wasting time in their scanline rendering.

Another trend in isosurface visualization investigates the use of points as rendering primitives instead of triangles when multiple facets are projected to a single pixel. Due to their structural simplicity, less computational effort has to be spent per primitive which reduces significantly computations. In addition, no connectivity or topological information is required which makes the geometric processing of highly complex 3Dmodels more flexible. The Dividing Cubes algorithm [START_REF] Cline | Two algorithms for three-dimensional reconstruction of tomograms[END_REF] is among the proposed point-based alternatives to solving problems related to MC-like methods. However, it has some drawbacks related mainly to its time complexity and the excessive number of points used in surface representation. Other uses of point-based primitives for isosurface representation and fast rendering appeared in [START_REF] Co | Iso-splatting: A Point-based Alternative to Isosurface Visualization[END_REF][START_REF] Livnat | Interactive point-based isosurface extraction[END_REF][START_REF] Zhang | Interactive point-based isosurface exploration and high-quality rendering[END_REF][START_REF] Guanfeng | Interactive Exploration of Remote Isosurfaces with Point-Based Non-Photorealistic Rendering[END_REF]. In our paper, we are interested by this family of isosurface point-based representation approaches, in which we propose a much more efficient discretization scheme for fast generation of most classically used point-based primitives. Our fast incremental point generation technique is then based on a Bresenham-like algorithm [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF] and its 3D extension [START_REF] Liu | Three-dimentional extension of Bresenhams algorithm and its application in straight-line interpolation[END_REF], and leads to correct and topologically consistent isosurfaces while avoid using excessive number of points in their representation.

Notions and Problem Definition

Notations and Preliminaries

The set Z n (resp. R n ) of points with integer (resp. real) coordinates is called the discrete (resp. continuous) space and its elements are called discrete (resp. Euclidean) points. Throughout we will assume that n ≤ 3, and we denote by upper-case letters e.g, A,B,C... (resp. calligraphic upper-case letters e.g, A,B,C... ) the points of the discrete space Z n (resp. the Euclidean space R n ). By [a..b], we denote the closed discrete interval in Z between a and b and by [a, b] the continuous closed interval in R between a and b.

For n = 3, a discrete (resp. continuous) domain D (resp. D) is a subset of Z 3 (resp. R 3 ) given by D = [1..X] × [1..Y ] × [1..Z] (resp. D = [1, X] × [1, Y ] × [1, Z]
). D can be seen as a 3D regular grid of size X × Y × Z. Each element in the grid is called a vertex. An edge is a line joining two neighboring vertices, i.e. vertices whose coordinates differ only by one unity, along one of the three coordinate axes. A 3D cell is the axes-aligned unit cube in the grid bounded by eight vertices, and which can be defined by

C = [i, i + 1] × [j, j + 1] × [k, k + 1] with i, j, k ∈ Z 3 .
When the 3D cell is subdivided -along the three coordinate axes -by a × b × c with a, b, c ∈ Z 3 , we then talk about a sub-grid of size a × b × c. A 3D sub-cell is the resulting axes-aligned unit cube in the sub-grid bounded by eight vertices.

We assume that our input data consists of a 3D regular grid dataset. This grid can be represented with a domain D given by D

= [1..X] × [1..Y ] × [1..Z].
The 3D image denoted by I can be represented as a function from D to [0, l], where [0, l] is the depth of the image, and which can be defined for every point P in D by :

I : D → [0, l] P → I(P )

Introduction to Object Digitization

To obtain more details about the following notions, the reader is referred to [START_REF] Kong | Digital topology: introduction and survey[END_REF]. Below we consider a discrete (resp. Euclidean) point as an element of Z 3 (resp. R 3 ), and a discrete (resp. Euclidean) object as a set of discrete (resp. Euclidean) points. Let P be a point in a discrete domain D. The coordinates of a point P are denoted by the tuple (x, y, z). A voxel is a unit cube the center of which is P . This voxel is denoted by V(P ).

The voxels V(P 1 ) and V(P 2 ), having their center points at P 1 (x 1 , y 1 , z 1 ) and P 2 (x 2 , y 2 , z 2 ) respectively, are said to be 26

-neighbors if max(|x 2 -x 1 |, |y 2 - y 1 |, |z 2 -z 1 |) = 1. if V(P 1
) and V(P 2 ) are 26-neighbors, they are also said to be 18-neighbors if

|x 2 -x 1 |+|y 2 - y 1 | + |z 2 -z 1 | ≤ 2.
Similarly, if V(P 1 ) and V(P 2 ) are 18neighbors, they are also said to be 6neighbors if |x 2 -x 1 |+|y 2 -y 1 |+|z 2 -z 1 | = 1. Intuitively, these three neighbouring schemes can be described as a set of voxels which share one face for 6-neighbors, share at least one edge for 18neighbors, and share at least one vertex for 26neighbors. A kpath, with k ∈ {6, 18, 26}, in a discrete object (A) is a sequence of discrete points from (A) such that every two consecutive points are

k -neighbors. (A) is called k -connected if there is a k -path connecting any two points of (A). A k - component is a maximal k -connected subset of (A).
The problem of discretization is to identify discrete points (mesh approach) or what cell (paving approach) belong to the discrete equivalent of an Euclidean object. Several discretization schemes have been proposed. One of the discretization schemes used for digitizing curves -open or closed-is the Grid Intersect Quantization model, denoted GIQ. Whenever an edge -in the working lattice-intersects the curve, then the GIQ model chooses the closest lattice point among the two ones connecting the intersected edge. In fact, this process is similar to the Bresenham algorithm used for line drawing. For closed curves, the scheme Object Boundary Quantization, denoted OBQ (respectively Background Boundary Quantization, denoted BBQ), in which all what is considered as a part of the discrete object are extreme points of each intersected edge and which are located in the interior (respectively exterior or background) of that object. Figure 1.(a), (b), and (c) illustrate these three quantization schemes for 2D curves. See [START_REF] Kiryati | Digital representation schemes for 3D curves[END_REF] for a survey on digitization schemes.

In paving approaches, the discretization scheme -in which the resulting discrete object consists of every cell intersected by the continuous object -is called supercover. This scheme has the advantage of being simple. However its drawback is the possible presence of what are called "bubbles" when a continuous curve passes exactly through a point of the lattice. This is because all the cells sharing this point will be selected (see Figure 1.(d) for an example of a 2D line). Although the probability that a continuous curve passes through a point of the lattice is very low, some cares can be taken to resolve this problem. The resulting discrete scheme is called standard [START_REF] Andres | Discrete linear objects in dimension n: the standard model[END_REF].

Basics of The Trilinear Interpolant and isosurfacing

The isosurfacing problem can be then rewritten as follows: Given the scalar field values T (x, y, z) of the trilinear function T (x, y, z), at the vertices of axis-aligned cube

C = [i, i + 1] × [j, j + 1] × [k, k + 1] with (i, j, k) ∈ Z 3
, determine the set of points P of the isosurface corresponding to isovalue σ where P(x, y, z) = {(x, y, z) ∈ C|T (x, y, z) = σ}

Without loss of generality, we transform the cubic cell C into the unit cube U = [0, 1]×[0, 1]×[0, 1] for convenience, so that trilinear interpolant will be expressed as follows: Given the grid cell of Figure 2 which is composed of eight grid vertices with the scalar values T (0, 0, 0), T (1, 0, 0), ..., T (1, 1, 1). Based on the scalar values at these vertices, the trilinear scalar function T on R 3 is given by:

T (x, y, z) = a+bx+cy+dz +exy+fxz+gyz+hxyz (1)
Where

a = T (0, 0, 0) b = T (1, 0, 0) -T (0, 0, 0) c = T (0, 1, 0) -T (0, 0, 0) d = T (0, 0, 1) -T (0, 0, 0) e = T (0, 0, 0) -T (1, 0, 0) -T (0, 1, 0) + T (1, 1, 0) f = T (0, 0, 0) -T (1, 0, 0) -T (0, 0, 1) + T (1, 0, 1) g = T (0, 0, 0) -T (0, 1, 0) -T (0, 0, 1) + T (0, 1, 1) h = T (1, 0, 0) -T (0, 0, 0) + T (0, 1, 0) -T (1, 1, 0) + T (0, 0, 1) -T (1, 0, 1) -T (0, 1, 1) + T (1, 1, 1)
The approximated continuous surface S (T ,σ) which is consistent to T for a given isovalue σ within a unit cube in R 3 can be defined by the zero set of a geometrically represented implicit function F given by:

F(x, y, z) = T (x, y, z) -σ = 0 (2)
Recall that by definition a cubic cell is being active if it is intersecting the surface. In a discrete domain D, this may occur if σ falls between the minimum and the maximum of all scalar values at the eight vertices of this cubic cell.

Problem Statement

Let I be our 3D input image defined on a discrete grid domain D, and F the implicit trilinear interpolation function defined on a continuous domain D. F coincides with I on the points of D. For all other points P of D, F(P) is given by trilinear interpolating values of I at the eight vertices of the 3D cell containing P. This allows to define F as a continuous surface, denoted S, on all points of D. The main problem of our work consists in development of a rapid voxelization algorithm for extracting an isosurface S, defined by an isovalue σ, from the 3D input dataset. The voxelized isosurface generated by our algorithm, denoted S, consists of a set of discrete points rather than polygons. However, the point-based representation used in our approach should guarantee producing correct and topologically consistent voxelized isosurfaces with respect to S.

Searching Strategies in Point-Based Isosurface Extraction

Several approaches were developed for accelerating and optimizing isosurface extraction, when dealing with polygonalbased techniques. To obtain more details about these approaches, the reader is referred to [START_REF] Newman | A survey of the marching cubes algorithm[END_REF]. We need to mention here that accelerating the search of active cells in the input 3D image is beyond the scope of our work.

In fact, we are looking to optimize and accelerate the searching process of isosurface points in the interior of active cells. Before proceeding to the presentation of our incremental approach, we are going first to present briefly the two main different possible strategies that can be used for isosurface points extraction; the exhaustive and the graph-based scan strategies.

Exhaustive Scan Search

This is the simple and basic strategy applied in the Dividing Cubes. The dividing cubes subdivides each cell intersected by the isosurface into o×p×q sub-cells along the X, Y , and Z directions respectively, such that the sub-cells project onto a single pixel on the image plane.

The scalar field values at the vertices of every sub-cell are obtained off-line by a trilinear interpolation of the values located at the eight vertices of the cell being subdivided. The same process is applied for interpolating the gradient values, that can be used as approximations of the surface normals. Thereafter, the searching of the active sub-cells is performed by applying an exhaustive search; i.e scan the o × p × q new 3D volume along the three axis. The computational time is however costly as the sampling resolution gets finer. The algorithm's complexity it exhibits is of O(n 3 ), when o, p, and q are equal to n. A statistical study on the number of points extracted as a function of the subdivision parameters is given in [START_REF] Boumghar | Complexity of discrete surfaces in the dividing-cubes algorithm[END_REF]. An optimization of the execution time of the "dividing-cubes" algorithm was proposed in [START_REF] Kebaili | Optimal blocs subdivision in the "Dividing-cubes" algorithm. Application to the 3D medical imagery[END_REF] by subdividing the input 3D image into blocks of voxels in order to reduce the time complexity by eliminating the blocks that do not contain active sub-cells.

Graph-based Search

The previous exhaustive strategy can be improved by applying a graph-based searching technique. In this way, the exhaustive scan process is avoided by propagating in the neighbourhood of every active sub-cell. However, in this approach, a starting point must be found for every connected component of the isosurface. Since every connected component must at least intersect one edge of the cell, an initial preprocessing test is made through the edges to identify those points. The isosurface is then tracked through the resulting 3D lattice in the interior of every active cell. The tracking process is performed by an iterative examination of the 6connected neighborhood of every active sub-cell, noting the active ones which have already been traversed. Flag labels are associated with every sub-cell, which is useful to check if it is already traversed. This process continues until all the active sub-cells have been traversed. We note here that the average complexity it exhibits is of O(n 2 ), except for the O(n 3 ) flags initialization. We will show in section4.3.6 an optimization that allows to reach an O(n 2 ) complexity even for this step.

Incremental Search

To optimize more the extraction process, we developed an incremental search method. It is based on the incremental point generation of the discrete implicit surface in the trilinear interpolant. This approach is based on a Bresenham-like algorithm adapted to 3D hyperbolic surface generation. The main advantages of this method are fast surface generation, thinner surface discretization very closer to GIQ, and much less operation spent per primitive. In the following, we give a detailed description of this approach. Before proceeding we need to define our discrete surface and how it is structured into different zones. First, the generated discrete surface S can be considered as the union of the three subsets S x , S y , and S z (S = S x ∪ S y ∪ S z ). Each of these subsets might be itself composed of one or several surface components belonging to the same zone type. The number of connected components in each zone depends on the shape of the trilinear function in the active 3D cell. We define S x as the zone type which consists of the set of points P (x, y, z) such that argmax • Determine the closest x for every (y, z) in zone S x .

• Determine the closest y for every (x, z) in zone S y .

• Determine the closest z for every (x, y) in zone S z .

Discretization Algorithm

First we present the main steps of our voxel traversing algorithm, then more details will be given in the next subsections. As we will show later, the algorithm generates discrete points of the trilinear interpolant by selecting one of the two possible points after a one-step increment along the axis to which the part of surface is most parallel. This axis is changed adaptively as the orientation of the surface is changed based on the normal associated with each point. In the following, every discrete point generated corresponds to the discrete point more closest to the exact surface.

The main steps of the discretization algorithm can be summarized as follows: First, with linear interpolation along edges we locate a starting point on one of the edges of an active cell or on the neighbourhood of a previously extracted surface part. Then, we select one of the three axis which is the most parallel to the surface's normal; the other two axes belong to the plane which is the most parallel to the surface. Once these axes are specified, we incrementally vary the two variables corresponding to the detected plane and calculate the remaining one (Bresenham-like method), followed by updating incrementally the surface's normal. After each increment step, we test whether we are still in the same zone of the surface; i.e check whether the normal is still most parallel to the same axis. The increment process is repeated until all points of that zone are extracted. Finally, we switch to another zone and repeat the whole extraction process until all parts of the surface are fully extracted.

Algorithm Illustration

For more clear presentation of our method, first we show its behaviour in generating points for the three different zone types separately, and then we give its overall process for generating all surface points. Assuming that the surface's implicit function is F(x, y, z) = 0, then we refer by ∆ x , ∆ y , and ∆ z to the step amounts by which the three variables x, y, and z are incremented/decremented respectively.

The algorithm used for the generation of each surface zone is similar, however the role played by control variables in the main loops are interchanged. For example, the generation of one connected component of S x requires two nested loops on y and z coordinates, and generating at each step the corresponding x coordinate.

In the following, we will give our algorithm's illustration for the generation of S x , and it will be the same thing for S y and S z by just changing their axis variables accordingly. In addition, we refer by E x i -in S x -as an indicator of how much a point P i (x i , y i , z i ) -generated at iteration i-is far away from the targeted implicit surface F(x, y, z) = 0. E y i and E z i are defined in a similar way in S y and S z respectively. In the following we show that this indicator value, which is related to scalar values, is proportional to the amount of error related to distances.

Recall that the scalar value of the trilinear interpolant T is equal to σ on the isosurface, less than σ below the isosurface, and greater than σ above the iso-surface. Now, consider F(x, y, z) = T (x, y, z)σ = 0 to be the implicit function of the targeted isosurface. Based on equation [START_REF] Nielson | On Marching Cubes[END_REF], and if we let T (x ′ , y, z) = σ ′ , then

T (x ′ , y, z)-T (x, y, z) = T x (x ′ -x, y, z) = σ-σ ′ = ∆σ = E x Where T x (x, y, z) = bx + exy + fxz + hxyz Finally x ′ -x = E x T x (x ′ -x, y, z)
Similarly, when either variable y or z is changed we obtain y ′y = E y Ty(x,y ′ -y,z) and z ′z = E z Tz(x,y,z ′ -z) respectively. Since the change in the trilinear interpolant's scalar value is proportional to the change in point's position, therefore the absolute value of that error will indicate the proportion of how much a point is far way from the exact surface. We will be based mainly on this error to choose the point more closest to the targeted isosurface during the incremental traversal. Now let's illustrate our algorithm for generating points in S x . Remember that in S x , variables z and y are being changed in two nested loops and x is computed incrementally for each step. The algorithm stops in either direction if the tracking process either switches to another zone or reaches the boundaries of bounding cubic volume. Consider the starting point P 0 , along one of the edges of an active cell, of coordinates (x 0 , y 0 , z 0 ). Based on equation (2), the initial amount of indicating error made at that point from being on the exact surface is given by

E x 0 = a + bx 0 + cy 0 + dz 0 + ex 0 y 0 + fx 0 z 0 + gy 0 z 0 + hx 0 y 0 z 0 -σ (3)
And using partial derivative formula, the normal components at that point are

g x0 = b + ey 0 + fz 0 + hy 0 z 0 (4)
g y0 = c + ex 0 + gz 0 + hx 0 z 0 (5)
g z0 = d + fx 0 + gy 0 + hx 0 y 0 (6)
From P 0 , a one step increment along O y z is performed iteratively in two nested loops. After every increment step, the point more closest to the exact surface along the O x is chosen among the two possible ones ( x i , or x (i+∆x) ), where x i is the xcomponent of the point extracted in the previous iteration and ∆ x is the step increment which is equal to ±1. Now assume that incrementing along O z corresponds to the inner most loop of the algorithm. Therefore, z variable is incremented with ∆ z ; z (i+∆z) equals to .

z i +∆ z , where ∆ z equals ±1 depending on the direction of the traversal along (O z ); upward or downward. The traversal direction can be determined based on the sign of the normal components at the point. If the point extracted in the previous iteration was P (x i , y i , z i ), then in the current iteration we look for the nearest discrete point to the exact surface among the two possible ones which are P (x i , y i , z (i+∆z) ) or P (x (i+∆x) , y i , z (i+∆z) ). ∆ x equals ±1, and it can be determined by just testing the sign of the partial derivative of x with respect to z based on the implicit function theorem. Figure 3 shows the two possible cases that may happen for the z incrementation. The recurrence formulas for the errors of the two possible points at the current step are

E x (i+∆z) = E x i + δE z (7) E x (i+∆z)+∆x = E x (i+∆z) + δE x (8)
where

δE z = K + L × ∆ z (9)
δE x = M + N × ∆ x (10)
Where K, L,M , and N are coefficients that depend only on y. Since y is unchanged, these coefficients are calculated one time as function of y 0 just before running the loop process. Since E x (i+∆z) = F(x i , y 0 , z i+∆z ) and E x (i+∆z) = F(x i+∆x , y 0 , z i+∆z ), the previous coefficients can be calculated as follows:

K = (d × ∆ z ) + (g × y 0 × ∆ z ) (11) L = (f × ∆ z ) + (h × y 0 × ∆ z ) (12) M = (b × ∆ x ) + (e × y 0 × ∆ x ) (13)
N = (f × ∆ x ) + (h × y 0 × ∆ x ) (14)
Either P (x i , y 0 , z i+∆z ) or P (x i+∆x , y 0 , z i+∆z ) is chosen according to argmin(|E x i+∆z |, |E x (i+∆z)+∆x |). Once this point is deduced, the normal components are updated accordingly to check whether we are still in the same zone. This is performed using the following equations

g x(i+1) = g x(i) + α × ∆ z ( 15
)
g z (i+1) = g z (i) if x i+1 = x i g z (i) + α × ∆ x if x i+1 = x i + ∆ x ( 16
)
g y (i+1) = β + γ × ∆ z ( 17
)
where

α = f + h × y 0 ( 18
)
β = unchanged if x i+1 = x i β + e × ∆ x if x i+1 = x i + ∆ x (19) γ = unchanged if x i+1 = x i γ + h × ∆ x if x i+1 = x i + ∆ x ( 20
)
Where α is a constant calculated just before starting the iterative process. Once the normal components are updated and their new values still satisfy the current zone requirements, we pass to the next iteration. Otherwise, P (x i+1 , y i+1 , z i+1 ) generated at the current iteration (i + 1) is added to S x , and the same process is repeated for one step increment along (O y ) in the outer loop. The generation of discrete points in S x is terminated when that zone boundaries are reached along both (O y ) and (O z ).

Algorithm 1 summarizes the iterative process for generating points of S x along O z .

Algorithm 1 points of S x generation along O z

Input: The current active 3D cell C Output: The set of points P (x, y, z) added to the isosurface S Let P 0 (x 0 , y 0 , z 0 ) be a starting point in S x along one of the edges in

C Let i = i 0 while P i ∈ (S x ∩ C) do z i+1 = z i + ∆ z E x (i+∆ z ) = E x i + δE z E x (i+∆ z )+∆ x = E x (i+∆ z ) + δE x if (|E x (i+∆ z )+∆ x | < |E x (i+∆ z ) |) then x i+1 = x i + ∆ x end if Incrementally update normal components if (|g x | > |g y | > |g z |) then if P i+1 (x i+1 , y i+1 , z i+1 ) / ∈ S then S = S ∪ Pi+1(xi+1, yi+1, zi+1) end if end if i = i + 1 end while

Overall Algorithm Behaviour

From the starting points on the edges of the active 3D cell, all the corresponding portions of the discrete surface are extracted. This is considered as the first stage of the extraction algorithm in which its iterative process was explained previously. The second stage of the algorithm is performed in the same manner as the first one, except that the starting points are now the end boundaries of the different zones extracted in the previous stage which were already stored in memory. This is to guarantee that no hole appears in the discrete surface. The process is repeated until all points are extracted; i.e all zones meet to form the full different discrete connected components of the desired discrete isosurface in the interior of the cell.

Figure 4 shows the resulting subvoxelized isosurfaces for two different configuration examples of the trilinear interpolant. For both examples, The scalar values at the eight vertices of the 3D cell are given and the isovalue is chosen to be σ = 90. Points in zones S x , S y , and S z are shown in red, green, and blue colors respectively. Starting points are shown in yellow color. In the first example, and starting from the cell's edge, all surface components are generated during one stage. In the second example, surface components are generated by applying the incremental process for two stages. Notice that not all starting point candidates for this second stage will be used. In fact, if any one of these points is generated from another one, then this later is omitted for being a starting point candidate.

Speed Up Hashing Functions

In equation (1), if any two variables are fixed and the third one is computed then this latter, if it exists, is always unique. In other words, if all points of the surface are projected on any of the three planes (O xy , O x z, or O y z) then no more than one point can be projected on the same location. This is true if we consider the euclidean space. In discrete space, this may not be the case especially for points of the surface that are close to an asymptotic plane.

To avoid processing every point more than one time, and therefore speed up the extraction process, three simple hashing functions were used which are h x (i, j, k) = (j, k), h y (i, j, k) = (i, k), and h z (i, j, k) = (i, j). Their corresponding Hash tables contain a linked list of the 3 rd component for points having the same projection on their corresponding plane. Each point is hashed to only one table according to its greatest normal component. When two or all normal components are equal, an arbitrary choice of any table is taken in order to avoid any collision that may happen. Experimental studies showed us that the length of the linked list in every hash table is almost one element, except for special configurations where this length is two, but not more, for just few surface points.

The advantage of these hashing functions is to use just (3×n 2 ) Booleans instead of (n 3 ) that help to check points traceability during the incremental search. This reduces the complexity of their initialization process from O(n 3 ) to O(n 2 ).

Voxelized Isosurface Correctness and Consistency Theoretical Analysis

Correctness and Consistency

As explained before, topological correctness implies that the discrete approximation should be homeomorphic to the implicit surface of the trilinear interpolant (i.e, it accurately matches its behavior). Consistency if satisfied if no holes are generated in the voxelized approximation.

To verify the accuracy of the extracted voxelized isosurface, we are going to prove that all points generated incrementally belong to the super-cover of the continuous surface S. In other words, we want to prove that the points generated incrementally are a subset of the points generated by the DC algorithm. Based on a recursive approach, we are going to show that if a subvoxel extracted at iteration "i" is being intersected with the trilinear interpolant implicit function, then the subvoxel extracted at the next iteration "i + 1" would be intersected too. Furthermore and without loss of generality, we limit our proof for one slice across a certain surface zone which can be extended to the whole surface by just using appropriate variables across every volume slice within any surface zone.

Before proceeding further, let's first consider that we are in the deepest loop of S x extraction process. Meaning that y variable is kept unchanged, z variable is incremented and the closest x variable, denoted X, is determined. Given F(x, y, z) = 0 as an implicit surface function, this function can be rewritten in the form x = F y (z). Then we refer by E i = X i -F y (z i ) which represents the horizontal error along O x between a generated point P (X i , y, z i ) and the point P (x i , y, z i ) on the exact surface. In addition, the sub-voxel containing P is intersected with the exact surface if and only if

|E i | ≤ 1 2 .
In the following, F ′ y (z) refers to the derivative of F y (z) with respect to z. The curve is most parallel to the O z axis if |F ′ y (z)| ≤ 1 Lemma 1

If |E i | ≤ 1 2 , and |F ′ y (z i + t)| ≤ 1 for all t ∈ [0, 1], then |E i+1 | ≤ 1 2 .
Proof Without loss of generality, assume that P (x i , y, z i ) is the point generated at iteration i and ∆ z = 1. ∆ z is the step amount by which z i is incremented, meaning

z i+1 = z i + ∆ z .
Since F y (z) is continuous and differentiable over [z i , z i+1 ], then applying intermediate values theorem we get

F y (z i+1 ) -F y (z i ) = F y (z i + t) 1 0 (21) = 1 0 F ′ y (z i + t)dt (22)
≤ 1 0 |F ′ y (z i + t)|dt (23)
≤ 1 0 1dt (24)
≤ 1 (25)

F y (z i+1 ) ≤ F y (z i ) + 1 (26)
Now let P (X i , y, z i ) and P (X i+1 , y, z i+1 ) be the discrete points generated at iterations i and i+1 respectively. In other side, let P (x i , y, z i ) and P (x i+1 , y, z i+1 ) be their corresponding points laying on the surface S.

According to the horizontal error defined above, we have the following two statements:

E i = X i -x i = X i -F y (z i ) E i+1 = X i+1 -x i+1 = X i+1 -F y (z i+1 )
Based on these two statements and using equation ( 26) we can write:

E i+1 ≤ X i+1 -(F y (z i ) + 1) (27)
E i+1 ≤ X i + ∆ x -(F y (z i ) + 1) (28) E i+1 ≤ (X i -F y (z i ) + (∆ x -1)) (29)
|E i+1 | ≤ |X i -F y (z i )| + |∆ x -1| (30)
|E i+1 | ≤ |E i | + |∆ x -1| (31)
|E i+1 | ≤ 1 2 + |∆ x -1|; as |E i | ≤ 1 2 (32) Since |∆ x -1| ≤ 1 then, equation (32) implies that |E i+1 | ≤ 3 2 . |E i+1 | ≤ 3
2 means that we have two possible cases for x i+1 which is either x i when kept unchanged or x i + 1 when incremented by one. Since the smallest one which is chosen every iteration, this yields that |E i+1 | ≤ 1 2 and this completes our proof.

Clearly, this proof shows that the set of incrementally extracted points are a subset of the supercover of the surface S. Furthermore, the way in which surface points are tracked along the three coordinate axes makes every resulting discrete surface component free from holes.

Sampling Issues For Topology Preserving

In the full trilinear function T over R 3 , saddle points are approached by contour surfaces from all directions [START_REF] Carr | Subdivision Analysis of the Trilinear Interpolant[END_REF]. If the analysis range is limited to a bounding cuboid in R 3 , there exist two kinds of saddles, face and body. Face saddle occurs on a face having two diagonally opposed vertices with different signs. Body saddles in the interior of the cuboid are points where the gradient of T vanishes. There are exactly 0, 1, or 2 body saddle points with distinct isovalues [START_REF] Nielson | On Marching Cubes[END_REF]. One of the properties of both face and body saddles is that the axis-parallel lines through them have constant isovalue [START_REF] Carr | Subdivision Analysis of the Trilinear Interpolant[END_REF], and isosurfaces only change overall topology at these saddle points. There are five possible cases for the number of face and body saddles that may exist, and the most complicated case occurs when there are six face saddles [START_REF] Bong-Soo | Topology Preserving Tetrahedral Decomposition Applied To Trilinear Interval Volume Tetrahedrization[END_REF].

The voxelized isosurface is correct and topologically consistent if it accurately matches the behaviour of the trilinear interpolant function without missing any point. This means that for topology preserving purposes, our discretization scheme should ensure that both surfaces must be having always the same number of continuous connected components and the same number of tunnels. Which implies that the accuracy of our approximate representation depends on the rate of sampling (i.e, the resolution of subdivision). A sufficient discretization sampling resolution should be specified in order to guarantee that all separated components should not meet, and tunnels should not be filled. With respect to this, and discarding the degenerate cases, a sufficient condition for topology preservation is to make sure that voxels containing saddle points and some of their neighbours are not active. For face saddles, voxels at saddle points with two of their 6connected neighbors but along different axes, should not be active. This is to make sure separating different components on the face. The same thing for body saddle where voxels at the body saddles with their 6-connected neighbours should not be active in order to keep all components that may exist separated and any tunnel not to be filled. A simple checker subroutine can be added to verify the previous conditions and to tell whether the chosen subdivision factors fulfil those conditions or we need to sample further.

Incremental Extraction Time Optimization

For further extraction time reduction in the incremental generation method, our algorithm was further optimized by generating first all zones boundaries then filling rapidly the interior of every zone. The prior knowledge of all zones limits makes the generation of points at the interior of the zones more straightforward. This is because the normal change checking is no more necessary which reduces much more the extraction time.

Zone Boundaries Extraction

The boundaries of all zones are first extracted on the faces of the cube then in its interior. On the faces, limit points are extracted incrementally starting from a point of intersection between the isosurface and a given edge of the cube. Starting from the end points of the discrete curves extracted previously on the faces, limit points separating surface zones in the interior of the cube are extracted in a zig-zag sequence by searching in the 26neighbour of already extracted point belonging to one zone in order to find the next limit point belonging to the adjacent zone. Already traversed branches are skipped.

Zone Interior Filling

Once all zone boundaries are extracted, the remaining parts of the isosurface, i.e zone interior filling, is processed based on a scan-line polygon fill like algorithm. The scan extrema points are the limit points found in the previous stage. The type, either lower or upper, of every extrema point is also deduced during their extraction by checking the sign of the normal components at that point with the type of its adjacent zone. All surface portions of the same zone type are filled by first sorting extrema points, laying on the same row, in the plane to which the surface is most aligned. Then incrementally generating the corresponding hyperbolic arc that falls between its corresponding lower and upper endpoints. The sorting process is not a bottleneck while trying to optimize our algorithm. This due to the limited number of extrema points on each plane's row. In most of the time there are either two or four extrema points. The main advantage of our optimization is to avoid calculating the normal components as described in the proposed incremental search method discussed previously.

Discussion

As seen from equations ( 7) to [START_REF] Bresenham | Algorithm for computer control of a digital plotter[END_REF], our incremental algorithm uses only six arithmetic operations for the computation of each point. This has to be compared to the twentythree arithmetic operations required by a direct implementation of the implicit surface equation [START_REF] Sreevalsan-Nair | Using Ray Intersection for Dual Isosurfacing[END_REF]. The computation of the surface's normal in our method requires one test and f our arithmetic operations when the coordinate of the generated point, along the axis most perpendicular to the surface, doesn't change. However it requires one test and eight arithmetic operations when this coordinate is incremented/decremented Fig. 5 The discrete approximation of all possible topological configurations of the trilinear interpolant using the same notation as in [START_REF] Lopes | Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing[END_REF] for their labelling.

(equations ( 15) to ( 20)). In case the surface's normal is required for further processing, as in surface shading, this later is already obtained by our method with less number of operations as compared to the twentyone additional operations required when computing it directly as in equations ( 4) to [START_REF] Cline | Two algorithms for three-dimensional reconstruction of tomograms[END_REF].

With the optimization performed (zone boundaries extracted first), the need for the normal information is no more required. In this case, the extraction process consists of just three steps. The first step is for extracting zone boundaries as shown above. The second step is for sorting extrema points on the plane to which the surface is most aligned. The last step is for the incremental scan hyperbolic curves filling between lower and upper endpoints (zone interior filling). Since the normal at each point is no more computed here, the algorithm complexity is further reduced as compared to algorithm 1.

Experimental Results and Analysis

In Figure 5, we present the obtained results for all possible configurations of the trilinear interpolant function including tunnel cases and multiple isosurface components per cell. Notice that the same notation as in [START_REF] Lopes | Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurfacing[END_REF] for labelling the cases was adopted here. For visual clarity, the different zones in every connected component are visualized with different colors by involving the free open source C++ DGtal Library [START_REF] Dgtal | Discrete Geometry Tools and Algorithms[END_REF] in our design environment. The three colors correspond to the three surface zones S x , S y and S z . In Figure 7, we show some isosurface examples of real CT medical datasets extracted with an isovalue of 80 and a sampling resolution of 5 × 5 × 5. The DC algorithm was applied on the first row and with our incremental approach on the second row. Isosurfaces of isovalue 80 and sampling resolution of 5×5×5 extracted from the (8bit, 256 × 256 × 256) Aneurism Dataset, downloaded from volvis.org, are shown in Figure 8. As we can see from these figures, there is no visual difference between images extracted by DC and our incremental approach, whereas our representation requires less points and shorter time for surface extraction as will be seen in the next figures and tables.

In Figure 9 the results of optimizing the incremental method are presented. Four different configurations have been chosen and shown in every column. Each row represents one step of the optimization process. The first row represents the results of the zones boundaries extraction stage for the four different configurations. In the second row we show all zone limits with their projection in their corresponding most aligned plane for every zone type. The last row shows the results of the surface filling process.

In Figure 11, a comparative study for the running time of the three isosurfacing methods (DC, graphbased, and our incremental method with normal checking) is given. As examples, we have chosen the three most complex configurations of figure 5 for our tests; i.e configs. 13.4, 13.5.1, and 13.5.2. All computation results were performed on a Mobile DualCore Intel Core i5-460M, 2.533 GHZ provided with an NVIDIA GeForce GT 425M. In this figure, we plotted curves for the running time (T ) for the different methods as a function of the input sampling resolution N along each axis. Therefore, the resulting sub-voxelized resolution in the input active 3D cell is of (N × N × N ). Notice that given running time consists of the time required for points generation with normal checking. The running time required for flag initialization is not considered.

To trace the complexities of these running times, mainly for graph-based and our incremental approach, we plotted in figure 10 the curves for T N 2 versus the input sampling resolution along each axis N ; where T is the computational time already shown in figure 11. We can see clearly that, both the graph-based and our incremental methods are changing quadratically but our incremental approach is faster; i.e ( T N 2 ≈ K 1 ) for graph search and ( T N 2 ≈ K 2 ) for incremental search where K 1 , and K 2 are positive constants with K 2 < K 1 .

Table 1, 2, and 3 show a comparison between the three methods -DC, graph-based, and our incremental approach-in terms of the total number of points generated within the previous three test configurations. We can see that the number of points generated by our incremental algorithm is significantly reduced as compared to the other two methods especially for high sampling resolutions. In fact the resulting discrete surface, produced by both DC and graph-based, corresponds to the super-cover of the continuous trilinear interpolant surface. Therefore the obtained discrete surfaces are 6connected, in this case. In counter part, the discrete surface produced by our incremental approach corresponds to the GID discretization scheme described previously. Therefore it is an 18-connected discrete surface. Since the number of points in a 6-connected surface is higher than the number of points in an 18-connected one, the surface generated by our method is thinner than the one generated by both DC and graph-based approaches. The results show that the number of points In Figure 12, we give some configuration examples in which we show how the discrete topology changes as a function of the discretization sampling resolution. As an example, the first column of this figure represents the case where the topology is not preserved due to the low sampling resolution (separated components are touching and tunnels are filled). For more illustration about what we have discussed in section 5.2, inactive sub-cells covering surface saddle points are shown in yellow. Here we can see that the state, active or inactive, of these sub-cells covering saddle points and their 6connected neighbours have an effect on the correctness of the resulting voxelized surface. In terms of implementation, the topological correctness of the resulting discrete isosurface can be checked with a simple iterative algorithm.

Conclusion and Future Work

In this work, we have presented a point-based isosurface extraction method based on an incremental voxelization approach. As compared to the DC and graph-based approaches, our method produces thinner discrete isosurfaces in a fast and optimized way. Furthermore, our voxelization technique guarantees obtaining correct and topologically consistent isosurfaces. Our method has been successfully tested on all different topological configurations of the trilinear interpolant. The complexity of our incremental algorithm is significantly reduced as compared to either an exhaustive scan or to a graphbased approach.

We have proved that the discrete surface generated by our incremental algorithm is a subset of the surface generated by DC algorithm and which is consistent with the trilinear interpolant. In addition, we have discussed the problem about the choice of the sufficient discretiza-tion sampling resolution which ensures the topological correctness of the voxelized surface around critical points.

The use of the trilinear interpolant for isosurface approximation is mainly due to its simplicity. However this type of interpolant generates G 0 continuous surfaces. For that reason, we plan to smoothly blend the generated discrete patches to make them globally G 1 continuous without changing their topologies. One way to solve this problem is by applying an appropriate re-parametrization of the domain of the scalar field as done in [START_REF] Theisel | Exact Isosurfaces for Marching Cubes[END_REF] followed by a discretization scheme similar to ours. Further direction includes developing a similar strategy for higher order interpolant, such as tricubic spline, in order to obtain G 1 continuous surfaces with low computational time. Finally, a more important computational time reduction can be achieved via a GP U implementation. One possible way to achieve this is by finding starting points laying on surface face curves and then extracting different parts of the surface in parallel.

Acknowledgement

These results were obtained during an eighteen months internship in LIRIS Laboratory, Lyon2 University of France. Therefore we would like to express our gratitude to all members of the LIRIS-M2DisCo team for their valuable feedback and guidance that helped us significantly throughout this work. The work was funded by Algerian Ministry of Higher Education and Research.

Fig. 1

 1 Fig. 1 (a), (b), and (c) are the GIQ, OBQ, and BBQ methods for 2D curve representation respectively. (d) A straight line segment and its super-cover (dashed), which contains a "bubble".

Fig. 2

 2 Fig. 2 Grid cell notation.

  (x,y,z) (|g x |, |g y |, |g z |) = |g x |. We define S y and S z in a similar way, i.e argmax (x,y,z) (|g x |, |g y |, |g z |) = |g y | for S y and argmax (x,y,z) (|g x |, |g y |, |g z |) = |g z | for S z . Based on the previous definition and for efficient surface discretization, the following three possible cases considered during the incremental traversal:

Fig. 3

 3 Fig. 3 The two possible cases that may happen along O x after one step increment of z. (a) case when ∂x ∂z > 0. (b) case when ∂x ∂z < 0. In both cases ∂x ∂z = -g z g x

Fig. 4

 4 Fig. 4 (a) One grid cell synthetic dataset. (b) First stage of extraction. (c) Second stage of extraction.

Figure 6 .

 6 Figure 6.(a) represents a synthetic grid dataset of size 3 × 3 × 3, with random scalar values between 0 and 255, extracted with isovalue value 150 and a sampling resolution of 50 × 50 × 50. Figure 6.(b) is (8bit, 41 × 41 × 41) marschnerlobb dataset extracted with isovalue 192 and sampling resolution of 5 × 5 × 5.Figure

Figure

  Figure 6.(a) represents a synthetic grid dataset of size 3 × 3 × 3, with random scalar values between 0 and 255, extracted with isovalue value 150 and a sampling resolution of 50 × 50 × 50. Figure 6.(b) is (8bit, 41 × 41 × 41) marschnerlobb dataset extracted with isovalue 192 and sampling resolution of 5 × 5 × 5.Figure

Fig. 6

 6 Fig. 6 (a) 3 × 3 × 3 synthetic random scalar values grid dataset. (b) marschnerlobb Dataset. (c) A part of the CT scan of an engine block.

Fig. 7

 7 Fig. 7 Examples of a real CT medical datasets. (a) Extracted with DC algorithm. (b)Extracted with our incremental algorithm.

Fig. 8

 8 Fig. 8 Example of the Aneurism Dataset. (a) Extracted with DC algorithm. (b) Extracted with our incremental algorithm.

Fig. 9

 9 Fig. 9 Zones boundaries extraction and surface filling (a) Config. 1 (b) Config. 9 (c) Config. 13.5.1 (d) Config. 13.5.2

Fig. 10 Fig. 11

 1011 Fig. 10 Time Complexity Analysis. (a) Config.13.4. (b) Config.13.5.1. (c) Config.13.5.2.

Fig. 12

 12 Fig. 12 The effect of the discretization sampling resolution on the correcteness of the discrete surface topology.(a) Config. 6.1.1 (b) Config. 4.1.1 (c) Config. 4.1.2

Table 1

 1 Comparison of the three extraction methods in terms of the number of extracted points for config. 13.5.2

	50 × 50 × 50 100 × 100 × 100 150 × 150 × 150 200 × 200 × 200 250 × 250 × 250 300 × 300 × 300	(6538pts) (26123pts) (58785pts) (104499pts) (163268pts) (235109pts)	(4409pts) (17287pts) (39391pts) (70008pts) (109329pts) (157581pts)
	extracted by our algorithm is almost half the number
	of points produced by DC and graph-based techniques.

Reso.(N × N × N ) DC/Graph.(♯pts) Our Algo.(♯pts)

Table 2

 2 Comparison of the three extraction methods in terms of the number of extracted points for config. 13.5.1

	Reso.(N × N × N ) DC/Graph.(♯pts) Our Algo.(♯pts)
	50 × 50 × 50 100 × 100 × 100 150 × 150 × 150 200 × 200 × 200 250 × 250 × 250 300 × 300 × 300	(6505pts) (26047pts) (58623pts) (104229pts) (162839pts) (234472pts)	(4186pts) (16743pts) (37626pts) (66689pts) (104067pts) (149580pts)

Table 3

 3 Comparaison of the three extraction methods in terms of the number of extracted points for config. 13.4 Reso.(N × N × N ) DC/Graph.(♯pts) Our Algo.(♯pts)

	50 × 50 × 50 100 × 100 × 100 150 × 150 × 150 200 × 200 × 200 250 × 250 × 250 300 × 300 × 300	(7177pts) (28709pts) (64526pts) (114761pts) (179372pts) (258142pts)	(5813pts) (22947pts) (51579pts) (90907pts) (142095pts) (204500pts)

Author Biographies (50 -100 words)