Tuesday, 10 February 2015

Studying Imagery from the Solar Dynamics Observatory with The Matlab Image Processing Toolbox


In this posting we present an overview if image processing techniques  using the matlab image processing toolbox. We apply these skills to a study of imagery generated by the NASA Solar Dynamics Observatory (SDO). SDO has provided for 5 years now high resolution imagery of our nearest star with a cadence of around 1 second. The video below illustrates the wonders captured by SDO and some of the wildly different eruptive events which can be produced by the sun. These events include solar flares, events where there is an additional ejection of solar material called a coronal mass ejection (CME). The event below also shows more complex moving structures in association with changes in magnetic field lines that loop up into the Sun's atmosphere, the corona. On July 19, 2012, an eruption occurred on the sun that produced all three. A moderately powerful solar flare exploded on the Sun's lower right hand limb, sending out light and radiation. Next came a CME, which shot off to the right out into space. And then, the Sun treated viewers to one of its dazzling magnetic displays -- a phenomenon known as coronal rain.

 The event studied was an M-9.3 flare (just short of an X class, the largest) from an active region right at the Sun's edge (Mar. 12-13, 2014).  The analysis focuses on a group of loops we  examine the properties of the loops over a sequence of times before and after the flare event.

Image processing Tools

 Solar physicists are provided with an armoury of analytical tools, these include SolarSoft (IDL)SunPY (Python)MyExperiment (see references 15 and 16) and the Heliophysics integrated laboratory  (see references 15 and 16).  An excellent overview of solar coronal loops is given in reference 22.

The Matlab image processing toolbox provides a range of functions for
Much of this functionality can take advantage of the performance of graphical processing units for processing imagery. The documentation provided with Matlab and online is excellent, there are many helpful examples. There are also a lot of tutorials which are referenced below (see ref.  1-4). Once data had been imported (from SDO) the  imread function was used to load data into Matlab, it takes as input a filename and returns a handle to the loaded image, a processed image can then be written using the imwrite function. Example syntax is below

imwrite(trans_image, looprotimfile);

There is a range of tools for displaying and exploring images. The simplest is the imshow function. The image tool is particularly useful for anlysing gray scale images it is started using the imtool command. which takes as input the handle to the image. Example syntax is below


imshowpair is used to visualise the difference between two images. There are many steps which can be taken to prepare our selected imagery for analysis. The simplest operations include the ability to select regions of interest, adjust the position and orientation of the image. It is also necessary to modify the colour and contrast within the image so as to enable a more effective study of the features of interest. The geometric transformation tools we used were the imtranslate, imrotate and imcrop functions. These were used to translate rotate and select a region of interest of the image respectively. Each of the functions takes a handle to an input image and returns a handle to the output image. They also take as input the parameters required to translate, rotae and crop the image. Examples are shown below.

imtranslate(li2, [-translation, 0],'FillValues',255);
imcrop(li2, [ xmin ymin width height]);

Matlab provides a range of tools for image enhancement from the contrast adjustment subset we used the histeq function to enhance the image contrast. Before doing histogram equalisation we converted images to grey scale using the rgb2gray function from the graphics toolbox. Basic examples are shown below;


The imsharpen function is a useful tool for bringing out features.  This uses the unsharp masking technique to enhance the image. This may sharpen any edges in the image a Gaussian lowpass filter can be tuned by specifying its diameter in terms of the number of pixels. 

To examine and make measurements from the the imagery we used the imshow and imtool functions. Although we can measure pixel positions and properties with image show the imtool provides tools for measuring the distance between pixels in the image viewer.

The final tools we used for processing were the plot, line and improfile tools. The plot and line graphics functions were used to mark and highlight features of interest. We used the improfile tool to provide pixel-value cross-sections along given line segments. Examples of using these are given below



Matlab also provides the impixel and imhist functions to provide pixel data and present information on the image characteristics. Although we've briefly described a limited subset further details are in the online documentation and in the tutorials  in reference 1-4.

Coronal Loops

