Oriol Pont
email: oriol.pont@inria.fr

Antonio Turiel

Hussein Yahia

An Optimized Algorithm for the Evaluation of Local Singularity Exponents in Digital Signals

Recent works show that the determination of singularity exponents in images can be useful to assess their information content, and in some cases they can cast additional information about underlying physical processes. However, the concept of singularity exponent is associated to differential calculus and thus cannot be easily translated to a digital context, even using wavelets. In this work we show that a recently patented algorithm allows obtaining precise, meaningful values of singularity exponents at every point in the image by the use of a discretized combinatorial mask, which is an extension of a particular wavelet basis. This mask is defined under the hypothesis that singularity exponents are a measure not only of the degree of regularity of the image, but also of the reconstructibility of a signal from their points.

Introduction

Since the introduction of wavelet theory, it has been recognized that the calculation of local singularity exponents from digital signals can be used to codify them in a more compact way [START_REF] Mallat | Wavelet transform maxima and multiscale edges[END_REF][START_REF] Mallat | Singularity detection and processing with wavelets[END_REF]. The early studies carried out over turbulent flows and other systems proved that the singularity exponents at the top points in a Wavelet Transform Modulus Maxima (WTMM) line can be easily calculated [START_REF] Mallat | Characterization of signals from multiscale edges[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Muzy | Wavelets and multifractal formalism for singular signals: Application to turbulence data[END_REF]. However, the extension of this methodology to any point at resolution scale was far from simple, and for those points the WTMM method become convoluted and rather imprecise [START_REF] Struzik | Determining local singularity strengths and their spectra with the wavelet transform[END_REF], even for the mere assessment of the statistical properties of the signal [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF].

A different approach using numerical determinations of the local gradient modulus convolved with wavelets, even with positive multiscaling bases, showed more stable results over discretized signals [START_REF] Turiel | The self-similarity properties of natural images resemble those of turbulent flows[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF], leading to accurate exponents values and fine spatial resolution [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF]. This new approach to singularity analysis is in the basis of the so-called Microcanonical Multiscale Formalism (MMF) [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF], which has been shown to be useful to assess physical properties of turbulent flows and other multiscale systems [START_REF] Turiel | Multifractal method for the instantaneous evaluation of the stream function in geophysical flows[END_REF][START_REF] Isern-Fontanet | Microcanonical multifractal formalism: application to the estimation of ocean surface velocities[END_REF][START_REF] Turiel | Tracking oceanic currents by singularity analysis of micro-wave sea surface temperature images[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF][START_REF] Turiel | Role of multifractal sources in the analysis of stock market time series[END_REF][START_REF] Pont | Application of the microcanonical multifractal formalism to monofractal systems[END_REF][START_REF] Pont | Empirical evidences of a common multifractal signature in economic, biological and physical systems[END_REF][START_REF] Pont | Description, modeling and forecasting of data with optimal wavelets[END_REF]. It has been hence demonstrated that obtaining the singularity exponents at each point of a signal can reveal many useful information not only about the image, but also about the physical processes giving raise to it.

However, the performance of MMF-based singularity analysis (namely, the quality of the exponents calculated with this approach) depends on the properties of the multiscaling function used [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF], the best multiscaling functions being positive functions with fast enough decay [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF]. Although for statistical analysis almost any scaling function gives the same results [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF], a precise geometrical determination of the underlying patterns, edges and textures requires a very fine tuning of the multiscaling function.

In this paper we discuss a method for singularity analysis introduced in a recent patent [START_REF] Turiel | Method and system for the singularity analysis of digital signals, patent registered under number[END_REF]. This method allows high-performance determination of the singularity exponents, with a method which is fast, computationally cheap, stable, accurate and provides fine resolution. Its definition relies in the connection of singularity exponents with the concept of image reconstruction from the Most Singular Component (MSC), as presented in [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF]. The paper is structured as follows: in the next Section we introduce the basics of singularity analysis, while Section 3 explains the key concept of reconstruction from the MSC. In Section 4 the conditions to define a UPM-based measure are discussed, while Section 5 gives the settings for the calculus on reduced neighborhoods. Finally, in Section 6 our method is presented and some results shown. The last Section, Section 7, presents the conclusions of our work.

Definition of Singularity Analysis

Singularity analysis is a term referring to different meanings in mathematical analysis (e.g. the studies of singularities of differentiable functions); in the present work we focus on its meaning in the theory of complex systems. Using singularity analysis we intend to describe and characterize the local behavior of a R m -valued function f (x) defined on R d around each one of its domain points x according to the so-called singularity exponent, Hölder exponent [START_REF] Jaffard | Multifractal formalism for functions. I. Results valid for all functions[END_REF] or Hurst exponent [START_REF] Simonsen | Determination of the hurst exponent by use of wavelet transforms[END_REF][START_REF] Jones | Wavelet packet computation of the hurst exponent[END_REF]. If the signal behaves at point x according to the following limiting behavior:

f (x + r) -f (x) = α(x)r h(x) + o(r h(x) ) (r → 0) (1)
then h(x) is the singularity exponent at x: small displacements around x lead to function increment which scale as powers of the displacement modulus r = r . A strictly n-times derivable function obviously leads to a Hölder exponent h(x) = n, and so this formulation allows to generalize the concept of integer differentiability to real differentiability. To complete the transposition, a slightly more exigent formulation of eq. ( 1) is required, namely we should assume that there exists a (1, 1) continuous tensor from R d to R m , α(x), such that

f (x + r) -f (x) = α(x)|r r h(x)-1 + o(r h(x) ) (r → 0) (2)
(denoting α(x)|r the standard duality bracket for (1, 1) tensors). When such representation is possible, the exponent h(x) is called the Hurst exponent of the function; if only eq. ( 1) can be applied we will prefer to speak about Hölder exponents.

Assessment based on Hölder exponents can only be applied to very specific signals; in general, the presence of long range correlations and the effects of noise and discretization would preclude a direct evaluation of the scaling exponent [START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF].

A more general framework is given by singularity exponents, which are defined using gradient-based measures [START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF]. Given a signal s(x), defined in R d and with values in R, we can define its associated gradient measure μ by its density dμ(x), which is given by:

dμ(x) = ∇s (x) dx ( 3
)
This measure is by definition absolutely continuous with respect to Lebesgue measure. Hence, the measure of any Borelian A is given by:

μ(A) = A dx ∇s (x) ( 4 )
Gradient measures also allow to characterize the local singularity of any point, and in a direction-independent manner. Following eq. ( 1), let us consider a function f (x) with a Hölder exponent h(x) + 1 at a point x (notice the shift +1 introduced for later convenience). Let B r (x) be the ball (using an arbitrary norm in R d ) of radius r centered around x. So, we obtain [START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF]:

μ (B r (x)) = α(x) r d+h(x) + o r d+h(x) (r → 0) ( 5
)
where d is the dimension of the domain space (d = 2 in images). The introduction of gradient measures is convenient, as measures can also be wavelet-projected to obtain smooth interpolations from discretized data. Given a wavelet Ψ , we define [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] the wavelet projection of the measure μ at the point x and scale r, denoted by T Ψ μ(x, r), as:

T Ψ μ(x, r) = R d dμ(x ) 1 r d Ψ x -x r (6) i.e. operator T Ψ is a map from the set M of σ-finite measures on R d to the set of functions R d × R + → R.
If the signal possesses a singularity exponent at the point x according to eq. ( 5), then the wavelet projections allow to infer this same exponent, as they verify [START_REF] Daubechies | Ten lectures on wavelets[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF]:

T Ψ μ(x, r) = α Ψ (x) r h(x) + o r h(x) (r → 0) (7)
The main advantage of using measures over discretized data is that any function Ψ can be used to evaluate singularity exponents using eq. ( 7), even positive functions [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF]. As discussed in [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF], the resolution capability of a wavelet depends on the number of zero-crossings it has, which is increased in higher-order wavelets but is minimum for positive wavelets. So, gradient measures improve the spatial resolution of singularity exponents [START_REF] Turiel | Relevance of multifractal textures in static images[END_REF].

MSC and Its Connection with Reconstruction

An important ingredient in the construction of wavelets with optimized resolution capability is the concept of reconstruction of signals from partial information about its gradient. The theoretical and practical implementation of this reconstruction algorithm was first introduced in [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] (see discussion there).

We consider signals having a singularity exponent h(x) at each point [START_REF] Parisi | On the singularity structure of fully developed turbulence[END_REF][START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF] and for which these exponents are organized forming sets with a particular multiscale structure [START_REF] Falconer | Fractal Geometry: Mathematical Foundations and Applications[END_REF]. That is, the values of the singularity exponents do not take arbitrary values but must be organized so that they define a hierarchy of multiscale geometrical structures "matching" and realizing closely the cascading properties of some random variables associated to the macroscopic description of the system under study [START_REF] She | Universal scaling laws in fully developed turbulence[END_REF]. Due to the difficulties of classical methodologies to assign a precise value of singularity exponent h(x) to each point, all the characterizations of this hierarchy that have been tried up to now are merely statistical. In [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] a new question was posed: if the hierarchy truly exists in complex signals, can they be reconstructed starting from the vertex of this hierarchy? For multifractals, the set associated to the vertex is well-known, at least from the theoretical point of view: it is the so-called Most Singular Component (MSC), which is the set comprising the points with most singular (i.e., most negative) values of h(x) [START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF][START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF].

The thesis in [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] is that the MSC contains enough information to fully reconstruct the signal (in that reference, the reconstruction of images is analyzed, although the formulas are valid for any number of dimensions). As we are working with gradient measures, the data to be retained at the MSC is the gradient of the signal. So, it was hypothesized that there exists an universal operator to reconstruct signals starting from the values of the signal, and leading to a reconstruction algorithm consistent with the known statistical invariances of turbulence and multiscale signals [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF].The algorithm was required to be deterministic, linear, translational invariant, isotropic and leading to the known power-spectrum shape. Under these requirements, it turned out that there exists, if any, only one possible operator to reconstruct signals from the gradient on the MSC. Let us first define a convenient notation for the starting data. For a given multiscale signal s let us denote by

F ∞ the MSC, that is to say F ∞ is the set of points x such that h(x) ∈]h ∞ -Δ, h ∞ + Δ[ with h ∞
being the minimum value of all h(x) over the finite domain of the discrete signal, and Δ a threshold parameter; we will define the essential gradient of s, ∇ F ∞ s, as follows:

∇ F ∞ s(x) = ∇s(x) δ F ∞ (x) ( 8 )
where δ F ∞ is a delta distribution associated to the continuum of the F ∞ , homogeneous in (Hausdorff) topological dimension to a repartition in between dimensions d -1 and d: it assigns uniform weight to the points on the MSC F ∞ and vanishes outside the MSC. According to this notation, the reconstruction formula [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] reads:

s(x) = (g • ∇ F ∞ s)(x) ( 9
)
where the symbol • means convolution dot-product of vectors and the vector field g is the universal reconstruction kernel, which can be easily expressed in Fourier space, namely:

ĝ(k) = i k k 2 (10)
and i = √ -1 is the imaginary unit. So defined, the reconstruction kernel g is a kind of inverse gradient operator. There is always a set F ∞ from which reconstruction is perfect, the whole domain: If F ∞ = R d , then eq. ( 9) reduces to the trivial identity ∇ F ∞ s = ∇s. But from eq. ( 9) is not evident if there exists a smaller set F ∞ ⊂ R d such that reconstruction is also perfect. Following the derivation in [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF], we can conclude that any set F leading to a perfect reconstruction must verify:

div ∇ F c s = 0 (11)
where F c is the complementary set of F . As the divergence operator is local and the formula above is linear, the decision to include or not a point can be taken on the basis of any neighborhood around that point. The points that must always be included to obtain a perfect reconstruction are hence those with values that cannot be predicted just knowing the values in their surroundings; they are hence called unpredictable points (in opposition to the other points, which are predictable). Predictability is a subject at the core of the analysis of complex systems and signals [START_REF] Boffetta | Predictability: a way to characterize complexity[END_REF][START_REF] Aurell | Predictability in the large: an extension of the concept of lyapunov exponent[END_REF], where for instance Lyapunov exponents and Kolmogorov-Sinai entropy are known measures of information growth in a dynamical system. The set F u formed by the collection of all the unpredictable points is what we will call the Unpredictable Points Manifold (UPM) and it is, by definition, the smallest set for which eq. ( 9) lead to a perfect reconstruction. The hypothesis in [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF] is that F u = F ∞ , a conjecture which is at the base of the framework of reconstructible systems. What is evident from experiences is that the MSC leads to good reconstructions (see discussion in [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. Part I: Singularity analysis[END_REF]).

General Conditions to Define UPM-Measures

We now step forward to generalize the concept of gradient measure introduced in the previous sections to the novel concept of UPM-measure. The basic requirements to define a singular positive UPM-measure μ are:

i) It is concerned with the local singular behavior of functions.

ii) It leads to a MSC as close to the UPM as possible.

