A spectral optical flow method for determining velocities from digital imagery

by   Neal Hurlburt, et al.

We present a method for determining surface flows from solar images based upon optical flow techniques. We apply the method to sets of images obtained by a variety of solar imagers to assess its performance. The opflow3d procedure is shown to extract accurate velocity estimates when provided perfect test data and quickly generates results consistent with completely distinct methods when applied on global scales. We also validate it in detail by comparing it to an established method when applied to high-resolution datasets and find that it provides comparable results without the need to tune, filter or otherwise preprocess the images before its application.



There are no comments yet.


page 4

page 9


Optical Flow on Evolving Surfaces with an Application to the Analysis of 4D Microscopy Data

We extend the concept of optical flow to a dynamic non-Euclidean setting...

Optical Flow on Evolving Surfaces with Space and Time Regularisation

We extend the concept of optical flow with spatiotemporal regularisation...

Velocity variations at Columbia Glacier captured by particle filtering of oblique time-lapse images

We develop a probabilistic method for tracking glacier surface motion ba...

Solar Image Restoration with the Cycle-GAN Based on Multi-Fractal Properties of Texture Features

Texture is one of the most obvious characteristics in solar images and i...

Attacking Optical Flow

Deep neural nets achieve state-of-the-art performance on the problem of ...

Where computer vision can aid physics: dynamic cloud motion forecasting from satellite images

This paper describes a new algorithm for solar energy forecasting from a...

Multidimensional Digital Filters for Point-Target Detection in Cluttered Infrared Scenes

A 3-D spatiotemporal prediction-error filter (PEF), is used to enhance f...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The most powerful events observed in the solar system are the result of convective dynamics in the outer layers of the Sun. Researchers trying to understand these dynamics need tools to measure and study these turbulent motions and their interaction with the local magnetic field. Several methods for deducing the flows in the solar photosphere (the visible surface of the sun) have been developed over the years. Many of these have been based upon some form of local correlation tracking (LCT) TTS86 ; N87 where pairs of successive images of the photosphere are broken into sub-images which are then shifted relative to each other to find an optimal relative shift. This shift is then associated with the local mean velocity of the flow. These methods appear to work well on high-resolution images collected from both ground-based B_al88 and space-based S_al88 ; TTS89 observatories. There have been a few attempts to assess and compare the performance of such methods over the past twenty years including those of Hurlburt et al. H95 , Welsch et al.W07 , Chae and Sakurai CS2008 and most recently, Verma, Steffen and Denker Verma13 . These have used a variety of data input for their assessments, from those derived from simulations of MHD flows in the photosphere, to creating image sets from solar images by applying known distortions, to direct comparisons with real solar data. All tested methods gave consistent results. However the best method was somewhat dependent on the test and the choice of the respective parameters for each method.

Alternative methods have been developed in other fields with similar goals. Machine vision researchers developed optical flow methods for deducing the relative motions of objects in digital images B2011 ; S2014 and atmospheric scientists developed various methods for deducing winds on Earth from cloud motions in satellite images L2013 ; C1999 . Experimentalists in other branches of fluid dynamics, including blood flow measurements deduced from x-ray images N2012 and flows with various tracer particles in suspensions W2001 , have explored similar methods. One difference between the typical machine vision problem and that of deducing fluid velocities is that the former seek motions of discrete objects while the latter seeks motions of continuous flows. Applying optical flow methods such as those described by Jähne J93 on high-resolution solar images typically underestimate the velocity of the flows. In part this is due to the simple spatial averaging used in deriving the velocity.

Hurlburt H99 presented a method which does not make use of such spatial averaging. Instead the flows are assumed to be smooth and continuous – being represented by a truncated Fourier series. Here we present a detailed description of the method and assess it using previously developed tests and compare it to local correlation tracking methods.

2 Method

The basic assumption is that the visible pattern observed in the fluid, as measured by the local intensity will be advected by the velocity field and hence should satisfy the equation


In the case of solar physics the pattern is typically formed by convective motions in the photosphere, which are clearly visible in white-light images and which appear to be advected by larger scale flow fields. Images are collected frequently relative to the flow speeds, such that the displacement of the pattern between any sequential pair of images is less than a pixel. The problem is to determine from a time sequence of two dimensional images in the presence of measurement noise and other “noise” sources, such as the acoustic oscillations present in the solar atmosphere or missing frames due to data dropouts. Using equation (1) we can seek the best fit velocity field using least squares. First we form the merit function of the fit for the full dataset .


