A NOISE ADAPTIVE FUZZY EQUALIZATION METHOD FOR PROCESSING SOLAR EXTREME ULTRAVIOLET IMAGES

Launched in 2010, the Atmospheric Imaging Assembly (AIA) aboard the Solar Dynamics Observatory (SDO) produces high resolution and high dynamic range images of the corona in the extreme ultraviolet (EUV) with unprecedented spatial resolution. To take full advantage of the information imbedded in these images, two types of products are necessary: (1) calibrated images necessary for the quantitative analysis of the coronal plasma and (2) images reflecting the instrumental resolution of approximately 1'', ideally suited for the qualitative studies of coronal structures created by the solar magnetic field. The second type of images is the focus of this paper.

The dynamic range of AIA images is up to 14 bits per pixel (16,384 brightness levels). There is typically an extreme contrast between bright features such as loops in active regions and a dark background which contains low contrast structures created by relatively weak magnetic fields. This contrast is so high that no monotonous pixel value transform y = t(x) (where x is the input and y the output pixel value) is able to produce an image suitable for computer display or for printing (typically 8 bit per pixel, i.e., 256 brightness levels), namely an image which shows all features with a contrast sufficient for human vision. That is why a compromise is usually made, whereby the processed images contain saturated parts in order to expose the dark, low contrast, features.

A classical method for the visualization of fine low contrast structures is based on attenuating low spatial frequencies by means of techniques based on the Fourier transform (or other orthogonal transforms) or convolution (Pratt 2004). If we use the convolution of the original image A with the kernel C for that purpose, the resulting image B may be written in the form

$$\begin{equation} B = A*C = A*(D-L) = A - A*L = A - M \end{equation} \tag{ 1 }$$

where * denotes convolution, D denotes a Dirac kernel (i.e., a kernel with the property A*D = A), L denotes the kernel of a low-pass filter. The image M is usually called the unsharp mask. (That is why the method is often called unsharp masking.) As the highest contrast is typically in the low spatial frequencies, the resulting image B has usually a significantly lower dynamic range than the original image A, and is therefore more suited for visualization. Unfortunately, this is not the case for AIA images, where the high contrast features are typically in the high spatial frequencies (for example, bright loops). That is why this classical method is not very effective.

A number of other techniques have been recently developed for processing high dynamic range images obtained with coronagraphs (SOHO LASCO, STEREO SECCHI) or during total solar eclipses (Morgan et al. 2006; Druckmüller et al. 2006; Druckmüllerová et al. 2011). These methods give very good results for images of the Thomson scattered K-corona emission and the emission line corona (Pasachoff et al. 2007, 2008; Habbal et al. 2010, 2011). They are based on the assumption that both types of emission have an extreme brightness gradient in the radial direction, as measured from the Sun center. Removing this gradient significantly decreases the dynamic range of the images. However, these methods are usable only for processing the parts of SDO AIA images, which extend off the limb, where the assumption of a steep radial gradient is valid.

There exists another class of methods, which is based on the adaptive histogram modification, in particular, the adaptive histogram equalization (Pratt 2004). These methods solve the problem of high dynamic range and extreme contrast very well but they suffer from two critical problems. The first problem is the appearance of artifacts caused by the shape of the neighborhood (usually a square) from which the histogram is computed. These artifacts may influence coronal structures significantly, leading to false interpretation. The second problem is the extreme amplification of noise in the low contrast parts of the image, resulting in the loss of faint details.

The Noise Adaptive Fuzzy Equalization method (NAFE), which is described in this paper, was inspired by both the adaptive histogram equalization and unsharp masking, and is in some sense a combination of both methods. The method does not suffer from the problems typical of these two methods when applied separately. As we show in what follows, the resulting images have a dynamic range which is suitable for viewing, they are free of processing artifacts and the noise may be well controlled.

Let A denote the original input image. In this image the pixel values have a linear dependence on the intensity of the coronal emission. Let B denote the resulting image. This image is a linear combination of two images: T_γ(A) and E_{N, σ}(A), such that

