Merge branch 'master' into FilFinder3D

Author: e-koch

.github/workflows/main.yml

@@ -21,49 +21,49 @@ jobs:
             toxenv: py39-test
 
           - os: ubuntu-latest
-            python-version: 3.8
-            name: Python 3.8 with minimal dependencies
-            toxenv: py38-test
+            python-version: '3.10'
+            name: Python 3.10 with minimal dependencies
+            toxenv: py310-test
 
           - os: ubuntu-latest
-            python-version: 3.7
-            name: Python 3.7 with minimal dependencies
-            toxenv: py37-test
+            python-version: 3.11
+            name: Python 3.11 with minimal dependencies
+            toxenv: py311-test
 
           - os: ubuntu-latest
             python-version: 3.9
             name: Python 3.9 with minimal and dev dependencies
             toxenv: py39-test-dev
 
           - os: ubuntu-latest
-            python-version: 3.9
-            name: Python 3.9, all dependencies, and dev versions of key dependencies
-            toxenv: py39-test-dev
+            python-version: 3.11
+            name: Python 3.11, all dependencies, and dev versions of key dependencies
+            toxenv: py311-test-dev
 
           - os: ubuntu-latest
-            python-version: 3.9
-            name: Python 3.9, all dependencies, and coverage
-            toxenv: py39-test-all-cov
+            python-version: 3.11
+            name: Python 3.11, all dependencies, and coverage
+            toxenv: py311-test-all-cov
 
           - os: macos-latest
-            python-version: 3.9
-            name: Python 3.9 with all dependencies on MacOS X
-            toxenv: py39-test-all-dev
+            python-version: 3.11
+            name: Python 3.11 with all dependencies on MacOS X
+            toxenv: py311-test-all-dev
 
           # - os: windows-latest
           #   python-version: 3.7
           #   name: Python 3.7, with all dependencies on Windows
           #   toxenv: py37-test-all-dev
 
           - os: ubuntu-latest
-            python-version: 3.9
+            python-version: 3.11
             name: Documentation
             toxenv: build_docs
 
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v2
+      uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install testing dependencies
@@ -72,6 +72,6 @@ jobs:
       run: tox -v -e ${{ matrix.toxenv }}
     - name: Upload coverage to codecov
       if: ${{ contains(matrix.toxenv,'-cov') }}
-      uses: codecov/codecov-action@v1.0.13
+      uses: codecov/codecov-action@v3
       with:
         file: ./coverage.xml

.github/workflows/publish.yml

@@ -7,11 +7,11 @@ jobs:
     name: Build source distribution
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: actions/setup-python@v2
+      - uses: actions/checkout@v3
+      - uses: actions/setup-python@v4
         name: Install Python
         with:
-          python-version: '3.9'
+          python-version: '3.11'
       - name: Install build
         run: python -m pip install build
       - name: Build sdist
@@ -26,11 +26,11 @@ jobs:
     runs-on: ubuntu-latest
     if: github.event_name == 'push' && startsWith(github.event.ref, 'refs/tags/v')
     steps:
-      - uses: actions/download-artifact@v2
+      - uses: actions/download-artifact@v4
         with:
           name: artifact
           path: dist
-      - uses: pypa/gh-action-pypi-publish@master
+      - uses: pypa/gh-action-pypi-publish@release/v1
         with:
           user: __token__
           password: ${{ secrets.pypi_password }}

.readthedocs.yml

@@ -1,14 +1,14 @@
 version: 2
 
 build:
-  image: latest
+  os: "ubuntu-20.04"
+  tools:
+    python: "3.12"
 
 # Install regular dependencies.
 python:
-  version: 3.7
   install:
     - method: pip
       path: .
       extra_requirements:
         - docs
-

fil_finder/__init__.py

@@ -2,7 +2,6 @@
 
 from ._astropy_init import __version__, test
 
-from .filfind_class import fil_finder_2D
 from .filfinder2D import FilFinder2D
 from .filfinderPPV import FilFinderPPV
 from .filfinderPPP import FilFinderPPP

