def convolve_map(m, k, cpix, threshold=0.001, imin=0, imax=None, wmap=None):
    """
    Perform an energy-dependent convolution on a sequence of 2-D spatial maps.

    Parameters
    ----------

    m : `~numpy.ndarray`
       3-D map containing a sequence of 2-D spatial maps.  First
       dimension should be energy.

    k : `~numpy.ndarray`
       3-D map containing a sequence of convolution kernels (PSF) for
       each slice in m.  This map should have the same dimension as m.

    cpix : list
       Indices of kernel reference pixel in the two spatial dimensions.

    threshold : float
       Kernel amplitude 

    imin : int
       Minimum index in energy dimension.

    imax : int
       Maximum index in energy dimension.

    wmap :  `~numpy.ndarray`
       3-D map containing a sequence of 2-D spatial maps of weights.  First
       dimension should be energy. This map should have the same dimension as m.

    """
    islice = slice(imin, imax)

    o = np.zeros(m[islice, ...].shape)
    ix = int(cpix[0])
    iy = int(cpix[1])

    # Loop over energy
    for i in range(m[islice, ...].shape[0]):

        ks = k[islice, ...][i, ...]
        ms = m[islice, ...][i, ...]

        mx = ks[ix, :] > ks[ix, iy] * threshold
        my = ks[:, iy] > ks[ix, iy] * threshold

        nx = int(max(3, np.round(np.sum(mx) / 2.)))
        ny = int(max(3, np.round(np.sum(my) / 2.)))

        # Ensure that there is an odd number of pixels in the kernel
        # array
        if ix + nx + 1 >= ms.shape[0] or ix - nx < 0:
            nx -= 1
            ny -= 1

        sx = slice(ix - nx, ix + nx + 1)
        sy = slice(iy - ny, iy + ny + 1)

        ks = ks[sx, sy]

#        origin = [0, 0]
#        if ks.shape[0] % 2 == 0: origin[0] += 1
#        if ks.shape[1] % 2 == 0: origin[1] += 1
#        o[i,...] = ndimage.convolve(ms, ks, mode='constant',
#                                     origin=origin, cval=0.0)

        o[i, ...] = scipy.signal.fftconvolve(ms, ks, mode='same')
        if wmap is not None:
            o[i, ...] *= wmap[islice, ...][i, ...]

    return o