[PIPE2D-218] Discuss FiberTrace's architecture Created: 13/Jul/17  Updated: 18/Aug/17  Resolved: 05/Aug/17

Status: Done
Project: DRP 2-D Pipeline
Component/s: None
Affects Version/s: None
Fix Version/s: None

Type: Task Priority: Normal
Reporter: rhl Assignee: rhl
Resolution: Done Votes: 0
Labels: None
Remaining Estimate: Not Specified
Time Spent: Not Specified
Original Estimate: Not Specified

Issue Links:
Relates
relates to PIPE2D-219 FiberTrace.createTrace() can create a... Done
relates to PIPE2D-235 Change sliced FiberTrace Image/Profil... Done
relates to PIPE2D-236 Get rid of non-public functions in th... Done
relates to PIPE2D-240 Move FiberTrace::extractFromProfile a... Done

 Description   

Having just been debugging some issues with FiberTrace I got the impression that it could be made significantly simpler and more robust.

We should discuss this.



 Comments   
Comment by rhl [ 27/Jul/17 ]

Some of the comments about unused code will be fixed on PIPE2D-207 (or, in one case, PIPE2D-221)

Here's a first dump of my thoughts.

The FiberTrace class describes the profiles and positions of spectra on a detector.

Creating a FiberTraceSet is done as follows:

    ftffc = drpStella.FiberTraceFunctionFindingControl()     // terrible name!
    quartzImage = butler.get("calexp", dataId).getMaskedImage()

    fiberTraceSet = drpStella.findAndTraceApertures(quartzImage, ftffc)

    fiberTraceSet.sortTracesByXCenter()

    for i in range(fiberTraceSet.size()):
        fiberTraceSet.getFiberTrace(i).setITrace(i)

    fiberTraceSet.setFiberTraceProfileFittingControl(fiberTraceProfileFittingControl)

    if doAll:
        inFiberTraceSet.calcProfileAllTraces()
    else:
        for i in inTraceNumbers:
	    inFiberTraceSet.getFiberTrace(i).calcProfile()

Extracting a spectrum then consists of:

   # Read a model of the mapping from pixels to traces
   xCenters, wavelengths, traceIds = readWavelengthFile(wLenFile)

   # assign trace number to fiberTraceSet based on model
   drpStella.assignITrace(fiberTraceSet, traceIds, xCenters)
   #
   spectrumSet = drpStella.SpectrumSet()

   scienceImage = butler.get("calexp", dataId).getMaskedImage()

   for i in traceIds:
      fiberTrace = fiberTraceSet.getFiberTrace(i)

      # Set pixels in FiberTrace from scienceImage
      fiberTrace.createTrace(scienceImage)
      
      # Extract spectrum from profile
      spectrum = fiberTrace.extractFromProfile()
      spectrumSet.addSpectrum(spectrum)

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

I'd expect this to look like:

    xCenters, wavelengths, fiberIds = readWavelengthFile(wLenFile)   # we'd actually read the DetectorMap
    detectorMap = DetectorMap(xCenters, wavelengths, fiberIds, ...)

    quartzImage = butler.get("calexp", dataId).getMaskedImage()

    ftffc = drpStella.FiberTraceFunctionFindingControl()
    ftpfc = fiberTraceProfileFittingControl

    fiberTraceSet = drpStella.findAndTraceApertures(quartzImage, detectorMap, ftffc, ftpfc)

which would now have profiles set.  Extracting a spectrum would then consist of:

    scienceImage = butler.get("calexp", dataId).getMaskedImage()
    spectrumSet = drpStella.extractSpectrumSet(scienceImage, fiberTraceSet)

Note in particular that we do not set the "trace" images in the FiberTrace; in fact, I don't think that
we need them at all.  We do need code to return a MaskedImage of the Profile (in the CCD coordinate system,
with XY0 set and mask bits to indicate the trace).

The FiberTrace does, of course, need to know all about its profile but I think that the current code is more
complex than is really necessary.

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

