Tomographic reconstruction of tokamak plasma light emission from single image using wavelet-vaguelette decomposition

Images acquired by cameras installed in tokamaks are difficult to interpret because the three-dimensional structure of the plasma is flattened in a non-trivial way. Nevertheless, taking advantage of the slow variation of the fluctuations along magnetic field lines, the optical transformation may be approximated by a generalized Abel transform, for which we propose an inversion technique based on the wavelet-vaguelette decomposition. After validation of the new method using an academic test case and numerical data obtained with the Tokam 2D code, we present an application to an experimental movie obtained in the tokamak Tore Supra. A comparison with a classical regularization technique for ill-posed inverse problems, the singular value decomposition, allows us to assess the efficiency. The superiority of the wavelet-vaguelette technique is reflected in preserving local features, such as blobs and fronts, in the denoised emissivity map.


Introduction
In the edge of tokamaks, cross-field particle and heat transport have been attributed in a large part to convection by flows, either due to sustained turbulence or on the occasion of violent events such as edge-localized modes [1]. Not only do the large particle and heat fluxes leaving the core in this manner deteriorate confinement, but the associated vessel erosion also impairs the awaited viability of long lasting discharges. It is thus of primary importance to describe and understand flows in the edge region, and in particular their spatial and temporal organization. Since the edge plasma is at a relatively low temperature, it can under certain conditions emit visible light, which may then be collected using an optical device to obtain a two-dimensional image [2]. Using a camera with a fast enough sampling rate (tens to hundreds of kHz), it is now possible to resolve in time part of the turbulent dynamics.
One of the current limitations of such optical diagnostics is that the received flux cannot be directly related to the volumic emissivity of the plasma, because the photons collected by each pixel on the camera sensor have been emitted all along a corresponding ray, rather than out of a single point in space. Nonetheless, impressive progress was achieved by qualitative analysis of sequences of images obtained from cameras as have been installed on several devices in recent years [3][4][5][6]. This has been possible because the dominant structures in tokamak edge turbulence happen to be field-aligned filaments, commonly known as 'blobs', that have a higher density than their surroundings, and whose structure varies much slower along magnetic field lines than in the orthogonal directions. Thanks to this remarkable property, it is possible to detect individual blobs on the projected movie and to analyse their behaviour. However, to describe blob trajectories more accurately, and also to obtain a picture of the plasma which goes beyond the detection of individual blobs, a way must be found to invert the optical projection effect and to reconstruct the emissivity in three dimensions. When it is acceptable to use concurrently many devices with different orientations in the experimental setup, relatively classical tomographic methods can be successfully applied to achieve this goal (see, e.g., [7] where several photomultiplier tubes are combined in this way). But in most cases, the restricted number-typically one-of available cameras imposes a completely different approach.
Fortunately, it has been shown that slow variation along magnetic field lines is not only a property of blobs, but of all the fluctuating quantities participating in the edge transport process, such as density, momentum and temperature. As we propose to demonstrate further down, this approximate symmetry of the plasma may be exploited to extract threedimensional information out of individual camera movies. A similar but geometrically simpler situation arises when taking pictures of axisymmetric objects which radiate volumically, as do some nebulae [8], or mechanical devices probed with x-rays [9]. When the object is viewed under a small enough angle, the mathematical transform which relates the emission of each point to the intensity of the image is called the Abel transform, a Volterra integral operator of the first kind [10]. In the case of the tokamak geometry, the symmetry is helical rather than axial, and the angle of view is relatively large since the camera cannot be placed much further away from the plasma than the confinement vessel. The problem is thus more challenging and, to our knowledge, has not been treated mathematically or even numerically in the past before, probably because the helical geometry is rather specific to tokamaks and was not encountered in other contexts. Note that related ideas were applied to analyse bolometer measurements from the tokamak DIII-D [11], but under the simpler assumption that the radiation coming from the core of the plasma was constant over whole magnetic surfaces rather than field lines.
Once the link between the three-dimensional radiation and the two-dimensional image is understood, the reconstruction becomes an inverse problem which has a formal solution under the assumed symmetry, but is ill-posed in the presence of noise. A suitable numerical inversion procedure needs to be devised in order to avoid amplification of experimental errors. Many classical approaches, as for example leastsquares iteration [12] or singular value decomposition (SVD) [13] are described in textbooks [14]. Although they regularize the problem by damping the modes of the inverse transform which would lead to amplification of the noise, they typically do not take advantage of the spatial localization of coherent structures present in the plasma, a shortcoming which limits their performance.
To better tackle inverse problems in the presence of localized structures, Donoho [15,16] has proposed to turn to the wavelet-vaguelette decomposition (WVD) which had been introduced by Tchamitchian [17]. The rationale behind this choice is that wavelet decompositions keep track of both scale and space localization, as illustrated, for example, by the extraction of coherent bursts in turbulent edge plasma [18]. WVD was successfully applied to tomographic reconstruction in the context of medical imaging [19][20][21], where the relevant mathematical model is the Radon transform, and also for solving partial differential equations [22,23]. Except for the advantage of spatial localization, WVD has many similarities with the SVD. An intuitive way of understanding the difference between WVD and SVD in this context is that the former emphasizes the space-scale structure of the field to be reconstructed, while the latter aspires to an optimal representation of the optical apparatus.
The purpose of this work is to provide numerical evidence that the method can also be applied successfully to the generalized Abel transform in the context of tokamak edge turbulence imaging, in both academic and realistic situations, and to illustrate the efficiency of the new method by comparing it with a naive least-squares approach, and to a more competitive SVD method. First, we introduce a generalized definition of the Abel transform, and propose a way to solve the associated inverse problem by WVD. We then show a few academic examples of inversion, validate the procedure by applying it to independently generated images, and finally, apply it to a movie acquired in the Tore Supra tokamak.

