Fix kaiser bessel filter and taper (#48)

HISTORY.rst

@@ -4,6 +4,7 @@ History
 
 0.1.3 (2018-03-28)
 ------------------
+* Fix Kaiser Bessel filter and taper (:pr:`48`)
 * Stokes/Correlation conversion (:pr:`41`)
 * Fix gridding examples (:pr:`43`)
 * Add simple dask gridder example (:pr:`42`)

africanus/filters/conv_filters.py

@@ -8,7 +8,7 @@
 
 import numpy as np
 
-from .kaiser_bessel_filter import (kaiser_bessel,
+from .kaiser_bessel_filter import (kaiser_bessel_with_sinc,
                                    estimate_kaiser_bessel_beta)
 
 ConvolutionFilter = collections.namedtuple("ConvolutionFilter",
@@ -47,6 +47,10 @@
 """
 
 
+class AsymmetricKernel(Exception):
+    pass
+
+
 def convolution_filter(half_support, oversampling_factor,
                        filter_type, **kwargs):
     r"""
@@ -68,6 +72,9 @@ def convolution_filter(half_support, oversampling_factor,
     beta : float, optional
         Beta shape parameter for
         `Kaiser Bessel <kaiser-bessel-filter_>`_ filters.
+    normalise : {True, False}
+        Normalise the filter by the it's volume.
+        Defaults to ``True``.
 
     Returns
     -------
@@ -78,6 +85,8 @@ def convolution_filter(half_support, oversampling_factor,
     full_sup = full_sup_wo_padding + 2  # + padding
     no_taps = full_sup + (full_sup - 1) * (oversampling_factor - 1)
 
+    normalise = kwargs.pop("normalise", True)
+
     taps = np.arange(no_taps) / oversampling_factor - full_sup // 2
 
     if filter_type == 'sinc':
@@ -89,16 +98,19 @@ def convolution_filter(half_support, oversampling_factor,
         except KeyError:
             beta = estimate_kaiser_bessel_beta(full_sup)
 
-        filter_taps = kaiser_bessel(taps, full_sup, beta)
-
-        # Normalise by integrating
-        filter_taps /= np.trapz(filter_taps, taps)
+        # Compute Kaiser Bessel and multiply in the sinc
+        filter_taps = kaiser_bessel_with_sinc(taps, full_sup,
+                                              oversampling_factor, beta,
+                                              normalise=normalise)
     else:
-        raise ValueError("Expected one of {'kaiser-bessel'}")
+        raise ValueError("Expected one of {'kaiser-bessel', 'sinc'}")
 
     # Expand filter taps to 2D
     filter_taps = np.outer(filter_taps, filter_taps)
 
+    if not np.all(filter_taps == filter_taps.T):
+        raise AsymmetricKernel("Kernel is asymmetric")
+
     return ConvolutionFilter(half_support, oversampling_factor,
                              full_sup_wo_padding, full_sup,
                              no_taps, filter_taps)

africanus/filters/filter_tapers.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 
-from .kaiser_bessel_filter import (kaiser_bessel_fourier,
+from .kaiser_bessel_filter import (kaiser_bessel_with_sinc,
                                    estimate_kaiser_bessel_beta)
 
 
@@ -28,23 +28,51 @@ def taper(filter_type, ny, nx, conv_filter, **kwargs):
     :class:`numpy.ndarray`
         Taper of shape :code:`(ny, nx)`
     """
+    cf = conv_filter
+
     if filter_type == "sinc":
-        raise NotImplementedError("Please implement sinc support")
+        return np.ones((ny, nx))
     elif filter_type == "kaiser-bessel":
         try:
             beta = kwargs.pop('beta')
         except KeyError:
-            beta = estimate_kaiser_bessel_beta(conv_filter.full_sup)
+            beta = estimate_kaiser_bessel_beta(cf.full_sup)
+
+        # What would Andre Offringa do?
+        # He would compute the numeric solution
+        taps = np.arange(cf.no_taps) / cf.oversample - cf.full_sup // 2
+        kb = kaiser_bessel_with_sinc(taps, cf.full_sup, cf.oversample, beta)
+        kbshift = np.fft.fftshift(kb)
+
+        width = nx * cf.oversample
+        height = ny * cf.oversample
+
+        # Put the first and last halves of the shifted Kaiser Bessel
+        # at each end of the output buffer, then FFT
+        buf = np.zeros(width, dtype=kb.dtype)
+        buf[:kbshift.size // 2] = kbshift[:kbshift.size // 2]
+        buf[-kbshift.size // 2:] = kbshift[-kbshift.size // 2:]
+        x = np.fft.ifft(buf).real
+
+        buf = np.zeros(height, dtype=kb.dtype)
+        buf[:kbshift.size // 2] = kbshift[:kbshift.size // 2]
+        buf[-kbshift.size // 2:] = kbshift[-kbshift.size // 2:]
+        y = np.fft.ifft(buf).real
+
+        # First quarter of the taper
+        quarter = y[:ny // 2, None] * x[None, :nx // 2]
 
-        py = np.arange(ny) / ny - 0.5
-        y = kaiser_bessel_fourier(py, conv_filter.full_sup, beta)
-        y *= np.sinc(y / conv_filter.oversample)
+        # Create the taper by copying
+        # the quarter into the appropriate bits
+        taper = np.empty((ny, nx), dtype=kb.dtype)
+        taper[:ny // 2, :nx // 2] = quarter[::-1, ::-1]
+        taper[ny // 2:, :nx // 2] = quarter[:, ::-1]
+        taper[:ny // 2, nx // 2:] = quarter[::-1, :]
+        taper[ny // 2:, nx // 2:] = quarter[:, :]
 
-        px = np.arange(nx) / nx - 0.5
-        x = kaiser_bessel_fourier(px, conv_filter.full_sup, beta)
-        x *= np.sinc(x / conv_filter.oversample)
+        # Normalise by oversampling factor
+        taper *= cf.oversample**2
 
-        # Produce 2D taper
-        return y[:, None] * x[None, :]
+        return taper
     else:
         raise ValueError("Invalid filter_type '%s'" % filter_type)

africanus/filters/kaiser_bessel_filter.py

@@ -8,84 +8,125 @@
 import numpy as np
 
 
-def estimate_kaiser_bessel_beta(full_support):
+def estimate_kaiser_bessel_beta(W):
     r"""
-    Estimate the kaiser bessel beta using he following heuristic,
-    with :math:`W` denoting ``full_support``:
+    Estimate the kaiser bessel beta using the following heuristic:
 
     .. math::
 
-        \beta = 1.2 \pi \sqrt{0.25 \text{ W }^2 - 1.0 }
+        \beta = 2.34 \times W
+
+    Derived from `Nonuniform fast Fourier transforms
+    using min-max interpolation
+    <nufft-min-max-ref_>`_.
+
+    .. _nufft-min-max-ref: https://ieeexplore.ieee.org/document/1166689/
 
     Parameters
     ----------
-    full_support : int
-        Full support of the filter
+    W : int
+        Width of the filter
 
     Returns
     -------
     float
-        kaiser Bessel beta shape parameter
+        Kaiser Bessel beta shape parameter
     """
-
-    # NOTE(bmerry)
-    # Puts the first null of the taper function
-    # at the edge of the image
-    beta = np.pi * np.sqrt(0.25 * full_support**2 - 1.0)
-    # Move the null outside the image,
-    # to avoid numerical instabilities.
-    # This will cause a small amount of aliasing at the edges,
-    # which ideally should be handled by clipping the image.
-    return beta * 1.2
+    return 2.34*W
 
 
-def kaiser_bessel(taps, full_support, beta):
+def kaiser_bessel(u, W, beta):
     r"""
-    Compute a 1D Kaiser Bessel filter.
+    Compute a 1D Kaiser Bessel filter as defined
+    in `Selection of a Convolution Function
+    for Fourier Inversion Using Gridding
+    <kbref_>`_.
+
+    .. _kbref: https://ieeexplore.ieee.org/document/97598/
 
     Parameters
     ----------
-    taps : :class:`numpy.ndarray`
-        Filter position
-    full_support : int
-        Full support of the filter
+    u : :class:`numpy.ndarray`
+        Filter positions
+    W : int
+        Width of the filter
     beta : float, optional
         Kaiser Bessel shape parameter
 
     Returns
     -------
     :class:`numpy.ndarray`
-        Kaiser Bessel filter of shape :code:`(taps,)`
+        Kaiser Bessel filter with the same shape as `u`
     """
 
     # Sanity check
-    M = full_support
-    hM = M // 2
-    assert np.all(-hM <= taps) & np.all(taps <= hM)
+    hW = W // 2
+    assert np.all(-hW <= u) & np.all(u <= hW)
 
-    param = 1 - (2 * taps / M)**2
+    param = 1 - (2 * u / W)**2
     param[param < 0] = 0  # Zero negative numbers
     return np.i0(beta * np.sqrt(param)) / np.i0(beta)
 
 
-def kaiser_bessel_fourier(positions, full_support, beta):
+def kaiser_bessel_with_sinc(u, W, oversample, beta, normalise=True):
+    """
+    Produces a filter composed of Kaiser Bessel multiplied by a sinc.
+
+    Accounts for the oversampling factor, as well as normalising the filter.
+
+    Parameters
+    ----------
+    u : :class:`numpy.ndarray`
+        Filter positions
+    W : int
+        Width of the filter
+    oversample : int
+        Oversampling factor
+    beta : float
+        Kaiser Bessel shape parameter
+    normalise : optional, {True, False}
+        True if the filter should be normalised
+
+    Returns
+    -------
+    :class:`numpy.ndarray`
+        Filter with the same shape as `u`
+    """
+    kb = kaiser_bessel(u, W, beta)
+    kb *= oversample
+    kb *= np.sinc(u / oversample)
+
+    if normalise:
+        kb /= np.trapz(kb, u)
+
+    return kb
+
+
+def kaiser_bessel_fourier(x, W, beta):
     r"""
     Computes the Fourier Transform of a 1D Kaiser Bessel filter.
+    as defined in `Selection of a Convolution Function
+    for Fourier Inversion Using Gridding
+    <kbref_>`_.
+
+    .. _kbref: https://ieeexplore.ieee.org/document/97598/
+
 
     Parameters
     ----------
-    positions : :class:`numpy.ndarray`
+    x : :class:`numpy.ndarray`
         Filter positions
-    full_support : int
-        Full support of the associated convolution filter.
+    W : int
+        Width of the filter.
     beta : float
         Kaiser bessel shape parameter
 
     Returns
     -------
     :class:`numpy.ndarray`
-        Array of shape :code:`(positions,)
+        Fourier Transform of the Kaiser Bessel,
+        with the same shape as `x`.
     """
-    alpha = beta / np.pi
-    inner = np.lib.scimath.sqrt((full_support * positions)**2 - alpha * alpha)
-    return full_support * np.sinc(inner).real / np.i0(beta)
+    term = (np.pi*W*x)**2 - beta**2
+    val = np.lib.scimath.sqrt(term).real
+    return np.sin(val)/val

africanus/filters/plot_filter.py

@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+
+import argparse
+import logging
+
+from ..filters import convolution_filter
+from ..util.cmdline import parse_python_assigns
+from ..util.requirements import requires_optional
+
+import numpy as np
+
+try:
+    import matplotlib.pyplot as plt
+    from matplotlib import cm
+    from matplotlib.ticker import LinearLocator, FormatStrFormatter
+    from mpl_toolkits.mplot3d import Axes3D  # noqa
+except ImportError:
+    pass
+
+
+def create_parser():
+    p = argparse.ArgumentParser()
+    p.add_argument("filter", choices=['kaiser-bessel', 'sinc'],
+                   default='kaiser-bessel')
+    p.add_argument("-hs", "--half-support", default=3, type=int)
+    p.add_argument("-os", "--oversample", default=63, type=int)
+    p.add_argument("-n", "--normalise", dest="normalise", action="store_true",
+                   help="Normalise filter by it's volume")
+    p.add_argument("--no-normalise", dest="normalise", action="store_false",
+                   help="Do not normalise the filter by it's volume")
+    p.add_argument("-k", "--kwargs", default="", type=parse_python_assigns,
+                   help="Extra keywords arguments used to create the filter. "
+                        "For example 'beta=2.3' to specify a beta shape "
+                        "parameter for the Kaiser Bessel")
+
+    return p
+
+
+@requires_optional('matplotlib.pyplot', 'mpl_toolkits.mplot3d')
+def _plot_filter(data):
+    X, Y = np.mgrid[-1:1:1j*data.shape[0], -1:1:1j*data.shape[1]]
+
+    fig = plt.figure()
+    ax = fig.gca(projection='3d')
+
+    # Plot the surface.
+    surf = ax.plot_surface(X, Y, data, cmap=cm.coolwarm,
+                           linewidth=0, antialiased=False)
+
+    zmax = data.max()
+
+    ax.plot([0, 0], [0, 0], zs=[0, 1.2*zmax], color='black')
+    ax.plot([0, 0], [-1, 1], zs=[zmax, zmax], color='black')
+    ax.plot([-1, 1], [0, 0], zs=[zmax, zmax], color='black')
+
+    # Customize the z axis.
+    ax.zaxis.set_major_locator(LinearLocator(10))
+    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
+
+    # Add a color bar which maps values to colors.
+    fig.colorbar(surf, shrink=0.5, aspect=5)
+
+    plt.show()
+
+
+def main():
+    logging.basicConfig(level=logging.INFO, format="%(message)s")
+
+    args = create_parser().parse_args()
+
+    logging.info("Creating %s filter with half support of %d "
+                 "and %d oversampling" %
+                 (args.filter, args.half_support, args.oversample))
+
+    if len(args.kwargs) > 0:
+        logging.info("Extra keywords %s" % args.kwargs)
+
+    conv_filter = convolution_filter(args.half_support,
+                                     args.oversample,
+                                     args.filter,
+                                     normalise=args.normalise,
+                                     **args.kwargs)
+
+    _plot_filter(np.abs(conv_filter.filter_taps))