In some sense, UPM-measures are gradient measures which also take into account the degree of predictability of points according to eq. [START_REF] Mallat | Characterization of signals from multiscale edges[END_REF]. So that, they take the form of a gradient (in the sense of finite difference over discretized signals) but with penalty terms associated to the lack of predictability. The best way to keep on working around singularities is to define UPM-measures as vectorial wavelet projections of standard gradient measures. So, the UPM-measure is a carefully designed vectorial wavelet projection of the gradient measure so that it penalizes unpredictability.

In our method, in contrast with standard singularity analysis, we will not perform many wavelet projections of the UPM measure in order to extract the singularity exponents by means of a log-log regression applied to eq. ( 7). Waveletprojecting the measure at several scales is costly in computer time and only serves to enhance the resolution of less singular structures at the cost of coarsening most singular ones (see a discussion on this in [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF]). But as we are mainly interested in the most singular structures, it is hence harmful to our interests to project across multiple scales. Instead, we will make use of point estimates [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: a comparative study[END_REF][START_REF] Pont | Application of the microcanonical multifractal formalism to monofractal systems[END_REF] of the singularity exponents, namely:

h(x) = log (T Ψ μ(x, r 0 )/ T Ψ μ(•, r 0 ) ) log r 0 + o 1 log r 0 (12)
where T Ψ μ(•, r 0 ) is the average value of the wavelet projection over the whole signal and serves to diminish the relative amplitude of the o 1 log r0 correction. When applying eq. ( 12) we will need that r 0 is small enough to neglect this correction. The scale r 0 will be defined as the smallest accessible one, that is, the pixel scale. We conventionally assign a Lebesgue measure of 1 to the whole space domain, so for a N × M image the value of r 0 is fixed to r 0 = 1 √ NM , so in general we need that images are large enough to make the first term in the right hand side of eq. ( 12) a good approximation of the singularity exponent. In practical terms, this implies a resolution around 100 × 100 pixels or larger.

Calculus on Reduced Neighborhoods: Cross Fourier Transform

In order to assess the degree of predictability of a given point, we will apply the reconstruction formula, eq. ( 9), for the smallest possible neighbor of a point, namely its 2d nearest neighbors in 2d connexity neighborhoods. In 2D (d = 2) this consists of 4 neighbors, that with the point altogether form a cross. For any quantity p(x) we will represent the neighborhood of any point x 0 by a 5component vector comprising this point and its 4 nearest neighbors, following the indexing convention established in Figure 1. So, the central point will be assigned the index 0, the point at its right will be indexed 1, the one on the left is indexed 2, that on top is indexed as 3 and the one on the bottom is indexed 4. So, we convert the neighbor in the vector (p 0 , p 1 , p 2 , p 3 , p 4 ). The notion of predictable point easily extends to any number of dimensions, regardless of the number of components of the neighbor vector, which grows as the dimension d increases.

We could apply the harmonics of the standard Fourier Transform on the discrete signal (p 0 , p 1 , p 2 , p 3 , p 4 ), but this is not a good idea. The harmonics of the standard Fourier Transform (i.e., (e 2ikπ/n ) k ) depend on the size of the embedding space, so they would lead to a dimension-dependent measure of the predictability . To overcome this problem, we note that relative to the center of the cross, the position of the other points correspond to displacements of ±1 (in pixel units) either in the x-direction or in the y-direction. So that, to define a special type of Fourier transform specialized to this formation, the basic Nyquist frequency in each direction is 2π/3. Consequently we introduce j = e 2πi/3 = cos(2π/3) + i sin(2π/3) = -

1 2 + i √ 3 2 , j = j 2 (13)
We define the direct Cross Fourier Transform of any 5-vector p = (p 0 , p 1 , p 2 , p 3 , p 4 ) as the complex 5-vector p = (p 0 , p1 , p2 , p3 , p4 ) obtained according to the following formula:

p = F p ( 14
)
where F is the following 5 × 5 complex matrix:

F = 1 3 ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 1 1 1 1 1 j j 1 1 1 j j 1 1 1 1 1 j j 1 1 1 j j ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (15)
This matrix represents the linear combination of the harmonics associated to the displacements in the cross and is designed to represent with the maximum fidelity the composition at the center of the cross, starting from the nearest points. The inverse of this matrix is:

F -1 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ -1 1 1 1 1 1 j j 0 0 1 j j 0 0 1 0 0 j j 1 0 0 j j ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (16)
We need to define surrogates of the gradient and the reconstruction formula restricted to the cross neighborhood, in order to evaluate in a fast way the degree of predictability of the central point. For that reason, we will construct appropriate implementations of the gradient and of the gradient reconstruction formula, based on the Cross Fourier Transform.

The Cross Gradient Operator is the operator (∂ x , ∂ y ) = F -1 • ( ∂x , ∂y ) • F. In Fourier space the operator acts by simply multiplying any function by the functions ∂x and ∂y to obtain the x and the y coordinate, respectively. The function ∂x is defined as:

∂x = (0, i √ 3, -i √ 3, 0, 0) (17)
and analogously we have:

∂y = (0, 0, 0, i √ 3, -i √ 3) (18)
The Cross Reconstruction Operator is one of the inverses of the Cross Gradient Operator. As the gradient operator eliminates any constant summed up to each component of the 5-vector representing the neighborhood, the reconstruction is defined up to a constant shift; our implementation of the cross reconstruction operator is such that the 5-vector has zero mean, 5 i=1 s i = 0. For that reason, signals should have the mean subtracted before applying these two operators (see below).

The Cross Reconstruction is the operator R = F -1 • R • F. In Fourier space R has two functional components, R = ( Rx , Ry ); the operator acts as the sum of the product of each component with the corresponding component (x and y) of the gradient on which it is operated. The component Rx is defined as:

Rx = (0, -i/ √ 3, i/ √ 3, 0, 0) (19)
and analogously for Ry ,

Ry = (0, 0, 0, -i/ √ 3, i/ √ 3) (20)
The Cross Gradient and the Cross Reconstruction are the two basic algorithms for the design of the UPM-measures. They can be simplified to a 5 × 5 matricial form, for faster numeric implementation.

Local Correlation Singularity Measure

The Local Correlation Singularity Measure is designed to measure the unpredictability of a given point, just quantifying the difference on the actual value of the detrended (i.e., after subtracting the mean) signal at a given point and the inferred one from their four neighbors. It is defined algorithmically as follows:

Goal: To evaluate T Ψ lcsm μ(x 0 , r 0 ) at a given point x 0 . Algorithm 1. The neighborhood of x 0 is converted into a 5-vector s = (s 0 , s 1 , s 2 , s 3 , s 4 ) according to the scheme in Figure 1. 2. The vector is conveniently detrended: we first obtain S = 1 3 5 i=1 s i , and we define the detrended vector, p = (p 0 , p 1 , p 2 , p 3 , p 4 ) as:

p 0 = s 0 + S ; p i = s 0 -S , i = 1, . . . , 4
3. We apply the Cross Gradient Operator to p, so we obtain the vector g x and g y . 4. We keep the value of the first components of these two vectors for a later use, A x = g x,0 , A y = g y,0 . 5. We set these two components to zero, g x,0 = g y,0 = 0. 6. We apply the Cross Reconstruction Operator to the resulting vectors g x and g y , to obtain the reconstructed signal r. 7. We apply once more the Cross Gradient Operator onto r to obtain ρ x and ρ y . 8. We define the Local Correlation Singularity Measure as the modulus of the difference of the cross gradients at the center of the cross, namely:

T Ψ lcsm μ(x 0 , r 0 ) = (A x -ρ x,0 ) 2 + (A y -ρ y,0 ) 2
In fact, this last step means to keep the modulus of a vector-valued wavelet projection, but to simplify notation we leave it as is. 9. The singularity exponent h(x 0 ) is then obtained in application of eq. ( 12).

In Figure 2 we show an example of the application of the singularity analysis based on the Local Correlation Singularity Measure.

Conclusions

The accurate estimation of singularity exponents in multiscale systems allows characterizing their relevant features and identifying their information content. This is particularly important for the case of digital images, where the degree of singularity is directly related to the distribution of information, and so its knowledge can be used for compact coding or reconstructing from the Most Singular Component (MSC). However, digital images are discretized and this fact is an important obstacle for precisely retrieving its singularity exponents: even when using wavelet projections, most standard wavelet bases only give average results.

In this article, we have presented a recent algorithm that allows obtaining the singularity exponents of an image at every point. The singularity exponents are extracted in a precise and meaningful way, by means of a discretized combinatorial mask. This mask is constructed by considering the singularity exponent of a given point as both a measure of the singularity/regularity degree and a measure of the unpredictability of that point. The result is a discretized, numerical extension of a particular wavelet basis.

We have presented and discussed the method for singularity analysis noted as "Local Correlation Singularity Measure" in patent [START_REF] Turiel | Method and system for the singularity analysis of digital signals, patent registered under number[END_REF]. This method attains at the same time good quality and spatial resolution in the estimation of singularity exponents. Additionally, the reconstruction from the MSC is of high quality. As an illustration, we have shown the singularity exponents from the Local Correlation Singularity Measure and the reconstruction from their MSC for Lena's image, for which the reconstruction quality is of 24.5 dB. The prospects of this method includes image compression [START_REF] Turiel | Reconstructing images from their most singular fractal manifold[END_REF], assessment of streamlines in turbulent flows [START_REF] Turiel | Tracking oceanic currents by singularity analysis of micro-wave sea surface temperature images[END_REF][START_REF] Turiel | The multifractal structure of satellite temperature images can be used to obtain global maps of ocean currents[END_REF], detection of convection meteorological systems [START_REF] Turiel | Multiscale techniques for the detection of precipitation using thermal IR satellite images[END_REF] or detection of investment cycles in stock market series [START_REF] Turiel | Role of multifractal sources in the analysis of stock market time series[END_REF], among others.

The presented methodology for 2D images can be easily generalized to any dimensionality. In addition, it is possible to define other UPM-measures other than the Local Correlation Singularity Measure, that can give better performance in certain cases. All these additional developments and their respective applications will be the object of future communications.

Fig. 1 .

 1 Fig. 1. Schematic representation of the indexing of the points in the 2D cross

Fig. 2 .

 2 Fig. 2. Top:Left: Original Lena image; Right: Singularity exponents estimated using the Local Correlation Singularity Measure; they are represented using a inverse grayscale palette (the brightest the smaller, so more singular). Bottom: Left: MSC, defined as {h < -0.5}; it comprises 30% of the points of the image ; Right: Reconstruction. Some details are missing due to the lack of capability of the method to capture every UPM point; however, the reconstruction is of high quality (24.5 dB).


and research institutions in France or abroad, or from public or private research centers.

A geometric convergence theory of the preconditioned steepest descent iteration

A geometric convergence theory of the preconditioned steepest descent iteration List of authors: Klaus Neymeyr [START_REF] Neymeyr | A geometric convergence of the preconditioned steepest descent iteration[END_REF] The topic of this talk is a new convergence analysis for the preconditioned steepest descent (PSD) iteration to solve (elliptic operator) eigenvalue problems. The PSD iteration combines the Rayleigh-Ritz procedure and the preconditioned gradient iteration for optimal convergence acceleration. Sharp non-asymptotic convergence estimates are derived by means of a geometric theory.

Up to now sharp convergence estimates are only known for the basic fixed-step size preconditioned gradient iteration (also called preconditioned inverse iteration). The new result substantiates that preconditioned eigensolvers are robust and stable solvers which allow to compute (e.g. for partial differential operator eigenproblems and with multigrid preconditioning) a fixed number of the smallest eigenvalues and the corresponding eigenfunctions with costs that increase only linearly in the number of unknowns.

Institut für Mathematik, Universität Rostock, Ulmenstraße 69, D-18057 Rostock.


J Marín

A Block ILU Preconditioner for Computational Electromagnetics Applications

R. Bru 2 N. Malla 3 In this work we consider the solution of the linear systems arising from electromagentism applications by preconditioned Krylov subspace methods [START_REF] Saad | Methods for Sparse Linear Systems[END_REF]. Simulation of electromagnetic wave propagation phenomena requires the numerical solution of Maxwell's equations which it is often done by means of integral equation methods. The discretization of the integral equations with the boundary element method results in dense linear systems of the form,

Ax = b (1)
where A has complex entries and it is quite challenging to solve. Sparse approximate inverse preconditioners based on Frobenius norm minimization [START_REF] Grote | Parallel preconditioning with sparse approximate inverses[END_REF] are quite effective for this type of problems. In [START_REF] Carpentieri | Sparse preconditioners for dense linear systems, from electromagnetics applications[END_REF] a number of preconditioners are compared and it is shown that the SPAI preconditioner performs the best. It is also shown that ILU-type preconditioners fail to produce good convergence results. In this work we experiment with ILU-type preconditioners. Following the work presented in [START_REF] Benzi | A two-level parallel preconditioner based on sparse approximate inverses[END_REF], we use graph partioning techniques in order to split the adjacency graph of a sparsified matrix  in several parts of almost equal size and with the number of edge cuts as small as possible. The matrix  is obtained from A after discarding small entries which correspond to far field interactions. Nodes which are connected with cut edges are removed from the subgraphs and put in a separator set. By numbering the nodes by subgraphs and taking the separator set last, the matrix is permuted into block form,

P T ÂP =        Â1 B1 Â2 B2 . . . . . . Âp Bp Ĉ1 Ĉ2 . . . Ĉp ÂS       
This permuted matrix can be decomposed as,

P T ÂP =        Â1 B1 Â2 B2 . . . . . . Âp Bp Ĉ1 Ĉ2 . . . Ĉp ÂS        =        L1 L2 . . . . . . Lp F1 F2 . . . Fp LS               Û1 Ê1 Û2 Ê2 . . . . . . Ûp Êp . . . ÛS       

1

  Institut de Matemàtica Multidisciplinar, Universitat Politècnica de València 2 Institut de Matemàtica Multidisciplinar, Universitat Politècnica de València 3 Institut de Matemàtica Multidisciplinar, Universitat Politècnica de València By computing approximate LU factors for the matrices Âi and Ŝ with an appropiate ILU technique, a two-level ILU precontioner for the system (1) is obtained. The results of the numerical experiments show that the block ILU precontioner proposed is competitive.

	where Âi = Li Ûi , Êi = L-1 i Bi , Fi = Ĉi	Û -1 i	and LS , ÛS are the triangular factors of the Schur
	complement matrix		p
	Ŝ = ÂS -	Ĉi	Â-1 i Bi .
			i=1
			1


D Neklyudov

V Tcheverda

I Silvestrov

A Helmholtz Iterative Solver with No Use of Finite-Difference Approximations List of authors

1. The numerical solution of the Helmholtz equation in 3D heterogeneous media is of great importance for theory and practice of seismic waves' propagation. The common approach to its solution lies in the finite-difference approximation with subsequent application of some iterative techniques for solving a system of linear algebraic equations ( [START_REF] Erlangga | Advances in iterative methods and preconditioners for the Helmholtz equation[END_REF], [START_REF] Plessix | A Helmholtz iterative solver for 3D seismic-imaging problems[END_REF]). One of the main trouble on this way is indefiniteness of the matrix which seriously hampers convergence of iterations. Therefore, the development of effective preconditioners is crucial here. In this work, we construct the right preconditioner as an inverse of the Helmholtz operator for vertically heterogeneous background with complex refractive index.

Let us consider Helmholtz equation with 3D heterogeneous wave propagation velocity c( x):

L < u >:= ∆ + ω 2 c 2 (x) u( x; x s ; ω) = f ( x s ; ω) (1)
satisfying the vanishing attenuation principle ( [START_REF] Vainberg | Principles of radiation, vanishing attenuation and limit amplitude in the general theory of partial differential equations Usp[END_REF]):

u( x; ω) = lim ε→0 u( x; ω + iε); u( x; ω + iε) → 0 f or x 2 + y 2 + z 2 → ∞
Let us introduce now another Helmholtz operator:

L 0 := ∆ + (1 + iβ) ω 2 c 2 0 (z)
with function c 0 (z) and attenuation β chosen in a specific way, and decompose the 3D heterogeneous Helmholtz operator L as follows:

L = L 0 -δL = ∆ + (1 + iβ) ω 2 c 2 0 (z) -(1 + iβ) ω 2 c 2 0 (z) - ω 2 c 2 ( x)
Using L -1 0 as a right preconditioner we come to the following equation:

(I -δLL -1 0 ) < u >= f, with u = L -1 0 < u >
3. The key point of any iterative technique is efficient implementation of the matrix-by-vector multiplication. The matrix we deal with comes from approximation of operators L -1 0 and δL. The latter one is pointwise multiplication, so let us deal with the first. In order perform Krylov's iterations one needs to compute the result of its action onto arbitrary 3D function, that is one should resolve the following Helmholtz equation:

∆ + (1 + iβ) ω 2 c 2 0 (z) u tmp = f input ( x)
After 2D Fourier transform with respect to lateral coordinates this equation is transformed into a two-parameter family of ordinary differential equations:

d 2 U (z; k x ; k y ) dz 2 + (1 + iβ) ω 2 c 2 0 (z) -k 2 x -k 2 y U (z; k x ; k y ) = F (z; k x ; k y ) (2)
To avoid artificial reflections at the bottom interface the transparent boundary condition (TBC) are introduced there. For the top one there is a choice -to consider it as a free surface or to apply there TBC as well.

In order to resolve the system (2) the 1D background velocity c 0 (z) is treated as a piecewise constant function. Then coefficients of ODE ( 2) may be considered as constant within each interval and its solution there is represented as a superposition of up-and downgoing waves:

U j (z) = B j exp(-i • K j • z) + C j exp(i • K j • z) + G j (z)
with G j being some partial solution of heterogeneous ODE. Constants B j and C j are computed from conditions at the interaces:

[U ]| z=z j = [ dU dz ]| z=z j = 0 (3)
which lead to a system of linear algebraic equations with tri-diagonal matrix. The desired solution is recovered from the computed coefficients in an evident manner. To come back to lateral coordinates, one needs to apply an inverse Fourier transform and finalize computation of the vector u tmp ( x). The last step of this matrix-by-vector product is pointwise multiplication by the given 3D function:

u output ( x) = f input ( x) -[δL]( x) • u tmp ( x) (4)
4. The proposed approach was tested on the velocity model presented in Fig. 1a (Gullfaks field, see ( [START_REF] Stovas | A sensitivity study based on 2D synthetic data from Gullfaks Field, using PP and PS time-lapse stacks for fluid-pressure discrimination[END_REF])). The wavefield was computed by the presented above iterative approach for Ricker source function with dominant frequency of 50 Hz. Computations were done for 601 temporal frequencies uniformly placed from 0.1 Hz to 90 Hz. Dependence of the number of iterations of the frequency one can see in Fig. 1b. The model is discretized with with spatial steps h x = h z = 5m, that is 6 points per dominant wavelength. Snapshot corresponding to t = 1.4s is presented in Fig. 1c. As one can see, there is no visible manifestations of numerical dispersion.

Figure 1 :

 1 Figure 1: a) Gullfaks velocity model. b) Number of iterations vs frequency c) Snapshot for t=1.4 s

