Multi-Scale Gaussian Normalization for Solar Image Processing

Extreme ultra-violet images of the corona contain information over a wide range of spatial scales, and different structures such as active regions, quiet Sun, and filament channels contain information at very different brightness regimes. Processing of these images is important to reveal information, often hidden within the data, without introducing artefacts or bias. It is also important that any process be computationally efficient, particularly given the fine spatial and temporal resolution of Atmospheric Imaging Assembly on the Solar Dynamics Observatory (AIA/SDO), and consideration of future higher resolution observations. A very efficient process is described here, which is based on localised normalising of the data at many different spatial scales. The method reveals information at the finest scales whilst maintaining enough of the larger-scale information to provide context. It also intrinsically flattens noisy regions and can reveal structure in off-limb regions out to the edge of the field of view. We also applied the method successfully to a white-light coronagraph observation. Electronic Supplementary Material The online version of this article (doi:10.1007/s11207-014-0523-9) contains supplementary material, which is available to authorized users.


Introduction
Extreme ultra-violet (EUV) observations currently provide the most important source of information on the low solar corona. As new EUV instruments are developed, the temporal, spatial, and spectral resolution becomes ever finer, giving new insight into the coronal and chromospheric structure and dynamics. The Atmospheric Imaging Assembly (AIA: Lemen et al., 2011) onboard the Solar Dynamics Observatory (SDO: Pesnell, Thompson, and Chamberlin, 2012) provides very fine temporal and spatial resolution of the Sun at multiple wavelengths and is having a strong impact on the field. Even as the community develop methods to digest the huge volume of data from AIA/SDO, new instruments are planned and tested with even finer resolution (e.g. High-Resolution Coronal Imager (Hi-C), see Cirtain et al., 2013).
Despite the development of automated detection tools for EUV observations (e.g. Martens et al., 2012), most scientific works begin from visual inspection of images. In particular, the volume of data provided by AIA/SDO is so high that low-resolution images are often used as a starting point to find features of interest. The higher-resolution data are then used for further analysis. Image processing is therefore an important step in analysing the data. The scientific return can be improved by the application of processing that better reveals features in the data, particularly in the early stages of analysis where visual inspection is most important.
A common approach to process EUV images is simply to display the square root (or a gamma-curve transformation), or alternatively, the logarithm, of the original pixel values. This is a quick and easy way of reducing the dominance of the image contrast range by a few small bright regions. To reveal dynamic features, time-differencing is used, where a previous image is subtracted from the current image. Such simple processes are commonly used not because of the quality of the output, but because of their simplicity. More advanced image processing methods are not so commonly used owing to their complexity and computational expense.
Wavelet-based techniques have been used by Stenborg and Cobelli (2003) and Stenborg, Vourlidas, and Howard (2008) to greatly improve the visual information available from the Extreme Ultraviolet Imaging Telescope (EIT: Delaboudiniere et al., 1995) aboard the Solar and Heliospheric Observatory (SOHO) and the Extreme UltraViolet Imagers (EUVI: Howard et al., 2008) onboard the Solar Terrestial Relations Observatory (STEREO: Kaiser et al., 2008). The technique involves the decomposition of images into different spatial scales and the filtering/enhancements of features at the multiple resolutions. The technique is computationally expensive but gives very good results. A less sophisticated, yet very efficient technique to reveal features above the limb is based on techniques originally developed for coronagraphs (Morgan, Habbal, and Woo, 2006;Druckmüllerová, Morgan, and Habbal, 2011). Most relevant to this work is the recently developed Noise Adaptive Fuzzy Equalization (NAFE) method (Druckmüller, 2013). The method is inspired by adaptive histogram equalisation, where local statistics govern the output value of a pixel. The method uses a fuzzy membership function of a Gaussian-weighted local set to enhance structural detail, preserve contextual detail, and to reduce noise. The NAFE method results in very clear images without artefacts and with excellent noise reduction. Its one downside is computational expense.
We summarise the problems in visualising features in EUV images (Section 2), introduce the new method (Section 3), apply the method to several images from various instruments (Section 4), and close with a brief summary (Section 5).

