How to get area of image with boundary conditions?

Question

Many functions in scipy.ndimage accept an optional mode=nearest|wrap|reflect|constant argument which determines how to handle cases in which the function needs some data from outside the area of the image (padding). The padding is handled internally by NI_ExtendLine() in native code.

Instead of running a ndimage function on padded data, I would like to just get the padded data using the same choice of padding modes as ndimage uses.

Here is an example (for mode=nearest only, assumes 2d image):

"""
Get padded data.  Returns numpy array with shape (y1-y0, x1-x0, ...)
Any of x0, x1, y0, y1 may be outside of the image
"""
def get(img, y0, y1, x0, x1, mode="nearest"):
    out_img = numpy.zeros((y1-y0, x1-x0))
    for y in range(y0, y1):
        for x in range(x0, x1):
            yc = numpy.clip(y, 0, img.shape[0])
            xc = numpy.clip(x, 0, img.shape[1])
            out_img[y-y0, x-x0] = img[yc, xc]
    return out_img

This does the right thing, but is slow since it iterates one pixel at a time.

What is the best (fastest, clearest, most pythonic) way to do this?

Answer 1

def get(img, y0, y1, x0, x1, mode="nearest"):
    xs, ys = numpy.mgrid[y0:y1, x0:x1]
    height, width = img.shape

    if mode == "nearest":
        xs = numpy.clip(xs, 0, height-1)
        ys = numpy.clip(ys, 0, width-1)

    elif mode == "wrap":
        xs = xs % height
        ys = ys % width

    elif mode == "reflect":
        maxh = height-1
        maxw = width-1

        # An unobvious way of performing reflecting modulo
        # You should comment this
        xs = numpy.absolute((xs + maxh) % (2 * maxh) - maxh)
        ys = numpy.absolute((ys + maxw) % (2 * maxw) - maxw)

    elif mode == "constant":
        output = numpy.empty((y1-y0, x1-x0))
        output.fill(0) # WHAT THE CONSTANT IS

        # LOADS of bounds checks and restrictions
        # You should comment this
        target_section = output[max(0, -y0):min(y1-y0, -y0+height), max(0, -x0):min(x1-x0, -x0+width)]
        new_fill = img[max(0, y0):min(height, y1), max(0, x0):min(width, x1)]

        # Crop both the sections, so that they're both the size of the smallest
        # Use more lines; I'm too lazy right now
        target_section[:new_fill.shape[0], :new_fill.shape[1]] = new_fill[:target_section.shape[0], :target_section.shape[1]]

        return output

    else:
        raise NotImplementedError("Unknown mode")

    return img[xs, ys]