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 `__. .. code:: python %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. .. code:: python from fil_finder import FilFinder2D from astropy.io import fits hdu = fits.open("../examples/filaments_updatedhdr.fits")[0] fil = FilFinder2D(hdu) .. code:: python # HDU data as an array arr = hdu.data hdr = hdu.header fil = FilFinder2D(arr) .. parsed-literal:: /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: .. code:: python fil = FilFinder2D(arr, header=hdr) If `spectral-cube `__ is installed, the ``Projection`` or ``Slice`` classes can also be passed to ``FilFinder2D``: .. code:: python from spectral_cube import Projection proj = Projection.from_hdu(hdu) fil = FilFinder2D(proj) .. parsed-literal:: 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``: .. code:: python fil = FilFinder2D(hdu, distance=250 * u.pc) **Angular Scale** -- If no header information is given, the pixel-to-angular conversion can be given: .. code:: python fil = FilFinder2D(arr, ang_scale=0.1 * u.deg) .. parsed-literal:: /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: .. code:: python fil = FilFinder2D(hdu, beamwidth=10 * u.arcsec) **Custom Filament Masks** -- If you have a pre-computed filament mask, the mask array can be passed: .. code:: python # 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: .. code:: python 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: .. code:: python fil = FilFinder2D(hdu, distance=260 * u.pc) The beamwidth is :math:`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: .. code:: python fil.preprocess_image(flatten_percent=95) .. code:: python 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() .. image:: tutorial_files/tutorial_26_0.png 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 :math:`\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: .. code:: python 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. .. code:: python 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: .. code:: python fil.create_mask(verbose=True) .. image:: tutorial_files/tutorial_32_0.png 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``: .. code:: python fil.create_mask(verbose=True, border_masking=False) .. image:: tutorial_files/tutorial_34_0.png 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: .. code:: python fil.create_mask(verbose=True, border_masking=False, size_thresh=400 * u.pix**2) .. image:: tutorial_files/tutorial_36_0.png 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: .. code:: python # 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') .. parsed-literal:: .. image:: tutorial_files/tutorial_38_1.png The dark regions will be excluded from the final mask. The filament mask with the threshold is then: .. code:: python fil.create_mask(verbose=True, border_masking=False, size_thresh=400 * u.pix**2, glob_thresh=0.0267) .. image:: tutorial_files/tutorial_40_0.png 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: .. code:: python fil.create_mask(use_existing_mask=True) .. parsed-literal:: /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: .. code:: python fil.medskel(verbose=True) .. image:: tutorial_files/tutorial_44_0.png 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: .. code:: python fil.analyze_skeletons() plt.imshow(fil.skeleton, origin='lower') .. parsed-literal:: .. image:: tutorial_files/tutorial_46_1.png 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``. .. code:: python fil.analyze_skeletons(branch_thresh=40 * u.pix, prune_criteria='length') plt.imshow(fil.skeleton, origin='lower') plt.contour(fil.skeleton_longpath, colors='r') .. parsed-literal:: .. image:: tutorial_files/tutorial_48_1.png 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: .. code:: python fil.analyze_skeletons(branch_thresh=400 * u.pix, prune_criteria='length') plt.imshow(fil.skeleton, origin='lower') plt.contour(fil.skeleton_longpath, colors='r') .. parsed-literal:: .. image:: tutorial_files/tutorial_50_1.png 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: .. code:: python 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') .. parsed-literal:: .. image:: tutorial_files/tutorial_52_1.png Another new feature of ``FilFinder2D`` is that each filament has its own analysis class defined in ``fil.filaments``: .. code:: python fil.filaments .. parsed-literal:: [, , , , , , , , , , , ] 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: .. code:: python fil1 = fil.filaments[0] fil1.skeleton_analysis(fil.image, verbose=True, branch_thresh=5 * u.pix, prune_criteria='length') .. image:: tutorial_files/tutorial_56_0.png .. image:: tutorial_files/tutorial_56_1.png .. image:: tutorial_files/tutorial_56_2.png 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: .. code:: python fil.lengths() .. math:: [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: .. code:: python fil.lengths(u.deg) .. math:: [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}} .. code:: python fil.lengths(u.pc) .. math:: [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: .. code:: python fil.branch_properties.keys() .. parsed-literal:: dict_keys(['length', 'intensity', 'pixels', 'number']) .. code:: python fil.branch_properties['number'] .. parsed-literal:: array([36, 15, 3, 7, 11, 10, 1, 8, 1, 8, 1, 3]) .. code:: python fil.branch_properties['length'][0] .. math:: [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: .. code:: python fil.branch_lengths(u.pix)[0] .. math:: [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} .. code:: python fil.branch_lengths(u.pc)[0] .. math:: [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: .. code:: python fil.exec_rht() fil1.plot_rht_distrib() .. image:: tutorial_files/tutorial_70_0.png 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 :math:`2\theta`. Since there is no preferred direction, :math:`0` and :math:`\pi` are equivalent direction for a filament, and so the distribution is defined over :math:`\theta \in [-\pi/2, \pi/2)`. Plotting the distribution as :math:`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 :math:`\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 :math:`\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. .. code:: python fil.orientation .. math:: [-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} .. code:: python fil.curvature .. math:: [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``: .. code:: python 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: .. code:: python _ = plt.hist(fil.orientation_branches[0].value, bins=10) plt.xlabel("Orientation (rad)") .. parsed-literal:: Text(0.5,0,'Orientation (rad)') .. image:: tutorial_files/tutorial_77_1.png .. code:: python 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)") .. parsed-literal:: Text(0.5,0,'Orientation (rad)') .. image:: tutorial_files/tutorial_78_1.png 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 :math:`\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: .. code:: python fil.find_widths(max_dist=0.2 * u.pc) fil1.plot_radial_profile(xunit=u.pc) .. parsed-literal:: /home/eric/Dropbox/code_development/filaments/build/lib.linux-x86_64-3.6/fil_finder/filament.py:926: 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" .. image:: tutorial_files/tutorial_81_1.png 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 ``xunit``\ s: .. code:: python fil1.plot_radial_profile(xunit=u.pix) .. image:: tutorial_files/tutorial_83_0.png 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: .. code:: python fil.widths() .. parsed-literal:: (, ) These widths can be returned in other units as well: .. code:: python fil.widths(u.pc) .. parsed-literal:: (, ) The 6th filament has a much larger width, and its uncertainty is very large. We can look at this radial profile more closely: .. code:: python fil.filaments[6].plot_radial_profile(xunit=u.pc) .. image:: tutorial_files/tutorial_89_0.png 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 `__: .. code:: python fil.width_fits(xunit=u.pc) .. parsed-literal:: 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] .. raw:: html 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: .. code:: python fil.total_intensity() .. math:: [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: .. code:: python fil.total_intensity(bkg_subtract=True) .. math:: [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: .. code:: python fil.median_brightness() .. parsed-literal:: 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: .. code:: python fil_mod = fil.filament_model() plt.imshow(fil_mod) plt.colorbar() .. parsed-literal:: .. image:: tutorial_files/tutorial_99_1.png 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: .. code:: python fil.covering_fraction() .. parsed-literal:: 0.5562639744273953 The same keywords given to ``FilFinder2D.filament_model`` can be passed here. The values aligned along the longest path are returned with: .. code:: python profs = fil.ridge_profiles() plt.subplot(211) plt.plot(profs[0]) plt.subplot(212) plt.plot(profs[1]) plt.tight_layout() .. image:: tutorial_files/tutorial_103_0.png 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: .. code:: python fil.output_table(xunit=u.pc) .. parsed-literal:: 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] .. raw:: html 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: .. code:: python fil.output_table(xunit=u.pc, world_coord=True) .. parsed-literal:: 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] .. raw:: html 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: .. code:: python branch_tables = fil.branch_tables() branch_tables[0] .. raw:: html 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: .. code:: python branch_tables = fil.branch_tables(include_rht=True) branch_tables[0] .. raw:: html 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: .. code:: python 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: .. code:: python 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.