$$\begin{equation} B = (1-w) T_\gamma (A) + w E_{N,\sigma }(A). \end{equation} \tag{ 2 }$$

T_γ(A) denotes the so-called gamma transform of image A, which is produced by the pixel value transform

$$\begin{equation} y = b_0 + (b_1 - b_0) \left(\frac{x - a_0}{a_1 - a_0}\right) ^{\gamma ^{-1}}. \end{equation} \tag{ 3 }$$

The constant a₀ is the input value representing 0 emission intensity. The constant a₁ is the maximum valid input value. Constants b₀, b₁ are minimum and maximum output values, respectively. The standard γ value for a PC computer display is γ = 2.2. Lower values for γ give darker images, and higher values give brighter images. Image T_γ(A) in formula (2) represents the "unmodified" image upon which only a nonlinear pixel value transform was applied. E_{N, σ}(A) is an image created by the NAFE method, to be described in what follows. The constant w will be called the NAFE weight. Typical values for w lie in the interval 〈0.05, 0.3〉. Constant w = 0 gives an image B without any enhancement. The specification of a value for w enables one to control the enhancement of the structures.

The NAFE method produces a pixel value transform which is not constant for the whole image like a gamma transform. It is different for every pixel and is dependent on the neighborhood of the pixel. An inspiration for this method was the histogram equalization method which is based on the idea that the optimal brightness level distribution in the image is the uniform one, i.e., all pixel values in the resulting image are used with the same probability.

Let us assume, for simplicity, that pixel values in image A are realizations of a continuous random variable V with a distribution function F(x). The pixel value transform y = mF(x) creates the image B pixel values which are realizations of a random variable U with a uniform distribution on the interval 〈0, m〉 (Hogg et al. 2005). In reality, the random variable V is discrete, which is why U is only an approximation of a uniform distribution. The images to be processed typically have several thousand discrete pixel values. Therefore the approximation is very good and the difference between a continuous and a discrete variable is unimportant. Because we do not know the random variable V, but only its realization (the image), the distribution function F(x) must be estimated by the normalized cumulative histogram of the image.

The adaptive histogram equalization (Pratt 2004) is based on the idea that the cumulative histogram is not computed for image A as a whole, but only for the processed pixel neighborhood. In this case the pixel value transform function in not a compromise derived from the whole image but it is optimally set according to the pixel neighborhood properties. The typical size of the neighborhood is 10⁴ or more pixels, which makes the algorithm extremely time consuming. That is why the square neighborhood is usually used, which enables the use of a recurrent algorithm for histogram computation.

The result of the adaptive histogram equalization is highly dependent on the size and shape of the neighborhood and the resulting image may contain severe artifacts. The reason may be explained by means of an analogy with convolution—formula (1). Let us denote the square neighborhood by $N_{i,j}^n$, where [i, j] are coordinates of the pixel, and the odd integer n is the width of the neighborhood in pixels. $N_{i,j}^n$ plays the analogous role of the convolution kernel L in which all elements are identical. Such a kernel is in principle a very bad low pass filter with a complicated Fourier spectrum which is direction-dependent. That is why kernels with much better properties are used—for example, a Gaussian kernel L. A significant improvement of the adaptive equalization may be achieved by using the fuzzy neighborhood $\widetilde{L}_{i,j}^n$ instead of the neighborhood $N_{i,j}^n$.

Let us denote by L = [l_{k, l}] the matrix of size n × n (in analogy with L in Equation (1)) elements which are in the interval 〈0, 1〉 and

$$\begin{equation} k, l \in \left\lbrace \frac{1-n}{2}, \dots, 0, \dots, \frac{n-1}{2} \right\rbrace. \end{equation} \tag{ 4 }$$