fil_finder/conftest.py

@@ -3,7 +3,8 @@
 # no matter how it is invoked within the source tree.
 from __future__ import print_function, absolute_import, division
 
-import pytest
+import os
+from setuptools._distutils.version import LooseVersion
 
 from astropy.version import version as astropy_version
 
@@ -13,8 +14,21 @@
 else:
     from pytest_astropy_header.display import PYTEST_HEADER_MODULES, TESTED_VERSIONS
 
+import pytest
+from fil_finder.tests.testing_utils import generate_filament_model
+
+
 def pytest_configure(config):
 
     config.option.astropy_header = True
 
     PYTEST_HEADER_MODULES['Astropy'] = 'astropy'
+
+
+@pytest.fixture
+def simple_filament_model():
+
+    mod = generate_filament_model(return_hdu=True, pad_size=31, shape=150,
+                                  width=10., background=0.1)[0]
+
+    yield mod

fil_finder/filament.py

@@ -16,9 +16,9 @@
 else:
     import cPickle as pickle
 
-from .length import (init_lengths, main_length, make_final_skeletons,
-                     pre_graph, longest_path, prune_graph, all_shortest_paths)
-from .pixel_ident import pix_identify
+from .length import (init_lengths, main_length,
+                     pre_graph, longest_path, prune_graph)
+from .pixel_ident import pix_identify, make_final_skeletons
 from .utilities import pad_image, in_ipynb, red_chisq
 from .base_conversions import UnitConverter
 from .rollinghough import rht
@@ -200,7 +200,8 @@ def skeleton(self, pad_size=0, corner_pix=None, out_type='all'):
     def skeleton_analysis(self, image, verbose=False, save_png=False,
                           save_name=None, prune_criteria='all',
                           relintens_thresh=0.2, max_prune_iter=10,
-                          branch_thresh=0 * u.pix):
+                          branch_thresh=0 * u.pix,
+                          return_self=False):
         '''
         Run the skeleton analysis.
 
@@ -230,6 +231,9 @@ def skeleton_analysis(self, image, verbose=False, save_png=False,
             Maximum number of pruning iterations to apply.
         branch_thresh : `~astropy.units.Quantity`, optional
             Minimum length for a branch to be eligible to be pruned.
+        return_self : bool, optional
+            Return the Filament2D object after skeleton analysis. This is needed
+            for parallel processing in FilFinder2D.
         '''
 
         # NOTE:
@@ -406,6 +410,9 @@ def skeleton_analysis(self, image, verbose=False, save_png=False,
              'number': branch_properties['number'][0],
              'pixels': branch_properties['pixels'][0]}
 
+        if return_self:
+            return self
+
     @property
     def branch_properties(self):
         '''
@@ -497,7 +504,8 @@ def plot_graph(self, save_name=None, layout_func=nx.spring_layout):
         # Add in the ipynb checker
 
     def rht_analysis(self, radius=10 * u.pix, ntheta=180,
-                     background_percentile=25):
+                     background_percentile=25,
+                     return_self=False):
         '''
         Use the RHT to find the filament orientation and dispersion of the
         longest path.
