Antonin Gilles

Patrick Gioia
email: patrick.gioia@b-com.com.

Rémi Cozot
email: remi.cozot@b-com.com.

Luce Morin
email: luce.morin@b-com.com.

Complex modulation computer-generated hologram

Keywords: Computer-Generated Hologram, Color holography, Three-dimensional imaging

Introduction

Holography is often considered as the most promising 3D visualization technology, since it can provide the most authentic three-dimensional illusion to the naked eye [START_REF] Schnars | Digital Holography: Digital Hologram Recording, Numerical Reconstruction, and Related Techniques[END_REF]. Over the past decades, several methods have been proposed to generate Computer-Generated Holograms (CGH) of synthetic or existing scenes by simulating the propagation of light scattered by the scene towards the hologram plane. One of the most challenging feature in CGH techniques is the ability to properly take into account occlusions between objects in a scene. Occlusion processing in CGH techniques, also called light shielding, is the counterpart of hidden surface removal in Computer Graphics (CG). However, since CGH provides motion parallax, the occlusion effect in CGH implies that the visibility of objects changes according to the movement of the viewer. Two approaches are commonly used for CGH computation: the point-source and the wave-field approaches.

The point-source approach samples 3D scenes by a collection of self-luminous points, and calculates light propagation from the scene towards the hologram plane as the sum of spherical waves scattered by each point. Light shielding is usually treated as a visibility test to check the existence of obstacles between an object point and each sampling point in the hologram plane [START_REF] Underkoffler | Occlusion processing and smooth surface shading for fully computed synthetic holography[END_REF][START_REF] Zhang | Full parallax three-dimensional computer generated hologram with occlusion effect using ray casting technique[END_REF]. This approach is very flexible and does not impose any restriction on the scene geometry. However, its computational com-plexity is very high since it requires one calculation per point of the scene per pixel of the hologram. In order to reduce the computational complexity, several methods have been proposed, including pre-computed look-up tables [START_REF] Lucente | Interactive computation of holograms using a look-up table[END_REF][START_REF] Kim | Effective generation of digital holograms of threedimensional objects using a novel look-up table method[END_REF], difference and recurrence formulas [START_REF] Yoshikawa | Fast computation of Fresnel holograms employing difference[END_REF][START_REF] Matsushima | Recurrence Formulas for Fast Creation of Synthetic Three-Dimensional Holograms[END_REF], image holograms [START_REF] Yamaguchi | Real-time image plane full-color and full-parallax holographic video display system[END_REF], wavefront recording planes [START_REF] Shimobaba | Rapid calculation algorithm of Fresnel computergenerated-hologram using look-up table and wavefront-recording plane methods for three-dimensional display[END_REF][START_REF] Weng | Generation of real-time large computer generated hologram using wavefront recording method[END_REF], and using GPU hardware [START_REF] Ahrenberg | Computer generated holography using parallel commodity graphics hardware[END_REF][START_REF] Masuda | Computer generated holography using a graphics processing unit[END_REF] or special purpose hardware [START_REF] Ichihashi | HORN-6 special-purpose clustered computing system for electroholography[END_REF]. Additionally, it is also possible to reduce the number of visibility tests needed for light shielding by grouping sampling points in the hologram plane and/or grouping scene points [START_REF] Chen | Computer generated hologram from point cloud using graphics processor[END_REF][START_REF] Hanák | Acceleration of detail driven method for hologram generation[END_REF].

In the wave-field approach [START_REF] Lohmann | Three-dimensional properties of wave-fields[END_REF][START_REF] Bayraktar | Method to calculate the far field of three-dimensional objects for computer-generated holography[END_REF][START_REF] Oi | Electronic generation of holograms by using depth maps of real scenes[END_REF], 3D scenes are sliced into depth layers parallel to the hologram and light scattered by the scene is propagated from one layer to another using diffraction formulas such as the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. Light shielding is performed at each diffraction step by multiplying the light field by the binary silhouette mask of the current layer. The computation of the Angular Spectrum involves the use of the Fast Fourier Transform (FFT) algorithm twice, and is therefore more timeconsuming than the computation of the spherical light wave scattered by a single point. However, complex waves scattered by scene points located within a single layer are calculated all at once using the Angular Spectrum. Therefore, this approach is more efficient than the point-source approach when the number of scene points within each layer is sufficiently important. However, when the scene geometry contains complex shapes, a large number of layers containing only one or a few points are needed to sample it, reducing the efficiency of the wave-field approach. In order to reduce the computational complexity of this approach, a technique that uses color-space conversion has been proposed [START_REF] Shimobaba | Acceleration of color computergenerated hologram from three-dimensional scenes with texture and depth information[END_REF].

These two approaches complement one another: while the point-source approach is very efficient when the scene contains complex shapes, such as trees or human bodies, the wave-field approach is more efficient when objects in the scene consist of large planar surfaces, such as roads or buildings. However, most real 3D scenes contain both complex shapes and planar surfaces. As a consequence, these two approaches are rarely fully efficient for computing CGH of real scenes in their entirety. In order to overcome these limitations, we propose a new CGH computation method with occlusion effect based on a fast hybrid point-source/wave-field approach. Whereas previously proposed methods tried to reduce the computational complexity of the point-source or the wave-field approaches independently, the proposed method uses the two approaches together and therefore takes advantages from both of them. Section 2 gives an overview of the method, Section 3 and Section 4 present respectively light propagation and light shielding techniques used by the method, and Section 5 gives the detailed algorithm of the method. Finally, Section 6 gives experimental results.

Overview of the method

Figure 1 shows the scene geometry and coordinate system used by the proposed method. The coordinate system is defined by (x, y, z) so that the hologram lies on the (x, y, 0) plane. The 3D scene is sliced into a set of N z depth layers parallel to the hologram plane and located between z min and z max . The distance between each layer is given by d z . N z is set such that the separation between two consecutive layers remains invisible [START_REF] Bayraktar | Method to calculate the far field of three-dimensional objects for computer-generated holography[END_REF]. The layers are numbered from 0 to N z -1, from the farthest to the nearest to the hologram plane. The layers and the hologram are sampled on a regular 2D grid of resolution N × N with sampling pitch p.

The layers are classified into two categories: the wave-field layers and the point-source layers. This classification is performed depending on the number of scene points within each layer: if the number of scene points M d within layer d exceeds a threshold value M d,max , this layer is considered as a wavefield layer, otherwise it is considered as a point-source layer. The farthest and nearest layers are always considered as wave-field layers.

Each layer d operates as a surface source of light which emits a complex wave o d given by

o d (x, y) = A d (x, y) exp [jφ d (x, y)] , (1)
where A d (x, y) is the amplitude of the (x, y) point within layer d, calculated using illumination formulas, and φ d (x, y) is its phase, set to a random value in order to render a diffusive scene. Light scattered by the scene is numerically propagated and shielded from one layer to another using the recurrence formula

u 0 (x, y) = o 0 (x, y), u d (x, y) = o d (x, y) + O d {P dz {u d-1 (x, y)}} , (2)
where u d is the total complex wave scattered by layer d, operator O d stands for light shielding by layer d, and operator P z stands for the numerical propagation of light between two parallel planes separated by a distance z. These numerical operations are performed using either a point-source or a wave-field approach depending on which category layer d belongs to. Finally, the complex wave scattered by layer N z -1 is numerically propagated to the hologram plane in order to obtain the final CGH, according to

H(x, y) = P zmin {u Nz-1 (x, y)} . (3)
3 Light propagation

Wave-field layers

The propagation of complex waves scattered by wavefield layers is numerically computed using the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF], which expresses light propagation between two parallel planes separated by a distance z as

P w z {u d (x, y)} = F -1 F {u d (x, y)} × e -j2π √ λ -2 -f 2 x -f 2 y z , ( 4
)
where λ is the wavelength of light, f x and f y are the spatial frequencies, and F is the forward Fourier Transform, which can be computed using the Fast Fourier Transform algorithm (FFT). The calculation of the right-hand side member of Eq. ( 4) has for computational complexity O(N 2 log(N )).

Point-source layers

Scene points located within point-source layers are considered as self-luminous points. The complex wave scattered by a point source k at coordinates (x k , y k , z) is given by the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF] as

w k (x, y) = u d (x k , y k )F -1 e -j2π √ λ -2 -f 2 x -f 2 y z ⊗ δ(x -x k , y -y k ), (5)
where u d (x k , y k ) is the complex amplitude of the point and ⊗ is the convolution operator. In order to speed up the computation, we use a pre-calculated LUT, as proposed in [START_REF] Kim | Effective generation of digital holograms of threedimensional objects using a novel look-up table method[END_REF]. The LUT is pre-computed as

T (x, y, z) = F -1 e -j2π √ λ -2 -f 2 x -f 2 y z h(x, y, z), (6)
h being a window function used to restrict the region of contribution of a given point source, equal to one within the region of contribution of the point and zero elsewhere [START_REF] Lucente | Interactive computation of holograms using a look-up table[END_REF]. Then, light propagation between two parallel planes separated by a distance z is computed simply by addressing this pre-calculated LUT, as

P s z {u d (x, y)} = M -1 k=0 u d (x k , y k )T (x -x k , y -y k , z), (7)
where M is the number of self-luminous points within the source plane. The calculation of the right-hand side member of Eq. ( 7) has for computational complexity O(N 2 M ).

4 Light shielding

Wave-field layers

Light shielding by wave-field layers is performed using a binary silhouette mask function, as proposed in [START_REF] Lohmann | Three-dimensional properties of wave-fields[END_REF].

An occluding obstacle is located within layer d. The complex wave incident on layer d is given by

u ′ d (x, y) = P dz {u d-1 (x, y)} , (8)
where u d-1 is the complex wave scattered by layer d -1. Part of u ′ d is shielded by the occluding scene points, and its amplitude vanishes in the area of the obstacle. This is expressed by multiplying u ′ d by a binary mask function m d that has value zero inside the obstacle and one outside. Operator O w d is thus given by

O w d {u ′ d (x, y)} = m d (x, y)u ′ d (x, y). (9)
If the scene points within layer d also emit a wave o d , The complex wave incident on layer d + 1 is given by

