Yohann Béarzi

Julie Digne Raphaëlle

Chaine Liris

Wavejets: A Local Frequency Framework for Shape Details Amplification

Keywords: (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Curve, surface, solid and object representations

Detail enhancement is a well-studied area of 3D rendering and image processing, which has few equivalents for 3D shape processing. To enhance details, one needs an efficient analysis tool to express the local surface dynamics. We introduce Wavejets, a new function basis for locally decomposing a shape expressed over the local tangent plane, by considering both angular oscillations of the surface around each point and a radial polynomial. We link the Wavejets coefficients to surface derivatives and give theoretical guarantees for their precision and stability with respect to an approximate tangent plane. The coefficients can be used for shape details amplification, to enhance, invert or distort them, by operating either on the surface point positions or on the normals. From a practical point of view, we derive an efficient way of estimating Wavejets on point sets and demonstrate experimentally the amplification results with respect to noise or basis truncation.

Introduction

Many shape processing methods, whether they target shape segmentation, shape denoising or shape editing, rely heavily on surface derivatives estimates. Surface derivatives are indeed useful to estimate important shape features such as normals or curvatures. The signal processing viewpoint is slightly different: instead of analyzing signal derivatives, signals are often processed by using a frequency analysis and by devising filters operating on the Fourier coefficients. In this paper, we propose to bring together these two trends in a new function basis taking into account both the local surface derivatives and the angular oscillations of the surface around each point of the surface. We consider surfaces as smooth manifolds that can locally be expressed as a height field over a planar parameterization. In this setting, we locally analyze the surface by extending the osculating Jets [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] in order to take into account both the local angular oscillation frequencies, and their evolution with respect to an increasing radius. This formulation, which we term Wavejets, gives valuable information on the shape by emphasizing that the local behavior of the surface along an arbitrary radial direction from p p p is a polynomial function of the distance to p p p in the parameterization plane. We propose to compute the Wavejets using the tangent plane for parameterization purpose and we demonstrate some theoretical properties of the Wavejets. In particular, we quantify the stability of the decomposition for a small deviation of the parameterization plane with respect to the tangent plane. This result is useful for estimating Wavejets on surfaces described by point sets. Furthermore, Wavejets coefficients can be used to compute in-teresting indicators of differential volumes that can be an alternative to using curvatures and further surface derivatives. Those indicators are used to devise efficient surface filters. We demonstrate two applications of these filters working directly on point sets: position filtering for detail modification (such as detail amplification or detail inversion) and normal modification. Figure 1 shows an example of detail amplification on the Armadillo shape. To summarize, our contributions are:

• A local frequency framework for representing a surface, whose precision and stability are proven. • A practical method for computing Wavejets for point set surfaces

• Efficient shape processing filters for detail amplification, inversion and warping using Wavejets coefficients.

Related work

Surface derivatives estimation. Computing surface derivatives has raised a lot of work in the geometry processing community. [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] introduced the notion of osculating n-jets to analyze a surface around a point. A n-jet is a truncated Taylor expansion used to locally estimate a smooth surface. Each monomial coefficient of a Taylor expansion is directly linked to a high order derivative of the surface, which yields information about the normal, local curvature, or higher order differential quantities. By fitting an osculating jet of order K to a set of surface points in a neighborhood of radius r, the precision of all k order derivatives is o(r K-k ). [START_REF] Tarrou | Invariant smoothness measures for surfaces[END_REF] first introduced a frequency interpretation of the local variations on a surface. Given the principal curvatures at one point, a second order smoothness measure is defined as the integral over a circle of a second order polynomial. A similar process is proposed to define a third order smoothness measure. Following this work, [START_REF] Joshi | An intuitive explanation of third-order surface behavior[END_REF] proposed to interpret third order derivatives of a surface as Fourier series coefficients of the height function above a circle in the tangent plane. Interpreting high order derivatives as Fourier coefficients is interesting since it makes the choice of the origin vector of the local parameterization plane irrelevant.

In many cases, however, surfaces are known only through a set of discrete and potentially noisy measures obtained by a 3D acquisition system. Local surface derivatives are directly impacted by this noise since the noise is amplified through an explicit derivation process. To alleviate this effect, one can compute some local differential quantities by integration. Integral invariants build on this principle to define principal curvatures based descriptors [MCH * 06, PWY * 07, PWHY09]. Integral invariants can be computed directly on meshes or on grids, or on point sets using local surface regression or interpolation. Similarly, [START_REF] Digne | Numerical analysis of differential operators on raw point clouds[END_REF] used covariance analysis to compute principal directions and curvatures of point sets in an asymptotically consistent way. On meshes, the derivatives can be estimated by considering the provided connectivity [START_REF] Meyer | Discrete Differential-Geometry Operators for Triangulated 2-Manifolds[END_REF]WMKG07]. Finally, it is also possible to use metrics derived from heat diffusion to give new expressions of surface derivatives [LSW09].

Signal Processing-like approaches. Many approaches have tackled the problem of surface filtering by mimicking standard signal processing algorithms [PG01]. [Tau95] noticed that Fourier basis functions are eigenfunctions of the Laplacian operator on R 2 . Thus, by building a Laplacian operator over a surface and extracting its eigenfunctions, the projection of the point coordinates functions on such a basis gives a spectral decomposition which can be used for designing low pass filters [Tau95, TZG96]. Pauly proposed to use iterative Laplacian smoothing to separate the signal into high and low frequency information and designs Fourier-like filters in this setting [PKG06]. Defining a Laplacian operator on a mesh has been widely discussed. [VL08] proposed a formulation yielding an orthogonal basis: the Manifold Harmonics Basis. Spectral processing can be directly performed in this basis as if it was a Fourier basis, hence low-pass, high-pass or high-boost filtering has a straightforward implementation, provided the shape eigenvectors are known. Spherical Harmonics were also explored as a way to get a spectral decomposition of a shape parameterized on a sphere in a rotation-invariant way [KFR03]. Spherical harmonics are restricted to genus 0 shapes and are a global basis of the shape, but they can be used to design high-or low-pass filters [ZBS04] or for detail transfer [MCAG08]. In this paper, we set up a framework for the local decomposition of a surface using a new basis. This basis is tailored for analyzing the surface around one point, and for computing differential quantities at that point. Decomposing the signal on this meaningful function basis is the core of Wavejets. Wavejets are therefore related to Zernike basis functions [Zer34], an orthogonal basis, introduced for optical lenses analysis and used often for processing images (e.g. [KH90]) or shape retrieval [NK03]. Although a linear relationship exists between Wavejets and Zernike polynomials, the latter do not give direct access to differential quantities, which are provided by our Wavejets basis by construction.

Detail Exaggeration. While the literature for detail exaggeration in image or video processing is large (e.g. [LTF * 05, DMIF15]) it has been far less studied for surfaces. Existing methods can be sorted out in two categories: the first one modifies the rendering of a shape while the second modifies the shape through its normals or point positions to enhance the shape details. Modifying the rendering of a shape to enhance the details can be considered as a non-photorealistic rendering method in the special case where the purpose is to amplifly the details. Rusinkiewicz et al. [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF] propose to modify the shader in a multiscale way and merge this multiscale information with the classic shading color. Detail Amplification through rendering has been also studied via view-dependent feature computation [VPB * 09]. The second category of approaches explicitly modifies the shape information in a viewpoint-independent manner. In image processing the equivalent operation is done by the so-called High Boost Filter, also known as Unsharp Masking [PRM00, MLLY91, RSKMY96], which moves each pixel value in a direction opposite to the one of the Laplacian smoothing, to sharpen the contrasts. This idea dates back to Gabor [Gab65, LFB94], and is known to produce artifacts if the filter is iterated too many times. This filter can be trivially adapted to point sets, and we will compare our method to it. Another way of amplifying the visualization of details is to modify only the normals at the points. Cignoni et al. [START_REF] Cignoni | A simple normal enhancement technique for interactive non-photorealistic renderings[END_REF] propose to modify the normals iteratively in a direction opposite to the mean normal of the neighbors. This process can seen as a high boost filter applied to the normals. Finally, as a side-benefit of Algebraic Point Set Surfaces rendering [GGG08], a curvature parameter controls the detail enhancement or inversion, a feature that can also be obtained in our Wavejets framework.

Wavejets

Definition

Osculating jets, originally introduced in [CP03] are high order polynomials that correspond to local truncated Taylor expansions of height field surfaces, thus providing interesting differential properties, such as the normal and principal curvatures, for smooth surfaces. We introduce a different kind of jets, named Wavejets, that Preprint Eurographics 2018

φ 0,0 φ 2,0 |φ 2,2 | |φ 3,1 | |φ 3,3 | φ 4,0 |φ 4,2 | |φ 4,4 |
Figure 2: Computation of some values of φ k,n (p p p) for a radius R φ corresponding to the height of the eyes of the statue. If k and n do not share the same parity, coefficient φ k,n is zero. For k ≥ 0 and n = 0, only the magnitude of φ k,n is displayed since the phase depends on the origin direction of the angles. For k ≥ 0, φ k,0 is real, does not depend on the origin direction, and can be negative. Note that φ 1,1 (p p p) = 0 if the parameterization plane is the tangent plane at p p p to the surface. φ 0,0 (p p p) measures the offset of the input point p p p to the surface which can be nonzero in the practical case where the Wavejets order is too low to catch every local variation around p p p.

retain all nice properties of the osculating Jets but also provide a natural interpretation in terms of local angular oscillations.

Let S be a smooth surface and p p p a point on S. As a consequence, the surface can be locally parameterized as a height field f (x, y) inside a neighborhood of radius R φ on a plane P(p p p) passing through p p p. The neighborhood of p p p can be expressed as a Taylor Expansion:

f (x, y) = ∞ ∑ k=0 k ∑ j=0 f x k-j y j (0, 0) (k -j)! j! x k-j y j (1)
where f x k-j y j = ∂ k f ∂x k-j ∂y j . Restricting to a circle of radius r (r < R φ ) centered at p p p in P(p p p), and expressing f with respect to an angle θ yields a periodic function θ → f (r, θ). θ is measured with respect to an arbitrary origin direction for the phases in P(p p p).

Using polar coordinates (r, θ) with (x, y) = (r cos θ, r sin θ) in equation 1 and Euler's formulas to express cos(nθ) and sin(nθ) as polynomials of cos θ and sin θ, one can show that

f (r, θ) = ∞ ∑ k=0 k ∑ n=-k r k φ k,n e i i inθ = ∞ ∑ n=-∞ ∞ ∑ k=|n| r k φ k,n e i i inθ
(2)