which we seek to minimize. Here the sums are taken over the discrete pixel and frame coordinates for and If the velocity field is a continuous field, we can express the fit velocity as a Fourier series


where and

are unit vectors,

and are complex amplitudes and and are the number of Fourier modes retained in the expansion. Substituting this into equation (2) and differentiating with respect to and results in the system of equations


This can be reorganized to form the complex system of equations




These matrices consist of discrete Fourier transforms of the various, time-averaged products of the spatial and temporal derivatives of the image. The matrices of the complex linear system (

6) for the spectral amplitudes , can be combined to form a single hermitian matrix of size

This method has been implemented in an IDL 111Trademark, Exelis. routine opflow3d and is available as part of the SolarSoft environment Freeland . The time derivative is evaluated using finite differencing between sequential images while the spatial derivatives is evaluated using 4th order finite differences on the average of the two images used for the time derivative. The matrices are then computed for the entire time-space cube and the system is solved.

Solving the system using direct methods can quickly become expensive, scaling as . The method requires 70 seconds on an 2013-vintage iMac workstation (with 3.4GHz processor and 32GB of memory) to obtain the solutions for on a image. This is partially offset by the fact that the method requires no preprocessing or filtering and that it fits flows over many instances in time in one go. The matrices in (7) can also be reused in subsequent calculations with minimal additional expense. Performance could be further optimized by taking advantage of the matrix structure in equation (7) which has a blocked-Toeplitz-Toeplitz-block (BTTB) structure C1999 .

3 Evaluation and comparisons

Figure 1: Comparison between a derived velocity field and input field used to distort a solar image (background). The input (derived) velocity field is displayed with black (white) arrows whose areas are proportional to the magnitude of the local velocity. The relative scale of the white arrows has been reduced slightly for aid in comparison.

As a first test of the method, we take the simulated observations developed by Hurlburt et al. H95 . Using a sixth-order accurate numerical scheme HR00 they took a single intensity image of solar granulation and evolved it with equation (1) with a steady velocity field. They then degraded and resampled the image to represent the expected resolution of the Michelson Doppler Imager (MDI) on the Solar and Heliospheric Observatory MDI95 . Since there is no source of noise and the imposed flows are themselves based on Fourier modes, we expect and observe that the opflow3d method can recover the flow with a high accuracy. The results for a case where on a pixel image is displayed in figure 1, along with the known input velocty field and a sample image. The two sets of arrows, corresponding to the known (black) and derived (white) velocity fields are almost perfectly correlated, both in direction and magnitude.

Figure 2: The differential rotation of the sun as determined by applying the spectral optical flow technique on one hour of full-disk MDI continuum images (solid) and published best fit from Doppler measurements (dotted). The former was calculated with =4, using sixty pixel images. The two agree within the noise level of the supergranular flow field.

With this basic validation of the method on perfect data, we turn to comparisons with other methods and operate on real solar data and then consider a detailed error analysis. First we compare the results of applying this method to other large-scale, full-disk measurements. Figure 2 displays the zonal (E-W) component of the solar velocity along the central meridian of the Sun as a function of solar latitude derived from one hour of MDI MDI95 data. We include the corresponding measure based upon fits to Doppler measurements S92 . Aside from the departures induced by sampling errors of the supergranular flows in this short time, the two curves agree very well.

4 Error Assessment

There are several factors that may contribute to errors in the velocity estimate provided by opflow3d. These can be broken into three classes: systematic errors introduced by the choice of velocity representation, errors due to image quality and errors introduced by physical effects in the solar atmosphere. As a first step, we take what is currently the most consistent and stable images available, using data obtained by the Helioseismic and Magnetic Imager (HMI) on SDO Scherrer12 . We then subject them to a variety of controlled tests to address the first two classes of error. This is the best-case scenario for studies of the solar photosphere: a stable imager observing ”quiet” sun (where magnetic effects are negligible). The following section considers a more complex situation of observing magnetic regions with a less-stable imager.

One hour of Level-1 HMI continuum images consisting of 80 individual frames were used for this study beginning at 2010-10-26T08:29:00. A set of 1024x1024 pixel sub-images were extracted centered on a coronal hole identified in the Heliophysics Event Knowledgebase Hurlburt12 222SOL2010-10-25T23:00:08L032C113 (herein referred to as C2010).

Figure 3: The horizontal velocity profile through the middle of the field of view. The magnitude of the simulated flow is displayed in black and the fitted velocity , and error, in red and green respectively. It is clear that error drops rapidly away from the boundaries, within half a wavelength of the truncated mode.