Notations
We start with magnetic coordinates ( , θ, ϕ), where is a flux coordinate, ϕ is the toroidal angle and θ is an appropriate poloidal coordinate such that the magnetic field lines satisfy the equations where q is the safety factor. In the following, we take for granted that a relationship has been established between those magnetic coordinates ( , θ, ϕ) and the fixed cylindrical coordinates (ρ, z, ϕ), where ρ is the distance to the axis of symmetry of the tokamak, and z is the vertical coordinate with respect to the horizontal midplane. In general, such a relation will be obtained numerically using a magnetic reconstruction code. Here, for demonstration purposes, we consider simplified analytical formulae, which apply to the special case where magnetic surfaces have circular cross sections, and can then be indexed by their radius = r: Here, the parameter γ controls the distance between the centre of the magnetic surfaces and the magnetic axis, which allows modelling of the Shafranov shift effect if necessary. The axially symmetric case is approached in the limit q → +∞. The magnetic configuration is schematized in figure 1. The camera is modelled by a screen located at a distance d from a vanishingly small diaphragm C whose coordinates are denoted (ρ C , z C , ϕ C ). The rays going through C are parametrized by two angles: α is the angle with respect to the half-plane ϕ = 0 and β is the angle with respect to the plane z = 0. We call H the point on a ray which has the smallest ρ, and ρ H the distance of closest approach of the ray to the z-axis. See figure 2 for a schematic view of the setup. Elementary geometry allows us to express α and β as functions of ρ H and z H as follows: and, denoting by s the arc length along a ray, with s = 0 at H , we obtain a parametric representation of the ray in cylindrical coordinates:

Helical Abel transform
Let S 0 ( , θ, ϕ, ) be the emissivity of the plasma at a point M characterized by its field line coordinates ( , θ, ϕ) and in the direction defined by the unit vector . By definition, S 0 ( , θ, ϕ, )d 3 Md 2 is the power radiated by an infinitesimal volume of the plasma located around M, towards the directions contained in an infinitesimal cone around . From the start, we assume that the radiation is isotropic, so that S 0 does not depend on . We leave aside the effect of the spectral dependence of the camera sensitivity within the excited frequency range, which could introduce complications since it comes as an additional convolution operator on top of the geometrical scrambling. The validity of this hypothesis should be checked case by case, and we will briefly come back to it in the experimental section below. Moreover, we consider that the plasma is optically thin or transparent, so that the flux density received by the camera screen around a point P = (x, y) is given, up to a dimensional constant depending on d and the opening of the diaphragm, by the integral of the volume emissivity along the unique ray going through P and C: where are the coordinates in the image plane, and is the position of the camera on the ray. We denote by K the operator such that I 0 = KS 0 . We are interested in the action of K on special kinds of emissivity fields, namely those that vary slowly along magnetic field lines as defined by equation (1). Due to the fact that the safety factor q can be irrational, it would not be realistic to assume that S 0 is constant along field lines since by continuity that would, in most cases, imply that S 0 is constant over whole magnetic surfaces, which is not observed in practice. In the following, we shall thus assume that S 0 is constant on any connected portion of a field line visible in the camera field. In practice, the fluctuations along the magnetic field lines should satisfy k 1 R for this procedure to work. The restriction of K to helically symmetric fields now associates a twodimensional image I (x, y) to any given poloidal cross section of the emissivity, S 0 ( , θ ). The classical Abel transform, which we have mentioned in the introduction, is approximately recovered when q 1 and ρ C R. In practice, to compute I 0 , (2a) and (2b) need to be inverted and (5) can then readily be applied. In our case we do the inversion analytically, using equation (2a) which defines r as the only positive root of a quadratic polynomial. In the generic case, no analytical representation of the field lines is available, and this inversion will have to be done point-wise using a lookup table. To discretize the operator K, we restrict (r, θ ) to a rectangular domain [r min , r max ]×[θ min , θ max ] and use a regular Cartesian grid (r i , θ j ) with i = 1, . . . , N r and j = 1, . . . , N θ . We assume that S 0 vanishes for r > r max and for r < r min , and is periodic with the period θ max − θ min in the θ direction. Similarly, (x, y) is restricted to [x min , x max ] × [y min , y max ] and discretized using a grid (x i , y j ) with i = 1, . . . , N x and j = 1, . . . , N y . The integral in (5) is approximated using a quadrature rule, the method of rectangles with 1024 points on the interval [−s max , s max ], where s max is defined by in order to cover the whole interval where the ray intersects the smallest vertical cylinder containing the plasma. From now on, we consider that the operator K has been discretized in this manner, but we keep the same notation for simplicity. K thus transforms arrays defined on the poloidal grid (r i , θ j ) into arrays defined on the image grid (x i , y j ).

