LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
Namespaces | Classes | Typedefs | Functions
lsst::meas::algorithms Namespace Reference

Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations. More...

Namespaces

 accumulator_mean_stack
 
 astrometrySourceSelector
 
 brightStarStamps
 
 coaddPsf
 
 convertRefcatManager
 
 convertReferenceCatalog
 
 detection
 
 dynamicDetection
 
 findCosmicRaysConfig
 
 flaggedSourceSelector
 
 gaussianPsfFactory
 
 htmIndexer
 
 indexerRegistry
 
 ingestIndexReferenceTask
 
 installGaussianPsf
 
 interp
 
 loadIndexedReferenceObjects
 
 loadReferenceObjects
 
 makeCoaddApCorrMap
 
 makePsfCandidates
 
 matcherSourceSelector
 
 measureApCorr
 
 objectSizeStarSelector
 
 pcaPsfDeterminer
 
 psfCandidate
 
 psfDeterminer
 
 psfSelectionFromMatchList
 
 readFitsCatalogTask
 
 readTextCatalogTask
 
 reserveSourcesTask
 
 simple_curve
 
 skyObjects
 
 sourceSelector
 
 stamps
 
 starSelector
 
 subtractBackground
 
 testUtils
 
 utils
 
 version
 

Classes

struct  CoaddBoundedFieldElement
 Struct used to hold one Exposure's data in a CoaddBoundedField. More...
 
class  CoaddBoundedField
 
class  CoaddPsfControl
 
class  CoaddPsf
 CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images. More...
 
class  DoubleGaussianPsf
 Represent a Psf as a circularly symmetrical double Gaussian. More...
 
class  ExposurePatch
 A convenience container for the exposure, peak and footprint that will be measured. More...
 
class  PsfImagePca
 
class  ImagePsf
 An intermediate base class for Psfs that use an image representation. More...
 
class  Defect
 Encapsulate information about a bad portion of a detector. More...
 
class  KernelPsf
 A Psf defined by a Kernel. More...
 
struct  KernelPsfPersistenceHelper
 A read-only singleton struct containing the schema and key used in persistence for KernelPsf. More...
 
class  KernelPsfFactory
 A PersistableFactory for KernelPsf and its subclasses. More...
 
class  PcaPsf
 Represent a PSF as a linear combination of PCA (== Karhunen-Loeve) basis functions. More...
 
class  PsfCandidate
 Class stored in SpatialCells for spatial Psf fitting. More...
 
class  ImagePsfTrampoline
 "Trampoline" for ImagePsf to let it be used as a base class in Python. More...
 
class  SingleGaussianPsf
 Represent a PSF as a circularly symmetrical Gaussian. More...
 
class  WarpedPsf
 A Psf class that maps an arbitrary Psf through a coordinate transformation. More...
 
class  evalChi2Visitor
 A class to pass around to all our PsfCandidates to evaluate the PSF fit's X^2. More...
 
class  MinimizeChi2
 

Typedefs

typedef std::vector< Defect::Ptr >::const_iterator DefectCIter
 

Functions

std::shared_ptr< afw::image::TransmissionCurve const > makeCoaddTransmissionCurve (std::shared_ptr< afw::geom::SkyWcs const > coaddWcs, afw::table::ExposureCatalog const &inputSensors)
 Create a TransmissionCurve that represents the effective throughput on a coadd. More...
 
template<typename MaskedImageT >
std::vector< std::shared_ptr< afw::detection::Footprint > > findCosmicRays (MaskedImageT &mimage, afw::detection::Psf const &psf, double const bkgd, daf::base::PropertySet const &ps, bool const keep)
 Find cosmic rays in an Image, and mask and remove them. More...
 
template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > makeExposurePatch (std::shared_ptr< ExposureT const > exp, std::shared_ptr< afw::detection::Footprint const > foot, geom::Point2D const &center)
 Factory function for ExposurePatch. More...
 
template<typename ExposureT >
std::shared_ptr< ExposurePatch< ExposureT > > makeExposurePatch (std::shared_ptr< ExposureT const > exp, afw::detection::Footprint const &standardFoot, geom::Point2D const &standardCenter, afw::geom::SkyWcs const &standardWcs)
 
template<typename MaskedImageT >
void interpolateOverDefects (MaskedImageT &mimage, afw::detection::Psf const &, std::vector< Defect::Ptr > &_badList, double fallbackValue, bool useFallbackValueAtEdge)
 Process a set of known bad pixels in an image. More...
 
template<typename PixelT >
std::shared_ptr< PsfCandidate< PixelT > > makePsfCandidate (std::shared_ptr< afw::table::SourceRecord > const &source, std::shared_ptr< afw::image::Exposure< PixelT >> image)
 Return a PsfCandidate of the right sort. More...
 
template<typename PixelT >
std::pair< std::shared_ptr< afw::math::LinearCombinationKernel >, std::vector< double > > createKernelFromPsfCandidates (afw::math::SpatialCellSet const &psfCells, lsst::geom::Extent2I const &dims, lsst::geom::Point2I const &xy0, int const nEigenComponents, int const spatialOrder, int const ksize, int const nStarPerCell, bool const constantWeight, int const border)
 Return a Kernel pointer and a list of eigenvalues resulting from analysing the provided SpatialCellSet. More...
 
template<typename PixelT >
int countPsfCandidates (afw::math::SpatialCellSet const &psfCells, int const nStarPerCell)
 Count the number of candidates in use. More...
 
template<typename PixelT >
std::pair< bool, double > fitSpatialKernelFromPsfCandidates (afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, int const nStarPerCell, double const tolerance, double const lambda)
 Fit spatial kernel using full-nonlinear optimization to estimate candidate amplitudes. More...
 
template<typename PixelT >
std::pair< bool, double > fitSpatialKernelFromPsfCandidates (afw::math::Kernel *kernel, afw::math::SpatialCellSet const &psfCells, bool const doNonLinearFit, int const nStarPerCell=-1, double const tolerance=1e-5, double const lambda=0.0)
 
template<typename ImageT >
double subtractPsf (afw::detection::Psf const &psf, ImageT *data, double x, double y, double psfFlux=std::numeric_limits< double >::quiet_NaN())
 