The fuzzy neighborhood $\widetilde{L}_{i,j}^n$ is a fuzzy set (Zadeh 1965) with support $N_{i,j}^n$ and membership function $\mu _{i,j}^n: N_{i,j}^n \rightarrow \langle 0,1 \rangle$ where the membership grade of pixel a_{i + k, j + l} to the fuzzy neighborhood is defined as

$$\begin{equation} \mu _{i,j}^n(i+k,j+l) = l_{k,l}. \end{equation} \tag{ 5 }$$

Now let us define the fuzzy histogram of $\widetilde{L}_{i,j}^n$ as

$$\begin{equation} h_{i,j}^n(x) = \sum _ {k=\frac{1-n}{2}}^{\frac{n-1}{2}} \sum _ {l=\frac{1-n}{2}}^{\frac{n-1}{2}} l_{k,l} \, \delta _{x, a_{i+k,j+l}} \end{equation} \tag{ 6 }$$

where δ denotes the Kronecker delta. Then we define the cumulative fuzzy histogram

$$\begin{equation} H_{i,j}^n(x) = \sum _{m=a_0}^x h_{i,j}^n(m) \end{equation} \tag{ 7 }$$

and the normalized cumulative fuzzy histogram

$$\begin{equation} C_{i,j}^n(x) = H_{i,j}^n(a_1)^{-1} H_{i,j}^n(x) \end{equation} \tag{ 8 }$$

where a₀, a₁ are minimum and maximum pixel values in $N_{i,j}^n$. Finally, we define the fuzzy equalizing function

$$\begin{equation} f_{i,j}^n(x) = b_0 + (b_1 - b_0) C_{i,j}^n(x) \end{equation} \tag{ 9 }$$

where b₀, b₁ are minimum and maximum output pixel values. This function is, unlike formula (3), different for every pixel a_{i, j}, and the output pixel q_{i, j} is computed according to the formula

$$\begin{equation} q_{i,j} = f_{i,j}^n(a_{i,j}). \end{equation} \tag{ 10 }$$

The use of the fuzzy equalizing function solves the majority of the problems of the classical histogram equalization. However, one serious problem persists, namely the extreme amplification of additive noise in areas with very low contrast, which results in the loss of faint low contrast details. If the neighborhood $N_{i,j}^n$ consists of noise only, the full dynamic range of the output image is used for noise display. The solution is to add more (artificial) additive noise to the input image. Of course, additive noise is not added to the image itself but to the image that is used for histogram computing. This will decrease the contrast of the (natural) original additive noise in the resulting image. This solution is not very practical, because it requires us to work with two copies of the input image and the additive noise must be added to one of them by the noise generator.

Because the original noise in the image and the artificial noise are stochastically independent, the equivalent solution is to compute the convolution of $C_{i,j}^n(x)$ with the probability density function of noise which we would like to add to the image (Hogg et al. 2005). Let us suppose that the added noise has a Gaussian distribution with mean value μ = 0, standard deviation σ and probability density function G_σ(x). Then we define the noise adaptive fuzzy equalizing function

$$\begin{equation} g_{i,j}^n(x) = b_0 + (b_1 - b_0) C_{i,j}^n(x)*G_\sigma (x). \end{equation} \tag{ 11 }$$

This function is used for the creation of image E_{N, σ}(A) in formula (2). The convolution of the normalized cumulative fuzzy histogram with a Gaussian kernel has a significant influence only in the pixel neighborhood dominated by noise in which the image has very low contrast. On the other hand, the influence is negligible in the contrasty parts of the image. The optimal value of σ is usually in the interval 〈2σ_A, 12σ_A〉 where σ_A is the standard deviation of additive noise in the input image A. The higher the value of σ is, the lower is the amount of noise in the low contrast parts of the image. It is not correct to suppress the noise too much because the low contrast details will be lost. The noise must be clearly visible but its intensity must not mask the visibility of fine details. Suitable values for SDO AIA images are found to be about σ = 10 for an extreme enhancement, σ = 20 as a standard value and σ = 50 for low-noise output images (see Figure 1).