Wavelet-vaguelette decomposition
To reconstruct S 0 from the observation I 0 , K needs to be inverted. The restriction of (r, θ ) to a rectangular domain, which we have imposed for numerical reasons, implies that depending on the values of r min , r max , θ min , θ max , part of the information may not be visible to the camera, so that K is not invertible in general. A rigorous mathematical study of the invertibility of K is out of the scope of this paper, and in the following we will be content to consider the pseudoinverse. The issue which, however, cannot be left aside is that the measured I is corrupted by noise. Assuming for simplicity that this noise is additive, Gaussian and white, we thus observe in practice and regularization is required to avoid amplification of the noise and thus obtain a reasonable estimate for S 0 from the noisy data. We propose to use the WVD [15,16] method, which is now briefly recalled. In the following, the L 2 scalar products are denoted by brackets, · | · , and the associated norms by · , in both the (r, θ) and (x, y) planes. When needed, the integrals are discretized using the method of rectangles over the Cartesian grids defined above. Consider a periodized orthogonal wavelet basis (ψ λ ) λ∈ associated with a multiresolution analysis on the domain [r min , r max ] × [θ min , θ max ]. The elements of the set used to index the wavelets are of the form (j, µ, i), where j indicates the scale of the wavelet, µ its direction of oscillation, and where the two integer components of i specify the position of the wavelet (details may be found in textbooks, e.g. [24]). The corresponding vaguelette families (ξ λ ) λ∈ and (χ λ ) λ∈ are then defined by where K * is the adjoint of K, and the constants κ λ are chosen in order to impose ξ λ = 1 for all λ. Therewith χ λ and ξ λ are defined in the image plane (x, y), while ψ λ is defined in the poloidal plane (r, θ). From their definitions ((10a) and (10b)) it follows that the families (ξ λ ) and (χ λ ) satisfy the biorthogonality relations: We refer the reader to the literature for mathematical background on vaguelettes [17,25,26]. The general idea of WVD to solve the inverse problem in a stable way is to take advantage of the formal inversion formula derived from (10a) and (11), which describes the emissivity S 0 in the poloidal plane in terms of the vaguelette coefficients I | ξ λ of the corresponding image I 0 = KS 0 . To understand how this can be done, first note that according to (9) the vaguelette coefficients of the noisy signal can be put in the following form: Now, since the noise is Gaussian and white, and since ξ λ = 1, the coefficients W | ξ λ are Gaussian, identically distributed random variables. The variance of the noise is thus spread over many vaguelette coefficients. In contrast, the variance of the signal is likely to be concentrated on a few coefficients corresponding to marked features in the image, blobs for example. A natural way of removing the noise is therefore to retain only the largest among those coefficients, which is precisely the approach adopted by WVD. Indeed, the WVD-reconstructed emissivity S R is defined similarly to (12) but retaining only the vaguelette coefficients whose magnitude is larger than a certain threshold :  The value of the threshold needs to be estimated from the data, since the optimal threshold to use depends on the level of noise which is unknown a priori. In order to do this, we propose to use an iterative algorithm [27,28], which we briefly recall now. Once all vaguelettes coefficients I | ξ λ of the image have been computed, a first estimate for the threshold is defined by where c is a dimensionless constant of order unity which controls the denoising sensitivity. Note that we normally prefer the notation q, as in 'quantile', instead of c (see [28]), but in the tokamak context q could be confused with the safety factor and has to be avoided. 0 is thus proportional to the standard deviation of all wavelet coefficients, which is a classical way of eliminating outliers in a statistical sequence. Then, successive approximations of are obtained by iterating the sequence until convergence. This complicated formula simply means that is in the end proportional to the standard deviation of the coefficients contained in the interval [−c , c ]. The choice of c should be made depending on the application one has in mind for the reconstructed emissivity. If c is very large, the field will be blurred, but if c is very small, some artefacts due to the noise will persist. We have found that for qualitative image analysis, c = 4 is a good compromise. This value has been consistently used to obtain all of the results that follow.
To implement formula (13) in practice, the vaguelettes χ λ are first computed by applying K to the wavelets, and the ξ λ are then obtained by inverting the linear system (11). To do this, several methods are possible. We have found that leastsquares iteration is inefficient due to slow convergence, and have therefore chosen to apply a sparse LU decomposition to an augmented regularized system [29]. This also has the benefit of keeping memory requirements relatively low (see table 2). Overall, the algorithm thus consists of the following steps: (i) explicitly construct a sparse matrix representation of K using the discretization of equation (5), (ii) obtain the χ λ from (10a), (iii) compute the ξ λ from (11), using sparse LU decomposition and the PETSC sparse parallel matrix library [30], (iv) determine the vaguelette coefficients I | ξ λ by sparse matrix multiplication, (v) reconstruct a denoised poloidal emissivity map S R by fast inverse wavelet transform according to (13), (vi) optionally, apply K to obtain a denoised image I R = KS R .
For given plasma and camera setups, steps (i)-(iii) need to be performed only once, while steps (iv)-(vi) must be repeated for each movie frame. Note that, in practice, the linear system inversion (iii) is by far the most expensive step in this algorithm. In the following we take the most regular Coiflets [31] with filters of length 6 as the wavelet family. The wavelets in this family have two vanishing moments, which is sufficient to ensure that the associated WVD is well defined, while at the same time the shortness of the wavelet filter ensures that the matrices remain sparse enough, so that the problem is computationally tractable. Note, however, that the fast wavelet transform algorithm [24] cannot be used when computing the vaguelette coefficients in step (iv). Indeed, the operator K not being translation invariant, there is no filter bank corresponding to the associated vaguelettes.

