docs: def _build_epsf_step(self, stars, epsf=None):
          """
          A single iteration of improving an ePSF.
  
          Parameters
          ----------
          stars : `EPSFStars` object
              The stars used to build the ePSF.
  
          epsf : `EPSFModel` object, optional
              The initial ePSF model.  If not input, then the ePSF will be
              built from scratch.
  
          Returns
          -------
          epsf : `EPSFModel` object
              The updated ePSF.
          """
  
          if len(stars) < 1:
              raise ValueError('stars must contain at least one EPSFStar or '
                               'LinkedEPSFStar object.')
  
          if epsf is None:
              # create an initial ePSF (array of zeros)
              epsf = self._create_initial_epsf(stars)
          else:
              # improve the input ePSF
              epsf = copy.deepcopy(epsf)
  
          # compute a 3D stack of 2D residual images
          residuals = self._resample_residuals(stars, epsf)
  
          self._residuals.append(residuals)
  
          # compute the sigma-clipped median along the 3D stack
          with warnings.catch_warnings():
              warnings.simplefilter('ignore', category=RuntimeWarning)
              warnings.simplefilter('ignore', category=AstropyUserWarning)
              residuals = self.sigclip(residuals, axis=0, masked=False,
                                       return_bounds=False)
              if HAS_BOTTLENECK:
                  residuals = bottleneck.nanmedian(residuals, axis=0)
              else:
                  residuals = np.nanmedian(residuals, axis=0)
  
          self._residuals_sigclip.append(residuals)
  
          # interpolate any missing data (np.nan)
          mask = ~np.isfinite(residuals)
          if np.any(mask):
              residuals = _interpolate_missing_data(residuals, mask,
                                                    method='cubic')
  
              # fill any remaining nans (outer points) with zeros
              residuals[~np.isfinite(residuals)] = 0.
  
          self._residuals_interp.append(residuals)
  
          # add the residuals to the previous ePSF image
          new_epsf = epsf.normalized_data + residuals
  
          # smooth the ePSF
          new_epsf = self._smooth_epsf(new_epsf)
  
          # recenter the ePSF
          new_epsf = self._recenter_epsf(new_epsf, epsf,
                                         centroid_func=self.recentering_func,
                                         box_size=self.recentering_boxsize,
                                         maxiters=self.recentering_maxiters,
                                         center_accuracy=1.0e-4)
  
          # normalize the ePSF data
          new_epsf /= np.sum(new_epsf, dtype=np.float64)
  
          # return the new ePSF object
          xcenter = (new_epsf.shape[1] - 1) / 2.
          ycenter = (new_epsf.shape[0] - 1) / 2.
  
          epsf_new = EPSFModel(data=new_epsf, origin=(xcenter, ycenter),
                               normalize=False, oversampling=epsf.oversampling)
  
          return epsf_new