The use of a truncated Fourier representation for the velocity field is a common practice in fluid dynamic investigations. However it has two well-known issues that must be considered: it imposes a periodic structure to the flows and may introduce aliasing or other errors due to truncation. To assess these effects, we use the first frame from the C2010 and generate ten artificial images from it using equation (1) with a known, hexagonal flow field, in a similar approach to that of Hurlburt H95 .

Our implementation uses a fast Fourier transform and can be sensitive to discontinuities at frame boundaries. While the Fourier representation forces the velocity to match across opposing boundaries, the application of a least squares fit works to confine such effects to near the edges. The resulting spikes in the residual decay rapidly away from frame boundaries (Figure 3). Thus, avoiding pixels near the boundaries is the practical means for avoiding these errors.

Aspects of image quality that may influence our fit include large displacements arising from insufficient sampling rates (which may introduce ambiguity into the possible solution), image noise and missing data. To assess the impact of image quality, we adapt the approach from the previous section. The first frame of the C2010 is advected to generate ten frames with a known velocity field composed of 22 by 22 randomized Fourier modes. A set of these data cubes is generated with a range of magnitudes for . We find the accuracy of the fit starts to break down when the RMS velocity exceeds 0.4 pixels/frame. This result applies to the images with repeating, high-frequency patterns such as solar granulation seen in C2010. In contrast, larger-scale patterns could allow larger unambiguous motions between frames.

The effect of truncation can be assessed by evaluating this same test case with an increasingly large number of modes () in equation (3), from to . We find that the overall trend is constant, but with higher values being more sensitive to the effects of noise as the effective sample size of the fit decreases. Thus the selection of the number of modes should take this trade off into account. A rule of thumb would be that where is the number of features (e.g. granules) required to span the image.

To assess the impact of noise on velocity estimate, we first measure the inherent noise in the synthetic datacube used above by setting in equation (2) to give counts per frame.

is used here to represent the original variation that is reduced by fitting V, and because it scales with image contrast. Next, we generate a sequence of datacubes by adding increasing levels of Gaussian white noise. The datacube with added noise standard deviation

has velocity estimate . By comparing with , we can gain some insight into the sensitivity of the method to noise. We find that for a maximum relative error of 1%, the maximum additional noise must have . Comparing this value of to the inherent signal noise , demonstrates that high noise images can still yield reliable velocity estimates. This robustness against measurement noise most likely results from averaging over many pixels.

Finally there are features in the solar atmosphere that may impact the performance of any method of velocity estimation. These include the presence of strong acoustic modes (known as five-minute oscillations) which generate a relatively-smooth, but random intensity fluctuation in solar images; Limb-darkening, which introduces fixed, large-scale intensity gradients due to line-of-sight effects; and strong magnetic features that may distort the intensity patterns in non-obvious ways. In exploring these cases, we cannot compare our results to a known solution. Instead we must make statistical inferences.

If we seek flows that persist on time scales significantly larger than five minutes, we can assess the effect of acoustic oscillations. We take to be the velocity fit for frames starting at frame . We examine the convergence of , as increases. Using C2010, we divide the 80 frames into sets 1 and 41, consisting of the first and last sets of 40 frames. Since the images in this case may possess an overall motion akin to camera motion (say ), we first subtract such motions from to produce the the Euclidean metric pixels/frame. Similarly, the distance between and is the Euclidean metric ppf. This distance is a rough measure of precision. Thus, with , considerable convergence is apparent.

To examine the rate of convergence, we can next estimate for . We then compute the distance . This distance, averaged over all , is 0.19 ppf. Clearly a two frame estimate is poor. Comparison of the 2 frame and 40 frame distances, 0.19 ppf and 0.019 ppf, imply that the rate of convergence is between and , which is expected for traveling wave patterns in acoustic oscillations. This is consistent with previous studies of solar flows using LCT methods N87 .

The contrast in the image varies smoothly across each frame in our sample due to the limb-darkening effect of the solar atmosphere. Such gradients can cause problems for some methods. However in our case, it only changes the weights of the sum of the squares in equation (2), so that the fit is only slightly affected. The same logic shows the method is insensitive to image-quality problems such as missing frames.

5 Comparison to other methods

Figure 4: A comparison of results of the method used by (a) V&D (from their Figure 2c) and (b) opflow3d for the same sample of Hinode/SOT data shows detailed agreement. Here we display the magnitude of the two velocity fields using approximately the same color map and scaling from black/dark blue to yellow/red. The flow velocities exhibit the same pattern of outward moat flows around sunspots and inflows around plage.