Alternative methods
For comparison, we have also implemented two other, more classical methods of inversion. First, to illustrate the difficulty of inverting the operator K in a stable way in the presence of noise, we shall consider the least-squares solution: which is easily computed using the LSQR iterative algorithm [12]. Strong noise amplification is anticipated to corrupt S LS and prevent its use as a meaningful estimate of S 0 . A classical way of mitigating the amplification of noise is to replace K in (16) by a lower rank approximation obtained by truncating the singular value decomposition (SVD) of K. Let us briefly recall how this is done. Denote N = min (N r N θ , N x N y ). By the SVD theorem, there exist orthogonal families (u i ) 1 i N and (v i ) 1 i N respectively in the source and target spaces of the discretized K, and a sequence of positive real numbers (η i ) 1 i N such that K admits the expansion The (η i ) are unique and are called the singular values of K, while the (u i ) and (v i ) are, respectively, called left and right singular vectors of K. The properties of the SVD [13] ensure that the least-squares solution (16) is also given by the formula Note the formal similarity of (18) with the WVD-inversion formula (12). As for the WVD, the coefficients I | u i appearing in (18) can be split into a contribution coming from S 0 , and a contribution coming from the noise. A classical procedure used to reduce the effect of the noise is thus to replace (18) by a truncated reconstruction formula: where N 0 N is adjusted in order to obtain a qualitatively good reconstruction S SVD . Note that the retained SVD modes are not determined by the data like in the wavelet case, but rather by the global cut-off parameter N 0 , for which no wellestablished automatic selection procedure exists. We shall study in a specific example the effect of varying N 0 .

Academic example
In this section we report the results of a simple test case, using the geometric parameters listed in the first column of table 1. For simplicity there is no Shafranov shift (γ = 0), and the magnetic lines are horizontal (q = ∞).
To give an idea of the type of functions we are working with, we first provide a representation of a wavelet and of the corresponding vaguelettes ξ λ and χ λ (figure 3). As expressed by (10a), the vaguelette χ λ (figure 3, right) is simply the image of the wavelet ψ λ (figure 3, left) as seen by the camera. In contrast, ξ λ is obtained by solving a linear system, which explains why it features more oscillations on the image (figure 3, middle). Despite these oscillations, ξ λ appears relatively well localized, which is essential to the success of the method, and also to its computational efficiency since it ensures sufficient sparsity (see table 2).
To test WVD-reconstruction, we now apply it to a simple academic test case. We start with an emissivity map which takes the constant value 1 in a toric shell extending from r = 0.47 to r = 0.73, and vanishes elsewhere (figure 4, left). First, we apply the operator K to obtain a synthetic camera image (figure 4, middle). It is interesting to remark the critical curves that appear in the image plane due to the integration along lines of sight intersecting the toric shell. Then, we perturb the image with a Gaussian white noise having a standard deviation σ = 0.5 ( figure 4, right), and apply WVD-reconstruction to reconstruct the emissivity field from the noisy image. The result (figure 5, left) preserves the main features of the input ( figure 4, left), although the consequences of the degradation, such as spurious peaks and oscillations, are visible. In particular, due to the oscillating nature of the wavelets, the approximation of the emissivity which is obtained takes negative values close to the discontinuities, which is why we have extended the colour scale to the interval [−0.2, 1.2] instead of just [0, 1]. As an interesting qualitative test for the method, we also show the denoised image I R = KS R (figure 5, right), where the improvement over the noisy observation ( figure 4, right) is more apparent. Note that when performing an analogous numerical experiment with σ = 0, i.e. in the absence of noise, the inversion is exact up to numerical truncation errors, and the inverted poloidal emissivity map (not shown) cannot be distinguished by the naked eye from figure 4 (left).
The same test case is then used to assess the advantages of WVD compared with other methods. In figure 6(a), the least-squares reconstruction S LS is shown, along with its corresponding camera image, KS LS . Note that S LS also corresponds to an SVD reconstruction with N 0 = N. As a consequence of noise amplification, the initial toric shell is barely visible in the inverted emissivity. By eliminating high order SVD modes (figures 6(b)-(d)), the noise is gradually suppressed, but at the same time the relevant information gets smeared out, and the radial localization of the steps in the emissivity is not well preserved. To obtain an acceptable compromise, one must go through the difficult process of finding the optimal value for the number N 0 of retained modes, and even once this has been done, the result remains qualitatively inferior to the WVD-inverted emissivity.
Admittedly, the discontinuities present in the emissivity make the test case presented in this section rather extreme. However, one is tempted to extrapolate that the same kind of qualitative advantage enjoyed by the WVD over the SVD will hold when more realistic data containing sharp features are considered.

