Aude Bolopion
email: aude.bolopion@femto-st.fr

Zhenjiang Ni

Joël Agnus
email: joel.agnus@ens2m.fr

Ryad Benosman

Stéphane Régnier

Stable Haptic Feedback based on a Dynamic Vision Sensor for Microrobotics

This work presents a stable vision based haptic feedback for micromanipulation using both an asynchronous Address Event Representation (AER) silicon retina and a conventional frame-based camera. At this scale, most of the grippers used to manipulate objects lack of force sensing. High frequency vision detection thus provides a sound solution to get information about the position of the object and the tool to provide virtual haptic guides. Artificial retinas present high update rates, which enables to address one of the major challenge of haptic feedback teleoperation systems, namely stability. However static objects are not detected. The haptic feedback is thus based on an asynchronous silicon retina to provide a high update rate of moving objects and a framebased camera to retrieve the position of the target object. This approach is validated by pick-and-place of microspheres (diameter: around 50 micrometers) using a piezoelectric microgripper. The displacement of the tool, as well as the opening and closing of the gripper are controlled by the haptic device. Haptic virtual guides are transmitted to users to assist them in the different steps of the pick-and-place task: a virtual stiffness ensures the correct alignment of the tool with respect to the object, a repulsive haptic force enables users to monitor the gripping step, and operators are assisted while picking and placing the object.

I. INTRODUCTION

Micromanipulation has received an increasing interest in the last few years, since the ability to handle objects with a size below 100 micrometers will open a wide range of applications, from electronics to biological fields. Several manipulation techniques have been developed, based on rolling, pushing-pulling or picking-and-placing techniques [START_REF] Haliyo | Mumad, the adhesion based dynamic micro-manipulation[END_REF], [START_REF] Haliyo | Controlled rolling of microobjects for autonomous manipulation[END_REF], [START_REF] Xie | Three-dimensional automated micromanipulation using a nanotip gripper with multi-feedback[END_REF]. However many issues must still be overcome to get intuitive manipulation platforms. In particular, only skilled operators can perform manual micromanipulation due to the scale difference, the fragility of both the tools and the objects, the complexity of the force fields and the high sensitivity of the systems to environmental conditions. One solution is to perform automated tasks. Good results have been obtained (i.e. high throughputs and accurate positioning) on objects with a size of hundreds of micrometers [START_REF] Tamadazte | A direct visual servoing scheme for automatic nanopositioning[END_REF]. This manipulation mode fits the needs of repeated tasks on large number of objects. However, operator's knowledge is not exploited, users cannot interact while the automated task is being performed and they do not receive enough feedback to learn from the manipulation. The use of assistance for manipulation appears as a promising alternative, which benefits from the interactivity of manual manipulation, while facilitating it. Haptic feedback is one of the solution to intuitively transmit information to operators, as it provides indications that users are used to receive when they perform manipulations at macroscales [START_REF] Ferreira | Virtual reality and haptics for nanorobotics[END_REF], [START_REF] Bukusoglu | Haptic manipulation of microspheres using optical tweezers under the guidance of artificial force fields[END_REF].