u ′ d+1 (x, y) = P dz {o d (x, y) + m d (x, y)u ′ d (x, y)} (10) = P dz {u d (x, y)} .
While perfectly adapted to wave-field layers, this light shielding technique is not suited to point-source layers. Indeed, at each diffraction step, the complex wave scattered by the scene spreads gradually from one layer to another, and u ′ d may spread out on a large number of samples. Thus, even if layer d contains only a few scene points, u d may have a large number of non-zero values, which are therefore considered as a point sources by P s z . Since the computational complexity of P s z is dependent on the number of point sources within the source layer, this light shielding technique is highly inefficient for pointsource layers. In order to compute light shielding by point-source layers efficiently, we use a binary aperture function.

Point-source layers

Light shielding using a binary aperture function has been proposed by Matsushima et al. [START_REF] Matsushima | Silhouette method for hidden surface removal in computer holography and its acceleration using the switch-back technique[END_REF]. However, to the best of our knowledge, this technique has never been applied to point-source approach neither to CGH computation based on a layered model of the scene.

This time, we use an aperture function defined as

a d (x, y) = 1 -m d (x, y). (11)
By substituting Eq. ( 11) in Eq. ( 10), u ′ d+1 becomes

u ′ d+1 (x, y) =P dz o d (x, y) + (1 -a d (x, y))P dz {u d-1 (x, y)} =P 2dz u d-1 (x, y) + P dz o d (x, y) -a d (x, y)P dz {u d-1 (x, y)} =P 2dz {u d-1 (x, y)} + P dz {û d (x, y)} , (12)
where the propagation P dz {P dz {u d-1 (x, y)}} is equivalent to P 2dz {u d-1 (x, y)}, according to the definition of P z . The computation of u ′ d+1 using the binary aperture involves using P z three times, whereas we only need to use it twice using the binary mask. However, P dz {u d-1 (x, y)} needs to be calculated only within the region defined by a d , which corresponds to the coordinates of the scene points located within layer d. Therefore, the number of non-zero values of ûd remains the same as the number of scene points. Using the binary aperture function has the great advantage of not increasing the number of point sources for which the complex wave has to be computed at each diffraction step. Operator O s d is thus given by

O s d {u ′ d (x, y)} = -a d (x, y)u ′ d (x, y). ( 13
)
5 Detailed algorithm waves scattered by layers located between d and the last wave-field layer l. The propagation of these complex waves towards layer d is computed using P s z and shielded using O s d . If the scene points within layer d also emit a wave o d , the total complex wave scattered by this layer is given by

Point-source layers

u d (x, y) = o d (x, y) -a d (x, y) d-1 i=l P s (d-i)dz {u i (x, y)} . ( 14
)
Since a d has value one only at the coordinates of the scene points located within layer d, Eq. ( 14) needs to be calculated within this region only. The computational complexity of this step is thus O(M d M ), where M is the total number of scene points within all the layers located between d and

l: M = d-1 i=l M i .

Wave-field layers

The complex wave incident on each wave-field layer d is given by the sum of the complex wave scattered by the last wave-field layer l and complex waves scattered by point-source layers located between d and l. The propagation of the complex wave scattered by the last wave-field layer l is computed using P w z , according to

u w d (x, y) = P w (d-l)dz {u l (x, y)} . ( 15
)
The computational complexity of this step is O(N 2 log(N )). By contrast, the propagation of the complex waves scattered by point-source layers is computed using P s z , according to

u s d (x, y) = d-1 i=l+1 P s (d-i)dz {u i (x, y)} . (16)
The computational complexity of this step is O(N 2 M ′ ), where M ′ is the total number of scene points within the point-source layers located between

d and l: M ′ = d-1 i=l+1 M i .
If the scene points within layer d also emit a wave o d , the total complex wave u d scattered by this layer is given by

u d (x, y) = o d (x, y) + m d (x, y) (u w d (x, y) + u s d (x, y)) . (17)

Layers classification

The first step to implement the proposed method is to determine the value of M d,max . We call t s d the time needed to compute light propagation and light shielding for layer d using the point-source approach, d and the time needed to compute the propagation of complex wave scattered by layer d using P s z :

t s d = αM d M + βM d N 2 , ( 18
)
where α and β are constant coefficients. t w d is given by

t w d = γN 2 log(N ), (19)
where γ is a constant. We found the numerical values for the coefficients α, β and γ using the Gnuplot implementation of the nonlinear least-squares Levenberg-Marquardt algorithm [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF]: α = 6, 37.10 -8 , β = 1, 07.10 -7 , γ = 5, 35.10 -8 .

In order to maximize the efficiency of the proposed method, M d,max must be set such that

t s d = t w d ⇔ M d,max = γN 2 log(N ) αM + βN 2 . ( 21
)
6 Experimental results

The proposed method was implemented in C++/CUDA on a PC system employing an Intel Core i7-4930K CPU operating at 3.40 GHz, a main memory of 16 GB as well as three GPUs NVIDIA GeForce GTX 780Ti. For the experiments, we used three different test scenes (cf. Figure 3). The scenes are sliced into a set of N z depth layers parallel to hologram plane and located between z min = 0 and z max = 2cm. N z is set to 512, 1024 and 2048 for Dice1, Dice2 and City, respectively. The layers and the hologram are sampled on a regular 2D grid of resolution 4096 × 4096 with a sampling pitch p = 8.1µm. We compared the proposed method with GPU implementations of two other methods:

(1) the point-source method without occlusion effect proposed in [START_REF] Kim | Effective generation of digital holograms of threedimensional objects using a novel look-up table method[END_REF], and (2) the wave-field method with occlusion effect proposed in [START_REF] Lohmann | Three-dimensional properties of wave-fields[END_REF]. We adapted both methods to produce colorful complex modulation CGH.

Figure 4 shows the scene images numerically reconstructed from the CGH patterns of the three test scenes generated by the point-source method (first row), the wave-field method (second row), and the proposed method (third row). As seen in first row, the point-source method does not take into account occlusions in the scene, and the objects appear therefore semi-transparent, strongly limiting the realism of the displayed image. On the other hand, occlu- sions between objects in the scene are accurately reproduced using the wave-field and proposed methods without producing any visible artifact (rows 2 and 3).

As shown in Table 1, using the wave-field method, the average CGH calculation time per layer remains constant for the three test scenes. This means that the total calculation time increases linearly with the number of layers and is not dependent on the total number of scene points. By contrast, using the pointsource method, the average CGH calculation time per scene point remains constant for the three test scenes. This means that the total calculation time increases linearly with the number of scene points and is not dependent on the number of layers. As a result, the point-source method is less efficient than the wavefield method when the number of points per layer is high, as in scene Dice1, and is more efficient when the average number of points within each layer is low, as in scenes Dice2 and City.

As shown in Table 1, the total calculation time using the proposed method is dependent both on the number of layers and on the number of scene points. By combining these two approaches, the proposed method takes advantages from both of them and is therefore always more efficient. Using the proposed method, the total CGH calculation time has been decreased by 26.0%, 37.4% and 28.7% for Dice1, Dice2 and City, respectively, compared to the wave-field method, and by 31.2%, 9.89% and 25.5% for Dice1, Dice2 and City, respectively, compared to the point-source method. These experimental results confirm the performance superiority of the proposed method over the point-source and wave-field methods in terms of computation time.

Conclusion

In this paper, we proposed a fast CGH computation method with occlusion effect based on a hybrid point-source/wave-field approach. The 3D scene is sliced into several depth layers parallel to the hologram plane. Then, light scattered by the scene is propagated and shielded from one layer to another using either a point-source or a wave-field approach according to a threshold criterion on the number of points within the layer. Finally, we compute light propagation from the nearest layer to the hologram plane in order to obtain the final CGH.

Experimental results reveal that our method accurately takes into account occlusions between objects in the scene without producing any visible artifact. Furthermore, the CGH calculation time has been reduced up to 31, 2% and 37.4% compared to the pointsource and wave-field methods, respectively. This confirms the performance superiority of our method over the point-source and wave-field approaches in terms of computation time.

Figure 1 :

 1 Figure 1: Scene geometry and coordinate system used by the proposed method.

Figure 2 Figure 2 :

 22 Figure2shows the overall block-diagram of the proposed method. The complex wave incident on each point-source layer d is given by the sum of complex

Figure 3 :

 3 Figure 3: Computer graphics images of the three test scenes used for the experiments.

Figure 4 :

 4 Figure 4: Numerical reconstructions of the CGHs generated by the point-source method (first row), the wave-field method (second row), and the proposed (third row).

Table 1 :

 1 CGH computation times for the three test scenes using the wave-field, point-source and proposed methods.

	Scene	Number of layers	Number of scene points	Method	Total	Calculation time Per layer Per point
				Wave-field method	6.77s (100%)	13.22ms	0.53µs
	Dice1	512	12, 820, 048	Point-source method 7.28s (108%)	14.22ms	0.57µs
				Proposed method	5.01s (74.0%)	9.79ms	0.39µs
				Wave-field method	13.97s (100%) 13.64ms	0.83µs
	Dice2	1024	16, 736, 640	Point-source method 9.71s (69.5%)	9.48ms	0.58µs
				Proposed method	8.75s (62.6%)	8.54ms	0.52µs
				Wave-field method	26.88s (100%) 13.12ms	0.56µs
	City	2048	47, 648, 608	Point-source method 25.74s (95.8%) 12.57ms	0.54µs
				Proposed method	19.17s (71.3%) 9.36ms	0.40µs

* This work has been achieved within the Institute of Research and Technology b<>com, dedicated to digital technologies. It has been funded by the French government through the National Research Agency (ANR) Investment referenced ANR-A0-AIRT-07.


Antonin Gilles

Patrick Gioia
email: patrick.gioia@b-com.com.

Rémi Cozot
email: remi.cozot@b-com.com.

Luce Morin
email: luce.morin@b-com.com.

Fast generation of complex

Keywords: Computer-Generated Hologram, Color holography, Video holography, Three-dimensional imaging

come

INTRODUCTION

Holography is often considered as the most promising 3D visualization technology, since it can provide the most realistic and natural three-dimensional illusion to the naked eye. Indeed, it provides complete human depth cues without the need for special viewing devices and without causing eye-strain [START_REF] Schnars | Digital Holography: Digital Hologram Recording, Numerical Reconstruction, and Related Techniques[END_REF]. Holograms can be generated by computer calculation techniques. These techniques calculate Computer-Generated Holograms (CGH) of synthetic or existing scenes by simulating the propagation of light scattered by the scene towards the hologram plane. Thanks to its attractive visualization properties, CGH may have application in the field of videoconferencing or telepresence systems. However, because of the huge calculation burden, real-time CGH computation is still very challenging. In this paper, we propose a new method which aims at reducing the CGH calculation time for 3D video in order to get closer to real-time computation.