Test with Tokam data
To validate the method further, we apply it to camera images that are generated artificially but using a different method than the one presented above. In order to come closer to the experimental situations we are going to face, we take as input emissivity field the fourth power of the ion density in a computation of edge plasma fluctuations obtained with the Tokam 2D code [32], shown in figure 7(a). The artificial camera image ( figure 7(b)) was obtained by accumulating projections of successive poloidal cross sections through S using Matlab, i.e. the operator K is discretized by splitting the integral (5) into a sum over ϕ instead of a sum over s as we did in the rest of this work. This difference of discretization is a way to test the robustness of our method. The geometric parameters for this test case are shown in the second column of table 1.
For numerical reasons we chose to perform the reconstruction on a grid of size only 32 × 32, whereas the original data were given on a much finer grid of size 512 × 512.
The emissivity map reconstructed by WVD is shown in figure 7(c). For better comparison, we show in figure 7(d) the input data downsampled on a grid of the same size. We see that the radial position of maximum average emissivity close to r = 0.4 is well captured, as well as some of the main blobs, like the one located close to θ = π 6 . Note that the limited resolution of the reconstruction only allows us to detect the presence of this blob, but does not reveal its fine structure. By applying the operator K again to this inverted emissivity map, we obtain an image (figure 7(e)) that is visually similar to the one we started from ( figure 7(b)). This supports the fact that our reconstruction, despite its low resolution, has faithfully extracted a large part of the information which was present in the image ( figure 7(b)).

Application to an experimental movie
In this section we present first results obtained by applying the WVD-reconstruction method to an experimental movie acquired during during the Tore Supra discharge TS42967, of which we first briefly describe the experimental conditions.
During the discharge, a regime of so-called fulldetachment was reached and stabilized over several seconds thanks to a feed-back control of gas puffing on the radiated power fraction (100% of radiated power). In that regime, the edge is strongly cooled down and radiates all the power losses from the plasma core: in the scrape-off layer (SOL), the electron temperature drops to a few eV and steep pressure and temperature gradients (similar to usual SOL profiles) establish at mid-radius (r/a ∼ 0.5) in the closed field line region. The spatial distribution of the radiations was deduced from three diagnostics. From bolometers, the 2D map (r, θ ) of radiation emissivity in the range 200-0.2 nm was reconstructed by tomography on 48 line-integrated channels. In that fully detached regime, the radiations consisted of a uniform shell of several centimeters width at r/a ∼ 0.5 ( figure 8, left). A UV spectrometer with a sweeping line of integration allowed the tomographic reconstruction of the radiations at given UV wavelengths. Emissions from D I (121.6 nm) and C IV (28.9 nm) species also formed a uniform shell at roughly the same position ( figure 8, left).
The strongly radiative toroidal plasma shell is also directly visible on the movies taken by the camera that we consider in this work (figure 8, right), which were acquired with an exposure time of 20 µs and a frequency of 40 kHz. These parameters are marginally sufficient to resolve part of the edge turbulence dynamics. The camera was oriented according to the geometric parameters provided in the last column of table 1, which were in fact estimated a posteriori from the movie itself using a key-point detection method based on some visible features of the vessel. During the parameter fitting process, discrepancies between the experimental movie and our simple camera model were noted, which corresponded to an estimated RMS error of 7 pixels in terms of pixel displacement likely due to optical distortions. In this experiment, the emission spectrum is dominated by the Hα line, around which the fluctuations of the camera sensitivity are of the order of 5%, which we have chosen to neglect in this first analysis.
Thanks to the sufficient resolution in space and time of the movies, previous analyses focusing on the region where Figure 6. Least-squares (a) and SVD (b)-(d) inversion results for the denoising test with a uniform radiating shell from r = 0.47 to r = 0.73, and for varying N 0 , see equation (19). The SVD reconstruction with N = N 0 = 1024 is equal to the least-squares solution (a), equation (16). The colour pictures correspond to S SVD , while the black and white pictures represent KS SVD . White pixels correspond to values falling outside of the colour scale boundaries. the line of sight is tangent to a magnetic surface have allowed extraction of structure advection velocities which are in very good agreement with direct measurement using Doppler backscattering [33]. Here, we apply the tomographic method to reconstruct the emissivity field over the whole sector of the poloidal plane visible to the camera, which is hard to analyse using classical techniques due to the effect of the optical projection.
To illustrate and discuss the results of the method, four consecutive movie frames are considered here (figure 9, left to right). Note that the movie was first cropped to focus on the most active region. In addition, a very simple preprocessing was applied to the movie beforehand, in order to remove vertically propagating horizontal bands (as seen in figure 8, right) which were likely due to an electronic artefact. Moreover, the time average of the whole movie was subtracted from each frame, which helps us to decrease the effect of reflection on the chamber wall. The algorithm is then applied directly to the fluctuations in the signal instead of the full signal ( figure 9, top row).
The inverted maps of emissivity fluctuations are shown in the second row of figure 9. The obtained emissivity fields are quite smooth except for some very intense and localized artifacts, which we interpret as leftovers from the noise. Since those fall out of the colour scale, they can be identified by looking for white pixels in the image. Examination of the four frames reveals that structures tend to propagate counterclockwise. In the last row of figure 9, we show the artificial movie frames obtained by applying K to the inverted emissivity map. The main features that were visible by eye in the original movie frames (figure 9, top row) are strongly enhanced in the denoised one, while the noise has been reduced to a very low level.
From figure 9 it appears that the most intense positive fluctuation, found around z = 0 and r = 0.4, is propagating both poloidally and radially. By measuring the distance travelled between the first and last frames, rough estimates for the poloidal and radial propagation velocities of, respectively, 660 m s −1 and 130 m s −1 are obtained.