Several teleoperated micro-and nanoscale tasks are reported in the literature, from feeling the substrate topology [START_REF] Hollis | Toward a tele-nanorobotic manipulation system with atomic scale force feedback and motion resolution[END_REF] to inducing a modification of the sample through indentations (e.g. direct patterning on a substrate [START_REF] Li | Videolized" atomic force microscopy for interactive nanomanipulation and nanoassembly[END_REF]). Tasks involving pushing/pulling or cutting objects are also presented [START_REF] Guthold | Investigation and modification of molecular structures with the nanoManipulator[END_REF], [START_REF] Sitti | Teleoperated touch feedback from the surfaces at the nanoscale: modeling and experiments[END_REF], [START_REF] Onal | Teleoperated 3-D force feedback from the nanoscale with an atomic force microscope[END_REF]. In addition to transmitting micro-and nanoscale interactions, haptic feedback is used as an enhancement for user assistance, with virtual guides for pushing and pick-andplace by adhesion tasks [START_REF] Ammi | Robotic assisted micromanipulation system using virtual fixtures and metaphors[END_REF]. In this case, haptic feedback is used to keep the user's motion on a specified path. However these examples are mostly limited to proof of concepts, and none of them can be extended to industrial applications.

To get an efficient and versatile teleoperation system for microscale applications a major issue must be faced: the lack of position and force feedback [START_REF] Chaillet | Microrobotics for Micromanipulation[END_REF]. Sensors have been developed [START_REF] Fahlbusch | Force sensing in microrobotic systemsan overview[END_REF], [START_REF] Beyeler | A six-axis MEMS forcetorque sensor with micro-Newton and nano-Newtonmeter resolution[END_REF], but their integration into the dedicated tools induces an important increase in the complexity and the cost of their fabrication. In particular, even if some microgrippers offer sensing capabilities at the expense of complex designs [START_REF] Duc | Electrothermal microgripper with large jaw displacement and integrated force sensors[END_REF], [START_REF] Kim | Micronewton force-controlled manipulation of biomaterials using a monolithic MEMS microgripper with two-axis force feedback[END_REF] most of them still lack force measurements [START_REF] Lopez-Walle | Principle of a submerged freeze gripper for microassembly[END_REF], [START_REF] Andersen | Multimodal electrothermal silicon microgrippers for nanotube manipulation[END_REF].

To avoid the complexity of sensors' integration, vision is a promising solution [START_REF] Chawda | Vision-Based force sensing for nanomanipulation[END_REF], [START_REF] Reddy | Miniature compliant grippers with Vision-Based force sensing[END_REF]. However, to ensure a stable haptic coupling, a frequency of 1 kHz is needed for the sensor output [START_REF] Colgate | Passivity of a class of sampled-data systems: Application to haptic interfaces[END_REF]. This is even more critical in case of micromanipulation since the dynamic of the objects is important due to their low inertia. Conventional frame-based cameras' lack of dynamic information and their redundancies set an important limit to potential micromanipulations. Eventbased computer vision based on Address Event Representation (AER) provides a sound solution to high-speed vision problems [START_REF] Delbruck | Frame-free dynamic digital vision[END_REF]. This newly developed discipline is motivated by mimicking biological visual systems [START_REF] Delbruck | Activity-driven, event-based vision sensors[END_REF]. The Dynamic Vision Sensor (DVS) silicon retina reacts to changes of contrast that are converted in a stream of asynchronous timestamped events [START_REF] Lichtsteiner | A 128*128 120 db 15us latency asynchronous temporal contrast vision sensor[END_REF]. The reduction of redundant information makes this technique promising for high-speed tracking.

This paper presents a haptic feedback teleoperation system for microscale applications based on visual sensing obtained from an asynchronous Address Event Representation (AER) silicon retina and a conventional frame-based camera. The Dynamic Vision Sensor enables to record high frequency phenomena, and to provide feedback at a high sampling frequency, but it gives no indication about immobile objects. The position of static objects are provided by a conventional frame-based camera. To validate this approach, a teleoperated pick-and-place task of spheres (diameter around 50 µm) using a piezoelectric gripper is presented. Haptic feedback directly estimated from the output of both the event-based retina and the conventional frame-based camera is provided to assist users during the manipulation. The first pick-and-place manipulation with 3D haptic feedback using a microgripper is successfully achieved. This work is based on several previous studies on the stability of haptic coupling schemes for applications at microscales [START_REF] Bolopion | Analysis of stability and transparency for nanoscale force feedback in bilateral coupling[END_REF], on the definition of virtual guides for pick-and-place operation of microspheres based on two atomic force microscopy cantilevers [START_REF] Bolopion | Haptic teleoperation for 3D microassembly of spherical objects[END_REF], on the use of visual feedback coming from a scanning electron microscope for teleoperation [START_REF] Bolopion | Remote microscale teleoperation through virtual reality and haptic feedback[END_REF] and on the use of dynamic vision sensors [START_REF] Ni | Asynchronous event based high speed vision for micro-particles tracking[END_REF]. However, it is the first time, to our knowledge, that dynamic vision sensors are used to provide stable haptic feedback to assist users while performing a pick-and-place operation. This approach will enable complex teleoperated assembly tasks.

This paper is organized as follows. The experimental setup is presented in Section II. Section III details the teleoperation system. The proposed approach is validated by experimental results presented in Section IV. Finally, Section V concludes the paper.

II. SETUP

A. Micromanipulation tool

The micromanipulation end-effector developed in FEMTO-ST Institute [START_REF] Lit | A four-degree-offreedom microprehensile microrobot on chip[END_REF] is a microgripper with tightened fingers that proves to be efficient and intuitive for microhandling problems (Fig. 1). This microgripper, named MMOC (Microprehensile Microrobot On Chip), has 4 degrees of freedom and is able to grasp, hold and release submillimetric-sized objects up to several tens of micrometers. The Fig. 1(d) and 1(e) show respectively the manipulation of a 20 µm silicate grain and 100 µm glass sphere.

Each finger is actuated by a piezoelectric cantilever with local electrodes, called duo-bimorph (Fig. 2). As mentioned on the cross section view, one finger has four electrodes referred to a central ground and two voltages (V 1 = V z -V y and V 2 = V z +V y ) are necessary to impose the displacements y and z based on the deflections of the piezoelectric beam. This configuration offers number of capabilities: open-andclose motion as well as up-and-down motion. Indeed, each finger is able to move independently from the others in two orthogonal directions (right picture of Fig. 2). The upand down motion can be useful for fine motion, for release strategies of objects by crossing the fingers or for insertion of microparts. It is also particularly convenient to align the finger-tips in case of misalignment after microfabrication and assembly process of the microgripper. Specific silicon end-effectors able to grasp objects from 10 to 100 micrometers have been made (Fig. 1(a) and 1(b)). They are designed to minimize the sticking effects between the end-effector and the objects, to facilitate the release step. They are microfabricated on SOI wafer with a thickness of 10 µm, mainly using DRIE process [START_REF] Agnus | Silicon endeffectors for microgripping tasks[END_REF]. The gap and the alignment of the tips are assured by breakable-parts. These parts are removed after the pair of fingers have been glued on the spatulas of the actuators (Fig. 1(a) and 1(b)).

Such a microgripper presents a typical stroke of open/close motion and up/down motion of respectively 320 and 200 micrometers at the end of the finger tips for ±100V and a blocking force of respectively 80 and 30 mN at the end of the actuator for 100 V. The gripper is controlled by commands sent from a personal computer to a high voltage interface (four channels of ±150V) via RS232 link.

B. Vision system

As shown in Fig. 3, the observed scene is monitored by two optical sensors, that record the same view. The light beam is divided into two optical paths, and redirected to an asynchronous event-based silicon retina (DVS, prototype from ETHZ) and a conventional frame-based camera (GigE vision, Basler). The scene recorded by both sensors is magnified by a 10× objective (Olympus).

The DVS used in this work is a silicon retina with 128 × 128 pixels [START_REF] Lichtsteiner | A 128*128 120 db 15us latency asynchronous temporal contrast vision sensor[END_REF]. It mimics the biological retinas and outputs compressed digital data in the form of events, removing redundancy and increasing dynamic range compared with conventional imagers. A review of the existing similar sensors are surveyed in [START_REF] Delbruck | Activity-driven, event-based vision sensors[END_REF]. The DVS's output consists of asynchronous address-events that signal scene luminance changes. Each pixel is independent and detects changes in log intensity larger than a threshold since the last event is emitted. When the change exceeds a set threshold, a +1 or -1 event is generated the pixel depending on whether the log intensity is increased or decreased. The timing of events can be conveyed with a very accurate temporal resolution of 1 µs. The events are then transmitted to a computer using a standard USB connection. The advantage of this sensor over conventional clocked cameras is that only moving objects produce data, thus increases significantly the processing speed. The reactivity and the stability of the haptic system can therefore be enhanced.

C. Teleoperation system

The gripper is mounted on a 3-axis motorized micromanipulator (MP285, Sutter Instrument) 1 to allow an accurate positioning with respect to the substrate. The manipulator used relies on stepper motors with a step size of 0.040 µm. It is a cable-driven system with cross roller bearings, and it has a sub micrometer resolution and a travel range of 25 mm. The manipulator was originally controlled through a serial port; however, to increase the communication frequency, a joystick is emulated by programming the manipulator's parallel port using a PCI6259 National Instrument acquisition card.

Both the positioning of the gripper and its actuation (opening and closing operations) are controlled through an Omega haptic device 2 . This device has 3 degrees of freedom for both displacement and force feedback. Haptic feedback is computed based on the vision detection. Forces higher than 5 N are saturated to avoid any damage to the interface.

A single PC (Intel Xeon core, 2.93 GHz) operating under Windows 7 runs the threads corresponding to the gripper, the micromanipulator, the vision detection, and the haptic feedback.

III. HAPTIC FEEDBACK TELEOPERATED SYSTEM

A. Manipulation modes

The coupling between the haptic device and the microgripper is represented in Fig. 3. Using the haptic interface, the 1 Sutter Instrument, http://www.sutter.com/MP 285 2 Force Dimension, http://www.forcedimension.com user controls the displacement of the micromanipulator as well as the opening and the closing of the gripper. Scaling factors α d and α oc are introduced to convert the position of the haptic handle to the variables used to control the positioning and the actuation of the gripper. To ensure ease of manipulation, different modes have been defined:

• planar displacement: the operator controls the displacement of the gripper in a plane parallel to the substrate, the (x, y) plane, • vertical displacement: the operator controls the displacement of the gripper along the vertical direction. It enables to lift-off the object, and to place it down, • gripper control: the operator controls both the opening and closing of the gripper and its position along the y-axis to align the gripper with the middle line of the object. The operator selects the appropriate manipulation mode on a graphical user interface developed in C++.

B. Determination of the relative position of the object and the tool

The haptic feedback is based on the distances between the gripper and the object. To determine these distances, both the DVS silicon retina and the frame-based camera are used.

The DVS sensor is installed to provide a fast vision detection. To track the gripper, the most recent events within a time period ∆t (called "active events" hereafter) are continuously stored in memory. Following the Event based Iterative Closest Point (EICP) algorithm, a model composed of points of the contour of the microgripper is required. An active event is matched with a model point by computing the minimal Euclidean distance between the event's position and all the points of the model. After the matching step, the rigid body transformation composed of rotations and translations is estimated by minimizing a mean square cost function. Readers interested in the algorithm details can refer to [START_REF] Besl | A method for registration of 3-d shapes[END_REF]. In our manipulation, the active events' decay time ∆t is tuned around 10 -15 ms. The processing is event driven, and its The calibration between the classical image and the DVS accumulation map (inset). Six points (crosses) have been to calculate the homography transform H. δ xl and δ xr describe the distance between the gripper's fingers and the sphere on the left side and the right side, δ x is the distance between the two fingers and δ y is the distance between the center of the fingers and the center of the sphere. update frequency has a mean value of 4 kHz. The algorithm is implemented in Java under JAER open-source project [33].

The frame-based camera serves as a complement to DVS silicon retina solely for the static object detection. The focal planes of both the DVS (128 × 128 pixels) and the frame camera (659 × 494 pixels) are related by a homography transform as both observe the same 2D plane [START_REF] Hartley | Multiple View Geometry in computer vision[END_REF]. The homography is estimated off-line by extracting from both sensors' the coordinates of six corner points of the gripper fingers and linking them to the actual metric of the gripper's points in the scene (see Fig. 4) [START_REF] Ni | Asynchronous event-based visual shape tracking for stable haptic feedback in microrobotics[END_REF]. During the manipulation, the circle corresponding to the sphere to be manipulated is detected using a Hough transform through conventional camera output. The distance between the gripper fingers δ x is estimated from DVS output. By combining the processing of the two sensors, if an object is detected between the two fingers, the relative finger-object distance on the left and the right sides δ xl , δ xr and the distance between the center of the sphere and the gripper δ y (Fig. 4) is calculated instead. These various distances are used to compute the haptic feedback.

C. Haptic feedback

To increase the success rate of the pick-and-place operation two criteria should be met: the sphere should be grasped on its middle line, and the grasping force should be enough to lift the sphere but controlled to avoid any damages to the object. The haptic feedback must assist the user for these two operations. The operator receives a force feedback F h T = F oc 0 F y based on the distance between the gripper and the sphere δ x , δ xl , δ xr and δ y determined from vision. The haptic force F y is transmitted on the vertical axis of the haptic device to ensure a correspondence between the visual feedback (on the screen the gripper is moving in the vertical plane) and the haptic feedback.

To help the user align the gripper with respect to the middle line of the sphere, a haptic force corresponding to a spring of stiffness k between the position of the gripper and the sphere is provided:

F y = -kδ y (1)
where δ y is the distance between the center of the gripper and the center of sphere along the y-axis (see Fig. 4). A haptic feedback F oc is provided so that the user can monitor the grasping force. Contrary to what is commonly presented in the literature [START_REF] Cappelleri | Two-dimensional, visionbased µN force sensor for microrobotics[END_REF], we are not interested here in computing the exact efforts applied on the object but only in deriving information to assist the user while performing a given task. The calibration process, which enables to relate the tool deformations to the applied force, is thus unnecessary. While closing the gripper, the user has to counteract a haptic force F oc :

F oc =    F max e -δ f x 2 α if not in the contact zone F contact if in the contact zone (2)
where δ f x is the free space between the two gripper's fingers. If the sphere is situated between the fingers δ f x = δ xl + δ xr , which corresponds to the sum of the distances between each of the fingers and the sphere; otherwise, δ f x = δ x . F max is the maximum force that can be transmitted to the user when the gripper is close to the sphere but has not entered the contact zone yet. α is a constant chosen to tune the decrease of the haptic force as the distance between the two fingers increases. F contact is the force sent while the gripper is grasping the sphere. The step between F max and F contact must be high enough to indicate clearly the contact between the sphere and the gripper. The contact zone is reached if δ xl and δ xr are less than a given distance δ 1 = 3 µm (which corresponds to 6% of the sphere diameter). The gripper will then reach the non contact zone if δ xl and δ xr are greater than δ 2 = 7.5 µm. This hysteresis avoids undesirable frequent transitions between contact and non-contact modes because of noises and tracking errors. The hysteresis values δ 1 and δ 2 are chosen according to our experience to achieve a comfortable user sensation. The force step that is sent to the user when contact is detected is filtered to avoid brutal force changes. Even if the user does not receive the maximum force feedback at the instant of contact, he/she can distinctly feel the increase in the force, and infers that contact happened.

IV. EXPERIMENTAL RESULTS

Experiments are performed to validate the integration of event-based vision sensor for stable haptic feedback. The microspheres are glass beads of 30 ∼ 50 micrometers of diameter from Polysciences, Inc3 . A Gel-Pack substrate has been used to provide enough adhesion to prevent the sphere from sticking to the gripper. The experiment consists in the following tasks: positioning the gripper with respect to the sphere (in plane displacements), grasping it, taking it off, moving it, placing it down, and releasing it. The evolution of the haptic force F y and F oc that help users monitoring the operation is illustrated in Fig. 5 for both grasping (a) and releasing stages (b). In Fig. 5(a), the haptic force along the y-axis assists the user in aligning the gripper with respect to the sphere. Users control the position of the gripper along the y-axis. A haptic force F y that corresponds to a virtual stiffness between the center of the gripper fingers and the center of the sphere (Eq. ( 1)) is transmitted to the operator. At the beginning of the experiment, the gripper is misaligned, and the user feels an attractive force that pulls him/her to the correct position. 13 s, the gripper is correctly aligned so that haptic feedback drops to zero. For time inferior than 11.8 s, the user closes the gripper towards the sphere. As the free space between the gripper's fingers and the sphere decreases, the operator has to counteract an increasing haptic force F oc (Eq. ( 2)). At 11.8 s, the gripper enters the contact zone, and the user feels a sudden increase of the haptic force. The sphere is grasped.

The user then begins the pick and place operation. During the lifting and the placing operations, a constant repulsive haptic force field, set to 2 N, is provided to avoid any involuntary contact with the substrate. To perform the vertical displacements a position scaling factor α z d = 12.5 × 10 3 is used. When the sphere has been lifted above the substrate at a desired height, the user can move it freely in the (x, y) plane parallel to the substrate. During this procedure, the haptic feedback is turned off. To move the sphere in the plane parallel to the substrate, the parameters along x, y axes are set to be α x d = 4.0 × 10 3 , α y d = 2.5 × 10 3 . Different factors are used along the three axes of the micromanipulator to achieve easy positioning. They are set according to the user's comfort of manipulation. Factors that are too high will lead to very small displacements and a time consuming experiment, whereas low factors will not enable a precise positioning. When the sphere is at the correct location, the user can place it down.

The last step is to release the sphere (see Fig. 5(b)). The plot starts at t = 44 s approximately since it took around 30 s to the user to lift the object, move it, and place it down. This duration depends on the distance between the starting and final positions, as well as on the user expertise with haptic devices. At the beginning of the operation, as the gripper contacts the sphere, the operator can feel a constant force (5 N) that assist him/her to open the gripper. The user releases the sphere at time 45.3 s. It can be noted that the curve does not reach zero as the force is still helping the user to open the gripper and avoids unexpected closing. Since the gripper and the sphere are aligned during releasing, F y is around zero.

During all the grasping-releasing and pick-and-place operation, the user receives haptic feedback based on vision detection that helps him/her to achieve comfortable and easy sphere manipulation. The dynamic vision sensor provides a high-frequency feedback that enables to provide a stable haptic system. Users successfully performed a 3D teleoperated manipulation on micron-sized objects. This system will surely benefit teleoperated or automated microassembly and opens new perspectives for complex micromanipulations. In the future a vertical haptic feedback should be provided to assist the user while positioning the gripper along the vertical axis during the grasping and the placing steps. To determine the distance between the gripper and the substrate the focusing and defocusing information might be used.

V. CONCLUSION

To overcome the lack of sensing capabilities at microscale, a vision-based system is presented. To ensure a stable haptic feedback, the frequency of the vision detection must be higher than 1 kHz at this scale as low inertia induces high dynamic phenomena. This is ensured by the output of a DVS sensor that enables to track the tool at more than 4 kHz. This feedback is combined with the output of a classical framebased camera, used to derive information about static parts of the scene, and in particular the position of the object that must be manipulated. This approach is validated on a pickand-place experiment of glass spheres with a diameter of around 50 µm using a piezoelectric gripper. Haptic feedback provided to assist users while aligning the gripper with respect to the tool, grasping it, picking and placing it and releasing it. This work will enable to provide stable haptic feedback to perform complex assembly tasks even with sensor-deprived systems.

Future works include extending this approach to other applications, involving different objects or tools. Automated tasks could also benefit from this high frequency vision detection.

VI. ACKNOWLEDGMENT

This work was supported by the French National Agency of Research, through the NANOROL project (contract ANR-07-ROBO-0003) and the NANOROBUST project (contract ANR 2011 NANO 006).

Fig. 1 .

 1 Fig. 1. FEMTO-STs microgripper. (a) and (b): end effectors mounted on breakable parts for easier mounting, (c): MMOC gripper, (d) and (e): handling of micronsized objects.

Fig. 2 .

 2 Fig. 2. Piezoelectric actuators on which the fingers (Fig. 1 (a)-(b)) are glued. Each finger is actuated by the piezoelectric actuator with four electrodes to enable four degrees of freedom of mobilities.

Fig. 3 .

 3 Fig.3. Haptic coupling scheme. The user sets the position of the gripper and the actuation of the gripper's fingers by controlling the position of the haptic interface P h . Scaling factors α d and α oc convert the position of the haptic handle to commands to control the gripper. The operator receives haptic feedback through the device. The haptic force is based on the distance between the gripper and the sphere, determined from vision algorithms.

  Fig.[START_REF] Tamadazte | A direct visual servoing scheme for automatic nanopositioning[END_REF].The calibration between the classical image and the DVS accumulation map (inset). Six points (crosses) have been to calculate the homography transform H. δ xl and δ xr describe the distance between the gripper's fingers and the sphere on the left side and the right side, δ x is the distance between the two fingers and δ y is the distance between the center of the fingers and the center of the sphere.

Fig. 5 .

 5 Fig.5. The haptic force during the experiment. (a) The haptic force F y assists the user in aligning the gripper with respect to the middle line of the sphere during grasping operation. Equation (1) is used with the following parameters: k = 50000 N • m -1 . The displacement scaling along the y axis is set to α y d = 2.5 × 10 3 . Haptic force F oc used to control the grasping force is computed using Equation (2) with the following coefficients: F contact = 5 N (set to the maximum admissible force of the haptic interface), F max = 2 N, α = 1.44 • 10 -8 . The opening/closing scaling factor is set to α oc = 1.8 × 10 -3 m.V -1 . (b) The same types of forces F y and F oc are illustrated for releasing the microsphere. Parameters are set the same as above.

Institut des Systèmes Intelligents et de Robotique, Université Pierre et Marie Curie, CNRS UMR 7222, 4 Place Jussieu, 75005 Paris, France.{ni, regnier}@isir.upmc.fr

Vision Institute, Université Pierre et Marie Curie, UMR S968 Inserm, UPMC, CNRS UMR 7210, CHNO des Quinze-Vingts, 17 rue moreau, France. ryad.benosman@upmc.fr

http://www.polysciences.com/


N Verrier
email: verrier@coria.fr

S Coëtmellec
email: coetmellec@coria.fr

M Brunel

D Lebrun

Digital in-line holography in thick optical systems: application to visualization in pipes

In this paper we apply digital in-line holography to image opaque objects through a thick plano-concave pipe. Opaque fibers and opaque particles are considered. Analytical expression of the intensity distribution in the CCD sensor plane is derived using generalized Fresnel transform. The proposed model has the ability to deal with various pipe shape and thickness and compensates for the lack of versatility of classical DIH models. Holograms obtained with a 12 mm thick plano-concave pipe are then reconstructed using fractional Fourier transform (FRFT). This method allows us to get rid of astigmatism. Numerical and experimental results are presented.

Introduction

Digital in-line holography (DIH) is a recognized optical technique for flow measurements in transparent liquid media. This method is widely used in various domains such as fluid mechanics [START_REF] Pan | Digital holography of particle fields: reconstruction by use of complex amplitude[END_REF][START_REF] Malek | Digital in-line holography for threedimensional-two-components particle tracking velocimetry[END_REF] where the flow is seeded with small particles, or in the biological microscopic imaging field [START_REF] Malkiel | The three-dimensional flow field generated by a feeding calanoid copepod measured using digital holography[END_REF][START_REF] Garcia-Sucerquia | Digital in-line holographic microscopy[END_REF].

Nevertheless, a transparent pipe can be viewed as an optical system that can introduce aberrations such as astigmatism [START_REF] Vikram | Particle field holography[END_REF]. These aberrations make difficult to image objects in the pipe. Compensation for these unwanted effects has been widely investigated.

De Nicola et al. [START_REF] De Nicola | Correct-image reconstruction in the presence of severe anamorphism by means of digital holography[END_REF] proposed a numerical method to compensate for anamorphism. By modifying the chirp function and the spatial frequency term of the Fresnel integral the authors managed to compensate for a severe anamorphism brought by a reflexion diffraction grating. Numerical methods for astigmatism compensation have been proposed in Refs. [START_REF] De Nicola | Correct-image reconstruction in the presence of severe anamorphism by means of digital holography[END_REF][START_REF] Grilli | Whole optical wavefields reconstruction by Digital Holography[END_REF][START_REF] De Nicola | Wave front reconstruction of Fresnel off-axis holograms with compensation of aberrations by means of phase-shifting digital holography[END_REF]. Here the authors used a modified chirp function where two different propagation distances are considered. Astigmatism can also be compensated for by optical means. For instance, the authors of [START_REF] Crane | Far-field holography of ampule contaminants[END_REF] used index matching to minimize astigmatism: the thin-lens like effect of an ampule filled with contaminants was investigated to determine index matching parameters for astigmatism compensation.

However, the pipe used to carry a flow must be considered as a thick and cylindrical optical system. Recently, an analytical solution of the scalar diffraction produced by an opaque disk illuminated by an elliptical, astigmatic and Gaussian beam under Fresnel approximation has been proposed [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]. In this publication, the astigmatism was introduced and controlled by using a thin plano-convex cylindrical lens. Using fractional Fourier transformation (FRFT), authors managed to retrieve a correct image of an opaque disk from astigmatic hologram. It was shown that FRFT is therefore well suited for these studies where classical methods such as Fresnel transform [START_REF] Onural | Digital decoding of in-line holograms[END_REF] or Wavelet transform [START_REF]Diffraction from a wavelet point of view[END_REF] fail.

The astigmatism introduced by thin lenses is thus well controlled. However, these DIH models can not be applied to hugger pipes. As a matter of fact, the optical thickness of these has to be taken into account.

In this paper, the aim is to apply DIH to thick pipe systems. In the first part of the paper, the intensity distribution, in the CCD plane, of the diffracted field produced by an opaque object located in a pipe is calculated. A method of Gaussian functions superposition is used to describe the object function. In the second part, definition of the FRFT is recalled and we demonstrate its ability to reconstruct holograms recorded in such systems. Finally simulations and experiments are performed to illustrate our results.

In-line Holography through pipes

Holography aims to record, on a CCD camera without objective lens, the intensity distribution of the diffraction pattern of an object illuminated by a monochromatic continuous wave [START_REF] Nicolas | Digital in-line holography with a subpicosecond laser beam[END_REF]. The numerical and experimental set-up is represented in Fig. [START_REF] Pan | Digital holography of particle fields: reconstruction by use of complex amplitude[END_REF]. The incident Gaussian beam propagates in free space over a distance z p and illuminates the pipe. Here, the pipe is modeled as two thick lenses. Their thickness is denoted by e. The opaque object is located between these two lenses: at a distance δ from the first thick lens, and at a distance z i from the second lens. The CCD sensor is located at the distance z next the pipe and records the intensity distribution of the diffraction pattern.

2.A. Intensity distribution in the CCD sensor plane

In this part, we consider the propagation of a Gaussian beam through our optical system. In the beam waist plane, located at a distance -z p from the pipe, the Gaussian beam, denoted G, is defined by:

G (µ, ν) = exp - µ 2 + ν 2 w 2 , ( 1
)
where ω is the waist width and (µ, ν) are the coordinates in the beam waist plane.

The propagation of the Gaussian beam through the pipe, to the CCD sensor is decomposed into two steps. The first step is the propagation of the illuminating beam (i.e G (µ, ν)) from the beam waist plane to the front of the opaque object. Using the generalized Huygens-Fresnel integral [START_REF] Palma | Extension of the Fresnel transform to ABCD systems[END_REF][START_REF] Lambert | Linear systems approach to simulation of optical diffraction[END_REF][START_REF] Yura | Optical beam wave propagation through complex optical systems[END_REF],

the amplitude of the field in the pipe, denoted G 1 , is defined by:

G 1 (ξ, η) = exp i 2π λ E 1 iλ B x 1 B y 1 R 2 G (µ, ν) exp i π λB x 1 A x 1 µ 2 -2ξµ + D x 1 ξ 2 × exp i π λB y 1 A y 1 ν 2 -2ην + D y 1 η 2 dµdν, (2)
1 are given by M x,y 1 matrices of Eq. (38) in App. (A). To obtain the relation of the intensity distribution, the linear canonical transformation is used. For the basic definitions and properties of the linear canonical transformation, we refer to [START_REF] Ozaktas | The Fractional Fourier Transform: with Applications in Optics and Signal Processing[END_REF]. This transformation is a parameterized integral operator of parameters A, B, C and D. Each parameter represents a coefficient of a transfer matrix denoted M (Cf appendix A). In our case the Fresnel transform is considered [START_REF] Goodman | Introduction to Fourier Optics[END_REF][START_REF] Siegman | Lasers[END_REF].

It should be noted that the values of the coefficients are different in both ξ and η direction. This is due to the cylindrical geometry of the pipe. The distance E 1 = z p + e + δ corresponds to the waistobject distance. Using the same formalism, the amplitude in the CCD plane, denoted G 2 (x, y), can be written as:

G 2 (x, y) = exp i 2π λ E 2 iλ B x 2 B y 2 R 2 G 1 (ξ, η) [1 -T (ξ, η)] exp i π λB x 2 A x 2 ξ 2 -2xξ + D x 2 x 2 × exp i π λB y 2 A y 2 η 2 -2yη + D y 2 y 2 dξdη, (3)
where the parameters A x,y 2 , B x,y 2 , D x,y 2 are given by M x,y 2 matrices of Eq. (39). The distance E 2 = z i + e + z is the distance between the object and the CCD sensor.

The spatial transmittance of the opaque 2D-object is defined by [1 -T (ξ, η)]. Here T (ξ, η) can be expressed as a superposition of Gaussian functions [START_REF] Wen | Gaussian beam functions as a base function set for acoustical field calculations[END_REF][START_REF] Wen | A diffraction beam expressed as the superposition of Gaussian beams[END_REF], such as:

T (r) = N k=1 A k exp -r T R T P k Rr , (4)
with

R =     cos θ sin θ -sin θ cos θ     , (5)
and

P k =     B k a 2 0 0 B k b 2     . (6)
This expression permits to deal with non symmetrical optical systems and elliptical opaque objects [START_REF] Zheng | Analytical expression of elliptical Gaussian beams through nonsymmetric systems with an elliptical aperture[END_REF][START_REF] Du | Propagation of decentered elliptical Gaussian beams in apertured and nonsymmetrical optical systems[END_REF][START_REF] Du | Propagation of elliptical Gaussian beams in apertured and misaligned optical systems[END_REF][START_REF] Du | Propagation of elliptical Gaussian beams modulated by an elliptical annular aperture[END_REF]. This Gaussian decomposition is very convenient; as a matter of fact it allows us to establish analytical expression of the diffracted pattern in the CCD plane. It also allows us to simulate holograms of 3D opaque objects whose 2D projection is elliptic or circular [START_REF] Slimani | Near-field Lorenz-Mie theory and its application to microholography[END_REF] (e.g. spheroids, fibers ...).

The object parameters are illustrated on Fig. [START_REF] Malek | Digital in-line holography for threedimensional-two-components particle tracking velocimetry[END_REF]. The coefficients a and b are the object radii within η and ξ axis respectively, θ is the angle between elliptical aperture principal axis and ξ axis.

The A k and B k coefficients are determined by numerical resolution of the Kirchhoff equation [START_REF] Wen | A diffraction beam expressed as the superposition of Gaussian beams[END_REF].

Let R ell = b/a representing the particle ellipticity. Considering R ell = 1 lead to the simulation of a circular particle, whereas R ell = 1 is used to simulate elliptical particles. The particular case of the opaque fiber, parallel to the pipe axis, is obtained when R ell → 0.

In further developments, θ will be considered to be π/2. With this assumption T (ξ, η) becomes:

T (ξ, η) = N k=1 A k exp - B k b 2 ξ 2 + R 2 ell η 2 . (7)
To simulate our opaque objects, N is fixed to 10. From Eq. ( 3), G 2 (x, y) is split into two integrals denoted R (x, y) for the reference beam and O (x, y) for the object beam so that :

G 2 (x, y) = exp i 2π λ E 2 iλ B x 2 B y 2 [R (x, y) -O (x, y)] , (8)
with

R (x, y) = R 2 G 1 (ξ, η) exp i π λB x 2 A x 2 ξ 2 -2xξ + D x 2 x 2 × exp i π λB x 2 A y 2 η 2 -2yη + D y 2 y 2 dξdη, (9)
and

O (x, y) = R 2 G 1 (ξ, η) T (ξ, η) exp i π λB x 2 A x 2 ξ 2 -2xξ + D x 2 x 2 × exp i π λB y 2 A y 2 η 2 -2yη + D y 2 y 2 dξdη. ( 10
)
The functions R (x, y) and O (x, y) are respectively given by:

R (x, y) = exp i 2π λ E 1 iλ B x 1 B y 1 K x 1 K y 1 K x 2 K y 2 × exp - π λ N x B x 2 x 2 + N y B y 2 y 2 exp i π λ x 2 + M y B y 2 y 2 , (11)
and

O (x, y) = exp i 2π λ E 1 iλ B x 1 B y 1 K x 1 K y 1 exp i π λ D x 2 B x 2 x 2 + D y 2 B y 2 y 2 N k=1 A k K xeq 2 K yeq 2 × exp - π λ N xeq B x 2 x 2 + N yeq B y 2 y 2 exp i π λ M xeq B x 2 x 2 + M yeq B y 2 y 2 . ( 12
)
Values of the different parameters of Eqs. ( 11) and ( 12) are defined in App. (B).

After theoretical developments the intensity distribution, denoted I (x, y), recorded by the CCD sensor is:

I (x, y) = G 2 (x, y) G 2 (x, y) = 1 λ 2 B x 2 B y 2 |R| 2 -2ℜ RO + |O| 2 , (13)
where the upper bar denotes the complex conjugate and ℜ the real part. The square modulus |R| ( 4) was shifted by 75 m along x-axis so that numerical and experimental results can be compared. This figure shows the good agreement between numerical and experimental data.

2.B. Case of thin lenses

We have proposed a theoretical model allowing us to deal with propagation of a laser beam through a thick pipe. In a former study, the effect of cylindrical lenses on the diffraction pattern of a particle has been investigated in details, leading to the expression of the intensity distribution in the CCD plane [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]. Using this approach we have a great opportunity to confirm our pipe model results.

In the following developments we aim to compare our thick and cylindrical pipe model with a thin cylindrical lenses approach. The phase transformation due to a thin lens is [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF][START_REF] Goodman | Introduction to Fourier Optics[END_REF]:

Φ (x l , y l ) = exp -i π λ x 2 l f x + y 2 l f y . (14)
Here (x l , y l ) are the coordinates in the lens plane and f x,y are the focal lengths of the lens in both directions. Using Huygens-Fresnel integral [START_REF] Goodman | Introduction to Fourier Optics[END_REF] the intensity distribution in the sensor plane is:

I thin (x, y) = exp i 2π λ E 2 iλE 2 [R thin (x, y) -O thin (x, y)] 2 . ( 15
)
Here, R thin and O thin can be expressed in a quite similar form than amplitude distributions given in Eqs. [START_REF] Onural | Digital decoding of in-line holograms[END_REF] and [START_REF]Diffraction from a wavelet point of view[END_REF].

The transverse intensity distribution profile presented in Fig. [START_REF] De Nicola | Correct-image reconstruction in the presence of severe anamorphism by means of digital holography[END_REF] shows:

lim e→0 I (x, y) = I thin (x, y) , (16)
remembering that e represents the glass thickness. This term is contained in matrices M x,y 1 and M x,y 2 (see Eqs. ( 38) and ( 39)).

For this example: 2b = 51.8 m, a → ∞, z p = 325 mm, n 2 = 1, δ = 18 mm and z = 23 mm. The comparison with other parameter values leads to the same conclusion. As a result, we can consider that our results are consistent with those of Ref. [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] and that our model is versatile enough to deal with various pipe shapes and thickness.

In this section a numerical model allowing to treat thick optical systems has been presented. We now aim to reconstruct the image of the object by means of FRFT from the calculated intensity distribution of the diffraction pattern (Eq. ( 13)). In section 3, mathematical definition of the FRFT is recalled leading to the reconstruction of holograms recorded in pipe systems.

Fractional Fourier transformation analysis of in-line holograms

3.A. Two-dimensional Fractional Fourier transformation

The FRFT is a generalization of the classical Fourier transform. This integral operator has numerous applications in signal processing [START_REF] Ozaktas | The Fractional Fourier Transform: with Applications in Optics and Signal Processing[END_REF]. Its mathematical expression is the following [START_REF] Mcbride | On Namias's fractional Fourier transforms[END_REF][START_REF] Namias | The fractional order Fourier transform and its application to quantum mechanics[END_REF][START_REF] Lohmann | Image rotation, Wigner rotation, and the fractional Fourier transform[END_REF]: the FRFT of order a x = (2α x )/π and a y = (2α y )/π (for x and y cross section respectively), with 0 ≤ |α x | ≤ π/2 and 0 ≤ |α y | ≤ π/2, of a two dimensional function I (x, y) is

F αx,αy [I(x, y)](x a , y a ) = R 2 N αx (x, x a ) N αy (y, y a )I(x, y) dx dy, (17)
where the kernel of the fractional operator is defined by

N αp (x, x a ) = C(α p ) exp iπ x 2 + x 2 a s 2 p tan α p exp - i2πx a x s 2 p sin α p , (18)
and

C(α p ) = exp(-i( π 4 sign(sin α p ) - αp 2 )) |s 2 p sin α p | 1/2 . ( 19
)
Here p = x, y. Generally, the parameter s p is considered as a normalization constant. In our case, its value is defined from the experimental set-up according to [START_REF] Mas | Fast numerical calculation of Fresnel patterns in convergent systems[END_REF]:

s 2 p = N px p • δ 2 p .
The number of samples is N px p in both intensity distribution I (x, y) and fractional domain. δ p is the sampling period along the two axes of the image. C(α p ), which is a function of the fractional order, insures the energy conservation law to be valid. Discretization of the FRFT kernel is performed using an orthogonal projection method proposed by Pei et al. [START_REF] Pei | Discrete Fractional Fourier Transform Based on Orthogonal Projections[END_REF].

3.B. Optimal fractional orders to refocus numerically over the object

In our case, the numerical reconstruction can be considered as a numerical refocusing over the object [START_REF] Verrier | Digital in-line holography with an elliptical, astigmatic Gaussian beam: wide-angle reconstruction[END_REF]. To do that, the quadratic phase, denoted ϕ = arg RO , and contained in the term 2ℜ RO of Eq. ( 13) must be evaluated. This term is composed of a linear chirp modulated by a sum of complex Gaussian functions. Information about the distance between the CCD sensor and the object is carried by the linear chirp (B 2x,y ) whereas information about object size is carried by the modulation (M x,y eq and N x,y eq ). The optimal reconstruction of the image consists in compensating the quadratic phase terms of the intensity distribution [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF].

The quadratic phase can be determined from Eqs. ( 11) and ( 12)

ϕ = π λ M x -D x 2 B x 2 x 2 + M y -D y 2 B y 2 y 2 , (20)
thus ℜ RO can be written as:

ℜ RO = RO cos (iϕ) . (21)
The quadratic phase term contained in the FRFT kernel, denoted ϕ a , is given by

ϕ a = π cot α x s 2 x x 2 + cot α y s 2 y y 2 . ( 22
)
The image reconstruction is obtained by applying the FRFT to the intensity distribution of Eq. ( 13):

F αx,αy [I(x, y)] ∝ F αx,αy [|R| 2 + |O| 2 ] -2F αx,αy [ RO cos ϕ]. (23)
The terms |R| 2 and |O| 2 contain no linear chirps, thus they do not influence the optimal fractional orders of reconstruction to be determined. Only the second term is useful for image reconstruction.

By noting that 2 cos ϕ = exp (-iϕ) + exp (iϕ), Eq. ( 23) becomes:

F αx,αy [I(x, y)] ∝ F αx,αy [|R| 2 + |O| 2 ] -C (α x ) C (α y ) R 2 RO exp [i (ϕ a -ϕ)] exp -i2π x a x s 2 x sin α x + y a y s 2 y sin α y dxdy -C (α x ) C (α y ) R 2 RO exp [i (ϕ a + ϕ)] exp -i2π x a x s 2 x sin α x + y a y s 2 y sin α y dxdy. ( 24
)
The best hologram reconstruction is reached when one of the quadratic phase term is brought to zero. Thus reconstruction is performed if

ϕ a ± ϕ = 0. ( 25
)
The optimal fractional orders α opt x and α opt y are defined from Eqs. ( 20), ( 22) and ( 25) and take the values:

α opt x = arctan ∓ B x 2 λ s 2 x (M x -D x 2 ) , α opt y = arctan ∓ B y 2 λ s 2 y (M y -D y 2 ) . (26)
It should be noted that we have checked that if we consider M 2x,y to be free space propagation, leading to B x,y 2 = z i + e + z and D x,y 2 = 1, the expression of the optimal fractional order of Ref. [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] are recovered.

3.C. Numerical simulations

Experimental set-up of Fig.

The reconstruction of the fiber with FRFT can be seen on Fig. [START_REF] Grilli | Whole optical wavefields reconstruction by Digital Holography[END_REF]. As we can see on this figure, the reconstructed image is disturbed by background fringes. This effect is due to the in-line configuration and is commonly known as the twin image effect [START_REF] Gabor | Microscopy by reconstructed wave-fronts[END_REF].

3.D. Experimental results

Since now, we have presented a method to simulate holograms in thick optical systems. Reconstruction has been successfully performed thanks to FRFT. In order to validate theoretical developments and simulations, a glass pipe with curvatures Rx 1 , Ry 27). To perform FRFT on this experimental image, FRFT orders have been adjusted to obtain the best image of the fiber i. e the best contrast between the reconstructed fiber and the background. Thus doing, the accuracy on the fractional order is approximately 10 -2 .