Our test case is that of processing a sequence of images of coronal loops in the solar atmosphere. Coronal loops are structures observed in the lower corona and transition region of the solar atmosphere. One of the earliest recorded observations of such structures was by  Angelo Secchi an Italian priest and astronomer. From observation's of the solar eclipses in Spain in 1860  he  presented sketches of coronal loop structures in his 1870 publication Le Soleil: Exposé des Principales Découvertes Modernes (The Sun: Presentation of the Major Modern Discoveries). He organised a further eclipse viewing trip in 1870. One of his sketches is shown below. Secchi's sketch also provides evidence for spicule structures observed on the solar limb. Further images of the solar corona from eclipse viewings are available from reference 12.

The coronal loops are actually spectacular evidence for the magnetic structures in the corona. Loops can extend upto 700Mm in length these structures contain plasma which is heated and rise along the magnetic field lines before falling back again with speeds of up to 100km/s. Loop temperatures can be between 2-3MK and the densities in the range 10^14-10^15 m^-3. The ends of the loops are located in regions of strong magnetic fields near the edge of active regions. Individual loops my last for about a day but a loop system may exist for a few solar rotations. More detailed studies reveal that these structures range in complexity consisting of many strands which can be braided and twisted. With the advent of more sophisticated solar telescopes, coronagraphs and space borne telescopes images of these structures are readily available for example the opening video highlighted coronal rain captured by SDO.

Space based observations of solar phenomena (using missions such as SDO, TRACE etc..) show frequent catastrophic cooling and evacuation of quiescent solar coronal loops over active regions. Analysis of this imagery can show plasma from a temperature of a few million degrees down to less than 100 000 K. Radiative emission from the cooling loops develops initially near the loop tops; later, that cool plasma usually slides down on both sides of the loop.  The evacuation of a coronal loop often occurs after plasma high in the corona has cooled to transition-region or even chromospheric temperatures. The relatively cool material often forms clumps that move at speeds of up to 100 km s−1. The downward acceleration is no more than 80 m s−2, less than a third of the surface gravity (ref. 21). Coronal loops are predominantly heated at the footpoints but their  catastrophic cooling can lead to localized brightenings. Loop heating can lead to an evaporation of plasma into the coronal loop which then cools rapidly due to a loss of thermal equilibrium. The confined region of "condensed'' plasma subsequently falls down under the effect of gravity in the form of a cool plasma blob (ref .17). Thermal nonequilibrium may nonetheless play an important role in prominences and catastrophic cooling events (e.g., coronal rain) that occupy a small fraction of the coronal volume ( ref. 20).


Our study focuses on an M2.5 class flare event which occurred on 12th March 2014 (see NOAA), the event was recorded in the SDO event gallery. The M2.5   flare occurred in NOAA group number 11996  located at N14W76. The region we study is NOAA group 12003 located at N06W37, no flares were reported for this group on this day. The Atmospheric Imaging Assembly (AIA) aboard the Solar Dynamics Observatory has state-of-the-art spatial resolution and shows the most detailed images of coronal loops ever observed. The series of coronal filters peak at different temperatures, which span the range of active regions. These features represent a significant improvement over earlier coronal imagers. There are many sources for gathering and exploring solar data including.

SDO AIA/HMI data access
Virtual solar observatory (VSO)

For our study we used the SDO data browser. We gathered two sets of data corresponding to the AIA193 telescope and the AIA171 telescope and including the output from the Helioseismic and magnetic imager. For the AIA193 data set we collected 184 images starting at 00:40:07 on 12th March 2014 and finishing at 22:40:07 in 13th March 2014. For the AIA171 (with HMI data) we collected 121 images starting at 00:29:42 on 12th March 2014 and finishing at 22:43:11 on 13th march 2014. For both sets we collected images with a resolution of 4096x4096 pixels. It would have been helpful to gather data using the VSO and to gather FITS data froom SDO. Using the flexible image transport system (FITS) tool box available with Matlab it would be straightforward to extract and  use SDO metadata in our analysis.  The event may also be viewed using helioviewer (2014-03-12T22:36:37) this was an easy way to view the metadata from SDO.

A video of the region we studied generated using helioviewer is shown below for the AIA193 detector.

The videos below were generated from the sequence of images used in our analysis the region of interest is just above the centre of the image. The brownish video below is that for AIA193, the bright flash is the tell tale sign of a flare. The brightness of the flare causes very bright saturation and blooming above and below the flare region on the CCD detector and caused extended diffraction patterns to spread out across the SDO imager. The next (yellowish) video is for the AIA171, the red/green mottling shows the magnetic field strength and polarity from measured using the HMI.

