<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def region_size(im):
 r""" Replace each voxel with size of region to which it belongs Parameters im : ND-array Either a boolean image wtih ``True`` indicating the features of interest, in which case ``scipy.ndimage.label`` will be applied to find regions, or a greyscale image with integer values indicating regions. Returns ------- image : ND-array A copy of ``im`` with each voxel value indicating the size of the region to which it belongs. This is particularly useful for finding chord sizes on the image produced by ``apply_chords``. """

    if im.dtype == bool:
        im = spim.label(im)[0]
    counts = sp.bincount(im.flatten())
    counts[0] = 0
    chords = counts[im]
    return chords