The FiberTrace class describes the profiles and positions of spectra on a detector.
Creating a FiberTraceSet is done as follows:
ftffc = drpStella.FiberTraceFunctionFindingControl() 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:
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);
virtual ~FiberTrace() {}
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
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.
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?
void setITrace(const std::size_t iTrace);
Shouldn't be needed
std::size_t getITrace() const;
ITrace -> FiberId
std::size_t getWidth() const;
std::size_t getHeight() const;
What are these used for?
bool isTraceSet() const;
bool isProfileSet() const;
bool isFiberTraceProfileFittingControlSet() const;
All of these must always be set in the construction step.
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
PTR(SpectrumT) extractFromProfile();
No; this is a free function that takes an image and a fiber trace.
PTR(SpectrumT) extractSum();
No; this should be an option to the free function.
const PTR(const FiberTraceFunction) getFiberTraceFunction() const;
void setFiberTraceFunction(PTR(FiberTraceFunction const) fiberTraceFunction);
Why do you need these? Maybe the getter
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?
const ndarray::Array< float, 1, 1 > getXCenters() const;
OK
void setXCenters( ndarray::Array< float, 1, 1 > const& xCenters );
Why do you need this?
ndarray::Array< float, 2, 1 > getXCentersMeas() const;
void setXCentersMeas( ndarray::Array< float, 2, 1 > const& );
How is this different from getXCenters?
PTR(Image) getReconstructed2DSpectrum(SpectrumT const& spectrum) const;
Should be a MaskedImage, I think. In CCD coords, of course.
PTR(Image) getReconstructedBackground(SpectrumT const& backgroundSpectrum) const;
Why do you need this too?
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
ndarray::Array<float, 1, 1> getTraceCoefficients() const;
Why do we need these? I can see a function that returns the model
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.