Next Article in Journal
What Do Leaders Know?
Next Article in Special Issue
3D Reconstruction of Coronal Loops by the Principal Component Analysis
Previous Article in Journal
Time Evolution of Relative Entropies for Anomalous Diffusion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimization of Curvilinear Tracing Applied to Solar Physics and Biophysics

by
Markus J. Aschwanden
1,*,
Bart De Pontieu
1 and
Eugene A. Katrukha
2
1
Lockheed Martin Solar and Astrophysics Laboratory, Building 252, Organization A021S, 3251 Hanover Street, Palo Alto, CA 94304, USA
2
Cell Biology, Faculty of Science, Utrecht University, Padualaan 8, Utrecht, The Netherlands
*
Author to whom correspondence should be addressed.
Entropy 2013, 15(8), 3007-3030; https://doi.org/10.3390/e15083007
Submission received: 23 May 2013 / Revised: 4 July 2013 / Accepted: 18 July 2013 / Published: 26 July 2013
(This article belongs to the Special Issue Advanced Signal Processing in Heliospheric Physics)

Abstract

:
We developed an automated pattern recognition code that is particularly well suited to extract one-dimensional curvilinear features from two-dimensional digital images. A former version of this Oriented Coronal Curved Loop Tracing (OCCULT) code was applied to spacecraft images of magnetic loops in the solar corona, recorded with the NASA spacecraft, Transition Region And Coronal Explorer (TRACE), in extreme ultra-violet wavelengths. Here, we apply an advanced version of this code (OCCULT-2), also, to similar images from the Solar Dynamics Observatory (SDO), to chromospheric H-α images obtained with the Swedish Solar Telescope (SST) and to microscopy images of microtubule filaments in live cells in biophysics. We provide a full analytical description of the code, optimize the control parameters and compare the automated tracing with visual/manual methods. The traced structures differ by up to 16 orders of magnitude in size, which demonstrates the universality of the tracing algorithm.

1. Introduction

Image segmentation is an image processing method that subdivides an image into its constituent regions or objects, which can have the one-dimensional geometry of curvilinear (1D) segments or the two-dimensional (2D) geometry of (fractal) areas. Common techniques include point, line and edge detection, edge linking and boundary detection, Hough transform, thresholding, region-based segmentation, morphological watersheds, etc. (e.g., [1]). Since there exists no omni-potent automated pattern recognition code that works for all types of images equally well, we have to customize suitable algorithms for each data type individually by taking advantage of the particular geometry of the features of interest, using a priori information from the data. In this study, we optimize an automated pattern recognition code to extract magnetized loops from images of the solar corona with the aim of optimum completeness and fidelity. We will demonstrate that the same code works also equally well for microscopic images in biophysics. The particular geometric property of the extracted features is the relatively large curvature radius of coronal magnetic field lines, which generally do not have sharp kinks and corners, but exhibit continuity in the variation of the local curvature radius along their length. Using a related strategy of curvature constraints, coronal loops were extracted also with the directional 2D Morlet wavelet transform ([2]). In addition, solar coronal loops, as well as biological microtubule filaments, have a relatively small cross-section compared with their length, so that they can be treated as curvilinear 1D objects in a tracing method. Tracing of 1D structures with large curvature radii simplifies an automated algorithm enormously, compared with segmentation of 2D regions with arbitrary geometry and, possibly, fractal fine structure [3].
The content of this paper includes a brief description of the automated tracing code (Section 2), applications to images in solar physics (Section 3), to images in biophysics (Section 4), discussion and conclusions (Section 5) and a full analytical description of the code in Appendix A.

2. Description of the Automated Tracing Code

Early versions of Oriented-Connectivity Methods (OCM) applied to solar images were pioneered by Lee, Newman and Gary [4,5]. This code was applied to a Transition Region And Coronal Explorer (TRACE) image, and a total of 57 coronal loops were detected in a solar active region, which supposedly outline the dipolar magnetic field. The results of this code were compared among five numerical codes based on similar curvilinear loop-segmentation methods in a study of Aschwanden [6], in terms of the cumulative size distribution of loop lengths, the median and maximum detected loop lengths, the completeness and detection efficiency, accuracy and flux sensitivity. One of these codes was developed further, which we call Oriented Coronal Curved Loop Tracing (OCCULT) code, and was found to approach the quality of visually and manually traced loops, detecting a total of 272 loop structures in the same TRACE image [7]. In the study, here, we developed this code further and tested it in a larger parameter space and with different types of images. Technical details of the original OCCULT code are given in [7], a concise description of the advanced OCCULT-2 code is given in the following, while a comprehensive analytical description of OCCULT-2 is provided in Appendix A.
(1)
Background suppression: The median, z m e d , of an intensity image, z i j = I 0 ( x i , y i ) , is computed, and the low intensity values with z i j < z m i n are set to the base value, z m i n = z m e d × q m e d , with q m e d being a selectable control parameter, with a default value, q m e d = 1 . 0 , if applied, and q m e d = 0 , if ignored, while a range of 1≲ qmed ≲ 2.5 was found to be useful for noisy data. The median value, zmed, is a good estimate of the background (if the features of interest cover less than 50% of the image area) and can manually be adjusted with the multiplier, qmed, otherwise. The parts of the original image that have intensities below this base level are then rendered with a constant value, are noise-free and will automatically suppress any structure detection in the background below this base level. This new method is more flexible and efficient in suppressing faint background structures.
(2)
Highpass and lowpass filtering: A lowpass filter with a boxcar smoothing constant, n s m 1 , smooths out the data noise (e.g., photon noise with Poisson statistics in astrophysical and microscopy images; e.g., see [8]), while a highpass filter with a boxcar smoothing constant, n s m 2 , enhances the fine structure. The two combined filters represent a bandpass filter (with n s m 1 < n s m 2 ), defined by:
I f i l t e r ( x i , y j ) = s m o o t h [ I 0 ( x i , y j ) , n s m 1 ] - s m o o t h [ I 0 ( x i , y j ) , n s m 2 ]
For theoretical and experimental reasons, a filter combination of n s m 2 = n s m 1 + 2 yields the sharpest enhancement of the intermediate spatial scale and, thus, warrants optimum performance in tracing of curvilinear structures.
(3)
Initialization of loop structures: The code initializes the first structure to be traced from the position, ( x 0 a , y 0 a ) , with the maximum brightness or flux intensity, f 0 a = I 0 ( x 0 a , y 0 a ) , in the original image, I 0 ( x , y ) . Once the full loop has been traced, the area of the detected loop is erased to zero, and the next loop structure is initialized at the position, ( x 0 b , y 0 b ) , with the next flux maximum, f 0 b , in the residual image. The initialization of subsequent loops, ( f 0 c , f 0 d , . . . ) , is continued iteratively, until the residual image becomes entirely zeroed out, the increase of detected structures stagnates or a maximum loop number, N m a x , is reached.
(4)
Loop structure tracing: An initialized structure starting at its flux maximum position, ( x 0 , y 0 ) , is then traced in the forward direction to the first end point of the loop and, then, in the opposite direction from the original starting point to the second endpoint. The two bi-directional segments are then combined into a single uni-directional 1D path, s i = s ( x i , y i ) , i = 1 , . . . , n s . The step-wise tracing along a loop structure position, s i , is carried out by determining, first, the direction of the local ridge (defined by the azimuthal angle, α l , with respect to the x-axis) and, secondly, by determining the local curvature radius, r m . The curved segment that follows a local ridge closest is used as a second-order polynomial to extrapolate the traced loop segment by one incremental step (of Δ s = 1 pixel). This second-order guiding criterion represents an improvement over the first-order guiding criterion used in the previous OCCULT code [7]. The second-order guiding criterion is defined by the brightness distribution, f ( s ) = f [ x s e g ( s ) , y s e g ( s ) ] , along a loop segment with constant curvature radius and length, n s . If the segment follows an ideal ridge with a constant curvature radius and a constant brightness, the summed (or averaged) flux along the ridge segment has a maximum value, while it exhibits a minimum in the perpendicular direction to the ridge, where the brightness profile collapses to a δ-function.
(5)
Loop subtraction in residual image: Once a full loop structure has been traced, the loop area, I 0 ( x i ± w , y j ± w ) , i = 1 , . . . , n s , is set to zero within a half width of w = ( n s m 2 / 2 - 1 ) , so that the area of a former detected loop is not used in the detection of subsequent loops. However, crossing loops can still be connected over a gap.
The original automated loop tracing code (OCCULT) employed the following control parameters: the highpass filter boxcar, n s m 2 , the noise threshold level, N σ , the minimum curvature radius, r m i n , the moving segment length, n s , the directional angle range, Δ α , and a filling factor, q f i l l , for the guiding criterion. In the advanced code (OCCULT-2), we use only two free parameters: the lowpass filter, n s m 1 , and the minimum curvature radius, r m i n . In addition, noise treatment is handled by selecting a typical noise area in the image, as well as the control factor, q m e d , that ignores faint structures below the base level, z m i n = z m e d × q m e d . Thus, the new version of the code offers a simpler choice of control parameters for automated detection of structures.