A FiberTrace has a number of data members:
    std::vector<PTR(std::vector<float>)> _overSampledProfileFitXPerSwath;
    std::vector<PTR(std::vector<float>)> _overSampledProfileFitYPerSwath;
    std::vector<PTR(std::vector<float>)> _profileFittingInputXPerSwath;
    std::vector<PTR(std::vector<float>)> _profileFittingInputYPerSwath;
    std::vector<PTR(std::vector<float>)> _profileFittingInputXMeanPerSwath;
    std::vector<PTR(std::vector<float>)> _profileFittingInputYMeanPerSwath;
What are these?  2-D arrays?  why are they vectors of vectors, and what do they all mean?
    
    PTR(lsst::afw::image::MaskedImage<ImageT, MaskT, VarianceT>) _trace;
Why do we need this in the data structure?
    
    PTR(lsst::afw::image::Image<float>) _profile;
Why is this an image, not a MaskedImage?

    ndarray::Array<float, 2, 1> _xCentersMeas;
What are these?  Do they need to be in the class?
    ndarray::Array<float, 1, 1> _xCenters;
In what sense are these not measured?

    std::size_t _iTrace;
Is this a fiberId?

    bool _isTraceSet;
    bool _isProfileSet;
    bool _isFiberTraceProfileFittingControlSet;
This will go once the profile extraction is part of the trace generation

    PTR(const FiberTraceFunction) _fiberTraceFunction;
This isn't a function, it's a set of parameters.  What do the coefficients mean -- there seems to
be only one per trace.
    PTR(FiberTraceProfileFittingControl) _fiberTraceProfileFittingControl;
There seems to be lots of old code in here.  What's the telluric code doing -- we'll never turn it on.  If it simplifies the code we should delete all references to Piskunov.

And methods:
    // ctors/dtors
    explicit FiberTrace(std::size_t width = 0, std::size_t height = 0, std::size_t iTrace = 0);
Do we need a default ctor?
    explicit FiberTrace(PTR(const MaskedImageT) const& maskedImage,
                        PTR(const FiberTraceFunction) const& fiberTraceFunction, 
                        std::size_t iTrace=0);
    FiberTrace(FiberTrace<ImageT, MaskT, VarianceT> const& fiberTrace);
    FiberTrace(FiberTrace<ImageT, MaskT, VarianceT> & fiberTrace, bool const deep);
// why two, not a default argument (here or in the pybind11 bindings)?
// do we ever need to copy FiberTraces?

    virtual ~FiberTrace() {}

    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    // getters and setters
    // "trace"
    PTR(MaskedImageT) getTrace();
    const PTR(const MaskedImageT) getTrace() const;
    void setTrace(PTR(MaskedImageT) & trace);
    void createTrace(PTR(const MaskedImageT) const& maskedImage);
I think all of these can go.

    // "profile"
    PTR(lsst::afw::image::Image<float>) getProfile() const;
This needs to be in CCD coordinates, not internal "sliced" images.
    void setProfile( PTR(lsst::afw::image::Image<float>) const& profile);
When do we do this?
    
    // ID number of this trace (_iTrace) to this number
    void setITrace(const std::size_t iTrace);
Shouldn't be needed
    std::size_t getITrace() const;
ITrace -> FiberId

    // dimensions of this FiberTrace
    std::size_t getWidth() const;
    std::size_t getHeight() const;
What are these used for?

    // Check if things are set
    bool isTraceSet() const;
    bool isProfileSet() const;
    bool isFiberTraceProfileFittingControlSet() const;
All of these must always be set in the construction step.

    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

    // Extract the spectrum of this fiber trace using the _profile
    PTR(SpectrumT) extractFromProfile();
No;  this is a free function that takes an image and a fiber trace.
    
    // Simple Sum Extraction of this fiber trace
    PTR(SpectrumT) extractSum();
No; this should be an option to the free function.

    // FiberTraceFunction (const)
    const PTR(const FiberTraceFunction) getFiberTraceFunction() const;
    void setFiberTraceFunction(PTR(FiberTraceFunction const) fiberTraceFunction);
Why do you need these?  Maybe the getter

    // Return FiberTraceProfileFittingControl
    PTR(FiberTraceProfileFittingControl) getFiberTraceProfileFittingControl() const;
I'm not sure you need the getter
    void setFiberTraceProfileFittingControl(PTR(FiberTraceProfileFittingControl) const&);
    void setFiberTraceProfileFittingControl(PTR(FiberTraceProfileFittingControl const) const&);