Reconstruction is presented on Fig. [START_REF] De Nicola | Wave front reconstruction of Fresnel off-axis holograms with compensation of aberrations by means of phase-shifting digital holography[END_REF], estimated optimal fractional orders are a x = 0.33 and a y = 0.27. These values are very close to the theoretical ones.

This good accordance is confirmed by the curves of Fig. [START_REF] Crane | Far-field holography of ampule contaminants[END_REF]. Here theoretical values of a opt x (solid line) and a opt y (dashed line) are plotted versus the distance between the pipe and the CCD sensor.

Theoretical values corresponds to an air filled glass pipe, δ is fixed to 10 mm and z varies from 0 to 140 mm. Note that, when z = 0, a x,y = 0, this is due to the fact that, in this case, z represents the distance between the pipe and the CCD sensor instead of the distance between the object and the CCD sensor which is considered in [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]. Experiments performed with this configuration allow us to plot estimated values of a opt x (circles) and a opt y (diamonds) on the same graph. As shown by Fig.

( 9), estimated and theoretical values are closely linked. Theoretical representation hereby presented is therefore well adapted to these optical systems.

Until then, we have performed reconstruction with objects for which position is well known. We now apply our formalism to reconstruct hologram of latex beads. The intensity distribution of the diffracted field is represented on Fig. [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]. Here, 100 m in diameter latex beads are used, z p = 325 mm (i.e. propagation distance between the source and the pipe) and z = 23 mm (i.e. distance between the pipe and the CCD sensor). This study gives us the opportunity to compare reconstruction using our novel approach with reconstruction using a thin lens approach. As far as every parameter of the experimental set-up is known, excepting δ, the reconstruction process is the following: for δ ∈ [0, 36 mm] we calculate the optimal FRFT orders using Eq. ( 26). Doing so, we are able to refocus on the particles which are in the pipe.

