Venceslas Biri
email: biri@univ-mlv.fr

Gaspard Monge

Real Time Single Scattering Effects

Keywords: Single scattering, Real Time, Hardware rendering

Rendering mist, haze or fog remains a challenge in current computer graphics since it is intrinsically a 3D problem. While the attenuation caused by fog is easy to implement, single scattering effects such as glows and shafts of light, that increase considerably the realism, are harder to render in realtime. This paper addresses the rendering of such effects along with volumetric shadows induced by shadow casters in the participating media. Whereas techniques related to shadow maps have been explored when rendering with single scattering, this paper proposes a real-time algorithm using the philosophy of shadow volumes, including volumetric shadows. With a spatial coherence method, simple shaders and an intensive use of the stencil buffer, we render the shadow planes in a back to front order to obtain the correct volumetric shadows. Therefore our method is easy to integrate in a graphics engine using the shadow volume technique since it requires only a little additional texture memory and is implemented with simple shaders. Realistic images can be produced in real-time for usual graphic scenes and at a high level framerate for complex scenes, allowing changes in the properties of participating medium, animations of objects and even light sources movements.

I. INTRODUCTION

If light scattering has been intensively explored, from Blinn in 1982 [START_REF] Blinn | Light reflection functions for simulation of clouds and dusty surfaces[END_REF] to Sun et al. [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF], it is mainly because it greatly enhances the realism of virtual scenes and can greatly improve the graphic quality of computer games. Indeed, the light scattering occurs everywhere in a scene and therefore is intrinsically a 3D phenomenon. For example, look at the differences between the Figure 1.a rendered classically and the Figure 1.b where the contributions of a participating medium have been added. In the first picture, it is not clear where the lights are while they can be roughly located in the second. Nevertheless, the Figure 1.c shows that, in some cases, rendering only the participating medium is not enough. In order to increase realism, we must render the volumetric shadows. Volumetric shadows are the shadows that cast objects in the participating medium itself. They greatly contribute on the presence in a virtual world. Here, the representation of shadow volumes is necessary to obtain a realistic image and understand clearly that a light stands inside the amphora.

A simple and common way to model light scattering is to handle the attenuation due to participating medium and consider the multiple scattering of light as homogeneous and constant over the scene. This is the OpenGL fog model which is popular since it only needs a fog color and is handled automatically by graphic cards. Figure 1.a is a perfect example of this technique. Despite its great speed, this model poorly represents all the effects induced by light scattering. On the other hand, methods seeking for great realism have investigated the computation of the multiple scattering of light through the medium. But due to this goal, these methods are very slow since they need Monte Carlo or finite element techniques.

In this paper we investigate the complete single scattering illumination model. In this local illumination approach, only the first scattering of light in the medium is taken into account. The remaining scattering of light is considered as homogeneous and constant all over the scene, like in the OpenGL fog model. The originality of the model is the introduction of the indirect single scattering. Indeed, we will call direct single scattering the effect of the first scattering of light along view rays. The indirect single scattering is the same effect but along illumination rays of any point of scene. Our contributions in this paper are :

• Define an analytic and comprehensive formulation of light scattering along view rays and illumination rays.

Based on a angular formulation of the radiative transfer equation, we present a way to use precomputed 2D table to compute directly these contributions. • Integration of the volumetric shadows. We build a method based on the shadow volume technique and using spatial coherence strategy, allowing the rendering of volumetric shadows in the scene, especially discernable around light sources. • Hardware implementation. Except for the determination of object's silhouette, all the work is done in hardware. We store our precomputed 2D tables in textures and use simple shaders to render the illumination of objects and participating media. The rendering of volumetric shadows also involves an intensive use of the stencil buffer.

Our goal is to design an algorithm that can render accurately participating media, including effects like light beams in foggy or smoky scenes. Our method is not restricted to isotropic participating media which can be lit by one or several, static or moving, point light sources since no precomputation is done involving either lights or camera. Our technique produces high resolution images and takes into account volumetric shadows, cast by occluders contained in the media. With very few texture memory cost, but using intensively graphics hardware, our method can render images at a high frame rate and is real-time for classical graphics scene. Our method is also easy to implement in traditional graphics engines since it follows the same strategy than the shadow volume algorithm, and use only shaders and textures. Therefore, it is straightforward with our method to obtain animations where objects or even light sources can move.

II. RELATED WORK

The representation of participating media has been a real challenge for years. We can easily divide all these studies between the single and the multiple scattering methods. The last ones try to compute all light reflections and inter-reflections inside the medium, whatever the number of these reflections. Despite their realism, they suffer from excessive computation time due to the complexity of light exchanges occurring in these cases. Therefore these approachs are not suitable for our goal and we will focus on methods considering only the single scattering case.

These techniques [START_REF] Hoffman | Rendering Outdoor Light Scattering in Real Time[END_REF], [START_REF] Dobashi | Interactive Rendering Method for Displaying Shafts of Light[END_REF], [START_REF] Dobashi | Interactive Rendering of Atmospheric Scattering Effects Using Graphics Hardware[END_REF]- [START_REF] Max | Efficient Light Propagation for Multiple Anisotropic Volume Scattering[END_REF] approximate multiple reflections of light as a constant ambient term and consider only the first scattering of light ray in the direction of the camera. This assumption allows a direct rendering of the illumination of the medium which is more suitable for interactive rendering. Visualization is often done by ray tracing or ray marching. View rays are followed to gather the participating media contributions. Unfortunately, these methods [START_REF] Foster | Modeling the motion of a Hot, Turbulent Gas[END_REF], [START_REF] Jensen | Efficient Simulation of Light Transport in Scenes with Participating Media using Photon Maps[END_REF] are far from being real-time on a conventional desktop computer. With the growing capacities of graphics hardware, the real-time problem has been investigated.

Two approaches can be used to achieve this goal : volume rendering or direct representation. In order to add the volumetric shadows, the first approach will naturally use shadow maps techniques whereas the second is implicitly shadow volumes oriented [START_REF] Heidman | Real Shadows Real Time[END_REF]. Volume rendering is a classic solution to render a participating medium which is a volume de facto. Methods like [START_REF] Behrens | Adding Shadows to a Texture-based Volume Renderer[END_REF]- [START_REF] Nulkar | Splatting with Shadows[END_REF] represent density or illumination in voxels encoded into 2D or 3D textures. Accumulation techniques using textured slices or virtual planes are then used to display the result. These methods could produce nice images of clouds or gas, but, in addition to requiring a lot of texture memory, they are not suitable for shafts of light where sharp edges exist. Special methods are defined to render beams and shafts of light precisely and most of them [START_REF] Dobashi | A Simple, Efficient Method for Realistic Animation of Clouds[END_REF]- [START_REF] Lefebvre | Volumetric Lighting and Shadowing NV30 shader[END_REF] use volume rendering techniques along with sampling shadows in shadow maps. Unfortunately, they suffer from artifacts due to the sampling. Dobashi et al. [START_REF] Dobashi | Interactive Rendering of Atmospheric Scattering Effects Using Graphics Hardware[END_REF] present a very elegant solution to solve this problem using specialized adaptive sampling for shadows. They obtain an interactive rendering of participating media without aliasing or artifacts. However the image resolution remains low since the method is expensive in terms of fillrate. Moreover, the method works only with static lights due to the precomputation of shadow maps and only addresses the direct single scattering.

The algorithms belonging to the second approach compute directly on every point of the scene the contribution of the participating medium. It is well adapted to classical graphics engines since it consists in one more rendering of the scene. In this case, methods like [START_REF] Mech | Hardware-Accelerated Real-Time Rendering of Gaseous Phenomena[END_REF], [START_REF] Hoffman | Rendering Outdoor Light Scattering in Real Time[END_REF] use participating medium boundaries, or special virtual planes, combined with vertex and fragments shaders. Other methods focus on the rendering of the atmosphere [START_REF] O'neil | Accurate atmospheric scattering[END_REF]. Despite their speed and simplicity, all these methods consider only the direct single scattering. The most advanced method of this group is proposed by Sun et al. [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF] and is the first to consider the effects of the indirect single scattering. Unfortunately, it does not take shadows into account, even though it is real-time. Our work belongs to this group but is the only one of them to integrate direct single scattering, indirect single scattering and volumetric shadows. Compared to our own work presented in [START_REF] Biri | Real Time Rendering of Atmospheric Scattering and Volumetric Shadows[END_REF], we present here an hardware implementation along with the representation of the indirect single scattering which are major improvements.

III. THE SINGLE SCATTERING ILLUMINATION MODEL

In this section, we will present an analytic formulation of the single scattering illumination model. We will mainly focus on in-scattering since absorption, emission, multiple scattering and out-scattering are really simple to integrate. We start in the following developments with the angular formulation we presented in [START_REF] Biri | Real Time Rendering of Atmospheric Scattering and Volumetric Shadows[END_REF], [START_REF] Lecocq | Mathematical approximation for real-time rendering of participating media considering the luminous intensity distribution of light sources[END_REF], and that were finely used by [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF], work we review and extend here to obtain the complete illumination model of single scattering. Our developments consider only homogeneous participating media defined all over the scene and point light sources, assumed to be isotropic.

A. Single scattering along a ray

Considering a view ray immersed in a participating medium, the radiance L observed along u P from a point P can be written (see Figure 2) : L uP = L dss ( u P ) + e -ktOP L(P ) + L a e -ktOP + L oiss ( u P )