with φ k,n = ∑ k j=0 1 j!(k-j)! b(k, j, n) f x k-j y j (0, 0). b(k, j, n) is defined as follows : • b(k, j, n) = 0 if k and n do not have the same parity • b(k, j, n) = 1 2 k i i i j ∑ n-k 2 h=0 k -j h j n-k 2 -h
(-1) h otherwise (see the supplementary material for the full derivation).

This amounts to decomposing the function on a new function basis B k,n (r, θ) = r k e i i inθ .

In other words, while osculating jets provide arbitrary high order derivatives f x k-j y j (0, 0), our representation provides arbitrary high order Fourier coefficients φ k,n ( f ) which are a linear combination of high order derivatives, combining them in a certain way that favors the independence of the coefficients with respect to the origin direction for the phases. More precisely, a rotation of the origin vector in the parameterization plane P(p p p) induces a phase shift of the coefficients. Figure 2 shows the amplitude of the first Wavejets terms φ k,n . Each φ k,n is related to an order of radial derivation k and to a number of oscillations n. A Wavejet of order K is called a K-Wavejet.

Properties

Curvatures By explicitly writing the link between φ k,n and the derivatives of f , the mean curvature H(p p p) and the Gaussian curvature K(p p p) at p p p can be obtained easily.

φ 1,1 = φ * 1,-1 = 1 2 ( fx + i i i fy) φ 2,0 = 1 2 ( fxx + fyy) ; φ 2,2 = φ * 2,-2 = 1 4 ( fxx -fyy + i i i fxy) (3)
Since the Gaussian curvature K(p p p) can be expressed w.r.t. partial derivatives of f at p p p as

K(p p p) = fxx fyy-f 2 xy (1+ f 2 x + f 2 y ) 2 , we get: K(p p p) = 4φ 2 2,0 -16φ 2,-2 φ 2,2 1 + 4φ 1,-1 φ 1,1 2 (4)
Similarly, the mean curvature is expressed as

H(p p p) = (1+ f 2 x ) fxx+(1+ f 2 y ) fyy-2 fx fy fyy 2(1+ f 2 x + f 2 y ) 3 2
, yielding:

H(p p p) = 2φ 2,0 1 + 4φ 1,-1 φ 1,1 + 4φ 2,-2 φ 2 1,1 + 4φ 2,2 φ 2 1,-1 1 + 4φ 1,-1 φ 1,1 3 2
(5)

If P(p p p) = T (p p p), the tangent plane to S at p p p, then φ 1,1 = φ 1,-1 = 0, and :

K(p p p) = 4 φ 2 2,0 -4φ 2,-2 φ 2,2 , H(p p p) = 2φ 2,0 (6)
The principal directions can be found using Wavejets by considering the signal ∑ 2 n=-2 φ 2,n e i i inθ . This signal contains a constant component φ 2,0 and a component that oscillates two times and whose maximum is aligned with the first principal curvature direction (corresponding to the phase of φ 2,2 ). As a consequence, the principal curvatures κ 1 and κ 2 can also be recovered using Wavejets:

κ 1 = 2 φ 2,0 + φ 2,2 + φ 2,-2 and κ 2 = 2 φ 2,0 -φ 2,2 -φ 2,-2 (7)

Stability

We now turn to our main result for the stability of the Wavejets coefficients when the parameterization plane is close but different from the tangent plane. Let us call T (p p p) the true tangent plane and P(p p p) the chosen parameterization plane, also passing through p p p. Since p p p belongs to both planes, they intersect along a line T (p p p) ∩ P(p p p) of direction u. We consider the angle γ such that the rotation of axis (p p p, u) and angle γ transforms P(p p p) into T (p p p). Let us parameterize T (p p p) and P(p p p) so that a point of the surface has coordinates (x = r cos θ, y = r sin θ, h) over T (p p p) and (x = R cos Θ, y = R sin Θ, H) over P(p p p). Let us first assume that θ (resp. Θ) corresponds to the angular coordinate of a point on the surface with respect to u in T (p p p) (resp. with u in P(p p p)). In this setting, the surface Wavejets decomposition at point p p p writes ∑ ∞ k=0 ∑ n=k n=-k φ k,n r k e i i inθ over T (p p p) and ∑ ∞ k=0 ∑ n=+k n=-k Φ k,n R k e i i inΘ over P(p p p). Using the rotation of angle γ, we can express the Φ k,n coefficients with respect to the φ k,n .

We will state our main theorem in this particular setting of origin vector for the phases. To generalize the theorem to an arbitrary origin vector for the angular coordinate θ, recall that a rotation of angle µ of the origin vector in T (p p p) amounts to a phase shift µ. Thus, one can always change the origin vector, compute the Wavejets coefficients φ k,n and recover the Wavejets coefficients for origin direction u as φ k,n e i i inµ (similar formulas hold for Φ k,n and an origin vector change in P(p p p)).

Theorem 1 The coefficients Φ k,n w.r.t to P(p p p) can be expressed with respect to the coefficients φ k,n in the tangent plane T (p p p) as follows:

Φ 0,0 = 0 Φ 1,1 = Φ * 1,-1 = γ 2 e -i i i π 2 + o(γ) Φ k,n = φ k,n + γF(k, n) + o(γ) (8)
where F(k, n) is a function of the φ coefficients of order lower than k.

Proof: see the supplementary material.

Corollary 1 It follows from Theorem 1 that |Φ 1,1 | = 1 2 γ + o(γ) and arg(Φ 1,1 ) = π 2 + o(γ).
Thus if the rotation is small enough, the phase of Φ 1,1 shifted by π/2 in the plane P(p p p) corresponds to the axis of rotation u. Therefore, it is possible to correct the parameterization plane into the tangent plane by performing a rotation of P(p p p) along the axis u with rotation angle 2|Φ 1,1 |.

Proof: see the supplementary material.

Error correction

Corollary 2 One can recover the true coefficients φ k,n iteratively, starting from the lowest order coefficients as:

φ k,n = Φ k,n -γ k-2 ∑ j=1 s j,k,n + o(γ) (9) s j,k,n = ∑ p+m=n |p|≤k-j |m|≤ j φ k-j,p 2i i i (φ j+1,m+1 (m + j + 2) + φ j+1,m-1 (m -j -2)) (10)
In particular,

φ 2,0 = Φ 2,0 + o(γ), φ 2,2 = Φ 2,2 + o(γ), φ 2,-2 = Φ 2,-2 + o(γ)
, which means that the mean, Gaussian and principal curvatures are also stable in o(γ).

Proof: see the supplementary material.

Difference with Jets and Zernike

There exists a linear map between Wavejets φ k,n , Jets J k, j [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] and Zernike polynomials Z n k [Zer34]. This means that there is a way to compute any quantity equivalently from either representation as soon as the linear map is explicit. However this linear map is far from trivial and computing certain quantities will be easier using one or another function basis.

Jets give a direct expression of high-order cross derivatives of a surface. This representation is well-suited to tangent plane estimation, since the cross derivatives of order 1, J 1,0 and J 1,1 , are equal to 0 if the surface is parameterized with respect to the tangent plane. When fitting a jet to a set of surface points, the offset between the parameterization plane and the approximated surface is given directly by J 0,0 .

Zernike polynomials give a polar representation of the surface. Zernike basis is orthogonal, which makes it easy to estimate the coefficients on regular polar grids. However, the error between a parameterization plane and the tangent plane is hard to express using Zernike coefficients. This error can be found as a non-trivial linear combination of Z ±1 k . Similarly, the error offset obtained when fitting Zernike polynomials to a set of points is given by a non-trivial linear combination of Z 0 k .

Wavejets give a representation which is mid-way between Jets and Zernike polynomials. Wavejets coefficients explicitly hold information about the tangent plane. When used to fit a set of surface points, they also provide a direct information on the offset between the surface and the parameterization plane. Besides, they split the surface into components corresponding to different angular frequencies in a similar manner as Zernike decomposition. This angular separation holds important information about the surface, which we will use to design our filters. In particular, we will provide some integral invariants by integration over angle θ, which is straightforward using Wavejets or Zernike and difficult with the Jets. Finally, the error in the orientation of the tangent plane or the offset to the surface are easy to express using Wavejets and Jets but difficult using Zernike polynomials. Thus, Wavejets retain interesting properties from Jets and Zernike polynomials while avoiding some of their weaknesses.

Details filtering and enhancement

Principle

Let us consider the signed volume V V V (s) delimited by the surface and the parameterization plane in a small radius s < R φ around a point p p p. We can express V V V (s) as the sum of infinitesimal angular slices of the volume between the surface and the tangent plane: gent plane, for a radius varying between 0 and s, each point of the surface being scaled by the radius.