As a final test, we provide a detailed comparison of our method to that used a recent study by Verma & Decker Verma11 (hereinafter referred to as V&D). They conducted a thorough investigation of horizontal flow fields observed in Hinode G-band images using a local correlation tracking (LCT) approach that was used in November and Simon N87 .

Following V&D, we selected an hour-long set of G-band images collected by the Hinode Solar Optical Telescope (SOT, T08 ) on June 4, 2007 between 14:27 and 15:27UT. In that study, the authors first applied the standard calibrations to the images and then further pre-processed them by correcting for foreshortening, applying a rigid alignment between the images to remove spacecraft jitter and solar rotation, and then employed a subsonic filter to remove acoustic oscillations.

We also calibrate the images using the SolarSoft routine fg_prep with its default settings. However, we do not apply the other preprocessing steps. Other than foreshortening, those corrections effectively remove noise from the velocity signal that we are seeking, be it jitter from the spacecraft, bulk motion across the field of view or distracting intensity fluctuations. Since the opflow3d method has already been shown to address such noise sources, we rely on it alone to do so. In addition, we found seven of the 238 images in the sequence were missing: rather than attempt to correct for these, we left those images blank and left it to opflow3d to deal with the consequences.

Figures 4 and 5 display a comparison of applying the two methods to the same image set. The only free parameters for opflow3d are the number of modes used to fit the velocity field, and whether to use a direct or iterative solution method: we select a direct solution with 20 modes in each direction. This corresponds to an effective pixel size of about 2Mm when compared to the Gaussian FWHM used by V&D.

Figure 5: Histograms of the speeds found by (a) V&D (from their Figure 6) and (b) opflow3d for the same one-hour sample of Hinode/SOT data also show detailed agreement. The solid curves display the normalized histogram of the velocities for the two methods. (The other curves in (a) correspond to evaluations over longer time intervals from 2 to 16 hours.)

We correct for foreshortening after the fact by scaling the components of the velocity and display the magnitude of the resulting velocity field (less the average velocity over the frame) in Figure 4. Strong moat flows are visible in the lower left, as well as converging flows elsewhere in plage areas. Figure 5 displays a normalized histogram of flow speeds, which can be compared to figure 6 of V&D. In both cases the peak value across the field of view is around 0.3 km/s. We find the overall rms speed to be 0.46km/s, the median to be 0.42 km/s and the maximum to be 1.68; as compared to 0.44 km/s, 0.40 km/s and 1.95 km/s respectively. The fact that opflow3d

method retains slightly higher rms velocities while reducing the extremes suggests that it might both retain a higher resolving power while mitigating the influences of outliers.

6 Discussion

We have described a method for deriving flows from sets of images obtained by a variety of solar imagers. The opflow3d procedure has been shown to extract accurate velocity estimates when provided perfect test data and quickly generates results consistent with completely distinct methods when applied on global scales. We have also verified that it agrees in detail with an established method when applied to high-resolution datasets – and without the need to tune, filter or otherwise preprocess the images before its application. It is currently running as part of the HEK system to identify regions of solar eruptions Hurlburt15 from data collected by the Atmospheric Imaging Array on SDO Lemen12 .

Our method has been found to work well on other types of image data, including magnetograms, since the only assumptions made are that the motions displayed in them are reasonably smooth and persistent. It can also be combined with other image processing methods to extract motions of specific features within the field. For instance the motion of the two polarities (North/South) in magnetograms could be tracked by thresholding the images prior to using opflow3d. Similarly, particular scales could be extracted by using high- or low-pass filters.

With the basic approach established, there are several avenues for improvement. First, we could replace the model equation (1) with a more elaborate one, say one that solves the vertical component of the induction equation to extract velocities from sets of magnetograms. Second, we could provide a more elaborate fitting function, say one that permits a simple time dependence. Finally one can seek to optimize the method using more sophisticated tools of linear algebra. We will explore some of these options in future work.