The steps followed in our analysis procedure are as follows

Step 1.  Image translation

From our video generated using the output from helioviewer we observe the features moving from left to right as the sun rotates. To perform our analysis we translate the images to correct for this rotation so that features of interest are at a fixed location within each image. From the relatively short  helioviewer movie it appears that the rotation is constant for all locations on the solar disc. Since the sun is a gaseous plasma there is a variation of this rotation with the solar latitude. This is given by the relation.

  A= 14.713 ± 0.0491 °/d
B= -2.396 ± 0.188 °/d
C= -1.787 ± 0.253 °/d

phi is the solar latitude and omega is the solar rotation. Since the solar radius is 695.5Mm and the solar disc in our images is 3304 pixels in diameter, each pixel corresponds to 0.421Mm. Using the solar rotation the translation is approximately
\Delta is the resulting shift in pixels, P_s is the pixel sclae (0.421Mm) \omega is the solar ration in radian/s. t_i and t_0 are the times of the observation and the initial observation respectively. Aschwanden 2011 (ref. 6 ) provides details of this correction in their paper. Our implementation of this correction is in the Matlab script file translatelopsequence.m which we applied to the full sequence of images

Step 2. -Region selection and image cropping

Our region of study was NOAA group 12003 located at N06W37, no flares were reported for this group on this day. We set the challenge of seeing if we could detect changes in the group which occurred as a result of the flare event at NOAA group number 11996  located at N14W76. Inspection of group 12003 indicated that it may be sufficiently near and it had a number of features which we sufficiently stable and well defined to measure.  With the script  croploopsequence.m we applied the Matlab imcrop function to the sequence of translated images.

Step 4 - loop examination and comparison

Tracking the changes in a particular coronal loop is actually quite challenging. The eye is good at picking out features as they change within a sequence such as that within a video. However when we study individual frames and wish to pick out and detect the evolution of an individual loop this is challenging. The process of identifying the loops was achieved by altering the image contrast using the matlab histeq function, for each image set. We compared these equalised images with the original cropped image, with a check against the video frames and by marking off points on earlier and later frames in the image sequence. It was necessary to undertake this process iteratively sometimes returning to correct earlier time steps. In this way there is greater confidence that we are tracking a particular loop.   

For our preliminary examination of the loops we prepared the comparepairs.m Matlab script. This was used to experiment with some of the functions within the image processing toolbox and computer vision toolbox. This included tests of the image registration functions, image registration is the process of aligning images of the same scene. We also experimented with some of the edge detection techniques. These tests were undertaken using the processloopsequence.m, detectfeatures.m and houghlines1.m scripts. Given the dynamic nature of our coronal loop images it was quickly realised that automating a procedure for isolating and measuring specific loops according to a prescription is challenging. It was decided to manually measure a subset of the images. We selected 13 images centered around the flare event. 

For the AIA193 data set the event occured at 22:40:43 on 20140312, this corresponded to image 89 in the overall data set. Our subset was the following images
The fifth image in this subset corresponded to the flare event.
These images correspond to the following time in minutes
20.3   50.7   80.7   95.3  110.7  125.8  140.1  154.7 171.0  198.7  214.5  231.1  261.3

For the AIA171/HMI data set the image nearest to the event was taken at 22:44:23 on 20140312, this corresponded to image 54 in the overall data set. Our subset was the following images
The fifth image in this subset corresponded to the flare event.
These images correspond to the following time in minutes
24.4   54.0   84.4   99.4  114.4  129.4  144.4  159.4 174.4  202.4  217.4  232.4  262.2

For the indices above note that we did not select all the images for the AIA193 we selected those which were closest in time stamp to the AIA171 images

step 5 - loop measurements

For this study we measured the length, height and width of a well defined loop within the group and which remained sufficiently define over our selected  sequence of 13 images. We also measured the loop intensity. The measured parameters are illustrated by the diagrams a,b and c shown below.

