docs: def is_local_maximum(image, labels, footprint):
      '''Return a boolean array of points that are local maxima
      
      image - intensity image
      labels - find maxima only within labels. Zero is reserved for background.
      footprint - binary mask indicating the neighborhood to be examined
                  must be a matrix with odd dimensions, center is taken to
                  be the point in question.
      '''
      assert((np.all(footprint.shape) & 1) == 1)
      footprint = (footprint != 0)
      footprint_extent = (np.array(footprint.shape)-1) // 2
      if np.all(footprint_extent == 0):
          return labels > 0
      result = (labels > 0).copy()
      #
      # Create a labels matrix with zeros at the borders that might be
      # hit by the footprint.
      #
      big_labels = np.zeros(np.array(labels.shape) + footprint_extent*2,
                            labels.dtype)
      big_labels[[slice(fe,-fe) for fe in footprint_extent]] = labels
      #
      # Find the relative indexes of each footprint element
      #
      image_strides = np.array(image.strides) // image.dtype.itemsize
      big_strides = np.array(big_labels.strides) // big_labels.dtype.itemsize
      result_strides = np.array(result.strides) // result.dtype.itemsize
      footprint_offsets = np.mgrid[[slice(-fe,fe+1) for fe in footprint_extent]]
      footprint_offsets = footprint_offsets[:, footprint]
      #
      # Order by distance, low to high and get rid of center pt.
      #
      d = np.sum(footprint_offsets **2, 0)
      footprint_offsets, d = footprint_offsets[:, d > 0], d[d > 0]
      footprint_offsets = footprint_offsets[:, np.lexsort([d])]
      
      fp_image_offsets = np.sum(image_strides[:, np.newaxis] *
                                footprint_offsets, 0)
      fp_big_offsets = np.sum(big_strides[:, np.newaxis] *
                              footprint_offsets, 0)
      #
      # Get the index of each labeled pixel in the image and big_labels arrays
      #
      indexes = np.mgrid[[slice(0,x) for x in labels.shape]][:, labels > 0]
      image_indexes = np.sum(image_strides[:, np.newaxis] * indexes, 0)
      big_indexes = np.sum(big_strides[:, np.newaxis] * 
                           (indexes + footprint_extent[:, np.newaxis]), 0)
      result_indexes = np.sum(result_strides[:, np.newaxis] * indexes, 0)
      #
      # Now operate on the raveled images
      #
      big_labels_raveled = big_labels.ravel()
      image_raveled = image.ravel()
      result_raveled = result.ravel()
      #
      # A hit is a hit if the label at the offset matches the label at the pixel
      # and if the intensity at the pixel is greater or equal to the intensity
      # at the offset.
      #
      for fp_image_offset, fp_big_offset in zip(fp_image_offsets, fp_big_offsets):
          same_label = (big_labels_raveled[big_indexes + fp_big_offset] ==
                        big_labels_raveled[big_indexes])
          less_than = (image_raveled[image_indexes[same_label]] <
                       image_raveled[image_indexes[same_label]+ fp_image_offset])
          mask = ~same_label
          mask[same_label] = ~less_than
          result_raveled[result_indexes[~mask]] = False
          result_indexes = result_indexes[mask]
          big_indexes = big_indexes[mask]
          image_indexes = image_indexes[mask]
          
      return result