In Fig. [START_REF] De Nicola | Correct-image reconstruction in the presence of severe anamorphism by means of digital holography[END_REF], it is shown that when the pipe thickness e is close to zero, the two approaches are equivalent. Thus, to reconstruct the hologram of Fig. [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] with the thin lens approach, we only need to use the orders given in Eq. ( 26) with e = 0.

Comparison between the two reconstruction methods can be made with Fig. [START_REF] Onural | Digital decoding of in-line holograms[END_REF].

Here, the reconstruction is realized with the same value of δ in both cases. We can notice that refocusing on the particle is impossible with the thin lens approach, whereas our thick lens approach allows a good reconstruction of the image of particles.

Therefore, the thick lens formalism is well adapted to pipe flow studies: as we are able to refocus on particles in a pipe, metrologies of particles diameter is possible.

Conclusion

In this paper, an analytical expression of the scalar diffraction, under Fresnel approximations, in thick optical systems such as pipes is derived. As such optical systems reveals astigmatism, FRFT is a reliable tool to perform hologram reconstruction. Optimal fractional orders have been calculated leading to a satisfactory reconstruction of either numerical or experimental images. Numerical ex-periments have been performed showing a good agreement with experimental results. Comparison between thin lens and thick lens approaches has been performed showing that the thin lens formalism is not well adapted to the study of pipe flows.

