docs: 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]