def fit_labels_to_mask(label_image, mask):
    r"""
    Reduces a label images by overlaying it with a binary mask and assign the labels
    either to the mask or to the background. The resulting binary mask is the nearest
    expression the label image can form of the supplied binary mask.
    
    Parameters
    ----------
    label_image : array_like
        A nD label map.
    mask : array_like
        A mask image, i.e., a binary image with False for background and True for foreground.
        
    Returns
    -------
    best_fit : ndarray
        The best fit of the labels to the mask.
    
    Raises
    ------         
    ValueError
        If ``label_image`` and ``mask`` are not of the same shape.
    """
    label_image = scipy.asarray(label_image)
    mask = scipy.asarray(mask, dtype=scipy.bool_)

    if label_image.shape != mask.shape:
        raise ValueError('The input images must be of the same shape.')
    
    # prepare collection dictionaries
    labels = scipy.unique(label_image)
    collection = {}
    for label in labels:
        collection[label] = [0, 0, []]  # size, union, points
    
    # iterate over the label images pixels and collect position, size and union
    for x in range(label_image.shape[0]):
        for y in range(label_image.shape[1]):
            for z in range(label_image.shape[2]):
                entry = collection[label_image[x,y,z]]
                entry[0] += 1
                if mask[x,y,z]: entry[1] += 1
                entry[2].append((x,y,z))
                
    # select labels that are more than half in the mask
    for label in labels:
        if collection[label][0] / 2. >= collection[label][1]:
            del collection[label]
                
    # image_result = numpy.zeros_like(mask) this is eq. to mask.copy().fill(0), which directly applied does not allow access to the rows and colums: Why?
    image_result = mask.copy()
    image_result.fill(False)         

    # add labels to result mask
    for label, data in list(collection.items()):
        for point in data[2]:
            image_result[point] = True
            
    return image_result