template<typename Image >
std::pair< std::vector< double >, afw::math::KernelListfitKernelParamsToImage (afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
 Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary. More...
 
template<typename Image >
std::pair< std::shared_ptr< afw::math::Kernel >, std::pair< double, double > > fitKernelToImage (afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
 Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary. More...
 
 PYBIND11_MODULE (psfCandidate, mod)
 
geom::Box2I getOverallBBox (std::vector< std::shared_ptr< afw::image::Image< double >>> const &imgVector)
 
void addToImage (std::shared_ptr< afw::image::Image< double >> image, std::vector< std::shared_ptr< afw::image::Image< double >>> const &imgVector, std::vector< double > const &weightVector)
 
void setSpatialParameters (afw::math::Kernel *kernel, std::vector< double > const &coeffs)
 Fit a Kernel's spatial variability from a set of stars. More...
 
void setSpatialParameters (afw::math::Kernel *kernel, Eigen::VectorXd const &vec)
 Fit a Kernel's spatial variability from a set of stars. More...
 
template<typename MaskedImageT >
double subtractPsf (afw::detection::Psf const &psf, MaskedImageT *data, double x, double y, double psfFlux)
 Subtract a PSF from an image at a given position. More...
 

Detailed Description

Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.

Typedef Documentation

◆ DefectCIter

Definition at line 49 of file Interp.cc.

Function Documentation

◆ addToImage()

void lsst::meas::algorithms::addToImage ( std::shared_ptr< afw::image::Image< double >>  image,
std::vector< std::shared_ptr< afw::image::Image< double >>> const &  imgVector,
std::vector< double > const &  weightVector 
)

Definition at line 201 of file CoaddPsf.cc.

203  {
204  assert(imgVector.size() == weightVector.size());
205  for (unsigned int i = 0; i < imgVector.size(); i++) {
206  std::shared_ptr<afw::image::Image<double>> componentImg = imgVector[i];
207  double weight = weightVector[i];
208  double sum = ndarray::asEigenMatrix(componentImg->getArray()).sum();
209 
210  // Now get the portion of the component image which is appropriate to add
211  // If the default image size is used, the component is guaranteed to fit,
212  // but not if a size has been specified.
213  geom::Box2I cBBox = componentImg->getBBox();
214  geom::Box2I overlap(cBBox);
215  overlap.clip(image->getBBox());
216  // JFB: A subimage view of the image we want to add to, containing only the overlap region.
217  afw::image::Image<double> targetSubImage(*image, overlap);
218  // JFB: A subimage view of the image we want to add from, containing only the overlap region.
219  afw::image::Image<double> cSubImage(*componentImg, overlap);
220  targetSubImage.scaledPlus(weight / sum, cSubImage);
221  }
222 }
An integer coordinate rectangle.
Definition: Box.h:55
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
T size(T... args)
afw::table::Key< double > weight

◆ countPsfCandidates()

template<typename PixelT >
int lsst::meas::algorithms::countPsfCandidates ( afw::math::SpatialCellSet const &  psfCells,
int const  nStarPerCell = -1 
)

Count the number of candidates in use.

Definition at line 337 of file SpatialModelPsf.cc.

337  {
338  countVisitor<PixelT> counter;
339  psfCells.visitCandidates(&counter, nStarPerCell);
340 
341  return counter.getN();
342 }

◆ createKernelFromPsfCandidates()

template<typename PixelT >
std::pair< std::shared_ptr< afw::math::LinearCombinationKernel >, std::vector< double > > lsst::meas::algorithms::createKernelFromPsfCandidates ( afw::math::SpatialCellSet const &  psfCells,
lsst::geom::Extent2I const &  dims,
lsst::geom::Point2I const &  xy0,
int const  nEigenComponents,
int const  spatialOrder,
int const  ksize,
int const  nStarPerCell,
bool const  constantWeight,
int const  border 
)

Return a Kernel pointer and a list of eigenvalues resulting from analysing the provided SpatialCellSet.

The Kernel is a LinearCombinationKernel of the first nEigenComponents eigenImages

N.b. This is templated over the Pixel type of the science image

Parameters
psfCellsA SpatialCellSet containing PsfCandidates
dimsDimensions of image
xy0Origin of image
nEigenComponentsnumber of eigen components to keep; <= 0 => infty
spatialOrderOrder of spatial variation (cf. afw::math::PolynomialFunction2)
ksizeSize of generated Kernel images
nStarPerCellmax no. of stars per cell; <= 0 => infty
constantWeightshould each star have equal weight in the fit?
borderBorder size for background subtraction

Definition at line 193 of file SpatialModelPsf.cc.

203  {
204  typedef typename afw::image::Image<PixelT> ImageT;
205  typedef typename afw::image::MaskedImage<PixelT> MaskedImageT;
206 
207  //
208  // Set the sizes for PsfCandidates made from either Images or MaskedImages
209  //
210  // lsst::meas::algorithms::PsfCandidate<ImageT>::setWidth(ksize);
211  // lsst::meas::algorithms::PsfCandidate<ImageT>::setHeight(ksize);
212  // lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setWidth(ksize);
213  // lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setHeight(ksize);
216 
217  PsfImagePca<MaskedImageT> imagePca(constantWeight, border); // Here's the set of images we'll analyze
218 
219  {
220  SetPcaImageVisitor<PixelT> importStarVisitor(&imagePca);
221  bool const ignoreExceptions = true;
222  psfCells.visitCandidates(&importStarVisitor, nStarPerCell, ignoreExceptions);
223  }
224 
225  //
226  // Do a PCA decomposition of those PSF candidates.
227  //
228  // We have "gappy" data; in other words we don't want to include any pixels with INTRP set
229  //
230  int niter = 10; // number of iterations of updateBadPixels
231  double deltaLim = 10.0; // acceptable value of delta, the max change due to updateBadPixels
232  lsst::afw::image::MaskPixel const BAD = afw::image::Mask<>::getPlaneBitMask("BAD");
233  lsst::afw::image::MaskPixel const CR = afw::image::Mask<>::getPlaneBitMask("CR");
234  lsst::afw::image::MaskPixel const INTRP = afw::image::Mask<>::getPlaneBitMask("INTRP");
235 
236  for (int i = 0; i != niter; ++i) {
237  int const ncomp =
238  (i == 0) ? 0
239  : ((nEigenComponents == 0) ? imagePca.getEigenImages().size() : nEigenComponents);
240  double delta = imagePca.updateBadPixels(BAD | CR | INTRP, ncomp);
241  if (i > 0 && delta < deltaLim) {
242  break;
243  }
244 
245  imagePca.analyze();
246  }
247 
248  std::vector<std::shared_ptr<MaskedImageT>> eigenImages = imagePca.getEigenImages();
249  std::vector<double> eigenValues = imagePca.getEigenValues();
250  int const nEigen = static_cast<int>(eigenValues.size());
251 
252  int const ncomp = (nEigenComponents <= 0 || nEigen < nEigenComponents) ? nEigen : nEigenComponents;
253  //
254  // Set the background level of the components to 0.0 to avoid coupling variable background
255  // levels to the form of the Psf. More precisely, we calculate the mean of an outer "annulus"
256  // of width bkg_border
257  //
258  for (int k = 0; k != ncomp; ++k) {
259  ImageT const& im = *eigenImages[k]->getImage();
260 
261  int bkg_border = 2;
262  if (bkg_border > im.getWidth()) {
263  bkg_border = im.getWidth() / 2;
264  }
265  if (bkg_border > im.getHeight()) {
266  bkg_border = im.getHeight() / 2;
267  }
268 
269  double sum = 0;
270  // Bottom and Top borders
271  for (int i = 0; i != bkg_border; ++i) {
272  typename ImageT::const_x_iterator ptrB = im.row_begin(i),
273  ptrT = im.row_begin(im.getHeight() - i - 1);
274  for (int j = 0; j != im.getWidth(); ++j, ++ptrB, ++ptrT) {
275  sum += *ptrB + *ptrT;
276  }
277  }
278  for (int i = bkg_border; i < im.getHeight() - bkg_border; ++i) {
279  // Left and Right borders
280  typename ImageT::const_x_iterator ptrL = im.row_begin(i),
281  ptrR = im.row_begin(i) + im.getWidth() - bkg_border;
282  for (int j = 0; j != bkg_border; ++j, ++ptrL, ++ptrR) {
283  sum += *ptrL + *ptrR;
284  }
285  }
286  sum /= 2 * (bkg_border * im.getWidth() + bkg_border * (im.getHeight() - 2 * bkg_border));
287 
288  *eigenImages[k] -= sum;
289  }
290  //
291  // Now build our LinearCombinationKernel; build the lists of basis functions
292  // and spatial variation, then assemble the Kernel
293  //
294  afw::math::KernelList kernelList;
296  geom::Box2D const range = geom::Box2D(geom::Point2D(xy0), geom::Extent2D(dims));
297 
298  for (int i = 0; i != ncomp; ++i) {
299  {
300  // Enforce unit sum for kernel by construction
301  // Zeroth component has unit sum
302  // Other components have zero sum by normalising and then subtracting the zeroth component
303  ImageT& image = *eigenImages[i]->getImage();
304  double sum = std::accumulate(image.begin(true), image.end(true), 0.0);
305  if (i == 0) {
306  image /= sum;
307  } else {
308  for (typename ImageT::fast_iterator ptr0 = eigenImages[0]->getImage()->begin(true),
309  ptr1 = image.begin(true), end = image.end(true);
310  ptr1 != end; ++ptr0, ++ptr1) {
311  *ptr1 = *ptr1 / sum - *ptr0;
312  }
313  }
314  }
315 
316  kernelList.push_back(std::shared_ptr<afw::math::Kernel>(new afw::math::FixedKernel(
317  afw::image::Image<afw::math::Kernel::Pixel>(*eigenImages[i]->getImage(), true))));
318 
319  afw::math::Kernel::SpatialFunctionPtr
320  // spatialFunction(new afw::math::PolynomialFunction2<double>(spatialOrder));
321  spatialFunction(new afw::math::Chebyshev1Function2<double>(spatialOrder, range));
322  spatialFunction->setParameter(0, 1.0); // the constant term; all others are 0
323  spatialFunctionList.push_back(spatialFunction);
324  }
325 
327  new afw::math::LinearCombinationKernel(kernelList, spatialFunctionList));
328 
329  return std::make_pair(psf, eigenValues);
330 }
int end
T accumulate(T... args)
T begin(T... args)
static void setWidth(int width)
Set the width of the image that getImage should return.
Definition: SpatialCell.h:136
static void setHeight(int height)
Set the height of the image that getImage should return.
Definition: SpatialCell.h:141
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
T make_pair(T... args)
std::vector< std::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:462
T push_back(T... args)
Key< int > psf
Definition: Exposure.cc:65

◆ findCosmicRays()

template<typename MaskedImageT >
std::vector< std::shared_ptr< afw::detection::Footprint > > lsst::meas::algorithms::findCosmicRays ( MaskedImageT &  mimage,
afw::detection::Psf const &  psf,
double const  bkgd,
daf::base::PropertySet const &  ps,
bool const  keep 
)

Find cosmic rays in an Image, and mask and remove them.

Returns
vector of CR's Footprints

In this loop, we look for strings of CRpixels on the same row and adjoining columns; each of these becomes a Span with a unique ID.

In this loop, we look for strings of CRpixels on the same row and adjoining columns; each of these becomes a Span with a unique ID.

Parameters
mimageImage to search
psfthe Image's PSF
bkgdunsubtracted background of frame, DN
psPropertySet directing the behavior
keepif true, don't remove the CRs

Definition at line 315 of file CR.cc.

321  {
322  typedef typename MaskedImageT::Image ImageT;
323  typedef typename ImageT::Pixel ImagePixel;
324  typedef typename MaskedImageT::Mask::Pixel MaskPixel;
325 
326  // Parse the PropertySet
327  double const minSigma = ps.getAsDouble("minSigma"); // min sigma over sky in pixel for CR candidate
328  double const minDn = ps.getAsDouble("min_DN"); // min number of DN in an CRs
329  double const cond3Fac = ps.getAsDouble("cond3_fac"); // fiddle factor for condition #3
330  double const cond3Fac2 = ps.getAsDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
331  int const niteration = ps.getAsInt("niteration"); // Number of times to look for contaminated
332  // pixels near CRs
333  int const nCrPixelMax = ps.getAsInt("nCrPixelMax"); // maximum number of contaminated pixels
334  /*
335  * thresholds for 3rd condition
336  *
337  * Realise PSF at center of image
338  */
339  std::shared_ptr<afw::math::Kernel const> kernel = psf.getLocalKernel();
340  if (!kernel) {
341  throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
342  }
343  afw::detection::Psf::Image psfImage =
344  afw::detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
345  kernel->computeImage(psfImage, true);
346 
347  int const xc = kernel->getCtr().getX(); // center of PSF
348  int const yc = kernel->getCtr().getY();
349 
350  double const I0 = psfImage(xc, yc);
351  double const thresH =
352  cond3Fac2 * (0.5 * (psfImage(xc - 1, yc) + psfImage(xc + 1, yc))) / I0; // horizontal
353  double const thresV = cond3Fac2 * (0.5 * (psfImage(xc, yc - 1) + psfImage(xc, yc + 1))) / I0; // vertical
354  double const thresD = cond3Fac2 *
355  (0.25 * (psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
356  psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1))) /
357  I0; // diag
358  /*
359  * Setup desired mask planes
360  */
361  MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
362  MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
363  MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
364  MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels
365  MaskPixel const nodataBit = mimage.getMask()->getPlaneBitMask("NO_DATA"); // Non data pixels
366 
367  MaskPixel const badMask = (badBit | interpBit | saturBit | nodataBit); // naughty pixels
368  /*
369  * Go through the frame looking at each pixel (except the edge ones which we ignore)
370  */
371  int const ncol = mimage.getWidth();
372  int const nrow = mimage.getHeight();
373 
374  std::vector<CRPixel<ImagePixel>> crpixels; // storage for detected CR-contaminated pixels
375  typedef typename std::vector<CRPixel<ImagePixel>>::iterator crpixel_iter;
376  typedef typename std::vector<CRPixel<ImagePixel>>::reverse_iterator crpixel_riter;
377 
378  for (int j = 1; j < nrow - 1; ++j) {
379  typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data
380 
381  for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
382  ImagePixel corr = 0;
383  if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma, thresH, thresV, thresD, bkgd, cond3Fac)) {
384  continue;
385  }
386  /*
387  * condition #4
388  */
389  if (loc.mask() & badMask) {
390  continue;
391  }
392  if ((loc.mask(-1, 1) | loc.mask(0, 1) | loc.mask(1, 1) | loc.mask(-1, 0) | loc.mask(1, 0) |
393  loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) &
394  interpBit) {
395  continue;
396  }
397  /*
398  * OK, it's a CR
399  *
400  * replace CR-contaminated pixels with reasonable values as we go through
401  * image, which increases the detection rate
402  */
403  crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
404  loc.image() = corr; /* just a preliminary estimate */
405 
406  if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
407  reinstateCrPixels(mimage.getImage().get(), crpixels);
408 
409  throw LSST_EXCEPT(pex::exceptions::LengthError,
410  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
411  }
412  }
413  }
414  /*
415  * We've found them on a pixel-by-pixel basis, now merge those pixels
416  * into cosmic rays
417  */
418  std::vector<int> aliases; // aliases for initially disjoint parts of CRs
419  aliases.reserve(1 + crpixels.size() / 2); // initial size of aliases
420 
421  std::vector<afw::detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
422  spans.reserve(aliases.capacity()); // initial size of spans
423 
424  aliases.push_back(0); // 0 --> 0
425 
431  int ncr = 0; // number of detected cosmic rays
432  if (!crpixels.empty()) {
433  int id; // id number for a CR
434  int x0 = -1, x1 = -1, y = -1; // the beginning and end column, and row of this span in a CR
435 
436  // I am dummy
437  CRPixel<ImagePixel> dummy(0, -1, 0, -1);
438  crpixels.push_back(dummy);
439  // printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id,
440  // dummy.col, dummy.row, (double)dummy.val);
441  for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1; ++crp) {
442  // printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col,
443  // crp->row, (double)crp->val);
444 
445  if (crp->id < 0) { // not already assigned
446  crp->id = ++ncr; // a new CR
447  aliases.push_back(crp->id);
448  y = crp->row;
449  x0 = x1 = crp->col;
450  // printf(" Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
451  }
452  id = crp->id;
453  // printf(" Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row,
454  // crp[1].col);
455 
456  if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
457  // printf(" Adjoining! Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
458  crp[1].id = id;
459  ++x1;
460  } else {
461  assert(y >= 0 && x0 >= 0 && x1 >= 0);
462  spans.push_back(afw::detection::IdSpan::Ptr(new afw::detection::IdSpan(id, y, x0, x1)));
463  // printf(" Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
464  }
465  }
466  }
467 
468  // At the end of this loop, all crpixel entries have been assigned an ID,
469  // except for the "dummy" entry at the end of the array.
470  if (crpixels.size() > 0) {
471  for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
472  assert(cp->id >= 0);
473  assert(cp->col >= 0);
474  assert(cp->row >= 0);
475  }
476  // dummy:
477  assert(crpixels[crpixels.size() - 1].id == -1);
478  assert(crpixels[crpixels.size() - 1].col == 0);
479  assert(crpixels[crpixels.size() - 1].row == -1);
480  }
481 
482  for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
483  sp++) {
484  assert((*sp)->id >= 0);
485  assert((*sp)->y >= 0);
486  assert((*sp)->x0 >= 0);
487  assert((*sp)->x1 >= (*sp)->x0);
488  for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
489  assert((*sp2)->y >= (*sp)->y);
490  if ((*sp2)->y == (*sp)->y) {
491  assert((*sp2)->x0 > (*sp)->x1);
492  }
493  }
494  }
495 
496  /*
497  * See if spans touch each other
498  */
499  for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end;
500  ++sp) {
501  int const y = (*sp)->y;
502  int const x0 = (*sp)->x0;
503  int const x1 = (*sp)->x1;
504 
505  // this loop will probably run for only a few steps
506  for (std::vector<afw::detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
507  if ((*sp2)->y == y) {
508  // on this row (but not adjoining columns, since it would have been merged into this span);
509  // keep looking.
510  continue;
511  } else if ((*sp2)->y != (y + 1)) {
512  // sp2 is more than one row below; can't be connected.
513  break;
514  } else if ((*sp2)->x0 > (x1 + 1)) {
515  // sp2 is more than one column away to the right; can't be connected
516  break;
517  } else if ((*sp2)->x1 >= (x0 - 1)) {
518  // touches
519  int r1 = afw::detection::resolve_alias(aliases, (*sp)->id);
520  int r2 = afw::detection::resolve_alias(aliases, (*sp2)->id);
521  aliases[r1] = r2;
522  }
523  }
524  }
525 
526  /*
527  * Resolve aliases; first alias chains, then the IDs in the spans
528  */
529  for (unsigned int i = 0; i != spans.size(); ++i) {
530  spans[i]->id = afw::detection::resolve_alias(aliases, spans[i]->id);
531  }
532 
533  /*
534  * Sort spans by ID, so we can sweep through them once
535  */
536  if (spans.size() > 0) {
537  std::sort(spans.begin(), spans.end(), afw::detection::IdSpanCompar());
538  }
539 
540  /*
541  * Build Footprints from spans
542  */
544 
545  if (spans.size() > 0) {
546  int id = spans[0]->id;
547  unsigned int i0 = 0; // initial value of i
548  for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
549  if (i == spans.size() || spans[i]->id != id) {
551 
553  spanList.reserve(i - i0);
554  for (; i0 < i; ++i0) {
555  spanList.push_back(afw::geom::Span(spans[i0]->y, spans[i0]->x0, spans[i0]->x1));
556  }
557  cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(spanList)));
558  CRs.push_back(cr);
559  }
560 
561  if (i < spans.size()) {
562  id = spans[i]->id;
563  }
564  }
565  }
566 
567  reinstateCrPixels(mimage.getImage().get(), crpixels);
568  /*
569  * apply condition #1
570  */
571  CountsInCR<typename ImageT::Pixel> CountDN(bkgd);
572  for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator cr = CRs.begin(), end = CRs.end();
573  cr != end; ++cr) {
574  // find the sum of pixel values within the CR
575  (*cr)->getSpans()->applyFunctor(CountDN, *mimage.getImage());
576 
577  LOGL_DEBUG("TRACE4.lsst.algorithms.CR", "CR at (%d, %d) has %g DN", (*cr)->getBBox().getMinX(),
578  (*cr)->getBBox().getMinY(), CountDN.getCounts());
579  if (CountDN.getCounts() < minDn) { /* not bright enough */
580  LOGL_DEBUG("TRACE5.lsst.algorithms.CR", "Erasing CR");
581 
582  cr = CRs.erase(cr);
583  --cr; // back up to previous CR (we're going to increment it)
584  --end;
585  }
586  CountDN.reset();
587  }
588  ncr = CRs.size(); /* some may have been too faint */
589  /*
590  * We've found them all, time to kill them all
591  */
592  bool const debias_values = true;
593  bool grow = false;
594  LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing initial list of CRs");
595  removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
596 #if 0 // Useful to see phase 2 in display; debugging only
597  (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
598  mimage.getMask()->getPlaneBitMask("DETECTED"));
599 #endif
600  /*
601  * Now that we've removed them, go through image again, examining area around
602  * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
603  *
604  * We iterate niteration times; niter==1 was sufficient for SDSS data, but megacam
605  * CCDs are different -- who knows for other devices?
606  */
607  bool too_many_crs = false; // we've seen too many CR pixels
608  int nextra = 0; // number of pixels added to list of CRs
609  for (int i = 0; i != niteration && !too_many_crs; ++i) {
610  LOGL_DEBUG("TRACE1.lsst.algorithms.CR", "Starting iteration %d", i);
611  for (std::vector<std::shared_ptr<afw::detection::Footprint>>::iterator fiter = CRs.begin();
612  fiter != CRs.end(); fiter++) {
614  /*
615  * Are all those `CR' pixels interpolated? If so, don't grow it
616  */
617  {
618  // this work should be taken on in DM-9538
619  // std::shared_ptr<afw::detection::Footprint> om = footprintAndMask(cr, mimage.getMask(),
620  // interpBit);
621  // Placeholder until DM-9538 then remove empty footprint
622  auto om = std::make_shared<afw::detection::Footprint>();
623  int const npix = (om) ? om->getArea() : 0;
624 
625  if (static_cast<std::size_t>(npix) == cr->getArea()) {
626  continue;
627  }
628  }
629  /*
630  * No; some of the suspect pixels aren't interpolated
631  */
632  afw::detection::Footprint extra; // extra pixels added to cr
633  for (auto siter = cr->getSpans()->begin(); siter != cr->getSpans()->end(); siter++) {
634  auto const span = siter;
635 
636  /*
637  * Check the lines above and below the span. We're going to check a 3x3 region around
638  * the pixels, so we need a buffer around the edge. We check the pixels just to the
639  * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
640  * column direction, but only 1 in the row direction.
641  */
642  int const y = span->getY() - mimage.getY0();
643  if (y < 2 || y >= nrow - 2) {
644  continue;
645  }
646  int x0 = span->getX0() - mimage.getX0();
647  int x1 = span->getX1() - mimage.getX0();
648  x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
649  x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;
650 
651  checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
652  bkgd, 0, keep);
653  checkSpanForCRs(&extra, crpixels, y, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
654  bkgd, 0, keep);
655  checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage, minSigma / 2, thresH, thresV, thresD,
656  bkgd, 0, keep);
657  }
658 
659  if (extra.getSpans()->size() > 0) { // we added some pixels
660  if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
661  too_many_crs = true;
662  break;
663  }
664 
665  nextra += extra.getArea();
666 
667  std::vector<afw::geom::Span> tmpSpanList(cr->getSpans()->begin(), cr->getSpans()->end());
668  for (auto const &spn : (*extra.getSpans())) {
669  tmpSpanList.push_back(spn);
670  }
671  cr->setSpans(std::make_shared<afw::geom::SpanSet>(std::move(tmpSpanList)));
672  }
673  }
674 
675  if (nextra == 0) {
676  break;
677  }
678  }
679  /*
680  * mark those pixels as CRs
681  */
682  if (!too_many_crs) {
683  for (auto const &foot : CRs) {
684  foot->getSpans()->setMask(*mimage.getMask().get(), crBit);
685  }
686  }
687  /*
688  * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
689  * first value stored (hence the uses of rbegin/rend)
690  *
691  * We have to do this if we decide _not_ to remove certain CRs,
692  * for example those which lie next to saturated pixels
693  */
694  if (keep || too_many_crs) {
695  if (crpixels.size() > 0) {
696  int const imageX0 = mimage.getX0();
697  int const imageY0 = mimage.getY0();
698 
699  std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
700 
701  crpixel_riter rend = crpixels.rend();
702  for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
703  if (crp->row == -1)
704  // dummy; skip it.
705  continue;
706  mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
707  }
708  }
709  } else {
710  if (true || nextra > 0) {
711  grow = true;
712  LOGL_DEBUG("TRACE2.lsst.algorithms.CR", "Removing final list of CRs, grow = %d", grow);
713  removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
714  }
715  /*
716  * we interpolated over all CR pixels, so set the interp bits too
717  */
718  for (auto const &foot : CRs) {
719  foot->getSpans()->setMask(*mimage.getMask().get(), static_cast<MaskPixel>(crBit | interpBit));
720  }
721  }
722 
723  if (too_many_crs) { // we've cleaned up, so we can throw the exception
724  throw LSST_EXCEPT(pex::exceptions::LengthError,
725  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
726  }
727 
728  return CRs;
729 }
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
#define LOGL_DEBUG(logger, message...)
Log a debug-level message using a varargs/printf style interface.
Definition: Log.h:515
int y
Definition: SpanSet.cc:48
T capacity(T... args)
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T empty(T... args)
T end(T... args)
T erase(T... args)
T move(T... args)
int resolve_alias(const std::vector< int > &aliases, int id)
Follow a chain of aliases, returning the final resolved value.
Definition: CR.cc:96
std::int32_t MaskPixel
default type for Masks and MaskedImage Masks
lsst::afw::detection::Footprint Footprint
Definition: Source.h:57
float Pixel
Typedefs to be used for pixel values.
Definition: common.h:37
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
T rbegin(T... args)
T rend(T... args)
T reserve(T... args)
T sort(T... args)
int id
Definition: CR.cc:143
int row
Definition: CR.cc:145
int col
Definition: CR.cc:144