We used the  matlab histeq function and balanced the image contrast, it was then easier to identify and measure loop features. The images show curves plotted on the selected loop elements after the initial loop observations, these were the lines which were used to measure the loop intensities. The most effective method is to click and select points from features, for example we clicked points along the flux loop. For this step we used the plotimpair.m script, the measurements were  stored in the loopproperties_img2.m, for use by the later image processing routines. In order to ensure we were measuring the same and correct loop this process was repeated a number of times.

For the image below we have selected the image nearest to the flare event. For the AIA193 we can see clearly the flare event. The brightness of the flare causes saturation and blooming above and below the flare region on the CCD detector it also caused extended diffraction patterns to spread out across the SDO imager. The AIA171 iamge also includes the mottled patches from the HMI corresponding to the strength and polarity of the magnetic field.

The results from this measurement process are shown below. In the following images its important to note that the flare event occured at 110.7mins (for AIA193) and 114.4  (for AIA171)
Variation of Loop Height with Time (Flare event at 110.7 minutes)

Variation of Loop Length with Time (Flare event at 110.7 minutes)

Variation of Loop Width with Time (Flare event at 110.7 minutes)

In the above images it is encouraging that the measurements agree to within the errors we expect for this measurement process. There is a couple of points which are questionable indicating the possibility of incorrect lop/feature selection.

Step 6 Loop thermal analysis

The ability to measure and study thermal characteristics and the thermodynamics of features in the solar atmosphere is a highly valuable tool, this  has developed increasingly since the availability of multispectral, high resolution and high cadence solar imagery. In our final analysis step we used the technique of isothermal analysis to determine the temperature of the coronal loops. This analysis of the loops was computed using the matlab script, loopintensitycomputation.m. In the previous section we had attempted to identify and highlight the full profile of a particular coronal loop. Very often a section of the loop would overlap or intersect with a different loop. For the thermal analysis we used only well defined sections of the loops for the analysis. We edited and stored these sections in the script file loopproperties_partialloops_img2.m. Using these sections of loop we compute the pixel profiles for the gray scale image generated using the the Matlab rgb2gray function. The pixel profile was computed using the improfile function. Using our loop sections as input for the  improfile function we can use the resulting image profile to compute a pixel data number (DN), representing the intensity of the emission from the loop. We describe the procedure for computing the DN values below.

The result of this analysis lead to the following image, here it is important to note that the flare event occured at 110.7mins (for AIA193) and 114.4  (for AIA171)
Measured Loop Intensity and its Variation with Time for the AIA171 and AIA193 Images
In the above image we see that the loop intensity falls initially but that there is rise in intensity occurring after the flare event. To determine if this is at all significant we see if these values give temperature values which are consistent with typical solar coronal loops. To make sense of this it is necessary to recognise that the SDO AIA cameras uses different filters which have varying and known sensitives to radiative emissions of different wavelengths. AIA has seven EUV and three UV–visible channels. Six of the EUV channels are centered on a different spectral line of ionized iron and allow the construction of relatively narrowband temperature maps of the solar corona from 0.4 to 20 million degrees. The plot below (see ref. 12) shows the response curves for the different filters on AIA, these curves show how the image intensity varies with temperature. It can be seen how the different filters have different sensitivities for different temperatures. Using this information it is possible to use our measured intensity data to estimate temperature values.

Response Curves for the SDO AIA Cameras, Responses are Shown for the Different Filters (See ref. 12)
In the figure above the fluxes are given in units of data numbers per second (DN s−1) these counts  are normalised by detector the exposure time, which is different in each filter. Since different spectral lines have varying sensitivies for different temperatures this fact can be used to identify characteristics for different regions (i.e. different heights)  within the solar atmosphere (e.g. see solar facts1). The temperatures for the main spectral lines are given in the table below, the full information of these spectral lines for astrophysical plasmas is provided by the Chianti database