A. Appendix A : Transfer matrices of the optical system Each part of an optical system can be represented by a matrix. Using this principle, we can build a set of matrices corresponding to our experiments (Fig. ( 1)). After the beam shaping step, the beam propagates in free space over a distance z p . The associated M zp matrix is:

M zp =     1 z p 0 1     . (28)
Then the beam is refracted at the interface between free space and glass (n 1 = 1.5). The curvatures along both directions are Rx 1 and Ry 1 . M Rx1 , M Ry1 matrices along x and y axis are:

M Rx1 =     1 0 n0-n1 Rx1 1     , M R1y =     1 0 n0-n1 Ry1 1     . (29)
After refraction, the beam propagates in glass (n 1 =1.5). For this step:

M e =     1 e n1 0 1     , (30)
Next step is the refraction of the beam at the interface between glass and the medium inside the pipe (refractive index n 2 ). Curvature are given by Rx 2 and Ry 2 in x and y direction.

M Rx2 =     1 0 n1-n2 Rx2 1     , M Ry2 =     1 0 n1-n2 Ry2 1     , (31)
We are now in the pipe. To reach the object, we have to propagate over δ.

M δ =     1 δ n2 0 1     . (32)
Doing the same over z i permit to reach the output of the pipe:

M zi =     1 zi n2 0 1     . (33)
After a refraction at the interface (curvature Rx 3 and Ry 3 ) between the medium with refractive index n 2 and glass:

M Rx3 =     1 0 n2-n1 Rx3 1     , M Ry3 =     1 0 n2-n1 Ry3 1     , (34)
a propagation in glass (n 1 =1. 5)

M e =     1 e n1 0 1     , (35)
a refraction at the interface (curvature Rx 4 and Ry 4 ) between glass (n 1 =1.5) and free space:

M Rx4 =     1 0 n1-n0 Rx4 1     , M Ry4 =     1 0 n1-n0 Ry4 1     , (36)
and free space propagation over z:

M z =     1 z 0 1     , (37)
the whole ABCD system is described.

B. Appendix B : Amplitude distributions R (x, y) and O (x, y)

In Appendix A each part of our optical system has been represented with transfer matrices. Using this formalism allows us, under paraxial conditions, to deal with propagation of a Gaussian point source through the pipe.

Intensity distribution of the diffraction pattern in the CCD sensor plane is determined by considering two matrix systems: M x,y 1 and M x,y 2 .

M x,y 1 is composed of three steps : propagation in free space over z p , propagation through the first thick lens, propagation in a medium of refractive index n 2 over δ. It is characterized by two transfer matrices:

M x,y 1 = M δ × M x,y L1 × M zp =     A x,y 1 B x,y 1 C x,y 1 D x,y 1     . (38)
is also composed of three steps : propagation in a medium of refractive index n 2 over z i , propagation through the second thick lens, propagation in free space over δ. Transfer matrices for this system are:

M x,y 2 = M z × M x,y L2 × M zi =     A x,y 2 B x,y 2 C x,y 2 D x,y 2     . (39)
Thanks to Eqs. ( 29), [START_REF] Mas | Fast numerical calculation of Fresnel patterns in convergent systems[END_REF], [START_REF] Pei | Discrete Fractional Fourier Transform Based on Orthogonal Projections[END_REF] one can build the transfer matrices of the first thick lens:

M x L1 = M Rx2 × M e × M Rx1 , M y L1 = M Ry2 × M e × M Ry1 . (40)
By using the same method for the second lens, we obtain:

M x L2 = M Rx4 × M e × M Rx3 , M y L2 = M Ry4 × M e × M Ry3 . (41) B.A. Propagation through M x,y 1
After analytical developments of Eq. ( 2), the complex amplitude distribution in the object plane is:

G 1 (ξ, η) = exp i 2π λ E 1 iλ B x 1 B y 1 K x 1 K y 1 exp - ξ 2 ω 2 1x + η 2 ω 2 1y exp -i π λ ξ 2 R 1x + η 2 R 1y , (42)
with K x,y 1 given by:

K x,y 1 = πω 2 1 -iA x,y 1 πω 2 λB x,y 1 1/2 , (43)
ω 1x,y and R 1x,y are respectively the beam radii and the wavefront curvature in the particle plane.

Their mathematical expressions are:

ω 1x,y = λB x,y 1 πω 1 + A x,y 1 πω 2 λB x,y 1 2 1/2 , R 1x,y = - B x,y 1 D x,y 1 - A x,y 1 πω 2 λB x,y 1 2 1+ A x,y 1 πω 2 λB x,y 1 2 .
(44)

B.B. Propagation through M x,y 2
Propagation to the CCD sensor plane is calculated thanks to the generalized Huygens-Fresnel integral (see Eq. ( 3)).

B.B.1. Amplitude distribution R (x, y)
R is associated with the reference wave. After analytical developments of Eq. ( 9), the amplitude distribution R is:

R (x, y) = exp i 2π λ E 1 iλ B x 1 B y 1 K x 1 K y 1 K x 2 K y 2 × exp - π λ N x B x 2 x 2 + N y B y 2 y 2 exp i π λ M x B x 2 x 2 + M y B y 2 y 2 , (45)
where,

M x,y = D x,y 2 + πω 2 1x,y λB x,y 2 2 B x,y 2 R1 x,y -A x,y 2 1 + πω 2 1x,y λB x,y 2 2 B x,y 2 R1 x,y -A x,y 2 2 , N x,y = πω 2 1x,y λB x,y 2 1 + πω 2 1x,y λB x,y 2 2 B x,y 2 R1 x,y -A x,y 2 2 , (46)
and

K x,y 2 =    πω 2 1x,y 1 + i πω 2 1x,y λB x,y 2 B x,y 2 R1 x,y -A x,y 2    1/2 . ( 47
) B.B.2. Amplitude distribution O (x, y)
O is the amplitude of the diffracted wave. We define ω 1x,y eq and R 1x,y eq

1 ω 2 1x eq = 1 ω 2 1x + ℜ{B k } b 2 , 1 ω 2 1y eq = 1 ω 2 1y + R 2 ell ℜ{B k } b 2 , (48) and 1
R 1x eq = 1 R 1x + ℑ{B k }λ πb 2 , 1 R 1y eq = 1 R 1y + R 2 ell ℑ{B k }λ πb 2 , (49)
to simplify notations. It should be noted that ℜ and ℑ stand for real and imaginary part respectively.

Thus O becomes:

O (x, y) = exp i 2π λ E 1 iλ B x 1 B y 1 K x 1 K y 1 exp i π λ D x 2 B x 2 x 2 + D y 2 B y 2 y 2 N k=1 A k K