3. Application to Solar Physics

In the following three subsections, we process three different types of solar images, one from the TRACE spacecraft (Section 3.1), one from the SDO spacecraft (Section 3.2) and one from a ground-based solar telescope (Section 3.3). Some parameters of the analyzed images are listed in Table 1.
Table 1. Parameters of five analyzed images, including the range of brightness in the image ( z m i n , z m a x ), the minimum curvature radius, r m i n (in units of pixels), the bandpass filter ( n s m 1 , n s m 2 ), the total number of detected loops, N d e t ( L > 30 ) , the number of detected long loops, N d e t ( L > 70 ) , and the power-law slope of the cumulative length distribution, p L .
Table 1. Parameters of five analyzed images, including the range of brightness in the image ( z m i n , z m a x ), the minimum curvature radius, r m i n (in units of pixels), the bandpass filter ( n s m 1 , n s m 2 ), the total number of detected loops, N d e t ( L > 30 ) , the number of detected long loops, N d e t ( L > 70 ) , and the power-law slope of the cumulative length distribution, p L .
ImageBrightness rangeMinimum curvatureFilter rangeNumber of all loopsNumber of long loopsPower-law slope
z min , z max r min n sm 1 , n sm 2 N det ( L > 30 ) N det ( L > 70 ) p L
TRACE56, 2606305, 74371343.1
SDO/AIA28, 255309, 115031212.7
SST339, 916303, 517573763.9
Cell-HC416, 65535153, 5208513.1
Cell-LC289, 18699307, 9151392.7

3.1. TRACE Data

The first image we are processing has been used as a standard in a number of previous publications [4,5,6,7], which shows coronal loops in a dipolar active region that was recorded by the TRACE spacecraft on May 5, 1998, 22:21 UT, in the 171 Å wavelength (Figure 1, center). The original image has a size of 1 , 024 × 1 , 024 pixels, with a pixel size of 0.5 ( 360 km on the solar surface). The EUVbrightness or intensity in every pixel is quantified by a data number (DN), which has a minimum of 56 DN and a maximum of 2,606 DN in this image. We show the bandpass-filtered image (with a lowpass boxcar of n s m 1 = 5 pixels and a highpass boxcar of n s m 2 = 7 pixels) in Figure 2. The image shows at least four different textures (Figure 1, a–d): coronal loops (L: curvilinear features), so-called moss regions (M: high-contrast reticulated or spongy features in the center of the image, which represent the footpoints of hot coronal loops), transition region emission (T: low-contrast irregular features in the background) and faint emission areas. The faint image areas show, in addition, a ripple with diagonal stripes at a pedestal level of 57 ± 1 DN (Figure 2, bottom left panel, R) that results from some interference in the electronic readout, which can produce unwanted non-solar structures in the automated detection of curvilinear features.
Figure 1. A solar EUV image of an active region, recorded with the Transition Region And Coronal Explorer (TRACE) spacecraft on May 15, 1998, is shown with a colorscale that has the highest brightness in the white regions (center). In addition, we show ( 100 × 100 pixel) enlargements of four subregions with different textures, which contain (a) coronal loops, (c) electronic ripples, (b) chromospheric and transition region emissions and (d) moss regions with footpoints of hot coronal loops.
Figure 1. A solar EUV image of an active region, recorded with the Transition Region And Coronal Explorer (TRACE) spacecraft on May 15, 1998, is shown with a colorscale that has the highest brightness in the white regions (center). In addition, we show ( 100 × 100 pixel) enlargements of four subregions with different textures, which contain (a) coronal loops, (c) electronic ripples, (b) chromospheric and transition region emissions and (d) moss regions with footpoints of hot coronal loops.
Entropy 15 03007 g001
The automated detection of curvilinear features with the OCCULT-2 code yields a total of 437 loop structures with lengths of L 30 pixels (Figure 2), whereof the longer loops (with lengths of L ≳ 50) coincide well with the 210 visually/manually traced loops (see Figure A1 in [7]). The good agreement between the automated and visually detected loops can also be seen from the cumulative size distributions of loop lengths obtained with both methods (Figure 1, bottom right panel). The two methods detect N = 154 and N = 134 loops with a length of L ≥ 70 pixels, and both distributions have a power-law slope of αL ≈ 2.9 ± 0.1. The challenges of coronal loop detection in this image are confusion and interference from all three types of non-loop structures (Figure 1): chains of “dotted moss structures”, filamentary and specular transition region emission and the diagonal stripes of electronic ripple in the background, which become all comparable with the signal of loop structures once the image contrast is enhanced with a highpass filter.
Figure 2. A bandpass-filtered ( n s m 1 = 5 , n s m 2 = 7 ) version of the original image rendered in Figure 1 is shown (grey scale), with automated loop tracings overlaid (red curves). Cumulative size distributions, N ( > L ) , of loop lengths are also shown (bottom right panel), comparing the automated tracing (red distribution) with visually/manually traced loops (black distribution). The maximum lengths, L m (in pixels), are listed for the longest loops detected with each method.
Figure 2. A bandpass-filtered ( n s m 1 = 5 , n s m 2 = 7 ) version of the original image rendered in Figure 1 is shown (grey scale), with automated loop tracings overlaid (red curves). Cumulative size distributions, N ( > L ) , of loop lengths are also shown (bottom right panel), comparing the automated tracing (red distribution) with visually/manually traced loops (black distribution). The maximum lengths, L m (in pixels), are listed for the longest loops detected with each method.
Entropy 15 03007 g002
Figure 3. The optimization of the highpass filter constant, n s m 2 (y-axis), is shown for the number of detected loops (with lengths longer than 70 pixels) for all analyzed cases, N d e t ( L 70 ) (a–e). The optimization of detected loops as a function of the curvature radius, r m i n , is also shown (f–j).
Figure 3. The optimization of the highpass filter constant, n s m 2 (y-axis), is shown for the number of detected loops (with lengths longer than 70 pixels) for all analyzed cases, N d e t ( L 70 ) (a–e). The optimization of detected loops as a function of the curvature radius, r m i n , is also shown (f–j).
Entropy 15 03007 g003
The successful or false detection in an image, thus, depends most strongly on the chosen low and highpass filter constant, n s m 2 , and the background base level, z m i n = z m e d × q m e d , and, to a lesser degree, on the other control parameters. A good indicator of the completeness and efficiency of automated loop detection is the number of coherently detected long loops, say with a length above 70 pixels, here, N d e t ( L > 70 ) (which is marked with a dotted line in the size distribution in Figure 2, bottom right panel). In Figure 3a, we show how this detection efficiency, N d e t ( n s m 2 ) , varies as a function of the chosen control parameter, n s m 2 . The detected number has a maximum of N d e t = 134 at n s m i 1 = 5 and n s m 2 = 7 , which corresponds to a bandpass filter in the range of 5–7 pixels ( 2 . 5 - 3 . 5 or 1,800–2,500 km on the solar surface). It appears that this is the most typical cross-section width (FWHM) of coronal loops. Other statistical studies of coronal loops observed with TRACE yield similar values (e.g., FWHM = 1 , 420 ± 340 km [9]). Thus, if an image contains structures with a preferential cross-section width, the relevant cross-section range can be bracketed with a bandpass filter ( n s m 1 , n s m 2 ), providing useful a priori information for automated detection of curvilinear features. We conducted tests with all possible filter widths, n s m 1 = 1 , 3 , . . . , 21 and n s m 2 > n s m 1 , and found that the largest number of detected structures virtually always occurs at n s m 2 = n s m 1 + 2 , which can be explained also by the theoretical argument that the best signal-to-noise ratio is obtained for maximum smoothing of the highpass-filtered (unsharp mask) image. We vary also the minimum curvature radius, r m i n , and find a maximum of detected structures at r m i n 30 pixels (Figure 3f).

