S Juma

D Honoré

B Lecordier

X Paubel

A Cessou

INFUB -11th European Conference on Industrial Furnaces and Boilers, INFUB-11 On the effect of separated oxygen and carbon dioxide injections on the stabilisation of diluted oxyfuel flames

Keywords: oxyfuel combustion, diluted combustion, Flue Gas Recirculation, CCS 1

Oxyfuel combustion with exhaust gas recycle coupled with CO 2 capture and storage (CCS) is a promising way to meet low CO 2 emission standards in industrial facilities with limited economical impact. In such systems, the flame stability is very sensitive to the dilution by injection of exhaust gases, particularly in configurations where they are not premixed with the oxygen, making the design of the burners more complex but offering a larger operation flexibility and a better control of flame heat transfer. In order to study the strategies of injection of the oxygen and exhaust gases, this paper presents an experimental study of the aerodynamic mechanisms influencing the stabilization of CO 2 -diluted oxy-fuel flames, for four different configurations of 23 kW quadri-coaxial burners with separated injections for oxygen and carbon dioxide. The four burners have same axisymmetric geometry consisting in injections of methane in the center surrounded successively by a first oxygen (O 2i ) inner annular jet, the CO 2 co-flow and a second oxygen (O 2e ) outer annular jet. Dimensions of burners are chosen to keep constant CH 4 and O 2i injections and to be able to change independently CO 2 and O 2e velocities for constant thermal power, total equivalence ratio, oxygen repartition and dilution ratio. The interaction between combustion and the aerodynamic features is investigated by CH* chemiluminescence imaging and Particle Image Velocimetry (PIV). Mean tomographic images of the flame structure are obtained by Abel's inversion of averaged chemiluminescence images. PIV measurements are performed for two fields of view in order to obtain the global aerodynamic features of the turbulent oxyfuel flames and a more precise characterization in the vicinity of the burner exit. For the latter, the spatial resolution of the measurements is optimized by the development of a specific multi-step PIV processing. Low flow-velocity and high flow-velocity configurations are tested at maximum CO 2 dilution allowed by the burners. Different structures of flames are obtained with a long continuous annular shape or with local extinction for some operating conditions. The intensity of the mixing processes and the resulting stability of the flame depend largely on the shear constraints between CO 2 and O 2e jets. A better stabilization is found for low CO 2 velocity, which favors its mixing with outer oxygen annular jet prior to direct dilution of the flame. An increase of O 2e velocity further improves centrifugal entrainment of CO 2 and then reduces the radial stratification around the flame. These results obtained in a reference configuration are useful guides for the design of flexible and efficient oxy-fuel industrial burners for CCS units.

Introduction

The socioeconomic costs linked to CO 2 emissions in the atmosphere require the development of new technologies to control the environmental impact of combustion. Oxyfuel combustion produces fumes with high CO 2 content, ready for storage (CCS) or valorization (EOR) after purification at a lower cost compared to aircombustion products. It appears as a promising mean-term solution, as it is economically competitive and does not require new technological development [START_REF] Chapel | Recovery of CO 2 from Flue Gases: Commercial Trends[END_REF]- [START_REF] Andersson | Process evaluation of an 865 MW e lignite-fired O 2 /CO 2 power plant[END_REF].

The elimination of air nitrogen affects the thermal transfers when retrofitting industrial facilities originally designed for air-combustion. An exhaust gases recycle (EGR) loop must be set-up to replace the missing volume of nitrogen by exhaust gases, which allows the control of temperatures and thermal transfers as a function of the recycle rate. However, flame properties and stability are strongly influenced by the injection of exhaust gases with high CO 2 content, which inhibit the combustion [START_REF] Heil | Experimental investigation on the effect of O 2 and CO 2 on burning rates during oxyfuel combustion of methane[END_REF]- [START_REF] Oh | The effect of CO 2 addition on the flame behavior of a non-premixed oxy-methane jet in a lab-scale furnace[END_REF]. The flame stability is increasingly challenged as the recycle rate rises. Therefore, the design of new burners adapted to oxy-fuel combustion with EGR and providing stable flames at high recycle rates is a major issue for the development of this technology.

The trend is to use existing air-burners and replace the air by a mixture of oxygen and exhaust gases. This method requires the cleansing of the exhaust gases to remove the unburned compounds prior to the mixing with oxygen within the burner. This operation is costly, and could likely favor the introduction of airborne nitrogen within the recycled flue gas, which may lead to the formation of NOx within the flame.

These drawbacks led to considering a new burner design based on the separation of the oxygen from the exhaust gases within the burner. However, this method may lead to stability issues due to separate injection of CO 2 which may act as a barrier between the fuel and the oxygen through stratification. The control of mixing processes is crucial to produce a homogeneous mixture and a stable flame, which is why the design of separate injection burners should be carefully led. The objective of the present work is to experimentally characterized these mixing processes and the corresponding diluted oxyflames with multi-coaxial burner configurations. This study is performed on a laboratory-scale facility by Particle Image Velocimetry to determine velocity fields at two different spatial resolutions and chemiluminescence imaging to obtain flame topology.

Experimental setup

