A repartition structure for virtual sculpture of very large multiresolution volume datasets
X. Heurtebise and S. Thon
LSIS: Laboratoire des Sciences de l�Information et des Systèmes
IUT de l�Université de Provence, Rue Raoul Follereau, 13200 Arles � France
and
Contents
1.1.1 Multiresolution representations
2 Extended multiresolution model
4 Our collision detection algorithm
4.1 The complete hierarchical collision detection
4.2 Rough collision detection between two octrees
4.3 The choice of bounding volumes (BVs)
4.4 Fine collision detection between two wavelet enumerations
4.4.1 Simple fine hierarchical collision detection
4.4.2 Enhanced fine hierarchical collision detection
5.1 Steps of the addition/subtraction of matter
5.2 Enhancement of the addition/subtraction of matter
In a virtual sculpture project, we would like to sculpt in realtime 3D objects sampled in volume elements (voxels). The drawback of this kind of representation is that a very huge number of voxels is required to represent large and detailed objects. Consequently, the memory cost will be very large and the user/object interaction will be slowed down. We proposed a multiresolution representation of 3D objects thanks to combination of 3D wavelet transform and octree. Each 3D object is roughly sampled in an octree, where each leaf, called block, containing data is thinly sampled thanks to a 3D wavelet transform. This representation allows us to manage very large volume datasets, by reducing the memory cost and by adapting the processing and display times with a desired level of detail. In order to allow realtime performance during the sculpting process, in this paper, we propose for each block a repartition structure, which is an octree whose each node contains maximal and minimal density for each area of the 3D object. Moreover, we combine this structure with a multiresolution collision detection, to accelerate the sculpting process during the addition and subtraction of matter into the 3D object thanks to a tool, both with the same multiresolution representation.
In this paper, we present a multiresolution model based on 3D wavelets to represent a 3D object as a discrete set of volume elements (voxels). Such a discrete representation is of great use in:
� medical imaging: management of MRI or CTscan data;
� surgical simulations: bone surgery, dental surgery;
� scientific simulations: simulation of heterogeneous systems such as semiconductor device simulation, molecular dynamics, plasma physics, and fluid mechanics;
� virtual sculpture: easy simulation of sculpture operations such as addition or subtraction of material by simply adding or removing voxels.
The major issue of discrete representations of volumes is the huge memory cost of very large 3D objects. We proposed a multiresolution sculpture system based on a combination of 3D Haar wavelets and octree [1]. A major advantage of this model is that sculpted objects can then be used as new tools, because the same model is used for both objects and tools. However, during sculpture operations, when 3D objects and tools are very large, the number of voxels of tools in collision with the ones of objects can become very huge. Consequently, the sculpture operations become slower. In order to accelerate the sculpture operations, we propose in this paper several solutions:
� First, we propose a collision detection algorithm based on an octree in order to refine the collision detection between objects and tools. Both use the same model, based on a combination of 3D wavelets and octree.
� Second, we use a min/max octree, that we called repartition structure, where each node contains maximal and minimal density for each area of the associated object.
� Third, we speed up the collision detection algorithm thanks to the repartition structure during sculpture operations.
In this paper, we would like to sculpt in realtime 3D objects with tools, both using the same model, based on 3D wavelets. In order to accelerate the sculpture operations, we would like a multiresolution collision detection algorithm. We present multiresolution representation in section 1.1.1, in order to represent 3D objects and tools.� Then, we describe in section 1.1.2 existing methods of virtual sculpture. Finally, we present in section 1.1.3 existing collision detection algorithms.
In this paper, we tackle the problem of virtual sculpture of a very large 3D object with a tool, both represented with spatial enumerations. Such a spatial enumeration is a set of volume elements called voxels, obtained by sampling the volume of a 3D object. It can be seen as a 3D image composed of voxels, while a 2D image is a bidimensional array composed of pixels.
To make a spatial enumeration from a 3D object, several methods have already been suggested. The simplest way is a uniform discrete spatial enumeration, by regularly sampling the object into voxels with the same size. However, a major drawback of this representation is the large number of voxels needed to represent very large objects with detailed features (for example, a 3D image in 1024×1024×1024 has more than one billion voxels). This entails three main problems. The first one is the important memory cost to store this uniform spatial enumeration. The second one is that the display of these objects becomes slower. Finally, operations on these objects such as sculpture actions or displacements become less and less interactive.
To further improve the use of spatial enumeration, several methods of multiresolution representations have been proposed. Therefore, processing and display times are adapted with the desired level of detail. Among these methods, there are trees and 3D wavelet decomposition.
An octree can also be seen as a hierarchical representation of 3D object. The maximum level of subdivision of the octree defines the finest level of detail of a multiresolution representation. Boada et al. [2] define a section in an octree that determines the displayed nodes for a given level of detail. This method is extended to a �ntree� by Ferley et al. [3].
The second multirésolution methods use bounding volume trees, used in collision detection. These methods propose to modify properties of the voxels, such as the size with AABB (axisalignedboundingbox) method [4], the orientation with OBB (orientedboundingbox) method [5], or even the shape with sphere tree [6][7]. Thanks to these three methods, the object rendering is optimized because the original object shape can be approached with less voxels than with a simple uniform spatial enumeration, but to the cost of higher computation times.
The third multiresolution method uses wavelet transform. Wavelets are a mathematical tool for representing functions hierarchically. In our case, these functions are discrete 3D functions that define a set of voxels. Muraki [8] shows the use of 3D Haar wavelets [9] to represent a 3D object. Pinnamaneni et al. [10[ build a 3D Haar wavelet decomposition from a sequence of 1D Haar wavelet decomposition in each direction of the 3D voxels grid. Wavelet decomposition allows us to display a 3D object faster according to the level of detail. It also permits to drastically cut down the memory cost, because high compression ratio can be achieved on wavelets coefficients, especially if lossy compression schemes are used. Heurtebise [11] uses same wavelet decomposition into a sequence of 1D wavelet transform in each direction of the 3D voxels grid, but he studied several wavelets, such as Haar wavelet, orthogonal Daubechies wavelet, biorthogonal CohenDaubechiesFeauveau wavelet�
In a voxel representation, several kinds of information can be stored, such as density of matter, color, material, hardness, elasticity� According to the kind of information, the mode of representation may be different.
A first volumetric method has been implemented by Galyean and Hughes [12]. Their model is defined with a uniform discrete spatial enumeration that contains density values. 3D tools modify the discrete potentials describing the object. Basic operations like addition or subtraction, and several tool definitions (heat gun, sand paper or color modifier) are proposed.
Ayasse and Müller [13] perform sculpture operations by the use of CSG (Constructive Solid Geometry). Complex objects are created by successive modifications of the material with a tool according to simple operations such as difference, union or intersection. However, the object and the tool are represented by simple uniform spatial enumerations. Moreover, voxels are limited to binary values (full or empty). The authors propose to reduce the computation time for each sculpture operation by using only the effective voxels according to a given displacement of the tool. However, they do not use a multiresolution representation to improve the display performance.
Raffin et al. [14] perform sculpture operation by moving the matter, into a 3D object, with a tool. Both of them are represented by a set of voxels. The authors proposed a method of diffusion of the matter by decomposing it along each axis x, y and z according the normal of collision point between the tool and the matter to be sculpted. They proposed also an algorithm to distribute the matter to surface neighbors, using a plasticity value of the object matter, in order to perform a more realistic result. The main drawback of this method is the size of the tool in order to have realtime performances.
Dewaele and Cani [15] proposed a deformable model for virtual clay. Sculpture operations are addition, subtraction and deformation of matter through the interaction with rigid tools. Their deformations, both large and small scale, mimic the effects of tools on real clay. The authors proposed two steps. The first one is the processing of the influence of static tools on the matter, by determining the displacement vector for each voxel of the matter. The second one is the displacement of the matter from each voxel of the object to its neighbours.
Angelidis et al. [16] presented sweepers, a new class of space deformations suitable for interactive and intuitive virtual sculpture. When an artist moves a tool, it causes a deformation of the working shape along the path of the tool: the authors used simple path (translation, scaling or rotation). Tools are simply shapes, subsets of 3D space. An advantage is that he artist can use shapes already created as new tools to make more complex shapes. Furthermore, more complex deformations are achieved by using several tools simultaneously in the same region. For representing shapes, Angelidis et al. [16] presented a mesh reﬁnement and decimation algorithm that takes advantage of the definition of deformations.
In the Kizamu project, Perry and Frisken [16] use ADFs (Adaptively sampled Distance Fields) to model and to sculpt the material. A 3D object is sampled adaptively with a 3D grid according to the details of the object. Each grid cell contains a scalar specifying the minimum distance to the object shape. This distance is signed to distinguish the inside from the outside of the shape. Sculpture operations use CSG, but the Euclidean distance field used instead of density raises discontinuity problems and increases the update computation time.
Bærentzen and Christensen [18] propose the LevelSet method to deform the material. This method stores distance fields around the exterior of a 3D object. The single tool is a blob (a sphere) represented by an implicit function, which limits the sculpture capabilities.
To represent an object to be sculpted, Ferley et al. [3] also use distance fields, stored in a �ntree� hierarchical representation where the sampling rate depends on object�s details. The tool is limited to an ellipsoid defined by an implicit function discretized to perform sculpture actions on the object: this limitation is very restrictive if the user wishes to use more complex tool shapes to perform specific sculpture actions. Raffin et al. [19] propose a hierarchical model of virtual sculpture based on an octree [19]. But the tools, defined as voxels sets, remain parallel to the axis, which limits the orientation of the tool and the sculpture operations.
We proposed in [21] a multiresolution sculpture system based on 3D Haar wavelets. A major advantage of this model is that sculpted objects can then be used as new tools, because the same model is used for both objects and tools. Moreover, any orientation of tools [22] is possible, which does not limit the sculpture operations. Furthermore, we proposed in [1] an extension of this model,� that combines octree and wavelet transform (a 3D object is roughly sampled in an octree, where each leaf containing data is thinly sampled thanks to a 3D wavelet transform), to manage very large volume datasets.
In order to deform the matter with a tool, we need to know if tools collide the matter. Two approaches can be used to determine the collision between tools and the matter: static collision detection (only the positions of tools and the matter are known at distinct instant) and dynamic collision detection (the movement of tools and object are considered: the positions of the tools and the matter are known at all instants).
The obvious approaches to collision detection for multiple and/or very detailed objects are very slow. Indeed, if we want to know the potential collision between N objects, checking the collision between every object against every other object is too inefficient to be used, because the complexity is quadratic, in O(N^{2}). Moreover, if we want to know the potential collision between two objects with complex geometry, checking the collision between these objects each face against each other face, is itself quite slow.
Thus, considerable research has been applied to speed up the problem. Kitamura et al. [23], Hubbard [6], and O�Sullivan and Dingliana [24] used hybrid collision detection, with the following phases:
� Broad phase:
o Phase 1: progressive delimitation levels. This phase restricts the subspaces with a potential by using hierarchies of space subdivisions.
o Phase 2: accurate broad levels. In this phase, the approximate test is performed to identify interfering objects using a coarse representation of object shapes (such as bounding volumes).
� Narrow phase:
o Phase 3: progressive refinement levels. Hierarchical approximations are suitable to well determine the objectparts in potential collision.
o Phase 4: exact level. The tests use a tightly representation of object shapes to accurately identify any object parts, selected in the previous phase, that actually cause interference.
In the literature, several data structures are used to solve collision queries. We can classify these data structures in different categories depending on the criterion followed to model the workspace and objects.
For geometric models, collision, proximity and interference queries are computed by using the geometry of the two possible candidate objects. In the narrow phase at the exact level, the queries are formulated by using the geometry of objects. Dobkin and Kirkpatrick [25] proposed a collision detection algorithm between two polyhedrons in O(log^{2} v) where v is the average number of vertex of two polyhedrons. This method requires the convexity of polyhedrons, and only one collision point is given even if more collision points exist. Linear programming algorithms [26] allow checking collision between two convex polytopes if and only if there exists a separation plane between them. Feature based algorithms focus on the relationships between the sets of features (vertices, edges and faces) of the two polytopes. The algorithm of LinCanny [27] constructs the Voronoï Region (set of points closer to a feature than any other) of each feature. Then it computes the distance between the closest features of two polytopes to elucidate whether they collide. This algorithm takes O(f) time, where f is the number of features. The method takes advantage of coherence because closest features will not change significantly between two consecutive frames.� Simplex based algorithms treat the polytopes as the convexhull of a point set. Gilbert et al. [28] presented the GJK algorithm, which detects collisions and gives a measure of interpenetration. Cameron [29] developed the Enhanced GJK algorithm. Volume minimization based algorithms focus on the intersection volume of two objects. Faure et al. [30] proposed an imagebased method of the intersection volume between two polyhedra, using surface rasterization in three orthogonal directions. Moreover, the authors proposed the integration of pressure forces over the pixels of the intersection volume, to compute forces applied to the vertices of the polyhedra. Their method handles deformable and rigid objects without any precomputation, by combining the speed of surfacebased methods [31][32] with the robustness of distancebased methods [33][34].
For space partition trees, space decomposition techniques in a hierarchical ways are considered. In the broad phase at the progressive delimitation level, the queries are often formulated by using the space partition trees. However, these trees are also used to refine collision detection between two objects, during the narrow phase at the progressive refinement levels. Examples of space partition hierarchies include regular grids of voxels [30], octrees [6][23], BSPtrees [35], kDtrees and their extensions [36][37]. The advantage of regular grids of voxels is the computation time in O(1) for the collision detection, but the major drawback is the huge memory cost of this method, for very detailed workspace. The advantage of octrees is the reduced memory cost in relation to the previous method. The computation time is proportional to the number of levels of subdivision of the octrees. For BSPtrees and kDtrees [36], the efficiency depends on the size and/or the number of levels of subdivision of the tree. The main drawback is the build of the tree in O(N^{2}), where N is the number of objects.
For bounding volumes (BV), the objects are enclosed in volumes of simple geometry such as spheres, AABBs, OBBs, kDOPs and convex hulls. BVs are used during the broad phase at accurate broad levels, in order to check easily and quickly collision between two objects. However, BVs are often used in hierarchical volume representations, called bounding volume trees: sphere trees [6][7][37], AABBtrees [4], OBBtrees [5], kDOP trees [38] and convex hull trees [26]. BV trees are used in the narrow phase at progressive refinement levels. Klosowski [38] quantized the computation time for collision detection. They demonstrated that the thinness of BV influences the number of collision tests between each pair of BVs. Moreover, the simplicity of BV influences the efficiency of the intersection test between two BVs.
Many models have already been proposed in order to represent 3D objects as discrete sets of voxels. We have already proposed in previous papers a discrete representation to reduce the memory cost and display times, which become very huge when 3D objects are very large.
In this paper, we propose to solve the problem of computation times, which become very important, when 3D objects are very large and/or very detailed. For that, we propose the use of a min/max octree that we called �repartition structure�. Indeed, it allows us to easily and quickly determine the maximal and minimal densities of an area of a 3D object.� The main contribution of this paper is to use this structure to enhance the collision detection between a tool and the matter, during the sculpture operations: as the tool is moved by the user, its voxels are in different local frame than the sculpture�s one. Such operations of addition and subtraction of material are then simply performed by modifying density values of voxels.
The remainder of the paper is organized as follows. In section 2, we describe our extended multiresolution model based on 3D wavelets and octree, proposed in a previous paper [1]. In section 3, we describe our repartition structure. Then, in section 4, we propose a collision detection algorithm between objects and tools, both using the same multiresolution model. In section 5, we present dynamical and interactive modifications of 3D objects thanks to sculpting tools. Then, we conclude and present future works.
In a virtual sculpture project, we represent the 3D objects to be sculpted as a discrete set of voxels to easily handle subtraction, addition and displacement of matter by tools. Each voxel contains a density value coded in a byte (from 0 for an empty voxel to 255 for a full one). However, a uniform spatial enumeration is expensive in processing and display times.
Thus, we use our extended multiresolution model [1] based on a combination of octree [19] and 3D Haar wavelets [21]. This principle can be extended to other discrete wavelet transform [11], such as orthogonal wavelets, biorthogonal wavelets, interpolated wavelets�
Each object is roughly sampled in an octree where each leaf, called block, containing data is thinly sampled thanks to a 3D wavelet transform (Figure 1).
Fig. 1 Octree where (a) leaves are empty and (b) leaves contain wavelet enumerations
For each block, we use the hierarchical structure proposed by Pinnamaneni et al. [10]. Therefore, for each level of details, the 1D wavelet transformation is applied in x, y and zdirection successively (Figure 2). For each transformation step, we have a block �L� with lowresolution coefficients obtained by a lowpass filter, and a block �H� with detail coefficients obtained by a highpass filter.
�� Original volume����������������������� First run: xdirection���������������� Second� run: ydirection��������������� Third run: zdirection
Fig. 2. 3D Haar wavelet transform
The advantages of this extended model are multiple:
� It takes the advantages of the octree:
o The memory cost of an octree is smaller than the one of uniform enumeration.
o An octree can be seen as a hierarchical representation of a 3D object.
o An octree can be extended or reduced.
� It takes the advantages of wavelet enumeration:
o Each leaf of the octree can be decomposed in wavelet enumeration with a different maximum level of detail.
o The wavelet permits to compress efficiently the 3D objects.
o We can choose different level of detail to display each leaf of the octree to satisfy the conditions of realtime interaction between a user and the 3D object.
We display this discrete object with the Marching Cubes algorithm [39] that provides a smooth surface instead of a set of blocky voxels. During display, they take the advantage of the multiresolution nature of the model given by the 3D wavelets to display the more appropriate level of details according to the situation (distance between the object and the point of view, needed frame rate). However, the limitation of the current display method is that we can see some joint problems if the Marching Cubes algorithm [39] is independently applied on each leaf of the octree, because the links between adjacent blocks are not treated
Moreover, we implemented a data cache [21] that improves performances by storing in memory the useful levels of detail of each block, in order to avoid extracting the level each time we need to use it.
In this paper, we use this same model to represent both objects and tools. So, we will take the advantage of the multiresolution representation to accelerate the processing times, during the collision detection and the sculpting process.
Indeed, consider an area of the 3D object to be sculpted. If the density values at a given level of detail n allow us to deduce that no sculpture operation is possible in this area, we can stop the refinement processing steps of the collision detection at the level of detail n. However, the density values at the level of detail n do not allow us to define the interval of density values (from 0 to 255) at the rougher level of detail n � 1, except if we use the detail coefficients of the wavelet transform. So, in order to define the density repartition of an area of N voxels of the 3D object at a level of detail, we need to know the following data at a finer level of detail: the density repartition of this area (N/8 voxels), and all the detail coefficients (7N/8 detail coefficients) of the wavelet transform in this area. Consequently, in order to define the density repartition of an area of N voxels of the 3D object, at a level of detail, we need to know N coefficients, at a finer level of detail.
In order to reduce the processing time to define the density repartition of an area of a 3D object, we propose the use of a min/max octree, called repartition structure. This structure is an octree that stores the density repartition of this 3D object, at different levels of detail. Each node contains two values: the maximal and minimal densities of an area of the 3D object contained in this node. Figure 3 illustrates this structure in the 2D case thanks to a quadtree: this principle is the same in 3D thanks to an octree. As shown on figure 3, the quadtree is completely developed: the leaves are the areas, with a size of 2×2×2 voxels, and the root is the area, with a size of n×n×n voxels (with n = 2^{p} and p is the number of levels of detail of the 3D object).
Fig. 3. Repartition structure, in the 2D case: a quadtree (on right) contains the minimal (min(Bi))
and maximal (max(Bi)) values of the densities for each zones Bi of the 2D image (on left).
Moreover, for each nonterminal node (i.e. each node different from a leaf) of this octree, the minimal density (respectively the maximal density) is the smallest (respectively the greatest) of the minimal densities (respectively the maximal densities) of all the child nodes. For example of figure 3, the minimal density of the area B_{1} is defined by: min(B_{1}) = min( min(B_{5}), min(B_{6}), min(B_{9}), min(B_{10}) ).
As the octree is completely developed with a structure set by the size of the 3D object, we can represent this octree with a vector V, whose handling is easier, faster and cheaper in memory:
� The first element, at the position 0, is the root of the octree.
� The size of the vector V is N = (n^{3} � 1)/7, where n is the size of the 3D object.
� The size of the repartition structure is S = 2N ´ d, where d is the number of bytes to store a density value (for example, 1 for �char�, 2 for �short int� and 4 for �long int�, on 32bits operating system).
� The element, at the position i < INT( (N � 1)/8 ), has eight children from the position
j_{1} = 8i + 1 to the position j_{8} = 8i + 8. The operator �INT(x)� correspond to the integer part of �x�
� The element, at the position i > 0 has one parent at the position j = INT ( (i � 1)/8 ).
� The level of detail p of the octree contains the elements from the position
a_{p} = 1 + 8 + � + 8p � 1 to the position b_{p} = a_{p} + 8p � 1. Recursively, we have: a_{1} = 2; then if p > 0, b_{p} = 8a_{p}; then, if p > 1, a_{p} = b_{p � 1} + 1.
The octree structure allows us to move easily from a node to its parent or its children into the repartition structure, as each node is set. The build of the repartition structure is easy and starts from the leaves to the root. For each leaf of the octree, we set the minimal and maximal densities thanks to the 3D objects at the finest level of detail. Once all the leaves are filled, we update the repartition structure from the leaves to the root.
During the sculpture process (see section 5), the repartition structure will be updated according to the modified voxels of the 3D object. If one voxel is modified, the minimal or maximal density of the corresponding leaf of the octree will be modified, and we update the repartition structure from this leaf to the root of the octree. However, if several voxels are modified, the update of the octree will be made in one pass, rather than density by density.
In section 1.1.3, we have presented the classification given by O�Sullivan and Dingliana [24]. For the broad phase, there exist several performing algorithms to check quickly collision between N objects and to determine the candidate pairs of objects in collision. Among these algorithms, we will use the �sweep and prune� algorithm [40] that is very efficient.
For the narrow phase, we will use the hierarchical structure of extended multiresolution model of the 3D objects, in order to refine progressively the collision detection from the entire object (roughest level of detail) to the different voxels (finest level of detail).
For the narrow phase, the complete hierarchical collision detection based on our extended multiresolution model is performed in two steps, for each candidate pairs of objects in collision:
� a first step for rough collision detection between the octrees of two objects in collision: we obtain a list of candidate pairs of leaves, containing wavelet enumeration, in collision;
� a second step for fine collision detection between the wavelet enumerations of two objects in collision.
During the narrow phase, we refine the collision detection of each candidate pairs of objects in potential collision thanks to �sweep and prune� algorithm, by using the collision detection based on octrees [6][23].
For that, we refine the collision detection from the root of each object to their leaves. At each step, when two nodes collide each other, we refine the collision detection with the children of the smallest node. We stop the refinement when two nodes are leaves: if two nodes contain data, we refine the collision detection with the algorithms described in section 4.4.
The main drawback of this method is the necessity to go to the leaves in order to check collision between two wavelet enumerations for each object. Indeed, if two objects are very detailed, the computation time will be very high, and the collision detection can be noninteractive. In order to accelerate the collision detection, we propose to use a bounding volume of each node.
The choice of BVs for the detection collision is very crucial, because the computation time depends on this choice.
Each 3D object is sampled in a set of voxels that are cubes. If we check collision between two 3D objects, sampled in two sets of cubes, the simplest method of collision detection uses the OBB to represent each cube, in their respective orientation. However, the computation time of collision detection between two OBBs is greater than two AABBs or two spheres: 200 operations needed to perform the interference test between two OBBs versus 6 for AABBs and 8 for spheres. Therefore, we prefer to use AABB or sphere to represent each voxel in order to reduce the computation time of collision detection. If we use AABBs to represent each voxel of 3D objects, it is necessary to compute again the AABBs for each voxel when 3D object rotates. In order to allow the rotation of 3D objects, we prefer using spheres to represent each voxel during the collision detection.
However, if we want to reduce the most possible the computation time and to allow any orientation of 3D objects, we use the smallest AABB of the bounding sphere for each voxel during collision detection, in despite of the tightness of the detection. Indeed, for any orientation of 3D objects, the AABB of the bounding sphere for each voxel does not change. Moreover, the computation time of collision detection between two spheres is greater than two AABBs of two spheres: eight operations needed to perform the interference test between two spheres versus six for AABBs.
In order to refine the collision detection between two blocks during the narrow phase, we use two methods: the first one is based on octree method, and the second one is based on our repartition structure.
In order to refine the collision detection between two wavelet enumerations during the narrow phase, we propose a very easy method, based on the collision detection thanks to octree method, with following conditions:
� Each node of the octree becomes a voxel at a given level of detail of 3D wavelet transform.
� The eight children of a node of the octree become the eight underlying voxels at a more detailed level of detail of 3D wavelet transform.
� The root of the octree becomes the voxel at the roughest level of detail of 3D wavelet transform.
� A voxel whose density value is lower than threshold T, is considered as empty, else it is considered as full or partially full. It is necessary that threshold T is very small in relation to the maximal density.
The main drawback of this algorithm is the necessity to go to the finest level of detail of each object in order to check collision between two voxels for each object. Indeed, if two objects are very detailed, the computation time will be very high, and the collision detection can be noninteractive. In order to accelerate the collision detection, we propose to use not only a bounding volume of each voxel (described in previous section), but also the repartition structure (described in section 3).
We propose an enhanced hierarchical collision detection algorithm, based on our repartition structure.� Consider Ni the node of the repartition structure of the 3D object i, for each level of detail LOD_{i}. If the BVs of N_{1} and N_{2} have no intersection, there is no collision: refinement of collision detection is useless.
Consider the case of the BVs of N_{1} and N_{2} are in collision. If repartition structures give for N_{1} or/and N_{2} maximal densities that do not exceed the threshold T, there is no collision. On the other hand, if repartition structures give for N_{1} and N_{2} minimal densities that exceed the threshold T, there is necessarily collision for the voxels of each object at the finest level of detail, contained in the intersection of the BVs of N_{1} and N_{2}. In these two cases, refinement of collision detection is useless. Else, we refine the smallest voxel in size between N_{1} and N_{2}, if possible. However, if no refinement is possible for both N_{1} and N_{2}, so all the voxels of two 3D objects, at the finest level of detail, are in the intersection of the BVs of N_{1} and N_{2}.
�With this method, it is no longer useful to know the 3D object for each level of detail, but only the repartition structure associated to the 3D object, which is cheaper in memory than the 3D object and its levels of detail, or even the 3D wavelet transform. Indeed, the repartition structure has a memory cost of (n^{3} � 1)/7 where n×n×n is the size of the 3D object (with n = 2^{p} and p is the number of levels of detail of the 3D object). Meanwhile, all the levels of detail of the 3D object have a memory cost of (8n^{3} � 1)/7 and the 3D wavelet transform has a memory cost of n^{3}.
Our model can easily handle sculpture operations by simply adding or removing matter into a 3D object thanks to a tool, i.e. by modifying the density values of each voxel of the object to be sculpted. A major advantage of our method is that the tool used for virtual sculpture has the same representation than the matter. So, the user can create his or her own tools to sculpt another 3D object. The tools can take any orientation and any position in relation to the object [21].
During the sculpture operations, such as addition and subtraction of matter thanks to a tool, a first phase allows us to determine which voxels of the tool and the matter to be sculpted are in collision: this phase is called �collision detection phase�. The following operations are performed:
� In a first step, we perform the collision detection, by using the �sweep and prune� algorithm to determine which tools Ti collides the matter M to be sculpted.
� In a second step, we refine the collision detection between each tool Ti� and the matter M. We obtain a list L_{BLOCKS}[Ti,M] containing all the pairs of wavelet enumerations, called blocks, for tool Ti and matter M, in collision. During the sculpting process, there are two cases of sculpting operations, described in [1]:
o When we add matter thanks to the tool, we determine which leaves of the octree of the object intersect the tool, and we create new leaves with 3D Haar wavelet if necessary.
o When we subtract matter thanks to the tool, the matter contained in a leaf of the octree can be entirely deleted. Therefore, we can delete the wavelet data of this leaf, as well as the leaf itself and potentially many empty parent nodes.
Function FINE_COLL (voxels: V_{Ti}, V_{M}; LOD: lod_{Ti}, lod_{M}) �If V_{Ti} collides V_{M} �Then �� Variable V�_{M} �� If (Op = �subtraction of matter�) �� Then V�_{M} = V_{M} �� Else V�_{M} = Max density � V_{M} �� If (V_{Ti} ³ T and V�_{M} ³ T) then ����� If (lod_{Ti} = _{}0 and lod_{M} = _{}0) then �������� Add V_{Ti} to the list L_{VOXELS}[Ti], and V_{M} to the list L_{VOXELS}[M] ����� Else if (lod_{Ti} = _{}0 or (lod_{M} ≠ _{}0� and size(V_{M}) > size(V_{Ti}))) then �������� For each V_{M}[k] Ì V_{M}, with LOD lod_{M}1����������� Do FINE_COLL (V_{Ti}, V_{M}[k], lod_{Ti}, lod_{M}�1) ����� Else �� ������For each V_{Ti}[k] Ì V_{Ti}, with LOD lod_{Ti}�1������ Do FINE_COLL (V_{Ti}[k], V_{M}, lod_{Ti}�1, lod_{M}) 
Alg. 1. Simple collision detection algorithm between two blocks.
� In a third step, we refine the collision detection between each tool Ti� and the matter M by using the algorithm 1. We obtain two lists L_{VOXELS}[Ti] and L_{VOXELS}[M] of voxels for tool Ti and matter M, respectively. For this step, when we subtract the matter with a tool, we check collision between each voxel V_{Ti} of the tool Ti with the dual V'_{M} of the voxel V_{M} of the matter M, whose value is equal to the maximal density minus the current density.
Then, we perform a �sculpting phase� whose operations are the following:
� In a first step, we apply the method of discrete rotation, presented by Heurtebise and Thon [21], only for the voxels in the intersection area, which is the AABB of the list L_{VOXELS}[Ti] of voxels for tool Ti. We obtain a new list L_{VOXELS+ROTATION}[Ti] of voxels of the tool after the discrete rotation.
� In a second step, for each voxel V_{M} of the matter M, with V_{M} included in L_{M}, we find which voxels V_{Ti} of the tool Ti, with V_{Ti} included in L_{VOXELS+ROTATION}[Ti], intersect it. Then we compute the new value of V_{M} according to the selected sculpting mode and the filling percentage D of the voxel V_{M} by the voxels V_{Ti} of the tool Ti. This voxel value is clamped to the maximal value (respectively the minimal value) during the operation of addition (respectively subtraction) of matter, because we have density values in a given interval.
o In the �addition of matter� mode, the filling percentage D is given to the value of the voxel V_{M} of the matter M, only if the filling percentage D is greater than the value of the voxel V_{M}.
o In the �subtraction of matter� mode, the empty percentage D', equal to the difference between the maximal density and D, is given to the value of the voxel V_{M} of the matter M, only if the empty percentage D' is smaller than the value of the voxel V_{M}.
� In a third step, we modify the wavelet data according to this new voxel value. In order to accelerate this step, we can update the wavelet data after the modification of a set of voxels V_{M} of the matter, instead of voxel by voxel.
� In a fourth step, we update the repartition structure according to the modified voxels.
� In a last step, for each level of detail, we use the Marching Cubes algorithm to rebuild the triangulated surface, only for the modified parts of the 3D object, corresponding to the matter M, to improve the computation time.�����
Fig. 4. Addition and subtraction of matter to an object with a tool.
Function ENH_FINE_COLL (nodes: N_{tool}, N_{obj}; mode: Op) �If N_{tool} collides N_{obj} �Then �� Variable N_{max}, N_{min} �� If (Op = �subtraction of matter�) �� Then N_{max} = N_{obj,max}�� and�� N_{min} = N_{obj,min} �� Else N_{max} = Max density � N_{obj,min}�� and ��N_{min} = Max density � N_{obj,max} �� If (N_{tool,max} ³ T) and (N_{max} ³ T) then ����� If ((N_{tool} is leaf or N_{tool,min} ³ T) and (N_{obj} is leaf or N_{min} ³ T)) then �������� Add to the list C_{VOXELS}[tool,obj], ����� Else if (N_{obj} is leaf or N_{min} ³ T) then �������� For each child N_{tool}[i] of node N_{tool}���������� Do ENH_ FINE_COLL (N_{obj}, N_{tool}[i], Op) ����� Else if (N_{tool} is leaf or N_{tool,min} ³ T) then �������� For each child N_{obj}[i] of node N_{object}��������� Do ENH_ FINE_COLL (N_{obj}[i], N_{tool}, Op) ����� Else if (size(N_{tool}) > size(N_{obj})) then �������� For each child N_{tool}[i] of node N_{tool}���������� Do ENH_ FINE_COLL (N_{obj}, N_{tool}[i], Op) ����� Else �������� For each child N_{obj}[i] of node N_{obj}�������������� Do ENH_ FINE_COLL (N_{obj}[i], N_{tool}, Op) 
Alg. 2. Hierarchical collision detection algorithm between two blocks.
An enhancement of this sculpting method is the use of the repartition structure, in the algorithm 2, to apply sculpting operations not at the finest level of detail, but at a rougher level of detail, in order to reduce the computation times, if necessary. The interest is to permit to modify interactively the matter with a tool. In the second step of the collision detection phase, we stored in the list C_{VOXELS}[Ti,M] the couples of voxels (V_{Ti},V_{M}) or (V_{Ti},V'_{M}) of the tool Ti and the matter M, respectively, in collision at the finest level of detail of two objects.
The results given in this paper have been obtained on a PC with an Intel Core 2 Duo 3 GHz with 4 GB memory, a NVidia GeForce 9800 GTX+ with 1 GB video memory, and Windows Seven 32 bits.
To determine the performance of each algorithm, we set the size of an object in 512×512×512 voxels, where each block, containing wavelet enumerations, are in 32×32×32 voxels. Then for several sizes of tools, we perform two kinds of tests: addition and subtraction of matter.
For that, we create a filled (respectively empty) cube and a tool with a given shape (ex: sphere), for subtraction (respectively addition) of matter. Then, we move the tool inside the 3D object with a given step and we compute the average times per frame of sculpture operations. For each algorithm, the sculpture times are given in table 1.
In table 1, we compare our algorithms with a virtual sculpture method without collision detection algorithm to refine the sculpture operations. In order to compare the collision detection algorithms, the sculpture time excludes the display times, which depends on the level of detail for display.
Alg. 
Hierarchical 
Addition of matter 
Subtraction of matter 

32x32x32 
64x64x64 
128x128x128 
32x32x32 
64x64x64 
128x128x128 

#0 
without 
49.8 
233 
728 
49.5 
232 
729 
#1 
simple 
83.9 
433 
1353 
83.6 
435 
1359 
#2 
enhanced 
51.9 
194 
637 
83.6 
435 
1359 
Tab. 1. Average sculpture times (in milliseconds) per frame, for addition and subtraction of matter.
The table 1 shows that the sculpture times depend on the size of the tool and the collision detection algorithm. Indeed, if the size of the tool increases, the average number of voxels of the tool in collision with the matter increases, consequently the sculpture times increase.
Furthermore, the sculpture times, for the algorithm 1, are higher than for the algorithm without collision detection. Indeed, the algorithm 1 uses a collision detection based on an octree method, where the detection is made on all the levels of detail of the object and the tool. So, the collision times become so high that the sculpture times, for the algorithm 1, is higher than in the case where no collision detection is used.
However, if we use the algorithm 2, based on our repartition structure, the sculpture times become smaller than the ones with the algorithm 1 or the algorithm without collision detection. Indeed, it is not needed to detect the collision between a tool and the matter at the finest level of detail, because the algorithm 1 stops the refinement before the finest level of detail when minimal density is greater than a given threshold. Furthermore, as the repartition structure, based on an octree, is completely developed, the processing time for each displacement into the repartition structure is optimized.
On the other hand, the repartition structure associated to each block of a 3D object is cheaper in memory than the block of 3D object and its levels of detail, or even the 3D wavelet transform. Indeed, the repartition structure has a memory cost of (n^{3} � 1)/7 where n×n×n is the size of the 3D object (with n = 2^{p} and p is the number of levels of detail of the 3D object). All the levels of detail of the block of 3D object have a memory cost of (8n^{3} � 1)/7 and the 3D wavelet transform has a memory cost of n^{3}.
Consequently, the repartition structure allows us to reduce the memory cost and the sculpture times during the sculpture process, thanks to the use of detection collision algorithm based on the repartition structure.
We have presented in this paper a multiresolution model for virtual sculpture of 3D objects with tools, using 3D wavelet. The major drawbacks of our multiresolution model [1] are the problem of computation times and the limitation of the memory size, which become very important, when 3D objects are very large and/or very detailed. In order to solve these problems, we proposed in this paper the use of a min/max octree that allows us to easily and quickly know the maximal and minimal densities of an area of a 3D object. Moreover, during the sculpture operation, we do not need to store all the levels of detail of our multiresolution model thanks to this structure. Then, our main contribution was the use of this structure to enhance the collision detection between a tool and the matter to be sculpted, during the sculpture operations. Finally, we simply perform operations of addition and subtraction of material by modifying density values of voxels. To verify the applicability of our sculpting system, we have conducted many sculpting sessions, which have resulted in numerous interesting sculptures. Sculpture examples are shown on figure 5: Dental Scan [41] in 512´512´167 voxels, and microCT scan of a human lumbar vertebral body L3 [42] in 244´512´512 voxels. The size of the blocks is 32´32´32 voxels, as in [1].




Teeth before 
Teeth after 
Vetebra before 
Vertebra after 
Fig. 5. Teeth (on left), and human lumbar vertebral body L3 (on right), both sculpted with several tools (sphere and cylinder, in different size between 8´8´8 and 64´64´64), thanks to the two kinds of sculpture operation: addition and subtraction of matter. The collision detection algorithm used for these sculptures is the algorithm 2. The framerate during the sculpture operations and the display is upper than 5 FPS.
Many improvements of our sculpture system are possible, by investigating open issues such as dynamical collision detection, other sculpture operations, realistic deformation, memory cost or computation time.
Our collision detection algorithm does not allow us to know if a tool collides the matter between two distinct instants. We plan to extend our collision detection algorithm to the spatiotemporal collision detection, in order to have realistic sculpting operation.
We will also investigate other sculpture operations, such as the displacement of the matter, and we plan to take more advantage of the levels of detail of the 3D wavelet enumerations, in order to accelerate these sculpture operations. In order to obtain a realistic deformation of the matter by a tool, we plan to add other information into our model, such as viscosity or hardness.
[1] X. Heurtebise and S. Thon, Multiresolution representation and deformation of very large volume datasets based on Haar wavelets, In GMAI �08 Proceedings of the 2008 3^{rd} International Conference on Geometric Modeling and Imaging, IEEE Computer Society, pp.34�40, 2008.
[2] I. Boada, I. Navazo, and R. Scopigno. Multiresolution Volume Visualization with a TextureBased Octree, In The Visual Computer, Vol. 17(3), pp.185�197, 2001.
[8] S. Muraki, Approximation and rendering of volume data using wavelet transforms, In Proceedings of IEEE Visualization'92, pp.21�28, 1992.
[9] A. Haar, Zur Theorie der orthogonalen Funktionensysteme, Ph.D. thesis, Gottingen, Germany, 1909.
[10] P. Pinnamaneni, S. Saladi, and J. Meyer, 3D Haar Wavelet Transformation and TextureBased 3D Reconstruction of Biomedical Data Sets, In VIIP�01Proceedings of the International Association of Science and Technology for Development, Marbella, Spain, ACTA Press, pp.389�394, 2001.
[12] T. A. Galyean, and J. F. Hughes, Sculpting: an Interactive Volumetric Modeling Technique, ACM SIGGRAPH Computer Graphics, Vol. 25(4), pp.267�274, 1991.
[14] R. Raffin, G. Thibault, and G. Gesquière, Simple and efficient tools for virsculpt, In GRAPP�06 Proceedings of the first international conference on computer graphics theory and applications, Lisboa, Portugal, pp.436�440, 2006.
[16] A. Angelidis, G. Wyvill, and M.P. Cani Sweepers: Swept UserDefined Tools for Modeling by Deformation, In SMI�04 proceedings of 2004 International Conference on Shape Modeling and Applications, IEEE, Genoa, Italy, pp.63�73, 2004.
[17] R. N. Perry, S. F. Frisken, Kizamu: A System for Sculpting Digital Characters, In Proc. of ACM SIGGRAPH 2001, pp.47�56, 2001.
[19] J. A. Bærentzen, Octreebased volume sculpting, In Proc. of IEEE Visualization 1998, IEEE CS Press, pp.9�12, 1998.
[20] R. Raffin, G. Gesquière, E. Rémy, and S. Thon, VirSculpt: A Virtual Sculpting Environment, In GraphiCon�04 proceedings, Moscow, pp.184�187, 2004.
[21] X. Heurtebise, S. Thon, and G. Gesquière, Multiresolution representation and deformation of waveletbased 3D objects, In Short Communication proceedings of WSCG�06, Plzen, Czech Republic, 2006.
[23] Y. Kitamura, A. Takemura, and F. Kishino, Efficient collision detection among objects in arbitrary motion using multiple shape representations,. In Proc. of the 12th IAPR International Conference, Jerusalem, Israel, pp.390�396, 1994.
[24] C. O'Sullivan, and J. Dingliana, Realtime collision detection and response using spheretrees, In Proc. of Spring Conference on Computer Graphics, pp. 8392, 1999.
[25] D. P. Dobkin, and D. G. Kirkpatrick, Determining the Separation of Preprocessed Polyhedra  A Unified Approach. In ICALP�90 proceedings of the 17th International Colloquium on Automata, Languages and Programming, London, UK, pp.400�413, 1990.
[26] R. Seidel, Linear programming and convex hulls made easy, In Proc. of the 6th annual ACM symposium on Computational Geometry, New York, USA, pp.211�215, 1990.
[27] M.C. Lin, and J.F. Canny, A fast algorithm for incremental distance calculation, In Proc. of IEEE Int. Conf. on Robotics and Automation, Sacramento, CA, USA, pp.1008�1014, 1991.
[30] F. Faure, S. Barbier, J. Allard, and F. Falipou, Imagebased Collision Detection and Response between Arbitrary Volumetric Objects, In SCA'08 proceedings of the 2008 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, Dublin, Irland, pp.155�162, 2008.
[31] G. Hirota, S. Fisher, A. State, C. Lee, and H. Fuchs, An implicit finite element method for elastic solids in contact, In Proc. of Computer Animation 2001, Seoul, pp.136�146, 2001.
[36] M. Held, J. T. Klosowski, and J. S. B. Mitchell, Speed comparison of generalized bounding box hierarchies, Technical report, Department of Applied Math, State Univ. of New York at Stony Brook, 1995.
[37] S. Redon, Algorithms for interactive dynamics simulation of rigid bodies, Ph.D. thesis, Evry, France, 2002.
[39] W. E. Lorensen, and H. E. Cline, Marching Cubes: A High Resolution 3D Surface Construction Algorithm, In Computer Graphics, Vol. 21(4), pp.163�169, 1987.
[40] J. D. Cohen, M. C. Lin, D. Manocha, and M. K. Ponamgi, ICOLLIDE: an interactive and exact collision detection system for largescale environments, in SI3D�95 proceedings of the 1995 symposium on Interactive 3D graphics, pp.189�196, 1995.
[41] DICOM sample image sets (Osirix), http://pubimage.hcuge.ch:8080/
[42] G. Beller, M. Burkhart, D. Felsenberg, W. Gowin, H.C. Hege, B. Koller, S. Prohaska, P. I. Saparin, and J. S. Thomsen, Vertebral Body Data Set ESA2999L3, http://bone3d.zib.de/data/2005/ESA2999L3/