docs: def associate_by_distance(labels_a, labels_b, distance):
      '''Find the objects that are within a given distance of each other
      
      Given two labels matrices and a distance, find pairs of objects that
      are within the given distance of each other where the distance is
      the minimum distance between any point in the convex hull of the
      two objects.
      
      labels_a - first labels matrix
      labels_b - second labels matrix
      distance - distance to measure
      
      returns a n x 2 matrix where m[x,0] is the label number in labels1 and
      m[x,1] is the label number in labels2
      
      Algorithm for computing distance between convex polygons taken from
      Chin, "Optimal Algorithms for the Intersection and the Minimum Distance 
      Problems Between Planar Polygons", IEEE Transactions on Computers, 
      vol. C-32, # 12, December 1983
      '''
      if np.max(labels_a) == 0 or np.max(labels_b) == 0:
          return np.zeros((0,2),int)
      
      hull_a, point_counts_a = convex_hull(labels_a)
      hull_b, point_counts_b = convex_hull(labels_b)
      centers_a, radii_a = minimum_enclosing_circle(
          labels_a, hull_and_point_count = (hull_a, point_counts_a))
      centers_b, radii_b = minimum_enclosing_circle(
          labels_b, hull_and_point_count = (hull_b, point_counts_b))
      #
      # Make an indexer into the hull tables
      #
      indexer_a = np.cumsum(point_counts_a)
      indexer_a[1:] = indexer_a[:-1]
      indexer_a[0] = 0
      indexer_b = np.cumsum(point_counts_b)
      indexer_b[1:] = indexer_b[:-1]
      indexer_b[0] = 0
      #
      # Compute the distances between minimum enclosing circles =
      # distance - radius_a - radius_b
      #
      i,j = np.mgrid[0:len(radii_a),0:len(radii_b)]
      ab_distance = np.sqrt((centers_a[i,0]-centers_b[j,0])**2 +
                            (centers_a[i,1]-centers_b[j,1])**2)
      ab_distance_minus_radii = ab_distance - radii_a[i] - radii_b[j]
      # Account for roundoff error
      ab_distance_minus_radii -= np.sqrt(np.finfo(float).eps)
      #
      # Exclude from consideration ab_distance > distance and automatically
      # choose those whose centers are within the distance
      #
      ab_easy_wins = ab_distance <= distance
      ij_wins = np.dstack((hull_a[indexer_a[i[ab_easy_wins]],0], 
                           hull_b[indexer_b[j[ab_easy_wins]],0]))
      ij_wins.shape = ij_wins.shape[1:]
      
      ab_consider = (ab_distance_minus_radii <= distance) & (~ ab_easy_wins)
      ij_consider = np.dstack((i[ab_consider], j[ab_consider]))
      ij_consider.shape = ij_consider.shape[1:]
      if np.product(ij_consider.shape) == 0:
          return ij_wins
      if True:
          wins = []
          distance2 = distance**2
          for ii,jj in ij_consider:
              a = hull_a[indexer_a[ii]:indexer_a[ii]+point_counts_a[ii],1:]
              b = hull_b[indexer_b[jj]:indexer_b[jj]+point_counts_b[jj],1:]
              d = minimum_distance2(a,centers_a[ii,:],
                                    b,centers_b[jj,:])
              if d <= distance2:
                  wins.append((hull_a[indexer_a[ii],0],
                               hull_b[indexer_b[jj],0]))
          ij_wins = np.vstack((ij_wins, np.array(wins)))
          return ij_wins
      else:
          #
          # For each point in the hull, get the next point mod # of points in hull
          #
          hull_next_a = np.arange(hull_a.shape[0])+1
          hull_next_a[indexer_a+point_counts_a-1] = indexer_a
          hull_next_b = np.arange(hull_b.shape[0])+1
          hull_next_b[indexer_b+point_counts_b-1] = indexer_b
          #
          # Parallelize the algorithm for overlap
          #
          # For each pair of points i, i+1 mod n in the hull, and the test point t
          # the cross product of the vector from i to i+1 and the vector from i+1
          # to t should have the same sign.
          #
          next_b = hull_b[hull_next_b,1:]
          vector_b = hull_b[:,1:] - next_b
          #
          # For each i,j, we have to compare the centers_a against point_counts_b[j]
          # crosses.
          #
          b_len = point_counts_b[ij_consider[:,1]]
          b_index = np.cumsum(point_counts_b)
          b_elems = b_index[-1]
          b_index[1:] = b_index[:-1]
          b_index[0] = 0
          #
          # First create a vector that's b_elems long and every element contains an
          # index into the ij_consider vector. How we do this:
          # 1) mark the first element at a particular index by 1, all others = 0
          # 2) Erase the first 1
          # 3) Take the cumulative sum which will increment to 1 when it hits the
          #    first 1, again when it hits the second...etc.
          #
          b_indexer = np.zeros(b_elems, int)
          b_indexer[b_index[1:]] = 1
          b_indexer = np.cumsum(b_indexer)
          #
          # The sub-index is the index from 1 to n for each of the vertices
          # per b convex hull
          #
          b_sub_index = np.arange(b_elems) - b_index[b_indexer]
          #
          # For each element of b_indexer, get the i and j at that index
          #
          b_i = ij_consider[b_indexer,0]
          b_j = ij_consider[b_indexer,1]
          #
          # Compute the cross-products now
          #
          b_vector_b = vector_b[indexer_b[b_j]+b_sub_index,:]
          b_center_vector = (next_b[indexer_b[b_j]+b_sub_index,:] - 
                             centers_a[indexer_a[b_i]])
          cross = (b_vector_b[:,0] * b_center_vector[:,1] -
                   b_vector_b[:,1] * b_center_vector[:,0])
          hits = (all_true(cross > 0, b_index) | all_true(cross < 0, b_index))
          
          ij_wins = np.vstack((ij_wins, ij_consider[hits,:]))
          ij_consider = ij_consider[~hits,:]
          if ij_consider.shape[0] == 0:
              return ij_wins