def haarpsi_weight_map(img1, img2, axis):
    r"""Weighting map for directional features along an axis.

    Parameters
    ----------
    img1, img2 : array-like
        The images to compare. They must have equal shape.
    axis : {0, 1}
        Direction in which to look for edge similarities.

    Returns
    -------
    weight_map : `numpy.ndarray`
        The pointwise weight map. See Notes for details.

    Notes
    -----
    The pointwise weight map of associated with input images :math:`f_1, f_2`
    and axis :math:`k` is defined
    as

    .. math::
        \mathrm{W}_{f_1, f_2}^{(k)}(x) =
        \max \left\{
            \left|g_3^{(k)} \ast f_1 \right|(x),
            \left|g_3^{(k)} \ast f_2 \right|(x)
        \right\},

    see `[Rei+2016] <https://arxiv.org/abs/1607.06140>`_ equations (11)
    and (13).

    Here, :math:`g_3^{(k)}` is a Haar wavelet filter for scaling level 3
    that performs high-pass filtering in axis :math:`k` and low-pass
    filtering in the other axes. Such a filter can be computed as ::

        f_lo_level1 = [np.sqrt(2), np.sqrt(2)] # low-pass Haar filter
        f_hi_level1 = [-np.sqrt(2), np.sqrt(2)] # high-pass Haar filter
        f_lo_level3 = np.repeat(f_lo_level1, 4)
        f_hi_level3 = np.repeat(f_hi_level1, 4)

    References
    ----------
    [Rei+2016] Reisenhofer, R, Bosse, S, Kutyniok, G, and Wiegand, T.
    *A Haar Wavelet-Based Perceptual Similarity Index for Image Quality
    Assessment*. arXiv:1607.06140 [cs], Jul. 2016.
    """
    # TODO: generalize for nD
    impl = 'pyfftw' if PYFFTW_AVAILABLE else 'numpy'

    # Haar wavelet filters for level 3
    dec_lo_lvl3 = np.repeat([np.sqrt(2), np.sqrt(2)], 4)
    dec_hi_lvl3 = np.repeat([-np.sqrt(2), np.sqrt(2)], 4)

    if axis == 0:
        fh_lvl3 = dec_hi_lvl3
        fv_lvl3 = dec_lo_lvl3
    elif axis == 1:
        fh_lvl3 = dec_lo_lvl3
        fv_lvl3 = dec_hi_lvl3
    else:
        raise ValueError('`axis` out of the valid range 0 -> 1')

    # Filter with level 3 wavelet filter
    img1_lvl3 = filter_image_sep2d(img1, fh_lvl3, fv_lvl3, impl=impl)
    img2_lvl3 = filter_image_sep2d(img2, fh_lvl3, fv_lvl3, impl=impl)

    # Return the pointwise maximum of the filtered images
    np.abs(img1_lvl3, out=img1_lvl3)
    np.abs(img2_lvl3, out=img2_lvl3)

    return np.maximum(img1_lvl3, img2_lvl3)
r"""Weighting map for directional features along an axis.

    Parameters
    ----------
    img1, img2 : array-like
        The images to compare. They must have equal shape.
    axis : {0, 1}
        Direction in which to look for edge similarities.

    Returns
    -------
    weight_map : `numpy.ndarray`
        The pointwise weight map. See Notes for details.

    Notes
    -----
    The pointwise weight map of associated with input images :math:`f_1, f_2`
    and axis :math:`k` is defined
    as

    .. math::
        \mathrm{W}_{f_1, f_2}^{(k)}(x) =
        \max \left\{
            \left|g_3^{(k)} \ast f_1 \right|(x),
            \left|g_3^{(k)} \ast f_2 \right|(x)
        \right\},

    see `[Rei+2016] <https://arxiv.org/abs/1607.06140>`_ equations (11)
    and (13).

    Here, :math:`g_3^{(k)}` is a Haar wavelet filter for scaling level 3
    that performs high-pass filtering in axis :math:`k` and low-pass
    filtering in the other axes. Such a filter can be computed as ::

        f_lo_level1 = [np.sqrt(2), np.sqrt(2)] # low-pass Haar filter
        f_hi_level1 = [-np.sqrt(2), np.sqrt(2)] # high-pass Haar filter
        f_lo_level3 = np.repeat(f_lo_level1, 4)
        f_hi_level3 = np.repeat(f_hi_level1, 4)

    References
    ----------
    [Rei+2016] Reisenhofer, R, Bosse, S, Kutyniok, G, and Wiegand, T.
    *A Haar Wavelet-Based Perceptual Similarity Index for Image Quality
    Assessment*. arXiv:1607.06140 [cs], Jul. 2016.