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)
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)
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)
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)
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.
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)
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]:
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)
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)
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]:
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]:
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]:
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]:
Another new feature of FilFinder2D is that each filament has its own analysis class defined in fil.filaments:
In [146]:
fil.filaments
Out[146]:
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:
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]:
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]:
In [150]:
fil.lengths(u.pc)
Out[150]:
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]:
In [152]:
fil.branch_properties['number']
Out[152]:
In [153]:
fil.branch_properties['length'][0]
Out[153]:
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]:
In [155]:
fil.branch_lengths(u.pc)[0]
Out[155]:
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]:
In [158]:
fil.curvature
Out[158]:
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]:
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]:
No orientation is strongly preferred in the example data.
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)
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]:
These widths can be returned in other units as well:
In [165]:
fil.widths(u.pc)
Out[165]:
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)
Out[167]:
In [168]:
fil.total_intensity()
Out[168]:
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]:
The median brightness along the skeleton is calculated with:
In [170]:
fil.median_brightness()
Out[170]:
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]:
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]:
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.
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)
Out[174]:
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)
Out[175]:
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]:
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]:
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()
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.