<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_surface_areas(regions, voxel_size=1, strel=None):
 r""" Extracts the surface area of each region in a labeled image. Optionally, it can also find the the interfacial area between all adjoining regions. Parameters regions : ND-array An image of the pore space partitioned into individual pore regions. Note that zeros in the image will not be considered for area calculation. voxel_size : scalar The resolution of the image, expressed as the length of one side of a voxel, so the volume of a voxel would be **voxel_size**-cubed. The default is 1. strel : array_like The structuring element used to blur the region. If not provided, then a spherical element (or disk) with radius 1 is used. See the docstring for ``mesh_region`` for more details, as this argument is passed to there. Returns ------- result : list A list containing the surface area of each region, offset by 1, such that the surface area of region 1 is stored in element 0 of the list. """

    print('_'*60)
    print('Finding surface area of each region')
    im = regions.copy()
    # Get 'slices' into im for each pore region
    slices = spim.find_objects(im)
    # Initialize arrays
    Ps = sp.arange(1, sp.amax(im)+1)
    sa = sp.zeros_like(Ps, dtype=float)
    # Start extracting marching cube area from im
    for i in tqdm(Ps):
        reg = i - 1
        if slices[reg] is not None:
            s = extend_slice(slices[reg], im.shape)
            sub_im = im[s]
            mask_im = sub_im == i
            mesh = mesh_region(region=mask_im, strel=strel)
            sa[reg] = mesh_surface_area(mesh)
    result = sa * voxel_size**2
    return result