◆ fitKernelParamsToImage()

template<typename Image >
std::pair< std::vector< double >, afw::math::KernelList > lsst::meas::algorithms::fitKernelParamsToImage ( afw::math::LinearCombinationKernel const &  kernel,
Image const &  image,
geom::Point2D const &  pos 
)

Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

Returns
std::pair(coefficients, std::pair(kernels, center amplitude))
Parameters
kernelthe Kernel to fit
imagethe image to be fit
posthe position of the object

Definition at line 1037 of file SpatialModelPsf.cc.

1041  {
1042  typedef afw::image::Image<afw::math::Kernel::Pixel> KernelT;
1043 
1044  afw::math::KernelList kernels = kernel.getKernelList(); // the Kernels that kernel adds together
1045  int const nKernel = kernels.size();
1046 
1047  if (nKernel == 0) {
1048  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Your kernel must have at least one component");
1049  }
1050 
1051  /*
1052  * Go through all the kernels, get a copy centered at the desired sub-pixel position, and then
1053  * extract a subImage from the parent image at the same place
1054  */
1055  std::vector<std::shared_ptr<KernelT>> kernelImages = offsetKernel<KernelT>(kernel, pos[0], pos[1]);
1056  geom::BoxI bbox(kernelImages[0]->getBBox());
1057  Image const& subImage(Image(image, bbox, afw::image::PARENT, false)); // shallow copy
1058 
1059  /*
1060  * Solve the linear problem subImage = sum x_i K_i + epsilon; we solve this for x_i by constructing the
1061  * normal equations, A x = b
1062  */
1063  Eigen::MatrixXd A(nKernel, nKernel);
1064  Eigen::VectorXd b(nKernel);
1065 
1066  for (int i = 0; i != nKernel; ++i) {
1067  b(i) = afw::image::innerProduct(*kernelImages[i], *subImage.getImage());
1068 
1069  for (int j = i; j != nKernel; ++j) {
1070  A(i, j) = A(j, i) = afw::image::innerProduct(*kernelImages[i], *kernelImages[j]);
1071  }
1072  }
1073  Eigen::VectorXd x(nKernel);
1074 
1075  if (nKernel == 1) {
1076  x(0) = b(0) / A(0, 0);
1077  } else {
1078  x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
1079  }
1080 
1081  // the XY0() point of the shifted Kernel basis functions
1082  geom::Point2I const xy0 = kernelImages[0]->getXY0();
1083 
1084  afw::math::KernelList newKernels(nKernel);
1085  std::vector<double> params(nKernel);
1086  for (int i = 0; i != nKernel; ++i) {
1087  std::shared_ptr<afw::math::Kernel> newKernel(new afw::math::FixedKernel(*kernelImages[i]));
1088  newKernel->setCtr(xy0 + newKernel->getDimensions() / 2);
1089 
1090  params[i] = x[i];
1091  newKernels[i] = newKernel;
1092  }
1093 
1094  return std::make_pair(params, newKernels);
1095 }
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
double x
table::Key< int > b
Reports attempts to exceed implementation-defined length limits for some classes.
Definition: Runtime.h:76
double innerProduct(Image1T const &lhs, Image2T const &rhs, int const border=0)
Calculate the inner product of two images.
Definition: ImagePca.cc:414