2

  corresponds to the directly transmitted beam whereas |O| 2 is associated with the diffracted part of the beam.To illustrate Eq. (13), the optical set-up given on Fig. (1) is considered. We investigate the diffraction pattern obtained with an opaque fiber parallel to the pipe axis. From the Gaussian decomposition of the object function, opaque fibers, parallel to the pipe axis, can be obtained taking R ell → 0. The value of the opaque fiber diameter is 2b = 51.8 m and a → ∞. The glass-made (n 1 = 1.5) pipe is filled with water of refractive index n 2 = 1.33. The beam propagates in free space over z p = 325 mm. The dimensions δ and z are fixed to 18 mm and 23 mm respectively. The image of Fig. (3) illustrates the intensity distribution of the diffracted field. The simulation is carried out by calculating the gray level of each pixel using Eq. (13). The size of the hologram is 768×576 pixels, and pixel pitch is 11 m. To validate our simulation, an experiment, which consists in placing an opaque fiber in the pipe, was performed. The hologram of Fig (4) represents the intensity of the diffraction recorded with the CCD sensor using the same parameters than previously. This illustration reveals a good accordance between numerical and experimental diffraction patterns. To confirm this point, the transverse intensity profiles obtained from Fig. (3) and Fig. (4) are presented on Fig (5). Here the normalized intensity I * is plotted against the x-axis. The intensity profiles are calculated by cumulating the gray levels along the η-axis (direction of the fiber) over 50 rows, which are figured out by the rectangular selection on Fig. (4). The diffraction pattern of Fig.

( 1 )

 1 is used to perform simulations. An application of the model is illustrated on Fig.[START_REF] Malkiel | The three-dimensional flow field generated by a feeding calanoid copepod measured using digital holography[END_REF]. Recall that we consider an opaque fiber (R ell → 0) 2b = 51.8µm in width.The glass made pipe is filled with water (n 2 = 1.33). δ and z are set to 18 mm and 23 mm respectively. Rx 1 , Ry 1 , Ry 2 → ∞ and Rx 2 = 18 mm are the first thick lens curvatures along x-and y-axis. The curvature radii of the second thick lens have the same values but opposite signs. With these parameters and owing to Eq. (26) the optimal fractional orders within this configuration are a opt x = 0.337, a opt y = 0.273.

Fig. 1 .Fig. 2 .

 12 Fig. 1. a) Schematic representation of the optical set-up (not to scale). Definition of the

Fig. 3 .

 3 Fig. 3. Simulation of the diffraction pattern of a 51.8 µm opaque fiber, parallel to the axis

Fig. 4 .Fig. 5 . 1 -Fig. 6 .Fig. 7 .

 45167 Fig. 4. Experimental diffraction pattern of a 51.8 µm opaque fiber, parallel to the axis of a

Fig. 8 .Fig. 9 .

 89 Fig. 8. Reconstruction of the fiber image from the diffraction pattern of Fig. (4) by FRFT with a opt x = 0.33 and a opt y = 0.27.

Fig. 10 .

 10 Fig. 10. Experimental diffraction pattern of 100m latex beads. z p = 325 mm, z = 23 mm.

Fig. 11 .

 11 Fig. 11. Reconstruction of the latex beads image. (a) Reconstruction with thin lens param-

  1 , Ry 2 → ∞ and Rx 2 = 18 mm has been used for the experiments. Curvatures on the other side of the pipe are deduced from Rx 1 , Ry 1 , Rx 2 , Ry 2 by taking negative values.

The image of Fig.

[START_REF] Garcia-Sucerquia | Digital in-line holographic microscopy[END_REF]