Why do you need both of these?
    
    // fitted x-centers of the fiber trace
    const ndarray::Array< float, 1, 1 > getXCenters() const;
OK
    void setXCenters( ndarray::Array< float, 1, 1 > const& xCenters );
Why do you need this?

    // Return the measured x-centers of the fiber trace
    ndarray::Array< float, 2, 1 > getXCentersMeas() const;
    void setXCentersMeas( ndarray::Array< float, 2, 1 > const& );
How is this different from getXCenters?

    // Return an image containing the reconstructed 2D spectrum of the FiberTrace
    // @param spectrum : 1D spectrum to reconstruct the 2D image from
    PTR(Image) getReconstructed2DSpectrum(SpectrumT const& spectrum) const;
Should be a MaskedImage, I think.  In CCD coords, of course.

    // Return an image containing the reconstructed background of the FiberTrace
    // @param backgroundSpectrum : 1D spectrum to reconstruct the 2D image from
    PTR(Image) getReconstructedBackground(SpectrumT const& backgroundSpectrum) const;
Why do you need this too?

    // Return an image containing the reconstructed 2D spectrum + background of the FiberTrace
    // @param spectrum : 1D spectrum to use for the 2D reconstruction
    // @param background : 1D background spectrum to use for the 2D reconstruction
    PTR(Image) getReconstructed2DSpectrum(SpectrumT const& spectrum,
                                          SpectrumT const& background) const;
Why do you need this too?
    
    /**
     * @brief Calculate the spatial profile for the FiberTrace
     * Normally this would be a Flat FiberTrace, but in principle, if the spectrum
     * shows some kind of continuum, the spatial profile can still be calculated
     */
    void calcProfile();
What does this do?

    /**
     * @brief Helper function for calcProfile, calculates profile for a swath
     * A swath is approximately FiberTraceProfileFittingControl.swathWidth long
     * Each swath is overlapping the previous swath for half of the swath width
     * spectrum:
     * |-----------------------------------------------------------------
     * swaths:
     * |---------------|--------------|--------------|--------------|----
     *         |---------------|--------------|--------------|-----------
     * @param imageSwath : array containing the CCD image of the FiberTrace swath
     * @param maskSwath : array containing the mask of the FiberTrace swath
     * @param varianceSwath : array containing the variance of the FiberTrace swath
     * @param xCentersSwath : 1D array containing the x center positions for the swath
     * @param iSwath : number of swath
     */
    ndarray::Array<float, 2, 1> calcProfileSwath(ndarray::Array<ImageT const, 2, 1> const& imageSwath,
                                                 ndarray::Array<MaskT const, 2, 1> const& maskSwath,
                                                 ndarray::Array<VarianceT const, 2, 1> const& varianceSwath,
                                                 ndarray::Array<float const, 1, 1> const& xCentersSwath,
                                                 std::size_t const iSwath);
Why is this part of the public interface?  Those image swaths should be invisible to the user.

    /**
     * @brief Calculate boundaries for the swaths used for profile calculation
     * @param swathWidth_In : Approximate width for the swaths, will be adjusted
     * to fill the length of the FiberTrace with equally sized swaths
     * @return 2D array containing the pixel numbers for the start and the end
     * of each swath
     */
    ndarray::Array<std::size_t, 2, 1> calcSwathBoundY(const std::size_t swathWidth_In) const;
Ditto

    // Return the coefficients of the trace function of the xCenters
    ndarray::Array<float, 1, 1> getTraceCoefficients() const;
Why do we need these?  I can see a function that returns the model

    // Return a smart pointer to this FiberTrace
    PTR(FiberTrace) getPointer();
Why do we need this?

-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

  template<typename ImageT, typename VarianceT=lsst::afw::image::VariancePixel>
  FindCenterPositionsOneTraceResult findCenterPositionsOneTrace(
	PTR(lsst::afw::image::Image<ImageT>) & ccdImage,
        PTR(lsst::afw::image::Image<VarianceT>) & ccdImageVariance,
        PTR(const FiberTraceFunctionFindingControl) const& fiberTraceFunctionFindingControl);