Observations
An observation from 04 May 2005 00:00 UT by the 171 Å channel of AIA is used as a working example throughout this section. The data were first pre-processed using the standard SDO Solarsoft software 'aia_prep'. This observation is shown without further image processing in the left panel of Figure 1. The 171 Å channel of AIA has a narrow bandwidth that is dominated by emission from highly ionised iron that has a formation temperature that peaks at ∼ 0.8 MK. The intensity of this image is therefore a function of the emitting plasma temperature and density. The line-of-sight and optical thickness of the plasma also has an effect on the measured intensity, as well as small contributions from other weaker spectral lines that share the same bandpass. The left panel of Figure 1 illustrates the main challenge of revealing information in EUV images. For various physical reasons (mostly due to density), the display range is dominated by the large difference between dark regions and the small bright regions near, or at, the base of active regions. Most of the structural detail appears dark and is hidden. A quick and easy way of revealing some of this structure is by taking the square root of the image (processing of this type is generally known as a gamma transformation). This is shown in the right panel of Figure 1. Although a large improvement on the unprocessed image, there is still much hidden structure, and bright areas become 'washed out'. In particular, there is very little visible structure off the limb. The same is true of applying a logarithmic transformation, or any one-to-one mapping between input and output.
The dominance of the image contrast range by small bright regions can be shown quantitatively by selecting regions and plotting histograms of pixel values. This is shown in Figure 2 for four areas -quiet Sun, active region base, active region, and off-limb. Note that these are pixel values normalised by the exposure time. The top ∼ 75 % of brightness values is taken by the small very brightest region of the active region (cyan). This leaves the information from all other regions in the lower 25 %. In fact, the off-limb and quiet-Sun regions have pixel values restricted to below ∼ 300. Furthermore, the bright active regions contain a large number of pixels that share the low values of the quiet Sun. There are sharp boundaries between very dark and very bright regions. These characteristics become more extreme when flares occur.
Any image processing that attempts to reveal the information hidden in EUV images must deal with the wide range of brightness values, preserve the steep spatial gradients in brightness, and respect the fact that even the brightest regions contain dark patches. Another challenge is the existence of structures at many spatial scales. We cannot apply edge- detection or enhancement exclusively at small scales since this will enhance some features whilst other structures and the larger-scale context is lost. Furthermore, with full-resolution AIA images of 4096 × 4096 pixels taken every 11 s (and the anticipation of the resolution of future instruments), computational efficiency is very important.

