Merge pull request #299

Test image processors module

pyobs/images/processor.py

@@ -20,11 +20,9 @@ async def __call__(self, image: Image) -> Image:
         Returns:
             Processed image.
         """
-        ...
 
     async def reset(self) -> None:
         """Resets state of image processor"""
-        pass
 
 
 __all__ = ["ImageProcessor"]

pyobs/images/processors/_daobackgroundremover.py

@@ -0,0 +1,37 @@
+from typing import Tuple
+
+import numpy as np
+from astropy.stats import SigmaClip
+
+from pyobs.images import Image
+
+
+class _DaoBackgroundRemover:
+    def __init__(self, sigma: float, box_size: Tuple[int, int], filter_size: Tuple[int, int]):
+        from photutils.background import MedianBackground
+
+        self._sigma_clip = SigmaClip(sigma=sigma)
+        self._box_size = box_size
+        self._filter_size = filter_size
+
+        self._bkg_estimator = MedianBackground()
+
+    def __call__(self, image: Image) -> Image:
+        background = self._estimate_background(image)
+        return self._remove_background(image, background)
+
+    def _estimate_background(self, image: Image):
+        from photutils.background import Background2D
+
+        bkg = Background2D(
+            image.data, box_size=self._box_size, filter_size=self._filter_size, sigma_clip=self._sigma_clip,
+            bkg_estimator=self._bkg_estimator, mask=image.mask
+        )
+
+        return bkg.background
+
+    @staticmethod
+    def _remove_background(image: Image, background: np.ndarray) -> Image:
+        output_image = image.copy()
+        output_image.data = output_image.data - background
+        return output_image

pyobs/images/processors/detection/_pysep_stats_calculator.py

@@ -0,0 +1,99 @@
+import asyncio
+from copy import copy
+from functools import partial
+from typing import Optional
+
+import numpy as np
+
+from pyobs.images.processors.detection._source_catalog import _SourceCatalog
+
+
+class PySepStatsCalculator:
+    def __init__(self, catalog: _SourceCatalog, data: np.ndarray, mask: np.ndarray, gain: Optional[float]):
+        self._catalog = copy(catalog)
+        self._data = data
+        self._mask = mask
+        self._gain = gain
+
+    async def __call__(self, *args, **kwargs) -> _SourceCatalog:
+        self._calc_ellipticity()
+        self._calc_fwhm()
+        self._calc_kron_radius()
+
+        loop = asyncio.get_running_loop()
+        await loop.run_in_executor(None, partial(self._calc_flux))
+
+        self._calc_flux_radii()
+        self._calc_winpos()
+
+        return self._catalog
+
+    def _calc_ellipticity(self):
+        self._catalog.sources["ellipticity"] = 1.0 - (self._catalog.sources["b"] / self._catalog.sources["a"])
+
+    def _calc_fwhm(self):
+        fwhm = 2.0 * (np.log(2) * (self._catalog.sources["a"] ** 2.0 + self._catalog.sources["b"] ** 2.0)) ** 0.5
+        self._catalog.sources["fwhm"] = fwhm
+
+    def _calc_kron_radius(self):
+        import sep
+
+        kronrad, krflag = sep.kron_radius(
+            self._data,
+            self._catalog.sources["x"],
+            self._catalog.sources["y"],
+            self._catalog.sources["a"],
+            self._catalog.sources["b"],
+            self._catalog.sources["theta"],
+            6.0,
+        )
+        self._catalog.sources["flag"] |= krflag
+        self._catalog.sources["kronrad"] = kronrad
+
+    def _calc_flux(self):
+        import sep
+
+        flux, _, flag = sep.sum_ellipse(
+            self._data,
+            self._catalog.sources["x"],
+            self._catalog.sources["y"],
+            self._catalog.sources["a"],
+            self._catalog.sources["b"],
+            self._catalog.sources["theta"],
+            2.5 * self._catalog.sources["kronrad"],
+            subpix=5,
+            mask=self._mask,
+            gain=self._gain,
+        )
+
+        self._catalog.sources["flag"] |= flag
+        self._catalog.sources["flux"] = flux
+
+    def _calc_flux_radii(self):
+        import sep
+
+        flux_radii, flag = sep.flux_radius(
+            self._data,
+            self._catalog.sources["x"],
+            self._catalog.sources["y"],
+            6.0 * self._catalog.sources["a"],
+            [0.25, 0.5, 0.75],
+            normflux=self._catalog.sources["flux"],
+            subpix=5,
+        )
+
+        self._catalog.sources["flag"] |= flag
+        self._catalog.sources["fluxrad25"] = flux_radii[:, 0]
+        self._catalog.sources["fluxrad50"] = flux_radii[:, 1]
+        self._catalog.sources["fluxrad75"] = flux_radii[:, 2]
+
+    def _calc_winpos(self):
+        import sep
+
+        sig = 2.0 / 2.35 * self._catalog.sources["fluxrad50"]
+        xwin, ywin, flag = sep.winpos(self._data, self._catalog.sources["x"], self._catalog.sources["y"], sig)
+
+        self._catalog.sources["flag"] |= flag
+        self._catalog.sources["xwin"] = xwin
+        self._catalog.sources["ywin"] = ywin
+

pyobs/images/processors/detection/_source_catalog.py

@@ -0,0 +1,56 @@
+from copy import copy
+from typing import Optional, List
+
+import numpy as np
+import pandas as pd
+from astropy.coordinates import Angle
+from astropy.table import Table
+
+from pyobs.images import Image
+
+
+class _SourceCatalog:
+    def __init__(self, sources: pd.DataFrame):
+        self.sources = sources
+
+    @classmethod
+    def from_array(cls, sources: np.ndarray):
+        source_dataframe = pd.DataFrame(sources)
+        return cls(source_dataframe)
+
+    @classmethod
+    def from_table(cls, sources: Table):
+        sources.rename_column("xcentroid", "x")
+        sources.rename_column("ycentroid", "y")
+
+        source_dataframe = sources.to_pandas()
+        return cls(source_dataframe)
+
+    def filter_detection_flag(self):
+        if "flag" not in self.sources:
+            return
+
+        self.sources = self.sources[self.sources["flag"] < 8]
+
+    def wrap_rotation_angle_at_ninty_deg(self):
+        if "theta" not in self.sources:
+            return
+
+        self.sources["theta"] = np.arcsin(np.sin(self.sources["theta"]))
+
+    def rotation_angle_to_degree(self):
+        if "theta" not in self.sources:
+            return
+
+        self.sources["theta"] = np.degrees(self.sources["theta"])
+
+    def apply_fits_origin_convention(self):
+        self.sources["x"] += 1
+        self.sources["y"] += 1
+
+    def save_to_image(self, image: Image, keys: List[str]) -> Image:
+        cat = self.sources[keys]
+
+        output_image = copy(image)
+        output_image.catalog = Table.from_pandas(cat)
+        return output_image

pyobs/images/processors/detection/daophot.py

@@ -1,9 +1,15 @@
 import asyncio
-from typing import Tuple, Any
 import logging
+from typing import Tuple, Any
+
+import numpy as np
+from astropy.stats import sigma_clipped_stats
+from astropy.table import Table
 
-from .sourcedetection import SourceDetection
 from pyobs.images import Image
+from ._source_catalog import _SourceCatalog
+from .sourcedetection import SourceDetection
+from .._daobackgroundremover import _DaoBackgroundRemover
 
 log = logging.getLogger(__name__)
 
@@ -13,6 +19,8 @@ class DaophotSourceDetection(SourceDetection):
 
     __module__ = "pyobs.images.processors.detection"
 
+    _CATALOG_KEYS = ["x", "y", "flux", "peak"]
+
     def __init__(
         self,
         fwhm: float = 3.0,
@@ -26,7 +34,7 @@ def __init__(
 
         Args:
             fwhm: Full-width at half maximum for Gaussian kernel.
-            threshold: Threshold pixel value for detection.
+            threshold: Threshold pixel value for detection in standard deviations.
             bkg_sigma: Sigma for background kappa-sigma clipping.
             bkg_box_size: Box size for background estimation.
             bkg_filter_size: Filter size for background estimation.
@@ -36,9 +44,16 @@ def __init__(
         # store
         self.fwhm = fwhm
         self.threshold = threshold
-        self.bkg_sigma = bkg_sigma
-        self.bkg_box_size = bkg_box_size
-        self.bkg_filter_size = bkg_filter_size
+
+        self._background_remover = _DaoBackgroundRemover(bkg_sigma, bkg_box_size, bkg_filter_size)
+
+    async def _find_stars(self, data: np.ndarray, std: int) -> Table:
+        from photutils import DAOStarFinder
+
+        daofind = DAOStarFinder(fwhm=self.fwhm, threshold=self.threshold * std)
+
+        loop = asyncio.get_running_loop()
+        return await loop.run_in_executor(None, daofind, data)
 
     async def __call__(self, image: Image) -> Image:
         """Find stars in given image and append catalog.