3.2. SDO/AIA Data

The next image to which we apply our automated loop tracing code is from the Atmospheric Imager Assembly (AIA) on board the Solar Dynamics Observatory (SDO), which replaced the TRACE mission and has been operating since 2010 [10]. AIA has a similar spatial resolution (pixel size 0.6 ) as TRACE (pixel size 0.5 ), but covers the full Sun disk, with an image size of 4 , 096 × 4 , 096 pixels. Figure 4 shows a subimage with a size of 1 , 450 × 650 pixels, which contains a complex of two magnetically coupled active regions, observed on August 03, 2011, 01 UT, in the 171 Å wavelength. This image is currently the subject of nonlinear force-free magnetic modeling (Mark DeRosa, private communication, 2012), and, thus, requires automated loop tracing to constrain the coronal part of the magnetic field configuration [11]. Differences to the TRACE image are the higher sensitivity of the AIA telescopes, different exposure times, the availability of simultaneous images in eight other wavelengths, different image compression and no apparent electronic ripple in the CCDreadout, which all affect the automated detection of faint structures. Synthesized loop tracings from multiple wavelength filters have been proven to provide a more robust and representative subset of loop structures for magnetic modeling than loop tracings from a single filter image ([12]).
We vary the lowpass filter constant in the range of n s m 1 = 1 , . . . , 21 , set the highpass filter constant to n s m 2 = n s m 1 + 2 and find a maximum detection rate of N d e t ( L > 70 ) = 121 loop structures for n s m 1 = 9 and n s m 2 = 11 (Figure 3b). We vary also the minimum curvature radius in the range of r m i n = 10 , . . . , 100 pixels and find a maximum detection rate of N d e t ( L > 70 ) = 121 at r m i n = 30 pixels (Figure 3g).
Figure 4. Bandpass-filtered image of an active region complex observed with Atmospheric Imager Assembly (AIA)/Solar Dynamics Observatory (SDO) on 3 Aug. 2011, 01 UT, 171 Å shown as (a) an intensity image, as (b) a bandpass-filtered version with n s m 1 = 9 and n s m 2 = 11 and (c) overlaid with automatically traced loop structures , where the low-intensity values below the median of f = 75 DN (data number) are blocked out (grey areas).
Figure 4. Bandpass-filtered image of an active region complex observed with Atmospheric Imager Assembly (AIA)/Solar Dynamics Observatory (SDO) on 3 Aug. 2011, 01 UT, 171 Å shown as (a) an intensity image, as (b) a bandpass-filtered version with n s m 1 = 9 and n s m 2 = 11 and (c) overlaid with automatically traced loop structures , where the low-intensity values below the median of f = 75 DN (data number) are blocked out (grey areas).
Entropy 15 03007 g004
Figure 5. High-resolution image of the solar Active Region 10380, recorded on June 16, 2003, with the Swedish 1-m Solar Telescope (SST) on (a) La Palma Spain and (b) automated tracing of curvilinear structures with a lowpass filter of n s m 1 = 3 pixels, a highpass filter of n s m 2 = 5 pixels and a minimum curvature radius of r m i n = 30 pixels, tracing out 1,757 curvilinear segments.
Figure 5. High-resolution image of the solar Active Region 10380, recorded on June 16, 2003, with the Swedish 1-m Solar Telescope (SST) on (a) La Palma Spain and (b) automated tracing of curvilinear structures with a lowpass filter of n s m 1 = 3 pixels, a highpass filter of n s m 2 = 5 pixels and a minimum curvature radius of r m i n = 30 pixels, tracing out 1,757 curvilinear segments.
Entropy 15 03007 g005

3.3. SST Data

Now, we apply automated loop tracing to a solar image in a completely different wavelength, namely in the Hα 6,563 Å line, the first line of the Balmer series of hydrogen. Figure 5 shows such an image of the solar upper chromosphere, which displays chromospheric spicules in the right side of the picture and filamentary structures in the upper chromosphere and transition region (in altitudes ≈ 2,000–5,000 km above the solar surface) [14,15]. The image was taken with the Swedish 1-m Solar Telescope (SST) on La Palma, Spain, using a tunable filter, tuned to the blue-shifted line wing of the Hα 6,563 Å line. The spicules are jets of moving gas at a lower temperature than the million degree hot corona and flow upward from the chromosphere to the transition region with a speed of 15 km s - 1 [14]. The image has a size of 1 , 507 × 999 pixels ( 62 × 41 Mm on the solar surface), with a pixel size of 0 . 041 (or 30 km on the solar surface). This picture is particularly intriguing for automated tracing of curvilinear structures, because of the ubiquity and complexity of fine structure.
This SST image has such a high contrast so that there is no significant noise that affects curvilinear tracing. Thus, we set the background level to zero ( q m e d = 0 . 0 ) in the OCCULT-2 code. We vary the lowpass filter constant in the range of n s m 1 = 1 , . . . , 21 , set the highpass filter constant to n s m 2 = n s m 1 + 2 and vary the minimum curvature radius in the range of r m i n = 10 , . . . , 100 pixels. We find a maximum number of detected structures of N d e t ( L > 70 ) = 376 (with a length above 70 pixels) at n s m 1 = 3 and N s m 2 = 5 (Figure 3c) and r m i n = 40 pixels (Figure 3h). Extending to shorter segments with lengths of L > 30 pixels, the automated tracing code identifies a total of 1,757 curvilinear segments, which outline the patterns of the flow field in the upper chromosphere (Figure 5 b).