◆ fitKernelToImage()

template<typename Image >
std::pair< std::shared_ptr< afw::math::Kernel >, std::pair< double, double > > lsst::meas::algorithms::fitKernelToImage ( afw::math::LinearCombinationKernel const &  kernel,
Image const &  image,
geom::Point2D const &  pos 
)

Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

Returns
std::pair(best-fit kernel, std::pair(amp, chi^2))
Parameters
kernelthe Kernel to fit
imagethe image to be fit
posthe position of the object

Definition at line 1104 of file SpatialModelPsf.cc.

1108  {
1110  fitKernelParamsToImage(kernel, image, pos);
1111  std::vector<double> params = fit.first;
1112  afw::math::KernelList kernels = fit.second;
1113  int const nKernel = params.size();
1114  assert(kernels.size() == static_cast<unsigned int>(nKernel));
1115 
1116  double amp = 0.0;
1117  for (int i = 0; i != nKernel; ++i) {
1118  std::shared_ptr<afw::math::Kernel> base = kernels[i];
1119  std::shared_ptr<afw::math::FixedKernel> k = std::static_pointer_cast<afw::math::FixedKernel>(base);
1120  amp += params[i] * k->getSum();
1121  }
1122 
1123  std::shared_ptr<afw::math::Kernel> outputKernel(new afw::math::LinearCombinationKernel(kernels, params));
1124  double chisq = 0.0;
1125  outputKernel->setCtr(kernels[0]->getCtr());
1126 
1127  return std::make_pair(outputKernel, std::make_pair(amp, chisq));
1128 }
std::pair< std::vector< double >, afw::math::KernelList > fitKernelParamsToImage(afw::math::LinearCombinationKernel const &kernel, Image const &image, geom::Point2D const &pos)
Fit a LinearCombinationKernel to an Image, allowing the coefficients of the components to vary.