V V V (s) = 2π 0 A(θ,
A(θ, s) = s 0 ∞ ∑ k=0 k ∑ n=-k r k φ k,n e i i inθ rdr = ∞ ∑ n=-∞ s 0 ∞ ∑ k=|n| φ k,n r k+1 dr e i i inθ = ∞ ∑ n=-∞ an(s)e i i inθ (11)
Coefficients an(s) correspond to a Fourier decomposition of A(θ, s). By developping the corresponding integrals for the radius s, we obtain the following closed form:

an(s) = ∞ ∑ k=|n| φ k,n s k+2 k + 2 (12)
One can show that if n = 0, 2|a ±n |(s) is the amplitude of the oscillating function s 0 φn(r, θ)rdr, where φn(r, θ), as defined in Figure 3, is the surface restricted to frequencies ±n. For the special case n = 0, a 0 (s) is the static part of s 0 f (r, θ)rdr. Thus, ∀θ, a 0 (s) = s 0 φ 0 (r, θ)rdr (see Figure 3). Each coefficient an(s) has interesting properties regarding the local surface dynamics. 2πa 0 (s) is equal to the signed volume between the surface and the parameterization plane. Indeed,

2π 0 s 0 f (r, θ)rdrdθ = 2π 0 A(θ, s)dθ = 2πa 0 (s) (13)
Importantly enough, the an are intrinsic quantities of the local surface neighborhood, and as such do not depend on the initial parameterization choice.

Remark 1 Equation (13) relates to the Volume Descriptor Vs(p p p) introduced by [MCH * 06] as follows:

Vs(p p p) -2πa 0 ≈ 2 3 πs 3 . ( 14
)
Thus, a 0 (s) measures the local deviation of the surface with regards to the tangent plane at p p p. Local mean curvature is commonly used to reflect this local deviation, but it is meaningless in some cases, for example some points might have 0 mean curvature but nonzero higher orders derivatives and a 0 . Such is the case of (0, 0) for f (x, y) = x 4 + y 4 . At those points, it is necessary to look at higher order derivatives to reveal the local dynamics of the surface. a 0 (s) involves higher order φ k,0 with k > 0, whereas local mean curvature is only proportional to φ 2,0 . In addition, the mean curvature is a measure per point, whose estimation precision from point sets is theoretically controlled by the polynomial order, therefore it will capture extremely small variations. a 0 (s) will on the contrary capture variations at a scale controlled by s.

Similarly, a ±1 (s) locally measures a first order antisymmetric tendency of the surface to change its normal direction when one moves away from p p p, along the direction given by θ. Here, this antisymmetric change corresponds to the variation of the tangent plane to be added on top of the variation induced by the mean variation of the position (reflected by a 0 (s)). More precisely, the normal direction tends to change less when one moves away in the tangent plane in the direction orthogonal to the phase arg(a -1 (s)). Thus, arg(a -1 (s)) gives the direction in which the rotation of the tangent plane is maximal. In a nutshell, it measures around which direction and with what intensity the normal evolves. φ 3,±1 should theoretically be enough to give such information. However, similarly to φ 2,0 , φ 3,±1 is an infinitesimal measure whereas a ±1 (s) gives smoother information while still being able to catch high-order local variations thanks to the high-order estimation of the surface. Figure 3 illustrates the sensitivity of coefficients φ 2,0 and φ 3,±1 compared to a 0 (s) and a ±1 (s) in the 9-Wavejets decomposition of a surface. Using the an instead of the Wavejets coefficients themselves can also be seen as an extension of the use of averaged quantities to compute surface derivatives [MCH * 06, PWY * 07, DM14].

To design the detail enhancement filters, we use Wavejets computed w.r.t. the tangent plane of the surface, and we only rely on a 0 (s) and a ±1 (s). Higher order values a ±n (s) would carry further information of higher order regarding the enclosed area in the direction θ. For example, a ±2 (s) measures the tendency of the evolution of the local horse saddle shape, a ±3 (s) gives information about the evolution of the monkey saddle shape, and so on. The only limit to this computation is the precision for the coefficients when the surface is only known through a set of points and that Wavejets are fit, since, as demonstrated by [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF], the precision of order k derivatives is a o(R K-k Φ ). In such cases, coefficient a ±n (s) is less interesting as n grows.

Position enhancement filter

Algorithm 1 Position Enhancement procedure ENHANCEPOSITION(S, α 0 , K, R φ , s) for all (p p p, n n n) ∈ S do N ←NEIGHBORHOOD(S,p p p, R φ ) φ ←COMPUTEWAVEJETS(p p p, n n n, N , K, R φ ) a 0 (s) ← ∑ K k=2 φ k,0 s k+2 k+2
p p p ← p p pφ 0,0 + 2π(α 0 -1)a 0 (s) n n n

Preprint Eurographics 2018

As noted previously, a 0 (s) gives the local deviation of the surface. If the volume, estimated by 2πa 0 (s) (See Equation 13), was to be increased, local details would be increased as well. Let α 0 be the amplitude of the targeted detail enhancement. In order to increase the volume, a point p p p with normal n n n is moved to a new position p p p such that: p p p = p p p + φ 0,0 -2π(α 0 -1)a 0 (s) n n n (15)

If α 0 = 1, the motion is independent of a 0 (s) and the points are simply projected on the Wavejets surface. Indeed, when the Wavejets order is K < ∞, φ 0,0 measures the distance between p p p and the truncated Wavejets surface. If the Wavejets approximation is exact, φ 0,0 = 0. However, it is not the case in the presence of noise, on sharp edges, or in the approximation framework for point sets described in section 5. The Wavejets surface is however considered as the true underlying surface and each point p p p is moved to this surface by displacing it by φ 0,0 . This has the effect of smoothening edges and canceling a part of the input noise. If α 0 > 1, each point is first moved to the underlying surface and then moved proportionally to its underlying volume. Conversely, if α 0 < 0, it will tend to invert the details and create anti-details.

The computation of a 0 (s) must be performed by using the tangent plane for parameterizing the surface at one point. Therefore, the parameterization plane should be corrected beforehand if it does not exactly fit the tangent plane as explained in section 3. Algorithm 1 assumes that this parameterization plane update has been performed. Notice that if the surface is smooth, the tangent plane and the enclosed volumes a 0 (s) evolve continuously over the surface. Therefore, the surface evolves continuously through the enhancement filter. In practice, this process cannot be applied to the infinite set of points on the continuous surface, however the filter definition remains valid in the continuous setting.

Normal enhancement filter

We propose a filter enhancing the details by exaggerating only the dynamics of the normals. Recall that a ±1 (s) measures the local balance of the shape by identifying the orientation (its phase) and the intensity (its absolute value) of the antisymmetric evolution of the tangent plane. The normal enhancement procedure amounts to modifying φ 1,±1 proportionally to a ±1 (s) at each point of the surface. The corresponding false normal is then estimated as the normal to a plane obtained by rotation of the parameterization following Corollary 1. Notice that this false normal is not coherent with the real surface anymore. Given the desired detail normal evolution gain α 1 = α * -1 , the value of the new coefficient φ 1,±1 is defined as follows:

φ 1,±1 = -π(α ±1 -1)a ±1 (s) (16)
Since φ 1,1 = φ * 1,-1 , it is enough to compute either of these coefficients and deduce the false normal using Corollary 1.

Setting α = 1 leaves the surface unchanged. If α 1 = 2, normals are enhanced, increasing the contrasts. If α 1 = 0, the dynamic of the normals is totally compensated, and the surface looks smoothed out. If α 1 < 0, normals are modified and "anti-details" appear in the rendering. α 1 can also take imaginary values, which skew normals towards one direction as shown in figure 15. Note that, when ). Notice how φ 1,1 , which is 0 without noise, captures most of the noise when the normals are perturbed. φ 0,0 = 0 because the order of a 3-Wavejets is too low to catch all the details of this shape given a large radius, leaving details as residue in φ 0,0 .

applying the normal enhancement filter, the positions of the surface points are unchanged. However since rendering algorithms are more sensitive to normals that positions, it looks as if the positions had been modified. Algorithm 2 sums up the normal enhancement filter.

Algorithm 2 Normal Enhancement procedure ENHANCENORMAL(S, α 1 , K, R φ , s) for all (p p p, n n n) ∈ S do N ←NEIGHBORHOOD(S,p p p, R φ ) φ ←COMPUTEWAVEJETS(p p p, n n n, N , K, R φ )

a 1 (s) ← ∑ K k=3 φ k,1 s k+2 k+2 φ 1,1 ← -π(α 1 -1)a 1 (s)
Compute the false normal n n n given by φ 1,1

Application to point sets

Wavejet Decomposition Equation

Given a surface S that is only known through a set of measured points possibly spoiled by noise, we want to compute the Wavejets representation of the underlying surface up to a chosen order K, at an input point p p p. Let us assume that the surface is locally sufficiently smooth, i.e. C K in a neighborhood of radius R φ around p p p.

Our goal is to compute the φ k,n coefficients that best decompose the underlying surface on the basis functions B k,n (r, θ) = r k e i i inθ in the neighborhood N (p) of radius R φ of p. Let L denote the number of samples in this neighborhood, and let q l be one of these samples, with cylindrical coordinates (r l , θ l , z l ) w.r.t. an axis that corresponds to a rough approximation of the normal direction at point p. Then, the decomposition problem is formulated as finding φ k,n minimizing:

E(Φ) = L ∑ l=1 z l - K-1 ∑ k=0 k ∑ n=-k r l k e i i inθ l φ k,n 2 2 ( 17
)
Preprint Eurographics 2018

For clarity, we state the problem using the 2 norm even if it is unreliable if there are outliers. When dealing with noisy point sets or outliers we solve this minimization using an iteratively reweighted least squares procedure. The weighting scheme involves the use of a diagonal matrix of weights W that are used to leverage the importance of emerging outliers within the norm.

To reformulate this energy minimization, let us reorder the B basis functions B k,n into a vector V b , and the B unknown φ k,n into a vector Φ, such that the b th component of Φ corresponds to the coefficient of the b th basis function in the decomposition. Let k b and n b respectively denote the order of derivation, and the oscillation frequency of the b th basis function (0 ≤ b ≤ B -1). Using these notations, the energy to be minimized is the following:

E(Φ) = L ∑ l=1 z l - B ∑ b=1 r l k b e i i in b θ l Φ b 2 (18)
This amounts to the minimization of MΦ -Z 2 2 , where Z is a vector of size L containing the heights z l of neighbors q l and M is a matrix of size L × B such that:

M l,b = r k b l e i i in b θ l (19)
Minimizing MΦ -Z 2 2 , is done by a QR decomposition of M. Thus computing the Wavejets decomposition around a point p p p amounts to building matrices M and Z and performing the QR decomposition of M. Using a Cholesky decomposition instead of QR fails because MM * is often ill-conditioned.

Algorithm

In order to compute Wavejets in the tangent plane, one can compute a first estimate of Wavejets in a parameterization plane close enough to the tangent plane. The initial parameterization plane is obtained through a Principal Component Analysis. Importantly enough the orientation is necessary neither to compute the Wavejets decomposition nor for the filtering, we only need a local parameterization with respect to an approximate tangent plane. Then the parameterization plane is corrected into the tangent plane using Corollary 1 and the Wavejets coefficients themselves are corrected using Corollary 2.

Given a point set of N points, K the Wavejets order and L the number of neighbors, the complexity of the computation, using Equation 19 for all points is O(NLK 2 ). To be able to solve the equation we pick L ≈ K 2 , yielding a final complexity of O(NK 4 ). As a consequence, when K is big (i.e. 13 for example), the computation cost increases a lot. Once the Wavejets decomposition is computed, the filtering amounts to computing a sum of K terms for each point, because the filter only involves coefficients of frequency 0 or 1 (K/2 coefficients instead of K 2 ). Then each point is moved in its normal direction, which is constant in time, thus the filtering complexity reduces to O(NK).

Experiments and comparisons

Our algorithm was implemented in C++ using Nanoflann and Eigen libraries, and parallelized using OpenMP. Detail enhancements is performed by first computing the Wavejets decomposition and then

Original

Enhanced normals Enhanced positions First column: naive Unsharp Masking with a precise mean curvature computed using 8-Wavejets. Unsharp Masking is unstable when applied to a precisely computed curvature. Second column: our position-based detail enhancement algorithm (2-Wavejets). This filtering is equivalent to a coarse Unsharp Masking since a 0 is proportional to the mean curvature when K = 2. Third column: result for K = 8. Increasing neighborhood radius makes the procedure less sensitive to small local variations while still catching larger details with a high resolution.

