Merge 02ca0e9 into 97b2118

radio_beam/beam.py

@@ -11,7 +11,7 @@
 from astropy.convolution import Kernel2D
 from astropy.convolution.kernels import _round_up_to_odd_integer
 
-from .utils import deconvolve, convolve, RadioBeamDeprecationWarning
+from .utils import deconvolve, deconvolve_opt, convolve, RadioBeamDeprecationWarning
 
 # Conversion between a twod Gaussian FWHM**2 and effective area
 FWHM_TO_AREA = 2*np.pi/(8*np.log(2))
@@ -383,10 +383,17 @@ def deconvolve(self, other, failure_returns_pointlike=False):
         """
 
         new_major, new_minor, new_pa = \
-            deconvolve(self, other,
-                       failure_returns_pointlike=failure_returns_pointlike)
+            deconvolve_opt(self.to_header_keywords(), other.to_header_keywords(),
+                           failure_returns_pointlike=True)
+
+        # Keep the units from before
+        new_major = (new_major * u.deg).to(self.major.unit)
+        new_minor = (new_minor * u.deg).to(self.minor.unit)
+        new_pa = (new_pa * u.rad).to(self.pa.unit)
+
         return Beam(major=new_major, minor=new_minor, pa=new_pa)
 
+
     def __eq__(self, other):
 
         # Catch floating point issues

radio_beam/commonbeam.py

@@ -10,7 +10,7 @@
     HAS_SCIPY = False
 
 from .beam import Beam
-from .utils import BeamError, transform_ellipse
+from .utils import BeamError, transform_ellipse, deconvolve_opt
 
 __all__ = ['commonbeam', 'common_2beams', 'getMinVolEllipse',
            'common_manybeams_mve']
@@ -336,13 +336,34 @@ def fits_in_largest(beams, large_beam=None):
     if large_beam is None:
         large_beam = beams.largest_beam()
 
-    for beam in beams:
-        if large_beam == beam:
+    large_hdr_keywords = large_beam.to_header_keywords()
+
+    majors = beams.major.to(u.deg).value
+    minors = beams.minor.to(u.deg).value
+    pas = beams.pa.to(u.deg).value
+
+    for major, minor, pa in zip(majors, minors, pas):
+
+        equal = abs(large_hdr_keywords['BMAJ'] - major) < 1e-12
+        equal = equal and (abs(large_hdr_keywords['BMIN'] - minor) < 1e-12)
+
+        # Check if the beam is circular
+        iscircular = (major - minor) / minor < 1e-6
+
+        # position angle only matters if the beam is asymmetric
+        if not iscircular:
+            equal = equal and (abs(((large_hdr_keywords['BPA'] % np.pi) - (pa % np.pi))) < 1e-12)
+
+        if equal:
             continue
-        deconv_beam = large_beam.deconvolve(beam,
-                                            failure_returns_pointlike=True)
 
-        if not deconv_beam.isfinite:
+        out = deconvolve_opt(large_hdr_keywords,
+                             {'BMAJ': major,
+                              'BMIN': minor,
+                              'BPA': pa},
+                             failure_returns_pointlike=True)
+
+        if np.any([ax == 0. for ax in out[:2]]):
             return False
 
     return True

radio_beam/tests/test_beam.py

@@ -18,7 +18,7 @@
 except ImportError:
     HAS_CASA = False
 
-from ..utils import RadioBeamDeprecationWarning
+from ..utils import RadioBeamDeprecationWarning, BeamError
 
 
 data_dir = os.path.join(os.path.dirname(__file__), 'data')
@@ -280,7 +280,7 @@ def test_deconv_pointlike(major, minor, pa, return_pointlike):
     else:
         try:
             beam1.deconvolve(beam1, failure_returns_pointlike=False)
-        except ValueError:
+        except BeamError:
             pass
 
 

radio_beam/utils.py

@@ -1,8 +1,11 @@
 
+import math
 import numpy as np
 import astropy.units as u
 
 
+DEG2RAD = math.pi / 180.
+
 class BeamError(Exception):
     """docstring for BeamError"""
     pass
@@ -16,6 +19,99 @@ class RadioBeamDeprecationWarning(Warning):
     pass
 
 
+def deconvolve_opt(beamprops1, beamprops2, failure_returns_pointlike=False):
+    """
+    An optimized, non-Quantity version of beam deconvolution.
+
+    Because no unit conversions are handled, the inputs MUST be in degrees for the major, minor,
+    and position angle.
+
+    Parameters
+    ----------
+    beamprops1: dict
+        Dictionary with keys 'BMAJ', 'BMIN', and 'BPA' for the beam to deconvolve from. Can be
+        produced with `~radio_beam.Beam.to_fits_keywords`.
+    beamprops2: dict
+        Same as `beamprops1` for the second beam.
+    failure_returns_pointlike : bool, optional
+        Return a point beam (zero area) when deconvolution fails. If `False`,
+        this will instead raise a `~radio_beam.utils.BeamError` when deconvolution fails.
+
+
+    Returns
+    -------
+    new_major : float
+        Deconvolved major FWHM.
+    new_minor : float
+        Deconvolved minor FWHM.
+    new_pa : float
+        Deconvolved position angle.
+
+    """
+
+    maj1 = beamprops1['BMAJ']
+    min1 = beamprops1['BMIN']
+    pa1 = beamprops1['BPA'] * DEG2RAD
+
+    maj2 = beamprops2['BMAJ']
+    min2 = beamprops2['BMIN']
+    pa2 = beamprops2['BPA'] * DEG2RAD
+
+    alpha = ((maj1 * math.cos(pa1))**2 +
+             (min1 * math.sin(pa1))**2 -
+             (maj2 * math.cos(pa2))**2 -
+             (min2 * math.sin(pa2))**2)
+
+    beta = ((maj1 * math.sin(pa1))**2 +
+            (min1 * math.cos(pa1))**2 -
+            (maj2 * math.sin(pa2))**2 -
+            (min2 * math.cos(pa2))**2)
+
+    gamma = 2 * ((min1**2 - maj1**2) * math.sin(pa1) * math.cos(pa1) -
+                 (min2**2 - maj2**2) * math.sin(pa2) * math.cos(pa2))
+
+    s = alpha + beta
+    t = math.sqrt((alpha - beta)**2 + gamma**2)
+
+    # Deal with floating point issues
+    # This matches the arcsec**2 check for deconvolve below
+    # Difference is we keep things in deg^2 here
+    atol_t = np.finfo(np.float64).eps / 3600.**2
+
+    # To deconvolve, the beam must satisfy:
+    # alpha < 0
+    alpha_cond = alpha + np.finfo(np.float64).eps < 0
+
+    # beta < 0
+    beta_cond = beta + np.finfo(np.float64).eps < 0
+    # s < t
+    st_cond = s < t + atol_t
+
+    if alpha_cond or beta_cond or st_cond:
+        if failure_returns_pointlike:
+            return 0., 0., 0.
+        else:
+            raise BeamError("Beam could not be deconvolved")
+    else:
+        new_major = math.sqrt(0.5 * (s + t))
+        new_minor = math.sqrt(0.5 * (s - t))
+
+        # absolute tolerance needs to be <<1 microarcsec
+        atol = 1e-7 / 3600.
+        if (math.sqrt(abs(gamma) + abs(alpha - beta))) < atol:
+            new_pa = 0.0
+        else:
+            new_pa = 0.5 * math.atan2(-1. * gamma, alpha - beta)
+
+    # In the limiting case, the axes can be zero to within precision
+    # Add the precision level onto each axis so a deconvolvable beam
+    # is always has beam.isfinite == True
+    new_major += np.finfo(np.float64).eps
+    new_minor += np.finfo(np.float64).eps
+
+    return new_major, new_minor, new_pa
+
+
 def deconvolve(beam, other, failure_returns_pointlike=False):
     """
     Deconvolve a beam from another
@@ -78,6 +174,7 @@ def deconvolve(beam, other, failure_returns_pointlike=False):
 
     # Deal with floating point issues
     atol_t = np.finfo(t.dtype).eps * t.unit
+    print(f"deconvolve atol_t {atol_t}")
 
     # To deconvolve, the beam must satisfy:
     # alpha < 0
@@ -87,6 +184,9 @@ def deconvolve(beam, other, failure_returns_pointlike=False):
     # s < t
     st_cond = s < t + atol_t
 
+    print(f"deconvolve alpha {alpha} beta {beta} s {s} t {t}")
+    print(f"deconvolve alpha {alpha_cond} beta {beta_cond} st {st_cond}")
+
     if alpha_cond or beta_cond or st_cond:
         if failure_returns_pointlike:
             return 0 * maj1.unit, 0 * min1.unit, 0 * pa1.unit