4. Applications to Biophysics

Finally, we apply our automated tracing code to images obtained in cellular biophysics, in order to test the versatility and universality of the OCCULT-2 code. As an example, we chose microscopy images of microtubules, but the same approach could be applied to any other filament-like structures: intermediate or actin filaments, fibrin, etc. Microtubules are long stiff polymers that are part of cytoskeleton (internal cellular scaffolding). They are important for cell division, motility and organization ([16,17]). The entangled network of microtubules serves as a “railroad” system for delivery of cargoes by molecular motors and, also, as a stiff carcass, controlling cellular mechanics ([17,18]). It is a dynamic network adapting and changing in time in response to external cues. The automated extraction of the microtubules network’s configuration from microscopy images and movies can provide insights on the mechanisms of these changes.
Figure 6 shows two images of cells from two different cell lines with microtubules networks of different density. Cells belonging to the CHO cell line (Figure 6a, image Cell-HC) are small and have a sparse microtubule network. Cells from the U2OS line (Figure 6b, image Cell-LC) are larger and contain a more dense radial network. Since microtubules are transparent to the visible light, a fluorescent tag is used to observe them in the living cells. In the analyzed images (Figure 6), microtubules were labeled with a green fluorescent protein (GFP) ([19]) that adsorbs light at 490 nm and has an emission peak at the 510 nm wavelength. The images were acquired with a spinning disk confocal microscope using the corresponding GFP’s emission filter. The magnified image was projected on the 16-bit chip of camera with dimensions of 512 × 512 pixels, resulting in the final image pixel size of 66 nm. The width of microtubules filaments measured in the electron microscopy studies is about 25 nm ([20]). Due to the diffraction limit, the effective width of microtubules in the image is much larger and is defined by the microscopy setup and the wavelength used ([21]). In our case, it is approximately equal to a half of the emission wavelength (510 nm) that is close to the measured microtubule width of 4 . 7 ± 2 . 0 pixels. The two pictures shown in Figure 6a,b show cases of opposite contrast, one with high contrast (Figure 6a, image Cell-HC) and one with low contrast (Figure 6b, image Cell-LC). The difference in contrast is explained by a different amount of GFP fluorescent tag associated with microtubules. The automated tracing of the cell filaments provides their locations inside cell and curvature radii, from which flexural rigidity and mechanical stress can be inferred.
Figure 6. False-colored images of microtubule filaments in live cells: (a) Cell-HC) a high-contrast image of the CHO cell and (b Cell-LC)a low-contrast image of the U2OS cell . The size of the images corresponds to 34 μm. The automated curvilinear tracing of both images was carried out with the parameters: n s m 1 = 3 , n s m 2 = 5 and r m i n = 15 for Cell-HC and n s m 1 = 7 , n s m 2 = 9 and r m i n = 30 for Cell-LC.
Figure 6. False-colored images of microtubule filaments in live cells: (a) Cell-HC) a high-contrast image of the CHO cell and (b Cell-LC)a low-contrast image of the U2OS cell . The size of the images corresponds to 34 μm. The automated curvilinear tracing of both images was carried out with the parameters: n s m 1 = 3 , n s m 2 = 5 and r m i n = 15 for Cell-HC and n s m 1 = 7 , n s m 2 = 9 and r m i n = 30 for Cell-LC.
Entropy 15 03007 g006
We optimize the automated filament tracing by varying the bypass filter constants and find a maximum number of detected filaments (with lengths L > 70 pixels) with n s m 1 = 3 , n s m 2 = 5 , r m i n = 50 for the high-contrast cell image and n s m 1 = 7 , n s m 2 = 9 , r m i n = 40 for the low-contrast cell image (Figure 3). The low-contrast image has a higher degree of noise, and thus, the filter constants have to be adjusted to a larger value for the low-contrast cell image ( n s m 1 = 7 ) than for the high-contrast image ( n s m 1 = 3 ).

5. Discussion and Conclusions

5.1. Optimization of Automated curvilinear Tracing

The efficiency and accuracy of automated curvilinear tracing can be controlled by a number of tuning or control parameters. For the present OCCULT-2 code, we have three independent control parameters ( n s m 1 , r m i n , q m e d ). A prerequisite for the application of curvilinear tracing codes is the assumption that the structures of interest have a much smaller width, w, than their length, l, so that they can be represented by a 1D path, which is generally curved, possibly limited by a minimum curvature radius, r m i n . The larger the minimum curvature radius is, the less ambiguity there is for tracing of crossing 1D structures. In this study, we explored the parameter space of the control parameters, n s m 1 and r m i n , in order to optimize the performance of the automated tracing code.
The lowpass filter ( n s m 1 ) and highpass filter ( n s m 2 ) represent the brackets or scale range of a bandpass filter, bracketing the range of cross-sections of detected structures. We expect to find structures with the smallest width preferentially in images with a high signal-to-noise ratio, while noisy images require more smoothing to enhance a fine structure and, thus, tend to have larger widths, due to the smearing effect of the smoothing. We found the narrowest structures indeed in the two images with the highest contrast (i.e., the SST and Cell-HC image), where highpass filters with boxcars of n s m 2 = 5 pixels were used. In images with lower contrast, optimum performance occurred for highpass filters of n s m 1 = 7 , . . . , 11 pixels. In addition, we found that the smallest bandpass filters produce the sharpest structures and, thus, the highest detection rate of structures. Since a symmetric boxcar requires odd numbers, n s m = 1 , 3 , 5 , . . . , the smallest difference between a lowpass and a highpass filter is two, and thus, the optimum combination is expected to be n s m 2 = n s m 1 + 2 , which we indeed confirmed also experimentally.
For the minimum curvature radius, we found optimum performance for a typical range of r m i n 30 , . . , 50 pixels (Figure 3), which seems not to depend on the contrast of the image.
The control parameter, q m e d = 1 , suppresses faint structures in an image that are below the median value of the image brightness. A meaningful value for this parameter depends very much what fraction of the image contains bright structures of interest, which has to be decided, depending on the area ratio of structures of interest to non-relevant background area. If the background has comparable brightness to the structures of interest, a thresholded separation may be impossible, but may still succeed if the background has a different texture than the curvilinear features (see examples in Figure 1, where moss and transition region emission have a different texture than coronal loops, while background ripples have the same texture as straight loops and can only be separated by a threshold control parameter). Future efforts aim at synthesizing the loop tracings from multi-wavelength image sets, which are more robust and representative for magnetic modeling than single-wavelength images ([12]).
In conclusion, we recommend the following procedure to achieve optimum performance of the curvilinear tracing code (OCCULT-2): (1) Start with the following recommended control settings: n s m 1 = 1 , n s m 2 = 3 , r m i n = 30 , and q m e d = 1 . 0 ; (2) Vary the filter combination in the range of n s m 1 = 1 , 3 , . . . , 15 , while setting n s m 2 = n s m 1 + 2 , to find the maximum detection rate for a given loop length (e.g., here, we used L > 70 pixels); (3) Low-contrast images are likely to require higher highpass filter values, n s m 2 , than high-contrast images; (4) Vary the minimum curvature radius, r m i n , within some range to find the maximum detection rate of structures; (5) If the code yields a lot of random structures in obvious background areas, increase the base level factor, q m e d > 1 .
The software of the numerical code OCCULT-2 is publicly accessible in Interactive Data Language (IDL) in the Solar Software (SSW) package. A tutorial and example is accessible at the authors homepage ([13]).