α 0 = 3/R φ .
applying the detail enhancement filters. Table 1 gives the parameters and computation times for different models. The computational bottleneck lies in the Wavejets decomposition. Although the computation times are already good enough for processing common 3D point sets, the implementation could be made faster by porting the code to GPU.

In our implementation, for numerical reasons, the local surface is rescaled before computing the Wavejets decomposition so that

α ±1 = 2 α ±1 = ±2i i i α ±1 = -2 α ±1 = ∓2i i i α ±1 = 2 α ±1 = ±2i i i α ±1 = -2 α ±1 = ∓2i i i α ±1 = 0 α ±1 = 1 α ±1 = 2 α ±1 = 4 α ±1 = 0 α ±1 = 1 α ±1 = 2 α ±1 = 4
Figure 8: Influence of normal amplification gain α ±1 and of order K. Similarly to Figure 6, a high order yields a fine enhancement. The phase of α ±1 sets the orientation followed by normal amplification. When α ±1 = 0, the normals are blurred.

R φ = 1. Sometimes, locally, the surface does not project well on the tangent plane, this corresponds to locations where the surface cannot be expressed as a height field over the parameterization plane in a neighborhood of radius R φ , but it would project nicely if the radius was smaller. This happens for example on the fingers of Armadillo on Figure 13. In such situations, Wavejets tend to have high amplitudes and so do the an coefficients, leading to high amplitude motion. To alleviate this problem, we set a threshold on the motion amplitude. In our experiments this threshold is set to Rφ 2 . In all our experiments we set the radius s for computing a 0 and a 1 equals R φ .

The input to our algorithm is a point set with coarse normals which can be computed with a Principal Component Analysis. If the normal orientation is provided, this information is used to re- Table 1: Computation times for various point sets and Wavejets orders. N is the number of points (desktop computer with 2 Intel Xeon E5-2623 of 3.00GHz processors and 32GB RAM). A lower order induces a coarser amplification. When α 0 < 0, the shape tends to be carved in the orientation opposite to the details (visually similar to normal filtering with α ±1 < 0, Fig. 8).

α 0 = -2 α 0 = -1 α 0 = 0 α 0 = 1 α 0 = 2
strict the neighborhoods to points having normals deviating up to π 2 with regard to the normal of the center point. This is done to avoid capturing two different surface sheets in a neighborhood.

The stability of the Wavejets coefficients w.r.t. normal error is shown in Figure 4. As can be seen, adding noise on the input normal direction mostly impacts φ 1,±1 , while order 2 coefficients are very stable, which confirms the results of Theorem 1. In order to test the stability of the a 0 coefficients, the position enhancement filter was tested on a perfect sphere of radius 1 randomly sampled (600k points). The parameters for the filter were R φ = s = 0.03, K = 8 and α 0 = 100. The resulting shape is, unsurprisingly, a sphere with a larger radius. The radius of this enhanced sphere, computed by measuring the average distance to the center, is 1.07 with a standard deviation σ = 3.8.10 -5 , meaning that our computation is quite stable. The resulting sphere is smooth and does not exhibit any bumps. The theoretical equivalence between Wavejets and the Zernike polynomials was also tested numerically. Given a Wavejets and a Zernike decomposition of the same surface w.r.t. the same parameterization plane, the average distance between the 2 reconstructed surfaces is around 10 -15 , i.e. of the order of numerical precision.

Figure 5 shows the effects of normal and position based detail enhancement filters on the Bunny. The exaggerated features have a large scale, for example the borders of the ears and the fur features. On Figure 6, one can see the influence of the chosen order K on the detail amplification. As the neighborhood radius increases, high order Wavejets allow to still amplify fine details, whereas low Preprint Eurographics 2018

Shape

Ground truth Wavejets APSS, scale 2 APSS, scale 4

Figure 10: Estimated mean curvature on a shape with varying sparsity and curvature. The Wavejets order is set to 15 but it decreases in sparse areas to adapt to the low number of neighbors. We compare it with APSS [GGG08] (Meshlab implementation, spherical parameter set to 1). The equation of the shape is ρ(θ, ϕ) = 1 + 0.005 cos(100ϕ 2 /π) in spherical coordinates. In sparse areas (left part of the shape), our method tends to underperform (due to the necessary order change), similarly to APSS (scale 4). However, APSS (scale 2) performs well in sparse regions but is less accurate in dense regions than ours (right part of the shape).

Dense Sparse order Wavejets tend to blur the shape. Since high order Wavejets give similar results over radius variation, a large neighborhood can be chosen to enhance details while being resilient to position noise. This holds as long as the height field assumption is respected. Figure 8 shows the influence of K on normal filtering. The higher K, the more precise the output is. One can note that applying a normal enhancement procedure with α ±1 = 0 cancels local normal dynamics. After normal update, a ±1 is close to 0, which means that local normals dynamics become low. Setting α ±1 ∈ C locally twists normals (see Figure 15) and gives a twisted look to local features. The nose and lips of the mask statue are good examples of this twist (first and second rows of Figure 8). Similarly, Figure 9 shows the influence of K and α 0 on position enhancement. Setting α 0 < 0 tends to carve the shape in the opposite direction of the details. Note that α 0 = -2 gives similar visual results as α ±1 = -2 (see Figure 8 and9).

Input α 0 = 2 α ±1 = 0 α ±1 = 2
Figure 12: Influence of noise on detail enhancement: Noise was added to the positions of the input points. R φ corresponds to 4.5% of the shape diameter. Position filtering with α ±1 = 0 gives very stable results. Position enhancement has been performed using preprocessed normals filtered with α ±1 = 0 so the direction of points motion remains stable over the shape.

We compare our method to a baseline high boost filter, by computing the mean curvature using Wavejets, and moving each point in a direction opposed to the normal direction at a rate equal to the mean curvature (Figure 6, first column). This naive amplification is very sensitive to noise or to small local variations because of the high resolution of the mean curvature computed with a high order K. In figure 7, we show Unsharp Masking applied to the mean curvature computed with Algebraic Point Set Surfaces (APSS) [GGG08] for different parameters. For high scales, the filter tends to produce similar outputs as our method with a low order K. For a low scale, it tends to enhance very small variations. There is a range of scales in which it produces similar results as our implementation using a 0 or a low order K with a small radius (see Figures 6 and7). In Figure 10, we show how Wavejets handle high curvature variations and sparsity compared to APSS for computing mean curvature. We use the PCA regression plane as parameterization plane with R φ as the neighborhood radius. Figure 11 shows the Root Mean Square Error of computed mean curvatures. In dense regions, our estimate starts to fail after frequency 40. Our representation approximates the surface with high precision in dense regions, yielding lower RMSE than APSS. As sparsity increases, Wavejets order decreases, yielding a coarser surface approximation. As one can visualize in Figure 10, APSS performs better than our method in such cases with MLS scale set to 2.

Detail exaggeration filters are more likely to be sensitive to position noise. If noise is considered as a detail to amplify by the algorithm, the output will tend to be noisier than the input. Our algorithm is resilient until a certain amplitude of noise. If the neighborhood radius R φ is large enough, noisy fluctuations in the signal tend to compensate over the radial integration. On Figure 12, different amplitudes of artificial noise over a shape are shown. On Figure 13, the detail exaggeration filter on the position is applied to the

Original

Positions enhanced Normals enhanced Figure 14: On the left: input shape. On the top right: our position and normal enhancement method. On the bottom right: the result of a high-boost filter applied to manifold harmonics of the dino (no normal enhancement is defined using this method). Our position enhancement method highlights more local details as the Manifold-Harmonics based high-boost filter. The rendering of the second row is a capture of the MHB demo provided by the authors.

Armadillo shape. The details are particularly well enhanced by the Wavejets filter: creases are more distinctive and small details, such as the tortoise-shell texture on the legs, are amplified. The eyebrow and teeth are also more prominent.

Figure 14 shows the difference of the output of our position filter with a high-boost filter applied to manifold harmonics. One can see that our result enhances the details in a different manner: the behavior highlighted by our approach is more local and less influenced by global shape structure. Indeed, both methods are very different, since Wavejets filters operate only locally while Manifold Harmonics take the whole shape into consideration for building the function basis. Furthermore, manifold harmonics often require a mesh (although some point clouds generalization exist). The man-ifold harmonics basis does not permit to filter or twist the normals in such a simple way as our method does.

α 1 = 24 α 1 = 24e i i i π 4 α 1 = 24i i i α 1 = 24e i i i 3π 4 α 1 = -24e i i i 3π 4 α 1 = -24i i i α 1 = -24e i i i π 4 α 1 = -24 Original Ours (α 1 = 24) [CST05] [RBD06]
On Figure 16, we compare the result of our normal filter with [START_REF] Cignoni | A simple normal enhancement technique for interactive non-photorealistic renderings[END_REF] and [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF]. [START_REF] Cignoni | A simple normal enhancement technique for interactive non-photorealistic renderings[END_REF] has an effect on the normals that is similar to our algorithm. This algorithm takes as input a neighborhood radius, a number of iterations to perform, and a step. If the number of iterations is too high, the normal enhancement [START_REF] Cignoni | A simple normal enhancement technique for interactive non-photorealistic renderings[END_REF] becomes unstable. On the contrary, our algorithm can perform higher detail enhancement while remaining stable. [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF] has a very different approach that is not based on enhancing normals, but rather on a shading algorithm taking the light direction into account in a multiscale way. Our Enhancement filters can also be used to exaggerate narrow details as shown in Figures 17 and18 on archaeological artifacts. It makes details appear more clearly than on the original point set.

Parameters. The parameters for our algorithm are the following: the radius s for the computation of the volumes a 0 and a 1 and R φ , which we set equal in our tests, the order for the Taylor Expansion K, and amplification gains α 0 or α ±1 depending on which filter is chosen. Note that s = R φ should be chosen so that the number of neighbors L of each point is above

(K+1)(K+2) 2
, which is the number of Wavejets coefficients of order K. The neighborhood size L is an important parameter. If L < (K+1)(K+2) 2

, the system to solve is underdetermined. In such cases, the order K is locally decreased until the system can be safely solved. To do so, we add a parameter β > 1, controlling the decrease rate of the order: K is decreased until L < β (K+1)(K+2) 2

. We set β = 1.1 in our implementation. In presence of noise, β should be greater.