We are grateful to K. Schrijver, R. Shine, T. Tarbell, A. Title and T. Berger for their helpful conversations and comments. This work has been supported by the National Aeronautics and Space Administration through contracts NAS8-39746 and NAS5-3077 and by Lockheed-Martin Independent Research Funds.


  • (1) Baker, S., Scharstein, D., Lewis, J., Roth, S, Black, M. and Szeliski, R., A Database and Evaluation Methodology for Optical Flow , Int J Comput Vis 92, 1 31 (2011)
  • (2) Brandt, P. N., Scharmer, G. B., Ferguson, S., Shine, R. A., Tarbell, T. D., and Title, A. M., “Vortex flow in the solar photosphere,” Nature, 335, 238 (1988)
  • (3) Cochran, W., Plemmons, R. and Torgersen, T., Exploiting Toeplitz structure in atmospheric image restoration in Structured Matrices in Mathematics, Computer Science, and Engineering II , ed by V. Olshevsky, Proc AMS-IMS-SIAM (1999)
  • (4) Chae and Sakurai, 2008, ApJ 689, 593, ”A Test of Three Optical Flow Technique - LCT, DAVE, and NAVE”
  • (5) Freeland, S., Bentley, R., and Murdin, P.: 2000, Encyclopedia of Astronomy and Astrophysics.
  • (6) Hurlburt, N.E., Schrijver, C.L., Shine, R.A. and Title, A.M. “Simulated MDI Observations of Convection,” in The 4th SoHO Workshop: Helioseismology, (ESA: Noordwijk, 1995 ) 2,239.
  • (7) Hurlburt, N. E. 1999, SOHO-9 Workshop on Helioseismic Diagnostics of Solar Convection and Activity
  • (8) Hurlburt, N., 2015 Space Weather and Space Climate, submitted, arXiv:1504.03395 [astro-ph.SR]
  • (9) Hurlburt, N. E., and Rucklidge, A. M. 2000, MNRS, 314, 793
  • (10) Hurlburt, N., Cheung, M., Schrijver, C., et al. 2012, Solar Phys., 275, 67
  • (11) Jähne, B., Digital Image Processing, (Springer-Verlag:Berlin, 1993)
  • (12) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Solar Phys., 275, 17
  • (13) Lou, T., Lin, L. and Zhan, N., Discussion Measurement Models and Algorithms of the Wind Vector Field Based on Satellite Images , Applied Mathematics, 4, 122-126 , (2013)
  • (14) Negahdar, M. and Amini, A., A 3D optical flow technique based on mass conservation for deformable motion estimation from 4-D CT images of the lung , Proc. SPIE 8317, Medical Imaging 2012: Biomedical Applications in Molecular, Structural, and Functional Imaging (March 2012)
  • (15) November, L.J., Simon, G.W., Tarbell, T.D., Title, A.M. and Ferguson, S.H. “Large-scale horizontal flows from SOUP observations of solar granulation,” in Theoretical Problems in High Resolution Solar Physics, (NASA/GSFC, 1987) eds. G. Athay and D.S. Spicer, p121.
  • (16) Scherrer, P. H., Bogart, R. S., Bush, R. I., Hoeksema, J. T., Kosovichev, A. G., Schou, J., Rosenberg, W., Springer, L., Tarbell, T. D., Title, A., Wolfson, C. J., Zayer, I. & The MDI Engineering Team, “The Solar Oscillations Investigation - Michelson Doppler Imager,” Solar Phys.,162,129 (1995).
  • (17) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Solar Phys., 275, 207
  • (18) Simon, G. W., Title, A. M., Topka, K. P., Tarbell, T. D., Shine, R. A., Ferguson, S. H., Zirin, H. and the SOUP Team, “On the relation between photospheric flow fields and the magnetic field distribution on the solar surface,” Astrophys. J., 327, 964 (1988)
  • (19) Snodgrass, in The Solar Cycle, ASP Conf. 27, 220, (1992)
  • (20) Sun, D., Roth, S. and Black, M., ”A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them,” Int J Comput Vis 106, 115-137 (2014)
  • (21) Title, A.M., Tarbell, T.D., Acton, L., Duncan, D. and Simon, G.W., “White-light movies of the solar photosphere from the SOUP instrument on Spacelab,” Adv. Space Res. 6,253 (1986).
  • (22) Title, A. M., Tarbell, T. D., Topka, K. P., et al. 1989, Astrophys. J., 336, 475
  • (23) Tsuneta, S., Ichimoto, K., Katsukawa, Y., et al. 2008, Solar Phys., 249, 167
  • (24) Verma, M., Steffen, M., and Denker, C.: 2013, Astron. Astroph. 555, A136.
  • (25) Verma, M. and Denker, C.: 2011, Astron. Astroph. 529, A153.
  • (26) Welsch, B.T., Abbett, W.P., De Rosa, M.L., Fisher, G.H., Georgoulis, M.K., Kusano, K., Longcope, D.W., Ravindra, B., and Schuck, P.W.: 2007, Astrophys. J. 670, 1434.
  • (27) Wereley, S., Gui, L. and Meinhart, C., ”Flow Measurement Techniques for the Microfrontier” in AIAA 2001-0243 at the 39th Aerospace Sciences Meeting (2001)