1

  Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russia 2 Institute of Petroleum Geology and Geophysics SB RAS and Novosibirsk State University, Novosibirsk, Russia 3 Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russia

Acknowledgements This research is done in cooperation with Schlumberger Moscow Research and partially supported by RFBR grants 10-05-00233, 11-05-00238 and 11-05-00947. We also thank StatoilHydro for providing the velocity model.


N Vigneau

G Rabatel

P Roumet

M Ecarnot

Field hyperspectral imagery as a tool for plant monitoring: application to wheat nitrogen content

Keywords: Close-range hyperspectral images, reflectance, chemometrical model, nitrogen content, wheat

In the study presented here, we used a pushbroom CCD camera fitted on a motorised tractor rail to take hyperspectral images of wheat plots. Reflectance correction was performed using a ceramic plate as a reference. The hyperspectral images have been used to build nitrogen (N) content images for wheat plots using two Partial Least Square regression models. The first chemometrical model was calibrated on nitrogen contents (in percentage of dry matter, % DM) obtained in laboratory and flat leaf normalised spectra (5 latent variables, R² = 0.907, RMSEP = 0.326 N % DM); the second one on greenhouse pot plant leaf spectra (8 latent variables, R² = 0.919, RMSEP = 0.412 N % DM). Subsequently, they were directly applied to field plant spectra, providing images of nitrogen content. The range of nitrogen content values in these images was not satisfactory. Due to important spectral differences between greenhouse pot plant leaf spectra and field leaf spectra, a model calibrated directly on field leaves is necessary.