should be in anon namespace

  
  ndarray::Array<float, 1, 1> calculateXCenters(PTR(const ::pfs::drp::stella::FiberTraceFunction) const& fiberTraceFunctionIn,
                                                 std::size_t const& ccdHeightIn = 0,
                                                 std::size_t const& ccdWidthIn = 0);
  ndarray::Array<float, 1, 1> calculateXCenters(PTR(const ::pfs::drp::stella::FiberTraceFunction) const& fiberTraceFunctionIn,
                                                 ndarray::Array<float, 1, 1> const& yIn,
                                                 std::size_t const& ccdHeightIn = 0,
                                                 std::size_t const& ccdWidthIn = 0);
Why do we need both in the public API?  And why not a member function?

  template< typename ImageT, typename MaskT=lsst::afw::image::MaskPixel,
            typename VarianceT=lsst::afw::image::VariancePixel >
  PTR(FiberTrace< ImageT, MaskT, VarianceT >) makeNormFlatFiberTrace(
	PTR(const lsst::afw::image::MaskedImage< ImageT, MaskT, VarianceT >) const& maskedImage,
        PTR(const ::pfs::drp::stella::FiberTraceFunction) const& fiberTraceFunctionWide,
        PTR(const ::pfs::drp::stella::FiberTraceFunctionControl) const& fiberTraceFunctionControlNarrow,
        PTR(const ::pfs::drp::stella::FiberTraceProfileFittingControl) const& fiberTraceProfileFittingControl,
        ImageT minSNR = 100.,
        std::size_t iTrace = 0);
Unused

  template< typename ImageT, typename MaskT, typename VarianceT, typename T, typename U, int I >
  void assignITrace( FiberTraceSet< ImageT, MaskT, VarianceT > & fiberTraceSet,
                     ndarray::Array< T, 1, I > const& traceIds,
                     ndarray::Array< U, 1, I > const& xCenters );
Part of the construction (and supported by DetectorMap)

    void addFiberTraceToCcdArray( FiberTrace< ImageT, MaskT, VarianceT > const& fiberTrace,
                                  lsst::afw::image::Image< arrayT > const& fiberTraceRepresentation,
                                  ndarray::Array< ccdImageT, 2, dim > & ccdArray );
    void addFiberTraceToCcdImage( FiberTrace< ImageT, MaskT, VarianceT > const& fiberTrace,
                                  lsst::afw::image::Image< arrayT > const& fiberTraceRepresentation,
                                  lsst::afw::image::Image< ccdImageT > & ccdImage );
    void addArrayIntoArray( ndarray::Array< smallT, 2, I > const& smallArr,
                            ndarray::Array< std::size_t, 2, 1 > const& xMinMax,
                            std::size_t const& yMin,
                            ndarray::Array< bigT, 2, J > & bigArr );
We don't need these once we can return an image of the profile

    dataXY<CoordT> ccdToFiberTraceCoordinates(
        dataXY<CoordT> const& ccdCoordinates,
        pfs::drp::stella::FiberTrace<ImageT, MaskT, VarianceT> const& fiberTrace);
Unused

    dataXY<CoordT> fiberTraceCoordinatesRelativeTo(
        dataXY<CoordT> const& fiberTraceCoordinates,
        dataXY<CoordT> const& ccdCoordinatesCenter,
        pfs::drp::stella::FiberTrace<ImageT, MaskT, VarianceT> const& fiberTrace
    );
unused
}   

    template<typename T>
    const T* getRawPointer(const PTR(const T) & ptr);
Fortunately unused.  This is a Very bad idea.

     /**
      * @brief mark FiberTrace pixels in Mask image
      * @param fiberTrace : FiberTrace to mark in maskedImage's Mask
      * @param mask : mask to mark the FiberTrace in
      * @param value : value to Or into the FiberTrace mask
      */
     template< typename ImageT, typename MaskT, typename VarianceT >
     void markFiberTraceInMask( PTR( FiberTrace< ImageT, MaskT, VarianceT > ) const& fiberTrace,
                                PTR( lsst::afw::image::Mask< MaskT > ) const& mask,
                                MaskT value = 1);
Superceded by getting the MaskedImage back.
Comment by swinbank [ 05/Aug/17 ]

We discussed this 2017-07-28. Marking it as done.

Generated at Sat Feb 10 15:49:40 JST 2024 using Jira 8.3.4#803005-sha1:1f96e09b3c60279a408a2ae47be3c745f571388b.