◆ fitSpatialKernelFromPsfCandidates() [1/2]

template<typename PixelT >
std::pair< bool, double > lsst::meas::algorithms::fitSpatialKernelFromPsfCandidates ( afw::math::Kernel kernel,
afw::math::SpatialCellSet const &  psfCells,
int const  nStarPerCell = -1,
double const  tolerance = 1e-5,
double const  lambda = 0.0 
)

Fit spatial kernel using full-nonlinear optimization to estimate candidate amplitudes.

Parameters
kernelthe Kernel to fit
psfCellsA SpatialCellSet containing PsfCandidates
nStarPerCellmax no. of stars per cell; <= 0 => infty
toleranceTolerance; how close chi^2 should be to true minimum
lambdafloor for variance is lambda*data

Definition at line 592 of file SpatialModelPsf.cc.

598  {
599  int const nComponents = kernel->getNKernelParameters();
600  int const nSpatialParams = kernel->getNSpatialParameters();
601  //
602  // visitor that evaluates the chi^2 of the current fit
603  //
604  evalChi2Visitor<PixelT> getChi2(*kernel, lambda);
605  //
606  // We have to unpack the Kernel coefficients into a linear array, coeffs
607  //
608  std::vector<double> coeffs; // The coefficients we want to fit
609  coeffs.assign(nComponents * nSpatialParams, 0.0);
610 
611  std::vector<double> stepSize; // step sizes
612  stepSize.assign(nComponents * nSpatialParams, 100);
613  //
614  // Translate that into minuit's language
615  //
616  ROOT::Minuit2::MnUserParameters fitPar;
617  std::vector<std::string> paramNames;
618  paramNames.reserve(nComponents * nSpatialParams);
619 
620  for (int i = 0, c = 0; c != nComponents; ++c) {
621  coeffs[i] = 1; // the constant part of each spatial order
622  for (int s = 0; s != nSpatialParams; ++s, ++i) {
623  paramNames.push_back((boost::format("C%d:%d") % c % s).str());
624  fitPar.Add(paramNames[i].c_str(), coeffs[i], stepSize[i]);
625  }
626  }
627  fitPar.Fix("C0:0");
628  //
629  // Create the minuit object that knows how to minimise our functor
630  //
631  MinimizeChi2<PixelT> minimizerFunc(getChi2, kernel, psfCells, nStarPerCell, nComponents, nSpatialParams);
632 
633  double const errorDef = 1.0; // use +- 1sigma errors
634  minimizerFunc.setErrorDef(errorDef);
635  //
636  // tell minuit about it
637  //
638  ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
639  //
640  // And let it loose
641  //
642  int maxFnCalls = 0; // i.e. unlimited
643  ROOT::Minuit2::FunctionMinimum min =
644  migrad(maxFnCalls, tolerance / (1e-4 * errorDef)); // minuit uses 0.1*1e-3*tolerance*errorDef
645 
646  float minChi2 = min.Fval();
647  bool const isValid = min.IsValid() && std::isfinite(minChi2);
648 
649  if (true || isValid) { // calculate coeffs even in minuit is unhappy
650  for (int i = 0; i != nComponents * nSpatialParams; ++i) {
651  coeffs[i] = min.UserState().Value(i);
652  }
653 
654  setSpatialParameters(kernel, coeffs);
655  }
656 
657 #if 0 // Estimate errors; we don't really need this
658  ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
659  for (int i = 0, c = 0; c != nComponents; ++c) {
660  for (int s = 0; s != nSpatialParams; ++s, ++i) {
661  char const *name = paramNames[i].c_str();
662  printf("%s %g", name, min.UserState().Value(name));
663  if (isValid && !fitPar.Parameter(fitPar.Index(name)).IsFixed()) {
664  printf(" (%g+%g)\n", minos(i).first, minos(i).second);
665  }
666  printf("\n");
667  }
668  }
669 #endif
670  //
671  // One time more through the Candidates setting their chi^2 values. We'll
672  // do all the candidates this time, not just the first nStarPerCell
673  //
674  psfCells.visitAllCandidates(&getChi2, true);
675 
676  return std::make_pair(isValid, minChi2);
677 }
table::Key< std::string > name
Definition: Amplifier.cc:116
int min
T assign(T... args)
bool isValid
Definition: fits.cc:399
T printf(T... args)
T isfinite(T... args)
void setSpatialParameters(afw::math::Kernel *kernel, Eigen::VectorXd const &vec)
Fit a Kernel's spatial variability from a set of stars.

