FilFinder Tutorial

This tutorial demonstrates the FilFinder algorithm on a simulated data set. The updated algorithm from FilFinder2D is used here, which is valid for versions >1.5. This tutorial was tested with python 3.6.

The example data is included in the github repository here.


In [120]:
%matplotlib inline
import matplotlib.pyplot as plt
import astropy.units as u

# Optional settings for the plots. Comment out if needed.
import seaborn as sb
sb.set_context('poster')

import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (12., 9.6)

Input Data

There are two caveats to the input data:

1) All angular and physical conversions assume that pixels in the image can be treated as squares. FilFinder2D is not aware of any axis misalignments! If you're data does not have aligned celestial axes, we recommend reprojecting the data onto a square grid.

2) The beam size is characterized by the major axis, assuming a 2D Gaussian beam shape. If the beam size of your data is highly elliptical, it is recommended to convolve the data to a circular beam.

FilFinder2D accepts several input types, including a FITS HDU and numpy arrays.


In [121]:
from fil_finder import FilFinder2D
from astropy.io import fits

hdu = fits.open("../examples/filaments_updatedhdr.fits")[0]

fil = FilFinder2D(hdu)

In [122]:
# HDU data as an array
arr = hdu.data
hdr = hdu.header

fil = FilFinder2D(arr)


/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filfinder2D.py:138: UserWarning: No beam width given. Using 0 pixels.
  warnings.warn("No beam width given. Using 0 pixels.")

In this case, no WCS information is given and all results will be returned in pixel units. Angular units can be returned when the header is specified:


In [123]:
fil = FilFinder2D(arr, header=hdr)

If spectral-cube is installed, the Projection or Slice classes can also be passed to FilFinder2D:


In [124]:
from spectral_cube import Projection

proj = Projection.from_hdu(hdu)

fil = FilFinder2D(proj)


WARNING: A 'NAXIS1' keyword already exists in this header.  Inserting duplicate keyword. [astropy.io.fits.header]

Other Inputs to FilFinder2D:

Note that numerical inputs must be given as ~astropy.units.Quantity object with the appropriate unit.

Distance -- To facilitate conversions to physical units, a distance can be given to FilFinder2D:


In [125]:
fil = FilFinder2D(hdu, distance=250 * u.pc)

Angular Scale -- If no header information is given, the pixel-to-angular conversion can be given:


In [126]:
fil = FilFinder2D(arr, ang_scale=0.1 * u.deg)


