Optimization of Curvi-Linear Tracing Applied to Solar Physics and Biophysics

We developed an automated pattern recognition code that is particularly well suited to extract one-dimensional curvi-linear features from two-dimensional digital images. A former version of this {\sl Oriented Coronal CUrved Loop Tracing (OCCULT)} code was applied to spacecraft images of magnetic loops in the solar corona, recorded with the NASA spacecraft {\sl Transition Region And Coronal Explorer (TRACE)} in extreme ultra-violet wavelengths. Here we apply an advanced version of this code ({\sl OCCULT-2}) also to similar images from the {\sl Solar Dynamics Observatory (SDO)}, to chromospheric H-$\alpha$ images obtained with the {\sl 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.


15
Image segmentation is an image processing method that subdivides an image into its constituent 16 regions or objects, which can have the one-dimensional geometry of curvi-linear (1D) segments, or original image that have intensities below this base level, are then rendered with a constant value, 58 are noise-free, and will automatically suppress any structure detection in the background below 59 this base level. This new method is more flexible and efficient in suppressing faint background 60 structures. 61 2. Highpass and lowpass filtering: A lowpass filter with a boxcar smoothing constant n sm1 smoothes out the data noise (e.g., photon noise with Poisson statistics in astrophysical and microscopy images, e.g., see [24]), while a highpass filter with a boxcar smoothing constant n sm2 enhances the fine structure. The two combined filters represent a bandpass filter (with n sm1 < n sm2 ), defined by For theoretical and experimental reasons, a filter combination of n sm2 = n sm1 + 2 yields the  I 0 (x i ± w, y j ± w), i = 1, ..., n s , is set to zero within a half width of w = (n sm2 /2 − 1), so 88 that the area of a former detected loop is not used in the detection of subsequent loops. However, 89 crossing loops can still be connected over a gap. Table 1. Parameters of five analyzed images, including the range of brightness in the image (z min , z max ), the minimum curvature radius r min (in units of pixels), the bandpass filter (n sm1 , n sm2 ), the total number of detected loops N det (L > 30), the number of detected long loops N det (L > 70), and the powerlaw slope of the cumulative length distribution p L .      virtually always occurs at n sm2 = n sm1 + 2, which can be explained also by the theoretical argument that 143 the best signal-to-noise ratio is obtained for maximum smoothing of the highpass-filtered (unsharp mask) 144 image. We vary also the minimum curvature radius r min and find a maximum of detected structures at       rigidity and mechanical stress can be inferred. 215 We optimize the automated filament tracing by varying the bipass filter constants and find a maximum 216 number of detected filaments (with lengths L > 70 pixels) with n sm1 = 3, n sm2 = 5, r min = 50 for the 217 high-contrast cell image, and n sm1 = 7, n sm2 = 9, r min = 40 for the low-contrast cell image (Fig. 3).

218
The low-contrast image has a higher degree of noise, and thus the filter constants have to be adjusted to 219 a larger value for the low-contrast cell image (n sm1 = 7) than for the high-contrast image (n sm1 = 3).  The lowpass filter (n sm1 ) and highpass filter (n sm2 ) represent the brackets or scale range of a bandpass 231 filter, bracketing the range of cross-sections of detected structures. We expect to find structures with the 232 smallest width preferentially in images with a high signal-to-noise ratio, while noisy images require 233 more smoothing to enhance fine structure, and thus tend to have larger widths due to the smearing effect 234 of the smoothing. We found the narrowest structures indeed in the two images with the highest contrast 235 (i.e., the SST and Cell-HC image), where highpass filters with boxcars of n sm2 = 5 pixels were used.

236
In images with lower contrast, optimum performance occurred for highpass filters of n sm1 = 7, ..., 11 237 pixels. In addition, we found that the smallest bandpass filters produce the sharpest structures and thus 238 the highest detection rate of structures. Since a symmetric boxcar requires odd numbers n sm = 1, 3, 5, ..., 239 the smallest difference between a lowpass and a highpass filter is 2, and thus the optimum combination 240 is expected to be n sm2 = n sm1 + 2, which we indeed confirmed also experimentlly.

241
For the minimum curvature radius we found optimum performance for a typical range of r min ≈ 242 30, .., 50 pixels (Fig. 3), which seems not to depend on the contrast of the image. representative for magnetic modeling than single-wavelength images ([6]).

253
In conclusion, we recommend the following procedure to achieve optimum performance of the curvi- field models to traced coronal loops is able to discriminate between potential field and force-free field 278 models, as well as to quantify the free magnetic energy (that is released in solar flares) and Lorentz 279 forces, e.g., [13].  Fig. 6. In our experiment with a high-contrast (Fig. 6, top panel) and low-contrast image (Fig. 6, bottom   293 panel), we demonstrated that a bipass or bandpass filter with a bandpass factor of n sm2 = n sm1 + 2 294 enhances the structures to an optimum contrast for automated tracing. Moreover we found that the 295 bandpass filter of low-contrast images requires a larger width (n sm2 = 9 pixels in Fig. 6, bottom panel) 296 than in a high-contrast image (n sm2 = 5 pixels in Fig. 6, top panel).

297
In our optimization exercise we concentrated mostly on the completeness of detected long curvi-linear

311
The input is a simple 2-dimensional (2D) image z ij , 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 curvi-linear 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, with n s the number of loop points for each loop.

312
The goal of the algorithm OCCULT-2 is to retrieve as many curvi-linear 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 curvi-linear 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". So 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 5 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, for a symmetric bi-directional array, used in the search of the local direction α l , and 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,

316
(1) Image base level (q base ): 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 base above the brightness level of unwanted structures, which we parameterize with the factor q med in units of the median brightness level z med = median[z i,j ], so that the corrected brightness z ′ i,j in each pixel fulfills the condition For q med = 0 the image is unchanged, while the image z ′ i,j appears to be flat in the fainter half area of the 317 image, if set q med = 1.0. So, the value q med can be adjusted depending on the estimated fraction of the 318 image that is covered with structures of interest. For solar images, this feature offers a convenient way 319 to filter out coronal loops in active regions (which are bright) from unwanted structures in the Quiet Sun 320 (which are faint).

321
(2) Bandpass Filtering (n sm1 , n sm2 ): The tracing of curvi-linear 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 curvi-linear 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 sm1 and a highpass filter with a boxcar n sm2 , which filters out broad structures with widths n w > ∼ n sm2 (highpass filter), but smoothes out fine structure with a boxcar of n sm1 < n sm2 . We experimented with a large number of combinations (n sm1 , n sm2 ) and found that the optimum choice is, which follows the principle of maximum possible smoothing of fine structures with a given width n sm2 .

327
(3) Noise Threshold: Our algorithm starts with the brightest curvi-linear 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 n1 : i n2 , j n1 : j n2 ). In this noise area we determine the median brightness z noise med and define a noise threshold at the doubled value, using only the pixels with positive values in the noise area (in order to prevent a too low value of z thresh = 328 0 in the case when more than 50% of the noisy pixels are below the previously chosen base value z base ).

