Article 3D Reconstruction of Coronal Loops by the Principal Component Analysis

Knowing the three dimensional structure of plasma filaments in the uppermost part of the solar atmosphere, known as coronal loops, and especially their length, is an important parameter in the wave-based diagnostics of this part of the Sun. The combination of observations of the Sun from different points of observations in space, thanks to the most recent missions, including the Solar Dynamics Observatory (SDO) and the Solar TErrestrial RElations Observatory (STEREO), allows us to infer information about the geometrical shape of coronal loops in 3D space. Here, we propose a new method to reconstruct the loop shape starting from stereoscopically determined 3D points, which sample the loop length, by principal component analysis. This method is shown to retrieve in an easy way the main parameters that define the loop, e.g., the minor and major axes, the loop plane, the azimuthal and inclination angles, for the special case of a coplanar loop.


Introduction
The dynamics of the uppermost part of the atmosphere of the Sun, known as the solar corona, is one of the key elements of the magnetohydrodynamic (MHD) interaction of the Sun with the Earth and other planets and is important for basic plasma physics.Understanding how the solar corona is structured in three dimensional (3D) space is a challenging and not a trivial task.Excluding solar eclipses, which allow us to see the corona in visible light for only a few minutes, routine observations from telescopes or coronagraphs aboard spacecraft (such as the Solar and Heliospheric Observatory (SoHO), the Transition Region And Coronal Explorer (TRACE), the Solar TErrestrial RElations Observatory (STEREO), Hinode and, more recently, the Solar Dynamics Observatory (SDO) at extreme ultra-violet (EUV) and X-ray bandpasses are the main source of information about the complex dynamics of the corona.In these bands, the corona is usually an optically thin medium: an emitted photon can go through it without experiencing any noticeable absorption or scattering.This seriously affects the identification or tracking of important features (e.g., coronal loops and plumes) in imaging data obtained from a single point of observation.A crucial parameter in this sense is the column depth [1,2], which gives a measure of how much the intensity of a given structure is integrated along the line-of-sight (LOS).
The launch of STEREO in 2006 has opened the possibility to make stereoscopy, thanks to the twin spacecraft observing the Sun at about 1 AU, but with two different LOS [3].The STEREO mission's aim is to understand the 3D structure of coronal mass ejections (CMEs) and their impact on the solar system and Earth environment (see [4]), but it is suitable also for 3D reconstruction of the geometry of loops in coronal active regions (ARs) (for a review of the topic, see [5]).In particular, observations show that coronal loops are subject to transverse (kink) oscillations (e.g., [6][7][8][9][10][11]).These oscillations are a good tool for the diagnostics of the magnetic field inside the loop, which depends on the period of the kink oscillations, its length, the plasma density [12] and the loop's sub-resolution structuring [13].Therefore, a better estimate of the loop length could easily improve the inference of the field.The value of the field is one of the decisive parameters of AR for answering the enigmatic questions of solar and stellar physics, such as the mechanisms for coronal heating and flaring energy releases, coronal mass ejections and the solar and stellar wind acceleration.
A forward-modelling 3D loop reconstruction method has been developed by Verwichte [14,15] in the context of studying kink oscillations of coronal loops and reducing the error in measuring loop lengths.The method works as follows.In one view, the plane-of-sky loop coordinates are traced.A third coordinate that denotes the loop observed depth is added by assuming that the loop is planar (or non-planar, following a simple geometrical model).The problem is thus reduced to finding the value of one (or a few) parameters, such as the inclination angle of the loop plane with respect to the normal to the solar surface.For a range of inclination angles, the loop is forward-modelled to a new view, and the angle with the best fit is determined visually.In [14,15], a second viewpoint was absent.In the first study, a second viewpoint was created from the data from several days later, using the change of the LOS by the solar rotation.In many respects, this method is then akin to the dynamic stereoscopy method by Aschwanden [16].In the second study, the inclination angle is found by minimising the variation in curvature.Later, this method had been applied to loops seen simultaneously with SDO and STEREO [10,[17][18][19].This method has the advantage of not requiring accurate tracing of the loop from both views and does not impose a specific shape of the loop (e.g., semi-circular), except for planarity.Because of the difference in the spatial resolution between the Atmospheric Imaging Assembly (AIA) in SDO and the Extreme Ultra Violet Imager (EUVI) in STEREO, tracing a loop in the AIA view and forward-modelling to the EUVI view gives the best results.
In a strict sense, 3D reconstruction is performed by observations of a solar feature from at least two observers placed at different positions-a technique that is called "stereoscopy"-and allows for the determination of the 3D coordinates in the space of an observed point (for an explanation of the principles of stereoscopy, see [5,20,21]; a description is given at [22], The SolarSoftWare (SSW) package [23] provides useful tools for making 3D reconstruction, e.g., the procedures scc measure.pro[24] or sunloop.pro[25].The former includes a widget application (see Figure 1 for a typical graphical output) that enables the user to select with the cursor a solar feature of interest appearing in both perspectives from one of the STEREO spacecraft (but it is also possible to combine SDO and the one of STEREO spacecraft).By selecting a point in one image, the program displays a line in the other image representing the LOS from the first image (or the epipolar line).The user then selects the point along this line with the same feature.The 3D coordinates are then calculated as (Earth-based) Stonyhurst heliographic longitude and latitude, along with the radial distance in solar radii.The coordinates can be easily converted into the Heliocentric Earth Equatorial (HEEQ) coordinates for a Cartesian representation of the data (for a complete review of the coordinate systems, please refer to [26]).As we have seen, fitting a loop with a 1D curvilinear feature (circular or elliptical) depends on six parameters, namely: the heliographic longitude and latitude of the midpoint of the loop baseline on the solar surface, the baseline length between the two footpoints, the height of the loop, the inclination and the azimuth angle (see Section 3.4.4 of [27]).
A representation of the loop and its main parameters are shown in Figure 2. Recently, Gary et al. [28] proposed a fitting approach of loops by cubic Bézier curves, in order to distinguish between competing models for the magnetic field topology.In this paper, we propose a new method allowing us to get the 3D shape of a loop by stereoscopic measurements, using principal component (PC) analysis.It reduces the number of parameters, dependencies and correlations between variables, by calculating a new basis of vectors, which defines the reference system relative to the loop.In the next section, we give an overview of PC analysis in the frame of 3D loop reconstruction.In Section 4, we present some examples, discussing analysis and comparing results with previous methods, and conclusions are given in the last section.

Three-Dimensional Reconstruction by Principal Component Analysis
Given a set of variables, which can be interrelated, PC analysis allows us to transform this set into a new one of "uncorrelated" or independent variables, known as principal components, but preserving the variation of the data distribution.Principal component analysis, also know as minimum variance analysis, has been applied in several research fields, e.g., in the determination of magnetic field geometry by in situ multi-spacecraft measurements [29].Algebraically, PCs are particular linear combinations of the set of original variables, and geometrically, this is equivalent to the rotation of the original coordinate system into a new one, where the new axes represent the maximum variability of the dataset.Avoiding a rigorous mathematical treatment, which can be found in textbooks, such as Jolliffe [30], Johnson and Wichern [31], we give here a specific step-by-step description of how to fit coronal loops by PC analysis.
Given a set of N 3D data points x i = (x i , y i , z i ) (we use, for simplicity, the HEEQ system) that sample the loop shape and is represented by an N × 3 matrix X, we can build the sample covariance matrix, S, defined as: where σ 2 q,q = N i=1 (q i − qi )(q i − qi )), with q = {x, y, z} and q the mean, identified as the centre C = ( x, ȳ, z) of the loop.Now, we would like to put this in a new reference frame, whose axes maximise the variances, minimise the covariances of the points and have means equal to zero.This is equivalent to finding the eigenvalues and eigenvectors that diagonalise the matrix, S, in a matrix, Λ.Let us obtain the eigenvalues and sort them in the ascending order, as: The corresponding eigenvectors will be e n , e a , e b , whose components are relative to the original frame.Thus, the matrix: e n,x e n,y e n,z e a,x e a,y e a,z e b,x e b,y e b,z defines the basis change X = EX l + C, and diagonalises S as Λ = E T SE (since E is the basis change matrix, E T = E −1 ); x l,i with i = 1, ..., N, are the coordinates of the data points in the new reference frame (they are also called "scores" in PC analysis); the eigenvector, e n , locates the normal to the loop plane, e a is directed like the minor axis of the ellipse and e b , like the major one.The axes values are related to the variance as (for statistical reasons, when N is small, it can be replaced by N − 1): Indeed, since we have defined the covariance matrix, S, without dividing its elements by N, this should be done for the eigenvalues (the eigenvectors remain the same despite the normalisation of S).A further justification for the presence of the factor, √ 2, is given in Appendix A.1.The "fitted" loop can be traced starting from the parametric equations of an ellipse: with t = [0, 2π], then, transformed in HEEQ coordinates, x e , y e , z e , by the basis change with the matrix E.
Those points, whose squared distance satisfies the condition r 2 = x 2 e + y 2 e + z 2 e ≥ 1.0R , defines the loop curve.
We can look at how much of the total variance of the data points is counted by each single eigenvalue.Let Λ be the total variance defined as Λ = λ n + λ a + λ b .Then, the ratio λi = λ i /Λ (with i = n, a, b) gives an indication of how much of the "energy" is carried out by λ i .According to Table 1, most of the variance is stored in along the e a,b directions, suggesting that the points are mostly distributed in a 2D plane.This is equivalent to saying that we have reduced the dimensionality of the dataset, by projecting it into a subspace, that is the 2D plane of the loop.
Table 1.The table shows the main quantities used and inferred from principal component (PC) analysis: the first rows are related to the inputs and show the number of 3D points, N, and the heliographic coordinates of the supposed loop centre; then, we tabulated the total variance, Λ, and the percentage relative to each eigenvalue, λ.The last rows show the estimates of each analysed loop: the minor (a) and major radii (b) with their corresponding standard deviations, the standard deviation, σ n , of the points from the loop plane, the inclination and azimuth angles, the eccentricity and, finally, the loop length.We can obtain the inclination and azimuth angles from the orientation of the loop plane in space.The inclination angle is defined as the angle between the normal to the solar surface and the loop plane.By considering the vector, e n , normal to the loop plane, we can find the inclination angle as: where e r is the normal, directed in the radial directions.The azimuthal angle is defined as the angle formed by the intersection of the loop baseline with the east-west solar direction, represented by the longitudinal vector, e φ .Therefore, the azimuthal angle can be found as: The errors of the fitting can be estimated from the "scores" and their distance, d, from the fitted curve; it is possible to give an estimate of the standard deviation of the a and b direction from the sum of the squared distance (see A.2, for more details).Another possible way is to perform PC analysis for M samples of the 3D points, or to construct synthetic samples by randomly varying the 3D original measurements within their errors (private communication: Antonio Vecchio): the best values of the radii, inclinations and azimuth angles, as well as their standard deviations, will be inferred from their distributions.

Analysis and Discussion
In this section, we demonstrate some applications of PC analysis to specific coronal loops, with some comparisons with previous analysis.

The Event on 27 June 2007
The coronal loop in Figure 3 appeared on the east solar limb and was subject to kink oscillations, triggered by a flare.An analysis of the loop geometry is given in Verwichte et al. [14], using a forward modelling technique based upon the assumption of the circular shape, and in Aschwanden [5], with a fitting approach of the 3D stereoscopic points.The loop centre is not clearly identified, because of its proximity to the solar limb.The knowledge of the possible footpoints in Verwichte et al. [14], was inferred by looking at the same active region three days later, when it passed the east solar limb.
Here, we deduce it directly from the stereoscopic measurements by calculating the centre of "mass" of the distribution of points and projecting the highest point (that can be regarded as the "apex" of the loop) to the solar surface, crossing the centre.Following this idea, we extrapolate a heliographic longitude of about φ C = −87.76degand latitude θ C = −11.41deg for the candidate loop centre.The results of the PC analysis are given in the first column of Table 1.We got an inclination of θ = −2.85deg and an azimuth angle α = −37.21deg.The relative errors for the minor and major radii are σ a /a = 5.6% and σ b /b = 6.1%, respectively.The loop length is about 325 Mm, and it is close to the value estimated by Aschwanden [5], Verwichte et al. [14].Moreover, it exhibits a strong eccentricity (∼0.9), indicating an elliptical geometry rather than a circular one.
We overplotted the data points and the reconstructed loop over the AR seen three days later, on 30 June 2007.The loop seems to agree with the underlying topology of the AR, especially for the position of the footpoint labelled by green in Figure 4.

The Event on 21 January 2013
Here, we show a loop appearing in the quiet region on 21 January 2013, close to the western solar limb.The loop is seen well with SDO/AIA at 171, 193 Å, and visible in STEREO/EUVI-A, as well, but with a spatial resolution 2.5-times lower.We were able to take 3D points with scc measure: we worked on filtered images (we subtracted the frames with those 5 min previous) in order to highlight the loop shape (Figure 5-top).Since the footpoints are well identified, we determined the centre as the midpoint between footpoints.We took 10 independent measurements and considered their mean as the best value for the centre.The outcomes of PC analysis are listed in Table 1.

The Event on 22 January 2013
Here, we present the 3D reconstruction of coronal loops belonging to the active region AR NOAA (National Oceanic and Atmospheric Administration) 11654, studied by Anfinogentov et al. [32].The active region exhibits a large number of coronal loops.From SDO/AIA, the active region is partly off-limb, and loop footpoints are not visible.The same active region is seen well from STEREO-A: the footpoints are easily identified, and the loops seem to span for a wide angle across the surface normal.Therefore, intuitively, one could assume that the only parameter that is changing is the inclination angle, θ, while the other properties, like the loop lengths and the azimuth angles, α, should remain almost constant.We use running difference images, in order to identify the loop in both FOVs when we note a simultaneous intensity variation.Figures 6-8 show the reconstruction for some of these loops, which were located with good certainty.The parameters are listed in Table 1; for all three loops, note the consistency in the parameters, like the minor and major radii, the eccentricity, the azimuth angle and the loop length.Furthermore, one can see the discrepancy in the estimation of the inclination, suggesting that the analysed loops become gradually more perpendicular to the solar surface.The loop in Figure 9 is different, showing different geometrical values and a stronger inclination.Entropy 2013, 15

Conclusions
In this paper, we propose a new method for the determination of the 3D shape of coronal loops from stereoscopic measurements.The necessity to understand the loop shape is important in order to better estimate, for example, the loop length, which is a crucial parameter in inferring the magnetic field by coronal seismology (or, generally, to compare the geometry of loops, which are assumed to trace magnetic field lines in the low-β coronal plasma, with potential field models of the coronal magnetic field [28]).Furthermore, the loop shape determines the variation of the column depth along the loop in different phases of kink oscillations and, hence, the distribution and dynamics of the LOS-integrated intensity of the loop.The previous approach to 3D reconstruction is based on fitting stereoscopic 3D tie-points, taken along the loops.A good fitting requires a large number of points, but taking them could be a laborious task, because of high observational and instrumental noise and the effect of the LOS integration.Usually, by procedures, e.g., scc measure, the number of tie-points selected is about 10-20.The method used in Aschwanden et al. [21] and his subsequent papers is to interpolate the 3D points with a spline, in order to retrieve the loop shape.On the other hand, this operation could add further degrees of freedom (e.g., non-planarity), and the corresponding reconstructed loop could deviate from a regular shape, like a circle or ellipse (e.g., see Figure 9 of Aschwanden et al. [21], and Figure 7 of Aschwanden [5]), also because of the natural spread of the measurements around the true value.Moreover, fitting the 3D points with a 1D curvilinear feature requires the determination of the six free parameters: the inclination and azimuth angles, the length of the baseline, the loop centre position and the loop height.Conversely, all this information is implicitly stored in the distribution of a sample of 3D measurements, and it can be retrieved easily by the technique of principal components analysis.Determination of these parameters can be made by a simple step-by-step process.The technique we developed can be summarised in the followings steps: • measurement of a sample of 3D tie-points that outline the loop shape from two different LOSs; • calculation of the sample covariance matrix, S, of the 3D points; • determination of the three pairs of eigenvalue/eigenvectors, λ i /e i , which identify the new reference frame relative to the loop, sorted in ascending order; • the smallest eigenvalues and the corresponding eigenvector is associated with the normal to the loop plane; • the largest eigenvalue/eigenvector locates the major axis and the remaining, the minor axis, according to Equation (4); • determination of the inclination and azimuth angles according to Equations ( 6) and ( 7); • loop tracing with the parametric Equation ( 5) and transformation in the original coordinates system, HEEQ.
The determination of the set of three vectors, e n , e a , e b , which constitutes the new basis for the reference frame relative to the loop, depends uniquely on the choice of the loop centre (or, also, called the midpoint), which plays the role of the decisive parameter of the problem.Moreover, according to the example of the loop on 27 June 2007 (Figure 3), we have shown that the information about the midpoint can be deduced or guessed directly from the 3D measurements, overcoming the issue of the lack of information about the loop baseline, when it is too close to the solar limb.The advantage of this method is that it works for a reasonably small number of data points (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20), without requiring any intermediate step, like interpolation.The only assumption made in our analysis is the planarity of the loop, but possible speculations on it can be made by studying the distribution of the data points around the loop plane and searching for some possible patterns, like arcs or s-shapes, which could outline, respectively, the fundamental and second harmonic in kink oscillations.On the other hand, the relative deviation of the measurements from the best-fitted loop plane is less the 1%, probably too small in order to appreciate any non-planarity.Moreover, a further step forward for reasonable modelling of the magnetic field would be to fit the loop shape by a dipole or stretched dipole line, as suggested in Hu et al. [33].These issues, as well as refinements of the technique by PC analysis could be the object of future work.measures the squared distance from x to µ in the standard deviation units.For N observations of a vector, x, this can be represented by: 1 with µ being the expected value and Σ being the variance-covariance matrix of the x observations (the apex indicates the complex conjugate of (x − µ), and similarly in the next equations, as well).The above quantity for a multivariate distributions represents the contour of constant density probability at specific values k 2 , or the surface of an ellipsoid centred in µ.In our case, given the x i observations, the centre, C, and the sample covariance matrix, S, we have: because of the principal components reduction.Therefore, if we expand the right side of Equation ( 11), we get the equation of a 3D ellipsoid: Since the eigenvalues, λ n , do not significantly contribute to the total variance, they can be eliminated, and we get the equation of an ellipse in the canonical form, with k = 1, a =

A.2. Error Estimations
The errors in the fitting can be estimated from the "scores" or the distance from the fitted curve.Let x l = (x l,a , x l,b , x l,n ) be the data point in the loop reference frame and x e = (x e,a , x e,b , x e,n ), a point of the fitted curve (that should be determined).The squared distance, d 2 , will be: ) 2 (18) The standard deviation along a, b and n can be estimated as: N − 1 with j = (a, b, n) A.3.Error Estimate on the Loop Length The loop length can be estimated numerically given the radii, a and b, and N p grid points along the loop as: The uncertainty of the length can be expressed as (note that, in general, the inequality ln(x + y) ≤ ln(x) + ln(y) is true when x ≥ y/(y − 1)): For very good fittings, we obtained the uncertainties, σ a /a and σ b /b, of the order of 3-5%, so that the loop length error should be less than 8 %.

Figure 1 .
Figure 1.Graphical output of the procedure scc measure: the windows give the field of view (FOV) from STEREO-A ( right) and STEREO-B ( left) spacecraft for a coronal loop on 27 June 2007, studied by Aschwanden [5], Verwichte et al. [14].

Figure 2 .
Figure 2. Schematic representation of a coronal loop relative to the solar surface, with the baseline connecting the footpoints, the inclination angle relative to the normal surface and the azimuth angle relative to the east-west direction.Courtesy of Markus Aschwanden.
Inclination angle θ F o o tp o in t b a s e li n e V e r t ic a l Definition of Loop Parameters

Figure 3 .
Figure 3.A coronal loop on the east solar limb during a flare event on 27 June 2007.The event was studied by Aschwanden[5], Verwichte et al.[14].The top figures show the FOVs from STEREO-A (right) and STEREO-B (left).The yellow symbols are the 3D data points measured with the routine scc measure.The light-blue line is the fitted curved loop; at the base, the centre is in blue, and the footpoints inferred from the reconstruction are in red and green colours.The middle plots show the reconstructed loop against the yellow data points, for the different orientations of the Heliocentric Earth Equatorial (HEEQ) coordinates system.The coloured vectors represent the new set of three axes: in red and blue, the minor and major axis; in green, the normal to the loop plane.Indeed, the bottom figures show the loop in this new reference frame for different orientations.A similar description is applied for the Figures of the other analysed events.

Figure 4 .
Figure 4. Image from STEREO-B at 171 Å of the active region, which was the site of the analysed coronal loop, on 30 June 2007.The 3D points and the "fitted" loop are overplotted in order to visually prove the consistency of our results.Note that there is a good agreement with the position of the top-right footpoint with the open loops of the active regions.The red meridian marks the position of the solar limb on 27 June 2007.

Figure 5 .
Figure 5.The loop event observed with SDO/AIA (left) and STEREO-A (right) in a quiet region on 21 January 2013, in the Southern hemisphere.

Figure 6 .
Figure 6.Analysis of the loop (a) on 22 January 2013.

Figure 7 .
Figure 7. Analysis of the loop (b) on the 22 January 2013.

Figure 8 .
Figure 8. Analysis of the loop (c) on the 22 January 2013.

Figure 9 .
Figure 9. Analysis of the loop (d) on the 22 January 2013.

2 a + d 2 b
+ d 2 n = (x l,a − x e,a ) 2 + (x l,b − x e,b ) 2 + (x l,n − x e,n ) 2(14)Let us normalise the distance along the axes to the radii, a and b:d 2 a = a( xl,a − xe,a ) 2 d 2 b = b( xl,b − xe,b ) 2(15)Now, we have reduced the loop to a circumference of radius 1.The coordinates, xe,a and xe,b , will be, respectively, equal to the cosine and sine of the angle spotted by the point, ( xl,a , xl,b ), so that we rewrite as: xe,a = cos γ = ,b = sin γ = xl,b x2 l,a + x2 l,b(17)and hence: d 2 a = a xl,a (1 − 1 x2 l,a + x2 l,b ) 2 d 2 b = b xl,b (1 −