Here is the technique for computing the DN's from ref 12. For each of the loop sections we averaged the loop- and background-pixel values and calculated standard deviations.  A background region is adjacent to the loop under consideration. We subtracted the averaged background results froom the averaged pixel counts for each loop section. Averages were taken over the loop length or the background section. The background-subtracted averages and uncertainties were then normalized by the exposure time for the appropriate filter, this results in units of Data Numbers s−1. Fluxes are now given in units of data numbers per second (DN s−1) after normalizing the counts (DN) by the exposure time, which is different in each filter. Aschwanden (ref 6.)  uses spline fitting for flux profile averaging. For each selected loop they mark the start, midpoint, and end position of the loop spine and interpolate a 2D-spline function with steps of one pixel in order to obtain a smooth curvature. Our implementation is presented in the script file loopintensitycomputation.m the AIA filter characteristics  are stored in loopproperties_partialloops_img2.m. For an isothermal analysis of flux loops we can use the following relation

Flux_171 and Flux_193 are the measured fluxes for the 171 and 193 filter respectively. Resp_171 and Resp_193 are the detector responses for a given temperature for the 171 and 193 filter respectively. Armed with this we can use this ratio to estimate the isothermal loop temperature. Using the detector response values the response ratios for each temperature are shown in the plot below.
Isothermal analysis for the 171-to-193 ratio.
Using the ratio of the fluxes for each of our loops we can read the loop temperature from the above plot (since we know flux ration is the same as the response ratio). This procedure enabled us to produce the following temperature profile. In the resulting image, below, we note that the flare event occurred at 110.7mins (for AIA193) and 114.4  (for AIA171). There is a hint of a temperature rise following the event.
Our approach is a simplistic preliminary attempt and we refer the reader to the excellent work and reviews in references 7,8,9,12,13,22. In this area there is a controversy surrounding the challenge of computing the cross-field temperature distribution for a loop. Also, are loops isothermal or multithermal? Is the observed loop a single flux tube or a collection of tangled magnetic strands? Here, a loop is a discernable structure in an observation and a strand is a fundamental flux tube where the temperature and density are constant perpendicular to the magnetic field. A loop could be composed of a single strand or many unresolved strands. Resolving this controversy has important implications for the coronal heating problem. Results from isothermal analysis, which is the average of the true plasma distribution weighted by the instrument response functions, appears to be deceptively low. These results have potentially serious implications and show substantially smaller temperature gradients than predicted by standard models for loops in hydrodynamic equilibrium these have been used as strong evidence in support of footpoint heating models.

Using the technique of differential emission analysis  has shown that the loop plasma is indeed multithermal and not isothermal. The Differential Emission Measure (DEM) ξ(Te) is a measure of the amount of emitting plasma along the Line Of Sight (LOS) as a function of the electron temperature Te. The intrinsic difficulties involved in this inverse problem lead to many complications in its inference, making its interpretation ambiguous (see reference 12 and 13).


We have demonstrated the application of the matlab image processing toolbox for studying imagery of coronal loops produce using SDO AIA camera.

We studied a loop system in a region close to an M2,5 class flare. Results appear to indicate an increase in the loop intensity after the flare event. This may not be surprising given the  amount of energy liberated by the flare event a and the possibility of radiative heating of our coronal loops.

A study of a single group is a provider of an indication, further results would be required to clarify this conclusion. Measurements of loop systems well isolated from flare events would be required, also necessary is a measurement of loops near and at varying distances from flare events.

It is intriguing that we observe the initial intensity peak for the 171 and 193 these peaks are separated by almost 49 minutes. It is a risk that we have selected a different loop although positional data and the iterative process indicates we've reduced this risk. 

Although we have established a procedure using an iterative approach to make these measurements it is a time consuming process. It is difficult to apply the Matlab visualisation toolbox methods without significant effort and development. A more successful approach might be that of a citizen science approach for example see the zooniverse sunspot project or the zooniverse solarstorm watch. This method beneficial when computer visualisation algorithms find it difficult to pick out and match the details in the images we have studied. Approaches used could be machine learning such as a neural network presented with a training set of coronal loop imagery, the training set includes a sequence of cropped  (and suitably translated images) along with a desired output.

We refer readers to reference 22 for a detailed review of the physics of coronal loops.
Key issues for the future remain the loop fine structure and dynamics. We need to address what are the ultimate elementary loop components, and whether they are unique or determined by local conditions, i.e., what determines the section scale size. We need the highest possible spatial resolution, probably in different bands. Further investigation of the plasma fine thermal structure and dynamics requires also high spectral resolution. High resolution broad-band X-ray spectroscopy is foreseeable to probe the hot components, signature of impulsive heating. Also the investigation of temporal variations still deserves attention. If loops are really so dynamic and subject to a distribution of heating events, whatever it is, the light curves are very difficult to interpret and signatures of any possible small-scale events are confused in a storming activity and by the plasma inertia.