For this purpose, an experimental study is performed with four 23 kW coaxial burners with separated injections for O 2 and CO 2 , in a 1.20 meter-high and 40 cm diameter cylindrical lab-scale combustion facility [START_REF] Juma | Oxycombustion diluted with CO 2 : study of the strategy of injection of CO 2[END_REF]. The four quadri-coaxial (CH 4 /O 2i /CO 2 /O 2e ) burners are based on a similar design (Figure 1 The four different configurations of quadri-coaxial burners are tested on a combustion facility [START_REF] Paubel | Analyse expérimentale des oxy-flammes turbulentes non prémélangées de gaz à faible pouvoir calorifique[END_REF], consisting in a 1.20 m high and 0.40 m diameter cylindrical chamber. The thermal power is fixed at 23 kW and the global equivalence ratio φ at 0.98, guaranteeing approximately 2% O 2 within the flue gas.

The volumetric repartition of oxygen between the primary injection (O 2i ) and the secondary injection (O 2e ) is represented by parameter Ro. The dilution rate α is defined as the volumetric proportion of carbon dioxide in the "oxidant" O 2 +CO 2 . This study is held at Ro=10%, and at the maximum dilution rate allowed by the configuration -i.e. the maximum dilution rate for which a flame is stabilized.

The methane and the primary oxygen (O 2i ) diameters and flow velocities are kept constant for the four burners C1-O1, C1-O2, C2-O1 and C2-O2, respectively in the order of 60 m/s and 15 m/s. The CO 2 and O 2e flow velocities are varied through d 3 and d 4 , in order to identify aerodynamic effects favoring the stabilization of the flame for high dilution rates α. This induces a variation of the velocity ratio U O2 /U CO2 noted as r 2 . The pipes thickness is 0.5 mm to avoid major bluff-body effects. The combustion chamber wall is continuously cooled down by water circulating within a double stainlesssteel casing. The water flowrate, inlet and outlet temperatures are controlled for thermal balance. The chamber is equipped with quartz optical windows (240x120 mm²) for chemiluminescence and PIV imaging. Eight thermocouples (penetrating 20 mm inside the chamber) and 5 fluxmeters are set on the chamber wall at different heights from the burner. A pressure probe is placed in the bottom of the chamber. A Siemens analysis unit with Ultramat 23 and Oxymat 61 analyzers records the mean dry flue gas concentrations (O 2 , CO 2 , CO, NOx and CH 4 ) after probing and cooling with an M&C® hot sampling line and cooler.

Mean flame structures are obtained from CH* chemiluminescence imaging by using a CCD camera (Roper Scientific Coolsnap 12 bits -1300x1030 pixels) equipped with a 50 mm f/2.8 Nikon lens and a set of a BG12 band-pass and a GG375 high-pass filter. The axisymmetric configuration allows the reconstruction of the mean tomographic image by operating Abel's inversion on the obtained profiles [START_REF] Susset | Développement de traitements d'images pour l'étude de la stabilisation de flammes turbulentes non-prémélangées générées par des brûleurs industriels modèles[END_REF]. The inverted images are normalized with respect to the longest exposure time, and to the highest intensity obtained for all the images. The characterization of the velocity fields is achieved with a Particle Image Velocimetry (PIV) system [START_REF] Boutier | Vélocimétrie laser pour la mécanique des fluides[END_REF]. The flows are seeded with dry ZrO 2 particles with a nominal diameter of 2 µm. The laser used is a double-pulse 200 mJ-15 Hz Nd-YAG Quantel laser operating at 532 nm. For high-resolution measurements at the burner nozzle, a 50 mm cylindrical lens and a 500 mm spherical lens are combined to perform a few centimeters high and a few-hundreds micrometers thick laser sheet, with maximum energy focused on the burner axis. The bottom of the laser sheet is cut with a sharp edge to avoid reflections on the burner surface. The measurement field size is 18x18 mm².The collection system is composed of a LaVision Imager ProX 14 bits camera equipped with an Infinity K2 microscope with its standard lens combination. A 532 nm interference filter is positioned in front of it to reduce the flame emission signal. The synchronization between the laser and the camera is achieved thanks to a programmable timing unit (PTU) and LaVision Davis 7.2 software, which is used for the collection of 300 pairs of images. The collection is realized with a 2 µs-delay between the laser pulses. The PIV correlation calculations are realized with a micro-PIV algorithm, adapted to the strong velocity gradients and low seeding behind the burner lips. This algorithm is based on the averaging of correlations, and leads directly to the meanvelocity field, through processing of the raw particle images obtained after background image subtraction. Because of the preferential flow direction, a rectangular meshing of the image is preferred, and two passes are applied with successive interrogation windows of 16x64 pixels and 6x12 pixels (eq. 0.06x0.12 mm²) with 1/3 overlap and window deformation. Intermediary validation steps consist in elimination of spurious vectors with displacement higher than 20 pixels, and a median filter over 3x3 pixels boxes. Instantaneous velocity fields are then calculated through an iterative two-steps FFT algorithm, with second-order mesh deformation, and with two passes for calculation and the mean-velocity field as an initial guess for the calculation. Final validation steps consist in elimination of spurious vectors with displacement higher than 20 pixels and a normalized median filter, adapted to strong velocity gradients, applied over 3x55 pixels boxes with parameter n=2 [START_REF] Westerweel | Efficient detection of spurious vectors in particle image velocimetry data sets[END_REF].

A global view of the mixing mechanisms further from the burner nozzle is provided by velocity measurements over large fields of view of 150x150 mm² above the burner. To provide these measurements, a cylindrical lenses combination resulting in 45 mm-focal-length, and a 750 mm-focal-length spherical lens, are combined to perform a constant-height laser sheet. The 150 mm diameter of the last spherical lens of this periscope provides illumination for this large field of view. The camera is equipped with a Nikon 50 mm f/1.2 lens with an aperture of f/5.6 and a 532 nm interference filter. LaVision Davis 7.2 software is used for the collection and processing of the data. The PIV correlation calculations are realized with two passes in 32x32 pixels interrogation windows with 50% overlap and window deformation. The absence of peak-locking is checked, as well as the convergence of mean and rms values for the 300 images acquired. The post-processing includes elimination of spurious vectors by limiting the allowed vector range (|u|<80 m/s and |v|<20 m/s), a median filter, and validation with a peak ratio threshold Q>1.3. These processes are followed by vector interpolation. The velocity fields obtained are reliable for x>D, with a validation rate over 99% on average for all considered field and for all the burners.

The combination of the different fields of measurements provides reliable velocity data from x=0 to x=37D even outside the flows thanks to the recirculation of seeded combustion products naturally occurring in the combustion chamber. The turbulent structures are well-resolved in the near-field up to x=5D, but more poorlyresolved for the 150x150 mm² field of view. For the latter, reliable fluctuation information is considered to be retrieved from x=7D.

Stability diagrams

For these four burners, the stability diagram is established as a function of Ro and dilution ratio  (Figure 2 and Figure 3). The limit of stability is defined at the extinction of flame or large CO emissions ([CO] > 400 ppm in flue gas). Except for the limits of the stability diagrams, flue gas emissions for all stable flames show that the combustion is always complete in the flames with very low CO and NOx concentrations (< 10 ppm). Considering flames structures, four types of flame are identified depending on the operating conditions (Ro ; α), and are classified as type A, A', B and C. Burner C2-O2 at (Ro = 50% ;  = 63%) illustrates type A flame on Figure 4.a, exhibiting a long continuous flame attached to the burner. The white bright flame at the burner outlet is due to the local high concentration of oxygen, and the blue color of the flame is explained by the high amount of CO 2 . The distinctive feature for flame A (Figure 4-a) is the continuous decrease of CH* chemiluminescence that can be associated to a decrease of heat release from x=0 to the end of the visualization window.

Flame A' (Figure 4-b) is characterized by a noticeable separation of a primary combustion zone attached to the burner, and a secondary combustion zone, with the separation marked by a local minimum of heat release and thinner flame front. For burner C1-O1 at (Ro = 10% ;  = 30%), this local minimum is located at x=4D. Within the second combustion downstream x=4D, the heat release increases with a local maximum at x=12D before decreasing again downstream.

The flame from burner C2-O1 at (Ro = 5% ;  = 50%) is chosen to illustrate type B flame (Figure 4-c). Chemiluminescence mean image shows a primary combustion zone extending from x=0 (burner outlet) to x=10D. The absence of chemiluminescence signal from x=10D to x=28D reflects a long flame extinction, until the signal of chemiluminescence and the flame re-appear downstream x=28D. Type C flames therefore exhibit very distinctive features: a primary combustion zone and a secondary combustion zone clearly separated by a local extinction zone.

The type C flame (Figure 4-d) has a very specific shape and is characterized by a curved bulb-shaped flame. For C2-O1 burner at (Ro = 2% ;  = 25%), a primary combustion zone extends from x=0 to x=3D and expands towards the periphery, while a secondary combustion zone converges centripetally towards the longitudinal axis from x=3D to x=5D before expanding again.

All these flames exhibit a primary combustion zone, corresponding to the injection of oxygen around the methane, and a secondary combustion zone, associated to the mixing between the methane and the secondary oxygen surrounding the CO 2 injection. The distinctive features for each type are linked to the transition area between the primary flame and the secondary flame.

The goal of this study is to analyze the mixing mechanisms in order to explain the changes between the different flames types. The operating conditions for this detailed study is set at the maximum achievable dilution ratio α max for Ro = 10%. The corresponding mean velocities for the fluids for each burner at these conditions are presented in Table 2.

 (%) V CH4 (m/s) V O2i (m/s) V CO2 (m/s) V O2e (m/s) = C1-O1 10

Identification of mean flame structures

The objective is to identify the burner configuration which provides the more stable flame structure, and to identify the mixing mechanisms that pilot that structure. The mean flame structures are identified through CH* chemiluminescence for the four burners at (Ro=10% ; α max ). Figure 5 Burner C1-O1 is referred to as the "base case", with high CO 2 velocity and low O 2e velocity (Table 2). The flame obtained for these operating conditions is type B (Figure 5-a) with primary flame extending from x=0 to x=5D, and secondary flame base starting in the region of x=26D. The zone between x=5D and x=26D is a zone of local extinction in average. When O 2e velocity is higher, for burner C1-O2 (Figure 5-b), at same dilution rate α and same CO 2 velocity, (Table 2), the flame is type A' with local decrease of heat release between x=5D and x=15D. When CO 2 velocity is lower compared to the base case, for burner C2-O1, the flame is a type A' (Figure 5-c), with a minimum of heat release and thinner mean flame localized in the region between x=15D and x=20D. When CO 2 velocity is lower and O 2e velocity is higher compared to the base case, C2-O2 configuration produces a type A flame, for which heat release intensity decreases continuously from x=0, due to dilution of the flame by CO 2 and combustion products. The global stability of the flame appears even better than for the previous configurations. From these observations of mean flame structures, it seems the stability of the flame improves from base case a (C1-O1) to case d (C2-O2). As this improvement can be expected to be correlated to the mixing mechanisms, velocity fields measurements are performed by PIV (Figure 6) from the burner outlet. The velocity fields, superimposed with chemiluminescence contours of the flame, enable to locate the flame within the flow over a total length close to 40D. The primary flame is attached to the burner, between the methane and primary oxygen injections, and the secondary flame develops at the mixing interface between the central methane jet and the peripheral CO 2 /O 2e mixture. The distinctive morphology for each flame is piloted by the mixing processes in the region of the nozzle exit. The spatial resolution of PIV measurements in such large field of view is not sufficient for a detailed description of the aerodynamics at the vicinity of the burner, as evidenced on Figure 6, and more particularly for C1-O2 and C2-O2 flame, where the O 2e exit velocity which is greater than 60 m/s, cannot be measured because of the small size of the injectors compared to the camera magnification ratio. For this purpose, a focus is made on the nozzle exit with a dedicated macro camera lens and a specific PIV algorithm. The resulting near-burner velocity fields, providing high spatial resolution, are studied comparatively for each burner (Figure 7 and Figure 8).

Stabilization mechanisms for high CO 2 and low O 2e injections velocities (C1-O1 case)

For the base case C1-O1, the flame is type B with a large extinction between 5D and 34D. In the near-field of the injector, the flame is clearly located behind the methane injection lip (Figure 7-a) at the interface between methane and O 2i [START_REF] Juniper | Edge diffusion flame stabilization behind a step over liquid reactant[END_REF], [START_REF] Candel | Structure and dynamics of cryogenic flames at supercritical pressure[END_REF]. The development of the jets, characterized by the velocity fields, puts in evidence that the flow is mainly vertical and confined over the first diameters: the O 2e flows with low expansion towards the periphery. The coaxial configuration promotes radial stratification for the jets, which have to mix through radial transfer of momentum in order to stabilize the flame. This can be seen on the radial velocity fields presented Figure 8-a. The O 2i is confined at the interface between CO 2 and methane, and entrained by the latter, as evidenced by the negative radial velocities, forming a strong primary flame. From x=4D, the mixture formed by O 2i , CO 2 and O 2e converges towards the center, which is evidenced by the centripetal velocities at this point. The CO 2 enters rapidly the flame causing the local extinction typical of this B flame. External oxygen is entrained by the faster CO 2 flow, and meets the methane jet downstream causing reactivation of the combustion.

Effect of external oxygen velocity on flame stabilization (C1-O2 case)

For C1-O2 operating conditions, the O 2e is injected at higher velocity than for C1-O1 ones (Figure 7-b). The O 2e jet flows towards the methane and transfers momentum towards the periphery, entraining burnt gases (Figure 8-b). Carbon dioxide is entrained by the O 2e flow towards the methane due to a velocity ratio r 2 higher than 1, which explains the weakening of the combustion zone, due to the interaction of the flame, with CO 2 and fresh gases, typical of type A' flame. As the O 2e converges faster towards the center, the mixing with methane occurs upstream compared to base case. The higher radial velocities (Figure 8-b) due to greater velocity gradients and jet entrainment explain the more efficient mixing between O 2 and CH 4 .

The comparison between radial velocities fluctuations (Figure 9) highlights the presence of eddies which propagate from CO 2 /O 2e interface and O 2e /burnt gases interface towards the reaction zone. It shows a better mixing for C1-O2 configuration due to the high O 2e velocity, which favors the momentum transfer and the mixing of O 2 and CH 4 , thus enabling to obtain a type A' flame, with attachment of the secondary flame to the primary flame (case C1-O2) instead of a type B flame for C1-O1 operating conditions.

C1-O1 (Ro

Effect of carbon dioxide velocity on flame stabilization (C2-O1 case)

For C2-O1 flame, the CO 2 is injected at lower velocity than for C1-O1, which is induced by a wider injection surface, and for a close dilution ratio α=50%. The O 2e injection velocity is the same than the C1-O1 base case. The flame expands more towards the periphery than the base case, as shown by chemiluminescence contours (Figure 7-c). The confinement of the primary flame is broken as evidenced by the diverging streamlines close to the flame (Figure 8-c), which favors the expansion of the methane towards the periphery and mixing with the oxygen. This causes remoteness of CO 2 far from the reaction zone. Thus, with lower local dilution, a higher heat release is observed for the primary flame for this configuration compared to previous flames: the chemiluminescence contour for 10,000 a.u. extends up to x=4D, against approximately x=2.5D for burners C1-O1 and C1-O2.

The CO 2 stream is pinched between O 2i and O 2e (Figure 7-c) due to the local depression linked to the lower velocity of CO 2 , favoring the mixing of the reactants, and thereby inducing the attachment of the secondary flame to the primary flame. This configuration emphasizes the stabilizing effects of the decrease of CO 2 velocity, which proves to be more efficient for flame stability in our configuration than increasing the velocity of O 2e at high CO 2 velocity.

4.4.

Optimization of injections strategy for flame stabilization enhancement (C2-O2 case)

The configuration C2-O2 (Figure 7-d) combines low CO 2 -velocity and high O 2e velocity, leading to a large velocity ratio (r 2 = 7.7). The attached flame expands again towards the periphery and is not confined at methane-O 2i interface, as shown by the CO 2 streamlines which diverge from the center (Figure 8-d). The primary oxygen jet expands towards the periphery with higher radial velocity (up to 6 m/s) than for burner C2-O1 (up to 4 m/s) (Figure 8-c andd). This effect, combined with the expansion due to the flame drives CO 2 away from the flame, favoring better stability of the flame. The secondary oxygen flow converges towards inner CO 2 jet with higher centripetal radial velocity (down to -6 m/s) than for burner C2-O1 (down to -2 m/s). This results in a stronger pinching of the CO 2 jet between higher-velocity O 2i and O 2e flows, which is characteristic of the local pressure decrease caused by high r 2 [START_REF] Rehab | Flow regimes of large-velocity-ratio coaxial jets[END_REF]. This fast merging between the two oxygen flows shortens the CO 2 potential core length to approximately 3D, and then enhances the mixing between carbon dioxide and oxygen.

The streamlines from O 2i , CO 2 and O 2e outlets tighten as the flows develop, emphasizing the increase of turbulent mixing mechanisms, which is consistent with the high velocity gradients generating high shear stress. Ultimately, the flame is type A with secondary flame directly attached to the primary flame, with no decrease of the combustion intensity thanks to the injection of O 2e in the combustion zone, balancing the CO 2 dilution in the flame.

High r 2 velocity ratio linked to C2-O1 and C2-O2 configurations promote high shear strain, which results in an efficient mixing of O 2e with CO 2 as observed on planar turbulent kinetic energy fields (Figure 10). The increase of r 2 induces shortenings of CO 2 and O 2e potential cores: their respective lengths decrease from 5D and 1.5D for burner C2-O1, to 2.5D and 0.5D for burner C2-O2. Moreover, for the latter case, the high kinetic energy, due to the high velocity gradients, enhances downstream mixing between the central CH These results put in evidence that the stability of the flame improves when O 2e velocity is increased and CO 2 velocity is lowered, due to mixing mechanisms improvement through velocity gradients. The importance of parameter r 2 in the stabilization process is underlined. In our studied configurations, CO 2 acts as a barrier between CH 4 and O 2e which is the main source of oxidizer (Ro = 10%). The flame stability improves when radial stratification is lowered, due to the convergence of O 2e outer flow towards CH 4 central jet.

The influence of r 2 velocity ratio is also pointed out through the kinetic turbulent energy growth at the center of methane jet (Figure 11) for both C2-O1 and C2-O2 configurations, which are characteristic of an enhanced mixing. The methane potential core length is even shorter for burner C2-O2 (~11D) than for burner C2-O1 (~13D) which is characteristic of the faster convection by O 2e jet of the turbulent structures generated in the main shear layers and provides better stabilization for the flame. For both C2 configurations, linked to high r 2 due to the geometry of the injectors, the maximum attainable dilution rate is higher than for C1 configurations (Figure 12), where only lower r 2 are possible.

Conclusions

Four quadri-coaxial configurations of CH 4 /O 2 /CO 2 /O 2 burners, with separated injections for CO 2 and O 2 have been investigated. The separation of carbon dioxide and oxygen injections in oxy-burners with recycling of exhaust gases offers a unique flexibility in terms of fuel nature and CO 2 dilution but may imply a stratification of the flows which lead to flame instability in some configurations. The four different configurations give the opportunity to study the parameters controlling flame stability through modification of the mixing mechanisms of the reactants within the combustion chamber.

Due to the low tip thickness and large velocity gradients between the jets, the development of wake instabilities is moderate, and the mixing processes are mainly controlled by the development of turbulent structures at the interface between the jets. The intensity of the mixing processes and resulting stability of the flame are found to depend largely on shear constraints initiated by the velocity gradients between the co-annular flows of CO 2 and O 2e .

Two parameters are identified in order to control the mixing, and to balance the radial stratification of the flows and the destabilization of the reaction zone: the CO 2 velocity, and its comparison with O 2e velocity (r 2 ratio). In particular, it is shown that delaying the entrance of CO 2 within the flame (e.g. by lowering CO 2 velocity) and favoring the O 2e mixing with CO 2 and convergence of O 2e jet towards central fuel jet through high r 2 ratio and low CO 2 pressure tend to stabilize the flame and suppress the local extinction zones observed at high CO 2 velocity. The homogeneity of the mixing between O 2 and CO 2 is a necessary condition for flame stability, and to avoid stratification at the flame base when CO 2 dilution is increasing. These observations are useful guides for the design of flexible separated O 2 /exhaust gases injection burners at industrial scale. Indeed, even if the present study is led at 23 kW with simple pipe-in-pipe configurations and voluntarily few stabilizing artefacts, it puts in evidence that separated injections of CO 2 and oxygen are possible, but require more complex configurations for industrial applications in order to reach higher dilution rates. Industrial burners are generally designed for power ranges in the order of 100 kW to several MW, and involve systems of swirls, thick walls that maximize bluff body effect for stabilizing the flame and multiple oriented injections. Such configurations would favor the mixing of the reactants and reduce the radial stratification between CO 2 and O 2 flows. According to this study, a burner configuration providing thick walls or swirls, low CO 2 velocity with injectors directed towards the periphery and high velocity O 2 directed towards a central natural gas jet, would be adapted to withstand higher dilution rates at higher power. Tests in a pilot facility of such designed burners could validate these guidelines and may be an intermediate step towards scale-up at industrial level.

  ). They are composed of a 4 mm-diameter (D) central injection of methane surrounded by a first coaxial injection of oxygen (O 2i ). The next coaxial injection of carbon dioxide is surrounded by the last annular injection of secondary oxygen (O 2e ). Carbon dioxide is injected between CH 4 and O 2e representing 90% of the total oxygen amount in order to guarantee a thorough dilution of CO 2 within the flame. The separation of CO 2 and O 2e injections promotes radial stratification of reactants dilution, which may destabilize the flame. The presence of internal oxygen injection O 2i is chosen to generate a central low intensity oxyflame attached to the burner. The variation of flow velocities of CO 2 and O 2e through changes of injection diameters d 3 and d 4 allows the study of the effect of mixing processes on the flame stabilisation. Two diameters d 3 and two diameters d 4 are considered, leading to four different burners providing independently either high or low CO 2 velocity, and either high or low O 2e velocity.

Figure 1 -

 1 Figure 1 -General design of the burner.

Figure 2 -Figure 3 -

 23 Figure 2 -Stability diagram for burners C1-O1 and C1-O2

Figure 4 -

 4 Figure 4 -Photographs and Abel's inversions for mean chemiluminescence image for A (image a), type A' (image b), type B (image c) and type C (image d) flames. The burner longitudinal axis corresponds to r=0. The chemiluminescence signal is given in arbitrary units.

Figure 5 -

 5 Figure 5 -Mean tomographic CH*chemiluminescence images obtained after Abel's inversion or the four burners at their maximal dilution rate for Ro=10%.

Figure 6 -

 6 Figure 6 -Mean fields of velocity magnitude obtained through PIV. Mean contours of the flame obtained through Abel's inversion of mean CH* chemiluminescence image (at 300 and 5000 a.u.) are superimposed.

Figure 7 -Figure 8 -

 78 Figure 7 -Mean fields of velocity magnitude superimposed with mean velocity vectors -1/35 vector along x/D-axis and 1/5 vector along r/D-axis -obtained through PIV in the near-field of the injector. Mean contours of the flame obtained through Abel's inversion of mean CH* chemiluminescence image (intensity in arbitrary units) are superimposed.

Figure 9 -

 9 Figure 9 -Radial fluctuations v' fields superimposed to mean flame contours obtained through Abel's inversion of CH* chemiluminescence image (500, 5000 and 10000 a.u.) for burners C1-O1 and C1-O2.

Figure 10 -

 10 Figure 10 -Turbulent kinetic energy fields for C2-O1 and C2-O2 burners superimposed with mean flame contours.

Figure 11 -

 11 Figure 11 -Longitudinal evolutions of turbulent kinetic energy along the burner centerline for all configurations.

Figure 12 -

 12 Figure 12 -Stability diagram for the four burners as a function of r 2 and α. The lines represent the critical operating conditions (α max ; r 2max ) where limits of flame stability are experimentally reached.

Table 1 -Outlets dimensions of the four burners Burner d 1 (mm) d 2 (mm) d 3 (mm) d 4 (mm) C1-O1 4 6 9.46 11.8 C1-O2 14.0 C2-O1 18 19.8 C2-O2 21.2

 1

Table 2 -

 2 Mean fluid velocities at burner outlets at the study defined operating points.

	Burner	Ro
		(%)


Isabelle Sivignon
email: sivignon@lis.inpg.fr

Florent Dupont
email: fdupont@liris.cnrs.fr

Jean-Marc Chassery
email: chassery@lis.inpg.fr

Discrete Surfaces Segmentation into Discrete Planes

This paper is composed of two parts. In the first one, we present an analysis of existing discrete surface segmentation algorithms. We show that two classes of algorithms can actually be defined according to discrete surface and plane definitions. In the second part, we prove the link between the two classes presented. To do so, we propose a new labelling of the surface elements which leads to a segmentation algorithm of the first class that can be easily transformed into a segmentation algorithm of the second class.

Introduction

The segmentation of a discrete object surface into pieces of discrete planes is the first step of more global processes like polyhedrization of the surface [START_REF] Françon | Polyhedrization of the boundary of a voxel object[END_REF] or surface area estimation [START_REF] Klette | Digital planar segment based polyhedrization for surface area estimation[END_REF]. This problem can be stated as follows: given a discrete object, a pair of connectivities used for the object and the background, and a definition of the discrete surface, label all the surface elements of the object such that the elements with the same label belong to the same piece of discrete plane. Each set of surface elements with the same label is called a discrete face.

This statement raises the problem that both discrete surfaces and discrete planes have to be defined before the segmentation process. In a general way, there exist two families of discrete surfaces definition: either the voxels or the surfels are the surface elements. In the same way, two types of discrete planes are mainly used: the naive planes and the standard planes.

In the literature, those discrete surface and discrete plane definitions are all used in segmentation algorithms, but most of the time, few arguments are given to justify the use of one definition instead of another. Some authors even propose several algorithms using different planes and surfaces definitions, but did not explore the relation between the different algorithms [START_REF] Papier | Polyédrisation et visualisation d'objets discrets tridimensionnels[END_REF]. Indeed, by now, no work has been accomplished in order to study the link between, first, the discrete surface and discrete plane definitions, and second, the segmentation processes that use different kinds of discrete planes.

In Section 2, we define those discrete surfaces and point out the links between the surfaces and the different types of discrete planes. Then, we propose a short overview of existing segmentation algorithms, in relation with the surface and plane definitions given previously. Next, in Section 3, we recall a segmentation algorithm we proposed in a previous work. This algorithm uses naive planes over the surface voxels and is quite different from existing algorithms since the voxels' faces are labelled instead of the voxels themselves. We moreover show how the 6connectivity between the object voxels can be handled in this algorithm. Finally, in Section 4, we propose a transformation from this segmentation algorithm with naive planes into a segmentation algorithm with standard planes. This transformation shows that the two segmentation processes are in fact equivalent. Some image results are given before a few words of conclusion.

Different Discrete Planes for Different Kind of Surfaces

Two Types of Surfaces

The two existing approaches differ by the surface elements considered: the surface elements are either voxels (dimension 3) or voxels' faces (dimension 2).

In his PhD thesis, Perroton [START_REF] Perroton | Segmentation parallèle d'images volumiques[END_REF] calls those approaches homogeneous and heterogeneous respectively, since in the first one the surface elements are space elements, whereas in the second one the dimension of the surface elements is lower than the dimension of the discrete space.

Homogeneous Approach. In 1982, a definition of discrete surface in dimension n is proposed by Udupa et al. [START_REF] Udupa | Boundary detection in multidimension[END_REF], and a surface tracking algorithm is described. Nevertheless, no Jordan's theorem proving that this surface separates the object voxels from the background voxels has been shown with this definition. The year before, Morgenthaler et al. [START_REF] Morgenthaler | Surfaces in three-dimensional digital images[END_REF] gave a local definition of a discrete surface, and a global Jordan's theorem based on this definition is shown. Malgouyres also proposed [START_REF] Malgouyres | A new definition of surfaces of Z 3 . A new 3D Jordan theorem[END_REF] a local definition for 26-connected discrete objects with a 6connected background. A Jordan's theorem is also shown and the link between Morgenthaler's surfaces is stated [START_REF] Bertrand | Some topological properties of discrete surfaces[END_REF][START_REF] Malgouyres | A new definition of surfaces of Z 3 . A new 3D Jordan theorem[END_REF]. Another surface definition involving ndimensional surface elements has been recently proposed in [START_REF] Brimkov | Curves, surfaces and good pairs[END_REF], where separating pairs of connectivity are studied according to this definition.

Nevertheless, this approach causes some problems [START_REF] Perroton | Segmentation parallèle d'images volumiques[END_REF]: first, the surface of an object and the surface of the background are different, and second, for some objects containing a cavity, the two surfaces expected (two surfaces should separate the object from the two background components: the cavity, and the infinite background) may have voxels in common.

Heterogeneous Approach. This approach is based on the abstract cellular complex introduced by Kovalevsky [START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF] for discrete images topology. In dimension 3, the following cells are defined (see Figure 1 In 1992, Herman [START_REF] Herman | Discrete multidimenional jordan surfaces[END_REF] generalizes this definition to any dimension. This definition has moreover been used to write surface tracking algorithms: first, Artzy et al. [START_REF] Artzy | The theory, design, implementation and evaluation of a three-dimensional surface detection algorithm[END_REF] proposed a tracking algorithm for 18-connected objects; then, Gordon and Udupa [START_REF] Gordon | Fast surface tracking in three-dimensional binary images[END_REF] extend this algorithm for 6-connected objects improving it; in 1994, Udupa [START_REF] Udupa | Multidimensional digital boundaries[END_REF] adds a connectivity definition over surfels, and proposes a very efficient tracking algorithm; in [START_REF] Perroton | A new 26-connected objects surface tracking algorithm and its related pram version[END_REF], Perroton describes a parallel surface tracking algorithm for 26-connected objects.

In the following, we will see how this definition may be used in two different ways in the context of surface segmentation.

2.2

Which Planes for Which Surfaces?

In a segmentation process, the discrete points of the object surface are labelled according to the discrete face they belong to. In the literature, two main approaches have been developed: the discrete points are either the voxels or the pointels of the object. The definition of a discrete surface given in Definition 2 states that the surface is composed of all the surfels s = {v, v} such that v ∈ O and v / ∈ O. From this definition, we derive that the segmentation is done over the pointels of the object surface. Nevertheless, another point of view is given by the definition of the inner boundary of an object: Definition 3. Let O be a 3D discrete object. The inner boundary of O is defined as:

II(O) = {v | {v, v} is a surfel of O's surface}
The voxels of II(O) are the underlying voxels of the surface surfels. When a segmentation is done over the object voxels, the voxels of II(O) are usually labelled.

Let us recall the analytical definition of a discrete plane [START_REF] Réveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Andrès | Discrete analytical hyperplanes[END_REF]:

Definition 4.
A discrete plane of parameters (a, b, c, µ) (a, b, c and µ belong to Z and a, b and c are relatively prime) and thickness ω is the set of voxels (x, y, z) verifying the following inequalities:

0 ≤ ax + by + cz + µ < ω
The thickness parameter ω defines many plane classes, but two types of discrete planes are mainly used in the literature. When ω = max(|a|, |b|, |c|), the discrete plane is called naive, and when ω = |a| + |b| + |c|, the plane is called standard. Naive planes are the thinnest planes without 6-connected holes, and standard planes are the thinnest planes without tunnels [START_REF] Andrès | Discrete analytical hyperplanes[END_REF]. For further information about digital planarity, see for instance [START_REF] Brimkov | Digital Planarity -A Review[END_REF].

Table 1 presents the links between discrete planes and discrete surfaces. On one hand, if the discrete points are the voxels of the inner boundary, then naive planes are used for the segmentation. Indeed, the voxels that are 18-connected but not 6-connected with the background are needed for a segmentation with standard planes, but do not belong to the inner boundary of the object. On the first line of the table, a naive plane recognition has been done over the inner boundary of the object: the voxels of II(O) are labelled.

On the other hand, the pointels belonging to the object's surface are 6connected, which suits well to a segmentation with standard planes. On the second line of the table, a standard plane recognition has been done over the pointels of the surfel boundary of the object: a surfel is colored if and only if its four adjacent pointels belong to the discrete face.

Short Overview of Existing Algorithms

The first segmentation algorithm was proposed by Debled-Rennesson [START_REF] Debled-Rennesson | Etude et reconnaissance des droites et plans discrets[END_REF]. She applies the naive plane recognition algorithm she proposed on the inner boundary of a discrete object. However, her algorithm works only for objects having known symmetries. With a totally different idea, Vittone proposes in [START_REF] Vittone | Caractérisation et reconnaissance de droites et de plans en géométrie discrète[END_REF] a segmentation algorithm based on the fact that any naive plane can be decomposed into basic structures [START_REF] Vittone | n -m)-cubes and farey nets for naive planes understanding[END_REF]. But only a few preliminary results are given and this method has not been further extended by now. The following year, Vittone [START_REF] Vittone | Recognition of digital naive planes and polyhedrization[END_REF] wrote another segmentation algorithm using the naive plane recognition algorithm she proposed. This recognition algorithm is based on a dual transformation, which we also use in our work. This will be detailed in the next section. In [START_REF] Papier | Polyédrisation et visualisation d'objets discrets tridimensionnels[END_REF], Papier proposes another segmentation algorithm based on Debled-Rennesson naive plane recognition algorithm. However, his algorithm does not give satisfactory results since many small discrete faces (with a few voxels) are recognized. This observation led him to propose a segmentation algorithm of the surfel boundary with standard planes [START_REF] Françon | Polyhedrization of the boundary of a voxel object[END_REF]. The discrete faces recognized are forced to be connected and homeomorphic to disks in order to be used in a polyhedrization process. A Fourier-Motskin elimination algorithm is used to recognize the standard planes, which makes the overall algorithm complexity high. Moreover, many small discrete faces are still recognized by this algorithm, and some tricks like discrete face size limitations are proposed to better the results. Finally, Klette and Sun proposed in [START_REF] Klette | Digital planar segment based polyhedrization for surface area estimation[END_REF] a segmentation algorithm close to the one of Papier in order to estimate the object surface area. In this case, no particular constraint over the shape of the discrete faces is required, except the connectivity, which is implicitly obtained since a breadth-first tracking of the surfels is done.

For the algorithms of the first paragraph, a segmentation of the inner boundary is done, and the voxels are labelled. The algorithms of the second paragraph use a standard plane segmentation of the object pointels, and a surface surfel is labelled by a discrete face when its four adjacent pointels belong to the face.

Surfel Labelling Using Naive Planes

In this section, we propose a segmentation algorithm of the inner boundary of an object using naive planes. But instead of labelling the voxels of the inner boundary, we label only some of the boundary's surfels. Indeed, all the existing algorithms using naive planes label the voxels of the object inner boundary. Nevertheless, the inner boundary definition is derived from the surfel surface def-inition. Thus the surface elements are the surfels, and it seems more appropriate to label the surfels. In Section 4, we show how this segmentation of the inner boundary can be easily transformed into a segmentation into standard planes, which link the two segmentation approaches existing by now.

Since the final segmentation uses standard planes, we consider a 6connectivity relationship for the voxels of the object. Indeed, 18-and 26-connectivities do not enable the description of the object surfel surface as a 2dimensional combinatorial manifold, which is inconsistent with the use of standard planes [START_REF] Françon | Sur la topologie d'un plan arithmétique[END_REF].

Algorithm Description

Part of this algorithm is presented in [START_REF] Coeurjolly | Reversible discrete volume polyhedrization using marching cubes simplification[END_REF] and we first recall here the main points of this method.

The general idea of this algorithm is to get the following labelling of the surface surfels: if a surfel s = {v, v} of the surface (v belongs to the object and v does not) is labelled with the discrete face f , then there exists a Euclidean plane p crossing the segment [vv[. In other words, the discretization of the plane p with the OBQ discretization [START_REF] Groen | Freeman-coe probabilities of object boundary quantized contours[END_REF] scheme contains the voxels of f .

In the following, we denote by l one direction of the set {(1, 0, 0), (-1, 0, 0), (0, 1, 0), (0, -1, 0), (0, 0, 1), (0, 0, -1)}. Then for each voxel v of the inner boundary of an object, there exist at least one direction l such that the voxel v + l does not belong to the object. The naive plane recognition used is a directional recognition algorithm: Definition 5. Let V be a set of voxels. The directional recognition in direction l computes the set of Euclidean planes crossing all the segments

[vv[ such that v ∈ V and v = v + l.
This recognition is achieved using Vittone's naive plane recognition algorithm [START_REF] Vittone | Recognition of digital naive planes and polyhedrization[END_REF]. This algorithm is based on a dual transformation. Consider for instance the direction l = (0, 0, 1). Then the Euclidean planes of parameters (α 0 , β 0 , γ 0 ) which cross the segment [v, v+l[ where v = (x, y, z) are those fulfilling the double inequality 0 ≤ α 0 x + β 0 y + z + γ 0 < 1. Then, in the parameter space (α, β, γ), the set of Euclidean planes computed by the directional algorithm on a set V is a convex polyhedron defined as the intersection of all the double inequalities related to the voxels of V .

But this recognition algorithm enables to label only the surfels s = {v, v } where v = v+l if l is the direction of recognition (see Figure 2(a)). A second step of the algorithm is to label the inner surfels of the plane, which are the surfels sharing two edges with surfels already labelled by the directional algorithm (see Figure 2(b) and (c)). It is easy to prove that the Euclidean planes that cut the segments related to the two surfels defining an inner surfel also cut the segment related to this inner surfel.

After this short presentation of the recognition process used, let us see how this algorithm is applied on a surface. Since the quality of the segmentation is not the main point of this paper, and for the sake of clarity of the rest of this work, we do not present the algorithm described in [START_REF] Coeurjolly | Reversible discrete volume polyhedrization using marching cubes simplification[END_REF], but quickly explain a simple method. A seed voxel v 0 is chosen for the discrete face to recognize in direction l. Neighbour voxels are added one by one if they belong to the same naive plane as the voxels of f and if they fulfill the condition presented in the next paragraph. This process starts over with a new seed voxel until all the surfels are labelled.

6-Connected Objects

Since 6-connectivity is considered for the discrete object, a discrete face recognized on the object surface should not be split over two different 6-connected components. Since naive planes are 18-connected, we have to define the cases where two voxels are 18-connected in a 6-connected object. Examples of segmentation results using 6-connectivity are proposed in Figure 5. The recognition in direction (0, 0, 1) only has been computed on those examples. In (a), one unique plane is recognized; in (b), two planes are recognized since the object is composed of two 6-connected components; in (c), three planes are recognized, one for each 6-connected component.

l v v (a) v v (b) v v (c) v v (d) v v (e)

Results

Some results of the algorithm proposed in [START_REF] Coeurjolly | Reversible discrete volume polyhedrization using marching cubes simplification[END_REF] with the 6-connectivity constraint are presented in Figure 6. The results obtained with the simple algorithm presented in this paper simply contain more small faces. The first image (a) is an ellipsoid of parameters [START_REF] Vittone | Caractérisation et reconnaissance de droites et de plans en géométrie discrète[END_REF][START_REF] Réveillès | Géométrie discrète, calcul en nombres entiers et algorithmique[END_REF][START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF]

Segmentation Using Standard Planes

In this section, we show that the segmentation using naive planes over the inner boundary voxels presented in the previous section can be very easily transformed into a segmentation of the surface pointels into standard planes.

Transformation

Let us consider an infinite 3D discrete object defined as a half-space bounded by a naive plane of parameters (a, b, c, µ). Let us moreover suppose that the normal vector of this plane is directed towards the discrete object itself. Then we study the surfel surface of this object (see Figure 7(a)).

We use here the definition of standard plane proposed in [START_REF] Andrès | Defining discrete objects for polygonlization : the standard model[END_REF] and directly derived from the standard digitization scheme: Definition 8. Let P be a Euclidean plane defined by the equation ax + by + cz + µ = 0, with a,b,c and µ integers. Then the standard digitization of P is the set of voxels (x, y, z) such that: Note that the possible values of ax + by + cz + µ are centered around zero. The following theorem makes the link between the naive plane parameters and the pointels of the surface.

- |a| + |b| + |c| 2 ≤ ax + by + cz + µ < |a| + |b| + |c| 2
Theorem 1. Let P be a naive plane of parameters (a, b, c, µ) and S be the surfel surface described as above. Then the pointels belonging to S belong to the standard plane of parameters (a, b, c, µ ) in the discrete grid translated by t = ( 12 , 1 2 , 1 2 ), with:

µ = τ - |a| + |b| + |c| 2 -v.(a, b, c)
where

• τ = µ if (a > 0) or (a = 0 and b > 0) or (a = 0 and b = 0 and c > 0) µ + 1 otherwise • v = -1 2 (sgn(a), sgn(b), sgn(c)) + ( 1 2 , 1 2 , 1 2 ) (sgn(x) is equal to 1 if x > 0, -1 if x < 0, and 1 or -1 if x = 0).
Proof. This proof works by construction. Figure 8 illustrates the different steps of this construction from a naive 2D segment to a standard 2D segment. This example will be detailed in the next paragraph.

Since we transform a naive plane into a standard plane, the first thing is to thicken the naive plane P . This is equivalent to adding some voxels that fulfill larger inequalities (the thickness determines the bounds of the inequalities in a discrete plane definition). This is done in three steps:

• Since the inequalities defining a standard plane depend on the signs of the parameters, a preprocessing step is required when (a < 0) or (a = 0 and b < 0) or (a = 0 and b = 0 and c < 0). We use the fact that the two doubleinequalities 0 ≤ ax + by + cz + µ < ω and 0 < ax + by + cz + µ + 1 ≤ ω are equivalent when integer values are considered to do the parameters translation (a, b, c, µ) → (a, b, c, µ + 1) in the case of a non standard orientation. This operation simply swap the large and strict inequalities in the discrete plane definition. In the following, we suppose that the parameters are in a standard orientation since the two cases are similar. ) must be done.

All the voxels of the discrete object that are 18-connected or 26-connected with the background are added by those transformations, and the translation (a, b, c, µ) → (a, b, c, µ(+1) -|a|+|b|+|c|

2

) has been done. The next part of the transformation consists in moving the discrete grid such that the discrete points are now the pointels.

A one-to-one and onto transformation is defined between the voxels of the plane P after thickening and the pointels of the surface: the pointel p = (x - For each case, the pointels that belong to the surface are underlined. We note that for each voxel (x, y, z) of the thickened plane P , the pointel (x - 12 , y -pointels of each voxel belong to O's surface, and then four transformations are possible. This transformation defines the first part of the translation vector v of the theorem, and a simple translation of the grid of the vector t = (1 2 , 1 2 , 1 2 ) leads to the final result.

The following Corollary shows how to transform a segmentation with naive planes into a segmentation with standard planes: Corollary 1. Let f be a discrete face recognized on the inner boundary of a discrete object O with the algorithm presented in Section 3. If the OBQ discretization of the Euclidean plane of parameters (a, b, c, µ) contains the voxels of f , then the standard discretization of the plane (a, b, c, µ ) (where µ is defined as in Theorem 1) contains all the pointels belonging to the labelled surfels of f , in a discrete grid translated by t = ( 12 , 1 2 , 1 2 ).

The proof of this corollary is straightforward from Theorem 1.

Example

Figure 8 is an example of the transformation steps of Theorem 1 for a discrete segment. The first column represents the discrete segment, the second one represents the set of Euclidean solution lines in the parameter space, and the third one the parameters of the Euclidean lines related to the vertices of the convex polygon of the parameter space. On the first line, a segment of the naive line of parameters (2, -3, 3) is represented. On the second line, the naive segment has been thickened in order to get a standard segment. The parameters translation corresponding to this thickening is done. The third line represents the standard pointel segment derived from the standard voxel segment of the second line. A second translation of the Euclidean planes parameters is accomplished.

Results

Since the transformation presented above only concerns the parameters of the Euclidean planes recognized, visual results are hard to show. Nevertheless, Figure 9 presents a result of this transformation for a catenoid 1 . In (a), the surfel labelling using a naive plane segmentation over the inner boundary is represented: all the voxels containing surfels of the same color belong to the same discrete naive plane. In (b), the pointel labelling is represented after the transformation defined in Theorem 1: all the pointels with the same color belong to the same standard discrete plane. The same colors are used for the two images, so that the comparison is easier. Note that since each pointel belongs to several discrete faces, its color is the one of the first recognized face it belongs to. 2 ) • B = (2, -3, 9 2 ) • C = (1, -2, 7 2 ) • D = (3, -4, 9 2 ) Fig. 8. All the steps from a naive voxel segment to a standard pointel segment

Conclusion

In this paper, we showed that two classes of discrete surface segmentation algorithms can actually be considered: the first one consists in looking for naive planes over the surface inner boundary, and the second one uses standard planes over the surface surfels. We recalled the different existing segmentation algorithms highlighting the fact that they all belong to one of those two classes.

The second part of this work was to show that those two types of segmentation are equivalent, and we expressed the transformation that changes from one to the other. To do this study, we proposed a segmentation algorithm using naive planes and generating a surfel labelling for 6-connected discrete objects.

Since the two segmentations are equivalent, an interesting future work would be to propose a complete study on the quality of segmentation algorithms, pointing out the influence of the choice of the first voxels of the discrete faces for instance. Such a study is very useful since the quality of the surface segmentation induces part of the quality of the polyhedrization or surface area estimation based on it.

  (a)):voxel: cell of dimension 3 -surfel: cell of dimension 2 -linel: cell of dimension 1 -pointel: cell of dimension 0.

Fig. 1 .Definition 2 .

 12 Fig. 1. (a) Cellular decomposition of a voxel. (b) A surfel defined by two 6-connected voxels

Fig. 2 .

 2 Fig. 2. Labelling of the surfels of a piece of plane of parameters (1, 3, -5, 0) using the directional algorithm in direction (0, 0, 1) (from [24]): (a) just after the directional algorithm; (b) illustration of an inner surfel; (c) after the labelling of inner surfels

Definition 6 .Definition 7 .

 67 Let v and v be two voxels of a 6-connected object O. Then v and v are 18-neighbours in O if and only if one of the following two conditions is fulfilled (see Figure3): -v and v are 6-connected (they share a surfel) -v and v share a linel, and one out of the other two voxels sharing this linel belongs to O. Now consider a directional recognition in direction l. Two voxels v and v are l-neighbours in a 6-connected object O if only if they fulfill the following three conditions: -v and v belong to the inner boundary of O and the voxels v + l and v + l do not belong to O the projection of v and v along direction l are 4-connected v and v are 18-neighbours in O (see Definition 6).

Figure 4

 4 Figure 4 presents some examples of voxels that are or are not l-neighbours. Finally, during the segmentation process, a new voxel can be added to a discrete face f recognized in direction l if and only if it is l-neighbour with a voxel of f .

Fig. 3 .

 3 Fig. 3. The voxels v and v are 18-connected in a 6-connected object if and only if they are in configuration (a) or if v1 or v2 belongs to the object in configuration (b)

Fig. 4 .

 4 Fig.4. The two voxels v and v are l-neighbours (see l on the left) only in configuration (a) : in (b), the first condition of Definition 7 is not fulfilled since the surfel {v, v + l} does not belong to the surface; in (c), the second condition is not fulfilled since the projections of v and v along direction l are not 4-connected; in (d), the two voxels dot not fulfill the conditions of Definition 6; in (e), the two last conditions are fulfilled, but like in (b), the first one is not

  which is composed of 4190 surface surfels. 131 discrete faces are recognized on this surface, and 15 among them are composed of less than 5 surfels. The image in (b) is a microscopic scale view of part of a bone. Its surface is composed of 5298 surfels, and 201 discrete faces are recognized by the algorithm. 65 discrete faces are composed of less than 5 surfels.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Influence of the 6-connectivity on the segmentation computed

Fig. 7 .

 7 Fig. 7. (a) Discrete half-space bounded by the naive plane of parameters (2, 3, -5). The surfels of the surface are in grey. (b) Illustration of the bijective transformation between a voxel V and a pointel p for given plane parameters

1 2

 1 sgn(a), y -1 2 sgn(b), z -1 2 sgn(c)) is associated to the voxel V = (x, y, z). This transformation is of course bijective but we must prove that if V belongs to P , then p belongs to the surface of O. Let us consider a plane of parameters (a, b, c) with a > 0, b > 0 and c < 0. The other cases are symmetrical. Then the transformation maps the point V to the pointel p = (x -1 2 , y -1 2 , z + 1 2 ) (see Figure7(b)). Only three types of voxels may appear in the thickened discrete plane P : -: 6-connected with the background -: 18-connected with the background -: 26-connected with the background

•

  Discrete segmentSolution Euclidean planes in the parameter space (α, β)Parameters (a, b, µ) of the solution Euclidean lines ax+ by+ µ = 0 A = (2, -3, 3) • B = (2, -3, 4) • C = (1, -2, 3) • D = (3, -4,

Fig. 9 .

 9 Fig. 9. Segmentation of a discrete catenoid surface: (a) surfel labelling; (b) pointel labelling

Table 1 .

 1 Links between the discrete surfaces and the discrete planes used for the segmentation process

	Surface definition	Discrete points	Discrete planes	Example
	Inner boundary	Voxels	Naive planes
	Set of surfels	Pointels Standard planes

•

  The thickness of a standard plane is |a| + |b| + |c|. Since we supposed that the normal vector is directed towards the object, the voxels V = (x, y, z) fulfilling the condition max(|a|, |b|, |c|) ≤ ax + by + cz + µ < |a| + |b| + |c| must be added to thicken the plane inside the discrete object.• At this point, our plane is defined by the double inequality 0 ≤ ax + by + cz + µ < |a| + |b| + |c|. To be consistent with Definition 8, a last translation of the parameters of vector (0, 0, 0, -|a|+|b|+|c|

	2

, z + 1 2 ) always belongs to O's surface. For other parameters' signs a similar construction can be done, and finally, the transformation always maps a voxel of P into a pointel of the surface. If one of the parameters a, b or c is equal to zero, then two pointels of each voxel of P are always on O's surface. Then two bijective transformations are possible and correspond to the two possible values of sgn(x) when x = 0. When two out of the three parameters are equal to zero, then 4

This image has been generated using the volgen program available on http://www.cb.uu.se/tc18/code_data_set/Code/Volgen/


Ibrahim Ardi

Hervé Carfantan

Simon Lacroix

Antoine Monmayrant

Fast hyperspectral cube reconstruction for a double disperser imager

We consider the problem of hyperspectral cube reconstruction with a new controllable imaging system. The reconstruction with a small number of images acquired with different configurations of the imager avoids a complete scanning of the hyperspectral cube. We focus here on a quadratic penalty reconstruction approach, which provides a fast resolution thanks to the high sparsity of the involved matrices. While such a regularization is known to smooth the restored images, we propose to exploit the system capability to acquire the panchromatic image of the scene, to introduce prior information on the sharp edges of the image, leading to a fast and edge-preserved reconstruction of the image.

I. INTRODUCTION

Hyperspectral (HS) Imaging, that is, acquisition of images with a large number of narrow spectral bands, has applications in numerous contexts. Conventional HS imaging technologies proceed by spatial or spectral scanning: they are not instantaneous and are restricted to the acquisition of static scenes, where the imager is either static or follows precisely controlled motions in the case of spatial scanning.

To circumvent this limitation various snapshot devices were proposed [START_REF] Hagen | Review of snapshot spectral imaging technologies[END_REF]. Some require complex optical systems, which include, for example, a series of imagers mounted behind a series of dichroic mirrors or spectral filters. Others exploit reconstruction algorithms using a spatio-spectral mixture of the HS cube obtained through an optical system defined once and for all, at the design phase. In the latter case, the acquisitions of the three-dimensional structure of the cube made by a two-dimensional matrix, necessarily induces a loss of spatial and/or spectral information. The reconstruction algorithms then require assumptions on the scene, and are expensive in computing resources. In addition, these systems require an accurate spatio-spectral calibration.

In this article we are interested in the reconstruction of hyperspectral cubes using a recently proposed device [START_REF] Mcgregor | Adaptive hyperspectral imager: design, modeling, and control[END_REF]. This imager has a dual disperser configuration analogous to the one presented in [START_REF] Gehm | Singleshot compressive spectral imaging with a dual-disperser architecture[END_REF]. It does not suffer from any dependent wavelength shift effect thanks to co-location property: all the spectral components of a given spatial position (x,y) in the cube are imaged at the same point on the CCD. Besides, it is programmable via the control of a matrix of micro-mirrors (Digital Micromirror Device, DMD). This allows the definition of adaptive acquisition schemes, i.e. that sequentially adapt the DMD configurations to the spectral content of the observed scene. For these acquisition schemes to be relevant, the reconstruction step at each iteration needs to be as fast as possible, ideally in real time. To implement such schemes, it is thus necessary to be able to reconstruct an HS cube from a small number of acquisitions in different configurations, at a reduced computation cost.

Section II presents the device and its optical and numerical modeling. In section 3, we address the reconstruction of an HS cube using quadratic regularization, fast and adapted to the system properties. Due to the size of the sparse matrices involved in the reconstruction process, typically 30.10 6 ×30.10 6 for a 1Mpixel-image with 30 spectral bands, a straightforward inverse is not feasible. In section 4 we expose the implementation of an alternative iterative algorithm CGNE (Conjugate Gradient for Normal Equations) and the speed up gained owing to the characteristics of the device. Finally, we present some reconstruction results obtained by simulating the device to acquire data from a HS cube.

II. DEVICE MODELING

A. Acquisition device

The device described in [START_REF] Mcgregor | Adaptive hyperspectral imager: design, modeling, and control[END_REF], is composed of two symmetrical 4f -lines (assembly of two lenses and a diffraction grating), separated by a DMD placed in the plane of symmetry as shown in Fig. 1. Each of the micro-mirrors of the DMD can be configured to transmit or block the incoming light, thus performing a spatial filtering of the signal. The first 4f -line behaves as a spectrometer with a completely open slit that outputs a diffracted image of the observed scene, which is then spatially filtered accordingly to the micro-mirrors configurations of the DMD. The joint effect of the first 4f -line and the micro-mirrors, is thus equivalent to a spatio-spectral filter completely characterized by the DMD configuration. Finally, the second 4f -line compensates the diffraction introduced by the first one: it cancels the wavelength-dependent shift and untangles the spatial and spectral components reaching the CCD. So, in the particular case where all the DMD mirrors are open, the CCD captures a panchromatic image of the observed scene.

B. Optical modeling

This section describes the different transformations that the image undergoes in the optical system, by modeling the propagation along the instrument of the spectral density of the observed object S 0 (x 0 , y 0 , λ), where x 0 and y 0 represent the location of a pixel on the plane and λ is the wavelength.

We denote by S n (x n , y n , λ) the spectral density of the object at the n-th plane of the device. Owing to the device's geometry each (x, λ) line of the HS cube is independent of others. Therefore, in the following, the variable y will be omitted and the one-dimensional spectral density (along the xaxis) S n (x n , λ) located at the n-th plane, as shown in Fig. 1, will be considered instead. We denote T (x) the transmittance of the DMD:

T (x) K k=1 t k 11 ∆ (x -x k )
where x is the position on the DMD's plane, x k and the boolean variable t k represent respectively the location and the state of the k-th micro-mirror of the DMD, 11 ∆ is a gate function of support [-∆/2, ∆/2], where ∆ is the width of the micro-mirrors, and K is the number of micro-mirrors.

The spectral density S 1 of the scene before the DMD is:

S 1 (x 1 , λ) = δ(x 0 -[x 1 + α(λ -λ c )])S 0 (x 0 , λ) dx 0 = S 0 (x 1 + α(λ -λ c ), λ),
where the central wavelength λ c represents the wavelength of the light propagating along the system's optical axis and the coefficient α represents the spectral dispersion of the system along the x-axis. After passing through the DMD of transmittance T (x), the density becomes:

S 2 (x 2 , λ) = T (x 2 )S 1 (x 2 , λ) = K k=1 t k 11 ∆ (x 2 -x k )S 0 (x 2 + α(λ -λ c ), λ).
As the second 4f -line undoes the wavelength dependent shift, the density S 3 at the CCD can be written as follows:

S 3 (x 3 , λ) = δ(x 2 -[x 3 -α(λ -λ c )])S 2 (x 2 , λ) dx 2 = K k=1 t k 11 ∆ (x 3 -x k -α(λ -λ c ))S 0 (x 3 , λ).
Finally, by integrating S 3 spatially over the width ∆ d of the CCD pixels, and over the spectral domain, we retrieve the intensity received by the pixels. The intensity I(c) measured on the c-th pixels is therefore written, as a function to the observed object S 0 , in the form:

I 3 (x 3 ) = K k=1 t k 1 1 ∆ (x 3 -x k -α(λ-λ c ))S 0 (x 3 , λ)dλ (1)
I(c) = 1 1 ∆ d (x 3 -c)I 3 (x 3 ) dx 3 (2)
Equation ( 1) reflects the co-location property of the instrument, since the intensity I(x 3 ) depends exclusively on the contents of S 0 at the location x 3 . Note that when all micromirrors are transmissive (t k = 1, ∀k), the data corresponds to the integration over all wavelengths of the observed object's spectral density, thus to its panchromatic image.

C. Matrix modeling

Considering the linear relationship between the acquisition I(c) and the HS cube S 0 , a matrix representation can be established via discretization. Hence considering an approximation of S 0 (x 3 , λ) by a staircase function of two variables, constant spectrally and spatially along intervals of width ∆ and ∆ d respectively, one can write:

I(c) ≈ ∆∆ d α K k=1 t k S 0 (c, λ k ). (3)
For the sake of simplicity ∆ and ∆ d are considered to be equal in the following, without any loss of generality. For a given configuration of the DMD, relation 3 can be written in matrix form I = To, where I is a vector containing the value of all the C pixels of a CCD column, o is the discretized version of S 0 with the same spatial size as the CCD i.e. C positions, and W wavelengths. Due to the co-location property, T is highly sparse which will be shown to be valuable in the reconstruction process.

Our objective is to reconstruct the HS cube from N image acquisitions I (n) each done with a different DMD configuration T (n) . Fortunately the matrix representation for multiple acquisitions can be readily drawn from the single acquisition model. To extend the first model, we concatenate all the N acquisitions I (n) to form a vector d of dimension N C and likewise for the DMD configuration matrices T (n) to form a matrix H of size N C × CW :

d =   
I (1) . . .

I (N )    =   
T (1) . . .

T (N )    o = Ho. (4)
Finally, to consider the general reconstruction case, that is, for a multiple rows CCD, we simply reiterate the concatenation process to obtain the multi-acquisition model. Note that this extension is possible thanks to the the line independence property of the device. Thus for a R rows CCD, the dimensions of d, H and o are NRC, NRC × RCW and RCW respectively.

III. HYPERSPECTRAL IMAGE RECONSTRUCTION

A. Reconstruction using a quadratic penalization

The reconstruction problem is a linear inverse problem, which consists in recovering the ground truth object from the data. We seek to retrieve the HS cube o of the observed scene from the noisy data d with a small set of acquisitions N : N < W . As a consequence, the system (4) is underdetermined, and additional hypotheses are necessary to solve it. Regularization have proved to yield efficient solutions to such inverse problems.

A classical approach to regularization is to minimize a penalized cost function of the form:

o = arg min o ||d -Ho|| 2 + Ω(o), (5)
The first term of this cost function (which is related to the likelihood function in the case of additive white Gaussian noise) aims to reconstruct a solution compatible with the data, whereas the second term favors certain properties of the object. The choice of the function Ω is a cornerstone of the regularization process: it should be based on prior information on the observed object. Classically, such information can be modelled with a stochastic model taken into account in a Bayesian framework (see e.g. [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]). In image restoration or reconstruction, the challenge is to model piecewise smooth images, to preserve the edges of the images. An alternative is based on assumptions of sparsity of the object in a certain representation space (see e.g. [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]). Recent works based on the sparsity of the object in a wavelet basis has shown its effectiveness for the reconstruction of hyperspectral cubes from devices close to ours [START_REF] Arce | Compressive coded aperture spectral imaging: An introduction[END_REF]. Despite many research on efficient optimization algorithms specific to such problems, they always require more computational resources than using a simple quadratic function for Ω (historically known as Tikhonov regularization) which leads to a linear least-squares problem.

As we are interested in fast reconstruction algorithms, we focused on Tikhonov quadratic regularization. In the 3dimension case of hyperspectral cubes this leads to penalize quadratically the spatial and spectral variations of the cube which results in the penalization function Ω(o) = µ x ||D x o|| 2 +µ y ||D y o|| 2 +µ λ ||D λ o|| 2 , where D x , D y and D λ represent the finite differences along the spatial (x and y) and spectral (λ) dimensions respectively, and µ x , µ y and µ λ their associated regularization coefficients. In this case solution (5) admits an analytic expression:

o = (H t H + µ x D t x D x + µ y D t y D y + µ λ D t λ D λ ) -1 H t d (6
) However, this solution has a major flaw, related to the quadratic penalty: it smoothes the image at the edges, resulting in a loss of spatial and spectral resolutions. To overcome this defect, many works have proposed to replace quadratic regularization with edge-preserving non-quadratic regularization [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]. However, computing the associated solution requires more computation time. Some iterative algorithms designed to solve these problems, sometime called Majorize-Minimization algorithms, are based on computation of solutions equivalent to (6) at each iteration (see e.g. [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF]).

In our case, we can benefit from the properties of the imaging device to detect the edges on the panchromatic image of the observed scene, easy to acquire. A simple segmentation algorithm can locate the edges in this image, which can be accounted for in the quadratic regularization scheme. In practice, this is simply implemented by removing the rows of matrices D x and D y corresponding to actual spatial edges, thereby preserving the quadratic regularization only in homogeneous areas of the panchromatic image. Note that such a procedure amounts to relax the smoothness constraint near the edges and does not force edges for all spectral bands.

B. Implementation

Due to the very large dimensions of matrices, the direct matrix inversion involved in solution ( 6) is intractable, even with the highly sparse matrices H, D x , D y , and D λ . Consequently, we propose to use an iterative algorithm, re-writting equation ( 6) as the following normal equation:

A t A o = A t d (7) with A = [H, √ µ x D x , √ µ y D y , √ µ λ D λ ] t .
Thus the CGNE (Conjugate Gradient for Normal Equation) algorithm can be used to compute the solution with a very low computational cost. Each CGNE iteration require mainly to compute the image of vectors by A and A t which boils down to compute the images by matrices T (n) , D y , D x and D λ and by their transposes. The images of the direct and transposed finite differences operators are straightforward to compute as simple spatial or spectral first order difference of pixels in the HS cube. Thanks to the co-location property of the imaging device, the image of operator T (n) , which corresponds to solving the forward model of eq. ( 3) for each spatial pixel, can be computed as a weighted projection of the object (with 0 or 1 weights): this requires less than N WRC additions. The image of the adjoint operator T (n) t is even easier to compute as it corresponds to a kind of retro-projection, without any addition or products. The CGNE algorithm is guaranteed to converge in at most WRC iterations but, in practice, it converges in less than 110 iterations for a precision of 10 -6 .

IV. RESULTS AND DISCUSSION

In this section we present some reconstruction results from data simulated with a real data cube of dimensions 820×820×31 excerpted from [START_REF] Foster | Frequency of metamerism in natural scenes[END_REF], illustrated Fig. 2. The CCD acquisitions I (n) used for the reconstruction are simulated using equation [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]. For the sake of simplicity, we considered the spectral sampling of the cube as constant and equal to ∆/α, and equal widths for the DMD and the CCD. Hence each unitary spectral displacement ∆/α corresponds to a spatial shift of one pixel on the CCD and the size of the acquired images I (n) is 820 × 820. We exploit acquisitions performed with random DMD patterns where each micro-mirror state (transmission or rejection) t k follows a Bernoulli distribution of parameter P (t k = 1) = 10%, that is, on average 10% of mirrors are in a transmission state. A white Gaussian noise is added to the data, with a signal to noise ratio (SNR) of 20dB.

Some reconstruction results obtained using 5 acquisitions for random DMD patterns are shown Fig. 3. One example of such acquisition is shown in Fig. 3 (e). Note that this figure 3 only shows the reconstruction related to the "ROI1" of 300 × 150 × 31 pixels (see Fig. 2), while the reconstruction has been performed on the whole image. The regularization parameters are set to (µ x , µ y , µ λ ) = (5, 5, 0.5) in order to favor the spectral resolution compared to the spatial one. The spatial resolution has been improved by accounting for the edges, as explained in § III-A. The edges have been detected from the panchromatic image using a simple spatial gradient thresholding method.

The images at a given wavelength of Figs (b), (c), and (d) as well as the spatial/spectral sections of Figs (f), (g), and (h), help to evaluate the quality of the reconstructions, with and without accounting for the edges detected on the panchromatic image. The improvement given by the edge preserving approach in terms of spatial resolution is obvious, in particular in the neighborhood of the edges, comparing the (x, y) crosssections of Fig (d) to the one of Fig. (c) which is blurred. Such edge-preserving effect is also clearly visible on the (x, λ) cross-sections of Figs (h), and (g). However, the improvement in terms of spectral reconstruction quality is more difficult to appreciate from such figures. Nevertheless, such improvement is clearly visible Figs. (i) and (j), where the amplitude of the spectra is clearly better estimated in the neighborhood of the edges when the latter are taken into account. However, such a result inherently depend on the performance of the edge detection process since over uniform regions (far from the edges) the reconstruction is nearly the same for both methods.

As can be seen near x = 275 in Fig. (i) and near x = 190 in Fig. (j), where no edge have been detected, the edge-preserving method fails to improve the reconstruction result. Finally, the distribution of the reconstruction errors helps to compare the overall reconstruction results without or with taking into account the edges. Such a distribution computed on the entire scene (Fig. (k)) is slightly more peaky when the edges are taken into account. But the improvement is clearly more visible if these distributions are computed in a region with a lot of edges as "ROI2" (defined Fig. 2) as shown Fig. (l).

V. CONCLUSION

In this paper we have proposed a fast HS cube edgepreserving reconstruction method that takes advantage of the main characteristics of a double disperser hyperspectral imager controllable thanks to a DMD, namely the co-location property and the ability to acquire a panchromatic image. Using only a few acquisitions (compared to the number of wavelenghts), for random configurations of the DMD, we obtained very satisfactory reconstruction results, in particular when the edges detected from the panchromatic images are taken into account. Indeed, the original and the reconstructed HS cubes are very similar despite the lack of information given by the data. Note that the degradation of the resolution due to quadratic regularization have been partially compensated by the edge preserving approach.

The influence of the acquisitions parameters (number of acquisition, ratio of transmission on the DMD, SNR) has to be studied before using such a method on real data. This requires the definition of metrics to evaluate the reconstruction quality in particular in terms of spatial and spectral resolution.

Our ultimate goal is to exploit the instrument in an adaptive way, in which the DMD configuration of each acquisition is controlled in order to improve the information on the observed object. Therefore, the ability to quickly reconstruct the HS cube is important to define adaptive acquisition schemes.

Figure 1 .

 1 Figure 1. Principle of the device for one line (x, λ) of the HS cube. Lines are independent.

Figure 2 .

 2 Figure2. Illustration of the observed object. Up: RGB image with two regions of interest "ROI1" and "ROI2" denoted respectively by the green and yellow frames, and panchromatic image of the "ROI1" where the green and cyan segments represent respectively the vertical and horizontal edges. Down: RGB intensities along the pink line, and associated panchromatic intensities where the red dotted lines mark the positions of detected edges -the shaded region corresponds to the "ROI1".

Figure 3 .

 3 Figure 3. Simulation results. (a) RGB image of "ROI1". (b) (x, y) cross-section of "ROI1", for band = 10, and its reconstruction (c) without and (d) with accounting for edges. (e) One simulated CCD acquisition of "ROI1". (f) (x, λ) cross-section of "ROI1", for y = 360 (pink line of (a)), and its reconstruction (g) without and (h) with accounting for edges. (i) and (j) Cross-sections of the figures (f), (g) and (h), at band=10 and 25 respectively. (k) and (l) Distribution of the reconstruction error computed over the entire scene and restricted to "ROI2" respectively.

Acknowledgment: The authors thank the Toulouse University and Occitanie Région for funding the Sir4Hyadim Project.