Introduction

Vegetation monitoring is a major issue for agriculture, crop genetical selection, etc. Nitrogen content, water status, can be obtained by destructive laboratory measurements. However, nondestructive measurement methods are needed to follow a same plant over time. One of them is the use of spectroscopic measurement. Light hitting a leaf is partially reflected depending on optical properties of leaf constituents. Studying the reflected signal can give information about leaf physiological status [START_REF] Curran | Remote sensing of foliar chemistry[END_REF]). In agronomy, vegetation reflectance spectra are largely used and make agronomical data acquisition possible at plant or leaf scale. A large amount of wavelengths generating a continuous spectrum provides us with accurate information even though partially redundant. Spectroscopy, formerly reserved to laboratory samples has been extended to field spectroscopy from the fifties [START_REF] Milton | Progress in field spectroscopy[END_REF]) and has been generalised over the last years thanks to the development of robust field spectroscopy devices [START_REF] Milton | Progress in field spectroscopy[END_REF]). Different devices can be used to collect spectra on an organ like a leaf or on a whole plant. Spectral measurements are used to obtain information about leaf area [START_REF] Baret | Monitoring wheat canopies with a high spectral resolution radiometer[END_REF]) or chlorophyll concentration [START_REF] Hansen | Reflectance measurement of canopy biomass and nitrogen status in wheat crops using normalized difference vegetation indices and partial least squares regression[END_REF]), and to assess total dry matter [START_REF] Aparicio | Relationship between growth traits and spectral vegetation indices in durum wheat[END_REF]). At another scale, many hyperspectral airborne or satellite sensors have been developed these last years (Casi, Hyperion, Hymap, Aviris…). Hyperspectral imagery can help to assess withinfield yield variability thanks to vegetation indices [START_REF] Zarco-Tejada | Temporal and spatial relationships between within-field yield variability in cotton and high-spatial hyperspectral remote sensing imagery[END_REF]) or chlorophyll content by combining an indice-based approach and a radiative transfer model inversion [START_REF] Haboudane | Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture[END_REF]). [START_REF] Inoue | Synergy of remote sensing and modeling for estimating ecophysiological processes in plant production[END_REF]) reviewed different ways of combining methods to access agronomic data from remote sensing. Due to their elevation, satellite or airborne sensors have a wide swath and are well-suited to studies on a large area. Unfortunately, due to their poor spatial resolution, they conduct to mixed pixels and to difficulties to estimate quantitatively and very accurately agronomical parameters like nitrogen content. At the contrary field spectroscopy can provide "pure" spectra, i.e. concerning vegetation only, but for a smaller area and with a lesser representativeness of the study area. Sometimes, it is required to have both pure spectra and good area representativeness. It is the case for crop genetical selection for example. Geneticists need accurate estimation of agronomical parameters for a small homogeneous area (micro-plot). A new technology can be tested for this kind of study: close-range hyperspectral sensing. Its representativeness is good because it provides an image of the whole plot. However, it provides pure spectra because spatial resolution is very high. The study presented here evaluates the potential of such a technology to obtain agronomical information on wheat. The aim of this article was to study the feasibility and the difficulties involved in using of ground hyperspectral imagery for crop characterisation. For that purpose, several chemometrical models between reflectance spectra and leaf nitrogen content were build. Subsequently, their application was tested on field images of durum wheat in order to see whether leaf nitrogen content could be predicted at several phenological stages.