CGH computation techniques usually sample 3D scenes by a set of primitives and calculate light propagation as the sum of complex light waves scattered by each primitive. Two approaches for CGH computation are commonly used: the point-source approach, which samples scenes by a collection of points, and the wave-field approach, which samples scenes by a collection of planar segments.

In the point-source approach, 3D scenes are populated by a collection of self-luminous points. The complex waves scattered by the points are computed using the monochromatic spherical light wave equation and summed up together in order to get the wave scattered by the entire scene, as shown in Figure 1a. This approach presents the advantage of not imposing any restriction on the scene geometry. However, its computational complexity is very high since it requires one calculation per point of the scene per pixel of the hologram. Moreover, to produce shapes that appear solid and continuous, the scene needs to be sampled at very high densities, making the CGH computation prohibitively slow. In order to reduce the computational complexity, several methods have been proposed, including geometric symmetry [START_REF] Juárez-Pérez | Nonredundant calculations for creating digital Fresnel holograms[END_REF], difference and recurrence formulas [START_REF] Yoshikawa | Fast computation of Fresnel holograms employing difference[END_REF][START_REF] Matsushima | Recurrence Formulas for Fast Creation of Synthetic Three-Dimensional Holograms[END_REF], image holograms [START_REF] Yamaguchi | Real-time image plane full-color and full-parallax holographic video display system[END_REF][START_REF] Yoshikawa | Real-Time Generation of Full Color Image Hologram with Compact Distance Look-up Table[END_REF], wave-front recording planes [START_REF] Shimobaba | Simple and fast calculation algorithm for computer-generated hologram with wavefront recording plane[END_REF][START_REF] Shimobaba | Rapid calculation algorithm of Fresnel computer-generated-hologram using look-up table and wavefront-recording plane methods for three-dimensional display[END_REF][START_REF] Weng | Generation of real-time large computer generated hologram using wavefront recording method[END_REF], look-up tables [START_REF] Lucente | Interactive computation of holograms using a look-up table[END_REF][START_REF] Kim | Effective generation of digital holograms of three-dimensional objects using a novel look-up table method[END_REF], interframe and interline redundancy reduction [START_REF] Kim | Fast computation of hologram patterns of a 3d object using run-length encoding and novel look-up table methods[END_REF][START_REF] Kim | Accelerated computation of hologram patterns by use of interline redundancy of 3-D object images[END_REF], using GPU hardware [START_REF] Ahrenberg | Computer generated holography using parallel commodity graphics hardware[END_REF][START_REF] Masuda | Computer generated holography using a graphics processing unit[END_REF], and special purpose hardware [START_REF] Ito | Special-purpose computer HORN-5 for a real-time electroholography[END_REF][START_REF] Ichihashi | HORN-6 special-purpose clustered computing system for electroholography[END_REF].

In the wave-field approach, 3D scenes are populated by a collection of self-luminous planar segments. The complex waves scattered by the segments are computed using diffraction formulas such as the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF] and summed up together in order to get the wave scattered by the entire scene, as shown in Figure 1b. The Angular Spectrum can be calculated using the Fast Fourier Transform (FFT) algorithm twice. The computation of the complex wave scattered by each segment is therefore more time-consuming than the calculation of the spherical light wave scattered by a single point. However, complex waves scattered by scene points located within a single planar segment are calculated all at once using the Angular Spectrum. Therefore, this approach is more efficient than the point-source approach when objects in a scene consist of large planar segments containing many points. However, when the scene geometry contains complex shapes, a large number of small planar segments containing only one or a few points are needed to sample it, making the wave-field approach less efficient than the point-source approach. In order to reduce the CGH calculation time, several methods have been proposed, including the use of analytic expression of the angular spectrum [START_REF] Ahrenberg | Computer generated holograms from three dimensional meshes using an analytic light transport model[END_REF][START_REF] Kim | Mathematical modeling of triangle-mesh-modeled threedimensional surface objects for digital holography[END_REF][START_REF] Sakata | Fast computation method for a Fresnel hologram using threedimensional affine transformations in real space[END_REF][START_REF] Liu | A fast analytical algorithm for generating CGH of 3d scene[END_REF][START_REF] Liu | High-speed full analytical holographic computations for true-life scenes[END_REF] and color-space conversion [START_REF] Shimobaba | Acceleration of color computer-generated hologram from three-dimensional scenes with texture and depth information[END_REF][START_REF] Shimobaba | Calculation reduction method for color digital holography and computer-generated hologram using color space conversion[END_REF].

These two approaches complement one another: while the point-source approach is very efficient when the scene contains complex shapes, such as trees or human bodies, the wave-field approach is more efficient when objects in the scene consist of large planar surfaces, such as roads or buildings. However, most real 3D scenes contain both complex shapes and planar surfaces. As a consequence, these two approaches are rarely fully efficient for computing holograms of real scenes in their entirety. In order to overcome these limitations, a new hybrid point-source/wave-field approach has recently been proposed [START_REF] Gilles | Complex modulation Computer-Generated Hologram by a fast hybrid Point-source/Wave-field approach[END_REF]. In this approach, the scene is sliced into several depth layers parallel to the hologram plane. The complex wave scattered by each depth layer is then computed using either a wave-field or a point-source approach according to a threshold criterion on the number of points within the layer. Finally, the complex waves scattered by all the depth layers are summed up in order to obtain the final CGH. Although outperforming both point-source and wave-field methods without producing any visible artifact, this method requires the CGH to be entirely recomputed for each video frame, increasing the number of redundant calculations.

Indeed, there exist only slight changes between consecutive frames in 3D video, and the computation time of video holograms can be reduced by removing these temporal redundancies [START_REF] Plesniak | Incremental update of computer-generated holograms[END_REF][START_REF] Kim | Fast generation of three-dimensional video holograms by combined use of data compression and lookup table techniques[END_REF][START_REF] Kim | Fast generation of video holograms of three-dimensional moving objects using a motion compensation-based novel look-up table[END_REF][START_REF] Dong | MPEG-based novel look-up table for rapid generation of video holograms of fast-moving three-dimensional objects[END_REF]. In [START_REF] Plesniak | Incremental update of computer-generated holograms[END_REF], the author proposed an incremental method for making localized changes to a CGH in order to rapidly modify the scene image it reconstructs. At each simulation time step, the method tracks scene changes and updates the CGH patterns of affected scene points only. In [START_REF] Kim | Fast generation of three-dimensional video holograms by combined use of data compression and lookup table techniques[END_REF], the authors proposed to remove temporal redundancies in the input 3D video images using the differential pulse code modulation (DPCM) algorithm. Then, for each video frame, the CGH pattern of the changed part of the scene only is computed and added to the CGH of the previous frame in order to get the final CGH for the current frame. In [START_REF] Kim | Fast generation of video holograms of three-dimensional moving objects using a motion compensation-based novel look-up table[END_REF], the authors proposed a CGH method based on motion compensation technique. For each video frame, motion vectors of 3D moving objects are extracted and used in order to update the CGH pattern of the previous frame. In [START_REF] Dong | MPEG-based novel look-up table for rapid generation of video holograms of fast-moving three-dimensional objects[END_REF], the authors proposed a MPEG-based method for CGH computation. In this method, the input 3D video frames are grouped into sets of four, in which the first and remaining three frames become the reference (RF) and general frames (GFs). Then, each video frame is divided into blocks, from which motion vectors between the RF and the current frame are estimated. These motion vectors are then used in order to update the CGH pattern of the RF.

Since only affected regions of the CGH are updated at each video frame, these methods allow the computational burden to be dramatically reduced. However, all of these proposed methods are based on the point-source approach, and the possible exploitation of temporal redundancies using the hybrid point-source/wave-field approach has never been studied. Accordingly, in this paper, we propose a fast method for the computation of video holograms based on the hybrid approach by taking into account these temporal redundancies. Section 2 presents our method and Section 3 gives experimental results. Figure 2 shows the overall block-diagram of the proposed method. The input of the method is a 3D video consisting of views and depth maps pairs of resolution M x × M y . For each video frame i, a CGH H i of resolution N x × N y with pixel pitch p is computed.

PROPOSED METHOD

The CGH computation is done differently for the first video frame and for the following frames. For the first frame, H 0 is computed from the view and depth data using the hybrid point-source/wave-field approach. For the following frames, the view and depth data are first compared to those of the previous frame in order to track the differences. A criterion, defined in Section 2.6, determines whether these differences constitute an important part of the whole scene or not. When the changed part of the scene is important, H i is computed from the view and depth data of the current frame i using the hybrid point-source/wave-field approach. Otherwise, if only a little part of the scene has changed, the temporal redundancies can be removed in order to reduce the computational burden. Therefore, the method computes the CGH of the changed part of the scene only and sums it with the CGH of previous frame H i-1 in order to get H i . Finally, H i is transmitted to the CGH video output and stored in the previous frame buffer.

Scene geometry

Figure 3 shows the scene geometry and coordinate system used by the proposed method. The coordinate system is defined by (x, y, z) so that the hologram lies on the (x, y, 0) plane. For each video frame i, a 3D point cloud is extracted from the view and depth data, where each point is given an amplitude proportional to its corresponding pixel value in the view image. Since each depth map is encoded as an 8-bits gray level image, the extracted 3D point cloud is naturally sliced as a set of N z = 256 depth layers parallel to the hologram plane and regularly located between z min and z max .

The depth layers are classified into two categories: the wave-field layers and the point-source layers. This classification is performed depending on the number of scene points within each layer: if the number of scene points M i,d within layer d exceeds a threshold value M d,max (selection criterion which will be determined in Section 2.5), this layer is considered as a wave-field layer, otherwise it is considered as a point-source layer. Moreover, each depth layer d operates as a surface source of light which emits a complex wave o i,d given by

o i,d (x, y) = a i,d (x, y) exp [jφ d (x, y)] , (1)
where a i,d (x, y) is the amplitude of the (x, y) point within layer d for frame i, and φ d (x, y) is its phase, set to a random value in order to render a diffusive scene. Light scattered by each depth layer d is numerically propagated towards the hologram plane using either a point-source or a wave-field approach depending on which category layer d belongs to. Then, the complex waves scattered by all the depth layers are summed up in the hologram plane in order to obtain the final CGH H i , according to

