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