code: def porosity_profile(im, axis):
      r
      if axis >= im.ndim:
          raise Exception('axis out of range')
      im = np.atleast_3d(im)
      a = set(range(im.ndim)).difference(set([axis]))
      a1, a2 = a
      prof = np.sum(np.sum(im, axis=a2), axis=a1)/(im.shape[a2]*im.shape[a1])
      return prof*100
docs: r"""
      Returns a porosity profile along the specified axis
  
      Parameters
      ----------
      im : ND-array
          The volumetric image for which to calculate the porosity profile
      axis : int
          The axis (0, 1, or 2) along which to calculate the profile.  For
          instance, if `axis` is 0, then the porosity in each YZ plane is
          calculated and returned as 1D array with 1 value for each X position.
  
      Returns
      -------
      result : 1D-array
          A 1D-array of porosity along the specified axis