Iterate object tracing

doc/object_finding.rst

@@ -138,4 +138,31 @@ vs spatial position vector. The best way to choose these pixels is to run PypeIt
 without it set. Then run :ref:`pypeit_show_2dspec` to view the sky-subtracted
 image and decide which pixels to use for object finding. Then re-run ``PypeIt``.
 
+Object Tracing
+==============
+
+For automatically identified objects (i.e., not manual extractions),
+PypeIt improves the object trace by fitting the spatial position of the peak
+as a function of wavelength. In some situations, the object trace is poorly
+determined by the peak location, and the code will fail to trace the object
+correctly. For example, if the slit edges are not well-defined, the object's
+position relative to the slit edges is also poorly defined, and the object trace
+is difficult to determine. In these cases, the user can attempt to iterate
+the object tracing by using a relatively low order polynomial in combination
+with an iterative fit, as follows
+
+.. code-block:: ini
 
+    [reduce]
+        [[findobj]]
+            find_numiterfit = 10
+            trace_npoly = 4
+
+Note that the default value is typically ``trace_npoly=5``. Also note that the
+object tracing will break if successive iterations yields the same object trace
+to within 0.1 pixels. If you notice a relatively poor object trace, sometimes in
+combination with the object counts being masked, an iterative fit may help to
+resolve your problem. If, on the other hand, your object is relatively faint, you
+may benefit from using the trace of a standard star (this is the default behaviour),
+and you can provide a 1D spectrum of a previously reduced standard star with the
+``std_spec1d`` parameter.

doc/releases/1.17.2dev.rst

@@ -17,6 +17,8 @@ Functionality/Performance Improvements and Additions
 - The WCS for datacubes now adopts the convention of North
   is up and East is left. In previous version of PypeIt,
   East was right.
+- The object trace can now be iterated over with the `find_numiterfit`
+  parameter.
 
 Instrument-specific Updates
 ---------------------------

pypeit/coadd2d.py

@@ -1416,6 +1416,7 @@ def compute_offsets(self, offsets):
                     inmask=inmask[iexp,:,:], fwhm=self.par['reduce']['findobj']['find_fwhm'],
                     trim_edg=self.par['reduce']['findobj']['find_trim_edge'],
                     maxdev=self.par['reduce']['findobj']['find_maxdev'],
+                    numiterfit=self.par['reduce']['findobj']['find_numiterfit'],
                     ncoeff=3, snr_thresh=self.par['reduce']['findobj']['snr_thresh'], nperslit=1,
                     find_min_max=self.par['reduce']['findobj']['find_min_max'],
                     show_trace=self.debug_offsets, show_peaks=self.debug_offsets)

pypeit/core/findobj_skymask.py

