def overlapped_convolution(bin_template, bin_image,
                           tollerance=0.5, splits=(4, 2)):
    """
    As each of these images are hold only binary values, and RFFT2 works on
    float64 greyscale values, we can make the convolution more efficient by
    breaking the image up into :splits: sectons. Each one of these sections
    then has its greyscale value adjusted and then stacked.

    We then apply the convolution to this 'stack' of images, and adjust the
    resultant position matches.
    """
    th, tw = bin_template.shape
    ih, iw = bin_image.shape
    hs, ws = splits
    h = ih // hs
    w = iw // ws
    count = numpy.count_nonzero(bin_template)
    assert count > 0
    assert h >= th
    assert w >= tw

    yoffset = [(i * h, ((i + 1) * h) + (th - 1)) for i in range(hs)]
    xoffset = [(i * w, ((i + 1) * w) + (tw - 1)) for i in range(ws)]

    # image_stacks is Origin (x,y), array, z (height in stack)
    image_stacks = [((x1, y1), bin_image[y1:y2, x1:x2], float((count + 1) ** (num)))
                    for num, (x1, x2, y1, y2) in
                    enumerate((x1, x2, y1, y2) for (x1, x2)
                              in xoffset for (y1, y2) in yoffset)]

    pad_h = max(i.shape[0] for _, i, _ in image_stacks)
    pad_w = max(i.shape[1] for _, i, _ in image_stacks)

    # rfft metrics must be an even size - why ... maths?
    pad_w += pad_w % 2
    pad_h += pad_h % 2

    overlapped_image = sum_2d_images(pad_bin_image_to_shape(i, (pad_h, pad_w))
                                     * num for _, i, num in image_stacks)
    #print "Overlap splits %r, Image Size (%d,%d),
    #Overlapped Size (%d,%d)"  % (splits,iw,ih,pad_w,pad_h)

    # Calculate the convolution of the FFT's of the overlapped image & template
    convolution_freqs = (rfft2(overlapped_image) *
                         rfft2(bin_template[::-1, ::-1],
                               overlapped_image.shape))

    # Reverse the FFT to find the result overlapped image
    convolution_image = irfft2(convolution_freqs)
    # At this point, the maximum point in convolution_image should be the
    # bottom right (why?) of the area of greatest match

    results = set()
    for (x, y), _, num in image_stacks[::-1]:
        test = convolution_image / num
        filtered = ((test >= (count - tollerance)) &
                   (test <= (count + tollerance)))
        match_points = numpy.transpose(numpy.nonzero(filtered))  # bottom right
        for (fy, fx) in match_points:
            if fx < (tw - 1) or fy < (th - 1):
                continue
            results.add((x + fx - (tw - 1), y + fy - (th - 1)))
        convolution_image %= num
    return list(results)