5.2. Solar Applications

The automated tracing of curvilinear structures in solar physics has mostly been applied to extreme-ultraviolet or soft X-ray images, which show magnetized coronal loops that outline the coronal magnetic field, due to the low plasma-β parameter in the solar corona (which is defined as the ratio of the thermal to the magnetic pressure). Thus, coronal loops represent the perfect tracers of the otherwise invisible coronal magnetic field. Standard magnetic field models of the solar corona, or parts of it, such as sunspot regions and active regions, have been modeled by extrapolating a photospheric magnetogram, either with a potential field solution or a force-free solution of Maxwell’s equations. However, since the chromosphere was found not to satisfy the force-free condition [22], extrapolations of photospheric magnetograms do not exactly render the coronal magnetic field, while coronal loops outline the true coronal magnetic field. It is therefore desirable to trace such coronal loops in EUV and soft X-ray images and to use them to constrain a magnetic field solution. Such attempts have been performed with single EUV images, as well as with stereoscopic EUV image pairs [11]. Forward-fitting of theoretical magnetic field models to traced coronal loops is able to discriminate between potential field and force-free field models, as well as to quantify the free magnetic energy (that is released in solar flares) and Lorentz forces, e.g., [23]. curvilinear tracing of filamentary structures in the chromosphere and transition region (Figure 5) may also help to constrain the horizontal magnetic field components in the non-force-free regions, which has been used in pre-processing of solar force-free magnetic field extrapolations [24].
One fundamental limitation of automated coronal loop tracing is the confusion by background structures resulting from EUV emission from the transition region, which has generally a cooler temperature than the coronal loops. Future efforts may use multi-wavelength image data sets to discriminate EUV emission from the chromosphere, transition region and the corona by its temperature, using a a deconvolution of the multi-wavelength temperature filter response functions in terms of a differential emission measure (DEM) method.

5.3. Biological Applications

The method of curvilinear tracing is increasingly used in the analysis of biological and medical images, such as to characterize blood vessel tracking in retinal images [25,26,27], neurons [28], dendritic spines [29] or microtubule tracing in fluorescent and phase-contrast microscopy [30,31,32], as shown in Figure 6. In our experiment with a high-contrast (Figure 6 a,c) and low-contrast image (Figure 6b,d), we demonstrated that a bypass or bandpass filter with a bandpass factor of n s m 2 = n s m 1 + 2 enhances the structures to an optimum contrast for automated tracing. Moreover, we found that the bandpass filter of low-contrast images requires a larger width ( n s m 2 = 9 pixels in Figure 6d) than in a high-contrast image ( n s m 2 = 5 pixels in Figure 6c).
In our optimization exercise, we concentrated mostly on the completeness of detected long curvilinear segments, but future efforts may also consider the optimization of linking multiple segments that are interrupted with gaps or subject to crossings. The efficiency and reliability of automated curvilinear algorithms became more important with the massive increase of imaging data over the last decade, which exceeds our limited capabilities of visual inspection. Note that the algorithm used here tracks curvilinear structures that differ by 16 orders of magnitude in size, from 66 nm to 360 km (pixel size in images).

Acknowledgments

We thank Ilya Grigoriev and Anna Akhmanova for providing the images of labeled microtubules. We acknowledge also helpful discussions with Karel Schrijver, Mark DeRosa, Allen Gary, Jake Lee, Bernd Inhester, James McAteer, Peter Gallagher, Alex Young, Alex Engels, Patrick Shami, Narges Fathalian, Fatemeh Amirkhanlou, Hossein Safari and participants of the 3rd Solar Image Processing Workshop in Dublin, Ireland, September 6–8, 2006, the 4th Solar Image Processing Workshop in Baltimore, Maryland, October 26–30, 2008, the 5th Solar Imaging Processing Workshop in Les Diablerets, Switzerland, September 13–16, 2010, and the 6th Solar Imaging Processing Workshop at Montana State University, Bozeman, August 13–16, 2012. Part of the work was supported by the NASA TRACE contract (NAS5-38099) and the NASA SDO/AIA contract, NNG04EA00C.

Conflict of Interest

The authors declare no conflict of interest.

