code: def calc_surface_intensity(self, factor=10):
          
          # First we calculate the surface intensity at native resolution
          pixels = self.roi.pixels_interior
          nside_in = self.config['coords']['nside_pixel']
          surface_intensity = self.kernel.pdf(pixels.lon,pixels.lat)
  
          # Then we recalculate the surface intensity around the kernel
          # centroid at higher resolution
          for i in np.arange(1,5):
              # Select pixels within the region of interest
              nside_out = 2**i * nside_in
              radius = factor*np.degrees(hp.max_pixrad(nside_out))
              pix = ang2disc(nside_in,self.kernel.lon,self.kernel.lat,
                             radius,inclusive=True)
  
              # Select pix within the interior region of the ROI
              idx = ugali.utils.healpix.index_pix_in_pixels(pix,pixels)
              pix = pix[(idx >= 0)]; idx = idx[(idx >= 0)]
  
              # Reset the surface intensity for the subsampled pixels
              subpix = ugali.utils.healpix.ud_grade_ipix(pix,nside_in,nside_out)
              pix_lon,pix_lat = pix2ang(nside_out,subpix)
              surface_intensity[idx]=np.mean(self.kernel.pdf(pix_lon,pix_lat),axis=1)
  
          return surface_intensity
docs: Calculate the surface intensity for each pixel in the interior
          region of the ROI. Pixels are adaptively subsampled around the
          kernel centroid out to a radius of 'factor * max_pixrad'.
  
          Parameters:
          -----------
          factor : the radius of the oversample region in units of max_pixrad
  
          Returns:
          --------
          surface_intensity : the surface intensity at each pixel