Compute the slice representation of intersection of two arrays.

    Given the shapes of two arrays and a reference point ref, compute the
    intersection of the two arrays.
    It returns a tuple of slices, that can be passed to the two \
    images as indexes

    :param shape: the shape of the reference array
    :param ref: coordinates of the reference point in the first array system
    :param sshape: the shape of the second array
    :param: sref: coordinates of the reference point in the \
    second array system, the origin by default
    :type sref: sequence or None
    :return: two matching slices, corresponding to both arrays \
    or a tuple of Nones if they don't match
    :rtype: a tuple

    Example:

      >>> import numpy
      >>> im = numpy.zeros((1000, 1000))
      >>> sim = numpy.ones((40, 40))
      >>> i,j = subarray_match(im.shape, [20, 23], sim.shape)
      >>> im[i] = 2 * sim[j]

def subarray_match(shape, ref, sshape, sref=None):
    """Compute the slice representation of intersection of two arrays.

    Given the shapes of two arrays and a reference point ref, compute the
    intersection of the two arrays.
    It returns a tuple of slices, that can be passed to the two \
    images as indexes

    :param shape: the shape of the reference array
    :param ref: coordinates of the reference point in the first array system
    :param sshape: the shape of the second array
    :param: sref: coordinates of the reference point in the \
    second array system, the origin by default
    :type sref: sequence or None
    :return: two matching slices, corresponding to both arrays \
    or a tuple of Nones if they don't match
    :rtype: a tuple

    Example:

      >>> import numpy
      >>> im = numpy.zeros((1000, 1000))
      >>> sim = numpy.ones((40, 40))
      >>> i,j = subarray_match(im.shape, [20, 23], sim.shape)
      >>> im[i] = 2 * sim[j]

    """
    # Reference point in im
    ref1 = asarray(ref, dtype='int')

    if sref is not None:
        ref2 = asarray(sref, dtype='int')
    else:
        ref2 = zeros_like(ref1)

    offset = ref1 - ref2
    urc1 = minimum(offset + asarray(sshape) - 1, asarray(shape) - 1)
    blc1 = maximum(offset, 0)
    urc2 = urc1 - offset
    blc2 = blc1 - offset

    def valid_slice(b, u):
        if b >= u + 1:
            return None
        else:
            return slice(b, u + 1)

    f = tuple(valid_slice(b, u) for b, u in zip(blc1, urc1))
    s = tuple(valid_slice(b, u) for b, u in zip(blc2, urc2))

    if not all(f) or not all(s):
        return (None, None)

    return (f, s)