References

  1. Gonzales, R.C.; Woods, R.E. Digital Image Processing; Pearson Prentice Hall: Upper Saddle River, NJ, USA, 2008; pp. 689–794. [Google Scholar]
  2. Biskri, S.; Antoine, J.-P.; Inhester, B.; Mekideche, F. Extraction of solar coronal magnetic loops with the directional 2D morlet wavelet transform. Solar Phys. 2010, 262, 373–385. [Google Scholar] [CrossRef]
  3. McAteer, R.T.J.; Gallagher, P.T.; Ireland, J. Statistics of active region complexity: A large-scale fractal dimension survey. Astrophys. J. 2005, 631, 628–635. [Google Scholar] [CrossRef]
  4. Lee, J.K.; Newman, T.S.; Gary, G.A. Automated Detection of Solar Loops by the Oriented Connectivity Method. In Proceedings of the 17th International Conference on Pattern Recognition (ICPR), Cambridge, UK, 23–26 August 2004; p. 315.
  5. Lee, J.K.; Newman, T.S.; Gary, G.A. Oriented connectivity-based method for segmenting solar loops. Pattern Recognit. 2006, 39, 246–259. [Google Scholar] [CrossRef]
  6. Aschwanden, M.J.; Lee, J.K.; Gary, G.A.; Smith, M.; Inhester, B. Comparison of five numerical codes for automated tracing of coronal loops. Solar Phys. 2008, 248, 359–377. [Google Scholar] [CrossRef]
  7. Aschwanden, M.J. A code for automated tracing of coronal loops approaching visual perception. Solar Phys. 2010, 262, 399–423. [Google Scholar] [CrossRef]
  8. Pawley, J.B.; Masters, B.R. Handbook of biological confocal microscopy. J. Biomed. Opt. 2008, 13, 029902. [Google Scholar] [CrossRef]
  9. Aschwanden, M.J.; Nightingale, R.W. Elementary loop structures in the solar corona analyzed from TRACE triple-filter images. Astrophys. J. 2005, 633, 499–517. [Google Scholar] [CrossRef]
  10. Lemen, J.R.; Title, A.M.; Akin, D.J.; Boerner, P.F.; Chou, C.; Drake, J.F.; Duncan, D.W.; Edwards, C.G.; Friedlaender, F.M.; Heyman, G.F.; et al. The Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO). Solar Phys. 2012, 275, 17–40. [Google Scholar] [CrossRef]
  11. Aschwanden, M.J. Nonlinear force-free magnetic field fitting to coronal loops with and without stereoscopy. Astrophys. J. 2013. [Google Scholar] [CrossRef]
  12. Aschwanden, M.J.; Sun, Y. The Active Region 11158 during the 2011 February 15 X-Class Flare: Nonlinear Force-free Magnetic Fields and Free Energies Inferred from Automated Loop Tracing. 2013. in preparation. [Google Scholar]
  13. Aschwanden, M.J. Available online: http://lmsal.com/∼aschwand/software/ (accessed on 23 May 2013).
  14. De Pontieu, B.; Erdelyi, R.; James, S.P. Solar chromospheric spicules from the leakage of photospheric oscillations and flows. Nature 2004, 430, 536–539. [Google Scholar] [CrossRef] [PubMed]
  15. De Pontieu, B.; Carlsson, M.; Rouppe van der Voort, L.; Löfdahl, M.; van Noort, M.; Nordlund, Å.; Scharmer, G. Rapid temporal variability of faculae: High-resolution observations and modeling. Astrophys. J. 2006, 1405–1420. [Google Scholar] [CrossRef]
  16. Kaverina, I.; Straube, A. Regulation of cell migration by dynamic microtubules. Semin. Cell Dev. Biol. 2011, 22, 968–974. [Google Scholar] [CrossRef] [PubMed]
  17. De Forges, H.; Bouissou, A.; Perez, F. Interplay between microtubule dynamics and intracellular organization. Int. J. Biochem. Cell Biol. 2012, 44, 266–274. [Google Scholar] [CrossRef] [PubMed]
  18. Subramaniam, R.; Kapoor, T.M. Building complexity: Insights into self-organized assembly of microtubule-based architectures. Dev. Cell 2012, 22, 874–885. [Google Scholar] [CrossRef] [PubMed]
  19. Tsien, R.Y. The green fluorescent protein. Annu. Rev. Biochem. 1998, 67, 506–544. [Google Scholar] [CrossRef] [PubMed]
  20. Wade, R.H.; Chrétien, D. Cryoelectron microscopy of microtubules. J. Struct. Biol. 1993, 110, 1–27. [Google Scholar] [CrossRef] [PubMed]
  21. Abbe, E. Note on the proper definition of the amplifying power of a lens or a lens-system. J. R. Microsc. Soc. 1884, 4, 348–351. [Google Scholar] [CrossRef]
  22. Metcalf, T.R.; Jiao, L.; Uitenbroek, H.; McClymont, A.N.; Canfield, R.C. Is the solar chromospheric magnetic field force-free? Astrophys. J. 1995, 439, 474–481. [Google Scholar] [CrossRef]
  23. Gary, G.A. Evaluation of a selected case of the minimum dissipative rate method for non-force-free solar magnetic field extrapolation. Solar Phys. 2009, 257, 271–286. [Google Scholar] [CrossRef]
  24. Fuhrmann, M.; Seehafer, N.; Valori, G.; Wiegelmann, T. A comparison of preprocessing methods for solar force-free magnetic field extrapolation. Astron. Astrophys. 2011, 526, A70. [Google Scholar] [CrossRef]
  25. Jiang, Y.; Bainbridge-Smith, A.; Morris, A.B. Blood Vessel Tracking in Retinal Images. In Proceedings of the Image and Vision Computing, New Zealand, 5–7 December 2007; pp. 126–131.
  26. Martinez-Perez, M.E.; Hughes, A.D.; Stanton, A.V.; Thom, S.A.; Barath, A.A.; Parker, K.H. Segmentation of Retinal Blood Bessels Based on the Second Directional Derivative and Region Growing. In Proceedings of the International Conference on Image Processing (ICIP 99), Kobe, Japan, 24–28 October 1999; Volume 2, pp. 173–176.
  27. Tang, J.; Acton, S.T. Vessel Boundary Tracking for Intravital Microscopy via Multiscale Gradient Vector Flow Snakes. In Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing ICASSP 04, Singapore, 17–21 May 2004; Volume 51, pp. 316–324.
  28. Xiong, G.; Zhou, X.; Degterev, A.; Ji, L.; Wong, S.T. Automated neurite labeling and analysis in fluorescence microscopy images. Cytom. Part A 2006, 69A, 494–505. [Google Scholar] [CrossRef] [PubMed]
  29. Zhang, Y.; Zhou, X.; Witt, R.M.; Sabatini, B.L.; Adjeroh, D.; Wong, S.T. Dendritic spine detection using curvilinear structure detector and Ida classifier. NeuroImage 2007, 36, 346–360. [Google Scholar] [CrossRef] [PubMed]
  30. Brangwynne, C.P.; Koenderingk, G.H.; Barry, E.; Dogic, Z.; MacKintosh, F.C.; Weitz, D.A. Bending dynamics of fluctuating biopolymers probed by automated high-resolution filament tracking. Biphys. J. 2007, 93, 346–359. [Google Scholar] [CrossRef] [PubMed]
  31. Koulgi, P.; Sargin, M.E.; Rose, K.; Manjunath, B.S. Graphical Model-based Tracing of curvilinear Structures in Bio-image Sequences. In Proceedings of the IEEE International Conference on Pattern Recognition, Istanbul, Turkey, 23–26 August 2010; pp. 2596–2599.
  32. Sargin, M.E.; Altinok, A.; Rose, K.; Manjunath, B.S. Tracing curvilinear Structures in Live Cell Images. In Proceedings of the IEEE International Conference on Image Processing ICIP 07, San Antonio, TX, USA, 16–19 September 2007; Volume 6, pp. 285–288.

Appendix

A.1. Analytical Description of the OCCULT-2 Code

