I think it is worth it: the cython code computes the statistics for all classes in two array sweeps (one for the means and one for the variances). The only alternative I came up with is to iterate over each class (there may be hundreds of classes), create a mask of all voxels having that label and compute the statistics from it, something like this:

    for label in range(num_labels): #typically num_labels==256
        selected = (labels == label) #create the mask
        samplesize = selected.sum() #how many voxels have this label
        if samplesize == 0:#skip if empty
            continue
        means[label] = np.ma.masked_where(~selected, v).mean() 
        if samplesize == 1: #handle the special case of only one voxel 
            variances[label] = np.inf
        else:
            variances[label] = np.ma.masked_where(~selected, v).var()

Note that numpy returns 0 as the variance of a sample with only one element, but I am marking it as infinite. By any chance do you know a better way of computing masked statistics with multiple labels? it would be very useful! =)