(1) where L dss is the direct single scattering and L(P ) the radiance of P attenuated by both absorption and out-scattering. k t is the extinction coefficient, sum of the absorption coefficient and the diffusion coefficient, properties of the participating medium. The third term represents both emission and multiple scattering inducing a constant radiance L a along the ray. The last term is the object indirect single scattering. Now we focus on the radiance of P that can be written :

L(P ) = L iss (P ) + e -ktSP L d (P ) + L a e -ktSP (2)
L iss is the indirect single scattering received on point P and L d (P ) an usual direct illumination, like Phong model, that is attenuated by absorption and out-scattering. The third term is once again the contribution of emission and multiple scattering on the illumination ray. Then, only L dss and L iss remain unknown. They both are the contribution of the first light scattering along, respectively, a view ray and an illumination ray. First of all, we will investigate the general contribution of this first scattering along a ray before to see its application on view rays and illumination rays. The integral transfer equation, presented in [START_REF] Siegel | Thermal Radiation Heat Transfert[END_REF], gives the incoming radiance L ss in direction u P stopped at vertex P and seen from O (see Figure 2) :

L ss ( u P ) = k t Ω 4π d 0 I S (α + β)e -ktx e -ktr r 2 p(α + π 2 )dx (3)
In this equation x, r, d, h are geometrical factors and Ω is the albedo, i.e. the fraction between the diffuse coefficient and the extinction coefficient of the medium. p is the phase function expressing for any incoming direction the ratio of light following the ray direction. Finally, I s is the directional intensity of the source computed relatively to the reference direction ω ref . Equation (3) expresses the in-scattering which is responsible for the subtle effects of atmospheric scattering.

The previous equation can be simplified [START_REF] Biri | Real Time Rendering of Atmospheric Scattering and Volumetric Shadows[END_REF], [START_REF] Lecocq | Mathematical approximation for real-time rendering of participating media considering the luminous intensity distribution of light sources[END_REF] using the angle α to minimize the dependency of the integral to it. Using the variable change x = t + h. tan(α), t and h are constant along the ray, we can obtain (see the annexes) :

L ss (φ, d, l) = k t Ωe -ktt 4πh γ d γ0 I S (α+β)p(α+ π 2 )e -kth sin(α)+1 cos(α) dα (4)
where

γ 0 = - π 2 + φ = atan( -t h ) γ d = atan( d -lcos(φ) lsin(φ) ) = atan( d -t h )

B. The direct single scattering

In the case of view rays, the point O is the position of the camera. Therefore, the distance between O and S remains constant whatever the view ray considered. The only variables in the equation ( 4) are the angle φ and the distance d. For an isotropic light, this equation becomes :

L dss (φ, d) = k t ΩI S e -ktt 4πh γ d γ0 p(α+ π 2 )e -kth sin(α)+1 cos(α) dα (5)
If we denote :

Γ(λ, γ) = γ 0 p(α+ π 2 )e -ktλ sin(α)+1 cos(α) dα (6)
the direct single scattering could be written :

L dss (d, φ) = k 2 t ΩIS 4π e -k t t kth Γ(k t h, atan( d -t h )) -Γ(k t h, atan( -t h )) (7)
This particular formulation is well suited to optimize the shader that computes the direct single scattering. Note that the 2D function Γ is purely numerical and only depends on the shape of the phase function. It can therefore be precomputed and stored in a 2D table. Figure 3 shows two examples for isotropic phase function and mie-hazzy phase function.

C. The indirect single scattering

Inspired by [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF], we also compute the illumination of vertex P due to the light scattering all over the scene. The indirect single scattering can be written (see Figure 4) as follows :

L iss (P ) = 2π θ=0 π φ=0 L ss (φ, d φ,θ , l)f r ( u c , u v )( u c . n)dω (8)
The first difficulty is to obtain for each direction u c the distance d φ,θ . This distance has been arbitrary fixed at ∞ in [START_REF] Sun | A practical analytic single scattering model for real time rendering[END_REF] but this is a coarse approximation. Nevertheless, it is necessary if we want to avoid too much complexity. So we will make the same assumption and take, d φ,θ = ∞, ∀φ, θ.

First of all, consider the Lambertian case where the function f r is the diffuse reflectivity of the surface k d . We then have :

L iss (P ) = k d 2π θ=0 π φ=0 L ss (φ, ∞, l)( u c . n)sin(φ)dφdθ (9)
Thanks to the symmetry on θ, we can choose the spherical coordinates centered on P such that L ss (φ, ∞, l)cos(φ n )cos(φ)sin(φ)dφ

Now we can substitute L ss with equation ( 4) to obtain :

L iss (P ) = 2πk d cos(φ n ) π φ=0 k 2 t I S Ω 4πl g(l, φ)cos(φ)dφ with g(λ, φ) = 1 λ π 2 φ-π 2 p(α+ π 2 )e -λ cos(α-φ)+1 cos(α)
dα And finally we have :

L iss (P ) = k 2 t ΩI S 4π 2πk d cos(φ n ) k t l Γ L (k t l) (10)
with

Γ L (λ) = π 0 cos(φ) π 2 φ-π 2 p(α+ π 2 )e -λ cos(α-φ)+1 cos(α) dαdφ (11)
The Phong model is much more difficult and we need the reparametrization of Ramamoorthi et al. [START_REF] Ramamoorthi | Frequency space environment map rendering[END_REF]. It gives use the following equation :

L iss (P ) = k s 2π θ=0 π φ=0 L ss (φ ′ , ∞, l)cos n (φ)sin(φ)dφdθ
where n is the shininess and k s the specular reflectivity of the surface. φ and θ are defined relatively to the axe u r . φ ′ is the angle formed by vector u s and u c (see Figure 5). We can choose the base to have :

φ ′ = acos(cos(φ S )cos(φ) + sin(φ)cos(θ) 1 -cos 2 (φ S ))
Therefore, the equation will be :

L iss (P ) = k 2 t ΩI S 4π k s k t l Γ P (k t l) (12)
with

Γ P (φ ′ , λ)= 2π 0 π 0 π 2 φ ′ -π 2 p(α+ π 2 )e -λ cos(α-φ ′ )+1 cos(α)
cos n (φ)dα (13)

D. Shadows on view rays

We want to integrate now the effect of occlusions along any view rays. Equation ( 7) describes the particular case where the view ray remains totally lit. In order to integrate shadow volumes, we need to consider more general cases, illustrated in Figure 6. Indeed, the ray must be split into lit and shadowed parts. In this example, the medium contribution along the ray is split into two parts on OA and BP. Using the laplace formula, it is straightforward to see that equation [START_REF] Nulkar | Splatting with Shadows[END_REF] becomes :

L dss (d, φ) = k 2 t ΩI S 4π e -ktt [Γ(k t h, γ P ) -Γ(k t h, γ B ) +Γ(k t h, γ A ) -Γ(k t h, γ O )] (14)

IV. HARDWARE IMPLEMENTATION

A. Overview of our method

Our algorithm is easy to implement. We present here the 5 steps of this method and we will precise for each step if the computation is done by CPU or by GPU.

1) (CPU) The silhouettes of every moving shadow caster are computed. If the light source is moving, every silhouette needs to be recomputed. 2) (GPU) Scene is rendered using the conventional polygonal rendering method to obtain the direct illumination and the indirect single scattering. Surface shadows can be obtained using shadow planes algorithms [START_REF] Heidman | Real Shadows Real Time[END_REF], [START_REF] Everitt | Practical and Robust Shadow Volumes[END_REF].

The stencil buffer now contains lit areas of the scene. An ambient fog is added to take into account both absorption and multiple scattering. 3) (GPU) Scene is rendered once more time and direct single scattering is computed for each vertex of the scene. Depth test is set to equality. Only lit parts of the scene are rendered thanks to the stencil buffer. The depth test function accepts only planes that are closer to the camera. Front facing planes add their contribution when back facing planes subtract them. Stencil function is set to allow fragments if the stencil is equal to 1 for front facing planes and 0 for back facing ones. Front facing planes always decrement the stencil buffer and back facing ones always increment it. All stages have to be done for each light source. As in [START_REF] Everitt | Practical and Robust Shadow Volumes[END_REF], a initialization stage is done to obtain the ambiant lighting and the first depth map. Each stage is detailed in the following sections.

B. Computation of silhouettes (step 1)

In our algorithm, we select some objects to be shadow casters. Their silhouettes are easily computed by determining all edges of their mesh common to a front-facing triangle regarding the light position and one back facing it. Then all these edges are linked together if possible, and stored in a loop list. In order to obtain correct silhouettes, we need well designed closed triangular meshes (2-manifold) for which connectivity information are available. These conditions for the shadow casters are the ones indicated in [START_REF] Everitt | Practical and Robust Shadow Volumes[END_REF]. Shadow planes are infinite quads formed by a silhouette edge and the light position. They are constituted by the two edge's vertices and two other points, projections of the previous vertices to infinity toward direction : light position -vertex [START_REF] Heidman | Real Shadows Real Time[END_REF]. They are oriented toward the unshadowed area of the scene. As we need to compute the medium contribution on all shadow planes, it is wise to use shadow plane silhouettes rather than the shadow planes of all little triangles. Of course, if the light does not move, only moving shadow caster silhouettes have to be computed. Finally, if the input geometry is modified by graphics hardware, using displacement mapping for example, a solution to obtain silhouettes of all objects quickly and accurately can be found in [START_REF] Brabec | Shadow Volumes on Programmable Graphics Hardware[END_REF].

C. Rendering the scene (step 2 and 3)