We have used the matlab image processing toolbox to process a sequence of images of solar coronal loops taken by SDO. In hindsight it would have been better to use the VSO, the heliophysics integrated laboratory or myexperiment (ref. 15 and 16) the sunpy or or the SSW. Using this we can extract the image metadata. The matlab toolbox also provides FITS toolbox, this provides highlevel functions such as fitsread and fitsdisp which can be used to read a FITS data file and the metadata respectively.


  1. Getting Started with Image Processing Toolbox-2014b
  2. Image processing examples - matlab central file exchange
  3. Image Processing with MATLAB-Pascal Getreuer
  4. MATLAB 6.5 Image Processing Toolbox Tutorial
  5. Automated Temperature and Emission Measure Analysis of Coronal Loops and Active Regions Observed with AIA/SDO Markus J. Aschwanden, Solar Physics 2011
  6. Solar Corona Loop Studies with AIA: I. Cross-Sectional Temperature Structure Markus J. Aschwanden, Paul Boerner, arXiv.org > astro-ph > arXiv:1103.0228, (ADS 2011ApJ...732...81A)
  7. Coronal Loop Observations and Diagnostics
  8. Evidence for Nonuniform Heating of Coronal Loops Inferred from Multithread Modeling of TRACE Data Aschwanden et. al. 2000ApJ...541.1059A
  9. Time Variability of the ``Quiet'' Sun Observed with TRACE. II. Physical Parameters, Temperature Evolution, and Energetics of Extreme-Ultraviolet Nanoflares Aschwanden et. al.2000ApJ...535.1047A
  10. Solar corona images observed using solar eclipses
  11. A Multi-Wavelength View on Coronal Rain Muller 2004 (see also ADS Muller  20052005ESASP.600E..30M  )
  12. Isothermal and Multithermal Analysis of Coronal Loops Observed with AIA, Schmelz, 2011ApJ...731...49S
  13. On the Accuracy of the Differential Emission Measure Diagnostics of Solar Plasmas. Application to AIA/SDO. Part II: Multithermal plasmas    C. Guennou et al 2012 arXiv:1210.2302 [astro-ph.SR]
  14. JSOC SDO AIA Metdata keywords
  15. Workflows for Heliophysics, Journal of Grid Computing, Anja Le Blanc et al, September 2013, Volume 11, Issue 3, pp 481-503
  16.   PIERANTONI, G. et al. A WORKFLOW-ORIENTED APPROACH TO PROPAGATION MODELS IN HELIOPHYSICS. Computer Science, North America, 15, jul. 2014. Available at: <https://journals.agh.edu.pl/csci/article/view/882>
  17. High Speed Coronal Rain Muller et al. 2005A&A...436.1067M
  18. E. O'Shea, D. Banerjee, and J. G. Doyle,  Plasma condensation in coronal loops,  A&A 475, L25-L28 (2007) , DOI: 10.1051/0004-6361:20078617
  19. Ionson, J. A., Resonant absorption of Alfvenic surface waves and the heating of solar coronal loops, 1978ApJ...226..650I
  20. James A. Klimchuk et al., Can Thermal Nonequilibrium Explain Coronal Loops? 2010 ApJ 714 1239
  21. Schrijver, CarolusJ., Catastrophic cooling and high-speed downflow in quiescent solar coronal loops observed with TRACE, Solar Physics, February 2001, Volume 198, Issue 2, pp 325-345, DOI 10.1023/A:1005211925515
  22. Fabio Reale, Coronal Loops: Observations and Modeling of Confined Plasma, Living Rev. Solar Phys. 11 (2014), 4.  doi: 10.12942/lrsp-2014-4


  1. Hi,,

    Nice representation

    Thank you for submission.

    can you tell us : How long will PV modules last?

    Answer: Based on manufacturers' in-field experience and reliability testing, PV modules will probably last longer, and are more reliable than just about any other capital investment for your business.

    for more info contact us:

    Solar Inverter supplier in delhi