def grey_reconstruction(image, mask, footprint=None, offset=None):
    '''Perform a morphological reconstruction of the image
    
    grey_dilate the image, constraining each pixel to have a value that is
    at most that of the mask.
    image - the seed image
    mask - the mask, giving the maximum allowed value at each point
    footprint - a boolean array giving the neighborhood pixels to be used
                in the dilation. None = 8-connected
    
    The algorithm is taken from:
    Robinson, "Efficient morphological reconstruction: a downhill filter",
    Pattern Recognition Letters 25 (2004) 1759-1767
    '''
    assert tuple(image.shape) == tuple(mask.shape)
    assert np.all(image <= mask)
    if footprint is None:
        footprint = np.ones([3]*image.ndim, bool)
    else:
        footprint = footprint.copy()
        
    if offset is None:
        assert all([d % 2 == 1 for d in footprint.shape]),\
               "Footprint dimensions must all be odd"
        offset = np.array([d//2 for d in footprint.shape])
    # Cross out the center of the footprint
    footprint[[slice(d,d+1) for d in offset]] = False
    #
    # Construct an array that's padded on the edges so we can ignore boundaries
    # The array is a dstack of the image and the mask; this lets us interleave
    # image and mask pixels when sorting which makes list manipulations easier
    #
    padding = (np.array(footprint.shape)//2).astype(int)
    dims = np.zeros(image.ndim+1,int)
    dims[1:] = np.array(image.shape)+2*padding
    dims[0] = 2
    inside_slices = [slice(p,-p) for p in padding]
    values = np.ones(dims)*np.min(image)
    values[[0]+inside_slices] = image
    values[[1]+inside_slices] = mask
    #
    # Create a list of strides across the array to get the neighbors
    # within a flattened array
    #
    value_stride = np.array(values.strides[1:]) // values.dtype.itemsize
    image_stride = values.strides[0] // values.dtype.itemsize
    footprint_mgrid = np.mgrid[[slice(-o,d - o) 
                                for d,o in zip(footprint.shape,offset)]]
    footprint_offsets = footprint_mgrid[:,footprint].transpose()
    strides = np.array([np.sum(value_stride * footprint_offset)
                        for footprint_offset in footprint_offsets],
                       np.int32)
    values = values.flatten()
    value_sort = np.lexsort([-values]).astype(np.int32)
    #
    # Make a linked list of pixels sorted by value. -1 is the list terminator.
    #
    prev = -np.ones(len(values), np.int32)
    next = -np.ones(len(values), np.int32)
    prev[value_sort[1:]] = value_sort[:-1]
    next[value_sort[:-1]] = value_sort[1:]
    #
    # Create a rank-order value array so that the Cython inner-loop
    # can operate on a uniform data type
    #
    values, value_map = rank_order(values)
    current = value_sort[0]
    slow = False
    
    if slow:
        while current != -1:
            if current < image_stride:
                current_value = values[current]
                if current_value == 0:
                    break
                neighbors = strides+current
                for neighbor in neighbors:
                    if neighbor < 0:
                        raise IndexError("Index out of bounds: %d, current=%d, current_value=%d"%
                                         (neighbor,
                                          np.unravel_index(current, dims),
                                          current_value))
                    neighbor_value = values[neighbor]
                    # Only do neighbors less than the current value
                    if neighbor_value < current_value:
                        mask_value = values[neighbor + image_stride]
                        # Only do neighbors less than the mask value
                        if neighbor_value < mask_value:
                            # Raise the neighbor to the mask value if
                            # the mask is less than current
                            if mask_value < current_value:
                                link = neighbor + image_stride
                                values[neighbor] = mask_value
                            else:
                                link = current
                                values[neighbor] = current_value
                            # unlink the neighbor
                            nprev = prev[neighbor]
                            nnext = next[neighbor]
                            next[nprev] = nnext
                            if nnext != -1:
                                prev[nnext] = nprev
                            # link the neighbor after the link
                            next[neighbor] = next[link]
                            prev[neighbor] = link
                            prev[next[link]] = neighbor
                            next[link] = neighbor
            current = next[current]
    else:
        grey_reconstruction_loop(values, prev, next, strides, current, 
                                 image_stride)
    #
    # Reshape the values array to the shape of the padded image
    # and return the unpadded portion of that result
    #
    values = value_map[values[:image_stride]]
    values.shape = np.array(image.shape)+2*padding
    return values[inside_slices]