code: def region_surface_areas(regions, voxel_size=1, strel=None):
      r
      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
docs: 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.