Material and method

Hyperspectral image acquisition system

All hyperspectral images were acquired with a pushbroom CCD camera (HySpex VNIR 1600 -160, Norsk Elektro Optikk, Norway) fitted on a tractor-mounted motorised rail (see Figure 1). The camera spectral range was from 400 nm to 1000 nm divided in 160 bands (3.7 nm spectral resolution). The first image spatial dimension was determined by the 1600 across-track pixels of the CCD matrix and the second one came from the camera forward movement on the ramp. At 1 m above the vegetation and with a nadir sighting, the ground track was about 30 cm and the spatial resolution across track was 0.2 mm (the lens and the view angle were fixed). The spatial resolution along track was set to 0.5 mm. The integration time, i.e. the time duration during which sensor is storing light energy was fixed manually by the user depending on meteorological conditions (cloudy or sunny weather).

Reflectance correction

Images were first corrected in radiance using sensor characteristics (e.g. spectral sensitivity) provided by the manufacturer. Radiance is the product of the target reflectance, which is intrinsic information and the illumination during image acquisition, i.e. in our case solar lighting.

E R L ⋅ = (1)
Radiance can not be used directly because illumination depends on date and meteorological conditions. To obtain the variable of interest R, it is thus necessary to know the illumination. For that purpose, spectralon® (Labsphere, Inc., New Hampshire, USA) placed in the field of view of the sensor is commonly used because it is a perfect Lambertian diffuser. Another alternative is to use a reference whose spectral characteristics are known. Indeed, in given lighting conditions:

