docs: def laplacian_of_gaussian(image, mask, size, sigma):
      '''Perform the Laplacian of Gaussian transform on the image
  
      image - 2-d image array
      mask  - binary mask of significant pixels
      size  - length of side of square kernel to use
      sigma - standard deviation of the Gaussian
      '''
      half_size = size//2
      i,j = np.mgrid[-half_size:half_size+1,
                     -half_size:half_size+1].astype(float) / float(sigma)
      distance = (i**2 + j**2)/2
      gaussian = np.exp(-distance)
      #
      # Normalize the Gaussian
      #
      gaussian = gaussian / np.sum(gaussian)
  
      log = (distance - 1) * gaussian
      #
      # Normalize the kernel to have a sum of zero
      #
      log = log - np.mean(log)
      if mask is None:
          mask = np.ones(image.shape[:2], bool)
      masked_image = image.copy()
      masked_image[~mask] = 0
      output = convolve(masked_image, log, mode='constant', cval=0)
      #
      # Do the LoG of the inverse of the mask. This finds the magnitude of the
      # contribution of the masked pixels. We then fudge by multiplying by the
      # value at the pixel of interest - this effectively sets the value at a
      # masked pixel to that of the pixel of interest.
      #
      # It underestimates the LoG, that's not a terrible thing.
      #
      correction = convolve((~ mask).astype(float), log, mode='constant', cval = 1)
      output += correction * image
      output[~ mask] = image[~ mask]
      return output