Method -Multi-Scale Gaussian Normalisation
Barring spurious calibration errors, the brightness values in the original EUV image will be positive. Let B be the original image normalised by exposure time. k w is a 2D Gaussian kernel of width w pixels in the x and y image dimensions. Throughout this work, w denotes the one-sigma width of the Gaussian. A normalised image C can be computed by where For a given pixel, the numerator of Equation (1) effectively subtracts the local 'mean' (weighted by the Gaussian function centred on the pixel), and the denominator σ w is the local 'standard deviation', also weighted by the Gaussian kernel. Thus, simply speaking, C is locally normalised to a standard deviation of one and a mean of zero. An example of C, computed with w = 20, is shown in the left image of Figure 3.
This process of local normalisation is in essence similar to adaptive histogram equalisation, where each pixel is scaled according to the statistics of the local group of pixels. The main difference is that histogram equalisation enables precise control of the output pixel value range, whereas Gaussian normalisation does not. This problem is solved by means of a further transformation applied to image C, C = arctan(kC).
(3) Figure 3 Left -AIA image locally normalised by Equation (1) with w = 20, giving image C. Right -the same image further processed by Equation (3), giving image C .
Since C is distributed across both negative and positive values, the arctan function in Equation (3)  A set of n locally normalised and arctan-transformed images, C i , is created for n values of w i . Values used throughout this work are w = 1.25, 2.5, 5, 10, 20, 40. A global gammatransformed image, C g is also created by γ can be set at values of 2.5 to 4. Throughout this work, we use γ = 3.2, and this gives good results for all instruments and bandpasses we have processed. a 0 and a 1 set the minimum and maximum input values, respectively. For processing single images, they can be set to the minimum and maximum image values. For batch processing of a large number of observations, they can be set as the minimum and maximum values of the whole data set. For AIA/SDO observations, these values can be easily found at the start of a batch process by reading the appropriate values stored in the data headers (that is, without having to read in the full images). The final processed image, I is calculated by a weighted average of the C i using weights g i , and addition of C g weighted by a global weight h, The g i can be calculated by a study of images consisting of only normally distributed random noise. For narrow Gaussian kernels, the mean of the local standard deviations across the whole image σ w is an underestimate of the global standard deviation, and there is a larger variation in the local standard deviation. This is intuitive -the ratio of the mean of the local standard deviation to the global standard deviation approaches unity with a wider kernel, and there is less variation in the local standard deviation. This is shown in Figure 4. Therefore the weights g i are set at lower values for small w, and approach values near one as w 3, as listed in Figure 4. As described in Equation (5), the original (normalised) image C g is included to give contextual information of the largest scale structure. We used h = 0.7 in this work. In practice, the weights g i and h may be adjusted according to the desired output, and also according to the type of input image (e.g. wavelength or channel). For most purposes, the g i can be set equal for all scales so that a straightforward mean of the locally normalised C i is made.

Pseudocode for MGN
1. Replace spurious negative pixels with zero or local mean/median. 2. Create Gaussian kernel of width w i . Kernel elements should sum to unity. 3. Convolve image with kernel to create local mean image B ⊗ k w . 5. Calculate difference between image and the local mean image, square the difference, and convolve with kernel. Square-root the resulting image to give 'local standard deviation' image σ w (Equation (2)). 5. Calculate normalised image C i by subtracting the local mean image and dividing by the local standard deviation image (Equation (1)). Store the result. 6. Apply arctan transformation on C i to give C i . 7. Repeat 2 -6 with the different kernel widths w i . 8. Take mean, or weighted mean if preferred, of the C i to give a weighted mean locally normalised image. 9. Calculate a global gamma-transformed image C g by applying Equation (4). 10. Sum the weighted mean locally normalised image with the global normalised image C g , with appropriate weight h (as Equation (5)).
Note that considerable efficiency is gained by creating a one-dimensional Gaussian kernel for convolving first along the x-direction, then convolving the resulting image along the y-direction with a transposed kernel. This is far more efficient than convolving directly with a two-dimensional kernel. Compared with the NAFE (Druckmüller, 2013), or wavelet-based routines, the MGN method is a simple one, and as Gaussian smoothing is a standard practice employed by most programming languages, can readily be programmed in a few lines of code. Without great rigour, we tested the efficiency of the processing using the Interactive Data Language (IDL) on a MacBook Pro laptop with 8 Gb Ram and a 2.6 GHz Intel Core i7 processor. The computing time increases approximately linearly from ∼ 1 s for a 500 × 500 image to 10 s for a 2K × 2K image. A full AIA/SDO 4096 × 4096 image takes 40 s. This is extremely fast compared with other methods such as the NAFE (at least an order of magnitude faster), and does not require a high-performance computer. Figure 5 shows the result of applying MGN to the example AIA image of the previous sections. Structure is enhanced down to fine spatial scales, although large-scale context is preserved (that is, bright regions are still brighter than dark regions). Off-limb structure, where present, is enhanced without being swamped by noise. Structure is enhanced in the extremely bright region at the base of an active region (the cyan-bordered region in Figure 2), which originally appeared saturated. Structure can also be seen in dark regions immediately neighbouring the bright region. It is also very useful to be able to trace structure from origins on the disk to off-limb heights. In this example, there are loops and fan-like bright features that extend from the disk and past the limb. Movies created using this processing reveal dynamic features at very small scales such as flows within filament channels, flows along large coronal loops, small-scale activity within active regions, and movement in extended off-limb structures in response to dynamic events.