E R L et t et t ⋅ = arg arg (2) E R L ref ref ⋅ = (3) ref ref et t et t R L L R ⋅ = arg arg (4)
where R designs the reflectance, L the radiance, and ref the reference. In this study, we used a ceramic plate appropriate for an every day field use. ref R was obtained from laboratory measurements.

Experimental protocol

Flat leaves

During the 2009 growing season, many wheat leaves of four genotypes (ixos, primadur, neodur and lloyd) were cut, oven-dried and conserved in a cold room. Dried wheat leaves were then flat imaged (see Figure 2). In order to have an homogeneous background, we used the leaf-clip disc of a field spectrometer (FieldSpec®, Analytical Spectral Devices, Inc. (ASD), Boulder, Colorado, USA) designed to be totally black. Leaf nitrogen content was measured in laboratory (Dumas method). Hyperspectral images were corrected in reflectance. Then, a mean reflectance spectrum was calculated for each leaf. Thus we had for each leaf a reflectance spectrum and an actual nitrogen content value. The dataset contained 146 couples spectrum/nitrogen content. The nitrogen content range was 0.40 -3.78 N % DM.

2.3.2.Standing plants

In winter 2009 -2010, some plants were grown in greenhouse with two nitrogen treatments: with or without nitrogen supply. Four genotypes (lloyd, neodur, primadur and ixos) were imaged at five phenological stages (tillering, 2 nodes, flowering, 450 degree days after flowering and maturity) with three repetitions. For each image, two genotypes and three repetitions at one stage were imaged in order to have several pots in the same image and thus to to simulate field conditions (see Figure 3). On each plant the two upper leaves were marked with coloured plastic collar and thereby located on the images. After each image acquisition, leaves were cut and sent to the laboratory for destructive nitrogen content measurement. For each image, we drew regions of interest on marked leaves and obtained therefore their mean reflectance spectrum. Consequently, we had for each plant two couples (one for each of the two cut leaves) leaf reflectance spectrum / actual nitrogen content. The dataset contained 180 couples spectrum/nitrogen content. The nitrogen content range was 0.81 -5.52 N % DM.