◆ fitSpatialKernelFromPsfCandidates() [2/2]

template<typename PixelT >
std::pair< bool, double > lsst::meas::algorithms::fitSpatialKernelFromPsfCandidates ( afw::math::Kernel kernel,
afw::math::SpatialCellSet const &  psfCells,
bool const  doNonLinearFit,
int const  nStarPerCell = -1,
double const  tolerance = 1e-5,
double const  lambda = 0.0 
)
Parameters
kernelthe Kernel to fit
psfCellsA SpatialCellSet containing PsfCandidates
doNonLinearFitUse the full-up nonlinear fitter
nStarPerCellmax no. of stars per cell; <= 0 => infty
toleranceTolerance; how close chi^2 should be to true minimum
lambdafloor for variance is lambda*data

Definition at line 876 of file SpatialModelPsf.cc.

883  {
884  if (doNonLinearFit) {
885  return fitSpatialKernelFromPsfCandidates<PixelT>(kernel, psfCells, nStarPerCell, tolerance);
886  }
887 
888  double const tau = 0; // softening for errors
889 
890  afw::math::LinearCombinationKernel const* lcKernel =
891  dynamic_cast<afw::math::LinearCombinationKernel const*>(kernel);
892  if (!lcKernel) {
893  throw LSST_EXCEPT(
895  "Failed to cast Kernel to LinearCombinationKernel while building spatial PSF model");
896  }
897 #if 1
898  //
899  // Set the initial amplitudes of all our candidates
900  //
901  setAmplitudeVisitor<PixelT> setAmplitude;
902  psfCells.visitAllCandidates(&setAmplitude, true);
903 #endif
904  //
905  // visitor that fills out the A and b matrices (we'll solve A x = b for the coeffs, x)
906  //
907  FillABVisitor<PixelT> getAB(*lcKernel, tau);
908  //
909  // Actually visit all our candidates
910  //
911  psfCells.visitCandidates(&getAB, nStarPerCell, true);
912  //
913  // Extract A and b, and solve Ax = b
914  //
915  Eigen::MatrixXd const& A = getAB.getA();
916  Eigen::VectorXd const& b = getAB.getB();
917  Eigen::VectorXd x0(b.size()); // Solution to matrix problem
918 
919  switch (b.size()) {
920  case 0: // One candidate, no spatial variability
921  break;
922  case 1: // eigen can't/won't handle 1x1 matrices
923  x0(0) = b(0) / A(0, 0);
924  break;
925  default:
926  x0 = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
927  break;
928  }
929 #if 0
930  std::cout << "A " << A << std::endl;
931  std::cout << "b " << b.transpose() << std::endl;
932  std::cout << "x " << x.transpose() << std::endl;
933 
934  afw::image::Image<double> img(b.size(), b.size());
935  for (int j = 0; j < b.size(); ++j) {
936  for (int i = 0; i < b.size(); ++i) {
937  img(i, j) = A(i, j);
938  }
939  }
940  img.writeFits("a.fits");
941 
942  if (x.cols() >= 6) {
943  for (int i = 0; i != 6; ++i) {
944  double xcen = 25; double ycen = 35 + 35*i;
945  std::cout << "x, y " << xcen << " , " << ycen << " b "
946  << (x[3] + xcen*x[4] + ycen*x[5])/(x[0] + xcen*x[1] + ycen*x[2]) << std::endl;
947  }
948  }
949 #endif
950 
951  // Generate kernel parameters (including 0th component) from matrix solution
952  Eigen::VectorXd x(kernel->getNKernelParameters() * kernel->getNSpatialParameters()); // Kernel parameters
953  x(0) = 1.0;
954  std::fill(x.data() + 1, x.data() + kernel->getNSpatialParameters(), 0.0);
955  std::copy(x0.data(), x0.data() + x0.size(), x.data() + kernel->getNSpatialParameters());
956 
957  setSpatialParameters(kernel, x);
958  //
959  // One time more through the Candidates setting their chi^2 values. We'll
960  // do all the candidates this time, not just the first nStarPerCell
961  //
962  // visitor that evaluates the chi^2 of the current fit
963  //
964  evalChi2Visitor<PixelT> getChi2(*kernel, lambda);
965 
966  psfCells.visitAllCandidates(&getChi2, true);
967 
968  return std::make_pair(true, getChi2.getValue());
969 }
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
T copy(T... args)
T endl(T... args)
T fill(T... args)