@@ -49,51 +64,23 @@ async def __call__(self, image: Image) -> Image:
         Returns:
             Image with attached catalog.
         """
-        from astropy.stats import SigmaClip, sigma_clipped_stats
-        from photutils import Background2D, MedianBackground, DAOStarFinder
 
-        # get data
         if image.data is None:
             log.warning("No data found in image.")
             return image
-        data = image.data.astype(float).copy()
-
-        # estimate background
-        sigma_clip = SigmaClip(sigma=self.bkg_sigma)
-        bkg_estimator = MedianBackground()
-        bkg = Background2D(
-            data,
-            self.bkg_box_size,
-            filter_size=self.bkg_filter_size,
-            sigma_clip=sigma_clip,
-            bkg_estimator=bkg_estimator,
-            mask=image.mask,
-        )
-        data -= bkg.background
-
-        # do statistics
-        mean, median, std = sigma_clipped_stats(data, sigma=3.0)
-
-        # find stars
-        daofind = DAOStarFinder(fwhm=self.fwhm, threshold=self.threshold * std)
-        loop = asyncio.get_running_loop()
-        sources = await loop.run_in_executor(None, daofind, data - median)
 
-        # rename columns
-        sources.rename_column("xcentroid", "x")
-        sources.rename_column("ycentroid", "y")
+        background_corrected_image = self._background_remover(image)
+        background_corrected_data = background_corrected_image.data.astype(float)
 
-        # match fits conventions
-        sources["x"] += 1
-        sources["y"] += 1
+        _, median, std = sigma_clipped_stats(background_corrected_data, sigma=3.0)
 
-        # pick columns for catalog
-        cat = sources["x", "y", "flux", "peak"]
+        median_corrected_data = background_corrected_data - median
+        sources = await self._find_stars(median_corrected_data, std)
 
-        # copy image, set catalog and return it
-        img = image.copy()
-        img.catalog = cat
-        return img
+        sources_catalog = _SourceCatalog.from_table(sources)
+        sources_catalog.apply_fits_origin_convention()
+        output_image = sources_catalog.save_to_image(image, self._CATALOG_KEYS)
+        return output_image
 
 
 __all__ = ["DaophotSourceDetection"]