The Oriented Coronal Curved Loop Tracing (OCCULT-2) code version 2 is an improved version of the original OCCULT code described in Aschwanden (2010). The improvements include: (1) a “curved” guiding segment that is adjusted to the local curvature radius of a traced loop (representing the second-order term of a polynomial), rather than the linear (first-order polynomial) guiding segment used in the original version; (2) suppression of faint structures; (3) bypass filtering instead of highpass filtering; and (4) simplification of selectable free parameters.
The input is a simple two-dimensional (2D) image, z i j , with pixel numbers, i = 0 , . . . , n x - 1 , on the x-axis and j = 0 , . . . , n y - 1 on the y-axis, respectively. The output is a number of curvilinear structures (also called “loops”, for short), which are parameterized in terms of x- and y-coordinates, [ x ( s k ) , y ( s k ) ] , where the loop length coordinate, s k = 0 , . . . , n s , is given in steps of Δ s = 1 pixel, so that for all k = 0 , . . . , n s - 1 :
Δ s k = [ x ( s k - x ( s k - 1 ) ] 2 + [ y ( s k ) - y ( s k - 1 ) ] 2 = 1 , k = 0 , . . . , n s - 1
with n s the number of loop points for each loop.
The goal of the algorithm, OCCULT-2, is to retrieve as many curvilinear structures as possible, without picking up false signals of non-existing structures in the noise of the image, which has some probability to form chains of random points in a curved array configuration. The challenges are, therefore, to evaluate an optimum threshold level that separates existing loops from noise structures and to retrieve the real curvilinear structures as completely as possible, without subdividing them into partial loop segments. Our strategy to obtain a fast numeric code is to retrieve the loops in a one-dimensional search algorithm, because any two-dimensional concept has a computation time that grows with the square of the image size. The one-dimensional parameter space is essentially the loop length coordinate, s k , k = 0 , . . . , n s - 1 . In addition, we define two other independent parameters in each loop point, which can be considered as the first-order and second-order term of a polynomial, namely the local direction angle, α l , l = 0 , . . . , n α , and the curvature radius, r m , m = 0 , . . . , n r . The algorithm selects iteratively the brightest position ( x i , y j ) in the image and starts a bi-directional loop tracing, determining first the local direction, α l ( x i , y j ) , and curvature radius, r m ( x i , y j ) , at the starting point and continues tracing the loop within a small (guided) range of the local curvature radius, which is the principle of “orientation-guided tracing”. Therefore, we are dealing essentially with 1D-tracing in a 5D-parameter space ( x i , y j , s k , α l , r m ) , which we parameterize with the five indices ( i , j , k , l , m ) that have the index ranges, i = 0 , . . . , n x - 1 , j = 0 , . . . , n y - 1 , k = 0 , . . . , n s , l = 0 , . . . , n α , m = 0 , . . . , n r . Specifically, we define the arrays:
s k b i = Δ s ( k - n s 2 ) , k = 0 , . . . , n s - 1
for a symmetric bi-directional array, used in the search of the local direction, α l , and:
s k u n i = Δ s k , k = 0 , . . . , n s - 1
for a uni-directional array, used in the search for the curvature radius in the forward-direction of a traced structure. For the directional angle, α l , which is only determined at the starting point of the loop, we use a fixed array with a resolution of one degree (or π / 180 radian),
α l = π l n α , l = 0 , . . . , n α , n α = 180
For the curvature radii, r m , we use a reciprocal scaling in order to obtain a uniform distribution of directional angles at the end of a curvature segment:
r m = r m i n [ - 1 + 2 m / ( n r - 1 ) ] , n r = 30
This parameterization covers positive and negative curvature radii in the ranges of [ - , - r m i n ] and [ r m i n , + ] . The choice of an even number, n r , prevents the singularity, r m = ± .
Now, we describe the consecutive steps of the algorithm one by one, which follow more or less the same flow chart as depicted in Figure 1 of Aschwanden (2010).
(1) Image base level ( q b a s e ): If the image consists of high-contrast structures without unwanted secondary structures in the background, we do not have to worry about the image base level and can set it to the lowest value (which should be zero or positive in astrophysical images that record an intensity or brightness). However, if there are unwanted structures at a faint brightness level, we can just set the image base level z b a s e above the brightness level of unwanted structures, which we parameterize with the factor, q m e d , in units of the median brightness level, z m e d = m e d i a n [ z i , j ] :
z b a s e = z m e d × q b a s e
so that the corrected brightness, z i , j , in each pixel fulfills the condition:
z i , j z b a s e
For q m e d = 0 , the image is unchanged, while the image, z i , j , appears to be flat in the fainter half area of the image, if set q m e d = 1 . 0 . Therefore, the value, q m e d , can be adjusted depending on the estimated fraction of the image that is covered with structures of interest. For solar images, this feature offers a convenient way to filter out coronal loops in active regions (which are bright) from unwanted structures in the quiet Sun (which are faint).
(2) Bandpass filtering ( n s m 1 , n s m 2 ): The tracing of curvilinear structures is considerably eased by enhancing of fine structures within a chosen bandpass that corresponds to the typical width, n w , of structures of interest, which is typically a few pixels, and assuming that the curvilinear structures has a much longer length, n s , than width, i.e., n s n w . We accomplish the enhancement with a bandpass filter, which consists of a lowpass filter with a boxcar, n s m 1 , and a highpass filter with a boxcar, n s m 2 :
z i , j f i l t e r = s m o o t h [ z i , j , n s m 1 ] - s m o o t h [ z i , j , n s m 2 ]
which filters out broad structures with widths, n w nsm2 (highpass filter), but smooths out fine structure with a boxcar of nsm1 < nsm2. We experimented with a large number of combinations (nsm1, nsm2) and found that the optimum choice is:
n s m 2 = n s m 1 + 2
which follows the principle of maximum possible smoothing of fine structures with a given width, nsm2. The values for the smoothing with a symmetric boxcar has to be an odd integer, i.e., nsm1 = 1, 3, 5, ...), which implies nsm2 = 3, 5, 7, ..., where the lowest value, nsm1 = 1, corresponds to the original image without any smoothing. This experimentally tested relationship for the optimum choice of bandpass filters (nsm1, nsm2) reduces also the possible parameter space by one dimension, and thus, we have to search only for nsm1 = 1, 3, 5, ..., while using nsm2 = nsm1 + 2.
(3) Noise threshold: Our algorithm starts with the brightest curvilinear structure and proceeds to fainter structures, and thus, we have to find a stop criterion that halts the procedure when it reaches the level of data noise. Such a noise threshold has to be determined empirically for every image, since there are many sources of possible data noise. Testing many images with completely different data types, we found that a most reasonable threshold level can be determined by an interactive choice of a noise area in the image that contains typical data noise, but little structures of interest. Such an image area can be characterized by the pixel ranges ( i n 1 : i n 2 , j n 1 : j n 2 ) . In this noise area, we determine the median brightness, z m e d n o i s e , and define a noise threshold at the doubled value:
z t h r e s h = 2 × ( z m e d n o i s e > 0 )
using only the pixels with positive values in the noise area (in order to prevent a too low value of z t h r e s h = 0 in the case when more than 50% of the noisy pixels are below the previously chosen base value, z b a s e ). The rationale for a factor of two in the threshold level comes from the fact that the median separates out only half of the noisy pixels, while the double value would separate out all noisy pixels if the distribution of noisy pixel values follows a linear relationship.
(4) Start of curvilinear structures: We are ready now to trace the first curvilinear structure. We determine the location ( i 0 , j 0 ) of the absolute brightness maximum, z 0 , in the bandpass-filtered image, z i , j f i l t e r :
z 0 = z ( x 0 , y 0 ) = m a x [ z i , j f i l t e r ( x i , y j ) ]
The rationale for the choice of this starting point is the expectation to trace first the most significant structure in the bandpass-filtered image, which can then be continued by going to the next-significant structure, once the tracing of the first structure has been successfully completed and the corresponding loop area is eliminated in a residual image. Consequently, the maximum of the brightness in the residual image marks the second-brightest structure, and we can repeat the same procedure by tracing the next loop.
(5) Loop direction at starting point: The next element of the structure to be traced is the first-order term of a polynomial, the direction angle, α l , which is also the direction of a possible ridge that outlines the local segment of the structure. We determine this directional angle simply by measuring the flux averaged over a straight loop segment symmetrically placed over the starting point ( i 0 , j 0 ) and rotated over a full range of possible angles, α l , from 0 to 180 degrees. The x,y-coordinates of this linear segment are, with the array, s k b i , defined by Equation (5):
x k , l = i 0 + s k b i cos α l
y k , l = j 0 + s k b i sin α l
where the index, k, runs along the length, s k , of the segment and the index, l, denotes a particular angle, α l . Among the set of angular values, α l , we determine the maximum of the summed flux in each rotated segment:
z m a x ( α l ) = m a x l 1 n s k = 0 n s - 1 z f i l t e r ( x k , l , y k , l )
which yields the local direction, α m a x = α l ( l = l m a x ) . In the ideal case of a straight ridge with a constant value, z 0 , along the ridge segment with length n s pixels and zero-values outside, a value of z = z 0 is found, while the value for a segment in perpendicular direction to the ridge is much smaller, namely z = 1 / n s . This ridge criterion works even for a close succession of parallel ridges, in which case it is still a factor of two smaller than in the parallel direction, i.e., z = 1 / 2 .
(6) Local curvature radius: Next, we determine the second-order term of a polynomial that follows the loop to be traced. We define a directional angle, β 0 , that is in perpendicular direction to the loop and defines the direction where all centers of possible curvature radii are located (that intersect the structure at location, ( x 0 , y 0 ) ), see geometric definition of the angles, α and β, in Figure A1:
β 0 = α 0 + π 2
The location ( x c , y c ) of the curvature center with a minimum curvature radius, r m i n , is then found at:
x c = x 0 + r m i n cos β 0
y c = y c + r m i n sin β 0
The loci of all curvature centers ( x m , y m ) of a set of curvature radii, r m = r m i n / [ - 1 + 2 m / ( n r - 1 ) ] (Equation 6), is then found at:
x m = x 0 + ( x c - x 0 ) r m r m i n
y m = y 0 + ( y c - y 0 ) r m r m i n
Figure A1. Geometry of curvature radii centers ( x r m , y r m ) located on a line at angle β (dash-dotted line), perpendicular to the tangent at angle α (solid line), that intersects a curvilinear feature (thick solid curve) at position, ( x 0 , y 0 ) . The angle, γ, indicates the half angular range of the curved guiding segment (thick solid line).
Figure A1. Geometry of curvature radii centers ( x r m , y r m ) located on a line at angle β (dash-dotted line), perpendicular to the tangent at angle α (solid line), that intersects a curvilinear feature (thick solid curve) at position, ( x 0 , y 0 ) . The angle, γ, indicates the half angular range of the curved guiding segment (thick solid line).
Entropy 15 03007 g007
Since we want to follow a loop along a curved segment for every possible curvature radius, r m , we determine the coordinates for each segment point, s k . It is useful to define the angle, β m , of the line that connects a curvature center ( x m , y m ) with a curved segment point, s k :
β m = β 0 + σ d i r s k r m
where σ d i r = ± 1 has two opposite signs, depending on the forward or backward tracing of a loop. The x,y-coordinates ( x k m , y k m ) of the loop segment, s k , are then:
x k m = x m - r m cos ( β m )
y k m = y m - r m sin ( β m )
In order to determine the curvature radius, r m , that fits the local loop segment best, we search for the segment with the maximum flux along the curve with radius, r m :
z m a x = m a x m [ 1 n r k = 0 n s - 1 z i , j ( x k m , y k m ) ]
Since we now know the optimum curvature radius, r m , we can trace the loop incrementally by a step, Δ s , and extrapolate the position ( x k + 1 , y k + 1 ) and angle, α k + 1 , by:
α k + 1 = α k + σ d i r Δ s r m
α m i d = ( α k + α k + 1 ) 2
x k + 1 = x k + Δ s cos α m i d + ( π / 2 ) ( 1 + σ d i r )
y k + 1 = y k + Δ s sin α m i d + ( π / 2 ) ( 1 + σ d i r )
(7) Bidirectional tracing: Step (6) describes the extrapolation or tracing of the loop segment from position s k = ( x k , y k ) to s k + 1 ( x k + 1 , y k + 1 ) by an incremental length step, Δ s . This step is repeated, starting from some arbitrary starting point, s 0 , inside the segment, until the first endpoint, s n s 1 , where the traced curvilinear structure seems to end. The first endpoint generally is demarcated at a location where the bandpass-filtered image has zero or negative values, so one could just detect the first image pixel along the guiding segment, s k , that is non-positive. However, in order to allow for some minor gaps in noisy structures, it turned out to be more reliable to define a stop criterion when at least a few pixels have a nonpositive value, say n g a p = 3 pixels. An example of a traced loop is shown in Figure A2, which illustrates that a loop can reliably be traced even in the presence of secondary loops that intersect at small angles.
After completing the first half segment ( σ d i r = + 1 ), we repeat the same procedure from the midpoint ( x 0 , y 0 ) in the opposite direction ( σ d i r = - 1 ), until we stop at the second endpoint at s n s 2 . We combine, then, the two segments, [ s 0 , s n 1 ] and [ s 0 , s n 2 ] , by reversing one segment in order to maintain the same direction and concatenating the two equal-directed segments into a single loop structure with indices, [ s 0 , . . . , s n s ] , where the length is s n = ( s n 1 + s n 2 ) .
(8) Loop iteration: Steps (4) to (7) describe the tracing of a single loop, specified by the n L coordinates, ( x i , y i ) = [ x ( s k ) , y ( s k ) ] , k = 0 , . . . , n L - 1 . In order to proceed to the next loop, we want first to erase the area of this traced loop in the residual image, so that no previously traced loop segment is traced a second time in a consecutive loop. The loop is erased in the residual image simply by setting those image pixels to zero that have a distance of d n w = ( n s m 2 / 2 - 1 ) from the loop coordinates, i.e., z f i l t e r [ x i ± n w , y j ± n w ] = 0 . An example of an image area with dense coverage of loops is shown in Figure A3, where some loops are traced at flux levels close to the noise threshold.
Figure A2. Example of loop tracing in the pixel area [525:680, 453:607] of the image shown in Figure 2. Loop #19 is traced (blue crosses) over a length of 115 pixels (orange numbers), crossing another structure at a small angle. The curves at position 115 indicate the three curved segments that have been used in the tracing of the last loop point. The black contours indicate the bandpass-filtered difference image and the red contours indicate the previously traced and erased structures in the residual difference image.
Figure A2. Example of loop tracing in the pixel area [525:680, 453:607] of the image shown in Figure 2. Loop #19 is traced (blue crosses) over a length of 115 pixels (orange numbers), crossing another structure at a small angle. The curves at position 115 indicate the three curved segments that have been used in the tracing of the last loop point. The black contours indicate the bandpass-filtered difference image and the red contours indicate the previously traced and erased structures in the residual difference image.
Entropy 15 03007 g008
Figure A3. Example of loop tracings in area [300:550, 600:800] of the full image shown in Figure 2. The grey scale indicates the bandpass-filtered image, ( n s m 1 = 5 , n s m 2 = 7 ) , and the loop tracings are shown with red curves.
Figure A3. Example of loop tracings in area [300:550, 600:800] of the full image shown in Figure 2. The grey scale indicates the bandpass-filtered image, ( n s m 1 = 5 , n s m 2 = 7 ) , and the loop tracings are shown with red curves.
Entropy 15 03007 g009
(9) Loop parameters: After we described analytically the numerical code, we summarize the free parameters and the dependent parameters (that do not require any selection by the user of the code). The code has the following free or input parameters: the lowpass filter boxcar, n s m 1 , the minimum curvature radius, r m i n , the base level factor, q m e d , in units of image median brightness values, and the corner coordinates of the noise area, [ i n 1 : i n 2 , j n 1 : j n 2 ] . All other parameters are dependent, or a fixed constant, such as the tracing step, Δ s = 1 , the highpass filter boxcar, n s m 2 = n s m 1 + 2 , the half width of the erased loop cross-sections, n w = n n s m 2 / 2 - 1 , and the loop termination gap, n g a p = 3 .

Share and Cite

MDPI and ACS Style

Aschwanden, M.J.; De Pontieu, B.; Katrukha, E.A. Optimization of Curvilinear Tracing Applied to Solar Physics and Biophysics. Entropy 2013, 15, 3007-3030. https://doi.org/10.3390/e15083007

AMA Style

Aschwanden MJ, De Pontieu B, Katrukha EA. Optimization of Curvilinear Tracing Applied to Solar Physics and Biophysics. Entropy. 2013; 15(8):3007-3030. https://doi.org/10.3390/e15083007

Chicago/Turabian Style

Aschwanden, Markus J., Bart De Pontieu, and Eugene A. Katrukha. 2013. "Optimization of Curvilinear Tracing Applied to Solar Physics and Biophysics" Entropy 15, no. 8: 3007-3030. https://doi.org/10.3390/e15083007

APA Style

Aschwanden, M. J., De Pontieu, B., & Katrukha, E. A. (2013). Optimization of Curvilinear Tracing Applied to Solar Physics and Biophysics. Entropy, 15(8), 3007-3030. https://doi.org/10.3390/e15083007

Article Metrics

Back to TopTop