Limitations Our method is still computationally expensive. The computational bottleneck lies in the system solve in order to compute the Wavejets for each point, with a complexity of O(NK 4 ). The Wavejets order choice is thus critical for the computation time, which is coherent with the measures of Table 1. In the presence of fast enough varying features, stability of Wavejets coefficients decreases (see figure 11). Another limitation of our algorithm lies in the assumption that the surface can be expressed as a height field over a parameterization plane inside a neighborhood of given radius r constant over the surface. For some points, r might be too large for the assumption to hold and spurious artifacts can appear. To alleviate this effect, Hamdi-Cherif et al. [START_REF] Hamdi-Cherif A | Superresolution of point set surfaces using local similarities[END_REF] proposed to encode the surface as a height field over a quadric. This simple idea would however require using the Fourier transform on geodesic circles, which is a nontrivial adaptation. Another solution could be to adapt the radius to the local feature size. [MGB * 12] proposes an interesting approach to solve the problem of the right scale estimation that could be adapted to Wavejets as a future work. As an alternative, using multiscale Wavejets could alleviate the scale choice.

Conclusion

We introduced a new basis for decomposing locally a surface in a basis that emphasizes both angular oscillations and surface derivatives. This decomposition is proven to be stable with respect to a small rotation of the parameterization plane from the tangent plane. This allows to define and compute surface characteristic functions at a very high precision and build explicit shape enhancing filters, acting either on the normals or the point positions with a view of amplifying, modifying or inverting shape details. Wavejets can have a broader interest than this application, in a future work we plan to study its application to surface description for shape matching and shape registration.

Original

Enhanced shader [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF] [CST05] Ours (K = 5, α ±1 = 3)

Figure 16: Normal enhancement with different methods. Note that our algorithm has very similar results as [START_REF] Cignoni | A simple normal enhancement technique for interactive non-photorealistic renderings[END_REF] (parameters: are a step of 0.1 with 30 iterations). The main difference between the two procedures is that ours is not iterative, and might be more straightforward for parameters choice. [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF] gives better visual results for this case, but is a totally different method which relies on a mesh where we use unstructured point sets. Each point is rendered through a shading procedure taking the light direction into account while we only rely on the local geometry of the surface. The snapshot of [START_REF] Rusinkiewicz | Exaggerated shading for depicting shape and detail[END_REF] was done using Xshade, a software provided by the authors.

Original Normal enhanced Position enhanced

Figure 1 :

 1 Figure 1: Amplifying details of the Armadillo shape with a Wavejets filter.

Figure 3 :

 3 Figure 3: Wavejets decomposition around a point of a surface. Left: approximated 9-Wavejets surface. Let φ k,n (r, θ) = r k φ k,n e i i inθ + φ k,-n e -i i inθ and φn = ∑ ∞ k=0 φ k,n . The amplitude of the functions φn is more comparable to the amplitude of the input surface than the functions φ k,n .

Figure 4 :

 4 Figure4: Evolution of the coefficients of a 3-Wavejets with increasing Gaussian noise on the normal direction (First row: no noise, 2nd row: Gaussian noise of standard deviation π/15, 3rd row: Gaussian noise of standard deviation π/9). Notice how φ 1,1 , which is 0 without noise, captures most of the noise when the normals are perturbed. φ 0,0 = 0 because the order of a 3-Wavejets is too low to catch all the details of this shape given a large radius, leaving details as residue in φ 0,0 .

Figure 5 :Figure 6 :

 56 Figure 5: Normal and position enhancement on a bunny with 6-Wavejets. R φ is equal to 3% of the shape diameter, and α 0 = α ±1 = 2.

Figure 7 :

 7 Figure 7: Unsharp Masking applied with mean curvature computed with APSS [GGG08] (Meshlab implementation) for different Minimum Least Squares (MLS) scales. MLS spherical parameter is set to 1. Each point is moved in the normal direction with a magnitude of 3 times the mean curvature

Figure 9 :

 9 Figure 9: Influence of position amplification gain α 0 and order K.A lower order induces a coarser amplification. When α 0 < 0, the shape tends to be carved in the orientation opposite to the details (visually similar to normal filtering with α ±1 < 0, Fig.8).

Figure 11 :

 11 Figure 11: Root Mean Square Error (RMSE) of estimated mean curvature with Wavejets and APSS [GGG08] on the shape of Figure 10. The left chart shows the RMSE computed on the part of the shape where the number of neighbors is stricly greater than 150 so K = 15. The right chart shows the RMSE computed where the number of neighbors is between 17 and 31. Wavejets of order 4 are computed in such cases. R Φ is 2.5% of the diameter of the shape.

Figure 13 :

 13 Figure 13: Outputs of our procedures on an armadillo (Left:original; Middle: normal-based detail enhancement with K = 7, α ±1 = 3; Right: position-based detail enhancement with K = 6, α 0 = 2).

Figure 15 :

 15 Figure 15: Normal enhancement on a golf ball. First and Second row : normal amplification for different α 1 on a 9-Wavejet. Note that α -1 = α * 1 on every examples. Last row: comparison with [CST05] normal enhancement algorithm and with [RBD06] detail exaggerating shading.

Figure 17 :Figure 18 :

 1718 Figure 17: Applying order 7 (normal filter) and order 6 filters (position filter) to the Pyramid datasets with α 0 = α±1 = 2.

Preprint Eurographics 2018

Acknowledgements This work was supported by the ANR PAPS project, grant ANR-14-CE27-0003 of the French Agence Nationale de la Recherche.


Nicolás Bayona-Carrillo

Dr Nathalie Bozzolo

Dr Jean-Jacques Fundenberger

Bertrand Thomas

Patrick Camelin

Emmanuel Lenarduzzi

Dr Francis Wagner
email: francis.wagner@univ-metz.fr

Mr Nicolas Bayona- Carrillo

Mr Patrick Camelin

Dr Mr Bertrand Thomas

Effect of

reduced simultaneously. Tailoring of the mechanical properties in cold pilgered tubes of Titanium-base and Zirconium-base alloys is mainly achieved through control of the strain paths during the cold working stages. It is now widely accepted that forming operations with a high ratio of wall thickness to diameter reduction produce a preferential orientation of the basal poles nearly parallel to the radial direction of the tubes which benefits the final mechanical properties. [1,2] This fact has promoted the use of Kearns factors as a common industrial parameter to quantify the anisotropy of cold-pilgered tubes. [START_REF] Kearns | Bettis Atomic Power Laboratory[END_REF] Nevertheless, Ti-3Al-2.5V seamless tubes are mostly employed in the CWSR condition. In industrial practice, some recrystallization might take place during the final heat treatment, which has a major effect on adjusting the mechanical properties (as the yield stress, the mechanical strength or the total relative elongation) to the user requirements. The purpose of this work is to study and explain the impact of the parameters of the final heat treatment on the microstructure and crystallographic texture of the tubes, and consequently, on their final mechanical properties and anisotropy.

Experimental

The Ti-3Al-2.5V (ASTM grade 9) tubes were formed by several cold pilgering passes with intermediate heat treatments. After the final pilgering pass, which involved 66% equivalent strain, the cold worked tubes were annealed for 120 minutes at different temperatures listed in Table 1. In this work, the effect of the heat treatment temperature was characterized through the recrystallized volume fraction (Rx) at the end of the heat treatment. Then, the effect of recrystallization over the evolution of microstructure, texture, tensile properties and anisotropy was analyzed. [START_REF] Kearns | Bettis Atomic Power Laboratory[END_REF] calculated for: RD, the radial direction, TD, the tangential direction and ED, the elongation direction.

For texture characterization, the tubular samples were cut to appropriate size (2 cm long), and then chemically thinned to achieve the final thickness allowing elastic flattening. The macroscopic coordinate system is then: ED -elongation direction, TDtangential direction, and RDradial direction. All the textures were determined from 5 incomplete pole figures measured by X ray diffraction using Co-Kα radiation. The Orientation Density Function (ODF) was computed using the MTEX method proposed by Hielscher and Schaeben. [START_REF] Hielscher | [END_REF] Microstructures were characterized using Electron BackScatter Diffraction (EBSD) in particular to assess the recrystallized volume fraction. EBSD samples were cut along the longitudinal direction and through one of the tube diameters. The surface normal to the circumferential direction was then mechanically polished and finally electropolished using a 10% HClO 4 -90% CH 3 OH solution at 17 V and 10°C. EBSD mapping was performed with a step size of 0.2 µm in a JEOL 6500F FEGSEM equipped with a HKL EBSD system. To ensure the statistical relevance of the recrystallized volume fractions, a minimal area of 18000 µm² was scanned for each sample.

The following EBSD data processing procedure was designed using the Channel5 software. A grain was defined as a group of at least 5 neighboring indexed points having the same crystal orientation within a tolerated angular deviation of 5°. Also, any given grain was considered as recrystallized if its internal mean misorientation was less than or equal to 1.5°, which is the maximum measured value of intragranular mean misorientation for the fully recrystallized sample. The recrystallized volume fraction was then calculated as the fraction of points in an EBSD orientation map belonging to the recrystallized grains.

Tensile tests were conducted at ambient temperature with a cross-head speed maintained at 5 mm/min, so that the initial strain rate was 2.4 × 10 -4 s -1 . CSR measurements were carried out in order to evaluate the anisotropy of the tubes. The CSR is defined as the ratio of circumferential to radial strain after a given tensile strain (3.5%) as described in the SAE AS 4076 specification. [START_REF]Contractile strain ratio testing of titanium hydraulic tubing[END_REF]

Results and Discussion

The estimated recrystallized volume fractions obtained after the heat treatments are displayed in Table 1.

Microstructure Evolution. EBSD pattern quality maps (as evaluated here by the band contrast value provided by the Channel 5 software) allow one to visualize the microstructural changes that occur with recrystallization. Figure 1 shows the pattern quality maps for the 7% Rx, 46% Rx and fully recrystallized conditions. In the 7% Rx condition, the pattern quality is low because of the crystal distortions caused by the forming process. Nevertheless, a few recrystallized grains can be recognized. At the 46% Rx condition, some grains have nearly achieved the final grain size, while other grains are smaller. The size distribution may be local stored energy level. In the fully recrystallized condition, the microstructure is nevertheless quite homogeneous, with an average grain diameter of 4 m. Because of the reduced ductility, CSR measurement was not possible in the cold worked condition. The CSR values for the other conditions show an increase of 30% between the least (7%) and the greatest (100%) recrystallized specimens, as shown in Figure 3a.

The values of yield point, mechanical strength and relative total elongation, measured at different degrees of recrystallization, are plotted in Figure 3b. As expected, the mechanical strength and the yield point of the tubes decreased with increasing recrystallization. Both values were found to be strongly correlated with the recrystallized volume fraction, showing a very good reproducibility of the tensile flow stresses on tubes over tubes after same heat treatment. As for the total relative elongation, the results showed a greater dispersion, this might be explained by complications related to the tensile test; because no notch was made on the tubular samples, necking did not always occur in the gauge length measured by the extensometer. Nevertheless, the observed tendency is the usual increase in ductility which comes with recrystallization, as a consequence of the lower dislocation density. Texture evolution. Figure 4 shows that the textures of all investigated samples (cold worked and annealed) have the orientation density maxima at the same position. The basal poles are spread in the plane normal to the axial direction, with a higher density at 60° from the radial direction, and crystals preferentially exhibit a <10-10> direction aligned with the axial direction. Only the texture strength decreases when recrystallization progresses over 45% Rx. Rx, e) 90% Rx and f) fully recrystallized conditions. The same grey-scale applies for all pole figures.