The second and third steps of the algorithm compute the illumination on surfaces. This can be done using one or

D. Sorting the shadow planes (step 4)

Before rendering all shadow planes, we have to make sure that we will not render shadow planes, or part of them, that are themselves in shadow. If we do not care about this problem, it will create artifacts we call shadow in shadows, shadow planes that are in shadow must not be rendered. Then, we render the shadow planes, back-or front-facing, in a "back to front" order and use the stencil buffer to avoid the drawing of shadowed shadow planes.

E. Rendering the shadow planes (step 5)

It remains to render the shadow volumes, i.e. to take into account the equation ( 14). We still have the stencil we have obtained in the stage 2. Shadow planes are rendered in the order defined in the previous stage. The same pixel shader of Figure 7 is used for the shadow planes. The only difference is blending : front facing planes add their contribution when back facing planes subtract them.

Taking into account correctly the shadow in shadow problem requires the use of the stencil buffer. The key idea is to draw only the planes that constitute the boundary between lit volumes and shadowed ones. This can be done using our ordering and the value stored in the stencil. We define the stencil test such that front facing planes pass the stencil test if its value is one, representing shadowed area, and back facing ones passes if it equals zero, value representing lit area. Ideally the back (resp. front) facing quads should always add (resp. subtract) one to the stencil buffer if it passes depth test. Unfortunately, the stencil test is before the depth test so we have to do our own depth test in the pixel shader. Thus, we only draw shadow planes that make the shadow volume boundaries. The indicated strategy works if the camera is in the light. A slightly different strategy can be used when the camera is in shadow but the philosophy remains the same.

V. RESULTS

The previous algorithm has been implemented on a standard computer using a 2.6 GHz processor and an ATI 9800 PRO graphics card (which is an old card now!). All images and videos we present have a 800x600 resolution. First of all, we point out, in Figure 10 the influence of each part of our single scattering illumination model. In this example, we have not considered the multiple scattering to concentrate on the other effects. On the top row, the first column shows the directe illumination of surfaces. The second column adds the direct single scattering. This effect highlights slightly all lit surfaces. Therefore, if we don't consider shadow volume, discontinuity in the illumination can occur between lit and shadowed areas, like in the bottom of the jug. After the shadow volumes rendering, illustrated in the third column, the discontinuities disapear. All these pictures have been rendered with indirect single scattering, that introduces subtle modifications : dimming of the specular highlights on the vase and brightening of darker region.

We also present some snapshots of our animations. The first image in Figure 11.a. is a simple scene, where a pen is bumping in front of two lights. It illustrates a classical situation where well design 3D objects are moving and casting shadows. Here the shadow planes are really intricated. This scene is rendered at more than 30 fps, including direct single scattering, indirect single scattering and of course volumetric shadows rendering. The image in Figure 11.b is a snapshot of an animation presenting two point light sources, including one that moves, in a complex scene containing about 100 000 triangles. The last picture 11.c is another snapshot from an animation where camera moves and properties of the participating medium evolve. The table II

VI. CONCLUSION

We have presented in this paper a complete single scattering illumination model along with a new algorithm able to render the main part of this model. It considers a single participating medium recovering the scene and lit by one or several, eventually dynamic, point light sources. Our algorithm is fast enough to handle more than 30 frames per second for moderately complex scenes. It implements the direct single scattering, the indirect single scattering and, most of all, volumetric shadows along with surfacic shadows. Our method can be implemented in programmable graphics hardware, and compared to volume rendering approach does not need a lot of texture memory. Therefore it can easily be integrated in any graphics engine. We also have to point out that our algorithm does not create the aliasing effect we can have with volume rendering techniques thanks to the use of the exact shadow planes.

The perspectives are numerous but two developments seems promising : Relax assumptions : So far we consider only isotropic point light and homogeneous participating media. We have made satisfying preliminary work about animating several bounded participating media and it could be a way to handle heterogeneous medium. Directional light source are easier to handle than point light and will not be a problem. The real challenge is the integration of directionnal point light sources. Make it soft ... and quick : An other way to obtain more realistic image will be to consider soft shadows, for surfacic shadows but also for volumetric shadows. Fortunately, our work is well adapted to algorithm such as [START_REF] Assarson | A Geometry-based Soft Shadow Volume Algorithm using Graphics Hardware[END_REF] since it uses shadow volumes. Finaly, note that no optimization have been done to render the scenes. We believe that occlusion culling, frustum culling and clustering would greatly speed up our algorithm.

Fig. 1 .

 1 Fig. 1. The same scene lit a. (left) classically, b. (center) with single scattering and c. (right) with single scattering and volumetric shadows. Who could say precisely where are the lights in the left picture ? This is obvious in the right picture

Fig. 2 .

 2 Fig. 2. Notations for our model

Fig. 3 .Fig. 4 .

 34 Fig. 3. The Γ function for isotropic and hazzy phase function

Fig. 5 .Fig. 6 .

 56 Fig. 5. Notations for Phong illumination model

Fig. 8 .

 8 Fig. 8. Shader for L iss computation

Fig. 9 .

 9 Fig. 9. Use of the stencil buffer in the rendering of shadow planes

Fig. 10 .

 10 Fig. 10. Insight of the influence of the different contributions. The first column is the result after step 2, the second presents the addition of the direct single scattering and the third adds the volumetric shadows

  presents the FPS compared to the number of triangles of those scenes.

	scene fig. 11 .a	without with dss and iss shadows 31	with iss, dss and shadows 30	number triangles 14 785	of
	fig. 11 .b	10.5	10.2	110 014
	fig. 11 .c	7.6	7	136 266
			TABLE II
		FPS COMPARED TO NUMBER OF TRIANGLES

ACKNOWLEDGMENT

We thank Pascal Lecocq for the developments done in his thesis about this subject and wish him a excellent career in Spain. We also thank Leonie Van Royeen for her english review.

two rendering of the scene depending on the ability of the graphic cards to handle large shaders. Indeed, the number of instructions in the shader for both indirect single scattering and direct single scattering computations could exceed the capacity of the graphic card.

We show in Figure 7 and 8 the two GLSL fragment shaders used for steps 2 and 3 and which "implement" respectively equation ( 1) and [START_REF] Nulkar | Splatting with Shadows[END_REF]. These shaders, and the corresponding steps, can be grouped together if possible. Note that we use a shader for any kind of participating medium of the scene. The two shaders need four textures : a 1D texture used as a lookup table for arctangent function, a 2D texture representing the function Γ defined in [START_REF] Fedkiw | Visual Simulation of Smoke[END_REF], 1D texture to store the Lambertian function [START_REF] Dobashi | A Simple, Efficient Method for Realistic Animation of Clouds[END_REF] and finally one for the Phong function [START_REF] Everitt | A Fast Algorithm for Area Light Source Using Backprojection[END_REF]. Note that in equations ( 7), [START_REF] O'neil | Accurate atmospheric scattering[END_REF], [START_REF] Dobashi | Interactive Rendering Method for Displaying Shafts of Light[END_REF] the first fraction could be a uniform variable and the second fraction could be computed directly by the shaders. Moreover we have avoided the direct computation of atan, cos or sin functions in these shaders,


Biri Venceslas

Gaspard Monge

Higher Order Radiosity : Using Higher Order Methods in Computer Graphics

Keywords: Higher Order Methods, Dynamic Radiosity, Temporal Coherence

The light ambiance of a virtual scene is the key of realistic computer generated images. But the computation of the whole light exchanges between elements of a complex scene constitutes a real challenge in computer graphics. Usually, in real-time video games for example, the modelisation of these exchanges remains quite simple. And achieve a realistic representation of lighting requires to solve the global illumination problem which means to compute all light exchanges.

), the illumination function (the radiosity) is separated from the geometric surface [8] assuming the surface to be parametric. This talk will present how higher order functions bases are used in computer graphics and especially in radiosity, and why they are efficient. We will see also how they are used to handle dynamic moves in virtual scenes, taking profit of temporal coherence. Figure 1. Two scenes lighted with radiosity using higher order functions bases

The term " radiosity " gather algorithms that are able to solve the global illumination problem in the case where all light exchanges are diffuse [START_REF] Goral | Modeling the interaction of light between diffuse surfaces[END_REF]. Radiosity associates each surface to a function modeling the light energy emited per square meters in every direction. And this function is also called radiosity. Traditionnal methods [START_REF] Cohen | The Hemi-Cube : A Radiosity Solution for Complex Environments[END_REF][START_REF] Nishita | Continuous Tone Representation of Tree-Dimentional Object Taking Account of Shadows and Interreflection[END_REF] divide surface into small patches where the radiosity is considered constant. But this approach leads to inefficient algorithms due to the high number of patches. Several algorithms [START_REF] Cohen | A Progressive Refinement Approach for Fast Radiosity Image Generation[END_REF][START_REF] Wallace | A Ray Tracing Algorithm for Progressive Radiosity[END_REF][START_REF] Hanrahan | A Rapid Hierarchical Radiosity Algorithm[END_REF][START_REF] Lischinski | Combining Hierarchical Radiosity and Discontinuity Meshing[END_REF] allow quicker computations using progressive or hierarchical methods but they still consider constant radiosity functions.

