docs: def _recenter_epsf(self, epsf_data, epsf, centroid_func=centroid_com,
                         box_size=5, maxiters=20, center_accuracy=1.0e-4):
          """
          Calculate the center of the ePSF data and shift the data so the
          ePSF center is at the center of the ePSF data array.
  
          Parameters
          ----------
          epsf_data : 2D `~numpy.ndarray`
              A 2D array containing the ePSF image.
  
          epsf : `EPSFModel` object
              The ePSF model.
  
          centroid_func : callable, optional
              A callable object (e.g. function or class) that is used to
              calculate the centroid of a 2D array.  The callable must
              accept a 2D `~numpy.ndarray`, have a ``mask`` keyword and
              optionally an ``error`` keyword.  The callable object must
              return a tuple of two 1D `~numpy.ndarray`\\s, representing
              the x and y centroids.  The default is
              `~photutils.centroids.centroid_com`.
  
          recentering_boxsize : float or tuple of two floats, optional
              The size (in pixels) of the box used to calculate the
              centroid of the ePSF during each build iteration.  If a
              single integer number is provided, then a square box will be
              used.  If two values are provided, then they should be in
              ``(ny, nx)`` order.  The default is 5.
  
          maxiters : int, optional
              The maximum number of recentering iterations to perform.
              The default is 20.
  
          center_accuracy : float, optional
              The desired accuracy for the centers of stars.  The building
              iterations will stop if the center of the ePSF changes by
              less than ``center_accuracy`` pixels between iterations.
              The default is 1.0e-4.
  
          Returns
          -------
          result : 2D `~numpy.ndarray`
              The recentered ePSF data.
          """
  
          # Define an EPSFModel for the input data.  This EPSFModel will be
          # used to evaluate the model on a shifted pixel grid to place the
          # centroid at the array center.
          epsf = EPSFModel(data=epsf_data, origin=epsf.origin, normalize=False,
                           oversampling=epsf.oversampling)
  
          epsf.fill_value = 0.0
          xcenter, ycenter = epsf.origin
  
          dx_total = 0
          dy_total = 0
          y, x = np.indices(epsf_data.shape, dtype=np.float)
  
          iter_num = 0
          center_accuracy_sq = center_accuracy ** 2
          center_dist_sq = center_accuracy_sq + 1.e6
          center_dist_sq_prev = center_dist_sq + 1
          while (iter_num < maxiters and
                 center_dist_sq >= center_accuracy_sq):
  
              iter_num += 1
  
              # extract a cutout from the ePSF
              slices_large, slices_small = overlap_slices(epsf_data.shape,
                                                          box_size,
                                                          (ycenter, xcenter))
              epsf_cutout = epsf_data[slices_large]
              mask = ~np.isfinite(epsf_cutout)
  
              # find a new center position
              xcenter_new, ycenter_new = centroid_func(epsf_cutout, mask=mask)
              xcenter_new += slices_large[1].start
              ycenter_new += slices_large[0].start
  
              # calculate the shift
              dx = xcenter - xcenter_new
              dy = ycenter - ycenter_new
              center_dist_sq = dx**2 + dy**2
              if center_dist_sq >= center_dist_sq_prev:  # don't shift
                  break
              center_dist_sq_prev = center_dist_sq
  
              # Resample the ePSF data to a shifted grid to place the
              # centroid in the center of the central pixel.  The shift is
              # always performed on the input epsf_data.
              dx_total += dx    # accumulated shifts for the input epsf_data
              dy_total += dy
              epsf_data = epsf.evaluate(x=x, y=y, flux=1.0,
                                        x_0=xcenter + dx_total,
                                        y_0=ycenter + dy_total,
                                        use_oversampling=False)
  
          return epsf_data