Conclusion
We have proposed a new method for reconstructing the volumic light emissivity map of a tokamak plasma using a single camera. Our method relies on the hypothesis that the emissivity varies sufficiently slowly along magnetic field lines. We have demonstrated its feasibility using simple academic test cases, and validated its robustness by applying it to independently generated artificial movies based on numerical computations of plasma edge fluctuations by the Tokam 2D code. We have also highlighted the benefits of the nonlinear denoising capabilities of wavelets as compared with more classical least-squares approaches to ill-posed inverse problems, namely the singular value decomposition (SVD). Finally, we have shown how the method could be applied to an actual experimental movie and reveal the propagation of a structure in the (r, θ ) plane.  The technical tool underlying our approach is the waveletvaguelette decomposition [15] (WVD), which is an efficient way of solving ill-posed inverse problems in the presence of noise. Thanks to the localization of the wavelets, features such as blobs and fronts are preserved in the denoised emissivity map. We have seen that some artifacts persist in the denoised output, but there are hopes that the method could be improved in the future, for example by choosing the threshold in a more refined way. Compared with the classical SVD, WVD yields much better results and appears to be a very promising method Figure 9. WVD-inversion of four consecutive experimental movie frames (left to right) from Tore Supra discharge TS42967. Top row: noisy movie frames used as input for the WVD-algorithm. Middle row: reconstructed emissivity maps in the (r, θ ) plane obtained as a result of WVD-inversion. Bottom row: denoised artificial movie frames obtained by applying K to the reconstructed emissivity maps.
for future application to data from high speed imaging of tokamaks.
The next issue that should be readily addressed is the thorough experimental validation of the method, by comparison with other existing diagnostics, for example Doppler reflectometry [34]. Note that the possible effect of the spectral response of the camera should be carefully investigated in future studies, especially those involving more radiating species, and the possibility of using a filter should be considered. The space resolved emissivity maps that will then be obtained will make new quantities available to experimental investigation, which could facilitate the study of blob velocity, mode coupling, etc, as is already done in linear devices [35].