docs: def gabor(image, labels, frequency, theta):
      '''Gabor-filter the objects in an image
  
      image - 2-d grayscale image to filter
      labels - a similarly shaped labels matrix
      frequency - cycles per trip around the circle
      theta - angle of the filter. 0 to 2 pi
  
      Calculate the Gabor filter centered on the centroids of each object
      in the image. Summing the resulting image over the labels matrix will
      yield a texture measure per object.
      '''
      #
      # The code inscribes the X and Y position of each pixel relative to
      # the centroid of that pixel's object. After that, the Gabor filter
      # for the image can be calculated per-pixel and the image can be
      # multiplied by the filter to get the filtered image.
      #
      nobjects = np.max(labels)
      if nobjects == 0:
          return image
      centers = centers_of_labels(labels)
      areas = fix(scind.sum(np.ones(image.shape),labels, np.arange(nobjects, dtype=np.int32)+1))
      mask = labels > 0
      i,j = np.mgrid[0:image.shape[0],0:image.shape[1]].astype(float)
      i = i[mask]
      j = j[mask]
      image = image[mask]
      lm = labels[mask] - 1
      i -= centers[0,lm]
      j -= centers[1,lm]
      sigma = np.sqrt(areas/np.pi) / 3.0
      sigma = sigma[lm]
      g_exp = 1000.0/(2.0*np.pi*sigma**2) * np.exp(-(i**2 + j**2)/(2*sigma**2))
      g_angle = 2*np.pi/frequency*(i*np.cos(theta)+j*np.sin(theta))
      g_cos = g_exp * np.cos(g_angle)
      g_sin = g_exp * np.sin(g_angle)
      #
      # Normalize so that the sum of the filter over each object is zero
      # and so that there is no bias-value within each object.
      #
      g_cos_mean = fix(scind.mean(g_cos,lm, np.arange(nobjects)))
      i_mean = fix(scind.mean(image, lm, np.arange(nobjects)))
      i_norm = image - i_mean[lm]
      g_sin_mean = fix(scind.mean(g_sin,lm, np.arange(nobjects)))
      g_cos -= g_cos_mean[lm]
      g_sin -= g_sin_mean[lm]
      g = np.zeros(mask.shape,dtype=np.complex)
      g[mask] = i_norm *g_cos+i_norm * g_sin*1j
      return g