◆ getOverallBBox()

geom::Box2I lsst::meas::algorithms::getOverallBBox ( std::vector< std::shared_ptr< afw::image::Image< double >>> const &  imgVector)

Definition at line 188 of file CoaddPsf.cc.

188  {
190  // Calculate the box which will contain them all
191  for (unsigned int i = 0; i < imgVector.size(); i++) {
192  std::shared_ptr<afw::image::Image<double>> componentImg = imgVector[i];
193  geom::Box2I cBBox = componentImg->getBBox();
194  bbox.include(cBBox); // JFB: this works even on empty bboxes
195  }
196  return bbox;
197 }

◆ interpolateOverDefects()

template<typename MaskedImageT >
void lsst::meas::algorithms::interpolateOverDefects ( MaskedImageT &  mimage,
afw::detection::Psf const &  psf,
std::vector< Defect::Ptr > &  _badList,
double  fallbackValue = 0.0,
bool  useFallbackValueAtEdge = false 
)

Process a set of known bad pixels in an image.

Parameters
mimageImage to patch
psfthe Image's PSF
_badListList of Defects to patch
fallbackValueValue to fallback to if all else fails
useFallbackValueAtEdgeUse the fallback value at the image's edge?

Definition at line 2047 of file Interp.cc.

2052  {
2053  /*
2054  * Allow for image's origin
2055  */
2056  int const width = mimage.getWidth();
2057  int const height = mimage.getHeight();
2058 
2059  std::vector<Defect::Ptr> badList;
2060  badList.reserve(_badList.size());
2061  for (std::vector<Defect::Ptr>::iterator ptr = _badList.begin(), end = _badList.end(); ptr != end; ++ptr) {
2062  geom::BoxI bbox = (*ptr)->getBBox();
2063  bbox.shift(geom::ExtentI(-mimage.getX0(), -mimage.getY0())); // allow for image's origin
2064  geom::PointI min = bbox.getMin(), max = bbox.getMax();
2065  if (min.getX() >= width) {
2066  continue;
2067  } else if (min.getX() < 0) {
2068  if (max.getX() < 0) {
2069  continue;
2070  } else {
2071  min.setX(0);
2072  }
2073  }
2074 
2075  if (max.getX() < 0) {
2076  continue;
2077  } else if (max.getX() >= width) {
2078  max.setX(width - 1);
2079  }
2080 
2081  bbox = geom::BoxI(min, max);
2082  Defect::Ptr ndefect(new Defect(bbox));
2083  ndefect->classify((*ptr)->getPos(), (*ptr)->getType());
2084  badList.push_back(ndefect);
2085  }
2086 
2087  sort(badList.begin(), badList.end(), Sort_ByX0<Defect>());
2088  /*
2089  * Go through the frame looking at each pixel (except the edge ones which we ignore)
2090  */
2091  typename MaskedImageT::Mask::Pixel const interpBit =
2092  mimage.getMask()->getPlaneBitMask("INTRP"); // interp'd pixels
2093 
2094  constexpr int nUseInterp = 6; // no. of pixels to interpolate towards edge
2095  static_assert(nUseInterp < Defect::WIDE_DEFECT,
2096  "make sure that we can handle these defects using"
2097  "the full interpolation not edge code");
2098 
2099  for (int y = 0; y != height; y++) {
2100  std::vector<Defect::Ptr> badList1D = classify_defects(badList, y, width);
2101 
2102  do_defects(badList1D, y, *mimage.getImage(),
2104  useFallbackValueAtEdge, nUseInterp);
2105 
2106  do_defects(badList1D, y, *mimage.getMask(), interpBit, useFallbackValueAtEdge, nUseInterp);
2107 
2108  do_defects(badList1D, y, *mimage.getVariance(),
2110  useFallbackValueAtEdge, nUseInterp);
2111  }
2112 }
int max
uint64_t * ptr
Definition: RangeSet.cc:88
Box2I BoxI
Definition: Box.h:780

◆ makeCoaddTransmissionCurve()

std::shared_ptr< afw::image::TransmissionCurve const > lsst::meas::algorithms::makeCoaddTransmissionCurve ( std::shared_ptr< afw::geom::SkyWcs const >  coaddWcs,
afw::table::ExposureCatalog const &  inputSensors 
)

Create a TransmissionCurve that represents the effective throughput on a coadd.

Parameters
[in]coaddWcsWCS that relates the coadd coordinate system to the sky.
[in]inputSensorsA catalog containing the WCSs, bounding boxes and polygons, coaddition weights (in a field called 'weight'), and TransmissionCurves of the sensor-level images that went into the coadd.
Exceptions
NotFoundErrorThrown if the 'weight' field does not exist in the schema.
InvalidParameterErrorThrown if one or more inputs do not have a TransmissionCurve or a Wcs (ValidPolygons may be null to indicate no spatial restrictions other than the bounding box).
Returns
a new TransmissionCurve object.