@@ -512,6 +520,8 @@ def rht_analysis(self, radius=10 * u.pix, ntheta=180,
         background_percentile : float, optional
             Float between 0 and 100 that sets a background level for the RHT
             distribution before calculating orientation and curvature.
+        return_self : bool, optional
+            If True, return the `Filament2D` object. Defaults to False.
         '''
 
         if not hasattr(radius, 'unit'):
@@ -539,6 +549,9 @@ def rht_analysis(self, radius=10 * u.pix, ntheta=180,
         self._orientation_hist = [theta, R]
         self._orientation_quantiles = [twofive, sevenfive]
 
+        if return_self:
+            return self
+
     @property
     def orientation_hist(self):
         '''
@@ -603,9 +616,10 @@ def plot_rht_distrib(self, save_name=None):
         else:
             plt.show()
 
-    def rht_branch_analysis(self, radius=10 * u.pix, ntheta=180,
+    def rht_branch_analysis(self,radius=10 * u.pix, ntheta=180,
                             background_percentile=25,
-                            min_branch_length=3 * u.pix):
+                            min_branch_length=3 * u.pix,
+                            return_self=False):
         '''
         Use the RHT to find the filament orientation and dispersion of each
         branch in the filament.
@@ -623,6 +637,8 @@ def rht_branch_analysis(self, radius=10 * u.pix, ntheta=180,
         min_branch_length : `~astropy.units.Quantity`, optional
             Minimum length of a branch to run the RHT on. Branches that are
             too short will cause spikes along the axis angles or 45 deg. off.
+        return_self : bool, optional
+            Return the `Filament2D` object with the results.
         '''
 
         # Convert length cut to pixel units
@@ -680,6 +696,9 @@ def rht_branch_analysis(self, radius=10 * u.pix, ntheta=180,
         self._orientation_branches = np.array(means) * u.rad
         self._curvature_branches = np.array(iqrs) * u.rad
 
+        if return_self:
+            return self
+
     @property
     def orientation_branches(self):
         '''
@@ -706,6 +725,7 @@ def width_analysis(self, image, all_skeleton_array=None,
                        beamwidth=None,
                        fwhm_function=None,
                        chisq_max=10.,
+                       return_self=False,
                        **kwargs):
         '''
 
@@ -756,6 +776,8 @@ def width_analysis(self, image, all_skeleton_array=None,
         chisq_max : float, optional
             Enable the fail flag if the reduced chi-squared value is above
             this limit.
+        return_self : bool, optional
+            Return the `Filament2D` object if True. Defaults to False.
         kwargs : Passed to `~fil_finder.width.radial_profile`.
 
         '''
@@ -1018,6 +1040,9 @@ def width_analysis(self, image, all_skeleton_array=None,
 
         self._radprof_failflag = fail_flag
 
+        if return_self:
+            return self
+
     @property
     def radprof_fit_fail_flag(self):
         '''
@@ -1122,17 +1147,25 @@ def radprof_model(self):
         '''
         return self._radprof_model
 
-    def plot_radial_profile(self, save_name=None, xunit=u.pix,
-                            ax=None):
+    def plot_radial_profile(self,
+                            ax=None,
+                            save_name=None,
+                            show_plot=True,
+                            xunit=u.pix
+                            ):
         '''
         Plot the radial profile of the filament and the fitted model.
 
         Parameters
         ----------
-        xunit : `~astropy.units.Unit`, optional
-            Pixel, angular, or physical unit to convert to.
         ax : `~matplotlib.axes`, optional
             Use an existing set of axes to plot the profile.
+        save_name : str, optional
+            Name of saved plot. A plot is only saved if a name is given.        
+        show_plot : bool, optional
+            Display open figure.        
+        xunit : `~astropy.units.Unit`, optional
+            Pixel, angular, or physical unit to convert to.
         '''
 
         dist, radprof = self.radprofile
@@ -1164,8 +1197,12 @@ def plot_radial_profile(self, save_name=None, xunit=u.pix,
 
         if save_name is not None:
             plt.savefig(save_name)
-
-        plt.show()
+            if not show_plot:
+                plt.close()
+            
+        if show_plot:
+            plt.show()
+            
         if in_ipynb():
             plt.clf()
 
@@ -1459,7 +1496,10 @@ def branch_table(self, include_rht=False):
                     names=branch_data.keys())
         return tab
 
-    def save_fits(self, savename, image, pad_size=20 * u.pix, header=None,
+    def save_fits(self, savename, image,
+                  image_dict=None,
+                  pad_size=20 * u.pix,
+                  header=None,
                   model_kwargs={},
                   **kwargs):
         '''
@@ -1468,6 +1508,8 @@ def save_fits(self, savename, image, pad_size=20 * u.pix, header=None,
 
         Parameters
         ----------
+        savename : str
+            Filename to save to.
         image : `~numpy.ndarray` or `~astropy.units.Quantity`
             The image from which the filament was extracted.
         pad_size : `~astropy.units.Quantity`, optional
@@ -1509,24 +1551,57 @@ def save_fits(self, savename, image, pad_size=20 * u.pix, header=None,
             else:
                 header = fits.Header()
 
+        # Add the pixel extents into the header
+        from astropy.table import Table, Column
+
+        tab = Table()
+        tab.add_column(Column([self.pixel_extents[0][0],
+                               self.pixel_extents[1][0]],
+                              name='lower_coord'))
+        tab.add_column(Column([self.pixel_extents[0][1],
+                               self.pixel_extents[1][1]],
+                              name='upper_coord'))
+
+
         # Strip off units if the image is a Quantity
         if hasattr(input_image, 'unit'):
             input_image = input_image.value.copy()
 
         hdu = fits.PrimaryHDU(input_image, header)
+        hdu.name = 'IMAGE'
 
         skel_hdr = header.copy()
         skel_hdr['BUNIT'] = ("", "bool")
         skel_hdr['COMMENT'] = "Skeleton created by fil_finder on " + \
             time.strftime("%c")
 
         skel_hdu = fits.ImageHDU(skels.astype(int), skel_hdr)
+        skel_hdu.name = 'SKELETON'
 
         skel_lp_hdu = fits.ImageHDU(skels_lp.astype(int), skel_hdr)
+        skel_lp_hdu.name = 'SKELETON_LONGPATH'
 
         model_hdu = fits.ImageHDU(model, header)
+        model_hdu.name = 'MODEL'
+
+        tab_hdu = fits.table_to_hdu(tab)
+        tab_hdu.name = 'PIXEXTENTS'
+
+        hdulist = fits.HDUList([hdu, skel_hdu, skel_lp_hdu, model_hdu, tab_hdu])
+
+        # If image_dict is provided, save cutouts from the image list
+        if image_dict is not None:
+            for key in image_dict:
+                img = image_dict[key]
+                img = pad_image(img, self.pixel_extents, pad_size)
+                if img.shape != skels.shape:
+                    img = self.image_slicer(img, skels.shape,
+                                            pad_size=pad_size)
+
+                img_hdu = fits.ImageHDU(img, header)
+                img_hdu.name = key.upper()
 
-        hdulist = fits.HDUList([hdu, skel_hdu, skel_lp_hdu, model_hdu])
+                hdulist.append(img_hdu)
 
         hdulist.writeto(savename, **kwargs)
 

fil_finder/filfind_class.py

@@ -14,6 +14,7 @@
 import os
 import time
 import warnings
+import concurrent.futures
 
 from .pixel_ident import recombine_skeletons, isolateregions
 from .utilities import eight_con, round_to_odd, threshold_local, in_ipynb
@@ -498,7 +499,7 @@ def create_mask(self, glob_thresh=None, adapt_thresh=None,
             if in_ipynb():
                 p.clf()
 
-    def medskel(self, verbose=False, save_png=False):
+    def medskel(self, verbose=False, save_png=False, rng=None):
         '''
         This function performs the medial axis transform (skeletonization)
         on the mask. This is essentially a wrapper function of
@@ -516,6 +517,9 @@ def medskel(self, verbose=False, save_png=False):
             Enables plotting.
         save_png : bool, optional
             Saves the plot made in verbose mode. Disabled by default.
+        rng : numpy.random.RandomState or int, optional
+            Random number generator for reproducibility. Used for tie breaks in
+            the `medial_axis <https://scikit-image.org/docs/stable/api/skimage.morphology.html#skimage.morphology.medial_axis>`_ function.
 
         Attributes
         ----------
@@ -525,8 +529,17 @@ def medskel(self, verbose=False, save_png=False):
             The distance transform used to create the skeletons.
         '''
 
-        self.skeleton, self.medial_axis_distance = \
-            medial_axis(self.mask, return_distance=True)
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # The kwarg for rng has changed in different skimage versions.
+        try:
+            self.skeleton, self.medial_axis_distance = \
+                medial_axis(self.mask, return_distance=True, rng=rng)
+        except TypeError:
+            self.skeleton, self.medial_axis_distance = \
+                medial_axis(self.mask, return_distance=True, random_state=rng)
+
         self.medial_axis_distance = \
             self.medial_axis_distance * self.skeleton * u.pix
         # Delete connection smaller than 2 pixels wide. Such a small
@@ -551,11 +564,18 @@ def medskel(self, verbose=False, save_png=False):
             if in_ipynb():
                 p.clf()
 
-    def analyze_skeletons(self, prune_criteria='all', relintens_thresh=0.2,
-                          nbeam_lengths=5, branch_nbeam_lengths=3,
-                          skel_thresh=None, branch_thresh=None,
+    def analyze_skeletons(self,
+                          nthreads=1,
+                          prune_criteria='all',
+                          relintens_thresh=0.2,
+                          nbeam_lengths=5,
+                          branch_nbeam_lengths=3,
+                          skel_thresh=None,
+                          branch_thresh=None,
                           max_prune_iter=10,
-                          verbose=False, save_png=False, save_name=None):
+                          verbose=False,
+                          save_png=False,
+                          save_name=None):
         '''
 
         Prune skeleton structure and calculate the branch and longest-path
@@ -564,6 +584,8 @@ def analyze_skeletons(self, prune_criteria='all', relintens_thresh=0.2,
 
         Parameters
         ----------
+        nthreads : int, optional
+            Number of threads to use to parallelize the skeleton analysis.
         prune_criteria : {'all', 'intensity', 'length'}, optional
             Choose the property to base pruning on. 'all' requires that the
             branch fails to satisfy the length and relative intensity checks.
@@ -638,25 +660,26 @@ def analyze_skeletons(self, prune_criteria='all', relintens_thresh=0.2,
             # Relabel after deleting short skeletons.
             labels, num = nd.label(self.skeleton, eight_con())
 
+
         self.filaments = [Filament2D(np.where(labels == lab),
                                      converter=self.converter) for lab in
                           range(1, num + 1)]
 
-        self.number_of_filaments = num
-
         # Now loop over the skeleton analysis for each filament object
-        for n, fil in enumerate(self.filaments):
-            savename = "{0}_{1}".format(save_name, n)
-            if verbose:
-                print("Filament: %s / %s" % (n + 1, self.number_of_filaments))
+        with concurrent.futures.ProcessPoolExecutor(nthreads) as executor:
+            futures = [executor.submit(fil.skeleton_analysis, self.image,
+                                                             verbose=verbose,
+                                                             save_png=save_png,
+                                                             save_name=save_name,
+                                                             prune_criteria=prune_criteria,
+                                                             relintens_thresh=relintens_thresh,
+                                                             branch_thresh=self.branch_thresh,
+                                                             max_prune_iter=max_prune_iter,
+                                                             return_self=True)
+                       for fil in self.filaments]
+            self.filaments = [future.result() for future in futures]
 
-            fil.skeleton_analysis(self.image, verbose=verbose,
-                                  save_png=save_png,
-                                  save_name=savename,
-                                  prune_criteria=prune_criteria,
-                                  relintens_thresh=relintens_thresh,
-                                  branch_thresh=self.branch_thresh,
-                                  max_prune_iter=max_prune_iter)
+        self.number_of_filaments = num
 
         self.array_offsets = [fil.pixel_extents for fil in self.filaments]
 
@@ -749,7 +772,9 @@ def end_pts(self):
         '''
         return [fil.end_pts for fil in self.filaments]
 
-    def exec_rht(self, radius=10 * u.pix,
+    def exec_rht(self,
+                 nthreads=1,
+                 radius=10 * u.pix,
                  ntheta=180, background_percentile=25,
                  branches=False, min_branch_length=3 * u.pix,
                  verbose=False, save_png=False, save_name=None):
@@ -774,6 +799,8 @@ def exec_rht(self, radius=10 * u.pix,
 
         Parameters
         ----------
+        nthreads : int, optional
+            The number of threads to use.
         radius : int
             Sets the patch size that the RHT uses.
         ntheta : int, optional
@@ -792,6 +819,7 @@ def exec_rht(self, radius=10 * u.pix,
         save_name : str, optional
             Prefix for the saved plots.
 
+
         Attributes
         ----------
         rht_curvature : dict
@@ -813,23 +841,35 @@ def exec_rht(self, radius=10 * u.pix,
         if save_name is None:
             save_name = self.save_name
 
-        for n, fil in enumerate(self.filaments):
-            if verbose:
-                print("Filament: %s / %s" % (n + 1, self.number_of_filaments))
 
-            if branches:
-                fil.rht_branch_analysis(radius=radius,
+        if branches:
+            with concurrent.futures.ProcessPoolExecutor(nthreads) as executor:
+                futures = [executor.submit(fil.rht_branch_analysis,
+                                           radius=radius,
                                         ntheta=ntheta,
                                         background_percentile=background_percentile,
-                                        min_branch_length=min_branch_length)
+                                        min_branch_length=min_branch_length,
+                                        return_self=True)
+                        for fil in self.filaments]
+                self.filaments = [future.result() for future in futures]
 
-            else:
-                fil.rht_analysis(radius=radius, ntheta=ntheta,
-                                 background_percentile=background_percentile)
 
-                if verbose:
+        else:
+            with concurrent.futures.ProcessPoolExecutor(nthreads) as executor:
+                futures = [executor.submit(fil.rht_analysis,
+                                           radius=radius,
+                                        ntheta=ntheta,
+                                        background_percentile=background_percentile,
+                                        return_self=True)
+                        for fil in self.filaments]
+                self.filaments = [future.result() for future in futures]
+
+
+            if verbose:
+                for n, fil in enumerate(self.filaments):
+
                     if save_png:
-                        savename = "{0}_{1}_rht.png".format(save_name, n)
+                        save_name = "{0}_{1}_rht.png".format(save_name, n)
                     else:
                         save_name = None
                     fil.plot_rht_distrib(save_name=save_name)
@@ -886,7 +926,9 @@ def pre_recombine_mask_corners(self):
         '''
         return self._pre_recombine_mask_corners
 
-    def find_widths(self, max_dist=10 * u.pix,
+    def find_widths(self,
+                    nthreads=1,
+                    max_dist=10 * u.pix,
                     pad_to_distance=0 * u.pix,
                     fit_model='gaussian_bkg',
                     fitter=None,
@@ -915,12 +957,8 @@ def find_widths(self, max_dist=10 * u.pix,
 
         Parameters
         ----------
-        image : `~astropy.unit.Quantity` or `~numpy.ndarray`
-            The image from which the filament was extracted.
-        all_skeleton_array : np.ndarray
-            An array with the skeletons of other filaments. This is used to
-            avoid double-counting pixels in the radial profiles in nearby
-            filaments.
+        nthreads : int, optional
+            Number of threads to use.
         max_dist : `~astropy.units.Quantity`, optional
             Largest radius around the skeleton to create the profile from. This
             can be given in physical, angular, or physical units.
@@ -967,23 +1005,26 @@ def find_widths(self, max_dist=10 * u.pix,
         if save_name is None:
             save_name = self.save_name
 
-        for n, fil in enumerate(self.filaments):
+        with concurrent.futures.ProcessPoolExecutor(nthreads) as executor:
+            futures = [executor.submit(fil.width_analysis, self.image,
+                                       all_skeleton_array=self.skeleton,
+                                       max_dist=max_dist,
+                                       pad_to_distance=pad_to_distance,
+                                       fit_model=fit_model,
+                                       fitter=fitter, try_nonparam=try_nonparam,
+                                       use_longest_path=use_longest_path,
+                                       add_width_to_length=add_width_to_length,
+                                       deconvolve_width=deconvolve_width,
+                                       beamwidth=self.beamwidth,
+                                       fwhm_function=fwhm_function,
+                                       chisq_max=chisq_max,
+                                       return_self=True,
+                                       **kwargs)
+                       for fil in self.filaments]
+            self.filaments = [future.result() for future in futures]
 
-            if verbose:
-                print("Filament: %s / %s" % (n + 1, self.number_of_filaments))
-
-            fil.width_analysis(self.image, all_skeleton_array=self.skeleton,
-                               max_dist=max_dist,
-                               pad_to_distance=pad_to_distance,
-                               fit_model=fit_model,
-                               fitter=fitter, try_nonparam=try_nonparam,
-                               use_longest_path=use_longest_path,
-                               add_width_to_length=add_width_to_length,
-                               deconvolve_width=deconvolve_width,
-                               beamwidth=self.beamwidth,
-                               fwhm_function=fwhm_function,
-                               chisq_max=chisq_max,
-                               **kwargs)
+
+        for n, fil in enumerate(self.filaments):
 
             if verbose:
                 if save_png:
@@ -1409,7 +1450,10 @@ def save_fits(self, save_name=None,
         out_hdu.writeto("{0}_image_output.fits".format(save_name),
                         **kwargs)
 
-    def save_stamp_fits(self, save_name=None, pad_size=20 * u.pix,
+    def save_stamp_fits(self,
+                        image_dict=None,
+                        save_name=None,
+                        pad_size=20 * u.pix,
                         model_kwargs={},
                         **kwargs):
         '''
@@ -1421,6 +1465,10 @@ def save_stamp_fits(self, save_name=None, pad_size=20 * u.pix,
 
         Parameters
         ----------
+        image_dict : dict, optional
+            Dictionary of arrays to save matching the pixel extents of each filament.
+            The shape of each array *must* be the same shape as the original image
+            given to `~FilFinder2D`.
         save_name : str, optional
             The prefix for the saved file. If None, the save name specified
             when `~FilFinder2D` was first called.
@@ -1436,10 +1484,21 @@ def save_stamp_fits(self, save_name=None, pad_size=20 * u.pix,
         else:
             save_name = os.path.splitext(save_name)[0]
 
+        if image_dict is not None:
+            for ii, key in enumerate(image_dict):
+                this_image = image_dict[key]
+                if this_image.shape != self.image.shape:
+                    raise ValueError("All images in image_dict must be same shape as fil.image. "
+                                     f"For index {ii}, found shape {this_image.shape} not {self.image.shape}")
+
+
         for n, fil in enumerate(self.filaments):
 
-            savename = "{0}_stamp_{1}.fits".format(save_name, n)
+            savename = f"{save_name}_stamp_{n}.fits"
 
-            fil.save_fits(savename, self.image, pad_size=pad_size,
+            fil.save_fits(savename,
+                          self.image,
+                          image_dict=image_dict,
+                          pad_size=pad_size,
                           model_kwargs=model_kwargs,
                           **kwargs)

fil_finder/filfinder2D.py

@@ -1,8 +1,6 @@
 # Licensed under an MIT open source license - see LICENSE
 
 from .utilities import *
-from .pixel_ident import *
-
 
 import numpy as np
 import scipy.ndimage as nd
@@ -830,3 +828,23 @@ def get_weight(pat):
     long_path_length = max(all_weights)
 
     return long_path, long_path_length
+
+  
+def merge_nodes(node, G):
+    '''
+    Combine a node into its neighbors.
+    '''
+
+    neigb = list(G[node])
+
+    if len(neigb) != 2:
+        return G
+
+    new_weight = G[node][neigb[0]]['weight'] + \
+        G[node][neigb[1]]['weight']
+
+    G.remove_node(node)
+    G.add_edge(neigb[0], neigb[1], weight=new_weight)
+
+    return G
+

fil_finder/length.py

@@ -862,22 +862,3 @@ def is_tpoint(vallist):
             return True
 
     return False
-
-
-def merge_nodes(node, G):
-    '''
-    Combine a node into its neighbors.
-    '''
-
-    neigb = list(G[node])
-
-    if len(neigb) != 2:
-        return G
-
-    new_weight = G[node][neigb[0]]['weight'] + \
-        G[node][neigb[1]]['weight']
-
-    G.remove_node(node)
-    G.add_edge(neigb[0], neigb[1], weight=new_weight)
-
-    return G

fil_finder/pixel_ident.py

@@ -56,7 +56,7 @@ def test_radial_profile_output(theta):
     model, skeleton = generate_filament_model(width=10.0,
                                               amplitude=1.0, background=0.0)
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,
@@ -75,7 +75,7 @@ def test_radial_profile_cutoff(cutoff):
     model, skeleton = generate_filament_model(width=10.0,
                                               amplitude=1.0, background=0.0)
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,
@@ -92,7 +92,7 @@ def test_radial_profile_padding(padding, max_distance=20.0):
     model, skeleton = generate_filament_model(width=10.0,
                                               amplitude=1.0, background=0.0)
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,
@@ -116,7 +116,7 @@ def test_radial_profile_fail_pad(padding=30.0, max_distance=20.0):
     model, skeleton = generate_filament_model(width=10.0,
                                               amplitude=1.0, background=0.0)
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,
@@ -138,7 +138,7 @@ def test_radial_profile_autocut():
 
     # all_skeleton += np.roll(skeleton, -30, axis=0)
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,
@@ -174,7 +174,7 @@ def test_radial_profile_autocut_plateau():
                                                   amplitude=5.0,
                                                   background=0.0)[::-1]
 
-    dist_transform = nd.distance_transform_edt((~skeleton).astype(np.int))
+    dist_transform = nd.distance_transform_edt((~skeleton).astype(int))
 
     dist, radprof, weights, unbin_dist, unbin_radprof = \
         radial_profile(model, dist_transform, dist_transform,

fil_finder/tests/test_whole_simplefilament.py

@@ -52,20 +52,6 @@ def generate_filament_model(shape=100, pad_size=0,
 
     return fits.PrimaryHDU(filament, header), centers
 
-
-# def make_skeleton_mask(shape=(256, 256), nfils=1):
-#     '''
-#     Generate a skeleton mask defining positions of filaments.
-#     '''
-
-
-#     mask = np.zeros(shape, dtype=bool)
-
-#     for _ in range(nfils):
-
-#         # Sample random start and end points
-#         start = (np.random.randint(shape[0]))
-
 def create_image_header(pixel_scale, beamfwhm, imshape,
                         restfreq, bunit):
     '''

fil_finder/tests/test_widths.py

@@ -14,7 +14,7 @@ github_project = e-koch/FilFinder
 [options]
 zip_safe = False
 packages = find:
-python_requires = >=3.7
+python_requires = >=3.9
 setup_requires = setuptools_scm
 install_requires =
     astropy

fil_finder/tests/testing_utils.py

@@ -1,6 +1,6 @@
 [tox]
 envlist =
-    py{36,37,38,39}-test{,-all,-dev,-cov}
+    py{39,310,311}-test{,-all,-dev,-cov}
     build_docs
     codestyle
 requires =
@@ -28,10 +28,11 @@ deps =
 extras =
     test
     all: all
+    dev: dev
 commands =
     pip freeze
-    !cov: pytest --open-files --pyargs fil_finder {toxinidir}/docs {posargs}
-    cov: pytest --open-files --pyargs fil_finder {toxinidir}/docs --cov fil_finder --cov-config={toxinidir}/setup.cfg {posargs}
+    !cov: pytest --pyargs fil_finder {toxinidir}/docs {posargs}
+    cov: pytest --pyargs fil_finder {toxinidir}/docs --cov fil_finder --cov-config={toxinidir}/setup.cfg {posargs}
     cov: coverage xml -o {toxinidir}/coverage.xml
 
 [testenv:build_docs]