The texture index gives a description of the sharpness of the texture, it is defined as follows, [START_REF] Bunge | Texture analysis in materials science[END_REF] (

Where f(g) is the orientation density function and g the orientation. The evolution of the texture index during recrystallization is shown in Figure 3a and the calculated values are specified in Table 1. J increases at the very beginning of recrystallization, and then decreases steadily as recrystallization progresses. This result shows that J is a suitable parameter to evaluate the advancement of the recrystallization process in Ti-3Al-2.5V tubes, at least in the particular situation where the texture type does not change.

As illustrated by Figure 4b andc, the rolling texture is reinforced at 7% Rx and 21% Rx which is in agreement with the texture evolution observed in commercially pure titanium at the beginning of recrystallization. [START_REF] Wagner | [END_REF] Further recrystallization leads on the contrary to a texture strength decrease (Figure 4d, e and f), but the final texture remains close to the deformation texture. In particular, the well-known 30°<0001> rotation that describes the texture change which is classically observed during annealing of hcp metals, does not occur here. It has been shown that this texture evolution during annealing of strongly cold rolled low alloyed titanium sheets, leading to the improperly so-called "recrystallization texture", is actually mainly controlled by grain growth. [8] In the present study, grain growth was very limited under the applied conditions, which probably contributes to the texture stability. In fact, Ti-3Al-2.5V contains some dispersed β phase, stabilized by the vanadium content. The β phase is known to reduce the rate of grain growth in alpha or near-α titanium alloys, [9] which may in part explain the low average grain size obtained in the fully recrystallized condition.

The Kearns factors, which give a description of the amount of basal poles aligned in each direction, were calculated from the (0002) pole figures according to Anderson et al. [10] Those calculated values, shown in table 1, do not vary much, and therefore cannot by themselves explain the observed CSR increase.

The low values of the CSR at the initial condition are a consequence of the inclination of the basal pole towards the tangential direction. [1,2] During recrystallization, the evolution of the CSR responds to the same phenomenon observed for the texture index J, also plotted on Figure 3a. The weak value of J in the fully recrystallized condition describes the smoothness of the final texture, which is equivalent to saying that the material has lost much anisotropy, and, consequently CSR comes closer to the isotropic value of 1. This explains the measured tendency of the material towards an isotropic behavior as recrystallization advances.

Conclusions

1. The variations in the CSR taking place during recrystallization are the result of the weakening of the crystallographic texture, which can be quantified by the texture index J.

2. Since the texture changes occurring during recrystallization in alpha titanium alloys do not affect the orientation distribution of basal poles, Kearns factors are not suitable to describe the observed evolutions.

Recrystallization is associated with the appearance of yield point phenomena and

Lüders elongation in strongly annealed Ti-3Al-2.5V tubes.

4. Besides the reported yield point phenomena, the tendencies observed for tensile behavior respond to the usual consequences of recrystallization in alpha titanium alloys.

Fig. 1 .

 1 Fig. 1. EBSD Pattern quality maps of (a) 7% Rx, (b) 46% Rx and (c) fully recrystallized Ti-

Fig. 2 .

 2 Fig. 2. a) Examples of the true stress-true strain curves obtained for the different

Fig. 3 .

 3 Fig. 3. a) Evolution of the texture index and CSR during recrystallization. b) Elongation (El),

Fig. 4 .

 4 Fig. 4. And pole figures of the a) cold worked, b) 7% Rx, c) 21% Rx, d) 46%

Table 1 .

 1 Recrystallized volume fraction, Kearns factors, and texture index values obtainedafter final annealing for 2h at the displayed temperatures.

	Annealing	Recrystallized	f RD	f TD	f ED	Texture
	temperature	vol. fraction	[a]	[a]	[a]	index J
	[K]	[%]
	-	0 (CW)	0.4354 0.4405 0.1241 3.188
	838	7	0.4528 0.4337 0.1134 3.604
	863	21	0.4408 0.4693 0.0898 4.142
	883	46	0.4381 0.4410 0.1209 2.847
	933	90	0.4699 0.3914 0.1387 1.959
	993	100	0.4356 0.4074 0.1570 1.817
	[a] f RD, f TD and f ED are the Kearns factors,


Nadine Candoni

Zoubida Hammadi

Romain Grossier

Manuel Ildefonso

Shuheng Zhang

Roger Morin

Stéphane Veesler

Addressing the stochasticity of nucleation: practical approaches

Keywords: spontaneous nucleation, induced nucleation, microfluidics, confinement, external fields

This chapter presents different practical ways to address nucleation stochasticity. The methods use either statistical studies on spontaneous nucleation or local control of nucleation. Techniques developed in our laboratory are described: droplet-based microfluidics, micro-injectors in oil, external electrical or mechanical fields in confined systems. Results of nucleation kinetics obtained on various molecules, are presented in terms of metastable zone, critical supersaturation, nucleation rate, induction time, interfacial energy of the critical nucleus. These practical approaches show considerable potential to increase understanding and control of the nucleation mechanism.

Introduction

In crystallization of materials such as pharmaceuticals, proteins, or nanomaterials, the most important properties are crystal morphology, habit, size distribution and phases. The control of these properties demands an understanding of the mechanisms governing the successive steps of nucleation and growth. Although a great deal is known about crystal growth, considerably less is known about crystal nucleation because of the difficulty of directly observing the nuclei. Indeed, the first nuclei are defined by Kashiev as randomly generated nanoscopic new phases "that have the ability for irreversible overgrowth to macroscopic size" [START_REF] Kashchiev | Nucleation: basic theory with applications[END_REF]. Hence nucleation is a stochastic phenomenon and the probability of observing the first nuclei is very low, for spatial and temporal reasons.

In this chapter we present practical ways of addressing the stochasticity of nucleation [START_REF] Candoni | [END_REF]. The first idea consists in performing many spontaneous nucleation experiments, and hoping that the law of large numbers permits us to get a sharper picture of nucleation. Observation is enhanced by reducing the volume of crystallization [3], and hence the number of crystals. The techniques presented here are droplet-based microfluidics [4][5][6][7][8][9] and micro-injection in oil [10][11][12][13]. The second idea consists in inducing nucleation where and when we choose controlling nucleation, which is no longer spontaneous. In this case, the techniques combine confinement of nucleation with localized external fields, such as electrical [14][15] and mechanical fields [START_REF] Hammadi | Faraday Discussions[END_REF].

The purpose of this chapter is to show what these techniques currently offer in terms of nucleation kinetics measurements and detection of the critical nucleus. The main questions are: what are the factors that influence nucleation and how they are revealed by these techniques? Since we explore nucleation from a purely practical point of view, we do not consider any a priori assumption on the nucleation mechanism. We will present practical examples with mineral, organic anud biological molecules confirming that nucleation mechanisms are identical for all molecules [START_REF] Chernov | [END_REF].

Theoretical background

The following is a brief reminder of nucleation theory and experimental parameters, introducing all the properties accessible with our techniques.

Critical nucleus

Until recently, solution nucleation has been described solely by the classical nucleation theory (CNT), a theory derived from nucleation of droplets in the bulk of pure supersaturated vapors [START_REF] Gibbs | The Collected Works[END_REF]. It considers that once a cluster has reached a critical size, nucleation starts. Hence the first nucleus is called the critical nucleus. In this model, the critical nucleus appears by densification of molecules and nucleation of a solid phase at the same time. Therefore the transition from the liquid phase to the solid phase is described by one order parameter, i.e. density. According to CNT, the critical size r* is given by the Gibbs-Thomson equation:

equation (1)
With  (m 3 ) the volume of a molecule inside the crystal,  (J•m -2 ) the interfacial energy of the critical nucleus, k the Boltzmann constant (J.K -1 ), T the temperature (K) and  the supersaturation (C/C s ). C is the concentration of the solute in solution and C s the equilibrium concentration. Hence, if >1, the crystal grows, if <1 the crystal dissolves, and if =1 crystals and solution are at equilibrium.