Chemometrical model calibration

As we see below, several models were calibrated, using either flat leaf data (cf. § 2.3.1) or standing plant leaf data (cf. § 2.3.2). Model calibrations were computed with partial least square regression (PLS) using Matlab software (TheMathWorks, Natick, MA, USA) and our own functions. For each experiment, we split the dataset in a calibration set of two third of the samples and a test set of one third of the samples, having both the same distribution. We calibrated the model by cross-validation leave-one-out on the calibration set. The best calibration equation was selected on the basis of a large coefficient of multiple determinations (R²) and a low standard error of cross-validation (SECV). SECV is a measure of the difference between the actual and predicted values calculated over all cross-validation calibrations. Next the model was tested on the independent test set.

Results and discussion

Reflectance correction

Reflectance theoretically depends on both illumination and viewing configurations, this phenomenon being formalised by a bidirectional reflectance distribution function (BRDF). However, to a first approximation, a leaf can be considered as a Lambertian surface [START_REF] Chelle | Could plant leaves be treated as Lambertian surfaces in dense crop canopies to estimate light absorption?[END_REF]; [START_REF] Grant | Diffuse and specular characteristics of leaf reflectance[END_REF]), i.e. its BRDF is a constant. So flat leaves and inclinated leaves have the same reflectance at a nadir viewing. Nevertheless, due to their dissimilar orientation toward the sun, all leaves do not receive the same level of illumination. They do not receive either the same level than the reference ceramic. Each illumination level is linked to the cosine of the angle between the surface and the light incidence. Because this difference is independent of the wavelength, it can be introduced as a multiplicative factor. The reflectance for field leaves obtained with the correction described above (see § 2.2) is thus their actual reflectance up to a scalar factor. Such a multiplicative effect can be removed by spectrum normalisation. Therefore, in order to be applied on field leaves, model must use normalised spectra.

Flat leaf model

Model prediction quality

On Figure 4, we can see the results of calibration and test steps. The best model was obtained for normalised spectra and with 5 latent variables (LV). Model was considered as satisfactory even if error seems to slightly increase with nitrogen content.

Application to field images

The model was applied on field images. Figure 5 presents the results for the genotype primadur at the stem elongation stage. NDVI (Normalised Difference Vegetation Index [START_REF] Tucker | Red and Photographic Infrared Linear Combinations for Monitoring Vegetation[END_REF])) was computed on the hyperspectral image. A mask was created by thresholding this NDVI value and was applied on image to select vegetation pixels. Although we do not have ground truth for each pixel of the image, the model seems to overestimate the nitrogen content on field leaves, which is not in agreement with expert expectations. Moreover, in April (at stem elongation), nitrogen content obtained by other measurements (leaf-clip spectrometer for example) is around 3.5 -4 %. We can put forth some hypotheses to explain that. Firstly, the model has been calibrated on dried leaves. Normally water content does not influence reflectance spectra in the 400 -900 nm spectral range but drying leaves can generate polyphenols, disrupting model calibration. Secondly, the model was calibrated on flat leaves and did not take into account field scene effects like specular or multiple reflexions. That is why we decided to calibrate a model directly on leaves on standing plants. For practical reasons (notably growing season was ended), the new model was calibrated on leaves of plants grown in pots in a greenhouse.

3.3.Standing plant leaf model

3.3.1.Model prediction quality

On Figure 6, we can see the results of calibration and test steps. The best model was obtained for normalised spectra and with 8 latent variables. The model is not totally satisfactory because the set contains three big groups and not a continuous range. Nevertheless this figure shows that it is possible to calibrate a model directly on standing plant leaves.

3.3.2.Application to field images

In order to see if a model calibrated on leaf of standing plant (and thus taking into account field scene effects) was better than a model calibrated on flat leaves, we applied this new model on field images. Once again, nitrogen content for leaf was overestimated (see Figure 7). An in-depth study of the model coefficients showed that the model was calibrated principally on the visible domain and on the ratio between visible and near infra-red domains. Figure 8 shows that these spectral characteristics are precisely the ones where greenhouse standing plant leaves spectra and field leaves spectra differ. That is why model calibrated on standing plant leaves do not give good results with field leaves.

This difference can probably be explained biologically. Plants in greenhouse grow faster than field plants and often the plants are delicate and their leaves are thinner, that could lead to these differences in the reflectance spectra.

4.Conclusion

In this article, we studied a new kind of hyperspectral technology: a close-range field imaging system. We took images at 1 m above the canopy of wheat plots. We used a ceramic plate as a reference to correct images in reflectance, conducting to a reflectance up to a scalar factor for inclinated leaves. We calibrated chemometrical models on flat isolated leaves and on pot plant leaves. These experiments demonstrated that it is possible to calibrate a chemometrical model between reflectance spectra and nitrogen content with a close-range imagery system and that it could be applied to various genotypes. However, it has been shown that a model calibrated on flat leaves or even on standing greenhouse plant leaves can not be applied directly to field leaves. As a conclusion, close-range hyperspectral imagery still seems a suitable tool for nitrogen content monitoring of field leaves, but it requires imperatively that the model is calibrated on field leaves also. The next step will be to acquire such a calibration data set.

Figure 1 :

 1 Figure 1: Prototype developed by Cemagref.

Figure 2 :

 2 Figure 2: Experimental protocol: flat leaves

Figure

  Figure 3: Experimental protocol: standing plants

Figure 4 :

 4 Figure 4: Prediction quality of flat leaf chemometrical model for calibration set (blue stars) and test set (red circle). SEP means standard error of prediction.

Figure 5 :

 5 Figure 5: a. Colour image of plot at stem elongation stage (ceramic reference is visible on the left), b. Cartography of predicted nitrogen content obtained with the model calibrated on flat leaves.

Figure 6 :

 6 Figure 6: Prediction quality of standing leaf chemometrical model for calibration set (blue stars) and test set (red circle). SEP means standard error of prediction.

Figure 7 :

 7 Figure 7: a. Colour image of plot at stem elongation stage (ceramic reference is visible on the left), b. Cartography of predicted nitrogen content obtained with the model calibrated on standing plants.

Figure 8 :

 8 Figure 8: Normalised spectra comparison between field leaf spectrum (blue dash) and greenhouse leaf spectrum (solid red).


Juan Undurraga
email: jundurra@clinic.ub.es