@@ -168,7 +168,7 @@ def ech_findobj_ineach_order(
     det='DET01', inmask=None, std_trace=None, ncoeff=5, 
     hand_extract_dict=None,
     box_radius=2.0, fwhm=3.0,
-    use_user_fwhm=False, maxdev=2.0, nperorder=2,
+    use_user_fwhm=False, maxdev=2.0, nperorder=2, numiterfit=1,
     extract_maskwidth=3.0, snr_thresh=10.0,
     specobj_dict=None, trim_edg=(5,5),
     show_peaks=False, show_single_fits=False,
@@ -260,6 +260,8 @@ def ech_findobj_ineach_order(
         maxdev (:obj:`float`, optional):
             Maximum deviation of pixels from polynomial fit to trace
             used to reject bad pixels in trace fitting.
+        numiterfit (:obj:`int`, optional):
+            Number of iterations to use when fitting the object traces.
         nperorder (:obj:`int`, optional):
             Maximum number of objects allowed per order.  If there are more
             detections than this number, the code will select the ``nperorder``
@@ -331,7 +333,7 @@ def ech_findobj_ineach_order(
                 spec_min_max=spec_min_max[:,iord],
                 inmask=inmask_iord,std_trace=std_in, 
                 ncoeff=ncoeff, fwhm=fwhm, use_user_fwhm=use_user_fwhm, maxdev=maxdev,
-                hand_extract_dict=hand_extract_dict,  
+                numiterfit=numiterfit, hand_extract_dict=hand_extract_dict,
                 nperslit=nperorder, extract_maskwidth=extract_maskwidth,
                 snr_thresh=snr_thresh, trim_edg=trim_edg, 
                 boxcar_rad=box_radius/plate_scale_ord[iord],
@@ -1038,7 +1040,7 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order
                 nabove_min_snr=2, pca_explained_var=99.0, 
                 box_radius=2.0, fwhm=3.0,
                 use_user_fwhm=False, maxdev=2.0, 
-                nperorder=2,
+                nperorder=2, numiterfit=1,
                 extract_maskwidth=3.0, snr_thresh=10.0,
                 specobj_dict=None, trim_edg=(5,5),
                 show_peaks=False, show_fits=False, 
@@ -1181,6 +1183,8 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order
         maxdev (:obj:`float`, optional):
             Maximum deviation of pixels from polynomial fit to trace
             used to reject bad pixels in trace fitting.
+        numiterfit (:obj:`int`, optional):
+            Number of iterations to perform when building the trace fits.
         hand_extract_dict (:obj:`dict`, optional):
             Dictionary with info on manual extraction; see
             :class:`~pypeit.manual_extract.ManualExtractionObj`.
@@ -1290,6 +1294,7 @@ def ech_objfind(image, ivar, slitmask, slit_left, slit_righ, slit_spat_id, order
         use_user_fwhm=use_user_fwhm,
         nperorder=nperorder,
         maxdev=maxdev,
+        numiterfit=numiterfit,
         box_radius=box_radius,
         objfindQA_filename=objfindQA_filename,
         hand_extract_dict=manual_extract_dict)
@@ -1500,7 +1505,7 @@ def get_fwhm(fwhm_in, nsamp, smash_peakflux, spat_fracpos, flux_smash_smth):
 def objs_in_slit(image, ivar, thismask, slit_left, slit_righ, 
                  inmask=None, fwhm=3.0,
                  sigclip_smash=5.0, use_user_fwhm=False, boxcar_rad=7.,
-                 maxdev=2.0, spec_min_max=None, hand_extract_dict=None, std_trace=None,
+                 maxdev=2.0, numiterfit=1, spec_min_max=None, hand_extract_dict=None, std_trace=None,
                  ncoeff=5, nperslit=None, snr_thresh=10.0, trim_edg=(5,5),
                  extract_maskwidth=4.0, specobj_dict=None, find_min_max=None,
                  show_peaks=False, show_fits=False, show_trace=False,
@@ -1588,6 +1593,10 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ,
         maxdev (:obj:`float`, optional):
             Maximum deviation of pixels from polynomial fit to trace
             used to reject bad pixels in trace fitting.
+        numiterfit (:obj:`int`, optional):
+            Number of iterations to use in the iterative trace fitting.
+            The number of iterations will be reduced if the fit has converged
+            to a value less than 0.1 pixels over the entire trace.
         spec_min_max (:obj:`tuple`, optional):
             2-tuple defining the minimum and maximum pixel in the spectral
             direction with useable data for this slit/order.  If None, the
@@ -1904,8 +1913,9 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ,
     msgs.info("Automatic finding routine found {0:d} objects".format(len(sobjs)))
 
     # Fit the object traces
+    tolerance = 0.1  # Tolerance in pixels before the trace fit has converged
     if len(sobjs) > 0:
-        msgs.info('Fitting the object traces')
+        msgs.info('Fitting the traces')
         # Note the transpose is here to pass in the TRACE_SPAT correctly.
         xinit_fweight = np.copy(sobjs.TRACE_SPAT.T).astype(float)
         spec_mask = (spec_vec >= spec_min_max_out[0]) & (spec_vec <= spec_min_max_out[1])
@@ -1914,9 +1924,16 @@ def objs_in_slit(image, ivar, thismask, slit_left, slit_righ,
                                  trace_bpm=np.invert(trc_inmask), fwhm=fwhm, maxdev=maxdev,
                                  idx=sobjs.NAME, debug=show_fits)[0]
         xinit_gweight = np.copy(xfit_fweight)
-        xfit_gweight = fit_trace(image, xinit_gweight, ncoeff, bpm=np.invert(inmask), maxshift=1.,
-                                 trace_bpm=np.invert(trc_inmask), fwhm=fwhm, maxdev=maxdev,
-                                 weighting='gaussian', idx=sobjs.NAME, debug=show_fits)[0]
+        for nn in range(numiterfit):
+            xfit_gweight = fit_trace(image, xinit_gweight, ncoeff, bpm=np.invert(inmask), maxshift=1.,
+                                     trace_bpm=np.invert(trc_inmask), fwhm=fwhm, maxdev=maxdev,
+                                     weighting='gaussian', idx=sobjs.NAME, debug=show_fits)[0]
+            maxdev = np.max(np.abs(xfit_gweight-xinit_gweight))
+            xinit_gweight = np.copy(xfit_gweight)
+            if maxdev <= tolerance:
+                break
+        if maxdev > tolerance:
+            msgs.warn('Trace fitting did not converge.  Check the fits.')
 
         # assign the final trace
         for iobj in range(nobj_reg):

pypeit/find_objects.py

@@ -819,6 +819,7 @@ def find_objects_pypeline(self, image, ivar, std_trace=None,
                                 use_user_fwhm=self.par['reduce']['extraction']['use_user_fwhm'],
                                 boxcar_rad=self.par['reduce']['extraction']['boxcar_radius'] / self.get_platescale(),  #pixels
                                 maxdev=self.par['reduce']['findobj']['find_maxdev'],
+                                numiterfit=self.par['reduce']['findobj']['find_numiterfit'],
                                 find_min_max=self.par['reduce']['findobj']['find_min_max'],
                                 extract_maskwidth=self.par['reduce']['skysub']['local_maskwidth'],
                                 qa_title=qa_title, nperslit=maxnumber,
@@ -959,6 +960,7 @@ def find_objects_pypeline(self, image, ivar, std_trace=None,
             use_user_fwhm=self.par['reduce']['extraction']['use_user_fwhm'],
             fof_link = self.par['reduce']['findobj']['fof_link'],
             maxdev=self.par['reduce']['findobj']['find_maxdev'],
+            numiterfit=self.par['reduce']['findobj']['find_numiterfit'],
             nperorder=nperorder,
             max_snr=self.par['reduce']['findobj']['ech_find_max_snr'],
             min_snr=self.par['reduce']['findobj']['ech_find_min_snr'],

pypeit/par/pypeitpar.py

@@ -4146,7 +4146,7 @@ class FindObjPar(ParSet):
 
     def __init__(self, trace_npoly=None, snr_thresh=None, find_trim_edge=None,
                  find_maxdev=None, find_extrap_npoly=None, maxnumber_sci=None, maxnumber_std=None,
-                 find_fwhm=None, ech_find_max_snr=None, ech_find_min_snr=None,
+                 find_fwhm=None, ech_find_max_snr=None, ech_find_min_snr=None, find_numiterfit=None,
                  ech_find_nabove_min_snr=None, skip_second_find=None, skip_final_global=None,
                  skip_skysub=None, find_negative=None, find_min_max=None, std_spec1d=None,
                  use_std_trace=None, fof_link = None):
@@ -4206,6 +4206,11 @@ def __init__(self, trace_npoly=None, snr_thresh=None, find_trim_edge=None,
         dtypes['find_maxdev'] = [int, float]
         descr['find_maxdev'] = 'Maximum deviation of pixels from polynomial fit to trace used to reject bad pixels in trace fitting.'
 
+        defaults['find_numiterfit'] = 1
+        dtypes['find_numiterfit'] = int
+        descr['find_numiterfit'] = 'Number of iterations to perform on the trace fitting. If the trace fitting ' \
+                                   'has converged, the number of iterations will be less than this number.'
+
         defaults['find_fwhm'] = 5.0
         dtypes['find_fwhm'] = [int, float]
         descr['find_fwhm'] = 'Indicates roughly the fwhm of objects in pixels for object finding'
@@ -4304,7 +4309,7 @@ def from_dict(cls, cfg):
         # Basic keywords
         parkeys = ['trace_npoly', 'snr_thresh', 'find_trim_edge',
                    'find_extrap_npoly', 'maxnumber_sci', 'maxnumber_std',
-                   'find_maxdev', 'find_fwhm', 'ech_find_max_snr',
+                   'find_maxdev', 'find_numiterfit', 'find_fwhm', 'ech_find_max_snr',
                    'ech_find_min_snr', 'ech_find_nabove_min_snr', 'skip_second_find',
                    'skip_final_global', 'skip_skysub', 'find_negative', 'find_min_max',
                    'std_spec1d', 'use_std_trace', 'fof_link']