H i (x, y) = H {o i,d } 0≤d<Nz = Nz-1 d=0 M i,d >M d,max P w z d {o i,d } (x, y) + Nz-1 d=0 M i,d ≤M d,max P s z d {o i,d } (x, y) = H w i (x, y) + H s i (x, y), (2)
where operator H stands for the numerical CGH computation using the hybrid approach, and operators P w z d and P s z d stand for the numerical propagation of light from layer d at depth z d to the hologram plane using the wave-field and the point-source approaches, respectively.

Wave-field layers

When the number of points M i,d within layer d exceeds M d,max , light scattered by this layer is numerically propagated towards the hologram plane using a wave-field approach. To this end, we use the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF], which expresses light propagation between two parallel planes separated by a distance z d by

P w z d {o i,d } (x, y) = F -1 F {o i,d } e -j2πz d √ λ -2 -f 2 x -f 2 y (x, y), (3)
where λ is the wavelength of light, f x and f y are the spatial frequencies, and F is the forward Fourier Transform.

In order to avoid one FFT per layer and therefore to speed-up the computation, the algorithm sums up the complex waves scattered by each wave-field layer directly in the frequency domain, and then inverse Fourier transforms the result to get H w i , as proposed in [START_REF] Ahrenberg | Computer generated holograms from three dimensional meshes using an analytic light transport model[END_REF]:

Ĥw i (f x , f y ) = Nz-1 d=0 M i,d >M d,max F {o i,d } e -j2πz d √ λ -2 -f 2 x -f 2 y ( 4
)
H w i (x, y) = F -1 Ĥw i (x, y). (5)

Point-source layers

When the number of points M i,d within layer d is less than M d,max , light scattered by this layer is numerically propagated towards the hologram plane using a point-source approach. Scene points located within the layer are therefore considered as isolated light sources, and light propagation is computed as the sum of complex waves scattered by each point. The complex wave scattered by a point k located at coordinates (x k , y k ) within layer d is given by the Angular Spectrum [START_REF] Goodman | Introduction to Fourier Optics[END_REF] as

w k (x, y) = o i,d (x k , y k )F -1 e -j2πz d √ λ -2 -f 2 x -f 2 y ⊗ δ(x -x k , y -y k ), (6)
where o i,d (x k , y k ) is the complex amplitude of the point and ⊗ is the convolution operator. Since convolving a function with a position-shifted Dirac delta shifts it by the same amount, knowing the inverse Fourier transform term in Eq. ( 6) beforehand allows us to compute w k simply by scaling this term with o i,d (x k , y k ), followed by a shifting operation. In order to speed up the computation, we use a pre-calculated LUT T , as proposed in [START_REF] Kim | Effective generation of digital holograms of three-dimensional objects using a novel look-up table method[END_REF]. T is pre-computed as

T (x, y, z) = F -1 e -j2πz √ λ -2 -f 2 x -f 2 y h(x, y, z), (7)
h being a window function used to restrict the region of contribution of a given point, equal to one within the region of contribution of the point and zero elsewhere. This function limits the spatial frequencies of the complex wave to avoid aliasing in the CGH.

According to the Nyquist Sampling Theorem, the maximum spatial frequency that can be represented with a sampling pitch p is f max = (2p) -1 . The relation between the maximum spatial frequency f max and the maximum diffraction angle θ is given by the grating equation [START_REF] Goodman | Introduction to Fourier Optics[END_REF] as sin(θ) = λf max . Therefore, the region of contribution of a point at depth z is given by its maximum radius R max (z) by

R max (z) = z tan(θ) = z tan arcsin λ 2p , (8)
as shown in Figure 4. The window function h can thus be defined as Thanks to its geometrical symmetry and to the window function h defined by Eq. ( 9), the LUT doesn't have to be computed for every (x, y) values. Instead, and in order to limit its number of samples, T is pre-computed only within the upper quarter of the square circumscribing the disk defined by h. Therefore, the number of samples N T,z d of the LUT for depth z d is given by

h(x, y, z) = 1 if x 2 + y 2 < R max (z) 0 otherwise. ( 9
)
N T,z d = R max (z d ) p 2 = z d p tan arcsin λ 2p 2 . (10)
Then, the complex wave scattered by point-source layer d in the hologram plane can be obtained by simply addressing this pre-calculated LUT, as

P s z d {o i,d } (x, y) = M i,d -1 k=0 o i,d (x k , y k )T (|x -x k |, |y -y k |, z d ). ( 11
)
The complex wave scattered by all the point-source layers in the hologram plane is therefore given by

H s i (x, y) = Nz-1 d=0 M i,d ≤M d,max M i,d -1 k=0 o i,d (x k , y k )T (|x -x k |, |y -y k |, z d ). ( 12
)

Layers classification

The first step for implementing the proposed method is determining the value of M d,max . We call t w the time needed to compute the complex wave scattered by a layer at depth z d with M i,d luminous points using operator P w z d , and t s the time needed to compute it using operator P s z d . Execution times t w and t s are obviously dependent on the implementation of P w z d and P s z d , respectively.

Since the computation of P w z d involves one complex multiplication per hologram pixel and a Fourier transform, t w is dependent on the number of hologram pixels N = N x × N y . In our implementation, the Fourier transform of Eq. ( 3) is performed on the GPU using the CUDA cuFFT library, so t w is dependent on the implementation of this library. In order to find the expression of t w , we measured the CGH computation time of one thousand random scenes containing only one wave-field layer for various hologram resolutions. t w is experimentally found to be linearly dependent on the number of hologram pixels N , such that

t w (N ) = kN. ( 13
)
On the other hand, the computation of P s z d involves one complex multiplication per sample of the LUT per point within the layer, so t s is dependent on the number of samples N T,z d of the LUT for depth z d and on the number of points M i,d within the layer. In order to find the expression of t s , we measured the CGH computation time of one thousand random scenes containing only one point-source layer for various number of scene points and various number of samples of the LUT. t s is experimentally found to be expressed by Figure 5 shows the graph of t s depending on the number of samples N T,z d of the LUT for three different values of M i,d . As shown in this figure, t s increases quickly with small values of N T,z d and then increases more slowly for higher values of N T,z d . This is due to the fact that when the number of samples N T,z d of the LUT is large, there is more chances for the regions of contribution of the different scene points to overlap in the hologram plane.