Results
An example AIA movie is available online (aia_movie.gif -an animated gif that can be viewed with a web browser), which shows a large portion of the western hemisphere and off-limb region observed in the 171 Å channel during 18 January 2013. Activity and smallscale flows are clearly seen, alongside larger-scale eruptions, and the movement/expansion of loops off the limb out to the edge of the field of view.
Another online movie (comparison_movie.gif) compares the appearance of a sequence of AIA observations processed by the MGN (left panel), and by a simple gamma transformation (right panel). This sequence shows a flaring active region and the eruption of a CME observed in the 171 Å channel during 08 March 2011 in the South-West. The disk and offlimb structure is seen far more clearly in the MGN images. Individual loops are resolved, even in the complicated off-limb region to the North of the active region. In the original data, the sequence is dominated by the bright core of the active region, and when the flare peaks, by the small saturated regions centred on the flare. The MGN images show this activity clearly, and also succeed in maintaining a clear and stable view of surrounding structure throughout the movie, including the 'wobble' of the active-region loops in response to flaring activity. This enables a more comprehensive analysis of the event, showing connections between activity and movement of structure that can be hidden without MGN processing. Figure 6 shows a High resolution Coronal Imager (Hi-C) image processed with MGN. A movie showing the northernmost portion of this image is available online (hic_movie.gif). Various features such as magnetic braiding (Cirtain et al., 2013) are clearly seen, as are plasma streams along the dark filament channels, and small-scale movement within the magnetic 'moss'. Such features cannot be seen without appropriate processing. Figure 7 shows the application of MGN to an observation by the Sun Watcher using the APS (SWAP) instrument onboard the PRoject for OnBoard Autonomy 2 (PROBA2) satellite (Berghmans et al., 2006;Seaton et al., 2013). This instrument has a coarser resolution than AIA, but has an extended field of view. Figure 7 reveals an erupting filament out to the extremity of the field of view, and other quiescent structures to ∼ 1.5 R . Again, even low-signal structures are enhanced without too much amplification of noise.
The MGN processing is not limited to EUV observations. Figure 8 shows its application to a white-light coronagraph image of the extended inner corona by the Large-Angle and Spectrometric Coronagraph (LASCO) onboard SOHO. The left image shows an observation of 14 January 2011, with a long-term background subtracted to remove instrumental straylight and F-corona. A point filter has also been applied to reduce spurious bright pixels. The right image shows an MGN-processed image, with parameters k = 0.8, h = 0.9 and γ = 1. Smaller-scale features are enhanced, in particular faint plumes over the poles that are difficult to see in the original image. In this respect, the MGN is an improvement over the NRGF (Morgan, Habbal, and Woo, 2006), although the NRGF provides structural context closer to the true K-corona. It is difficult to make a qualitative comparison of images created using different processing. Different processes may serve different analysis purposes -the important function for the process described here is to reveal information hidden in the original images. A non-rigorous comparison 'by eye' of the same images processed using the MGN and of the NAFE (Druckmüller, 2013) suggests that the results are similar, with the NAFE giving slightly better clarity overall, and improved noise suppression -most obvious at large heights off the limb. For the MGN, noise reduction arises naturally from taking the (weighted) mean across many spatial scales. The noise suppression is not as effective as the NAFE, and the MGN method lacks control over the degree of noise suppression. The NAFE method is also better suited to reveal contrast in very bright regions.

Summary
The MGN method normalises an image by using the local mean and standard deviation calculated using a Gaussian-weighted sample of local pixels. This normalised image is transformed by the arctan function (similar to a gamma transformation). This is applied over several spatial scales, and the final image is a weighted combination of the normalised components. The results compare well with multi-resolution wavelet enhancement or the NAFE  procedure, but is far more computationally efficient. We hope that the MGN will become an established tool for researchers, offering a good compromise between computational time and clarity of the final images. The method is simple to implement, and the lead author is happy to provide the IDL code by email request.