def METHOD_NAME(lith_block, fault_block, n_faults, neighbors=8, noddy=False):
    """
    Args:
        lith_block (np.ndarray): Lithology block model
        fault_block (np.ndarray): Fault block model
        n_faults (int): Number of faults.
        neighbors (int, optional): Specifies the neighbor voxel connectivity taken into account for the topology
            analysis. Must be either 4 or 8 (default: 8)
        noddy (bool): If a noddy block is handed to the function, equalizes the results to be comparable with GemPy
    Returns:
        (np.ndarray): Model block with uniquely labeled regions.
    """
    lith_block = np.round(lith_block).astype(int)
    fault_block = np.round(fault_block).astype(int)
    # label the fault block for normalization (comparability of e.g. pynoddy and gempy models)
    fault_block = label(fault_block, neighbors=neighbors, background=9999)
    if noddy:
        # then this is a gempy model, numpy starts with 1
        lith_block[lith_block == 0] = int(np.max(lith_block) + 1)  # set the 0 to highest value + 1
        lith_block -= n_faults  # lower by n_faults to equal with pynoddy models
        # so the block starts at 1 and goes continuously to max
    ublock = (lith_block.max() + 1) * fault_block + lith_block
    labels_block, labels_n = label(ublock, neighbors=neighbors, return_num=True, background=9999)
    if 0 in np.unique(labels_block):
        labels_block += 1
    return labels_block