represents the intensity of the diffraction pattern recorded with a 768 × 576 px 2 CCD sensor with 11 µm pitch. The object is a 2b = 51.8µm opaque fiber, δ and z are approximately and respectively equal to 18 mm and 23 mm. Theoretical orders associated with this experiment are given by Eq. (

y 2 , (50)

Acknowledgments

These works come within the scope of Hydro-Testing Alliance (HTA) with the support of the European Commission's Sixth Framework Programme under DG Research.

List of Figures


N Verrier

S Cotmellec
email: coetmellec@coria.fr

M Brunel

D Lebrun

A J E M Janssen

Digital in-line holography with an elliptical, astigmatic Gaussian beam : wide-angle reconstruction

We demonstrate in this paper that the effect of object shift in an elliptical, astigmatic Gaussian beam does not affect the optimal fractional orders used to reconstruct the holographic image of a particle or another opaque object in the field. Simulations and experimental results are presented.

Introduction

Digital in-line holography (DIH) is widely used in microscopy for biological applications [START_REF] Xu | Tracking particles in four dimensions with in-line holographic microscopy[END_REF][START_REF] Dubois | Improved three-dimensional imaging with digital holography microscope using a partial spatial coherent source[END_REF], in 3D Holographic Particle Image Velocimetry (HPIV) for fluid mechanics studies [START_REF] Skarman | Simultaneous 3D-PIV and temperature measurement using a new CCD based holographic interferometer[END_REF] and in refractometry [START_REF] Sebesta | Object characterization with refractometric digital Fourier holography[END_REF]. In most theoretical DIH studies, the optical systems or the objects used, for example particles, are considered to be centered on the optical axis [START_REF] Schnars | Direct phase determination in hologram interferometry with use of digitally recorded holograms[END_REF][START_REF] Shen | Digital holography particle image velocimetry for the measurement of 3Dt-3c flows[END_REF]. Nevertheless, in many practical applications, the systems and objects are not necessarily centered. This is not a problem under plane wave illumination. But more and more studies involve a diverging beam, or astigmatic beams directly 1 obtained at the output of fibered laser diodes or astigmatic laser diodes. The position of the object in the beam becomes particularly relevant. To reconstruct the image of an object, a reconstruction parameter must be determined in order to find the best focus plane. For wavelet transformation, the parameter is the scale factor [START_REF]Diffraction from a wavelet point of view[END_REF]. For Fresnel transformation, the parameter is the distance z from the object to the quadratic sensor plane [START_REF] Onural | Digital decoding of in-line holograms[END_REF] and for the fractional Fourier transformation, the parameters are the fractional orders [START_REF] Fogret | Agreement of fractional Fourier optics with the Huygens Fresnel principle[END_REF][START_REF] Cotmellec | Characterization of diffraction patterns directly from in-line holograms with the fractional Fourier transform[END_REF]. It is considered here that the same parameter can be used for all positions of the object. Recently, an analytical solution of scalar diffraction of an elliptical and astigmatic Gaussian beam (EAGB) by a centered opaque disk under Fresnel approximation has been proposed. By using the fractional Fourier transformation, a good particle image reconstruction is obtained [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]. However, there are no recent publications in DIH theory demonstrating that reconstruction parameter is constant for all transverse positions of the object in the field of the beam. This is, however, particularly important if we want to develop wide-angle metrologies.

In this publication, the aim is to demonstrate that the same fractional orders can be used to reconstruct an image of the particle whatever its position in the field of the beam. We exhibit the effect of the Gaussian beam on the reconstructed image. In the first part of this publication, the model of the analytical solution of scalar diffraction of an EAGB by a centered opaque disk is revisited to take into account a decentered object. In the second part, the definition of the fractional Fourier transformation is recalled and this transformation is used to reconstruct the image of the particle. It is in this part that we demonstrate that the same orders can be used. Finally, we propose to illustrate our results by simulations and experimental results.

In-line Holography with an elliptic and astigmatic Gaussian beam

The basic idea in DIH is to record by a CCD camera the intensity distribution of the diffraction pattern of an object illuminated by a continuous or pulsed wave [START_REF] Nicolas | Digital in-line holography with a sub-picosecond laser beam[END_REF]. Figure [START_REF] Xu | Tracking particles in four dimensions with in-line holographic microscopy[END_REF] represents the numerical and experimental set-up where all parameters are identified. The incident Gaussian beam of diameter ω crosses a plano-convex cylindrical lens. This cylindrical lens acts only in the η-axis. Its focal length is equal to f η = 200mm. After propagation in free space, the elliptical and astigmatic Gaussian beam illuminates an opaque particle. In the plane of the object, the beam widths along the ξ-axis and η-axis are defined by ω ξ and ω η . The wavefront curvatures of the astigmatic beam are denoted by R q with q = ξ, η. The CCD camera is located at a distance z from the opaque particle.

The basic model to describe the intensity distribution of the diffracted beam by an object recorded by the CCD camera is the integral of Kirchhoff-Fresnel given by the scalar integral for the complex amplitude A in the quadratic sensor plane:

A = exp(i 2π λ z) iλz R 2 E T (ξ, η) exp iπ λz (ξ -x) 2 + (η -y) 2 dξdη, (1)
in which E T (ξ, η) is the product of the optical incident beam, denoted E(ξ, η), by the spatial transmittance of the shifted opaque 2D-object, denoted 1 -T (ξ -ξ 0 , η -η 0 ). The quadratic sensor records the intensity defined by |A| 2 . If we consider that the function E T is the product of an elliptic and astigmatic Gaussian beam by an opaque disk of diameter D, then

E T (ξ, η) = exp c ξ ξ 2 + c η η 2 =E(ξ,η) • [1 -T (ξ -ξ 0 , η -η 0 )] shif ted object . (2)
The complex coefficients c ξ and c η are

c ξ = - 1 ω 2 ξ -i π λR ξ , c η = - 1 ω 2 η -i π λR η , (3)
For an opaque disk centered at the origin O, the transmittance function T (ξ, η) in the object plane is:

T (ξ, η) =                        1, 0 < ξ 2 + η 2 < D/2, 1/2, 0 < ξ 2 + η 2 = D/2, 0, ξ 2 + η 2 > D/2 > 0. ( 4
)
¿From Eqs. ( 1) and ( 2), the expression for A(x, y) can split into two integral terms, denoted A 1 and

A(x, y) = exp(i 2π λ z) i [A 1 -A 2 ] . (5)
The expression of these two terms is given in the Appendix (A). The complex amplitude A 1 represents the propagation of the incident beam without diffraction by the particle and A 2 contains the diffraction by a pinhole of diameter D. After some development detailed in the Appendix (A), the expression of A 2 becomes:

A 2 = πD 2 λz exp [Φ(ξ 0 , η 0 )] • exp ic z (x -ξ 0 ) 2 + (y -η 0 ) 2 • ∞ k=0 (-i) k ε k T k (r, γ) cos(2kθ), (6)
with

T k (r, γ) = ∞ p=0 β 2k 2k+2p (δ)V 2k 2k+2p (r, γ), (7)
and

V 2k 2k+2p (r, γ) = exp(iγ/2) ∞ m=0 (2m + 1)i m jm(γ/2) • m+p l=max(0,m-2k-p,p-m) (-1) l ω ml J 2k+2l+1 (r) r . (8)
The coefficients ω ml are given explicitly in [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] and [START_REF] Janssen | On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus[END_REF]. Finally, the expression of A 1 contains only the characteristics of the incident beam and A 2 contains a shifted linear chirp function linked to the object shift in the EAGB. The previous function is modulated by a series of Bessel functions which constitutes the envelope of the amplitude distribution of A 2 .

A similar expression has been established in [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] for a particle located on the axis of the beam, i.e. ξ 0 = 0 and η 0 = 0. This last expression is more general and can be applied to a particle located everywhere in the field of the beam.

In our previous study, we demonstrated that a fractional-order Fourier transformation allows a particle located on the axis of the beam to be reconstructed. Unfortunately, the expression of the diffraction pattern (and particularly A 2 ) is so complex that it cannot be proved theoretically that a digital (or optical) reconstruction leads effectively to an opaque disk function. The demonstration is empirical. In addition, nothing proves that the fractional orders used for the reconstruction would be appropriate if the object were shifted transversely.

to the object. We will then demonstrate that the reconstruction process is successful whatever the transversal position of the particle in the field of the beam. In addition, we show that the same fractional orders have to be applied. We will thus be able to conclude that reconstruction is possible for a whole wide-field object or wide particle field in such a beam. Now, to obtain the desired accuracy, we should analyze the number of terms necessary in the series over k in Eq. ( 6) and over p in Eq. ( 7). The series which should be analyzed are in Eq.( 8).

The first upper bound that is relevant here is [14-16]

j m γ 2 ≤ 1 (2m + 1) 1/2 min 1, π 2 1/2 |γ/4| m m! . (9)
The coefficients ω ml verify ω ml ≥ 0 and l ω ml = 1, so the ω's are completely innocent. Using the bound ( 9) for the function j m γ 2

requires a rough estimate of the variable γ. In practical experiments, we have

D ≈ 10 -4 m, λ ≈ 10 -6 m, z ≈ 10 -1 m, R q ≈ 0.5 • 10 -1 then γ ≈ 0.785 • 10 -1 + 0.816 • 10 -3 i.
The right-hand side of Eq. ( 9) is less than 0.142 • 10 -1 for m ≥ 1. The second upper bound we use concerns the Bessel functions J 2k+2l+1 . From [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], 9.1.62 on p.362 and [START_REF] Janssen | On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus[END_REF] we obtain:

J 2k+2l+1 (r) r ≤ min 1, 1 2 |r| 2k+2l exp(|r| 2 ) (2k + 2l + 1)! . ( 10
)
The evaluation of the accuracy concerns the product of the left-hand side of Eqs. ( 9) and ( 10):

j m γ 2 J 2k+2(m+p)+1 (r) r . (11)
With the previous values, all the quantities in [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] are less than 0.142 • 10 -1 for all (k, p, m) ≥ 1.

Thus, we only consider the case where k = p = m = 0. In conclusion, the function V 0 0 (r, γ) is given with good accuracy by:

V 0 0 (r, γ) ≃ exp(iγ/2)j0(γ/2) J1(r) r . ( 12
)
Note that this result is true if the object is far from the beam waist. The amplitude A 2 then becomes

A 2 = πD 2 2λz β 0 0 (δ) exp [Φ(ξ 0 , η 0 )] exp ic z (x -ξ 0 ) 2 + (y -η 0 ) 2 • V 0 0 (r, γ). (13)

2.A. Intensity distribution of the diffraction pattern

The intensity distribution of the diffraction pattern in the quadratic sensor plane, denoted I, is evaluated from the Eqs. ( 5), ( 13) and (36) in the following way:

I = AA = [A 1 -A 2 ] A 1 -A 2 = |A 1 | 2 + |A 2 | 2 -2 ℜ A 1 A 2 , (14)
where the overhead bar denotes the complex conjugate, ℜ denotes the real part. Thus, the intensity distribution recorded by the CCD sensor is described by the Eq. ( 14). As one sees from the form of I, the first and second terms, i.e., |A 1 | 2 and |A 2 | 2 , do not generate interference fringes with a linear instantaneous frequency, denoted f i (x), in the CCD plane, i.e. [START_REF] Mecklenbräuker | The Wigner Distribution[END_REF]:

f i (x) = 1 2π ∂ϕ(x) ∂x = 1 2π ∂ arg(|A 1,2 |) ∂x = 0. ( 15
)
But the third term exhibits a phase which is composed of a constant and a linear instantaneous frequency. This fact is important because the fractional Fourier transform is an effective operator for analyzing a signal containing a linear instantaneous frequency (linear chirp functions). From Eq. ( 14) we write

A 1 A 2 = |A 1 A 2 | exp i arg A 1 A 2 , (16)
where arg

A 1 A 2 = φ -φ 0 with φ = c z x 2 (M ξ -1) + y 2 (M η -1) + 2c z (xξ 0 + yη 0 ) -arg J 1 (r) r , (17)
and

φ 0 = ℜ(γ) 2 + ℑ(Φ(ξ 0 , η 0 )) + arg j 0 γ 2 + arg β 0 0 (δ) -arg (K ξ K η ) , (18)
where ℑ represent the imaginary part of a complex number. The first term in (17) yields a quadratic phase and the second term yields a linear phase. Remember that the aim of the reconstruction by means of FRFT is precisely to analyze a linear chirp.

To give two different examples, it is necessary to fix the values of the parameters (ω ξ , ω η ), (R ξ , R η ) and (D, λ, z).

In the first case, the values are defined by (7mm, 1.75mm) for the beam waists, (-∞, -50mm)

for the wave's curvatures and the diameter D of the particle is equal to 150m and located at 120mm from the CCD sensor. The wavelength of the laser beam is 632.8nm. The distance between the cylindrical lens and the particle is δ = 250mm. The particle is shifted from the origin by (ξ 0 , η 0 ) = (0.5mm, 0.2mm). Figure [START_REF] Dubois | Improved three-dimensional imaging with digital holography microscope using a partial spatial coherent source[END_REF] illustrates the diffraction pattern which is recorded by the camera. Note that the shift of the diffraction pattern (x 0 , y 0 ) observed in the camera plane is not equal to the object shift (ξ 0 , η 0 ). If the particle is considered far from the beam waist then we have the formula:

y 0 = | ∆ | ±z | ∆ | η 0 . ( 19
)
The sign of z depends on the position of the particle compared with the position of the waist. If the particle is after the waist, the sign is positive. If it is in front of the waist, the sign is negative. Along

x-axis, the parameter ∆ is infinite so that in the plane of the camera x 0 = ξ 0 = 0.5mm and along y-axis, ∆ = 50mm thus y 0 = 0.68mm. Now, in the second case, the values are defined by (7mm, 1.75mm) for the beam waists, (-∞, 50mm) for the wave's curvatures and the diameter D of the particle is equal to 150m and located at z = 120mm from the CCD sensor. The distance between the cylindrical lens and the particle is δ = 150mm. The particle is shifted from the origin by (ξ 0 , η 0 ) = (0.5mm, 0.2mm). Figure

(3) illustrates the diffraction pattern which is recorded by the camera. In the plane of the camera,

x 0 = ξ 0 = 0.5mm, ∆ = 50mm thus, from Eq. ( 19), it leads to y 0 = -0.28mm. The diffraction pattern changes from elliptical fringes to hyperbolic fringes. These diffraction patterns will be used to reconstruct the image of the particle by FRFT.

Fractional Fourier transformation analysis of in-line holograms

3.A. Two-dimensional Fractional Fourier transformation

FRFT is an integral operator that has various application in signal and image processing. Its mathematical definition is given in Ref. [START_REF] Mcbride | On Namias's fractional Fourier transforms[END_REF][START_REF] Namias | The fractional order Fourier transform and its application to quantum mechanics[END_REF][START_REF] Lohmann | Image rotation, Wigner rotation, and the fractional Fourier transform[END_REF]. The two-dimensional fractional Fourier transformation of order a x for x-cross-section and a y for y-cross-section with 0

where the kernel of the fractional operator is defined by

N αp (x, x a ) = C(α p ) exp iπ x 2 + x 2 a s 2 p tan α p exp - i2πx a x s 2 p sin α p , (21)
and

C(α p ) = exp(-i( π 4 sign(sin α p ) - αp 2 )) |s 2 p sin α p | 1/2 . ( 22
)
Here p = x, y. Generally, the parameter s p is considered as a normalization constant. It can take any value. In our case, its value is defined from the experimental set-up according to [START_REF] Mas | Fast numerical calculation of Fresnel patterns in convergent systems[END_REF]

s 2 p = N p • δ 2 p . ( 23
)
This definition is presented in the Appendix C. N p is the number of samples along the x and y axes in both spatial I(x, y) and fractional domains. The constant δ p is the sampling period along the two previous axes of the image. In our case the number of samples and the sampling period are the same along both axes, so that the parameters s p are equal to s. The energy-conservation law is ensured by the coefficient C(α p ) which is a function of the fractional order. One of the most important FRFT feature is its ability to transform a linear chirp into a Dirac impulse. Let g(x) = exp iπχx 2 the chirp function to be analyzed. If we consider 1 χ = -tan α, and using the fact that:

lim ε→0 1 √ iπε exp - x 2 iε = δ(x), (24)
F α [g(x)] = δ(x a ). (25)
Finally, by choosing the adequate value of the fractional order, a pure linear chirp function can be transformed into a delta Dirac distribution.

3.B. Reconstruction: optimal fractional orders

To reconstruct the image of the particle, the FRFT of the diffraction pattern (Eq.( 14)) must be calculated:

F αx,αy [I] = F αx,αy |A 1 | 2 -F αx,αy 2|A 1 A 2 | cos(φ -φ 0 ) + F αx,αy |A 2 | 2 (26)
The terms |A 1 | 2 and |A 2 | 2 do not contain a linear chirp, so they do not have any effect on the optimal fractional order to be determined. But the second term, denoted S t contains a linear chirp.

It will be considered for the image reconstruction of the particle. By noting that 2 cos(φ -φ 0 ) = exp(-i(φ -φ 0 )) + exp(i(φ -φ 0 )), the second term of Eq. ( 26) becomes :

F αx,αy 2|A 1 A 2 | cos(φ -φ 0 ) = exp iπ x 2 a s 2 tan α x exp iπ y 2 a s 2 tan α y {I -+ I + } (27)
with

I ± = C(α x )C(α y ) R 2 A 1 A 2 exp [i (φ a ± (φ -φ 0 ))] exp - 2iπ s 2
x a x sin α x + y a y sin α y dxdy (28)

The quadratic phase term of the FRFT is denoted by φ a = π s 2 x 2 cot α x + y 2 cot α y . Let us recall that the FRFT allows us to analyze a linear chirp. Thus, if the fractional orders check the following conditions

π cot α opt x s 2 = c z (M x -1), π cot α opt y s 2 = c z (M y -1), (29)
then, the FRFT of I -is no more than a classical Fourier transformation according to

I -= χ • F J 1 (r) r • exp - π λz ρ T N ρ (u, v), ( 30
) with χ = πD 2 2λz C(α x )C(α y )K ξ K η exp [Φ(ξ 0 , η 0 )] β 0 0 (δ) exp i γ 2 j 0 γ 2 .
The operator F is the 2D-Fourier transformation. The spatial frequencies u and v are equal to:

u = x a s 2 sin(α opt x ) + c z ξ 0 π v = y a s 2 sin(α opt y ) + c z η 0 π (31)
As the Fourier transform of the product of two functions is equal to the convolution of their transforms then

I -= χ • F J 1 (r) r * F exp - π λz ρ T N ρ , (32)
Remember here that the variables ρ and r pertain to the same coordinates (x, y). With the shift theorem for the Fourier transform, the Hankel transform and the discontinuous Weber-Schafheitlin integral in [ [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]], 11.4.42 on p.487, we have:

F J 1 (r) r = 2π λz πD 2 exp [-i2π(uX 0 + vY 0 )] ×                        1, 0 < √ u 2 + v 2 < D/2 λz , 1/2, 0 < √ u 2 + v 2 = D/2 λz , 0, √ u 2 + v 2 ¿ D/2 λz > 0, (33)
with X 0 = ξ 0 (1 -ic ξ /c z ) and Y 0 = η 0 (1 -ic η /c z ). The function defined by the right-hand side of Eq. ( 33), with spatial coordinates (x a , y a ), has the aperture of the pinhole with diameter equal to the diameter D of the opaque particle. We have thus demonstrated our result: reconstruction with a FRFT leads exactly to the object. In addition, shifting the object does not modify the fractional order. This point is important because in the case of a particle field, a single fractional order couple along the x-axis and y-axis is necessary to reconstruct the particle image. If one wishes to determine the shift of the diffraction patterns in the (x a , y a )-plane, the coordinates that should be considered are:

s 2 (u - c z ξ 0 π ) tan α opt x , s 2 (v - c z η 0 π ) tan α opt y ( 34
)
This correction is necessary because the fractional Fourier transformation is not invariant by translation: a space-shift in the spatial domain will lead to both frequency-shift and phase-shift in the

F α [f (x -b)](u) = exp iπb 2 sin α cos α exp (-i2πbu sin α) F α [f (x)](u -b cos α) (35)
Note that the nature of the Gaussian beam implies that the reconstructed image of the object is convolved by a 2D Gaussian function. The sign of the object function (opaque object with a luminous background) is retrieved by applying an inversion of I -that is realized by the minus in front of the second term of Eq. (26).

3.C. Numerical experiments

The simulations of the particle image reconstruction are realized from the diffraction patterns illustrated by the Figs. FRFT is successful. Note that we have checked that the apertures resulting from Eq. ( 32) and the image of the reconstructed particle give the same diameter D.

3.D. Experimental results

As the reconstruction process is successful whatever the transversal position of the particle in the field of the beam, for the same fractional orders we can now carry out a reconstruction for a wide-field object.

The previous theoretical developments and numerical experiments have been tested by using an RS-3 standard reticle (Malvern Equipment). This reticle is an optical glass plate with a pattern of the word "ELECTRO" photographically deposited on the surface. The word "ELECTRO" spread over 6.5mm. The reticle is located at δ ≈ 156mm from the cylindrical lens (CL). The distance z between the CCD camera and the reticle is approximatively equal to 117mm. The opaque word "ELECTRO" is in front of the waist of the beam. The intensity distribution of the EAGB diffracted by the word "ELECTRO" is shown in Fig. [START_REF] Shen | Digital holography particle image velocimetry for the measurement of 3Dt-3c flows[END_REF]. The image of the word "ELECTRO" is reconstructed by the fractional Fourier transformation. Fractional orders have been adjusted to obtain the best contrast between the reconstructed "ELECTRO" word and the background. Thus doing, the approximate orders are a opt x = 0.505 and a opt y = -0.785. The image of Fig. [START_REF]Diffraction from a wavelet point of view[END_REF] shows that two previous orders allow all parts of the image of the object to be reconstructed. Note that the word is reversed along the y-axis. This is due to the Gouy phase shift of a Gaussian beam that propagates from -∞ to +∞ through the focus point. In this case, the phase shift predicted is π.

Conclusion

The effect of a transversal object shift in an elliptical, astigmatic Gaussian beam does not affect the optimal fractional orders required to reconstruct the image of a particle or any other opaque object in the field. In this publication an analytical model has been developed to prove this. In this development, an opaque particle is considered. Experimental result presented by the recording and reconstruction of the word "ELECTRO" show that our method works, even when the object is spread out over the whole image field.

A.A. Expression for the amplitude distribution A 1

The development of integral A 1 has been given in a previous paper by [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF] :

A 1 = K ξ K η exp - π λz ρ T N ρ exp i π λz ρ T M ρ (36)
where ρ T represents the vector (x y) and the factors K q with q = ξ, η in Eq.(36) are defined by

K q =   πω 2 q λz 1 + i πω 2 q λz z Rq -1   1/2 (37)
and the diagonal matrices N and M by

N =     N x 0 0 N y     , M =     M x 0 0 M y     , (38)
with

N q = π ω 2 q λz 1 + π 2 ω 4 q (λz) 2 z Rq -1 2 , M q = 1 + π 2 ω 4 q (λz) 2 z Rq -1 1 + π 2 ω 4 q (λz) 2 z Rq -1 2 . (39)

A.B. Expression for the amplitude distribution A 2

To develop the second integral of A 2 , involving the product of E(ξ, η) and T (ξ -ξ 0 , η -η 0 ), i.e.:

A 2 = exp iπ λz (x 2 + y 2 ) λz R 2 E(ξ, η)T (ξ -ξ 0 , η -η 0 ) exp iπ λz (ξ 2 + η 2 ) exp -i 2π λz (xξ + yη) dξdη, (40)
we first replace ξ by ξ + ξ 0 and η by η + η 0 and get:

A 2 = exp iπ λz (x 2 + y 2 ) λz D E(ξ + ξ 0 , η + η 0 ) exp iπ λz ((ξ + ξ 0 ) 2 + (η + η 0 ) 2 ) × exp -i 2π λz (x(ξ + ξ 0 ) + y(η + η 0 )) dξdη. (41)
The domain D must be defined as the disk with center 0 and diameter D. By considering that c z = π/(λz) and by restating A 2 in cylindrical coordinates as follows: ξ = Dσ cos(ϕ)/2 and η = Dσ sin(ϕ)/2 for the object plane, we obtain:

A 2 = D 2 4λz exp c ξ ξ 2 0 + c η η 2 0 + ic z (x -ξ 0 ) 2 + (y -η 0 ) 2 • 1 0 2π 0 exp iγσ 2 exp iδσ 2 cos(2ϕ) exp [iaσ cos ϕ + ibσ sin ϕ] σdσdϕ (42) with γ = D 2 4 c z -i D 2 8 (c ξ + c η ), δ = i D 2 8 (c η -c ξ ), a = Dc z [ξ 0 (1 -ic ξ /c z ) -x] , b = Dc z [η 0 (1 -ic η /c z ) -y] . (43)
By writing: The function T k (r, 2γ) is defined as:

a cos ϕ + b sin ϕ = r cos(ϕ -θ) (44
T k (r, γ) = ∞ p=0 β 2k 2k+2p (δ)V 2k 2k+2p (r, γ), (49)

2 )F 2 N

 22 ≤ |α x | ≤ π/2 and 0 ≤ |α y | ≤ π/2, respectively, of a 2D-function I(x, y) is defined as (with α p = apπ αx,αy [I(x, y)](x a , y a ) = R αx (x, x a ) N αy (y, y a )I(x, y) dx dy,

  (2) and (3). The diffraction patterns consist of a 512 × 512 array of 11m × 11m size pixels. Consider the diffraction pattern presented in Fig. (2) and produced by a particle of D = 150m in diameter located at z = 120mm. The optimal fractional orders obtained from Eq. (29) are a opt x = -0.564 and a opt y = -0.850. The image of the reconstructed image is shown in Fig. (4). In this representation, the squared modulus of the FRFT, i.e. |F αx,αy [I]| 2 , is taken. The shape of the particle image is not modified: the width of the Gaussian function in the Eq.(32) is greater than the diameter D of the particle (typically 679m along x-cross axis and 490m along y-cross axis). Now, the reconstruction of the particle image from the diffraction pattern illustrated by Fig. (3), is realized by a fractional Fourier transformation of optimal orders a opt x = -0.564 and a opt y = 0.664. Figure (5) illustrates the result of the reconstruction. In both cases (Fig.(4) and Fig.(5)) reconstruction by

  ) for which we have the condition a = r cos θ, b = r sin θ (45) with complex r and θ. This representation is discussed in some detail in Appendix B. By means of the following equalities in [ [17]], 9.1.41 -9.1.45:exp iδσ 2 cos (2ϕ + 2θ) = J 0 δσ 2 + 2 +∞ k=1 i k J k δσ 2 cos 2k(ϕ + θ), ) exp[ix cos θ]dθ = i n J n (x),(47)the expression of A 2 becomes:A 2 = πD 2 λz exp [Φ(ξ 0 , η 0 )] • exp ic z (x -ξ 0 ) 2 + (y -η 0 ) 2 • ∞ k=0 (-i) k ε k T k (r, γ) cos(2kθ), (48)with ε k = 1/2 if k = 0 and 1 otherwise. The parameter denoted Φ(ξ 0 , η 0 ) is equal to c ξ ξ 2 0 + c η η 2 0 .