While this theory has the advantage of simplicity, some discrepancies (in predicting nucleation frequencies) have been observed with experiments [3,[START_REF] Dixit | [END_REF][20][21][22]. Consequently, a more complicated two-step mechanism was proposed: first, formation of a dense phase of clusters on the model of a liquid-liquid phase separation and second, organization of these clusters into structured clusters [23][24][25]. This theory is based on two order parameters, density and structure [26]. Subsequently, a third nucleation theory was proposed linked to biomineralization, consisting of a nucleation via stable prenucleation clusters [27]. For reviews see 28, 29. Concerning experimental discrepancies, different explanations are proposed. For instance, Knezic et al. [30] suggested that macromolecules take longer to rearrange themselves into crystalline nuclei than smaller molecules (this rearrangement step is not directly taken into account in CNT). Hence, the limiting step, densification or structuration, will depend on the solute.

It must be pointed out that experimental observations and interpretations are influenced by the fact that the location, the nature and the lifetime of the dense phase and/or the critical nucleus are currently unknown.

Nucleation kinetics

In this chapter we are looking at primary nucleation, i.e. the appearance of crystals in a solution which contains no crystals of that phase. The nucleation kinetics is commonly described by the nucleation rate J, which is the number of nuclei that form in the supersaturated solution per unit of volume and per unit of time [START_REF] Kashchiev | Nucleation: basic theory with applications[END_REF][31][START_REF] Abraham | Homogeneous nucleation theory[END_REF][START_REF] Toschev | Homogeneous nucleation. Crystal growth: an introduction[END_REF]. Theoretically, the nucleation rate J is determined by CNT derived from the Gibbs treatment of droplet nucleation [START_REF] Gibbs | The Collected Works[END_REF]. The general equation for primary nucleation is expressed as equation (2) [START_REF] Boistelle | [END_REF]: equation [START_REF] Candoni | [END_REF] with K 0 the kinetic factor (m -3 •s -1 ) and ΔG* the activation free energy (J).

CNT assumes a spherical form for the critical nucleus as mentioned by Fletcher [35]. In the general case of heterogeneous or homogeneous nucleation, a constant shape for the cluster (here spherical) is assumed, which depends on the possible presence of foreign substances: foreign molecules or particles, bubbles, crystallizer walls, liquid-liquid or liquid-air interfaces. Hence the value of ΔG* depends on surface interactions, which are taken into account by introducing an effective interfacial energy  ef of the critical nucleus and J is given by equation (3) [4]:

equation (3)

Critical supersaturation

According to equation ( 3), there is a zone in the supersaturated region (>1) where no nucleation occurs for a given time in a given volume. This is the metastable zone [36], which is delimited by the critical supersaturation (CS) for a given temperature or the temperature of metastable zone (TMZ) for a given concentration (Fig. 1). In the nucleation zone or zone of spontaneous nucleation, nucleation (i.e. the appearance of crystals) is promoted; in the metastable zone, growth (i.e. the subsequent expansion of the crystal) is promoted.

Experimental considerations

In practice, the two usual methods used to determine the nucleation rate J are:

-Measurement of the induction time, i.e. the time it takes to form a detectable amount of the new crystalline phase after a given supersaturation is established. However nuclei must grow to a minimum size before they are observable; thus the induction time includes the nucleation time and the time required for growth to a detectable size. The supersaturation range over which the experiment can be performed is rather limited due to the metastable zone. When β is too small heterogeneous nucleation may replace homogeneous nucleation. When β is too large, induction times are difficult to measure because they are too short. In practice, we are able to measure the induction time in the vicinity of the metastable zone limit and we could reasonably assume that the growing period can be considered negligible compared to the nucleation period.

-The Double Pulse Technique (DPT), which permits nucleation and growth to be separated [37][38][39]. At the beginning of a run, the solution is loaded at a temperature chosen to prevent nucleation of crystals. Then the temperature is lowered (respectively raised if the solubility is reverse with T) to a selected T 1 at which nucleation occurs. After a period of Δt 1 (nucleation time), the temperature is raised (respectively lowered if the solubility is reverse with T) from the nucleation temperature T 1 to the growth temperature T 2 in the metastable zone. At T 2 , crystals already formed can grow to detectable dimensions. The nucleation rate, at a given concentration C and temperature, is given by the variation of the amount of crystals nucleated as a function of nucleation time.

Studying nucleation, an unpredictable phenomenon

As we have seen above, nucleation plays a fundamental role in the physical properties of crystals. However the unpredictability of the spatial and temporal location of the critical nucleus makes its study difficult. In this part we address the problem of the stochasticity of nucleation in two ways: first by carrying out statistical studies on single crystals spontaneously nucleated; second by inducing nucleation.

Statistical studies of spontaneous nucleation

For statistical studies, a large number of samples in the same condition are required. In large volumes, it is difficult to accurately detect crystals as they are formed, so sample volumes need to be limited. Therefore the idea is to reduce the volume of nucleation to the nL-pL range. Hence the number of crystals nucleating and the zone of observation are reduced, thus increasing the resolution of the detection. Moreover, this has the advantage of decreasing the quantity of molecules while speeding up heat and mass transfer [40]. Volume reduction to nL is realized using droplet-based microfluidics. Then pL-fL volumes are obtained by micro-injection methods.

Nucleation in nanoliter droplets

Droplet-based microfluidics permits series of nanolilter droplets to be generated by injecting the crystallization solution in a non-miscible continuous phase and without using surfactant. Droplets are mono-dispersed (in size, composition, frequency) nano-crystallizers isolated from each other and in which crystallization can occur independently. Hundreds of droplets can be generated per condition, leading to a large number of samples for statistical studies, while consuming only small quantities of material.

-Microfluidic devices: (i) molded in PDMS and adapted from Salmon [41][42][43] and (ii) composed of Teflon tubing with junctions in PEEK (Fig. 2a-b). The continuous phase is silicone or fluorinated oil according to the solvent of crystallization (discontinuous phase). The crystallization solvent is water in PDMS systems. This microfluidic system thus becomes compatible with almost all solvents, since Teflon tubing is used instead of PDMS chips (Fig. 2c). Crystals of various organic molecules obtained in aqueous and organic solvents are presented in figs 2d-2h [6]. The microfluidic set-up [7] in Fig. 3 allows us to determine the properties of nucleation described above theoretically: metastable zone width, nucleation rate J, effective interfacial energy  ef of the critical nucleus and polymorphism. -Metastable zone width: evaluated from measurements of the temperature of the metastable zone (TMZ) illustrated in Fig. 1. This is the maximum temperature above which the supersaturated solution does not lose its metastability over a given time [36]. For a given concentration, TMZ is determined by storing solutions at given temperatures for a given time (here 20h). During this time if T<TMZ, crystals nucleate and if T>TMZ, no crystal appears. Hence TMZ at which there is no crystal is extrapolated (interpolation) from a plot of the average number of crystals versus temperature. Moreover a rapid screening of the solubility diagram (black curve of Fig. 4) can be carried out as shown by Laval et al [41][42]. -Nucleation rate: determined by using the Double Pulse Technique (DPT), based on the phase diagram of figure 4. Crystals are nucleated at a T Nucleation outside the metastable zone (below TMZ) and for a given nucleation time t; then, the crystal are stored at T Growth in the metastable zone, where the crystals already formed grow to detectable dimensions and are counted. From the variation of the average number of crystals nucleated as a function of t, the steady-state primary nucleation rate J is determined for given temperature and concentration. Our results obtained in nucleation volumes of 250nL of lysozyme are summarized in Fig. 5.

Supersaturation

Fig. 5 Lysozyme nucleation rate vs supersaturation, () data from our study at 20°C [5] and () data from Galkin and Vekilov, at 12.6°C [21][22]. Copyright 2011 American Chemical Society.

As shown in Fig. 5, our kinetic data are in good agreement with the data of Galkin and Vekilov [21][22] obtained for the same NaCl concentration at 12.6°C, using droplets of 1µL suspended in Teflon wells. Moreover, the supersaturation range experimentally accessible when reducing the volume of nucleation to the nL-pL range is increased for kinetic [5] and thermodynamic [14] reasons (see part 3.2.).

-Effective interfacial energy  ef computed from eq.3 shows that primary nucleation is not influenced by the container material (Fig. 5). Indeed, the oil creates a "containerless" environment, as pointed out by Chayen [45]. However,  ef determined in fluorinated oil and in silicone oil (Tab. 1) demonstrate the influence of the droplet-oil interface in droplet-based microfluidics [4]. Ildefonso et al. 5 Vekilov group 21-22 Ildefonso et al. 4 Fraden group 46 -Polymorphism: droplet-based microfluidics is adapted for detection of new phases for statistical reasons. For instance, different polymorphs can nucleate, the most stable polymorph corresponding to the minimal free energy of the system. Moreover, according to Ostwald's rule of states [47], a chemical system does not tend directly towards the most stable state, but rather tends towards the closest metastable state. Hence performing many spontaneous nucleation experiments increase the probability of observing an event. In the case of lysozyme, we observed a metastable phase, the sea urchin-like phase, in 6 droplets out of 237 (Fig. 6) [8]. It was already known [48][49], but is not easy to observe in mL crystallizers. Moreover, because of mononuclear nucleation, the nucleated crystal of the metastable phase cannot easily follow a solution-mediated phase transition to a more stable phase. Thus small volumes of droplets seem to ''freeze'' the metastable phase [50]. Lastly, the higher supersaturations reached in microfluidics due to volume reduction increase the probability of nucleating a metastable phase [51].

Fig. 6 Droplets of lysozyme solutions (20 mg/mL, 0.7 M NaCl-pH4.5) in Teflon tubing of 500µm diameter, observed at 20°C after storage for 20h at 6°C: the tetragonal form top left and the sea urchin-like form bottom right [8].

In screening the crystallization conditions of rasburicase, the same experimental conditions produce different crystal habits and phases. Moreover, the two polymorphs known for rasburicase [52] appear simultaneously in droplets (Fig. 7) [9]. This concomitant nucleation occurs mainly because the difference between their solubilities is weak. To summarize, droplet-based microfluidics is an accurate method for studying nucleation kinetics due to statistical approach, exploring one experimental condition in hundreds of identical droplets. Moreover, as droplet volumes are nanoliter, this statistical approach requires only a few milligrams of molecules.

Nucleation in picoliter to femtoliter droplets

Small-droplet generation with micro-pipets or micro-injectors renders micro-droplets accessible, and makes the size range easier to control compared to microfluidics, which is limited by channel size (hardware limitation). However, statistical analysis cannot be carried out with single micro-droplet injection [40,[53][54].

-The fluidic device developed in our laboratory [10] (Fig. 8a) generates arrayed aqueous phase micro-droplets through a layer of inert oil (sessile geometry), without using surfactant. Up to thousands of micro-droplets can be generated, with volumes ranging from nanoliter to femtoliter (Fig. 8b), thus allowing statistical studies of nucleation. -The droplet-contraction method is used to generate supersaturation and thus to create crystals, as we can see with picoliter droplets of NaCl solution in Fig. 9a [START_REF] Hammadi | Faraday Discussions[END_REF]. Indeed, at this scale, droplets slowly evaporate, linearly decreasing the volume (Fig. 9b) [11]. Hence the concentration, which is given by the refractive index [12], increases until it exceeds the critical supersaturation, i.e. the supersaturation required to nucleate a new phase [55]. Then single crystals nucleate (Fig. 9c).

50µm

-The induction time can be measured with the droplet-contraction method because the growth rate is so fast that the time required for the newly-formed nuclei to grow to a detectable size is negligible with regard to the induction time [13]. Thus, the time when a detectable crystal is observed can be considered as the induction time. Moreover, the smaller the volume, the faster the evaporation rate and the faster the nucleation (Fig. 10) [START_REF] Hammadi | Faraday Discussions[END_REF].

Fig. 10 Nucleation probability P(t) vs time. P(t) represent the normalized fraction of nucleated microdroplets for droplet diameters of () 104pL, () 172pL and () 284pL [START_REF] Hammadi | Faraday Discussions[END_REF]. Reproduced from Ref [START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry To summarize, this technique of droplet micro-injection in oil combines the advantage of channeled microfluidic techniques, generating thousands of droplets, with the advantage of micropipette techniques: control over size, micro-droplet accessibility and contraction of droplets in oil to increase supersaturation. Both of these techniques involve easy-to-build devices and are easy to use. Moreover, it must be noted that they require no surfactant that might affect crystallization.

Influence of volume on nucleation

Reducing the volume of crystallization raises several questions: How do kinetic parameters such as the metastable zone and the critical supersaturation vary? Does CNT, which supposes an infinite reservoir of molecules, i.e. a constant supersaturation during the nucleation event, remain valid in small volumes? What are the limitations of small volumes and how small? -A kinetic limitation in small volumes arises from nucleation rate J (number of nucleus per unit of volume and per unit of time). J is an intensive property (does not vary with volume). However J is inversely proportional to volume and proportional to induction time. Hence the smaller the volume, the longer the induction time. For instance, an induction time of 1s in a volume of 1mL becomes 11 days in a volume of 1nL.

-A thermodynamic limitation appears because the critical size and energy barrier required for nucleation increase when the crystallization volume decreases [56]. When the volume decreases, the concentration in solution can no longer be considered constant during the nucleation process but is decreasing. In practice, this means that the reservoir of molecules is finite, in contradiction with CNT. Hence the critical supersaturation, where no nucleation can occur, increases with confinement (Fig. 11), widening the metastable zone [13,56]. Thus "scale-down" requires creation of sufficient supersaturation within each droplet to ensure nucleation, with a risk of unwanted nucleation before droplet formation. In practice there is a "thermodynamic" influence of volume on nucleation from picoliter range down. To summarize, two effects occur depending on droplet volume: 1) for nL range droplets, a kinetic effect increases the induction time; 2) for pL-fL range droplets, the same kinetic effect is accompanied by a thermodynamic effect. In both cases these effects widen the metastable zone, hence requiring higher supersaturation for nucleation. However these effects stabilize high-supersaturated solutions, permitting critical clusters and nucleation to be controlled.

Predictive study of induced nucleation

Nucleation can even be induced in the metastable zone, thus becoming predictive. In other words, we determine where and when it will occur. Nucleation is induced using an external field, the implications of which were highlighted by Voss [57], Oxtoby [58] and Revalor [59]. Two effects on the supersaturated solution are expected: molecular orientation and density fluctuation. As a general rule, any disturbance triggers nucleation once a highlysupersaturated metastable state has been achieved [60]. First, localization in time is achieved by triggering nucleation in a metastable solution with a localized mechanical contact or a localized electrical field. Then localization in space is achieved by either confining the nucleation volume or by controlling fluxes in non-convective gel medium.

Mechanical action and confinement by volume

-Coupling a mechanical action using a sharp tip with volume confinement by microinjector triggers nucleation, as previously shown. A structural transformation is thereby induced in a dense supersaturated droplet at precisely determined points and times (Fig. 12). Then the crystal grows at a rate greater than 200μm/s for the first 20ms [START_REF] Hammadi | Faraday Discussions[END_REF]. [START_REF] Hammadi | Faraday Discussions[END_REF]. Reproduced from Ref [START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry -The induction time is shorter than 5ms, corresponding to the time between tip contact and observation of nucleation (Fig. 12a-b). However the droplet is repeatedly tapped with the tip during the generation of supersaturation by water diffusion in order to launch nucleation. The crystal nucleated is rough and transition to a facetted crystal is observed in less than 1s, i.e. the transition between the nucleation form and the equilibrium form. In practice, the tip position being controlled with micromanipulators, the position of the critical nucleus can be determined with an accuracy of 16nm.

-This method is suitable for a binary crystallization medium, with a solute and a solvent. The solute concentrates during droplet contraction and reaches high supersaturation. Then the solute nucleates due to the mechanical action. When the crystallization medium is ternary containing a solvent, a solute and a crystallization agent, e.g. protein crystallization, both solute and crystallization agent concentrate during droplet contraction and reach high supersaturation. Hence for kinetic reasons [61], the crystallization agent often nucleates first, before solute nucleation. Therefore, supersaturation must be achieved via a method other than droplet contraction.

Electrical field and confinement by gel

-Coupling a localized electric field generated by a sharp-tip electrode with a gel eliminating convection in a supersaturated metastable solution increases supersaturation at the tip apex and induces nucleation. A large electrical field [START_REF] Gomer | Field emission and field ionization[END_REF] and large field gradients are created near the tip at low DC voltage, due to its nanometer size (Fig. 13). The electrical field induces high current density inside the solution close to the region of high curvature [14]. The resulting flux accumulates molecules at a precise point in the vicinity of the tip apex [15] and the concentration gradient in solution tends to be counteracted by convection. However, the solution being in a gel, convection is eliminated and the confinement is kept in the vicinity of the tip until the critical supersaturation is reached. Then the crystal nucleates near the electrode tip [START_REF] Hammadi | Faraday Discussions[END_REF].

-Induction time is measured by optical observation and by the evolution of the current in the solution. An example of nucleation induction in a supersaturated metastable solution of lysozyme in agarose gel is shown in Fig. 14. Fig. 14 Time sequence in a supersaturated metastable solution of lysozyme 20mg/mL (NaCl 0.7M). Nucleation is induced by an external localized electric field with a sharp-tip electrode and in agarose gel 1% gel [START_REF] Hammadi | Faraday Discussions[END_REF]. Reproduced from Ref [START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry

The optical observation of nucleation gives an induction time shorter than 600s. The nucleated crystal is rough, due to high local supersaturation encountered during nucleation and growth. Growth rate (between 600 and 10800s) is 15μm/h in agreement with the growth rate obtained by Durbin et al. [START_REF] Durbin | [END_REF] at high supersaturation for lysozyme. Compared to optical determination of the induction time, the measures of the current during the experiment are clearly more sensitive (Fig. 15) [START_REF] Hammadi | Faraday Discussions[END_REF]. Indeed the induction time is given by the modification in the slope at 300s.

Conclusion

Faced with the stochasticity of nucleation, different experimental approaches have been developed to study nucleation. In the first approach, nucleation is allowed to proceed spontaneously and the stochastic nature of results is addressed by statistical studies. Thus the volume is reduced to droplets of nL to fL volume that are generated in large numbers. Droplet-based microfluidics and micro-injectors in oil are used. These are compatible with many solvents and involve easy-to-use and material-saving techniques. In the second approach, nucleation is controlled by external electrical or mechanical fields in confined systems. Hence nucleation becomes predictive, triggering nucleation for molecules that are not easy to crystallize, nucleating one crystal per crystallizer and localizing nucleation in space and in time for measurements on the critical nucleus.

All of these methods approach nucleation kinetics through the induction time or the nucleation rate. Moreover, factors affecting the nucleation process (confinement, complexity of solution) are brought to light. These practical approaches show considerable potential to increase understanding and control of the nucleation mechanism.

Fig. 1

 1 Fig. 1 Phase diagram of supersaturation vs temperature, with the solubility curve (----) and the metastable zone (••••).

Fig. 2

 2 Fig. 2 (a) Photo and (b) scheme of a T-junction made of PEEK; (c) Photo of isolated droplets of 60nL volume of crystallization solution separated by fluorinated oil in Teflon tubing of 500µm diameter, (d) Lysozyme crystals in aqueous droplets, (e) Rasburicase crystal in a viscous medium, (f) Caffeine crystals in ethanol droplets, (g) Isonicotinamide crystal in an ethanol droplet and (h) Isonicotinamide crystals in nitrobenzene droplets in Teflon tubing of 500µm diameter [6]. Copyright 2012 American Chemical Society.

Fig. 3

 3 Fig. 3 The microfluidique set-up: (1) and (2) Thermostatted Syringe-pumps, (3) T-junction, (4) Thermostatted sample-holder, (5) Camera + Zoom on XYZ stage.

Fig. 4

 4 Fig.4 Metastable zone[5] (---) and solubility curve[44] (----) of tetragonal lysozyme (0.7M NaCl, pH=4.5), lines improve legibility. Points (x) and (O) indicate nucleation and growth conditions, respectively. Copyright 2011 American Chemical Society.

Fig. 7

 7 Fig. 7 Photos of crystals of rasburicase obtained in 65nL droplets: each line corresponds to the same crystallization conditions (a) 10% PEG, 10µg/µL rasburicase, at 5°C. (b) 10% PEG, 10µg/µL rasburicase, at 20°C. (c) 5% PEG, 10µg/µl rasburicase, at 5°C. (d) 7.5% PEG, 5µg/µL rasburicase, at 5°C. (e) 7.5% de PEG, 5µg/µL rasburicase, at 20°C [9]. Copyright 2015 American Chemical Society (paper submitted).

Fig. 6

 6 Fig.6and 7 highlight the importance of performing statistical studies through a large number of experiments per condition because of the stochasticity of nucleation. Statistically each crystallization condition should be tested from 50 to 100 times in order to decide whether the result (obtaining crystals, targeted polymorph(s)) is positive or negative. This is the main problem with the conventional crystallization robots, where typically a few hundred different conditions are tested only 1 to 5 times each, leading often to non-reproducible results.

Fig. 8

 8 Fig. 8 (a) The micro-injector in oil set-up: (1) injector, (2) glass micropipette, (3) XYZ miniature translation stages for injector and (4) microscope; (b) Array of monodisperse water droplets with volumes ranging from nanoliter to femtoliter [10]. (a) Copyright 2010 Elsevier and (b) with permission from [10]. Copyright 2011 AIP Publishing LLC

Fig. 9

 9 Fig. 9 (a) Array of droplets (~300pL) of NaCl solutions at =0.1 generated through a layer of inert DMS oil. (b) Contraction of droplets due to evaporation of aqueous phase in oil. (c) Nucleation of single crystals in every droplet [16].

Fig. 11

 11 Fig. 11 Critical supersaturation (below which no nucleation can occur) vs droplet size, the red curve for a single nucleation event and the blue curve for a second event: (a) for lysozyme in NaCl solution [56]; (b) for NaCl [13]. (a) Copyright 2009 American Chemical Society and (b) Copyright 2010 Elsevier

Fig. 12

 12 Fig. 12 Time sequence showing nucleation induced by a sharp tip touching a 95pL droplet of NaCl (β>1.24) and growth of NaCl in paraffin oil: (a) to (b) nucleation in less than 5ms,(c) to (i) growth[START_REF] Hammadi | Faraday Discussions[END_REF]. Reproduced from Ref[START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry

Fig. 13 (

 13 Fig.13(a) Temperature-controlled crystallization cell with 2 micromanipulators and (b) SEM image of a W-tip[START_REF] Candoni | [END_REF]. Reproduced from Ref[START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry

Fig. 15

 15 Fig.15 Evolution of current over time during the experiment, semi-log scale. Line is a guide for the eye. Reproduced from Ref[START_REF] Hammadi | Faraday Discussions[END_REF] with permission from The Royal Society of Chemistry

  Effective interfacial energy  ef of the critical nucleus determined for oils of different chemical natures, by different authors. The deviation in absolute is given in parentheses; errors are not available in 21, 22 and 46.

	device	PDMS	Teflon wells	Teflon	PDMS
	oil	silicone	silicone	fluorinated	fluorinated
	 ef (mJ/m 2 )	0.62 (0.13)	0.56	0.88 (0.05)	0.91
	Tab. 1