def filter_image_sep2d(image, fh, fv, impl='numpy', padding=None):
    """Filter an image with a separable filter.

    Parameters
    ----------
    image : 2D array-like
        The image to be filtered. It must have a real (vs. complex) dtype.
    fh, fv : 1D array-like
        Horizontal (axis 0) and vertical (axis 1) filters. Their sizes
        can be at most the image sizes in the respective axes.
    impl : {'numpy', 'pyfftw'}, optional
        FFT backend to use. The ``pyfftw`` backend requires the
        ``pyfftw`` package to be installed. It is usually significantly
        faster than the NumPy backend.
    padding : positive int, optional
        Amount of zeros added to the left and right of the image in all
        axes before FFT. This helps avoiding wraparound artifacts due to
        large boundary values.
        For ``None``, the padding is computed as ::

            padding = min(max(len(fh), len(fv)) - 1, 64)

        A padding of ``len(filt) - 1`` ensures that errors in FFT-based
        convolutions are small. At the same time, the padding should not
        be excessive to retain efficiency.

    Returns
    -------
    filtered : 2D `numpy.ndarray`
        The image filtered horizontally by ``fh`` and vertically by ``fv``.
        It has the same shape as ``image``, and its dtype is
        ``np.result_type(image, fh, fv)``.
    """
    # TODO: generalize for nD
    impl, impl_in = str(impl).lower(), impl
    if impl not in ('numpy', 'pyfftw'):
        raise ValueError('`impl` {!r} not understood'
                         ''.format(impl_in))

    image = np.asarray(image)
    if image.ndim != 2:
        raise ValueError('`image` must be 2-dimensional, got image with '
                         'ndim={}'.format(image.ndim))
    if image.size == 0:
        raise ValueError('`image` cannot have size 0')
    if not np.issubsctype(image.dtype, np.floating):
        image = image.astype(float)

    fh = np.asarray(fh).astype(image.dtype)
    if fh.ndim != 1:
        raise ValueError('`fh` must be one-dimensional')
    elif fh.size == 0:
        raise ValueError('`fh` cannot have size 0')
    elif fh.size > image.shape[0]:
        raise ValueError('`fh` can be at most `image.shape[0]`, got '
                         '{} > {}'.format(fh.size, image.shape[0]))

    fv = np.asarray(fv).astype(image.dtype)
    if fv.ndim != 1:
        raise ValueError('`fv` must be one-dimensional')
    elif fv.size == 0:
        raise ValueError('`fv` cannot have size 0')
    elif fv.size > image.shape[0]:
        raise ValueError('`fv` can be at most `image.shape[1]`, got '
                         '{} > {}'.format(fv.size, image.shape[1]))

    # Pad image with zeros
    if padding is None:
        padding = min(max(len(fh), len(fv)) - 1, 64)

    if padding != 0:
        image_padded = np.pad(image, padding, mode='constant')
    else:
        image_padded = image.copy() if impl == 'pyfftw' else image

    # Prepare filters for the convolution
    def prepare_for_fft(filt, n_new):
        """Return padded and shifted filter ready for FFT.

        The filter is padded with zeros to the new size, and then shifted
        such that such that the middle element of old filter, i.e., the
        one at index ``(len(filt) - 1) // 2`` ends up at index 0.
        """
        mid = (len(filt) - 1) // 2
        padded = np.zeros(n_new, dtype=filt.dtype)
        padded[:len(filt) - mid] = filt[mid:]
        padded[len(padded) - mid:] = filt[:mid]
        return padded

    fh = prepare_for_fft(fh, image_padded.shape[0])
    fv = prepare_for_fft(fv, image_padded.shape[1])

    # Perform the multiplication in Fourier space and apply inverse FFT
    if impl == 'numpy':
        image_ft = np.fft.rfftn(image_padded)
        fh_ft = np.fft.fft(fh)
        fv_ft = np.fft.rfft(fv)

        image_ft *= fh_ft[:, None]
        image_ft *= fv_ft[None, :]
        # Important to specify the shape since `irfftn` cannot know the
        # original shape
        conv = np.fft.irfftn(image_ft, s=image_padded.shape)
        if conv.dtype != image.dtype:
            conv = conv.astype(image.dtype)

    elif impl == 'pyfftw':
        if not PYFFTW_AVAILABLE:
            raise ValueError(
                '`pyfftw` package is not available; you need to install it '
                'to use the pyfftw backend')

        import pyfftw
        import multiprocessing

        # Generate output arrays, for half-complex transform of image and
        # vertical filter, and full FT of the horizontal filter
        out_img_shape = (image_padded.shape[0], image_padded.shape[1] // 2 + 1)
        out_img_dtype = np.result_type(image_padded, 1j)
        out_img = np.empty(out_img_shape, out_img_dtype)

        out_fh_shape = out_img_shape[0]
        out_fh_dtype = np.result_type(fh, 1j)
        fh_c = fh.astype(out_fh_dtype)  # need to make this a C2C trafo
        out_fh = np.empty(out_fh_shape, out_fh_dtype)

        out_fv_shape = out_img_shape[1]
        out_fv_dtype = np.result_type(fv, 1j)
        out_fv = np.empty(out_fv_shape, out_fv_dtype)

        # Perform the forward transforms of image and filters. We use
        # the `FFTW_ESTIMATE` flag to not allow the planner to destroy
        # the input.
        plan = pyfftw.FFTW(image_padded, out_img, axes=(0, 1),
                           direction='FFTW_FORWARD',
                           flags=['FFTW_ESTIMATE'],
                           threads=multiprocessing.cpu_count())
        plan(image_padded, out_img)

        plan = pyfftw.FFTW(fh_c, out_fh, axes=(0,),
                           direction='FFTW_FORWARD',
                           flags=['FFTW_ESTIMATE'],
                           threads=multiprocessing.cpu_count())
        plan(fh_c, out_fh)

        plan = pyfftw.FFTW(fv, out_fv, axes=(0,),
                           direction='FFTW_FORWARD',
                           flags=['FFTW_ESTIMATE'],
                           threads=multiprocessing.cpu_count())
        plan(fv, out_fv)

        # Fourier space multiplication
        out_img *= out_fh[:, None]
        out_img *= out_fv[None, :]

        # Inverse trafo
        conv = image_padded  # Overwrite
        plan = pyfftw.FFTW(out_img.copy(), conv, axes=(0, 1),
                           direction='FFTW_BACKWARD',
                           flags=['FFTW_ESTIMATE'],
                           threads=multiprocessing.cpu_count())
        plan(out_img, conv)

    else:
        raise ValueError('unsupported `impl` {!r}'.format(impl_in))

    if padding:
        return conv[padding:-padding, padding:-padding]
    else:
        return conv
def filter_image_sep2d ( image , fh , fv , impl = 'numpy' , padding = None ) : impl , impl_in = str ( impl ) . lower ( ) , impl if impl not in ( 'numpy' , 'pyfftw' ) : raise ValueError ( '`impl` {!r} not understood' '' . format ( impl_in ) ) image = np . asarray ( image ) if image . ndim != 2 : raise ValueError ( '`image` must be 2-dimensional, got image with ' 'ndim={}' . format ( image . ndim ) ) if image . size == 0 : raise ValueError ( '`image` cannot have size 0' ) if not np . issubsctype ( image . dtype , np . floating ) : image = image . astype ( float ) fh = np . asarray ( fh ) . astype ( image . dtype ) if fh . ndim != 1 : raise ValueError ( '`fh` must be one-dimensional' ) elif fh . size == 0 : raise ValueError ( '`fh` cannot have size 0' ) elif fh . size > image . shape [ 0 ] : raise ValueError ( '`fh` can be at most `image.shape[0]`, got ' '{} > {}' . format ( fh . size , image . shape [ 0 ] ) ) fv = np . asarray ( fv ) . astype ( image . dtype ) if fv . ndim != 1 : raise ValueError ( '`fv` must be one-dimensional' ) elif fv . size == 0 : raise ValueError ( '`fv` cannot have size 0' ) elif fv . size > image . shape [ 0 ] : raise ValueError ( '`fv` can be at most `image.shape[1]`, got ' '{} > {}' . format ( fv . size , image . shape [ 1 ] ) ) if padding is None : padding = min ( max ( len ( fh ) , len ( fv ) ) - 1 , 64 ) if padding != 0 : image_padded = np . pad ( image , padding , mode = 'constant' ) else : image_padded = image . copy ( ) if impl == 'pyfftw' else image def prepare_for_fft ( filt , n_new ) : mid = ( len ( filt ) - 1 ) // 2 padded = np . zeros ( n_new , dtype = filt . dtype ) padded [ : len ( filt ) - mid ] = filt [ mid : ] padded [ len ( padded ) - mid : ] = filt [ : mid ] return padded fh = prepare_for_fft ( fh , image_padded . shape [ 0 ] ) fv = prepare_for_fft ( fv , image_padded . shape [ 1 ] ) if impl == 'numpy' : image_ft = np . fft . rfftn ( image_padded ) fh_ft = np . fft . fft ( fh ) fv_ft = np . fft . rfft ( fv ) image_ft *= fh_ft [ : , None ] image_ft *= fv_ft [ None , : ] conv = np . fft . irfftn ( image_ft , s = image_padded . shape ) if conv . dtype != image . dtype : conv = conv . astype ( image . dtype ) elif impl == 'pyfftw' : if not PYFFTW_AVAILABLE : raise ValueError ( '`pyfftw` package is not available; you need to install it ' 'to use the pyfftw backend' ) import pyfftw import multiprocessing out_img_shape = ( image_padded . shape [ 0 ] , image_padded . shape [ 1 ] // 2 + 1 ) out_img_dtype = np . result_type ( image_padded , 1j ) out_img = np . empty ( out_img_shape , out_img_dtype ) out_fh_shape = out_img_shape [ 0 ] out_fh_dtype = np . result_type ( fh , 1j ) fh_c = fh . astype ( out_fh_dtype ) out_fh = np . empty ( out_fh_shape , out_fh_dtype ) out_fv_shape = out_img_shape [ 1 ] out_fv_dtype = np . result_type ( fv , 1j ) out_fv = np . empty ( out_fv_shape , out_fv_dtype ) plan = pyfftw . FFTW ( image_padded , out_img , axes = ( 0 , 1 ) , direction = 'FFTW_FORWARD' , flags = [ 'FFTW_ESTIMATE' ] , threads = multiprocessing . cpu_count ( ) ) plan ( image_padded , out_img ) plan = pyfftw . FFTW ( fh_c , out_fh , axes = ( 0 , ) , direction = 'FFTW_FORWARD' , flags = [ 'FFTW_ESTIMATE' ] , threads = multiprocessing . cpu_count ( ) ) plan ( fh_c , out_fh ) plan = pyfftw . FFTW ( fv , out_fv , axes = ( 0 , ) , direction = 'FFTW_FORWARD' , flags = [ 'FFTW_ESTIMATE' ] , threads = multiprocessing . cpu_count ( ) ) plan ( fv , out_fv ) out_img *= out_fh [ : , None ] out_img *= out_fv [ None , : ] conv = image_padded plan = pyfftw . FFTW ( out_img . copy ( ) , conv , axes = ( 0 , 1 ) , direction = 'FFTW_BACKWARD' , flags = [ 'FFTW_ESTIMATE' ] , threads = multiprocessing . cpu_count ( ) ) plan ( out_img , conv ) else : raise ValueError ( 'unsupported `impl` {!r}' . format ( impl_in ) ) if padding : return conv [ padding : - padding , padding : - padding ] else : return conv
Filter an image with a separable filter.

    Parameters
    ----------
    image : 2D array-like
        The image to be filtered. It must have a real (vs. complex) dtype.
    fh, fv : 1D array-like
        Horizontal (axis 0) and vertical (axis 1) filters. Their sizes
        can be at most the image sizes in the respective axes.
    impl : {'numpy', 'pyfftw'}, optional
        FFT backend to use. The ``pyfftw`` backend requires the
        ``pyfftw`` package to be installed. It is usually significantly
        faster than the NumPy backend.
    padding : positive int, optional
        Amount of zeros added to the left and right of the image in all
        axes before FFT. This helps avoiding wraparound artifacts due to
        large boundary values.
        For ``None``, the padding is computed as ::

            padding = min(max(len(fh), len(fv)) - 1, 64)

        A padding of ``len(filt) - 1`` ensures that errors in FFT-based
        convolutions are small. At the same time, the padding should not
        be excessive to retain efficiency.

    Returns
    -------
    filtered : 2D `numpy.ndarray`
        The image filtered horizontally by ``fh`` and vertically by ``fv``.
        It has the same shape as ``image``, and its dtype is
        ``np.result_type(image, fh, fv)``.