Zoom In Zoom Out Reset image size

The program NAFE Image Analyzer, for MS Windows XP and later versions written in Borland Delphi, was developed to implement the NAFE method.¹ As the NAFE method is extremely time consuming, parallel implementation was necessary. Processing of a full (4, 096 × 4, 096 pixels) SDO AIA image takes approximately 5.5 minutes on a PC with Intel Core i7-3820 working on 3.6 GHz (8 processor kernels). However, it is often the case that a region of interest on the Sun is much smaller than the whole image. So, for example, for a 1, 024 × 768 pixels region, the processing takes approximately 18 s. However, the processing time is also dependent on the image structure. Low contrast parts are processed faster than higher contrasted ones.

Matrix L = [l_{k, l}] of the size n × n, n = 129, was defined in the form

$$\begin{equation} l_{k,l} = d \sum _{m=1}^{12} c_m G_{\sigma _m}(k,l), \quad \sigma _m = 2^\frac{m}{2} \end{equation} \tag{ 12 }$$

where d, c₁, c₂, ..., c₁₂ are constants and

$$\begin{equation} G_{\sigma _m}(k,l) = (2\pi \sqrt{\sigma _m})^{-1} \; e^{-0.5(k^2+l^2)\sigma _m^{-2}}. \end{equation} \tag{ 13 }$$

The form of matrix L as a linear combination of Gaussian kernels $G_{\sigma _m}$ was chosen because the constants c₁, c₂, ..., c₁₂ give enough variability which is necessary for fine tuning the properties of the fuzzy neighborhood. The constant

$$\begin{equation} d = 2\pi \left(\sum _{m=1}^{12} c_m \sigma _m^{-1/2} \right)^{-1} \end{equation} \tag{ 14 }$$

ensures that maximum element l_{0, 0} = 1. After extensive testing it was found that the optimum setting is very near to c_m = 1 for all m = 1, 2, ..., 12. Therefore this setting is used as a default in the NAFE Image Analyzer.

The resulting image computed by means of the formula (2) is a linear combination of two images T_γ(A), E_{N, σ}(A) and both of them use the full range of pixel values 〈b₀, b₁〉. However, a linear combination of these two images has a lower dynamic range—usually about 70% of 〈b₀, b₁〉 interval is used for the NAFE weight w = 0.2. Therefore the 16 bits per pixel representation of image B is used (b₀ = 0, b₁ = 65, 535) and the final 8 bits per pixel image is obtained by a linear transform of image B. This transform ensures that the full range of pixel values is used in the final 8 bits per pixel image.

The NAFE method described above was extensively tested on the full wavelength range of 〈93, 335〉 of the SDO AIA level 1 data. The method yields the best results for the 171 Å (Fe ix) images (see Figure 2). It seems to be optimal for this wavelength channel. The method gives very good results for the 193 Å (Fe xii + Fe xxiv), 211 Å (Fe xiv), 304 Å (He ii) images as well. In general, it seems that application of the method for 93 Å, 131 Å, 335 Å images does not yield any meaningful results because of the intense noise in the dark parts of these images. Suppression of this noise also causes the loss of details.

Zoom In Zoom Out Reset image size

The NAFE Image Analyzer software has a built-in automatic setting of all processing parameters for all SDO AIA level 1 images. Therefore knowledge of the underlying mathematical theory presented in this paper is, in principle, not necessary for a meaningful application of the software.

The NAFE method described in this paper is a powerful tool for qualitative studies of the fine structures of the coronal emission at EUV wavelengths as captured in the high resolution SDO AIA images. The method yields artifact-free images in which the local contrast and noise may be well controlled. The main disadvantage of the method is the processing time. Parallel implementation of this method is appropriate. Implementation for MS Windows is available as a freeware.

Solar Dynamics Observatory (SDO) is a NASA project. This work was supported by Grant Agency of Brno University of Technology, project FSI-S-11-3. I highly appreciate suggestions and help from Shadia Habbal.