/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/base_conversions.py:55: UserWarning: Cannot find 'BMAJ' in the header. Try installing the `radio_beam` package for loading header information.
  warn("Cannot find 'BMAJ' in the header. Try installing"
/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/base_conversions.py:63: UserWarning: Cannot find 'BMIN' in the header. Assuming circular beam.
  warn("Cannot find 'BMIN' in the header. Assuming circular beam.")
/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/base_conversions.py:69: UserWarning: Cannot find 'BPA' in the header. Assuming PA of 0.
  warn("Cannot find 'BPA' in the header. Assuming PA of 0.")
/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filfinder2D.py:138: UserWarning: No beam width given. Using 0 pixels.
  warnings.warn("No beam width given. Using 0 pixels.")

Beamwidth -- If the major axis of the beam is contained in the header, it will be automatically read in. If that information is not in the header, the beam size can be passed separately:


In [127]:
fil = FilFinder2D(hdu, beamwidth=10 * u.arcsec)

Custom Filament Masks -- If you have a pre-computed filament mask, the mask array can be passed:


In [128]:
# Example custom mask
mask = hdu.data > 1.

fil = FilFinder2D(hdu, mask=mask)

The custom mask must have the same shape as the inputed image.

Save Name -- A prefix for saved plots and table can be given:


In [129]:
fil = FilFinder2D(hdu, save_name="FilFinder_Output")

For the purposes of this tutorial, we will assume that the data has WCS information and a well-defined distance:


In [130]:
fil = FilFinder2D(hdu, distance=260 * u.pc)

The beamwidth is $24''$ and is defined in the header.

Image Preprocessing

Prior to creating the mask, it can be helpful to first flatten the image of bright compact sources. FilFinder2D uses an arctan transform, where the data are first normalized by some percentile value of the data:


In [131]:
fil.preprocess_image(flatten_percent=95)

In [132]:
plt.subplot(121)
plt.imshow(fil.image.value, origin='lower')
plt.title("Image")
plt.subplot(122)
plt.imshow(fil.flat_img.value, origin='lower')
plt.title("Flattened Image")
plt.tight_layout()


If a percentile is not given, FilFinder2D.preprocess_image will try to fit a log-normal distribution to the data and will set the threshold at $\mu + 2 \sigma$. There are no checks for the quality of the fit. Use only if you are confident that the brightness distribution is close to a log-normal.

If you wish to run the masking procedure without flattening the image, use the command:


In [133]:
fil.preprocess_image(skip_flatten=True)

The original image will be set as fil.flat_img and used in the masking step.

For this example, we will use the prior flattened image.


In [134]:
fil.preprocess_image(flatten_percent=95)

Masking

Creating the filament mask is a complex process performed by FilFinder2D.create_mask. There are several parameters that set the masking behaviour.

If a FITS header and distance were provided at the beginning, FilFinder2D will use default guesses based on typical filament sizes from Herschel studies of the Gould Belt clouds (e.g., Koch & Rosolowsky 2015). These choices will not be the optimal settings in general, and we recommend trying different different parameter setting before using the resulting mask for the analysis.

If a distance was not provided, these parameters must be set. FilFinder2D will raise errors until the required parameters are given.

This simulated data set is an example where the default FilFinder2D settings do not provide an ideal filament mask:


In [135]:
fil.create_mask(verbose=True)


Most of the filamentary structure has been ignored by the mask. There are several parameters that can be set to improve the mask:

  • glob_thresh -- Set a minimum intensity for a pixel to be included in the mask. This is useful for removing noisy regions in the data from the mask. Must have the same units as fil.image.
  • adapt_thresh -- The width of the element used for the adaptive thresholding mask. This is primarily the step that picks out the filamentary structure. The element size should be similar to the width of the expected filamentary structure. The default here, when distance is provided, is 0.1 pc.
  • smooth_size -- It is often helpful to smooth the data before calculating the mask. By smoothing in small scales, small noise variations are removed resulting in a simpler skeleton structure. The default is set to 0.05 pc.
  • size_thresh -- The minimum number of pixels a region of the mask must have to be considered real. The default is set by assuming a minimum filament size to be an ellipse with a 0.1 pc width and length of 0.5 pc. Most data sets will require this parameter to be manually set, as is used below.
  • regrid -- If the pixel width of adapt_thresh is less than 40 pixels, the resulting mask may be fragmented due to pixelization. To increase the size of adapt_thresh, regrid interpolates the data onto a larger grid, calculates the mask on the larger grid, then interpolates the mask at the original image size.
  • border_masking -- Observational maps may not fill the entire image, and the edges of the mapped regions tend to be noisier. border_masking finds regions of NaNs along the edge of the map and tries to remove noisy regions near the edges. Its behaviour can be controlled using border_kwargs, where the size of a NaN region (size), size of the filter used to define noisy edges (filt_width), and the number of times to apply that filter (eros_iter) can be controlled.
  • fill_hole_size -- If there are holes within a skeleton, fill_hole_size can be used to fill in holes smaller than the given size.
  • use_existing_mask -- If you gave a user-defined mask when calling FilFinder2D, enable this parameter to skip making a new mask.

Varying a few of these parameters will produce a much improved mask. First, since the data go right to the edges of the image, we can disable border_masking:


In [136]:
fil.create_mask(verbose=True, border_masking=False)


The mask now extends right to the edge of the data. However, only one structure was retained. This occurs because size_thresh is too large. We can manually set the value in pixel units. The size must have units of area:


In [137]:
fil.create_mask(verbose=True, border_masking=False, size_thresh=400 * u.pix**2)


Much better! Most of the filamentary structure is now being included in the mask.

This simulated image does not have noise added in, however, most data sets will. The cell below demonstrates how to set glob_thresh to avoid noisy regions:


In [138]:
# Define the noise value. As a demonstration, say values below the 20th percentile here are dominated by noise
noise_level = np.percentile(fil.image, 20)
noise_level

plt.imshow(fil.image.value > noise_level, origin='lower')


Out[138]:
<matplotlib.image.AxesImage at 0x7f2c6841e908>

The dark regions will be excluded from the final mask. The filament mask with the threshold is then:


In [139]:
fil.create_mask(verbose=True, border_masking=False, size_thresh=400 * u.pix**2, glob_thresh=0.0267)


A few small region have been removed compared to the previous mask, but the structure is largely unchanged in this case.

This is a usable mask for the filament analysis. The effects of altering the other parameters are shown in Koch & Rosolowsky 2015.

Try varying each parameter to assess its affect on your data.

If you gave a user-defined mask at the beginning, run:


In [140]:
fil.create_mask(use_existing_mask=True)


/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filfinder2D.py:288: UserWarning: Using inputted mask. Skipping creation of anew mask.
  warnings.warn("Using inputted mask. Skipping creation of a"

Skeletonization

The next step is to reduce the mask into single-pixel-width skeleton objects. These skeletons will define the location of a filament and its path.

In FilFinder2D, the medial axis is defined to be the skeleton:


In [141]:
fil.medskel(verbose=True)


Skeletons: Pruning & Length

We are now prepared to analyze the filaments. The first analysis step includes two parts: pruning the skeleton structures and finding the filament lengths. The first part removes small/unimportant spurs on the skeletons. To ensure important portions of the skeleton are retained, however, both parts are performed together.

Each skeleton is converted into a graph object using the networkx package. We use the graph to find the longest path through the skeleton, which is used to define the structure's length. All branches in the skeleton away from this longest path are eligible to be pruned off.

This process is handled by the FilFinder2D.analyze_skeletons function. When using verbose=True, a ton of plots will get returned. To save you some scrolling, the verbose mode is highlighted for just one filament below.

With just the default settings:


In [142]:
fil.analyze_skeletons()

plt.imshow(fil.skeleton, origin='lower')


Out[142]:
<matplotlib.image.AxesImage at 0x7f2c6812c860>

The skeletons are largely the same, with only short branches removed.

The default settings use minimum skeleton and branch lengths based off of the beam size. To be kept, a branch must be at least three times the length of the beam and a skeleton must have a length of 5 times the beam. Practically, this will only remove very small features.

These parameters, and ones related to the pruning, can be manually set:

  • prune_criteria -- The criteria for removing a branch can be altered. The default ('all') uses a mix of the average intensity along the branch and its length. The length alone can be used for pruning with prune_criteria='length'. All branches below this length will be removed. Finally, only the intensity can be used for pruning (prune_criteria='intensity'). A branch is kept in this case by comparing the average intensity along the branch to the average over the whole filament. The critical fraction that determines whether a branch is important is set by relintens_thresh.
  • relintens_thresh -- Set the critical intensity comparison for intensity-based pruning.
  • nbeam_lengths -- Number of beam widths a skeleton must have to be considered a valid structure. Default is 5.
  • branch_nbeam_lengths -- Number of beam widths a branch must have to avoid pruning. Default is 3.
  • skel_thresh -- Minimum length for a skeleton. Overrides nbeam_lengths. Must have a unit of length.
  • branch_thresh -- Minimum length for a branch. Overrides branch_nbeam_lengths. Must have a unit of length.
  • max_prune_iter -- Number of pruning iterations. The default is 10, which works well for multiple data sets used in the testing process. A warning is returned if the maximum is reached. New in FilFinder2D!

Here we will highlight the effect of pruning away moderately long branches. Note that re-running FilFinder2D.analyze_skeletons will start on the output from the previous call, not that original skeleton from FilFinder2D.medskel.


In [143]:
fil.analyze_skeletons(branch_thresh=40 * u.pix, prune_criteria='length')

plt.imshow(fil.skeleton, origin='lower')
plt.contour(fil.skeleton_longpath, colors='r')


Out[143]:
<matplotlib.contour.QuadContourSet at 0x7f2c6938c588>

The structure have now been significantly pruned. The red contours highlight the longest paths through each skeleton.

If we continue to increase the branch threshold, the skeletons will converge to the longest path structures:


In [144]:
fil.analyze_skeletons(branch_thresh=400 * u.pix, prune_criteria='length')

plt.imshow(fil.skeleton, origin='lower')
plt.contour(fil.skeleton_longpath, colors='r')


Out[144]:
<matplotlib.contour.QuadContourSet at 0x7f2c686d03c8>

This is an extreme case of pruning and a significant amount of real structure was removed. We will return to a less pruned version to use for the rest of the tutorial:


In [145]:
fil.medskel(verbose=False)
fil.analyze_skeletons(branch_thresh=5 * u.pix, prune_criteria='length')

plt.imshow(fil.skeleton, origin='lower')
plt.contour(fil.skeleton_longpath, colors='r')


Out[145]:
<matplotlib.contour.QuadContourSet at 0x7f2c63f5a8d0>

Another new feature of FilFinder2D is that each filament has its own analysis class defined in fil.filaments:


In [146]:
fil.filaments


Out[146]:
[<fil_finder.filament.Filament2D at 0x7f2c68088908>,
 <fil_finder.filament.Filament2D at 0x7f2c68123eb8>,
 <fil_finder.filament.Filament2D at 0x7f2c68123358>,
 <fil_finder.filament.Filament2D at 0x7f2c681233c8>,
 <fil_finder.filament.Filament2D at 0x7f2c681232e8>,
 <fil_finder.filament.Filament2D at 0x7f2c70ef64a8>,
 <fil_finder.filament.Filament2D at 0x7f2c68b69358>,
 <fil_finder.filament.Filament2D at 0x7f2c683ede80>,
 <fil_finder.filament.Filament2D at 0x7f2c684698d0>,
 <fil_finder.filament.Filament2D at 0x7f2c68469cf8>,
 <fil_finder.filament.Filament2D at 0x7f2c68469eb8>,
 <fil_finder.filament.Filament2D at 0x7f2c68469390>]

This allows for each skeleton to be analyzed independently, in case your analysis requires fine-tuning.

A separate tutorial on the Filament2D class is available from the docs page. We will highlight some of the features here to show the plotting outputs. Each Filament2D class does not contain that entire image, however, to avoid making multiple copies of the data.

The first filament is quite large with a lot of structure. We can plot the output from FilFinder2D.analyze_skeletons for just one filament with:


In [147]:
fil1 = fil.filaments[0]
fil1.skeleton_analysis(fil.image, verbose=True, branch_thresh=5 * u.pix, prune_criteria='length')


Three plots are returned:

  • The labeled branch array (left) with intersection points removed and the equivalent graph structure (right).
  • The longest path through the skeleton (left) and the same labeled branch array (right) as above.
  • The final, pruned skeleton structure.

Only one set of plots is shown after the iterative pruning has been finished.

The lengths of the filament's longest paths are now calculated:


In [148]:
fil.lengths()


Out[148]:
$[453.75945,~270.73506,~159.46804,~127.63961,~156.30866,~138.61017,~46.112698,~165.75231,~32.313708,~115.56854,~72.083261,~65.698485] \; \mathrm{pix}$

The default output is in pixel units, but if the angular and physical scales are defined, they can be converted into other units:


In [149]:
fil.lengths(u.deg)


Out[149]:
$[0.75626575,~0.45122511,~0.26578006,~0.21273268,~0.26051443,~0.23101696,~0.076854497,~0.27625385,~0.053856181,~0.19261424,~0.12013877,~0.10949747] \; \mathrm{{}^{\circ}}$

In [150]:
fil.lengths(u.pc)


Out[150]:
$[3.4318251,~2.0475946,~1.2060717,~0.9653503,~1.182177,~1.0483217,~0.34875465,~1.2536002,~0.2443916,~0.87405568,~0.54517244,~0.49688378] \; \mathrm{pc}$

The properties of the branches are also saved in the FilFinder2D.branch_properties dictionary. This includes the length of each branch, the average intensity, the skeleton pixels of the branch, and the number of branches in each skeleton:


In [151]:
fil.branch_properties.keys()


Out[151]:
dict_keys(['length', 'intensity', 'pixels', 'number'])

In [152]:
fil.branch_properties['number']


Out[152]:
array([36, 15,  3,  7, 11, 10,  1,  8,  1,  8,  1,  3])

In [153]:
fil.branch_properties['length'][0]


Out[153]:
$[62.284271,~20.485281,~6.2426407,~13.656854,~7,~20.071068,~10.656854,~14.313708,~12.828427,~10.242641,~18.313708,~26.313708,~67.627417,~58.769553,~5.4142136,~15.485281,~15.071068,~12.242641,~25.727922,~16.313708,~7.8284271,~5.6568542,~7.2426407,~20.313708,~19.727922,~7.6568542,~30.970563,~56.112698,~14.899495,~23.142136,~8.2426407,~36.384776,~16.142136,~34.556349,~11.828427,~51.355339] \; \mathrm{pix}$

Note that the pixels are defined with respect to the cut-out structures in Filament2D. These offsets are contained in FilFinder2D.filament_extents. See the FilFinder2D tutorial for more information.

The branch lengths can also be returned with:


In [154]:
fil.branch_lengths(u.pix)[0]


Out[154]:
$[62.284271,~20.485281,~6.2426407,~13.656854,~7,~20.071068,~10.656854,~14.313708,~12.828427,~10.242641,~18.313708,~26.313708,~67.627417,~58.769553,~5.4142136,~15.485281,~15.071068,~12.242641,~25.727922,~16.313708,~7.8284271,~5.6568542,~7.2426407,~20.313708,~19.727922,~7.6568542,~30.970563,~56.112698,~14.899495,~23.142136,~8.2426407,~36.384776,~16.142136,~34.556349,~11.828427,~51.355339] \; \mathrm{pix}$

In [155]:
fil.branch_lengths(u.pc)[0]


Out[155]:
$[0.47106176,~0.1549321,~0.047213675,~0.10328806,~0.052941654,~0.15179936,~0.080598784,~0.10825591,~0.097022593,~0.077466048,~0.13850829,~0.19901304,~0.51147247,~0.44447962,~0.040948203,~0.11711663,~0.11398389,~0.092592235,~0.19458268,~0.1233821,~0.059207126,~0.042783317,~0.054776768,~0.15363448,~0.14920412,~0.057909504,~0.23423326,~0.42438558,~0.11268627,~0.17502613,~0.062339862,~0.27518146,~0.12208448,~0.2613529,~0.089459499,~0.38840523] \; \mathrm{pc}$

Curvature and Orientation

A filament's curvature and orientation are calculated using a modified version of the Rolling Hough Transform (RHT). This can be run either on the longest path skeletons or on individual branches.

The default setting is to run on the longest path skeletons:


In [156]:
fil.exec_rht()
fil1.plot_rht_distrib()


The RHT distribution is shown for the first skeleton, along with its longest path skeleton. The polar plot shows the distribution as a function of $2\theta$. Since there is no preferred direction, $0$ and $\pi$ are equivalent direction for a filament, and so the distribution is defined over $\theta \in [-\pi/2, \pi/2)$. Plotting the distribution as $2\theta$ makes it easier to visualize with discontinuities. The solid green line shows the mean orientation of the filament, and the curvature region is indicated by the angle between the dashed blue lines.

The RHT distribution is built by measuring the angular distribution in a circular region around each pixel in the skeleton, then accumulating the distribution over all pixels in the skeleton. There are three parameters that affect the distribution:

  • radius -- the radius of the circular region to use. The default is 10 pixels. The region must be large enough to avoid pixelization (causing spikes at 0, 45, and 90 deg) but small enough to only include the local filament direction.
  • ntheta -- The number of bins in $\theta$ to calculate the distribution at. Default is 180.
  • background_percentile -- The accumulation process used to create the distribution will create a constant background level over $\theta$. Peaks in the distribution are better characterized by removing this constant level. The default setting is to subtract the 25th percentile from the distribution.

The RHT returns the orientation and curvature of each filament. The orientation is defined as the circular mean and the curvature is the interquartile region about the mean. See the documentation for the definitions.


In [157]:
fil.orientation


Out[157]:
$[-0.09372123,~-0.4688631,~0.046228431,~0.47607807,~-0.56841827,~-0.96319154,~-0.57605962,~-1.0046175,~-0.028319612,~-0.41163082,~-0.74576315,~0.12105612] \; \mathrm{rad}$

In [158]:
fil.curvature


Out[158]:
$[0.58720166,~0.76920664,~0.50536306,~0.89124894,~0.77025275,~0.73970447,~0.44143597,~1.076277,~0.46711777,~0.73452577,~0.54890669,~0.55313247] \; \mathrm{rad}$

It can be more useful to run this analysis on individual branches to understand the distribution of orientation and curvature across the whole map. This can be performed by enabling branches=True:


In [159]:
fil.exec_rht(branches=True, min_branch_length=5 * u.pix)

There is no default plot setting in this case.

An additional parameter is enabled in this mode: min_branch_lengths. This avoids running the RHT on very short branches, where pixelization will lead to large spikes towards the axis directions.

The outputs are contained in FilFinder2D.orientation_branches and FilFinder2D.curvature_branches, which return a list of lists for each filament. These can be visualized as distributions:


In [160]:
_ = plt.hist(fil.orientation_branches[0].value, bins=10)
plt.xlabel("Orientation (rad)")


Out[160]:
Text(0.5,0,'Orientation (rad)')

In [161]:
all_orient = np.array([orient.value for fil_orient in fil.orientation_branches for orient in fil_orient])
# Short, excluded branches have NaNs
all_orient = all_orient[np.isfinite(all_orient)]

_ = plt.hist(all_orient, bins=10)
plt.xlabel("Orientation (rad)")


Out[161]:
Text(0.5,0,'Orientation (rad)')

No orientation is strongly preferred in the example data.

Radial Profiles and Widths

FilFinder2D finds filament widths by creating radial profiles centered on the skeleton. A simple model can then be fit to the radial profile to find the width.

There are several parameters used to control the creation of the radial profile:

  • max_dist -- The maximum radial distance to build the radial profile to. Must be given in units of length (pixel, degree, pc, etc...). The default is 10 pixels. In order to not bias the fit, the profile should extend far enough to adequately fit the background level.
  • pad_to_distance -- FilFinder only includes pixels in the radial profiles that are closest to the filament which avoids double-counting pixels. But if the filaments are closely packed together, this can severely limit the number of points used to make the profile. pad_to_distance forces all pixels within the given distance to be included in the profile. Must be given in length units and be less than max_dist.
  • use_longest_path -- Will use the longest path skeleton instead of the full skeleton. Default is False.
  • kwargs -- These are passed to the radial_profile function. Please see the documentation in the link for the different options. If an error about empty bins or an array with a shape of 0 is returned, try using auto_cut=False.

FilFinder supports 3 simple models for fitting the radial profiles: a Gaussian with a mean fixed to 0 and a constant background, the same Gaussian without a background, and a non-parametric method to estimate Gaussian widths. However, FilFinder uses the astropy.modeling package and will accept any 1D astropy model. For example, the radfil package has an astropy implementation of a Plummer model, which could be used here.

The parameters that control the fitting are:

  • fit_model -- The model to the profiles to. The defaults are gaussian_bkg, gaussian_nobkg, and nonparam. Otherwise, a 1D astropy model can be given, as discussed above.
  • fitter -- The fitter to use. See astropy.modeling.fitter. Defaults to a least-squares fitter.
  • try_nonparam -- If the fit to the model fails, the non-parametric method can be used instead. Default is True.
  • add_width_to_length -- The fitted FWHM can be added to the lengths (FilFinder2D.lengths), assuming that the skeleton's length was shortened by the width in the medial axis transform (FilFinder2D.medskel). The width will not be added if the fit was poor or highly unconstrained. Default is True.
  • deconvolve_width -- Subtract off the beam width when calculating the FWHM width. Default is True.
  • fwhm_function -- Pass a function that takes the fit_model and returns the FWHM and its uncertainty. If None is given, the FWHM is passed assuming a Gaussian profile.
  • chisq_max -- The critical reduced $\chi^2$ used to determine "bad" fits. The default is 10, and is entirely subjective. This seems to flag most bad fits, but the quality of the fits should always be visually checked.

With the default settings, a Gaussian with a constant background is fit to the profiles:


In [162]:
fil.find_widths(max_dist=0.2 * u.pc)
fil1.plot_radial_profile(xunit=u.pc)


/home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filament.py:927: UserWarning: Ignoring adding the width to the length because the fail flag was raised for the fit.
  warnings.warn("Ignoring adding the width to the length because"

The radial profile of the first filament is shown above. The binned radial profile is shown as black diamonds and the fit is shown with the red solid line.

The profile can be plotted with different xunits:


In [163]:
fil1.plot_radial_profile(xunit=u.pix)


Based on the warning above, at least one of the filament profile fits failed. We can look at the list of widths. FilFinder2D.widths() returns the FWHMs and their uncertainties:


In [164]:
fil.widths()


Out[164]:
(<Quantity [ 6.8908461 ,  4.3148627 ,  7.61875285,  5.73253972,  6.89457072,
             7.01281673, 43.24213164,  5.51177282,  8.63979524,  4.62431374,
             2.61879647,  5.11334422] pix>,
 <Quantity [ 0.21120144,  0.21981508,  0.25891586,  0.34575726,  0.1521915 ,
             0.57087374, 56.06296197,  0.18323136,  0.71417106,  0.42027546,
             0.39571059,  0.23258477] pix>)

These widths can be returned in other units as well:


In [165]:
fil.widths(u.pc)


Out[165]:
(<Quantity [0.05211611, 0.03263371, 0.05762134, 0.04335573, 0.05214428,
            0.05303859, 0.32704428, 0.04168605, 0.06534358, 0.03497412,
            0.0198062 , 0.0386727 ] pc>,
 <Quantity [0.00159734, 0.00166248, 0.0019582 , 0.00261499, 0.00115104,
            0.00431757, 0.42400942, 0.0013858 , 0.00540134, 0.00317858,
            0.0029928 , 0.00175906] pc>)

The 6th filament has a much larger width, and its uncertainty is very large. We can look at this radial profile more closely:


In [166]:
fil.filaments[6].plot_radial_profile(xunit=u.pc)


This is a faint feature near another, and the simple modeling has failed here. This is a case where fine-tuning may lead to a better result for certain filaments. See the Filament2D tutorial.

The fit results can be returned as an astropy table:


In [167]:
fil.width_fits(xunit=u.pc)


WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
Out[167]:
Table length=12
amplitude_0amplitude_0_errstddev_0stddev_0_erramplitude_1amplitude_1_errfwhmfwhm_errfail_flagmodel_type
pcpc
float64float64float64float64float64float64float64float64boolstr12
0.82663186429858810.0219085641566048423.38355777483813960.077567606223617970.070003010370571910.00668383566120528750.0521161128477118860.0015973362521571198Falsegaussian_bkg
1.6108891413215330.058239987652100952.4985814669369790.06845662834851640.101496946023255350.0112929165253715170.032633709699399960.0016624820013050273Falsegaussian_bkg
0.88674642961607710.0317938098345882943.65419149824438440.097349962145332650.15586020558325510.007850185128682180.057621339553875790.001958204832120392Falsegaussian_bkg
0.67151485116378430.0329015805386588262.9684377589379080.120413447445316150.11747743302849960.0091608069385849360.043355733443118290.0026149944928117643Falsegaussian_bkg
0.43534295389730920.0092042764814183523.38492580219815050.055902729829974550.025263666659104630.00223553440645752950.052144282479727330.0011510385041146299Falsegaussian_bkg
0.136535132476197640.0088534431439776053.42845234034586750.210580905555919930.046903771519111840.00380595340778858840.0530385881396014460.004317571408578928Falsegaussian_bkg
0.081398364901542660.0168410930782556530.139475842352161820.179294770629756460.041665405035018920.01833464391529560.32704428152714930.4240094190495937Truenonparam
0.213575915883271880.0053015602045502272.892051088736550.06297522337359680.056207419896808340.00131824034812510760.041686052740444440.0013857959315123815Falsegaussian_bkg
1.72929363629390490.168865793394973654.04312106970688840.27521588902838573-0.0105282189934598980.0330821245189028760.065343578593703850.005401342450237702Falsegaussian_bkg
1.60673535904577910.094761259356856162.5964907066197310.13498300783419730.209457581835425680.032079499453898730.0349741168810364450.003178582552243453Falsegaussian_bkg
0.038929220068031440.00286156556700467152.0303097125732020.092045312117276110.0290550985849673920.000332169175640486530.0198062023662416430.002992796127670664Falsegaussian_bkg
2.5447408559587630.082948212849974582.75690592768258420.077794503143038820.221283806317223830.0318516174623795640.0386727000824745250.001759060359710792Falsegaussian_bkg

This provides the fit results, the parameter errors, whether or not the fit failed, and the type of model used. The table can then be saved.

Other Filament Properties

With the width models, we can define other filament properties, such as the total intensity within the FWHM of the filament:


In [168]:
fil.total_intensity()


Out[168]:
$[7133.1509,~3538.8506,~1187.0249,~853.59998,~824.00342,~350.38428,~98.006683,~359.08759,~400.62476,~2129.0601,~11.555851,~1204.0597] \; \mathrm{K}$

If a background was fit in the model, the background level can be subtracted off. The index of the background parameter needs to be given. For the ''gaussian_bkg'', this is bkg_mod_index=2 and set as the default:


In [169]:
fil.total_intensity(bkg_subtract=True)


Out[169]:
$[6738.7539,~3378.1809,~990.95276,~738.23712,~782.19208,~266.42651,~21.717327,~291.18903,~404.1517,~1898.6567,~6.3840432,~1115.9888] \; \mathrm{K}$

The median brightness along the skeleton is calculated with:


In [170]:
fil.median_brightness()


Out[170]:
array([0.995523  , 2.0136652 , 1.088453  , 0.8972937 , 0.503168  ,
       0.1784499 , 0.05824468, 0.27832663, 0.9761157 , 1.8300676 ,
       0.06206793, 3.0347648 ], dtype=float32)

Based on the radial profile models, we can create an image based on the models:


In [171]:
fil_mod = fil.filament_model()
plt.imshow(fil_mod)
plt.colorbar()


Out[171]:
<matplotlib.colorbar.Colorbar at 0x7f2c63d729b0>

By default, the background level is subtracted off, as was used for fil.total_intensity. The maximum radius around each skeleton to evaluate the model can also be given with max_radius. The default is 3 times the FWHM, which should account for most of the model flux.

This model can be used to estimate the fraction of the total flux contained in the filamentary structure:


In [172]:
fil.covering_fraction()


Out[172]:
0.5562639744273953

The same keywords given to FilFinder2D.filament_model can be passed here.

The values aligned along the longest path are returned with:


In [173]:
profs = fil.ridge_profiles()
plt.subplot(211)
plt.plot(profs[0])
plt.subplot(212)
plt.plot(profs[1])
plt.tight_layout()


This can be useful for examining the distribution of cores along filaments.

Output Table & Images

FilFinder2D returns result tables as astropy tables. FilFinder2D.width_fits is highlighted above.

The width results and additional properties are returned with:


In [174]:
fil.output_table(xunit=u.pc)


WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
Out[174]:
Table length=12
lengthsbranchestotal_intensitymedian_brightnessX_posnY_posnamplitude_0amplitude_0_errstddev_0stddev_0_erramplitude_1amplitude_1_errfwhmfwhm_errfail_flagmodel_type
pcKpixpixpcpc
float64int64float64float32float64float64float64float64float64float64float64float64float64float64boolstr12
3.4839412361229956367133.150878906250.995523106.095.50.82663186429858810.0219085641566048423.38355777483813960.077567606223617970.070003010370571910.00668383566120528750.0521161128477118860.0015973362521571198Falsegaussian_bkg
2.080228297847961153538.85058593752.0136652188.045.01.6108891413215330.058239987652100952.4985814669369790.06845662834851640.101496946023255350.0112929165253715170.032633709699399960.0016624820013050273Falsegaussian_bkg
1.263693005023334831187.024902343751.088453142.0126.50.88674642961607710.0317938098345882943.65419149824438440.097349962145332650.15586020558325510.007850185128682180.057621339553875790.001958204832120392Falsegaussian_bkg
1.00870603098267567853.59997558593750.897293743.0123.00.67151485116378430.0329015805386588262.9684377589379080.120413447445316150.11747743302849960.0091608069385849360.043355733443118290.0026149944928117643Falsegaussian_bkg
1.234321265121755911824.003417968750.503168213.0120.50.43534295389730920.0092042764814183523.38492580219815050.055902729829974550.025263666659104630.00223553440645752950.052144282479727330.0011510385041146299Falsegaussian_bkg
1.10136027667576310350.384277343750.1784499218.0178.00.136535132476197640.0088534431439776053.42845234034586750.210580905555919930.046903771519111840.00380595340778858840.0530385881396014460.004317571408578928Falsegaussian_bkg
0.3487546458890681198.006683349609380.0582446851.0177.50.081398364901542660.0168410930782556530.139475842352161820.179294770629756460.041665405035018920.01833464391529560.32704428152714930.4240094190495937Truenonparam
1.295286248765398359.087585449218750.2783266323.0216.00.213575915883271880.0053015602045502272.892051088736550.06297522337359680.056207419896808340.00131824034812510760.041686052740444440.0013857959315123815Falsegaussian_bkg
0.309735174894607171400.6247558593750.976115783.0189.51.72929363629390490.168865793394973654.04312106970688840.27521588902838573-0.0105282189934598980.0330821245189028760.065343578593703850.005401342450237702Falsegaussian_bkg
0.90902980081439282129.060058593751.8300676185.0229.01.60673535904577910.094761259356856162.5964907066197310.13498300783419730.209457581835425680.032079499453898730.0349741168810364450.003178582552243453Falsegaussian_bkg
0.5649786406338143111.5558509826660160.06206793234.0227.00.038929220068031440.00286156556700467152.0303097125732020.092045312117276110.0290550985849673920.000332169175640486530.0198062023662416430.002992796127670664Falsegaussian_bkg
0.535556478610422831204.05969238281253.0347648156.0238.02.5447408559587630.082948212849974582.75690592768258420.077794503143038820.221283806317223830.0318516174623795640.0386727000824745250.001759060359710792Falsegaussian_bkg

This will include units if attached to the image or radial profile models.

The median positions can also be returned in world coordinates if WCS information was given:


In [175]:
fil.output_table(xunit=u.pc, world_coord=True)


WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match ( != pc).  Using pc for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match ( != K).  Using K for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_0_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'stddev_0_err' the 'unit' attribute does not match (pc != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
WARNING: MergeConflictWarning: In merged column 'amplitude_1_err' the 'unit' attribute does not match (K != ).  Using  for merged output [astropy.utils.metadata]
Out[175]:
Table length=12
lengthsbranchestotal_intensitymedian_brightnessRADecamplitude_0amplitude_0_errstddev_0stddev_0_erramplitude_1amplitude_1_errfwhmfwhm_errfail_flagmodel_type
pcKdegdegpcpc
float64int64float64float32float64float64float64float64float64float64float64float64float64float64boolstr12
3.4839412361229956367133.150878906250.9955230.05250000000105-0.03500000000070.82663186429858810.0219085641566048423.38355777483813960.077567606223617970.070003010370571910.00668383566120528750.0521161128477118860.0015973362521571198Falsegaussian_bkg
2.080228297847961153538.85058593752.01366520.13666666666940.10166666666871.6108891413215330.058239987652100952.4985814669369790.06845662834851640.101496946023255350.0112929165253715170.032633709699399960.0016624820013050273Falsegaussian_bkg
1.263693005023334831187.024902343751.0884530.000833333333350.0250000000004999970.88674642961607710.0317938098345882943.65419149824438440.097349962145332650.15586020558325510.007850185128682180.057621339553875790.001958204832120392Falsegaussian_bkg
1.00870603098267567853.59997558593750.89729370.0066666666668-0.14000000000280.67151485116378430.0329015805386588262.9684377589379080.120413447445316150.11747743302849960.0091608069385849360.043355733443118290.0026149944928117643Falsegaussian_bkg
1.234321265121755911824.003417968750.5031680.0108333333335499990.14333333333620.43534295389730920.0092042764814183523.38492580219815050.055902729829974550.025263666659104630.00223553440645752950.052144282479727330.0011510385041146299Falsegaussian_bkg
1.10136027667576310350.384277343750.1784499359.91499999999830.15166666666970.136535132476197640.0088534431439776053.42845234034586750.210580905555919930.046903771519111840.00380595340778858840.0530385881396014460.004317571408578928Falsegaussian_bkg
0.3487546458890681198.006683349609380.05824468359.91583333333165-0.126666666669199980.081398364901542660.0168410930782556530.139475842352161820.179294770629756460.041665405035018920.01833464391529560.32704428152714930.4240094190495937Truenonparam
1.295286248765398359.087585449218750.27832663359.8516666666637-0.173333333336799980.213575915883271880.0053015602045502272.892051088736550.06297522337359680.056207419896808340.00131824034812510760.041686052740444440.0013857959315123815Falsegaussian_bkg
0.309735174894607171400.6247558593750.9761157359.89583333333127-0.07333333333481.72929363629390490.168865793394973654.04312106970688840.27521588902838573-0.0105282189934598980.0330821245189028760.065343578593703850.005401342450237702Falsegaussian_bkg
0.90902980081439282129.060058593751.8300676359.82999999999660.09666666666861.60673535904577910.094761259356856162.5964907066197310.13498300783419730.209457581835425680.032079499453898730.0349741168810364450.003178582552243453Falsegaussian_bkg
0.5649786406338143111.5558509826660160.06206793359.833333333330.17833333333690.038929220068031440.00286156556700467152.0303097125732020.092045312117276110.0290550985849673920.000332169175640486530.0198062023662416430.002992796127670664Falsegaussian_bkg
0.535556478610422831204.05969238281253.0347648359.81499999999630.04833333333432.5447408559587630.082948212849974582.75690592768258420.077794503143038820.221283806317223830.0318516174623795640.0386727000824745250.001759060359710792Falsegaussian_bkg

A table for each of the branch properties of the filaments is returned with:


In [176]:
branch_tables = fil.branch_tables()
branch_tables[0]


Out[176]:
Table length=36
lengthintensity
pix
float64float32
62.28427124746193.9204443
20.485281374238571.1869783
6.2426406871192860.59490615
13.656854249492385.536005
7.00.9801212
20.0710678118654761.3630383
10.656854249492381.0734288
14.3137084989847610.56572586
12.828427124746192.260275
......
30.9705627484771432.326289
56.1126983722080940.77015895
14.8994949366116671.2642363
23.142135623730950.95123225
8.2426406871192860.5127995
36.3847763108502350.7552973
16.142135623730950.45914
34.556349186104040.34504774
11.828427124746190.3573772
51.355339059327380.3015641

If the RHT was run on branches, these data can also be added to the branch tables:


In [177]:
branch_tables = fil.branch_tables(include_rht=True)
branch_tables[0]


Out[177]:
Table length=36
lengthintensityorientationcurvature
pixradrad
float64float32float64float64
62.28427124746193.92044430.66059863599116160.9210741882576339
20.485281374238571.18697830.80137498937860440.923494292293757
6.2426406871192860.59490615-1.35083160090503380.43948899090739646
13.656854249492385.5360051.27521240615270750.5512162614209193
7.00.9801212-1.217704747523665e-160.320209242394435
20.0710678118654761.3630383-0.198719061950573460.468575647829656
10.656854249492381.07342881.49463577277433910.3891476232440563
14.3137084989847610.565725860.82559871617785240.46475729202323324
12.828427124746192.260275-0.156692319016146410.47639931090836724
............
30.9705627484771432.3262890.6302067516016830.8559321180406212
56.1126983722080940.77015895-0.084704450455261620.7149396089534249
14.8994949366116671.2642363-0.63743412247414090.44963034275443436
23.142135623730950.95123225-0.47490033664233770.5011542068110733
8.2426406871192860.51279950.21384665797580320.45210663309611
36.3847763108502350.75529730.20785632168662020.7597649003949792
16.142135623730950.45914-0.78266636157651130.4520951862182143
34.556349186104040.34504774-0.36715026357398991.4416318768939498
11.828427124746190.35737720.25364362496446790.45310231131019507
51.355339059327380.3015641-0.60666069522566310.532684628103508

These tables can be saved to a format supported by astropy tables.

Finally, the mask, skeletons, longest path skeletons, and the filament model can be saved as a FITS file:


In [178]:
fil.save_fits()


---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-178-328dcd15745e> in <module>()
----> 1 fil.save_fits()

~/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filfinder2D.py in save_fits(self, save_name, **kwargs)
   1333         out_hdu.append(model_hdu)
   1334 
-> 1335         out_hdu.writeto("{0}_image_output.fits".format(save_name))
   1336 
   1337     def save_stamp_fits(self, save_name=None, pad_size=20 * u.pix,

~/anaconda3/lib/python3.6/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    485                         # one with the name of the new argument to the function
    486                         kwargs[new_name[i]] = value
--> 487             return function(*args, **kwargs)
    488 
    489         return wrapper

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py in writeto(self, fileobj, output_verify, overwrite, checksum)
    865         # file object that's open to write only, or in append/update modes
    866         # but only if the file doesn't exist.
--> 867         fileobj = _File(fileobj, mode='ostream', overwrite=overwrite)
    868         hdulist = self.fromfile(fileobj)
    869         try:

~/anaconda3/lib/python3.6/site-packages/astropy/utils/decorators.py in wrapper(*args, **kwargs)
    485                         # one with the name of the new argument to the function
    486                         kwargs[new_name[i]] = value
--> 487             return function(*args, **kwargs)
    488 
    489         return wrapper

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in __init__(self, fileobj, mode, memmap, overwrite, cache)
    173             self._open_fileobj(fileobj, mode, overwrite)
    174         elif isinstance(fileobj, str):
--> 175             self._open_filename(fileobj, mode, overwrite)
    176         else:
    177             self._open_filelike(fileobj, mode, overwrite)

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in _open_filename(self, filename, mode, overwrite)
    515 
    516         if mode == 'ostream':
--> 517             self._overwrite_existing(overwrite, None, True)
    518 
    519         if os.path.exists(self.name):

~/anaconda3/lib/python3.6/site-packages/astropy/io/fits/file.py in _overwrite_existing(self, overwrite, fileobj, closed)
    412                     os.remove(self.name)
    413             else:
--> 414                 raise OSError("File {!r} already exists.".format(self.name))
    415 
    416     def _try_read_compressed(self, obj_or_name, magic, mode, ext=''):

OSError: File 'FilFinder_output_image_output.fits' already exists.

This will save the file with prefix given at the beginning. This can be changed by specifying save_name here. The keywords for FilFinder2D.filament_model can also be specified here.

The regions and stamps around each filament can also be saved with:


In [ ]:
fil.save_stamp_fits()

The same arguments for FilFinder2D.save_fits can be given, along with pad_size which sets how large of an area around each skeleton is included in the stamp. This will create a FITS file for each filament.