t s (M i,d , N T,z d ) = a (bN T,z d + c) 1 2 + dN T,z d + e M i,d . (14
We find the numerical values for the coefficients a, b, c, d, e and k in Eq. ( 13) and Eq. ( 14) using the Gnuplot implementation of the nonlinear least-squares Levenberg-Marquardt algorithm [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF]:

a = 2,
In order to maximize the efficiency of the proposed method, M d,max must be set such that

t s (M d,max , N T,z d ) = t w (N ) ( 16
)
⇔M d,max = kN a (bN T,z d + c) 1 2 + dN T,z d + e . ( 17
)

Incremental update

Similarly to 2D video, consecutive frames in the 3D video differ slightly from each other. This means that most of the scene points extracted from view and depth data remain the same from one frame to another. As a consequence, the CGH of the scene does not need to be entirely recomputed for each frame. Instead, the proposed method first tracks the scene changes between two consecutive frames in order to remove the temporal redundancies. Then, if only a little part of the scene has changed, the method computes the CGH of the changed part of the scene only and sums it with the CGH of the previous frame H i-1 in order to get H i . This allows the CGH computational burden to be dramatically reduced.

In order to track the differences between two consecutive video frames i -1 and i, the method computes the complex wave o ′ i,d scattered by the appearing and disappearing scene points within each depth layer d. We have

o ′ i,d (x, y) = o i,d (x, y) -o i-1,d (x, y) = (a i,d (x, y) -a i-1,d (x, y)) exp [jφ d (x, y)] , (18)
where a i,d (x, y) and a i-1,d (x, y) are the amplitude of the (x, y) point within layer d for frames i and i -1, respectively. The CGH of the changed parts of the scene is therefore given by

H ′ i (x, y) = H o ′ i,d 0≤d<Nz . (19)
The CGH H i can therefore be computed as

H i (x, y) = H o ′ i,d + o i-1,d 0≤d<Nz = H o ′ i,d 0≤d<Nz + H o i-1,d 0≤d<Nz = H ′ i (x, y) + H i-1 (x, y), (20)
where

H o ′ i,d + o i-1,d 0≤d<Nz = H o ′ i,d 0≤d<Nz + H o i-1,d 0≤d<Nz
, according to the definition of operator H.

Since there exist only slight changes in view and depth data between two consecutive frames of the 3D video, the number M ′ i,d of appearing and disappearing points between frames i -1 and i may be far less than the total number of scene points M i,d for frame i. Therefore, the computation of the CGH for the changed parts of the scene H ′ i may be faster than the computation of the CGH for the whole scene H i . The CGH computation time t i for frame i using Eq. ( 2) is given by

t i = Nz-1 d=0 M i,d >M d,max t w (N ) + Nz-1 d=0 M i,d ≤M d,max t s (M i,d , N T,z d ). ( 21
)
We call t i and t ′ i the computation times of H i and H ′ i using Eq. ( 2), respectively. The computation of H i using Eq. ( 23) is therefore faster than using Eq. (2) if

t ′ i < t i Nz-1 d=0 M ′ i,d >M d,max M d,max + Nz-1 d=0 M ′ i,d ≤M d,max M ′ i,d < Nz-1 d=0 M i,d >M d,max M d,max + Nz-1 d=0 M i,d ≤M d,max M i,d . (22)
For each video frame i, the CGH H i is therefore computed using the recurrence formula

H i (x, y) =    H o i,d 0≤d<Nz if i = 0 or t ′ i ≥ t i H o ′ i,d 0≤d<Nz + H i-1 (x, y) otherwise. ( 23
)

EXPERIMENTAL RESULTS AND DISCUSSION

The proposed method was implemented in C++/CUDA on a PC system employing an Intel Core i7-4930K CPU operating at 3.40 GHz, a main memory of 16 GB and an operating system of Microsoft Windows 8 as well as three GPUs NVIDIA GeForce GTX 780Ti.

Input 3D video

For the experiments, we generated with Unity an input 3D video containing 375 views and depth map pairs of resolution 512 × 512. From each view and depth map pair, a 3D point cloud is extracted, where each point is given an amplitude proportional to its corresponding pixel value in the view image. The extracted 3D point cloud is considered to be located between z min = 0 and z max = 2cm in front of the CGH plane. Finally, the CGH to be computed has a resolution of 2048 × 2048 with a sampling pitch p = 8, 1µm. : Scene points whose complex waves are updated at each frame. The scene points whose complex waves are computed by the proposed method using the wave-field approach are shown in blue. The scene points whose complex waves are computed using the point-source approach are shown in yellow.

Figure 6 shows the view and depth maps pairs for nine different frames of the video, which is composed of three parts. The first part (1 st to 125 th frames) shows a car turning around a building in a street. The second part (126 th to 250 th frames) shows the same scene with some cheering people. Finally, the third part (251 st to 375 th frames) shows the same scene from the driver's viewpoint. In the first and second parts of the video the camera is fixed, whereas it is moving with the car in the third part of the video. As a consequence, while consecutive frames in the first and second parts of the video differ slightly from each other, there exist large changes between consecutive frames in the third part of the video.

Reconstructed scene images

Figure 7 shows the scene points whose complex waves are updated at each frame. The scene points whose complex waves are computed using the wave-field approach are shown in blue, whereas those whose complex waves are computed using the point-source approach are shown in yellow. As seen in this figure, since there exist slight changes between consecutive frames in the first and second parts of the video, the proposed method updates the complex waves scattered by only a few scene points at each frame. On the other hand, since a large part of the scene changes between consecutive frames in the third part of the video, the proposed method recomputes the complex wave scattered by the whole scene at each frame. Moreover, due to the sudden changes of the scene between the first and second parts of the video, the proposed method also recomputes the complex wave scattered by the whole scene for the 126 th frame. Figure 8 shows the scene images numerically reconstructed from the CGH patterns of nine different frames. For the first and second parts of the video, the reconstructed images are focused on the facade of the building, and for the third part of the video, the reconstructed images are focused on the steering wheel. As shown in this figure, the proposed method does not produce any visible artifact, even at the boundaries between the two categories of scene points. In order to evaluate the objective quality of the reconstructed images compared to the hybrid method proposed in [START_REF] Gilles | Complex modulation Computer-Generated Hologram by a fast hybrid Point-source/Wave-field approach[END_REF], we used the Peak Signal-to-Noise Ratio (PSNR). The PSNR of the reconstructed images were found to be higher than 50dB for each video frame. These experimental results show that the CGH pattern is accurately updated at each video frame and that the proposed method does not reduce the quality of the reconstructed scene images compared to the hybrid method. Figure 9: Performance comparison between the wave-field method (in blue), the point-source method (in red), the hybrid method (in purple), the incremental method (in orange), and the proposed method (in green): (a) number of points whose complex waves are computed at each frame, (b) CGH computation time per frame.

In order to evaluate the performance of the proposed method in terms of computation time, we compare it with GPU implementations of four other methods: (1) the method proposed in [START_REF] Bayraktar | Method to calculate the far field of three-dimensional objects for computer-generated holography[END_REF], based on a wave-field approach, (2) the method proposed in [START_REF] Kim | Effective generation of digital holograms of three-dimensional objects using a novel look-up table method[END_REF], based on a point-source approach, (3) the method proposed in [START_REF] Kim | Fast generation of three-dimensional video holograms by combined use of data compression and lookup table techniques[END_REF], based on a point-source approach and incremental update of the CGH at each frame, and (4) the method proposed in [START_REF] Gilles | Complex modulation Computer-Generated Hologram by a fast hybrid Point-source/Wave-field approach[END_REF], based on a hybrid approach. In the following, these methods are called wave-field method, point-source method, incremental method and hybrid method, respectively. We adapted the first three methods in order to produce colorful complex modulation CGH.

Figure 9a shows the number of scene points whose complex waves are computed at each frame for all the compared methods. Since the wave-field, point-source and hybrid methods do not take into account temporal redundancies between consecutive frames in the 3D video, the number of calculated scene points remains the same for each video frame. By contrast, the incremental and proposed methods take into account those redundancies in order to reduce the number of calculated scene points. As a consequence, since consecutive frames in the first and second parts of the video differ slightly from each other, the number of calculated points is much lower for the incremental and proposed methods than for the wave-field, point-source and hybrid methods. Due to the cheering people appearing in the second part of the video, the number of calculated scene points is higher for these frames than for the first part of the video. Finally, since there exist large changes between consecutive frames in the third part of the video, the number of calculated scene points is exactly the same for all the compared methods.

Figure 9b shows the CGH computation time per frame for all the compared methods. The computational complexity of the wave-field method is dependent on the number of depth layers containing scene points. By contrast, the computational complexity of the point-source method is dependent on the number of scene points and on their location within the scene. Since consecutive frames in the first and second parts of the video differ slightly from each other, the number of non-empty depth layers and the location of scene points are almost the same from one frame to another. As a consequence, the CGH computation times of the wave-field and point-source methods remain almost constant for the first and second parts of the video, as shown in Figure 9b. On the contrary, since there exist large changes between consecutive frames in the third part of the video, the number of non-empty depth layers and the location of scene points change from one frame to another. Therefore, the CGH computation times of the wave-field and point-source methods also change from one frame to another.

Since it combines the wave-field and point-source approaches, the computational complexity of the hybrid method is dependent on the number of depth layers containing scene points, on the number of scene points and on their location within the scene. As a consequence, its CGH computation time varies in the same way as those of the wave-field and point-source methods, as shown in Figure 9b. However, by combining these two approaches, the hybrid method takes advantages from both of them and is therefore always more efficient than the wave-field and point-source methods.

In the first and second parts of the video, the number of calculated points of the incremental and proposed methods is much lower than that of the wave-field, point-source and hybrid methods. As a consequence, the incremental and proposed methods take much less time to calculate the CGH patterns than the wave-field, point-source and hybrid methods, as shown in Figure 9b. While the CGH computation times of the incremental and proposed methods are almost the same for the first part of the video, the proposed method is more efficient than the incremental method for the second part of the video. This is due to the fact that when the number of calculated scene points increases, the efficiency of the hybrid approach also increases [START_REF] Gilles | Complex modulation Computer-Generated Hologram by a fast hybrid Point-source/Wave-field approach[END_REF]. Finally, in the third part of the video, the number of calculated scene points is exactly the same for all the compared methods, and the CGH computation times of the incremental and proposed methods reach those of the point-source and hybrid methods, respectively. As a consequence, the proposed method is much more efficient than the incremental method for the third part of the video. 1: Comparison of the average number of calculated scene points and the average CGH computation time per frame between the hybrid, the incremental and the proposed methods.

Video part

Table 1 shows the average number of calculated scene points and the average CGH computation time per frame for the hybrid, incremental and proposed methods. As shown in this table, for the first part of the video, the average number of calculated scene points per frame is equal to 215293 for the hybrid method and to 4282 for the incremental and proposed methods. The average number of calculated scene points has therefore been reduced by 98, 01% using the incremental and proposed methods. As a consequence, the average CGH computation time per frame is equal to 174, 46ms, 16, 69ms and 16, 60ms for the hybrid, incremental and proposed methods, respectively. The CGH computation time per frame has therefore been reduced by 90, 43% and 90, 48% using the incremental and proposed methods, respectively.

Similarly, for the second part of the video, the average number of calculated scene points per frame is equal to 215293 for the hybrid method and to 32537 for the incremental and proposed methods. The average number of calculated scene points has therefore been reduced by 84, 9% using the incremental and proposed methods. As a consequence, the average CGH computation time per frame is equal to 188, 02ms, 96, 15ms and 75, 37ms for the hybrid, incremental and proposed methods, respectively. The CGH computation time per frame has therefore been reduced by 48, 9% and 59, 9% using the incremental and proposed methods, respectively. This confirms that the proposed method is more efficient than the incremental method for the second part of the video.

Finally, for the third part of the video, the average number of calculated scene points per frame is equal to 215296 for the hybrid, incremental and proposed methods. As a consequence, the average CGH computation time per frame is equal to 164, 57ms, 492, 80ms and 168, 07ms for the hybrid, incremental and proposed methods, respectively. The proposed method is therefore much more efficient than the incremental method when there exist important changes between consecutive frames of the input 3D video. The CGH computation time per frame has been increased by 2% using the proposed method compared to the hybrid method. This is due to the additional pre-processing stage used in order to track the differences from one frame to another.

These experimental results show that for the 3D video cases in which consecutive frames differ slightly from each other, the number of scene points to be calculated and the CGH computation time can be significantly reduced using the proposed method. Additionally, for the 3D video cases in which there exist large changes between consecutive frames, the proposed method is still more efficient than the methods based on point-source or wave-field approaches.

CONCLUSION

In this paper, we proposed a fast Computer-Generated Hologram (CGH) computation method based on a hybrid point-source/wave-field approach and temporal redundancy compression. Our algorithm consists of three steps. First, intensity and depth data of the current 3D video frame are extracted and compared with those of the previous frame in order to remove temporally redundant data. Then, the CGH pattern for this compressed frame is generated using the hybrid point-source/wave-field approach. The resulting CGH pattern is finally transmitted to the video output and stored in the previous frame buffer. Experimental results with three kinds of input 3D videos revealed that our method is able to produce colorful video holograms at interactive rates without producing any visible artifact.

Figure 1 :

 1 Figure 1: Two approaches are commonly used for CGH computation: (a) the point-source approach, which samples 3D scenes by a collection of self-luminous points and calculate light propagation as the sum of spherical waves scattered by each point, and (b) the wave-field approach, which samples 3D scenes by a collection of self-luminous planar segments and calculate light propagation using diffraction formulas.

Figure 2 :

 2 Figure 2: Overview of the proposed method.

Figure 3 :

 3 Figure 3: Scene geometry and coordinate system used by the proposed method.

Figure 4 :

 4 Figure 4: Region of contribution of a given point source.

  time (in ms) Number of pixels of the LUT N T, z d M d = 10k M d = 20k M d = 30k

Figure 5 :

 5 Figure 5: Time needed to compute the complex wave scattered by a point-source layer depending on the number of samples of the LUT and on the number of scene points.

  69.10 -5 b = 0, 34 c = 1, 11 d = 8, 49.10 -7 e = 2, 76.10 -4 k = 5, 54.10 -7

(a) 1 Figure 6 :

 16 Figure 6: View and depth maps pairs for nine different frames of the input 3D video.

(a) 1

 1 st frame (b) 31 st frame (c) 51 st frame (d) 126 th frame (e) 156 th frame (f) 176 th frame (g) 251 st frame (h) 281 st frame (i) 301 st frame

Figure 7

 7 Figure 7: Scene points whose complex waves are updated at each frame. The scene points whose complex waves are computed by the proposed method using the wave-field approach are shown in blue. The scene points whose complex waves are computed using the point-source approach are shown in yellow.

(a) 1 Figure 8 :

 18 Figure 8: Scene images numerically reconstructed from the CGH patterns of nine different frames.

2.1 Overview Current rame C anges rev ous rame - Computat on o t e C H or t e w ole scene Computat on o t e C H or t e c ange part o t e scene v eo utput u er pirst frame or ig changes ew changes

* This work has been achieved within the Institute of Research and Technology b<>com, dedicated to digital technologies. It has been funded by the French government through the National Research Agency (ANR) Investment referenced ANR-A0-AIRT-07.


Pedram Masoudi
email: pedram.masoudi@univ-rennes1.fr

Babak Nadjar Araabi
email: araabi@ut.ac.ir

Tahar Aïfa
email: tahar.aifa@univ-rennes1.fr

Hossein Memarian
email: memarian@ut.ac.ir

Clustering as an efficient tool for assessing fluid content and movability by resistivity logs

Keywords: water saturation in carbonate reservoirs, clustering permeable zones, clustering saturated zone, application of unsupervised learning, Gustafson-Kessel clustering

Neither core measurements nor well tests provide precise measurement of fluid contents; so, at the moment there is no validated saturation measurement in the oil industry. However, resistivity logs contain some valuable information about reservoir fluids, and classes of saturation based on the behaviour of resistivity data. In order to extract saturation value, clustering algorithms are proposed and tested here. Clustering is an unsupervised categorization method, which relies on natural groupings of real data, instead of predefined labels. Distance or dissimilarity plays an important role in the formation of clusters in an algorithm. The application of three clustering algorithms on prediction of water saturation is discussed here. It is shown that Fuzzy C-Means clustering divides data only according to saturation property; while the Gustafson-Kessel algorithm considers not only saturation but also permeability. The reason is that Gustafson-Kessel can detect linear patterns. This algorithm is introduced as the most appropriate clustering method for predicting permeability, and for understanding the reservoir quality of the formation under examination. Gath-Geva clustering did not provide as much information as Gustafson-Kessel. In addition to these three clustering algorithms, two other methods were likewise checked, but they did not provide acceptable interpretations, so not reported. Another achievement of this study is the introduction of a cluster label instead of a single value which is very unreliable for expressing saturation and permeability simultaneously. Predictions of these two petrophysical properties are provided merely by two conventional resistivity welllogs: deep and shallow. The output of cluster analysis is much closer to well-scale reservoir properties as compared with that of core-scale properties, since clustering algorithms are applied to logs that are volumetric recordings. The applicability and efficiency of the proposed methods are examined and verified through the use of well-logs of Sarvak Formation in an anticlinal oil field in the Abadan Plain, Iran.

INTRODUCTION

The importance of evaluating water saturation in petroleum exploration is taken for granted by petroleum geoscientists. The widest usage of water saturation is in net pay detection. Petrophysicists apply cut-off values to water saturation logs, to distinguish oil-bearing intervals in both horizontally [START_REF] Mostafazadeh | The Productivity Estimation of Designed Horizontal Oil and Gas Wells Before a Drilling Operation, Using Seismic and Petrophysical Parameters and Modeling[END_REF] and vertically drilled wells [2, 3, 4 and 5]. Due to the effect of water saturation on the plasto-elastic properties of rocks, accurate determination of this property is essential in geomechanics [START_REF] Amalokwu | Water saturation effects on elastic wave attenuation in porous rocks with aligned fractures[END_REF]. In addition, [START_REF] Tokhmechi | Fracture detection from water saturation log data using a Fourier-wavelet approach[END_REF] have used water saturation log in the process of fracture detection within wells [START_REF] Tokhmechi | Fracture detection from water saturation log data using a Fourier-wavelet approach[END_REF]. Production planning is another criterion which water saturation is essential for [START_REF] Lopes | Forced imbibition into a limestone: measuring P-wave velocity and water saturation dependence on injection rate[END_REF].

Still, empirical relations for estimating water saturation are the most widely used methodologies. Among them, Archie relations, developed by early Archie works [9, 10, 11 and 12] are the best-known. There are so many other empirical relations which have also been developed, such as Carman [START_REF] Carman | Fluid flow through granular beds[END_REF] and Timur [START_REF] Timur | An Investigation Of Permeability, Porosity, and Residual Water Saturation Relationships For Sandstone Reservoirs[END_REF]. For the comprehensive study of empirical methods in estimating water saturation, please refer to [START_REF] Tiab | Petrophysics: theory and practice of measuring reservoir rock and fluid transport properties[END_REF]. [START_REF] Worthington | Improved Quantification of Fit-for-Purpose Saturation Exponents[END_REF] offered a wise fitfor-purpose usage of Archie and non-Archie relations due to rock type and our purpose of estimating water saturation [START_REF] Worthington | Improved Quantification of Fit-for-Purpose Saturation Exponents[END_REF].

There are some other methodologies, based on supervised methods and trainable machines, which have recently been developed. To use these supervised methods, a reference criterion is required to train a machine. For example, Kadkhodaie-Ilkhchi et al. (2009) proposed a Committee Machine with Training Algorithm (CMTA), utilizing a combination of Levenberg-Marquardt, Bayesian regularization, gradient descent, one step secant and resilient back-propagation to train a neural network for predicting Normalized Oil Content (NOC) [START_REF] Kadkhodaie-Ilkhchi | A committee neural network for prediction of normalized oil content from well log data: An example from South Pars Gas Field, Persian Gulf[END_REF].

The Bayesian classifier can be applied, while using Gassmann equation as labels [START_REF] Mollajan | Zonal classification by pattern recognition methods: An example from Asmari Formation (Mansuri oil field, south of Iran)[END_REF]. Also, Radial Basis Function Neural Networks [START_REF] Mollajan | Estimation of water saturation from petrophysical logs using radial basis function neural network[END_REF] and support vector regression machine [START_REF] Mollajan | Prediction of Reservoir Water Saturation Using Support Vector Regression in an Iranian Carbonate Reservoir[END_REF] could be trained according to core measurements on saturation as reference criteria. It seems that these trainable (supervised) methodologies are more reliable as we have evidence for measuring precision of output. [START_REF] Barros | Determination of water saturation by angular competitive neural network[END_REF] predicted water saturation by angular competitive neural network. In this methodology there is no training, and the algorithm only searches for patterns of cross-plots to predict a saturation parameter. This algorithm is really applicable and useful in cases in which we lack labelled data as validation [START_REF] Barros | Determination of water saturation by angular competitive neural network[END_REF]. We consider the current manuscript as complementary to [START_REF] Barros | Determination of water saturation by angular competitive neural network[END_REF], since here, the proposed approach searches for patterns without considering labelled data for machine learning. In this study, a review is addressed on previous works regarding water saturation. Then, some clusteringmethodologies are introduced as possible solutions to the problem of saturation prediction. Finally, outputs of different cluster analysis methods in studying water saturation are discussed.

As mentioned above, the industrial importance of this work is its applicability in cases where there is no evidence for validating outputs. Still, there is no reliable continuous saturation measurement within drilled wells. It is noteworthy that the experience and wisdom of petrophysicists is essential for a successful interpretation, and cluster-based methodology is only a powerful tool in distinguishing similar data due to different criteria of similarity.

HYPOTHESES: WHY CLUSTER ANALYSIS FOR SATURATION PREDICTION?

As discussed previously, there are several saturation prediction methods. All of them provide an indirect procedure for predicting water content in each horizon, using well-log data. Nearly all saturation prediction methods on well-log data assign a single value (not a range) of water content to each horizon through wells. We believe that such a fixed saturation value is not a credible way of introducing this variable in oil wells due to three issues: (i) A fixed and single value is not an appropriate representative of the volume of investigation regarding recorded measurements within the well. In other words, the certainty and accuracy of a single value for expressing saturation of a volume is subject to doubt. Experts are keen to have a single value as the value of saturation; but, here a range of values, instead of a single value for saturation is recommended as to render the predictions more reliable. When we have vagueness in the parameter under study, we have to use a vague language too. In other words, at the moment, science and technology do not let us express water saturation by a fixed value. It is just an approximation without any error measurement! (ii) If we can train a machine in a specific environment, in this case a well which is drilled in a reservoir, it is more accurate to use the trained machine in evaluating water saturation instead of empirical universal equations. So, the use of empirical relations such as Archie is not recommended when we are able to use trainable and locally available machines. In addition, empirical relations have only been developed in specific environments, assumptions and considerations. For example, the Archie equation was primarily developed in sandy reservoirs, but if it is addressed in shale-contaminated reservoirs, some corrections have to be considered. Therefore, the Archie relation in carbonate reservoirs is only a rough approximation of water saturation.

(iii) Using supervised methods in saturation prediction is not accurate enough at the moment, because we do not have any reliable evidence or label for training the machine. Evidences for training have to be either core measurements or well-tests. Training a machine based on the label of water saturation is not an accurate job due to the following points:

(a) Whereas each core plug could be tested several times at several intervals (with a vertical resolution of about centimeter) for porosity measurements; it could only be used once for measuring water saturation. Therefore it is impossible through cores to get real, accurate data for saturation with high vertical resolution. Also, measuring water saturation in laboratory conditions results in different outputs, compared to the in-situ situation. Water saturation is too dependent on the place of analysis (because of environmental causes such as temperature, pressure, etc.); i.e. temperature and pressure of reservoir conditions differ greatly from laboratory conditions. As we recover the core from the reservoir, its pressure would be dropped, and the balance between fluid contents would be violated. Even in measuring porosity, in-situ and laboratory measurements differ somewhat due to stress relief, but this variation boosts in measuring saturation property, since dynamic processes take part.

(b) Surprisingly, well-tests cannot provide a practical label of saturation. Two issues are under doubt in well-test results as labels of supervised classification of water saturation: (b1) the very first issue is that wells are tested at intervals, not horizons. Therefore, when we do not have more than one test in a specific reservoir interval through a well, we would not have more than an average value for that specific interval, which is not enough for training a machine. In order to be able to train a machine in the highly heterogeneous conditions of the reservoirs, we require much more than 100 observations, while by well testing; we always obtain less than 10 observations, i.e. well test intervals. (b2) Besides, the nature of well-test results differs from the characters of water saturation, i.e. a welltest result is a combination of water saturation and reservoir permeability. The effect of permeability has to be removed from the testing results if we want to use the tests as a measurement of saturation. This incompatibility between well-test output and saturation value is shown in Figure 1. As it is obvious, well-test class numbers 2 and 3 (i.e. oil production lower and higher than 1500 barrel per day) show the reverse order of the expected pattern, compared to water saturation value [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]. Therefore, in predicting water saturation, we face a problem: there is not enough credible measurement for saturation. To cope with this situation, we believe that the assessment has to be done by an unsupervised method able to analyze the behavior of datasets according to the target (here, saturation). So, cluster analysis methods are recommended and investigated here as unsupervised clustering methods to evaluate the water saturation of reservoir rocks.

DATASETS

Datasets of this work belong to the Sarvak Formation in an anticlinal (north-south trend) oil field in the Abadan Plain, Iran. For the sake of confidentiality, the name of the field is not disclosed. There are six drilled wells for exploratory purposes, and well-log data is available through all the wells. But well-test and core data are not available at all. The summary of the available information is presented in Table 1.

From the standpoint of sequence stratigraphy, the Sarvak Formation is a shallowing upward carbonate reservoir, deposited on a shallow Upper Albian to Upper Turonian platform. It is conformably laid on the Kazhdumi Formation, and the Laffan Formation is overlaid on the Sarvak Formation above this sharp erosional unconformity [START_REF] Ghabeishavi | Microfacies and depositional environment of the Cenomanian of the Bangestan anticline, SW Iran[END_REF].

METHODOLOGY

A. Archie Equation

The famous empirical Archie equation was first developed in sandy reservoirs; however, it is used in carbonate rocks too. Two general forms of this relation are [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]:

The famous empirical Archie equation was first developed in sandy reservoirs; however, it is used in carbonate rocks too. Two general forms of this relation are [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]:

𝑅 𝑡 = 𝑅 𝑜 × 𝑆 𝑤 -𝑛 (1) 𝑆 𝑤 = ( 𝑅 𝑜 𝑅 𝑡 ) 1 𝑛 (2)
where, 𝑅 𝑡 is the recorded resistivity of rock, filled with fluids at reservoir conditions. 𝑅 𝑜 is the resistivity of sand when all the pores are filled with brine; 𝑆 𝑤 is the water saturation, in fraction form; and n is the saturation exponent. The saturation exponent for consolidated sands and clean unconsolidated sand appears to be close to 2. Hence, the equation for clean sand could be simply rewritten as:

𝑆 𝑤 = √ 𝑅 𝑜 𝑅 𝑡 (3)
𝑅 𝑡 is the resistivity measured by the logging instrument at each depth. Therefore, by having 𝑅 𝑜 , we can easily estimate water saturation by relation 3. Here, 𝑅 𝑜 is considered as the minimum value of 𝑅 𝑡 through the well. The logic behind choosing the minimum of 𝑅 𝑡 as the value of 𝑅 𝑜 could be stated through two facts: (i) Resistivity logs are mostly affected by reservoir fluids; therefore, they are measures of reservoir fluids. (ii) The resistivity of water is the lowest amongst other reservoir fluids.

A. Fuzzy Clustering

Finding available structures (clusters) within datasets is the main goal of clustering algorithms. The problem of finding several cluster centers that are appropriately representative of relevant classes of X, where X is a finite set of data [START_REF] Klir | Fuzzy Sets and Fuzzy Logic, Theory and Applications[END_REF]. The word "relevant classes" in the above sentence shows the impreciseness in defining clusters. In other words, the sentence mentioned is a general definition, and should be specified by each clustering algorithm. Hartigan and Wong, (1979) have defined k-means clustering as: "The aim of the k-means algorithm is to divide M points in N dimensions into K clusters so that the within-cluster sum of squares is minimized" [START_REF] Hartigan | Algorithm AS 136: A k-means clustering algorithm[END_REF].

In this study, three fuzzy clustering algorithms are applied to datasets to compare differences between the outputs of these algorithms. Also to find the most suitable one for the purpose of studying water saturation properties among Fuzzy C-Means (FCM), Gustafson-Kessel (GK) and Gath-Geva (GG) algorithms. The application of two other cluster analysis methods, Affinity Propagation (AP) and Spectral clustering, is also checked. The results of neither of the latter algorithms are included since their outputs were not well-interpretable.

npv=1 3 3 4 1 0 1 npv=2 3 0 1 1 0 0 npv=3 0 1 1 1 0 1 Petrophysical Well Logs Caliper (CALI)       Gamma Ray (GR)       Corrected Gamma Ray (CGR)       Sonic Log (ΔT)       Neutron Porosity (NPHI)       Bulk Density (RHOB)       Density Correction (DRHO)       Deep Laterolog Resistivity (LLD)       Shallow Laterolog Resistivity (LLS)       Microspherically Focused Log (MSFL)       Photoelectric Effect Log (PEF)       Core Tests Porosity       Permeability       Water Saturation Not Reliable Not Reliable Not Reliable Not Reliable Not Reliable  Fuzzy C-Means (FCM)
clustering is an extension of hard C-means (or k-means) clustering. The main difference between FCM and hard C-means is in a fuzzy weight, powered by a weighting exponent, which is multiplied by Euclidian distance values in cost function. FCM is not an optimum clustering algorithm for oriented or multi-size clusters. The FCM algorithm is introduced in four main steps in the literature [START_REF] Klir | Fuzzy Sets and Fuzzy Logic, Theory and Applications[END_REF][START_REF] Bezdek | FCM: The fuzzy cmeans clustering algorithm[END_REF]:

(1) Define fixed number of clusters, C, fuzziness, m, and norm type (i.e. distance or similarity function). Also randomly create similarity matrix, U, for all cluster centers.

(2) Compute the average value of similarities of all C clusters using the following equation:

𝑑(𝑥 𝑖 , 𝑤 𝑗 ) = ∑ 𝑢 𝑘𝑗 𝑚 ‖𝑥 𝑖 ,𝑤 𝑗 ‖ 2 𝑛 𝑘=1 ∑ 𝑢 𝑘𝑗 𝑚 𝑛 𝑘=1 ; 𝑗 ∈ [1, 𝑐]; 𝑖 ∈ [1, 𝑛] (4)
(3) Compute an updated membership matrix by:

𝑢 𝑘𝑗 = 1 ∑ [ 𝑑(𝑥 𝑖 ,𝑤 𝑙 ) 𝑑(𝑥 𝑖 ,𝑤 𝑗 ) ] 1 𝑚-1 𝑐 𝑗=1 ; 𝑙 ∈ [1, 𝑐]; 𝑖 ∈ [1, 𝑛] (5)
(4) Comparing Uk with Uk+1. If ‖𝑈 𝑘 , 𝑈 𝑘+1 ‖ ≤ 𝜀 , stop; otherwise set 𝑈 𝑘 = 𝑈 𝑘+1 , and go back to stage [START_REF] Worthington | The Role of Cut-offs in Integrated Reservoir Studies[END_REF].

The core difference between GK and FCM is in the use of the Mahalanobis distance, instead of the Euclidian one in cost function (norm in relation 4), which results in the recognition of oriented clusters [START_REF] Babuka | Improved covariance estimation for Gustafson-Kessel clustering[END_REF][START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF]. Therefore, GK tends to recognize oriented (linear at extreme mode) clusters. As in FCM, the identified clusters of GK are of the same size in terms of number of data, i.e. the number of data in each cluster is about the same.

The limitation of multi-size clusters is solved in the GG algorithm. The definition of cluster means and covariance matrix are derived from FCM and GK algorithms, respectively [START_REF] Gath | Unsupervised optimal fuzzy clustering[END_REF][START_REF] Abonyi | Modified Gath-Geva fuzzy clustering for identification of Takagi-Sugeno fuzzy models[END_REF]. The novelty of GG that results in solving the limitation of multi-size clusters is in defining distance based on decision surface. For example, a small blue cluster, a medium-sized green cluster and a large black cluster are detected by the GG algorithm in a single run (Figure 5c,d). The algorithm is much the same as FCM; and for detailed study, please refer to [START_REF] Abonyi | Modified Gath-Geva fuzzy clustering for identification of Takagi-Sugeno fuzzy models[END_REF].

INPUT SELECTION FOR CLUSTERING

Clustering algorithms do not consider any class label for data and only search the optimum clusters within input datasets. Therefore, selected features/datasets are very important in a clustering project. A logical interrelation between petrophysical parameters, inspired from [START_REF] Archie | Introduction to petrophysics of reservoir rocks[END_REF] is depicted in Figure 2.

Five features are related directly to saturation property: capillary pressure, permeability, fluid relief, fluid-solid interaction and electrical response. Four of the features are mostly categorized as dynamic properties, and usually we do not have any accurate estimation of these parameters, especially in exploratory stages. Only resistivity has a reliable and continuous recording throughout the reservoir. In the datasets of the current work, we have three resistivity log recordings: Deep Laterolog Resistivity (LLD), Shallow Laterolog Resistivity (LLS) and Microspherically Focused Log (MSFL). By trial and error, we found that removing MSFL provides better clustering results; therefore, two resistivity logs (LLS and LLD) are selected as input features of clustering algorithms. Cross plots of LLD to LLS show upper and lower bounds of data (Figure 3a), and the area of concentrated data (Figure 3b). fonts (modified after [START_REF] Archie | Introduction to petrophysics of reservoir rocks[END_REF]).

LLD TO LLS, A MEASURE OF PERMEABILITY

In Figure 3b, we have interpreted two arrows: (i) the direction of increasing water saturation, which is downward, and is taken for granted by researchers; (ii) the direction of increasing permeability, which is the main questionable issue of this figure. The interpretation (ii) is inspired from a previous publication, for which the ratio of LLD to LLS is used to relate invasion to permeability [START_REF] Sami | Permeability Prediction from Wireline Well Logs Using Fuzzy Logic and Discriminant Analysis[END_REF].

It is not far from our intuition of the reservoir conditions that the more permeability around the well, the more uniform the fluid distribution outward from the wellaxis. Here, LLS is representative of near-well conditions (including infiltration), compared to LLD which represents rather a larger zone around the wellbore. Therefore, in an extreme situation that recorded value of LLD is equal to LLS (i.e. ratio of LLD to LLS is one), the high permeability of the reservoir resulted in a mixture of reservoir fluids and drilling mud effectively around the well. In the other extreme case (i.e. LLD>>LLS), the mud cake and infiltration zone will be easily constructed due to the insufficiency of mixing drilling mud particles with reservoir fluids.

Besides all analytical descriptions, we have checked the correctness of this hypothesis (i.e. correlation between permeability and the ratio of LLD to LLS) in wells number one to five. This is fundamental to find the applicability of our analytics in real conditions and datasets. So, we plotted a logarithmic scale of permeability vs. the ratio of LLD to LLS in the five above-mentioned wells to investigate correlations (Figure 4). It is noteworthy that permeability is a lognormal property in most reservoir conditions.

In the visual assessment of Figure 4, the correlation between permeability and LLD:LLS, on a logarithmic scale could be easily found in wells two, three and five. There is no correlation in wells number one and four; however it does not mean that LLD:LLS is not a good measure of permeability in these two wells. This is because core permeability only shows permeability of intact rock on a small scale, without considering fractures and joints, whereas well-logs are volumetric responses and could be better representatives of real heterogeneous reservoir conditions. It is worth mentioning that from a structural viewpoint, well number four is drilled on a deviation of the axis of the anticline from its normal trend, i.e. the trend of the anticline changes from normal, which is NS, to another trend, which is NW-SE, at the emplacement of well number four. Therefore, lots of fractures and joint sets are probably to be found within well number four, thus increasing the permeability of the reservoir, while the permeability of the intact rock (cores) is unchanged. Intact rock properties are not controlled by structural geology. These properties are controlled by sedimentary processes and lithological properties. All in all, there is a correlation, however not that strong between LLD:LLS and permeability on the reservoir scale.

RESULTS: OUTPUT OF CLUSTERING METHODS

Clustering algorithms, FCM, GK and GG were run on the datasets (Figure 3). The number of clusters is set to be 3, 5, 10 and 20 to find the optimum number of clusters. The optimum number of clusters is found to be 5 for all mentioned algorithms. The outputs of all the fuzzy clustering methods are given in Figure 5.

Clusters of FCM are relatively of the same size, with no elongation. Clustering is done along the direction of "higher water saturation" (Figure 3b). Therefore, the blue cluster represents the lowest water saturation, (Figure 5a) while the red cluster shows the highest water content data. But as this algorithm is not able to detect lineation, it is unable to differentiate along the arrow of "Higher Perm and Larger Pore Throats" on Figure 3b. The green cluster shows two different trends: One trend belongs to high permeability values, the other represents the less permeable or shale contaminated part. This phenomenon is detected by GK clustering as this algorithm is able to distinguish elongated clusters. The blue cluster represents low water saturation and high permeability (Figure 5b), which is ideal to be known as net pay zone, i.e. a hydrocarbon production zone. The red data shows below oil water contact but still permeable, i.e. a water-producing zone. The area of green data is the least permeable part of the Sarvak Formation, either oil-bearing or water-bearing. The black and magenta parts have less permeability than red and blue clusters; and magenta contains less water, compared to black data. GK provided the best clustering due to petrophysical interpretations.

While GK clusters are relatively similar in size (number of data in each cluster), the GG clusters differ considerably in size (Figure 5c,d) with a net pay zone represented by the magenta colour (Figure 5c). Red, black, green and blue clusters have lower reservoir quality, respectively. Contrary to FCM and GK, clusters of GG do not show a specific pattern; therefore they cannot be as easily interpreted as FCM and GK clusters.

DISCUSSION

In the previous section, the outputs of three clustering methods on a two-dimensional dataset were tested. FCM was the most compatible tool for saturation prediction, while GK is the best for permeability study and GG could be used as a rough measure of productivity. The reason why the functionality of these clustering methods differ from each other could be seen at the bottom of Table 2.

FCM provides simple clusters of the same size, while GG provides elongated unequal clusters, which is difficult to interpret due to the complexity of the structure of clusters. On the other hand, GK provides elongated clusters (same size), but the structure of the clusters is not complex for interpretation. It is not difficult to imagine that when we have a change in the specifications of clusters, their physical meaning (characterization target) would be changed too. So, three issues have to be kept in mind for successful clustering: (i) our goal of study or characterization target; (ii) the relation or correlation of input dataset with target, which was discussed above; and (iii) the specifications of produced clusters. Among the three issues mentioned, first we define the purpose of our study; i.e. "What is to be measured or predicted?" According to this goal, the dataset is chosen; e.g. here we used resistivity logs for the sake of their correlation with saturation and reservoir fluids. It is noteworthy that sometimes the dataset is not optional, and we have to use some fixed parameters. But we are always able to choose among various clustering algorithms. To choose the most appropriate one, we can either go in deep through the mathematical basis of the algorithm, analytically; or we can simply check the specifications of their clusters as is provided in Table 2.

CONCLUSIONS

This paper shows the strength of clustering algorithms in interpreting resistivity well-logs to extract various inferences regarding water saturation, and other reservoir parameters just by means of two well-known resistivity well-logs. The number of publications on saturation prediction is much less than other petrophysical properties such as porosity and permeability, mostly due to a lack of real saturation data for the validation procedure. However, its importance is not less than other reservoir properties. Furthermore, we proposed clustering algorithms for water saturation prediction, because clustering methods are unsupervised methods that do not require real labeled data for training the machine.

Due to petrophysical interpretations, the FCM method provided separation of the space only by considering saturation property. However, GK clustering divided the datasets according to two properties: saturation and permeability (or pore throat size). GG could divide the space according to overall reservoir quality. The output of the Spectral algorithm was much the same as GG but more noise-contaminated. The affinity propagation clustering method did not provide well-interpretable outputs though we are not disappointed with this algorithm, and believe that further investigations have to be carried out on the application of clustering algorithms, on petrophysical datasets.

It is noteworthy that by means of unsupervised methods, we do not train the machine, considering cores, i.e. the minute scale. Therefore, the results would be free from the falsification resulting from the use of idealistic small scale measurements of cores, when interpreting highly heterogeneous reservoir conditions based on core analysis. The validation procedure is done qualitatively by checking the behaviour of clusters against the behaviour of core permeability. The ratio of LLD:LLS has established a link between the reservoir permeability and core permeability. Quantitative comparison is not carried out here because: (i) there are two different languages in expressing permeability values: core permeability, in milliDarcy, is a continuous value, while inferred well-log permeability is in discrete labels (i.e. clusters); (ii) There is no guarantee for compatibility of core-scale permeability with reservoir-scale permeability, inferred from well-logs. Hence, we lack a reliable base point for quantitative evaluation of water saturation and permeability in reservoir-scales; however, qualitative reasoning is possible and necessary as is done here.

Figure 1 .

 1 Figure 1. Comparing well-test results and water saturation (Archie) value in five wells of Sarvak Formation. Horizontal axis is a simple data numbering, which makes saturation values in ascending order. Well-test class value means water production (1), oil production with a rate of less than (2) or of higher than (3) 1500 barrel per day.

Figure 2 .

 2 Figure 2. Interrelations between petrophysical features. Water saturation and directly related features are highlighted by color and boldfonts (modified after[START_REF] Archie | Introduction to petrophysics of reservoir rocks[END_REF]).

Figure 3 .

 3 Figure 3. Cross plots of input features. Trends show interpretation of petrophysical interrelations. (a) Whole the data, and (b) focused on dense area. Arrows show interpretation of permeability, pore throat size and saturation on cross plot.

Figure 4 .

 4 Figure 4. Checking correlation of permeability (mD) and LLD (Ω.m) to LLS (Ω.m) in five wells, logarithmic scale.

Figure 5 .

 5 Figure 5. Outputs of fuzzy clustering algorithms. (a) Fuzzy C-Means clustering with 5 clusters. (b) Gustafson-Kessel with 5 clusters. (c) Gath-Geva with 5 clusters. (d) Gath-Geva with 5 clusters. Focused to show blue cluster.

TABLE I .

 I DATASET AVAILABLE FROM TWO WELLS OF THE OIL FIELD.

	"NPV" STANDS FOR "NET PAY VALUE". NPV=1 MEANS THAT THE
	INTERVAL OF TESTING IS NOT OIL PRODUCING, NPV=2 AND 3 MEAN
	THAT THE TESTED INTERVAL PRODUCES OIL LESS AND MORE THAN
	1500 BARREL PER DAY ( IN SI UNITE SYSTEM) RESPECTIVELY.
	AVAILABLE AND UNAVAILABLE INFORMATION ARE SHOWN BY
	SYMBOLS () AND (), RESPECTIVELY. FURTHER INVESTIGATION IS
	LIMITED TO AVAILABLE LOGS ONLY.
		Well 1	Well 2	Well 3	Well 4	Well 5	Well 6
	No. of Well	Test Intervals

TABLE II .

 II SUMMARY OF APPLICATION OF CLUSTERING METHODS IN RESERVOIR CHARACTERIZATION. SYMBOLS () AND () MEAN COMPATIBILITY AND INCOMPATIBILITY, RESPECTIVELY; E.G. FCM IS COMPATIBLE FOR PREDICTING WATER SATURATION BUT IT IS NOT COMPATIBLE FOR PERMEABILITY INFERENCE. "SIMPLICITY OF CLUSTERS" SHOWS WHETHER THE CLUSTERING ALGORITHM COULD BE INTERPRETED EASILY OR NOT; E.G. GG COULD NOT BE INTERPRETED SIMPLY; WHILE GK IS EASILY INTERPRETABLE. THE DOUBLE SYMBOL () MEANS MORE COMPATIBILITY, COMPARED TO THE SINGLE SYMBOL ().

				Fuzzy C-	Gustafson-	Gath-
				Means	Kessel	Geva
	characteriza	tion target	water saturation permeability productivity	  	  	  
	properties of	clusters	equity of cluster size elongation interpretability	  	  	  

ACKNOWLEDGMENT

This work has been supported by the Center for International Scientific Studies & Collaboration (CISSC) and French Embassy in Iran through PHC Gundishapur Program. It is our pleasure to acknowledge the Exploration Directorate of National Iranian Oil Company (NIOC) for providing data, and permitting the publication of scientific achievements.