Higher order functions bases have finally be used to represent the radiosity functions. Such higher order methods [START_REF] Zatz | Galerkin Radiosity: A Higher Order Solution Method for Global Illumination[END_REF][START_REF] Gortler | Wavelet Radiosity[END_REF][START_REF] Yizhou | Multiresolution B-Spline Radiosity[END_REF] lead to powerfull algorithms that are able, not only to compute efficiently the radiosity function [START_REF] Franc ¸ois Cuny | A novel approach makes higher order wavelets really efficient for radiosity[END_REF] of all surfaces, but also give us the possibility to handle dynamic scenes [START_REF] Biri | Dynamic Radiosity using Higher Order Functions Bases and Temporal Coherence[END_REF]. For N surfaces representing the scene, the radiosity in a point x for a surface i can be written [START_REF] Gortler | Wavelet Radiosity[END_REF] :

B i (x) = E i (x)+ρ i (x) N j=1 Sj K ij (s,t,u,v)B j (u,v)du dv (1)
where E i and ρ i are the exitance (self illumination) and the reflectivity of the surface. The kernel function K is the product of a geometric term, the form factor, the differential area A and a visibility term V :

K ij (s,t,u,v) = - ( n i . r)( n j . r) π r 4 V ij (s,t,u,v) A j (u,v) (2)
The kernel function can be analytically computed only for very specific positions between two surfaces and usually needs a numerical computation. Then, if each function is projected over an orthonormal set of base functions {Γ k (s, t), k = 1..N }, we can project the equation ( 1) over the k th function of Γ which leads to the equation :

b k i = e k i + ρ i j,l b l j K kl ij (3)
where

K kl ij = K ij Γ l |Γ k (4)
The b k i are the coefficients of the radiosity functions and so are the unknowns we seek. All these equations, for each surface i and function k can be summarize in a system that have to be solved in order to obtain our unknowns. And this is far from being easy since matrix coefficients are difficult to compute (this is a double surface integral).


Yukiko Kenmochi
email: y.kenmochi@esiee.fr

Yusuke Nomura

Local Configurations in Discrete Combinatorial Surfaces

Keywords: discrete surface, local configurations, polyhedral complex

Representing discrete objects by polyhedral complexes, we can define all conceivable topological characteristics of points in discrete objects, namely those of vertices of polyhedral complexes. Such a topological characteristic is determined for each point by observing a configuration of object points in the 3 × 3 × 3 local point set of its neighbors. We study a topological characteristic such that the point is in the boundary of a 3D polyhedral complex and the boundary forms a 2D combinatorial surface. By using the topological characteristic, we present an algorithm which examines whether the central point of a local point set is in a combinatorial surface, and show how many local point configurations exist in combinatorial surfaces in a 3D discrete space.

Introduction

Surfaces are often used in the three-dimensional image analysis, sometimes implicitly. For example, the method called active balloon uses a deformable surface for the three-dimensional image segmentation [START_REF] Cohen | On active contour models and balloons[END_REF]. For such image segmentation, first an initial surface is set and then its shape is changed by using geometric characteristics such as curvatures. In many cases, such a deformable surface is realized by a continuous surface which is interpolated from discrete points in a three-dimensional digital space, such as a spline surface, and its shape is changed by moving its control points in a continuous space [START_REF] Mcinerney | Deformable models in medical images analysis: a survey[END_REF]. If we need a result such as a segmented object in a discrete space after its continuous deformation, we therefore rediscretize the continuous surface to obtain a discrete object in a digital image. Obviously, continuous surface models are not efficient for the case such that both input and output are required to be discrete. It is more efficient to use a discrete surface model which allows us to make discrete deformation [START_REF] Esnard | Discrete deformable boundaries for 3D image segmentation[END_REF]. In this paper, we study discrete combinatorial surfaces and their local configurations, especially in the 3 × 3 × 3 point region. We insist that such local configurations determine local geometric parameters such as curvatures in a discrete space. Therefore, if we can obtain all local configurations appearing in discrete surfaces a priori, we can expect to obtain useful properties which may give us efficient algorithms for calculation of geometric parameters as well as those for surface deformation based on such parameters.

Let us consider a three-dimensional discrete space Z 3 , consisting of lattice points whose coordinates are all integers in a three-dimensional Euclidean space R 3 . In our previous work [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF], we presented a boundary extraction algorithm which provides a triangulation (or polygonization) of a set of boundary points given by

Br m (V) = {x ∈ V : N m (x) ∩ V = ∅} (1)
where V is the input, i.e., a discrete object in Z 3 and V is the complement. N m (x) is the m-neighborhood of a point x in Z 3 , defined by N m (x) = {y ∈ Z 3 : xy 2 ≤ t} where t = 1, 2, 3 for m = 6, 18, 26 respectively. Applying discrete polyhedral complexes [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF] based on combinatorial topology [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF] to object representation, we can obtain topologies for boundary points. With a help of the topologies, we found that boundary points in Br m (V) include not only surface points, i.e., points on 2-dimensional combinatorial manifolds, but also non-surface points, i.e., singular points, as shown in Fig. 1. In this paper, we use local topological notions similarly to our work [START_REF] Kenmochi | Discrete combinatorial geometry[END_REF] to discriminate surface Fig. 1. Examples of local configurations in a 3×3×3 point region so that the central point is considered to be a boundary point [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF]; a surface point (left), a surface point but not a simplicity surface point [START_REF] Couprie | Simplicity surface: a new definition of surfaces in Z 3[END_REF] (center), and a non-surface point, i.e., a singular point (right).

points from boundary points. Such notions enable us to present an algorithm to count the local point configurations appearing in discrete combinatorial surfaces for the 6-, 18-and 26-neighborhood systems; discrete combinatorial surfaces are defined for any m-neighborhood system, m = 6, 18, 26, such that their adjacent vertices are m-neighboring. Because there are a finite number of points in the local region, there must be a finite number of local point configurations in discrete combinatorial surfaces.

For the 6-neighborhood system, the definition of discrete combinatorial surfaces is given by Françon in [START_REF] Françon | Discrete combinatorial surfaces[END_REF] and he showed that there are 6 local configurations of discrete surfaces for the 6-neighborhood system as illustrated in Fig. 2. Note that the similar results are obtained by using different approaches, for example, in [START_REF] Imiya | The Euler characteristics of discrete objects and discrete quasi-objects[END_REF][START_REF] Jonker | The generation of N dimensional shape primitives[END_REF]. The discrete deformation model based on such discrete surface configurations for 6-neighborhood system is also presented in [START_REF] Esnard | Discrete deformable boundaries for 3D image segmentation[END_REF]. Moreover, Françon and Kenmochi et al. show that there are five local configurations which appear in discrete planes, illustrated as the five left configurations in Fig. 2 [9,13]. In other words, there is only one configuration (the most right one in Fig. 2) which does not appear in discrete planes but appears in discrete non-planar surfaces.

In [START_REF] Françon | Discrete combinatorial surfaces[END_REF], however, the 18-and 26-neighborhoods are not practically treated so that we do not see how to generate discrete combinatorial surfaces for the 18-and 26-neighborhood systems, even if the mathematical definition is given for any neighborhood system. Morgenthaler et al. defined discrete surfaces by using the point connectivity based on the Jordan surface theorem; any Jordan surface divides the space into two [START_REF] Morgenthaler | Surfaces in three-dimensional digital images[END_REF]. In [START_REF] Couprie | Simplicity surface: a new definition of surfaces in Z 3[END_REF], Couprie et al. pointed out that, for the 26-neighborhood system, Morgenthaler's discrete surfaces have only 13 local configurations while their discrete surfaces, called simplicity surfaces, have 736 configurations. However, we see that even simplicity surfaces do not give enough configurations if we would like to treat our boundary points. For example, we obtain a boundary point by applying our boundary extraction algorithm [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF] as illustrated in Fig. 1 (center) and we see that it is not considered to be a simplicity surface. Ciria et al. also presented a graphbased notion of discrete surfaces for the 26-neighborhood system, though the number of all local point configurations has not been studied yet [START_REF] Ciria | Local characterization of a maximum set of digital (26, 6)-surfaces[END_REF].