Fig. 1 . 20 x

 120 Fig. 1. Numerical and experimental optical set-up.

Fig. 2 . 21 x

 221 Fig. 2. Diffraction pattern with , ω ξ = 7mm, ω η = 1.75mm, R ξ = ∞, R η = -50mm, D = 150m, λ = 632.8nm, z = 120mm, δ = 250mm, ξ 0 = 0.5mm and η 0 = 0.2mm.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Diffraction pattern with , ω ξ = 7mm, ω η = 1.75mm, R ξ = ∞, R η = -50mm, D = 150m, λ = 632.8nm, z = 120mm, δ = 150mm, ξ 0 = 0.5mm and η 0 = 0.2mm.

Fig. 5 .=Fig. 6 .Fig. 7 .

 567 Fig. 5. Fractional Fourier transform of the diffraction pattern with a opt x = -0.564 and

Recall here that the expression of V 2k 2k+2p (r, γ) is:

The coefficients ω ml is given explicitly in [ [START_REF] Nicolas | Application of the fractional Fourier transformation to digital holography recorded by an elliptical, astigmatic Gaussian beam[END_REF]] and [ [START_REF] Janssen | On the computation of the Nijboer-Zernike aberration integrals at arbitrary defocus[END_REF]].

B. Appendix B : Elaboration of condition (45)

With u = a + ib, v = a -ib, we should find τ and ω = exp(iθ) such that

Assume that u = 0, v = 0, and write u = r exp(iα), v = s exp(iβ) with r, s > 0 and α, β ∈ R then we see that ω = exp(iθ) given by

satisfy the relations (51). This does not work in the case that u = 0 or v = 0. Indeed, when a = 1, b = i we get from (45) and cos 2 θ + sin 2 θ = 1 that r = a + ib = 0, i.e., r = 0.

C. Appendix C : Definition of s p

To determine the value of s p , it is necessary to write the definition of the one-dimensional fractional Fourier transformation in the particular case of α = π/2:

Its discrete version is

where δ x and δ xa are the sampling periods of I(x) and its transform. The sampling periods are equal to δ x , and N is the number of samples of I(x) and its transform. The relation (54) can be written

By identification of Eqs. ( 55) and (54), one obtains:

In the case of two-dimensional function one finally has s 2 p = N p • δ 2 p with p = ξ, η.

List of Figures