docs: def euler_number(labels, indexes=None):
      """Calculate the Euler number of each label
      
      labels - a label matrix
      indexes - the indexes of the labels to measure or None to
                treat the labels matrix as a binary matrix
      """
      if indexes is None:
          labels = labels != 0
          indexes = np.array([1],dtype=np.int32)
      elif getattr(indexes,'__getitem__',False):
          indexes = np.array(indexes,dtype=np.int32)
      else:
          indexes = np.array([indexes],dtype=np.int32)
      fix = fixup_scipy_ndimage_result
      #
      # The algorithm here is from the following reference:
      # S.B. Gray, "Local Properties of Binary Images in Two Dimensions",
      # IEEE Transactions on Computers, Vol c-20 # 5 p 551, May 1971
      #
      # The general idea is that crossings into objects can be measured locally
      # through counting numbers of patterns resulting in crossings. There
      # are three sets that are applicable in Euler Numbers:
      # Q1: 1 0  0 1  0 0  0 0 (or more simply, 1 bit per quad)
      #     0 0  0 0  1 0  0 1
      #
      # Q3: 0 1  1 0  1 1  1 1 (or 3 bits per quad)
      #     1 1  1 1  1 0  0 1
      #
      # QD: 1 0  0 1
      #     0 1  1 0
      #
      # and the Euler number = W of an object is
      #
      # 4W = n(Q1) - n(Q3) - 2n(QD) (equation 34)
      # W  = (n(Q1) - n(Q3) - 2n(QD))/4
      #
      # We shift the label matrix to make matrices, padded by zeros on the
      # sides for each of the four positions of the quad:
      # I00 I01
      # I10 I11
      # 
      # We can then assign each bitquad to a label based on the value
      # of the label at one of the "on" bits. For example, the first pattern
      # of Q1 has the label I00 because that bit is on. It's truth value is
      # I00 != I01 and I00 != I02 and I00 != I03.
      #
      I_shape = (labels.shape[0]+3,labels.shape[1]+3)
      I00 = np.zeros(I_shape,int)
      I01 = np.zeros(I_shape,int)
      I10 = np.zeros(I_shape,int)
      I11 = np.zeros(I_shape,int)
      slice_00 = [slice(1,labels.shape[0]+1),slice(1,labels.shape[1]+1)]
      slice_01 = [slice(1,labels.shape[0]+1),slice(0,labels.shape[1])]
      slice_10 = [slice(labels.shape[0]),slice(1,labels.shape[1]+1)]
      slice_11 = [slice(0,labels.shape[0]),slice(0,labels.shape[1])]
      I00[slice_00] = labels
      I01[slice_01] = labels
      I10[slice_10] = labels
      I11[slice_11] = labels
      #
      # There are 6 binary comparisons among the four bits
      #
      EQ00_01 = I00 == I01;              EQ01_00 = EQ00_01
      EQ00_10 = I00 == I10;              EQ10_00 = EQ00_10
      EQ00_11 = I00 == I11;              EQ11_00 = EQ00_11
      EQ01_10 = I01 == I10;              EQ10_01 = EQ01_10
      EQ01_11 = I01 == I11;              EQ11_01 = EQ01_11
      EQ10_11 = I10 == I11;              EQ11_10 = EQ10_11
      NE00_01 = np.logical_not(EQ00_01); NE01_00 = NE00_01
      NE00_10 = np.logical_not(EQ00_10); NE10_00 = NE00_10
      NE00_11 = np.logical_not(EQ00_11); NE11_00 = NE00_11
      NE01_10 = np.logical_not(EQ01_10); NE10_01 = NE01_10
      NE01_11 = np.logical_not(EQ01_11); NE11_01 = NE01_11
      NE10_11 = np.logical_not(EQ10_11); NE11_10 = NE10_11
      #
      # Q1: 1 0 
      #     0 0
      Q1_condition = (NE00_01 & NE00_10 & NE00_11).astype(int)
      #     0 1
      #     0 0
      Q1_condition[slice_00] += (NE01_00 & NE01_10 & NE01_11)[slice_01]
      #     0 0
      #     1 0
      Q1_condition[slice_00] += (NE10_00 & NE10_01 & NE10_11)[slice_10]
      #     0 0
      #     0 1
      Q1_condition[slice_00] += (NE11_00 & NE11_01 & NE11_10)[slice_11]
      Q1 = fix(scind.sum(Q1_condition, I00, indexes))
      #
      # Q3: 1 1
      #     1 0
      Q3_condition = (EQ00_10 & EQ00_01 & NE00_11).astype(int)
      #     0 1
      #     1 1
      Q3_condition[slice_00] += (NE11_00 & EQ11_10 & EQ11_01)[slice_11]
      #     1 0
      #     1 1
      Q3_condition += (NE00_01 & EQ00_10 & EQ00_11)
      #     1 1
      #     0 1
      Q3_condition += (NE00_10 & EQ00_01 & EQ00_11)
      Q3 = fix(scind.sum(Q3_condition, I00, indexes))
      # QD: 1 0
      #     0 1
      QD_condition = (NE00_01 & NE00_10 & EQ00_11).astype(int)
      #     0 1
      #     1 0
      QD_condition[slice_00] += (NE01_00 & NE01_11 & EQ01_10)[slice_01]
      QD  = fix(scind.sum(QD_condition, I00, indexes))
      W = (Q1 - Q3 - 2*QD).astype(float)/4.0
      if indexes is None:
          return W[0]
      return W