To tackle the problem above, we take the approach based on polyhedral complexes similarly to our previous work [START_REF] Kenmochi | Combinatorial topologies for discrete planes[END_REF][START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF]. The organization of our pa-Fig. 2. The 6 local configurations of discrete combinatorial surfaces for the 6-neighborhood system [START_REF] Françon | Discrete combinatorial surfaces[END_REF][START_REF] Imiya | The Euler characteristics of discrete objects and discrete quasi-objects[END_REF][START_REF] Jonker | The generation of N dimensional shape primitives[END_REF] where the five left ones appears in discrete planes [9,[START_REF] Kenmochi | Combinatorial topologies for discrete planes[END_REF].

per is as follows. In Section 2, we first define discrete polyhedral complexes [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF], based on combinatorial topology [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF], which gives a topology to a three-dimensional discrete space. In Section 3, we then study topological characteristics of discrete polyhedral complexes similarly to [START_REF] Kenmochi | Discrete combinatorial geometry[END_REF]. By using these topological characteristics, in Section 4, we consider all conceivable topological characteristics of object points and distinguish those of boundary points such that the boundary forms a discrete combinatorial surface. After presenting an algorithm to examine whether a local point set forms a discrete combinatorial surface, we count local point configurations in discrete combinatorial surfaces, called local surface configurations, for each neighborhood system in Section 5.

Finally, we show that such local surface configurations for the 6-neighborhood system are the same 6 ones as illustrated in Fig. 2 and derive new results for the 18-and 26-neighborhood systems. We also discuss the utilities of such study on local configurations in discrete surfaces for three-dimensional shape analysis in Section 6. This paper is an extended version of [START_REF] Kenmochi | Local point configurations of discrete combinatorial surfaces[END_REF].

Discrete polyhedral complexes

In combinatorial topology, any object in R 3 is represented by a set of simplexes [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF] or more generally polyhedra [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF]. An n-dimensional simplex is considered to be an n-dimensional element in R 3 where n can be from 0 to 3; 0-, 1-, 2-and 3-dimensional simplexes are defined as isolated points, line segments, triangles, and tetrahedra, respectively. A set of simplexes which are combined together without contradiction is called a complex [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF]. Replacing simplexes with convex polyhedra, we obtain polyhedral complexes [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Ziegler | Lectures on Polytopes[END_REF] instead of simplicial complexes.

In this paper, we construct convex polyhedra whose vertices are lattice points in Z 3 and adjacent vertices are m-neighboring for m = 6, 18, 26. Such convex polyhedra are called discrete convex polyhedra, and a discrete polyhedral complex is constructed as a set of discrete convex polyhedra combined together without contradiction [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF].

Convex polyhedra and polyhedral complexes in R 3

For the definitions of convex polyhedra and polyhedral complexes in R 3 , we follow the notions in [START_REF] Ziegler | Lectures on Polytopes[END_REF]. Similar notations are also seen in [START_REF] Alexandrov | Combinatorial Topology[END_REF][START_REF] Stillwell | Classical Topology and Combinatorial Group Theory[END_REF].

Definition 1 A convex polyhedron σ is the convex hull of a finite set of points in some R n .

The dimension of a convex polyhedron σ is the dimension of its hull. An n-dimensional convex polyhedron σ is abbreviated to an n-polyhedron. For instance, a point is a 0-polyhedron, a line segment is a 1-polyhedron, a triangle or a square is a 2-polyhedron, and a tetrahedron or a hexahedron is a 3polyhedron.

Definition 2 Let σ be an n-polyhedron. A linear inequality a • x ≤ z is said to be valid for σ if it is satisfied for all points x ∈ σ. A face of σ is defined by any set of the form

δ = σ ∩ {x ∈ R 3 : a • x = z} where a • x ≤ z is valid for σ.
For instance, a 3-polyhedron which is a tetrahedron has four 0-polyhedra, six 1-polyhedra and four 2-polyhedra for its faces. If an n-polyhedron τ is a face of σ, τ is called an n-face and such a binary relation is denoted by τ ≺ σ.

The point of a 0-polyhedron, the endpoints of a 1-polyhedron and the vertices of 2-and 3-polyhedra are called the vertices of each convex polyhedron.

Definition 3 A polyhedral complex C is a finite collection of convex polyhedra such that (1) the empty polyhedron is in C, ( 2
) if σ ∈ C and τ ≺ σ, then τ ∈ C, (3)
if σ, τ ∈ C, then the intersection σ ∩ τ is a common face of σ and τ .

The dimension of C is the largest dimension of a convex polyhedron in C.

We see that any C is a partially ordered set which can be identified with a topological space called a discrete space; the detail is found in Section 6 of Chapter 1 in [START_REF] Alexandrov | Combinatorial Topology[END_REF].

Discrete convex polyhedra

We consider discrete polyhedral complexes which are polyhedral complexes such that vertices of convex polyhedra are all lattice points in Z 3 and adjacent vertices are m-neighboring for m = 6, 18, 26. The constraints allow us to look for a discrete convex polyhedron which is not larger than the unit cubic region as follows. The explication is also found in [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF].

Let us consider all possible convex polyhedra in a unit cubic region such that the vertices of each convex polyhedron are vertices of a unit cube. A unit cube has eight lattice points for its vertices. For each lattice point we assign the value of either 1 or 0 and call the point a 1-or 0-point, respectively. We then obtain 2 8 different 1-point configurations in a unit cube and, up to rotations and symmetries, we can reduce the number of different configurations to 22. Note that, among the 22 configurations, there are 21 configurations including at least one 1-point. For each of the 21 configurations, we obtain a convex polyhedron such that vertices of the polyhedron are 1-points. We then classify the 21 convex polyhedra with the dimension of n = 0, 1, 2, 3 and with the mneighboring relations of adjacent vertices for m = 6, 18, 26 as shown in Table 1. We see that all 21 convex polyhedra are classified for m = 26 and parts of them are classified for m = 6, 18.

For any neighborhood system, 0-dimensional discrete convex polyhedra are isolated points. For 1-, 2-and 3-dimensional discrete convex polyhedra, in the case of the 6-neighborhood, the distance between any adjacent vertices is equal to 1, and in the case of the 26-neighborhood, it is equal to or less than √ 3. Hereafter, we abbreviate n-dimensional discrete convex polyhedra to n-dimensional discrete polyhedra or simply n-polyhedra.

Discrete polyhedral complexes

We construct a discrete polyhedral complex which is a finite collection of discrete convex polyhedra satisfying the three conditions in Definition 3 for each m-neighborhood system. In this paper, we need to construct a discrete polyhedral complex only in the local region such as 3 × 3 × 3. Even though we present an algorithm for construction of a discrete polyhedral complex in such a local region in this subsection, we remark that the same algorithm can be applied to three-dimensional binary images of various sizes. Mathematical details for obtaining discrete polyhedral complexes are found in [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF].

Let us consider a point x ∈ Z 3 and the 3 × 3 × 3 point region around x, that is N 26 (x). Given a subset V ⊆ N 26 (x), we explain how to construct a discrete polyhedral complex C m for m = 6, 18, 26 from V in this subsection.

We first consider the case of the 26-neighborhood, and then show how to modify our algorithm for the 18-neighborhood. For the 6-neighborhood, we refer to another method which is simple for constructing C 6 from V.

Hereafter, we abbreviate an n-dimensional discrete polyhedral complex to an n-dimensional discrete complex or simply an n-complex.

Discrete complex construction for the 26-neighborhood system

We first consider a 2 × 2 × 2 unit cubic region which is a set of eight lattice points for any x = (p, q, r) in Z 3 such that

D(x) = {p, p + 1} × {q, q + 1} × {r, r + 1}.
Obviously, D(x) includes x itself. Therefore, any 3×3×3 point region N 26 (x) is divided into eight unit cubic regions D(y) such that

N 26 (x) = y∈W(x) D(y)
where

W(x) = D(x -(1, 1, 1)).
All points in V (resp. the complement V) are 1-points (resp. 0-points). Setting the set of 1-points in a unit cube D(y) as

U(y) = D(y) ∩ V,
we first construct a discrete complex from U(y) in a unit cube D(y), denoted by cmp 26 (U(y)). We then obtain a discrete complex C 26 with respect to V, which is also denoted by cmp 26 (V), such that

C 26 = y∈W(x) cmp 26 (U(y)) (2)
for x is the central point in N 26 (x). In [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF], mathematical proof of ( 2) is given.

A discrete complex cmp 26 (U(y)) is constructed from U(y) by the following algorithm. Let P n be a set of n-polyhedra in cmp 26 (U(y)) for n = 0, 1, 2, 3.

Note that |A| denotes the number of elements of a set A and that V er(σ) denotes a set of 0-faces (or lattice points) of an n-polyhedron σ, namely vertices of σ.

Algorithm 1

Input: A 1-point set U(y) in a unit cubic region D(y).

Output: A discrete complex cmp 26 (U(y)). begin 1 if |U(y)| ≥ 5
, then create a 3-polyhedron as a convex hull of all points in U(y) and put it in P 3 ; 2 else if |U(y)| = 4 and the four points are not coplanar, then create a 3-polyhedron as a convex hull of all points in U(y) and put it in P 3 ; 3 else if |U(y)| = 4 and the four points are coplanar, then create a 2-polyhedron as a convex hull of all points in U(y) and put it in P 2 ; 4 else if |U(y)| = 3, then create a 2-polyhedron as a convex hull of all points in U(y) and put it in P 2 ; 5 else if |U(y)| = 2, then create a 1-polyhedron as a convex hull of all points in U(y) and put it in P 1 ; 6 for each 3-polyhedron σ ∈ P 3 , look for all 2-faces τ of σ such that V er(τ ) is a three-or four-point subset of V er(σ) where all points in V er(τ ) are coplanar and all points in V er(σ) \ V er(τ ) are in one side of the plane spanned by points in V er(τ ), and put all such τ in P 2 ; 7 for each 2-polyhedron σ ∈ P 2 , 7.

1 if |V er(σ)| = 3, consider all 1-faces τ of σ such that V er(τ ) is any two-point subset of V er(σ)
, and put all τ in P 1 ; 7.2 if |V er(σ)| = 4, look for all 1-faces τ of σ such that the distance between the two points of V er(τ ) is one, or such that the distance between the two points of V er(τ ) is √ 2 when there exist exactly two pairs of points whose distance is 1 in V er(σ), and put all such τ in P 1 ; 8 set P 0 = {{z} : z ∈ U(y)}; 9 obtain cmp 26 (U(y)) = n=0,1,2,3 P n . end Table 2 Five 1-point configurations of a unit cube where different discrete complexes are constructed for the 18-and 26-neighborhood systems.

N18 N26

a discrete complex in a unit cube In Steps 6 (resp. 7), we also store binary relations ≺ between 3-and 2polyhedra (resp. 2-and 1-polyhedra).

Discrete complex construction for the 18-neighborhood system

Let us consider the case of m = 18. The right column of Table 1 shows that there are 22 different 1-point configurations of a unit cube for discrete complex construction in the case that m = 26. It also shows that there are five 1-point configurations missing for the 18-neighborhood system in comparison with those for the 26-neighborhood system. This is because each of those five convex polyhedra includes a pair of adjacent vertices which are 26-neighboring but not 18-neighboring. Therefore, to construct a discrete complex in a unit cube whose 1-point configuration is one of those five for the 18-neighborhood system, we cannot apply Algorithm 1. Table 2 illustrates the discrete complexes expected for the 18-neighborhood, which are different from those for the 26-neighborhood.

In order to obtain each discrete complex in the first line of Table 2, we first decompose a unit cube into two cubes as shown in Table 3. We then apply Algorithm 1 for each unit cube except for the following exceptional cases.

Let us consider a unit cube D(y) which is the last case in Table 3. At an adjacent cube D(z) which shares four 1-points with D(y), we need to modify its polyhedral complex cmp 18 (U(z)). Figure 3 illustrates that if cmp 18 (U(z)) has a squared 2-polyhedron σ whose vertices are the four shared 1-points, we must replace σ with two triangle 2-polyhedra τ 1 , τ 2 and their common face τ 3 which is a 1-polyhedron, such that cmp 18 (U(z)) \ {σ} ∪ {τ 1 , τ 2 , τ 3 }. We remark that such a replacement does not destroy the topology of a discrete complex. 2 is decomposed into two unit cubes for discrete complex construction in the case of the 18-neighborhood system.

original cube cube decomposition

Table 4 A unit cube of the last line in Table 3 is decomposed into three unit cubes for discrete complex construction in the case of the 18-neighborhood system, if unit cubes are adjacent as illustrated in Fig. 4 (left).

original cube cube decomposition

If an adjacent cube D(z) also has the last configuration in Table 3, that is, both D(y) and D(z) have the last 1-point configuration in Table 3 and they share four 1-points, and if the distance between two 1-points which are not shared by D(y) and D(z) is √ 5, as illustrated in the left figure of Fig. 4, we modify the polyhedral complexes in both cubes, cmp 18 (U(y)) and cmp 18 (U(z)), as illustrated in the right figure of Fig. 4. For such a modification, we decompose each unit cube of D(y) and D(z) into three unit cubes as shown in Table 4 before applying Algorithm 1 at each unit cube.

We remark that the modification of polyhedral complexes as illustrated in Fig. 4. Modification of polyhedral complexes in two adjacent unit cubes such that each of them has the last 1-point configuration in Table 3, and they have four common 1-points and also have two non-common 1-points whose distance is √ 5, for the 18-neighborhood system. Fig. 4 causes the problem such that a discrete complex C 18 constructed from V ⊆ N 26 (x) does not depend only on the 3 × 3 × 3 point region, namely N 26 (x), when x is the leftmost (or rightmost) point in Fig. 4. We will discuss how we deal with this problem for topological classification of local point configurations at the end of Section 4.

Discrete complex construction for the 6-neighborhood system

In order to construct a discrete complex C 6 from a 1-point set V, we can have a geometric algorithm similar to Algorithm 1 which is applied at each unit cube D(y) to construct cmp 6 (U(y)) and then obtain C 6 by (2). However, it is more efficient to apply a method to construct a Khalimsky topology [START_REF] Khalimsky | Computer graphics and connected topologies on finite ordered sets[END_REF] or a partially order set associated to Z n [2,3]; our discrete polyhedral complex for the 6-neighborhood system is equivalent to the above partially order set and it gives a topology equivalent to a Khalimsky topology. In [START_REF] Bertrand | New notions for discrete topology[END_REF] it is mentioned that a three-dimensional array whose size is 5 × 5 × 5 is used for storing such a topology for a 3 × 3 × 3 local point set.

Properties of discrete polyhedral complexes

We will present several properties of discrete complexes which we will need in the following sections.

Definition 4 Let C be an n-complex. If we have at least one n-polyhedron σ ∈ C for every s-polyhedron τ ∈ C such that τ ≺ σ, C is said to be pure.

Figure 5 shows examples of pure and non-pure discrete complexes. The 3complex in Fig. 5 (a) is not pure because it includes 0-, 1-and 2-polyhedra which do not belong to any 3-polyhedron. If we remove these 0-, 1-and 2polyhedra from Fig. 5 (a), we obtain a pure 3-complex in Fig. 5 (b).

We need the following notion of combinatorial closure to obtain a discrete complex from a set of discrete polyhedra, and by using this notion we present the connectivity of a discrete complex. binatorial closure of B is defined as

Cl(B) = B ∪ {τ ∈ C : τ ≺ σ, σ ∈ B}.
Note that B may not be a complex while Cl(B) is always a complex, and B = Cl(B) if B is not a complex. Figure 6 shows an example of the closure of a subset for a given complex.

Definition 6 Let C be a discrete complex, and σ, τ be arbitrary elements in C. We say that C is connected, if we have a path σ = a 1 , a 2 , . . . , τ = a n which satisfies the following conditions:

(1) a i ∈ C for every i = 1, 2, . . . , n;

(2) Cl({a i }) ∩ Cl({a i+1 }) = ∅ for every i = 1, 2, . . . , n -1.

Topological point characterization on discrete complexes

Let Sk(C) be the union of all 0-polyhedra in a discrete complex C, namely the union of V er(σ) for all discrete convex polyhedra σ ∈ C, called the skeleton of C. Obviously,

Sk(C) = V if C is made from V.
The goal of this paper is to present an algorithm to verify whether each point in Sk(C) is considered to be on a discrete surface or not. For that goal, we study topological characterization of each point in Sk(C) by observing its local point configuration and investigate all topological characteristics which can be maintained by points in Sk(C). In Z 3 , we have discrete complexes whose dimensions can be from zero to three. Thus, we present topological characterization of discrete complexes for each dimension from one to three by using the notions of star and link [START_REF] Ziegler | Lectures on Polytopes[END_REF] similarly to the previous work [START_REF] Kenmochi | Discrete combinatorial geometry[END_REF]. We then show that there are 12 topological types of points in Sk(C). In the next sections, we will classify all points in Sk(C) by their topological characteristics and study the type of points in discrete surfaces.

In this section, we do not have to distinguish the three different neighborhood systems. Thus, we abbreviate a discrete complex C m for m = 6, 18, 28 simply to C.

Star and link

The star and the link are defined for each point in Sk(C) as follows.

Definition 7 For a discrete complex C, the star of a point x ∈ Sk(C) is defined such that

star(x) = {σ ∈ C : x ∈ σ}.
Definition 8 For a discrete complex C, the link of a point x ∈ Sk(C) is defined such that link(x) = Cl(star(x)) \ star(x).

The star and link with respect to C are denoted by star(x : C) and link(x : C), respectively, when we emphasize C.

The star is defined with respect to a point x in Sk(C) as the set of discrete convex polyhedra including x. For example, let us consider the star of a point

x in Fig. 7 (a). It includes, as shown in Fig. 7 (b), a pentahedral 3-polyhedron, four triangle 2-polyhedra, four line-segment 1-polyhedra, and a 0-polyhedron which is the point x itself. The link of x in Fig. 7 (a) is obtained such as a squared 2-polyhedron and its faces as shown in Fig. 7 (c). Similarly to discrete complexes, we define the dimension of star(x : C) as the maximum dimension of discrete convex polyhedra belonging to star(x : C) and denoted by dim(star(x : C)). Note that a star is not always a discrete complex because it may not satisfy the second condition in Definition 3.

Topological characteristics of stars in one dimension

Stars of points in Sk(C) are classified into the following three types: linear stars, semi-linear stars, and neither of them, when the dimensions of stars are one.

Definition 9 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is linear if link(x) consists of two 0-polyhedra. Definition 10 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is semi-linear if link(x) consists of one 0-polyhedron. Figure 8 illustrates points whose stars are linear, semi-linear and neither of them. We see that a point is an endpoint of a curve if its star is semi-linear, and an intermediate point of a curve if its star is linear. If the star of a point is neither linear nor semi-linear, it is an intersection of a curve. By using the above definitions, we define discrete curves in Z 3 .

Definition 11 Let C be a connected and pure 1-complex. We say that C is a discrete curve with endpoint if the star of every point in Sk(C) is either linear or semi-linear and there is at least one point whose star is semi-linear in Sk(C).

Definition 12 Let C be a connected and pure 1-complex. We say that C is a discrete closed curve if the star of every point in Sk(C) is linear.

Topological characteristics of stars in two dimensions

Stars of points in Sk(C) are classified into the following three types: cyclic stars, semi-cyclic stars, and neither of them, when the dimensions of stars are two.

Definition 13 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is cyclic if link(x) is a discrete closed curve. Definition 14 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is semi-cyclic if link(x) is a discrete curve with endpoint. Figure 9 illustrates points of stars which are cyclic, semi-cyclic and neither of them. We see that a point is an edge point of a surface if its star is semi-cyclic, and an interior point of a surface if its star is cyclic. If the star of a point is neither cyclic nor semi-cyclic, the point is at an intersection of surfaces. By using these characteristics, we define discrete surfaces in Z 3 .

Definition 15 Let C be a connected and pure 2-complex. We say that C is a discrete surface with edge if every point in Sk(C) has either a cyclic or semicyclic star, and there is at least one point whose star is semi-cyclic in Sk(C).

Definition 16 Let C be a connected and pure 2-complex. We say that C is a discrete closed surface if every point in Sk(C) has a cyclic star.

Definition 15 (resp. [START_REF] Khalimsky | Computer graphics and connected topologies on finite ordered sets[END_REF]) is similar to the definition of a two-dimensional combinatorial manifold with (resp. without) boundary in reference [START_REF] Françon | Discrete combinatorial surfaces[END_REF]. In fact, the notion of cyclic stars in Definition 13 corresponds to that of umbrellas in [START_REF] Françon | Discrete combinatorial surfaces[END_REF].

Topological characteristics of stars in three dimensions

Stars of points in Sk(C) are classified into the following three types: spherical stars, semi-spherical stars, and neither of them, when the dimensions of stars are three.

Definition 17 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is spherical if link(x) is a discrete closed surface.

We define semi-spherical stars by using the notion of combinatorial boundary.

Definition 18 Let C be a pure n-complex and H be the set of all (n -1)polyhedra in C each of which is a face of exactly one n-polyhedron in C. The combinatorial boundary of C is then defined as a pure (n -1)-complex such (a) (b) Fig. 10. An example of a point whose star is considered to be semi-spherical with our previous definition in [START_REF] Kenmochi | Discrete combinatorial geometry[END_REF], but not semi-spherical with that in this paper, illustrated as the white point in (a). A halftone region and two black bold closed curves in (b) show its link and its combinatorial boundary.

that ∂C = Cl(H).

We now see that the endpoints of a discrete curve and the edges of a discrete surface in Definitions 11 and 15 are their combinatorial boundaries. An example is shown in Fig. 10 (b). If we consider the link of the central white point in Fig. 10 (a) illustrated as a half-tone region in Fig. 10 (b), its combinatorial boundary is shown as black bold lines.

Definition 19 Let C be a discrete complex and x be a point in Sk(C). We say that star(x) is semi-spherical if link(x) is a discrete surface with edge and the edge, i.e., the combinatorial boundary ∂(link(x)) is a discrete closed curve.

In the previous work [START_REF] Kenmochi | Discrete combinatorial geometry[END_REF], semi-spherical stars are simply defined such that link(x) is a discrete surface with edges. However, we found a counter example such as a white central point in Fig. 10 (a); its link is a discrete surface with two boundaries, and its star should not be regarded as a semi-spherical star because it is not topologically equivalent to a semi-sphere. We therefore modify our definition of semi-spherical stars.

Figure 11 illustrates points whose stars are spherical, semi-spherical and neither of them. It also shows that a point whose star is spherical is an interior point in a 3-complex, a point whose star is semi-spherical is a boundary point of a 3-complex, and a point whose star is neither spherical nor semi-spherical is a singular point, i.e., an intersection point of the boundaries.

We present the following proposition which plays an important role in this paper.

Proposition 20 Let C be a pure 3-complex and x be a point in Sk(C). If star(x : C) is semi-spherical, then star(x : ∂C) is cyclic. (3)

For any σ ∈ H, from Definition 18, there is exactly one 2-polyhedron τ ∈ link(x : C) such that σ ≺ τ . Therefore, for any σ ∈ H, there is exactly one 3-polyhedron γ ∈ star(x : C) such that σ ≺ γ as well. We then easily see that σ becomes the common 1-face of two 2-faces τ 1 , τ 2 of γ. Because we know from the above that there is only one 2-face in link(x : C) for σ, set to be τ 1 , another 2-face τ 2 is not in link(x : C) but in star(x : C). Such τ 2 is a 2-face of exactly one 3-polyhedron in star(x : C), that is γ. Therefore, τ 2 ∈ star(x : ∂C), and thus σ ∈ Cl(star(x : ∂C)) for any σ ∈ H. Because σ / ∈ star(x : ∂C), we obtain

Cl(H) = link(x : ∂C). (4)
From ( 3) and (4), link(x : ∂C) becomes a discrete closed curve and consequently, star(x : ∂C) is cyclic from Definition 13.

Twelve point types by topological star characteristics

By using topological characterization of stars described above, we classify points x ∈ Sk(C) into twelve types each of which satisfies one of the following conditions. Note that each point in Sk(C) is classified into one of them.

type 0: dim(star(x)) = 0; type 1a: star(x) is linear; type 1b: star(x) is semi-linear; type 1c: dim(star(x)) = 1 and star(x) is neither linear nor semi-linear; type 2a: star(x) is cyclic; type 2b: star(x) is semi-cyclic; type 2c: dim(star(x)) = 2, Cl(star(x)) is pure and star(x) is neither cyclic nor semi-cyclic; type 2d: dim(star(x)) = 2 and Cl(star(x)) is not pure; type 3a: star(x) is spherical; type 3b: star(x) is semi-spherical; type 3c: dim(star(x)) = 3, Cl(star(x)) is pure and star(x) is neither spherical nor semi-spherical; type 3d: dim(star(x)) = 3 and Cl(star(x)) is not pure. Some examples of points whose types are 2d and 3d are shown in Fig. 12.

Table 5

Classification results of points from P1 to P41 in Figs. 8, 9, 11, 12 (P1-P5 in Fig. 8, P6-P18 in Fig. 9, P19-P28 in Fig. 11, P29-P41 in Fig. 12

Topological classification of local point configurations

As we mentioned above, any point is classified into one of the twelve types. We also see that these twelve types have a hierarchical structure as shown in Fig. 13. By using the hierarchical structure, we can easily obtain an algorithm to clas- sify a local 1-point set V with respect to a type of the star of the central point

x ∈ V:

(1) construct a discrete complex C m for m = 6, 18 or 26 from V;

(2) obtain star(x) in C m ;

(3) classify star(x) by its dimension; (4) if the dimension is more than zero, classify star(x) by its topological characteristics (including the purity of Cl(star(x)) for more than one dimension).

Note that x is always a 1-point for any 1-point configuration of V. In addition, a discrete complex C m for m = 6, 18, 26 is uniquely obtained from V as described in Section 2 mostly except for the case of m = 18 as mentioned at the end of Section 2.3.2. In order to simplify the situation for m = 18 and to obtain a unique type of the central point x for any input V, we therefore assume that 1-points exist only in N 26 (x) and there are only 0-points in the exterior. Consequently, if we have a unit cube which is the last case in Table 3 and whose four coplanar 1-points are located at the boundary of N 26 (x) and shared by an adjacent unit cube in the exterior of N 26 (x), the cube decomposition is simply considered as the last one of Table 3 but not as that of Table 4 to construct C 18 .

We apply the algorithm to every local 1-point configurations of V ⊆ N 26 (x) whose central point x is 1-point. The number of all possible 1-point configurations of V is 2 26 = 67108864, that is reduced to 1426144 up to rotations around the x-, y-and z-axes and symmetries with respect to the xy-, yz-, zx-planes.

Among them, we count the number of each type of local point configurations.

Table 6 shows the result for each neighborhood system.

Obviously, because of (5), they have the same type 3b, and have the same forms around the central points.

For m = 6, 18, the following equation does not always hold;

Sk(C m ) = Sk(Cl(star(x : C m ))), (6)
while it always holds for m = 26. For example, we see in Fig. 14 that ( 6) holds for C 6 (a) but does not for C ′ 6 (b). In order to avoid counting the local point configurations twice for C 6 and C ′ 6 in Fig. 14, we count different configurations of Sk(Cl(star(x : C m ))) instead of those of Sk(C m ). Such configurations are called star configurations and we obtain Table 7. Note that star(x : C m ) is not a discrete complex so that we make a smallest complex by using the closure function before taking its skeleton. We also mention that we count configurations of star(x : C m ) up to rotations and symmetries; thus no redundant configuration is contained in Table 7.

Clearly, the number of star configurations for each type in Table 7 is smaller than that of local point configurations in Table 6 for m = 6, 18. We also see that the total numbers of different star configurations for m = 6, 18 are much smaller than that for m = 26. Now let us also consider the cube decomposition of Table 4 for m = 18 to obtain the more precise results for Table 7. Experimentally, if we consider all possible cube decompositions such as Tables 3 and4, the results for m = 18 in Table 7 become slightly bigger numbers for all the 2-dimensional stars and the types 3c and 3d: 88 for type 2a, 347 for type 2b, 1309 for type 2c, 558 for type 2d, 58890 for type 3c, and 184393 for type 3d. For the other stars including type 3b which is of our interest in this paper, the numbers of configurations become the same as those in Table 7. This is because, as far as the case of type 3b be concerned, influences on star configurations may be made by modification of polyhedral complexes in Fig. 4 when the central point x is considered to be the leftmost (or rightmost) point in Fig. 4, and such a modification seems not to change Sk(Cl(star(x : C 18 ))) even if it may change C 18 itself.

Local configurations in discrete combinatorial surfaces

Our discrete combinatorial surfaces appear at the 2-dimensional combinatorial boundaries of 3-complexes, that is, ∂C where dim(C) = 3. Because star(x : ∂C) is cyclic if star(x : C) is semi-spherical from Proposition 20, we see that semi-spherical stars whose point type is 3b give all local configurations appearing in such discrete combinatorial surfaces.

Let us consider a set of boundary points Br m ′ (V) of (1) for m ′ = 6, 18, 26. In [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF], we see that Br m ′ (V) includes all types of points except for spherical points (type 3a) which are interior points of V and that there are relations between our discrete complexes C m constructed from V and Br m ′ (V) only for the pairs (m, m ′ ) = (6, 18), (6, 26), [START_REF] Morgenthaler | Surfaces in three-dimensional digital images[END_REF][START_REF] Couprie | Simplicity surface: a new definition of surfaces in Z 3[END_REF], (26, 6) [START_REF] Kenmochi | Combinatorial boundary of a 3D lattice point set[END_REF]. Those relations indicate that boundary points of Br m ′ (V) do not have always semi-spherical stars, but also the other stars such as one-, two-and three-dimensional stars except for spherical stars depending on their local point configurations. By using our topological classification of local point configurations in the previous section, we easily discriminate semi-spherical stars from the other stars on boundaries.

As shown in Table 7, the numbers of semi-spherical star configurations of Sk(Cl(star(x : C m ))) are still large, especially for m = 18, 26. What we are interested in, however, is point configurations of stars of x in combinatorial boundaries ∂C m but not in C m . For example, we have different discrete complexes C m and C ′ m such that Sk(C m ) ⊃ Sk(C ′ m ) and star(x : ∂C m ) = star(x : ∂C ′ m ) as illustrated in Fig. 15. Note that this does not occur for m = 6. However, we have the case that C 6 = C ′ 6 and star(x : ∂C 6 ) = star(x : ∂C ′ 6 ) as shown in Fig. 16. We can also say that two discrete surfaces in Fig. 16 have the same shape but do not have the same orientation if we consider that they have two sides, the inside and outside.

In 8. For m = 18, we first calculate the result in Table 8 considering a unique complex for each local point set by using only Table 3 as mentioned before. If we consider all possible polyhedral complexes for each local point set, including the case of Fig. 4, we obtain 1720 instead of 1412 in Table 8. We verified that our results for m = 6 are the same as those in reference [START_REF] Françon | Discrete combinatorial surfaces[END_REF], which are illustrated in Fig. 2.

For m = 6, we see that surface star configurations are exactly the same as cyclic star configurations (type 2a in Table 7). However, for m = 18, 26, they contain more configurations than cyclic star configurations. One of the reasons is that a cyclic star requires that both interior and exterior points which are separated by a discrete surface exist in N 26 (x). For m = 18, 26, there are some surface star configurations where there is no interior point as illustrated in Fig. 1 (center). that of simplicity surfaces [START_REF] Couprie | Simplicity surface: a new definition of surfaces in Z 3[END_REF].

References [9,[START_REF] Kenmochi | Combinatorial topologies for discrete planes[END_REF] show that there are 5 and 32 different configurations of stars which appear in discrete combinatorial planes for m = 6, 26, respectively. Such planar stars for m = 6 are shown as the five left configurations in Fig. 2. We also illustrate the 32 configurations of planar stars for m = 26 in Fig. 17. Note that oriented surfaces are considered in [9,13] so that 8 and 34 configurations are obtained for m = 6, 26. While there is only one non-planar star for m = 6, we see that, for m = 26, most of the 6028, namely 5994 surface stars are non-planar and they do not appear on discrete planes but appear on discrete non-planar surfaces.

Fig. 3 .

 3 Fig. 3. Additional polyhedral decomposition for the 18-neighborhood system.

Definition 5 Fig. 5 .

 55 Fig. 5. (a) Non-pure and (b) pure 3-complexes.

Fig. 6 .

 6 Fig. 6. (a) A 3-complex C, (b) its subset B consisting of three convex polyhedra whose dimensions are from one to three respectively, and (c) the closure of B in C.

Fig. 7 .

 7 Fig. 7. (a) A 3-complex C; (b) the star of x ∈ Sk(C); (c) the link of x.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Examples of points whose stars are linear, semi-linear and neither of them, illustrated as white, grey and black points, respectively.P6

Fig. 11 .

 11 Fig. 11. Examples of points whose stars are spherical, semi-spherical and neither of them, illustrated as white, grey and black points, respectively.

Fig. 12 .

 12 Fig. 12. Examples of points whose types are 2d (a) and 3d (b), shown as black points, in non-pure discrete complexes.

Fig. 13 .

 13 Fig. 13. Hierarchical point classification by topological characterization of stars in Z 3 .

Fig. 14 .

 14 Fig. 14. Examples of two different 3-complexes C 6 (a) and C ′ 6 (b) whose central points, illustrated as white points, have type 3b such that Sk(C 6 ) = Sk(C ′ 6 ) but Sk(Cl(star(x : C 6 ))) = Sk(Cl(star(x : C ′ 6 ))).

  order to avoid counting twice for such two discrete complexes C m and C ′ m in Figs. 15, 16 respectively, we consider point configurations of Sk(Cl(star(x : ∂C m ))) instead of those of Sk(Cl(star(x : C m ))). Such point configurations are called surface star configurations. The results of counting different surface star configurations up to rotations and symmetries for m = 6, 18, 26 are shown in Table

Fig. 15 .

 15 Fig. 15. Two different discrete complexes C 26 (a) and C ′ 26 (b) around the central points x such that star(x : ∂C 26 ) = star(x : ∂C ′ 26 ).

Fig. 16 .

 16 Fig. 16. Two different discrete complexes C 6 (a) and C ′ 6 (b) around the central points x such that star(x : ∂C 6 ) = star(x : ∂C ′ 6 ).

6 Conclusions

 6 Given a subset V ⊆ N 26 (x), we presented a method for classifying the central point x into one of the twelve types by the topological characterization of its star after obtaining a complex C m . Considering that boundary points having type 3b form discrete combinatorial surfaces as Proposition 20, we counted local configurations in discrete surfaces such as local point configurations whose central point has type 3b, namely Sk(C m ), and obtained 14031, 345997 and 290979 configurations for m = 6, 18, 26, respectively, up to rotations and symmetries. We also obtained 9, 102793 and 290979 semi-spherical star configurations, namely, Sk(Cl(star(x : C m ), and 6, 1412 and 6028 surface star configurations, namely, Sk(Cl(star(x : ∂C m ))), for m = 6, 18, 26, respectively. For m = 18, if we consider all possible polyhedral complexes C 18 which may break the uniqueness of complex construction for certain local point sets V (but still keep the finiteness), we obtain 1720 surface star configurations. The same surface star configurations for m = 6 are already presented in reference[START_REF] Françon | Discrete combinatorial surfaces[END_REF] and they are illustrated in Fig.2. We see that a boundary point illustrated as the central point in Fig.1(center) has a surface star configuration. This explains why our discrete surfaces have more configurations than

Fig. 17

 17 Fig. 17. The 32 local star configurations which appear in discrete combinatorial planes for the 26-neighborhood system.

  Fig. 17. The 32 local star configurations which appear in discrete combinatorial planes for the 26-neighborhood system.

Figure 18 Fig. 18 .Fig. 19 .

 181819 Fig. 18. Boundary points of three-dimensional digitized objects, such as a cube, a sphere, a one-sheet hyperboloid and a hyperbolic paraboloid, are classified into two types in the 26-neighborhood system: they are illustrated as white and black points if the stars are planar and non-planar, respectively.

Table 1

 1 All possible discrete convex polyhedra for the 6-, 18-and 26-neighborhood systems.

Table 3

 3 Each unit cube in Table

  respectively) into the twelve types from 0 to 3d (no point for type 0).

	classified points
	type 1a 4
	type 1b 1, 3, 5, 33, 34
	type 1c 2
	type 2a 10
	type 2b 6, 7, 8, 9, 11, 12, 15, 16, 17, 18, 29, 30, 31, 40, 41
	type 2c 13, 14
	type 2d 32
	type 3a 22
	type 3b 19, 20, 21, 23, 25, 26, 27, 28, 36, 37
	type 3c 24
	type 3d 35, 38, 39

Table 5

 5 shows the classification results of all points from P1 to P41, illustrated in Figs. 8, 9, 11, 12, into the twelve types.

Table 6

 6 The numbers of different local point configurations in a 3 × 3 × 3 region for each of the twelve types, with respect to the 6-, 18-and 26-neighborhood system, up to rotations and symmetries.

		6-neighborhood 18-neighborhood 26-neighborhood
	type 0	23 520	22	1
	type 1a	200 712	520	11
	type 1b	134 280	143	3
	type 1c	103 092	735	77
	type 2a	26 862	2 852	55
	type 2b	345 016	18 741	398
	type 2c	113 008	19 324	3 203
	type 2d	399 329	18 434	3 664
	type 3a	1	23 520	23 520
	type 3b	14 031	345 997	290 979
	type 3c	374	131 073	321 371
	type 3d	65 919	864 783	782 862
	total	1 426 144	1 426 144	1 426 144

Table 7

 7 The numbers of star configurations for each of the twelve types, with respect to the 6-, 18-and 26-neighborhood system, up to rotations and symmetries.

		6-neighborhood 18-neighborhood 26-neighborhood
	type 0	1	1	1
	type 1a	2	6	11
	type 1b	1	2	3
	type 1c	6	17	77
	type 2a	6	80	55
	type 2b	14	313	398
	type 2c	123	938	3 203
	type 2d	74	461	3 664
	type 3a	1	21 425	23 520
	type 3b	9	102 793	290 979
	type 3c	11	58 532	321 371
	type 3d	274	179 893	782 862
	total	522	364 461	1 426 144
	but
	Sk(Cl(star(x : C 6 ))) = Sk(Cl(star(x : C ′ 6 ))).

Table 8

 8 The numbers of surface star configurations for the 6-, 18-and 26-neighborhood systems.

		6-neighborhood 18-neighborhood 26-neighborhood
	# of surface stars	6	1 412	6 028

Acknowledgements. A part of this work was supported by Grant-in-Aid for Scientific Research of the Ministry of Education, Culture, Sports, Science and Technology of Japan under the contract of 16650040 and was done under the research collaboration framework of National Institute of Informatics in Japan.