329
The rationale for the factor 2 in the threshold level comes from the fact that the median separates out

332
(4) Start of Curvi-Linear Structures: We are ready now to trace the first curvi-linear structure. We determine the location (i 0 , j 0 ) of the absolute brightness maximum z 0 in the bandpass-filtered image The rationale for the choice of this starting point is the expectation to trace first the most significant 333 structure in the bandpass-filtered image, which can then be continued by going to the next-significant 334 structure, once the tracing of the first structure has been successfully completed and the corresponding 335 loop area is eliminated in a residual image. Consequently, the maximum of the brightness in the residual 336 image marks the second-brightest structure and we can repeat the same procedure by tracing the next 337 loop.

338
(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 bi k defined by Eq. (5), 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, which yields the local direction α max = α l (l = l max ). In the ideal case of a straight ridge with a constant 339 value z 0 along the ridge segment with length n s pixels and zero-values outside, a value of z = z 0 is 340 found, while the value for a segment in perpendicular direction to the ridge is much smaller, namely 341 z ⊥ = 1/n s . This ridge criterion works even for a close succession of parallel ridges, in which case it is 342 still a factor of 2 smaller than in parallel direction, i.e., z ⊥ = 1/2.

343
(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 Fig. 7, The location (x c , y c ) of the curvature center with a minimum curvature radius r min is then found at, y c = y c + r min sin β 0 .
The loci of all curvature centers (x m , y m ) of a set of curvature radii r m = r min /[−1 + 2m/(n r − 1)] (Eq. 6) is then found at,  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 , where σ dir = ±1 has two opposite signs, depending on the forward or backward tracing of a loop. The x,y-coordinates (x km , y km ) of the loop segment s k is then, 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 , Since we know now 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 x k+1 = x k + ∆s cos [α mid + (π/2)(1 + σ dir )] , y k+1 = y k + ∆s sin [α mid + (π/2)(1 + σ dir )] .
(28) (7) Bidirectional tracing: The step (6) describes the extrapolation or tracing of the loop segment from 344 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, 345 starting from some arbitrary starting point s 0 inside the segment until the first endpoint s n s1 , where the 346 traced curvi-linear structure seems to end. The first endpoint generally is demarcated at a location where 347 the bandpass-filtered image has zero or negative values, so one could just detect the first image pixel 348 along the guiding segment s k that is non-positive. However, in order to allow for some minor gaps in 349 noisy structures, it turned out to be more reliable to define a stop criterion when a least a few pixels 350 have a nonpositive value, say n gap = 3 pixels. An example of a traced loop is shown in Fig. 8, which   351 illustrates that a loop can reliably be traced even in the presence of secondary loops that intersect at small 352 angles.

353
After completing the first half segment (σ dir = +1), we repeat the same procedure from the midpoint 354 (x 0 , y 0 ) in opposite direction (σ dir = −1), until we stop at the second endpoint at s n s2 . We combine than where the length is s n = (s n 1 + s n2 ).  i.e., z f ilter [x i ± n w , y j ± n w ] = 0. An example of an image area with dense coverage of loops is shown 364 in Fig. 9, where some loops are traced at flux levels close to the noise threshold.

365
(9) Loop Parameters: After we described analytically the numerical code, we summarize the free 366 parameters and the dependent parameters (that do not require any selection by the user of the code).

367
The code has the following free or input parameters: the lowpass filter boxcar n sm1 , the minimum