PhD Md

Inmaculada Baeza

Marc Valentí

ML Lázaro

Sir:

The association between brain lesions and psychiatric disorders has long been known. A wide variety of these syndromes have been reported since Dr John Harlow's description of his patient Phineas Gage more than 150 years ago [START_REF] Harlow | Passage of an iron rod through the head[END_REF]. These syndromes have included personality changes, psychosis, obsessive-compulsive symptoms, depression, anxiety disorders, mania, cognitive deterioration and anorexia nervosa [START_REF] Bunevicius | Brain lesions manifesting as psychiatric disorders: eight cases[END_REF][START_REF] Madhusoodanan | Brain tumor and psychiatric manifestations: a case report and brief review[END_REF][START_REF] Madhusoodanan | Psychiatric manifestations of brain tumors: diagnostic implications[END_REF]. Most reports describe tumors of the frontal or temporal lobes, and very few discuss the psychiatric sequelae of tumors primarily affecting subcortical structures [START_REF] Martínez-Lage | Thalamic tumors: clinical presentation[END_REF][START_REF] Rushing | Germinoma: unusual imaging and pathological characteristics. Report of two cases[END_REF]. Moreover, there have been few reports in the pediatric literature concerning this phenomenon [START_REF] Craven | Pineal germinoma and psychosis[END_REF][START_REF] Gamazo-Garrán | Obsessive-compulsive disorder secondary to brain dysgerminoma in an adolescent boy: a positron emission tomography case report[END_REF][START_REF] Mordecai | Case study: suprasellar germinoma presenting with psychotic and obsessive-compulsive symptoms[END_REF].

We report on a well-adjusted 17-year-old adolescent with no prior psychiatric history who was brought by his parents to our psychiatric emergency room. He presented with a 3 month history of progressive social withdrawal, apathy, lack of drive and deterioration of academic achievement.

During the month prior to consultation, he also showed neglect in personal cleanliness and daytime somnolence. There was no significant medical or developmental history, and he was receiving no medication. No family history of mental illness was reported.

At clinical admission, the patient was conscious, attentive and well-oriented in terms of place and person but with respect to time he could only identify the year. He had short term memory loss, was perplexed, exhibited psychomotor inhibition and maintained poor eye contact. His affect was flattened and incongruent in quality, with unprovoked smiles. He was apathic and presented loss of drive. He showed poverty of content of speech, and spoke in an unfluent, low and unexpressive tone. Thought processes were retarded and overly concrete. He was not delusional, but described hearing noises and voices calling his name and sporadic visual sensations on the edge of his visual field which he described as "wing flutters". There were no suicidal thoughts or plans. He had daytime hypersomnolence. He demonstrated little insight into his condition but accepted hospitalization with indifference. The Positive and Negative Syndrome Scale (PANSS) was administered and the scores were 13 on PANSS positive, 36 on PANSS negative, 33 on PANSS general and 82 on PANSS total. Physical and neurological examination revealed a limitation in supraversion of the eyes, with no other significant abnormalities. Serum chemistry and complete blood count was normal. A brain magnetic resonance image (MRI) revealed a heterogeneous cerebral mass, with necrotic and hemorrhagic areas, infiltrating the thalamus, mesencephalon, right mammillary tubercle, interventricular septum, and corpus callosum with mild perilesional edema extending to both internal capsules. Secondary ventricular dilatation was also observed.

The patient was transferred to the neurosurgery unit where a cerebral biopsy and ventriculoperitoneal derivation was performed. The biopsy revealed an intracranial germinoma with syncitiotrophoblast cells. The patient had a new surgical intervention because of hydrocephalic signs and was subsequently transferred to the oncology department to receive chemotherapy. A week later, he presented a brain stem haemorrhage and died.

Discussion

This is a case of an adolescent male with negative and positive psychotic symptoms who had an intracranial germinoma. Germinomas, a type of germ cell tumor, reach a peak incidence in the second half of the second decade and affect males ten times as often as females [START_REF] Mordecai | Case study: suprasellar germinoma presenting with psychotic and obsessive-compulsive symptoms[END_REF]. Up to 35% of patients will be asymptomatic for more than six months [START_REF] Packer | Intracraneal germ cell tumors[END_REF]. In general, germinomas carry an excellent prognosis with most series suggesting five-year progression free survival rates and a cure in well over 90% of patients [START_REF] Packer | Intracraneal germ cell tumors[END_REF]. However, prognosis is highly dependent on the histological subtype of the tumor present. Syncytiotrophoblastic tumors tend to relate to less favorable outcome, although this has no been found to be the case in all series [START_REF] Packer | Intracraneal germ cell tumors[END_REF].

Our patient had prominent psychiatric symptoms that made us suspect a new-onset psychosis.

However, in the neurological examination, he showed a limitation in the supraversion of the eyes. This limitation is part of the dorsal midbrain syndrome described by Henri Parinaud in the late nineteenth century [START_REF] Pearce | Parinaud's syndrome[END_REF]. This syndrome usually results from damage to premotor pathways for binocular upward gaze in the dorsal mesencephalon [START_REF] Miller | Late improvement in upward gaze in a patient with hydrocephalus related Parinaud dorsal midbrain syndrome[END_REF]. It may be caused by extrinsic or intrinsic lesions, with hydrocephalus being a common aetiology [START_REF] Miller | Late improvement in upward gaze in a patient with hydrocephalus related Parinaud dorsal midbrain syndrome[END_REF]. Additionally, the patient could have had cerebellar pathways involvement, as suggested by the cerebellar cognitive affective syndrome [START_REF] Riva | The cerebellum contributes to higher functions during development: evidence from a series of children surgically treated for posterior fossa tumours[END_REF][START_REF] Schmahmann | The neuropsychiatry of the cerebelluminsights from the clinic[END_REF] he presented. This syndrome, described by Schamahmann and colleagues [START_REF] Schmahmann | The neuropsychiatry of the cerebelluminsights from the clinic[END_REF], supposes that the cerebellum participates in the organization of higher order function. Lesions that involve the cerebellum or its pathways can manifest as impairment of executive functions, difficulties with spatial cognition, personality change with blunting of affect or disinhibited and inappropriate behaviour and language deficits including agrammatism and dysprosodia.

Patients with brain tumors may present with a wide variety of psychiatric symptoms, and have minimal or no neurological symptoms or signs. These symptoms are very unspecific. A recent review of the literature concluded that neither tumor type nor location is associated with any particular psychiatric manifestation [START_REF] Madhusoodanan | Psychiatric manifestations of brain tumors: diagnostic implications[END_REF].

Early detection can have a direct influence on treatment options and outcome, so a high degree of awareness is recommended. We think it is also important to assess patients that present with a first episode of psychosis with neuroimaging tools to rule out the possibility of brain lesions.