Definition at line 257 of file CoaddTransmissionCurve.cc.

258  {
259  return std::make_shared<CoaddTransmissionCurve>(coaddWcs, inputSensors);
260 }
afw::table::Key< int > coaddWcs

◆ makeExposurePatch() [1/2]

template<typename ExposureT >
std::shared_ptr<ExposurePatch<ExposureT> > lsst::meas::algorithms::makeExposurePatch ( std::shared_ptr< ExposureT const >  exp,
afw::detection::Footprint const &  standardFoot,
geom::Point2D const &  standardCenter,
afw::geom::SkyWcs const &  standardWcs 
)

Definition at line 98 of file ExposurePatch.h.

99  {
100  return std::make_shared<ExposurePatch<ExposureT> >(exp, standardFoot, standardCenter, standardWcs);
101 }
T exp(T... args)

◆ makeExposurePatch() [2/2]

template<typename ExposureT >
std::shared_ptr<ExposurePatch<ExposureT> > lsst::meas::algorithms::makeExposurePatch ( std::shared_ptr< ExposureT const >  exp,
std::shared_ptr< afw::detection::Footprint const >  foot,
geom::Point2D const &  center 
)

Factory function for ExposurePatch.

Definition at line 92 of file ExposurePatch.h.

93  {
94  return std::make_shared<ExposurePatch<ExposureT> >(exp, foot, center);
95 }

◆ makePsfCandidate()

template<typename PixelT >
std::shared_ptr<PsfCandidate<PixelT> > lsst::meas::algorithms::makePsfCandidate ( std::shared_ptr< afw::table::SourceRecord > const &  source,
std::shared_ptr< afw::image::Exposure< PixelT >>  image 
)

Return a PsfCandidate of the right sort.

Cf. std::make_pair

Parameters
sourceThe detected Source
imageThe image wherein lies the object

Definition at line 176 of file PsfCandidate.h.

180  {
181  return std::make_shared<PsfCandidate<PixelT>>(source, image);
182 }
afw::table::Key< afw::table::Array< ImagePixelT > > image
const char * source()
Source function that allows astChannel to source from a Stream.
Definition: Stream.h:224

◆ PYBIND11_MODULE()

lsst::meas::algorithms::PYBIND11_MODULE ( psfCandidate  ,
mod   
)

Definition at line 77 of file psfCandidate.cc.

77  {
78  declarePsfCandidate<float>(mod, "F");
79 }

◆ setSpatialParameters() [1/2]

void lsst::meas::algorithms::setSpatialParameters ( afw::math::Kernel kernel,
Eigen::VectorXd const &  vec 
)

Fit a Kernel's spatial variability from a set of stars.

Definition at line 521 of file SpatialModelPsf.cc.

521  {
522  int const nComponents = kernel->getNKernelParameters();
523  int const nSpatialParams = kernel->getNSpatialParameters();
524 
525  assert(nComponents * nSpatialParams == vec.size());
526 
527  std::vector<std::vector<double>> kCoeffs; // coefficients rearranged for Kernel
528  kCoeffs.reserve(nComponents);
529  for (int i = 0; i != nComponents; ++i) {
530  std::vector<double> spatialCoeffs(nSpatialParams);
531  for (int j = 0; j != nSpatialParams; ++j) {
532  spatialCoeffs[j] = vec[i * nSpatialParams + j];
533  }
534  kCoeffs.push_back(spatialCoeffs);
535  }
536 
537  kernel->setSpatialParameters(kCoeffs);
538 }

◆ setSpatialParameters() [2/2]

void lsst::meas::algorithms::setSpatialParameters ( afw::math::Kernel kernel,
std::vector< double > const &  coeffs 
)

Fit a Kernel's spatial variability from a set of stars.

Definition at line 500 of file SpatialModelPsf.cc.

500  {
501  int const nComponents = kernel->getNKernelParameters();
502  int const nSpatialParams = kernel->getNSpatialParameters();
503 
504  assert(nComponents * nSpatialParams == static_cast<long>(coeffs.size()));
505 
506  std::vector<std::vector<double>> kCoeffs; // coefficients rearranged for Kernel
507  kCoeffs.reserve(nComponents);
508  for (int i = 0; i != nComponents; ++i) {
509  kCoeffs.push_back(std::vector<double>(nSpatialParams));
510  std::copy(coeffs.begin() + i * nSpatialParams, coeffs.begin() + (i + 1) * nSpatialParams,
511  kCoeffs[i].begin());
512  }
513 
514  kernel->setSpatialParameters(kCoeffs);
515 }

◆ subtractPsf() [1/2]

template<typename ImageT >
double lsst::meas::algorithms::subtractPsf ( afw::detection::Psf const &  psf,
ImageT *  data,
double  x,
double  y,
double  psfFlux = std::numeric_limits< double >::quiet_NaN() 
)

◆ subtractPsf() [2/2]

template<typename MaskedImageT >
double lsst::meas::algorithms::subtractPsf ( afw::detection::Psf const &  psf,
MaskedImageT *  data,
double  x,
double  y,
double  psfFlux 
)

Subtract a PSF from an image at a given position.

Parameters
psfthe PSF to subtract
dataImage to subtract PSF from
xcolumn position
yrow position
psfFluxobject's PSF flux (if not NaN)

Definition at line 976 of file SpatialModelPsf.cc.

981  {
982  if (std::isnan(x + y)) {
984  }
985 
986  //
987  // Get Psf candidate
988  //
990 
991  //
992  // Now find the proper sub-Image
993  //
994  geom::BoxI bbox = kImage->getBBox();
995 
997  new MaskedImageT(*data, bbox, afw::image::PARENT, false)); // shallow copy
998  //
999  // Now we've got both; find the PSF's amplitude
1000  //
1001  double lambda = 0.0; // floor for variance is lambda*data
1002  try {
1003  double chi2; // chi^2 for fit
1004  double amp; // estimate of amplitude of model at this point
1005 
1006  if (std::isnan(psfFlux)) {
1007  std::pair<double, double> result = fitKernel(*kImage, *subData, lambda, true);
1008  chi2 = result.first; // chi^2 for fit
1009  amp = result.second; // estimate of amplitude of model at this point
1010  } else {
1012  amp = psfFlux / afw::math::makeStatistics(*kImage, afw::math::SUM).getValue();
1013  }
1014  //
1015  // Convert kImage to the proper type so that I can subtract it.
1016  //
1018  new typename MaskedImageT::Image(*kImage, true)); // of data's type
1019 
1020  *kImageF *= amp;
1021  *subData->getImage() -= *kImageF;
1022 
1023  return chi2;
1024  } catch (lsst::pex::exceptions::RangeError& e) {
1025  LSST_EXCEPT_ADD(e, (boost::format("Object at (%.2f, %.2f)") % x % y).str());
1026  throw e;
1027  }
1028 }
py::object result
Definition: _schema.cc:429
char * data
Definition: BaseRecord.cc:61
#define LSST_EXCEPT_ADD(e, m)
Add the current location and a message to an existing exception before rethrowing it.
Definition: Exception.h:54
double getValue(Property const prop=NOTHING) const
Return the value of the desired property (if specified in the constructor)
Definition: Statistics.cc:1047
Reports when the result of an operation cannot be represented by the destination type.
Definition: Runtime.h:115
T isnan(T... args)
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359
@ SUM
find sum of pixels in the image
Definition: Statistics.h:77
T quiet_NaN(T... args)