lsst::meas::algorithms Namespace Reference

Namespaces
	debugger
	defects
	detection
	findCosmicRaysConfig
	gaussianPsfFactory
	interp
	loadReferenceObjects
	makeCoaddApCorrMap
	measureSourceUtils
	objectSizeStarSelector
	pcaPsfDeterminer
	psfDeterminerRegistry
	psfSelectionFromMatchList
	secondMomentStarSelector
	shapelet
	sizeMagnitudeStarSelectorFactory
	starSelectorRegistry
	testUtils
	utils
	version

Classes
struct	CoaddBoundedFieldElement
	Struct used to hold one Exposure's data in a CoaddBoundedField. More...
class	CoaddBoundedField
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
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	PsfAttributes
class	PsfCandidate
	Class stored in SpatialCells for spatial Psf fitting. More...
class	Shapelet
class	ShapeletInterpolation
class	LocalShapeletKernel
class	ShapeletKernel
class	ShapeletPsfCandidate
class	SingleGaussianPsf
	Represent a PSF as a circularly symmetrical double Gaussian. More...
class	SizeMagnitudeStarSelector
class	WarpedPsf
	A Psf class that maps an arbitrary Psf through a coordinate transformation. More...
class	ShapeletImpl
class	LoadCandidatesVisitor
class	ShapeletInterpolationImpl
class	ShapeletKernelFunction
class	ShapeletSpatialFunction
class	SizeMagnitudeStarSelectorImpl
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
template<typename MaskedImageT >
std::vector< boost::shared_ptr < lsst::afw::detection::Footprint > >	findCosmicRays (MaskedImageT &image, lsst::afw::detection::Psf const &psf, double const bkgd, lsst::pex::policy::Policy const &policy, bool const keep=false)
	Find cosmic rays in an Image, and mask and remove them. More...
template<typename ExposureT >
boost::shared_ptr < ExposurePatch< ExposureT > >	makeExposurePatch (boost::shared_ptr< ExposureT const > exp, boost::shared_ptr< afw::detection::Footprint const > foot, afw::geom::Point2D const &center)
	Factory function for ExposurePatch. More...
template<typename ExposureT >
boost::shared_ptr < ExposurePatch< ExposureT > >	makeExposurePatch (boost::shared_ptr< ExposureT const > exp, afw::detection::Footprint const &standardFoot, afw::geom::Point2D const &standardCenter, afw::image::Wcs const &standardWcs)
template<typename MaskedImageT >
void	interpolateOverDefects (MaskedImageT &mimage,lsst::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 >
boost::shared_ptr < PsfCandidate< PixelT > >	makePsfCandidate (boost::shared_ptr< afw::table::SourceRecord > const &source, boost::shared_ptr< afw::image::Exposure< PixelT > > image)
Eigen::Matrix2d	getJacobian (const lsst::afw::image::Wcs &wcs, const lsst::afw::geom::PointD &pos)
	a helper function to deal with the fact that Wcs doesn't directly return a Jacobian. More...
template<typename PixelT >
std::pair < lsst::afw::math::LinearCombinationKernel::Ptr, std::vector< double > >	createKernelFromPsfCandidates (lsst::afw::math::SpatialCellSet const &psfCells, lsst::afw::geom::Extent2I const &dims, lsst::afw::geom::Point2I const &xy0, int const nEigenComponents, int const spatialOrder, int const ksize, int const nStarPerCell=-1, bool const constantWeight=true, int const border=3)
template<typename PixelT >
int	countPsfCandidates (lsst::afw::math::SpatialCellSet const &psfCells, int const nStarPerCell=-1)
template<typename PixelT >
std::pair< bool, double >	fitSpatialKernelFromPsfCandidates (lsst::afw::math::Kernel *kernel, lsst::afw::math::SpatialCellSet const &psfCells, int const nStarPerCell=-1, double const tolerance=1e-5, double const lambda=0.0)
template<typename PixelT >
std::pair< bool, double >	fitSpatialKernelFromPsfCandidates (lsst::afw::math::Kernel *kernel, lsst::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 (lsst::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 > , lsst::afw::math::KernelList >	fitKernelParamsToImage (lsst::afw::math::LinearCombinationKernel const &kernel, Image const &image, lsst::afw::geom::Point2D const &pos)
template<typename Image >
std::pair < lsst::afw::math::Kernel::Ptr, std::pair< double, double > >	fitKernelToImage (lsst::afw::math::LinearCombinationKernel const &kernel, Image const &image, lsst::afw::geom::Point2D const &pos)
afw::geom::Box2I	getOverallBBox (std::vector< boost::shared_ptr< afw::image::Image< double > >> const &imgVector)
void	addToImage (boost::shared_ptr< afw::image::Image< double > > image, std::vector< boost::shared_ptr< afw::image::Image< double > >> const &imgVector, std::vector< double > const &weightVector)
template<typename MaskedImageT >
std::vector < detection::Footprint::Ptr >	findCosmicRays (MaskedImageT &mimage, detection::Psf const &psf, double const bkgd, lsst::pex::policy::Policy const &policy, bool const keep)
	Find cosmic rays in an Image, and mask and remove them. More...
template<class ShapeletPtr >
int	getImageSize (ShapeletPtr shapelet, int size)
std::vector < lsst::afw::math::Kernel::SpatialFunctionPtr >	buildSetOfSpatialFunctions (ShapeletInterpolation::ConstPtr interp)
template<typename PixelT >
std::pair < afwMath::LinearCombinationKernel::Ptr, std::vector< double > >	createKernelFromPsfCandidates (afwMath::SpatialCellSet const &psfCells, lsst::afw::geom::Extent2I const &dims, lsst::afw::geom::Point2I const &xy0, int const nEigenComponents, int const spatialOrder, int const ksize, int const nStarPerCell, bool const constantWeight, int const border)
void	setSpatialParameters (afwMath::Kernel *kernel, std::vector< double > const &coeffs)
void	setSpatialParameters (afwMath::Kernel *kernel, Eigen::VectorXd const &vec)
template<typename MaskedImageT >
double	subtractPsf (afwDetection::Psf const &psf, MaskedImageT *data, double x, double y, double psfFlux)
template<typename Image >
std::pair< std::vector< double > , afwMath::KernelList >	fitKernelParamsToImage (afwMath::LinearCombinationKernel const &kernel, Image const &image, afwGeom::Point2D const &pos)
template<typename Image >
std::pair < afwMath::Kernel::Ptr, std::pair< double, double > >	fitKernelToImage (afwMath::LinearCombinationKernel const &kernel, Image const &image, afwGeom::Point2D const &pos)

Detailed Description

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

Typedef Documentation

typedef std::vector<Defect::Ptr>::const_iterator lsst::meas::algorithms::DefectCIter

Definition at line 53 of file Interp.cc.

Function Documentation

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

Definition at line 198 of file CoaddPsf.cc.

  {
    assert(imgVector.size() == weightVector.size());
    for (unsigned int i = 0; i < imgVector.size(); i ++) {
        PTR(afw::image::Image<double>) componentImg = imgVector[i];
        double weight = weightVector[i];
        double sum = componentImg->getArray().asEigen().sum();

        // Now get the portion of the component image which is appropriate to add
        // If the default image size is used, the component is guaranteed to fit,
        // but not if a size has been specified.
        afw::geom::Box2I cBBox = componentImg->getBBox();
        afw::geom::Box2I overlap(cBBox);
        overlap.clip(image->getBBox());
        // JFB: A subimage view of the image we want to add to, containing only the overlap region.
        afw::image::Image<double> targetSubImage(*image, overlap);
        // JFB: A subimage view of the image we want to add from, containing only the overlap region.
        afw::image::Image<double> cSubImage(*componentImg, overlap);
        targetSubImage.scaledPlus(weight/sum, cSubImage);
    }
}

std::vector<lsst::afw::math::Kernel::SpatialFunctionPtr> lsst::meas::algorithms::buildSetOfSpatialFunctions

(

ShapeletInterpolation::ConstPtr

interp

)

Definition at line 167 of file ShapeletKernel.cc.

    {
        const int nParam = interp->getFitSize();
        typedef lsst::afw::math::Kernel::SpatialFunctionPtr SFPtr;
        std::vector<SFPtr> ret;
        ret.reserve(nParam);
        for(int i=0;i<nParam;++i) {
            ret.push_back(SFPtr(new ShapeletSpatialFunction(interp,i)));
        }
        return ret;
    }

template<typename PixelT >

int lsst::meas::algorithms::countPsfCandidates	(	afwMath::SpatialCellSet const &	psfCells,
		int const	nStarPerCell
	)

Count the number of candidates in use

Definition at line 352 of file SpatialModelPsf.cc.

{
    countVisitor<PixelT> counter;
    psfCells.visitCandidates(&counter, nStarPerCell);

    return counter.getN();    
}

template<typename PixelT >

std::pair<lsst::afw::math::LinearCombinationKernel::Ptr, std::vector<double> > lsst::meas::algorithms::createKernelFromPsfCandidates	(	afwMath::SpatialCellSet const &	psfCells,
		lsst::afw::geom::Extent2I const &	dims,
		lsst::afw::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::Ptr 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

psfCells	A SpatialCellSet containing PsfCandidates
dims	Dimensions of image
xy0	Origin of image
nEigenComponents	number of eigen components to keep; <= 0 => infty
spatialOrder	Order of spatial variation (cf. afw::math::PolynomialFunction2)
ksize	Size of generated Kernel images
nStarPerCell	max no. of stars per cell; <= 0 => infty
constantWeight	should each star have equal weight in the fit?
border	Border size for background subtraction

Definition at line 206 of file SpatialModelPsf.cc.

{
    typedef typename afwImage::Image<PixelT> ImageT;
    typedef typename afwImage::MaskedImage<PixelT> MaskedImageT;
    typedef typename afwImage::Exposure<PixelT> ExposureT;
    
    //
    // Set the sizes for PsfCandidates made from either Images or MaskedImages
    //
    //lsst::meas::algorithms::PsfCandidate<ImageT>::setWidth(ksize);
    //lsst::meas::algorithms::PsfCandidate<ImageT>::setHeight(ksize);
    //lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setWidth(ksize);
    //lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setHeight(ksize);
    lsst::meas::algorithms::PsfCandidate<PixelT>::setWidth(ksize);
    lsst::meas::algorithms::PsfCandidate<PixelT>::setHeight(ksize);

    
    PsfImagePca<MaskedImageT> imagePca(constantWeight, border); // Here's the set of images we'll analyze

    {
        SetPcaImageVisitor<PixelT> importStarVisitor(&imagePca);
        bool const ignoreExceptions = true;
        psfCells.visitCandidates(&importStarVisitor, nStarPerCell, ignoreExceptions);
    }

    //
    // Do a PCA decomposition of those PSF candidates.
    //
    // We have "gappy" data;  in other words we don't want to include any pixels with INTRP set
    //
    int niter = 10;                     // number of iterations of updateBadPixels
    double deltaLim = 10.0;             // acceptable value of delta, the max change due to updateBadPixels
    lsst::afw::image::MaskPixel const BAD = afwImage::Mask<>::getPlaneBitMask("BAD");
    lsst::afw::image::MaskPixel const CR = afwImage::Mask<>::getPlaneBitMask("CR");
    lsst::afw::image::MaskPixel const INTRP = afwImage::Mask<>::getPlaneBitMask("INTRP");
    
    for (int i = 0; i != niter; ++i) {
        int const ncomp = (i == 0) ? 0 :
            ((nEigenComponents == 0) ? imagePca.getEigenImages().size() : nEigenComponents);
        double delta = imagePca.updateBadPixels(BAD | CR | INTRP, ncomp);
        if (i > 0 && delta < deltaLim) {
            break;
        }
        
        imagePca.analyze();
    }
    
    std::vector<typename MaskedImageT::Ptr> eigenImages = imagePca.getEigenImages();
    std::vector<double> eigenValues = imagePca.getEigenValues();
    int const nEigen = static_cast<int>(eigenValues.size());
    
    int const ncomp = (nEigenComponents <= 0 || nEigen < nEigenComponents) ? nEigen : nEigenComponents;
    //
    // Set the background level of the components to 0.0 to avoid coupling variable background
    // levels to the form of the Psf.  More precisely, we calculate the mean of an outer "annulus"
    // of width bkg_border
    //
    for (int k = 0; k != ncomp; ++k) {
        ImageT const& im = *eigenImages[k]->getImage();

        int bkg_border = 2;
        if (bkg_border > im.getWidth()) {
            bkg_border = im.getWidth() / 2;
        }
        if (bkg_border > im.getHeight()) {
            bkg_border = im.getHeight() / 2;
        }

        double sum = 0;
        // Bottom and Top borders
        for (int i = 0; i != bkg_border; ++i) {
            typename ImageT::const_x_iterator
                ptrB = im.row_begin(i), ptrT = im.row_begin(im.getHeight() - i - 1);
            for (int j = 0; j != im.getWidth(); ++j, ++ptrB, ++ptrT) {
                sum += *ptrB + *ptrT;
            }
        }
        for (int i = bkg_border; i < im.getHeight() - bkg_border; ++i) {
            // Left and Right borders
            typename ImageT::const_x_iterator
                ptrL = im.row_begin(i), ptrR = im.row_begin(i) + im.getWidth() - bkg_border;
            for (int j = 0; j != bkg_border; ++j, ++ptrL, ++ptrR) {
                sum += *ptrL + *ptrR;
            }
        }
        sum /= 2*(bkg_border*im.getWidth() + bkg_border*(im.getHeight() - 2*bkg_border));

        *eigenImages[k] -= sum;
    }
    //
    // Now build our LinearCombinationKernel; build the lists of basis functions
    // and spatial variation, then assemble the Kernel
    //
    afwMath::KernelList  kernelList;
    std::vector<afwMath::Kernel::SpatialFunctionPtr> spatialFunctionList;
    afwGeom::Box2D const range = afwGeom::Box2D(afwGeom::Point2D(xy0), afwGeom::Extent2D(dims));

    for (int i = 0; i != ncomp; ++i) {
        {
            // Enforce unit sum for kernel by construction
            // Zeroth component has unit sum
            // Other components have zero sum by normalising and then subtracting the zeroth component
            ImageT& image = *eigenImages[i]->getImage();
            double sum = std::accumulate(image.begin(true), image.end(true), 0.0);
            if (i == 0) {
                image /= sum;
            } else {
                for (typename ImageT::fast_iterator ptr0 = eigenImages[0]->getImage()->begin(true),
                         ptr1 = image.begin(true), end = image.end(true); ptr1 != end; ++ptr0, ++ptr1) {
                    *ptr1 = *ptr1 / sum - *ptr0;
                }
            }
        }

        kernelList.push_back(afwMath::Kernel::Ptr(new afwMath::FixedKernel(
                                      afwImage::Image<afwMath::Kernel::Pixel>(*eigenImages[i]->getImage(),true)
                                                                          )));

        afwMath::Kernel::SpatialFunctionPtr
//            spatialFunction(new afwMath::PolynomialFunction2<double>(spatialOrder));
          spatialFunction(new afwMath::Chebyshev1Function2<double>(spatialOrder, range));
        spatialFunction->setParameter(0, 1.0); // the constant term; all others are 0
        spatialFunctionList.push_back(spatialFunction);
    }

    afwMath::LinearCombinationKernel::Ptr
        psf(new afwMath::LinearCombinationKernel(kernelList, spatialFunctionList));

    return std::make_pair(psf, eigenValues);
}

template<typename PixelT >

std::pair<afwMath::LinearCombinationKernel::Ptr, std::vector<double> > lsst::meas::algorithms::createKernelFromPsfCandidates	(	afwMath::SpatialCellSet const &	psfCells,
		lsst::afw::geom::Extent2I const &	dims,
		lsst::afw::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::Ptr 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

psfCells	A SpatialCellSet containing PsfCandidates
dims	Dimensions of image
xy0	Origin of image
nEigenComponents	number of eigen components to keep; <= 0 => infty
spatialOrder	Order of spatial variation (cf. afw::math::PolynomialFunction2)
ksize	Size of generated Kernel images
nStarPerCell	max no. of stars per cell; <= 0 => infty
constantWeight	should each star have equal weight in the fit?
border	Border size for background subtraction

Definition at line 206 of file SpatialModelPsf.cc.

{
    typedef typename afwImage::Image<PixelT> ImageT;
    typedef typename afwImage::MaskedImage<PixelT> MaskedImageT;
    typedef typename afwImage::Exposure<PixelT> ExposureT;
    
    //
    // Set the sizes for PsfCandidates made from either Images or MaskedImages
    //
    //lsst::meas::algorithms::PsfCandidate<ImageT>::setWidth(ksize);
    //lsst::meas::algorithms::PsfCandidate<ImageT>::setHeight(ksize);
    //lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setWidth(ksize);
    //lsst::meas::algorithms::PsfCandidate<MaskedImageT>::setHeight(ksize);
    lsst::meas::algorithms::PsfCandidate<PixelT>::setWidth(ksize);
    lsst::meas::algorithms::PsfCandidate<PixelT>::setHeight(ksize);

    
    PsfImagePca<MaskedImageT> imagePca(constantWeight, border); // Here's the set of images we'll analyze

    {
        SetPcaImageVisitor<PixelT> importStarVisitor(&imagePca);
        bool const ignoreExceptions = true;
        psfCells.visitCandidates(&importStarVisitor, nStarPerCell, ignoreExceptions);
    }

    //
    // Do a PCA decomposition of those PSF candidates.
    //
    // We have "gappy" data;  in other words we don't want to include any pixels with INTRP set
    //
    int niter = 10;                     // number of iterations of updateBadPixels
    double deltaLim = 10.0;             // acceptable value of delta, the max change due to updateBadPixels
    lsst::afw::image::MaskPixel const BAD = afwImage::Mask<>::getPlaneBitMask("BAD");
    lsst::afw::image::MaskPixel const CR = afwImage::Mask<>::getPlaneBitMask("CR");
    lsst::afw::image::MaskPixel const INTRP = afwImage::Mask<>::getPlaneBitMask("INTRP");
    
    for (int i = 0; i != niter; ++i) {
        int const ncomp = (i == 0) ? 0 :
            ((nEigenComponents == 0) ? imagePca.getEigenImages().size() : nEigenComponents);
        double delta = imagePca.updateBadPixels(BAD | CR | INTRP, ncomp);
        if (i > 0 && delta < deltaLim) {
            break;
        }
        
        imagePca.analyze();
    }
    
    std::vector<typename MaskedImageT::Ptr> eigenImages = imagePca.getEigenImages();
    std::vector<double> eigenValues = imagePca.getEigenValues();
    int const nEigen = static_cast<int>(eigenValues.size());
    
    int const ncomp = (nEigenComponents <= 0 || nEigen < nEigenComponents) ? nEigen : nEigenComponents;
    //
    // Set the background level of the components to 0.0 to avoid coupling variable background
    // levels to the form of the Psf.  More precisely, we calculate the mean of an outer "annulus"
    // of width bkg_border
    //
    for (int k = 0; k != ncomp; ++k) {
        ImageT const& im = *eigenImages[k]->getImage();

        int bkg_border = 2;
        if (bkg_border > im.getWidth()) {
            bkg_border = im.getWidth() / 2;
        }
        if (bkg_border > im.getHeight()) {
            bkg_border = im.getHeight() / 2;
        }

        double sum = 0;
        // Bottom and Top borders
        for (int i = 0; i != bkg_border; ++i) {
            typename ImageT::const_x_iterator
                ptrB = im.row_begin(i), ptrT = im.row_begin(im.getHeight() - i - 1);
            for (int j = 0; j != im.getWidth(); ++j, ++ptrB, ++ptrT) {
                sum += *ptrB + *ptrT;
            }
        }
        for (int i = bkg_border; i < im.getHeight() - bkg_border; ++i) {
            // Left and Right borders
            typename ImageT::const_x_iterator
                ptrL = im.row_begin(i), ptrR = im.row_begin(i) + im.getWidth() - bkg_border;
            for (int j = 0; j != bkg_border; ++j, ++ptrL, ++ptrR) {
                sum += *ptrL + *ptrR;
            }
        }
        sum /= 2*(bkg_border*im.getWidth() + bkg_border*(im.getHeight() - 2*bkg_border));

        *eigenImages[k] -= sum;
    }
    //
    // Now build our LinearCombinationKernel; build the lists of basis functions
    // and spatial variation, then assemble the Kernel
    //
    afwMath::KernelList  kernelList;
    std::vector<afwMath::Kernel::SpatialFunctionPtr> spatialFunctionList;
    afwGeom::Box2D const range = afwGeom::Box2D(afwGeom::Point2D(xy0), afwGeom::Extent2D(dims));

    for (int i = 0; i != ncomp; ++i) {
        {
            // Enforce unit sum for kernel by construction
            // Zeroth component has unit sum
            // Other components have zero sum by normalising and then subtracting the zeroth component
            ImageT& image = *eigenImages[i]->getImage();
            double sum = std::accumulate(image.begin(true), image.end(true), 0.0);
            if (i == 0) {
                image /= sum;
            } else {
                for (typename ImageT::fast_iterator ptr0 = eigenImages[0]->getImage()->begin(true),
                         ptr1 = image.begin(true), end = image.end(true); ptr1 != end; ++ptr0, ++ptr1) {
                    *ptr1 = *ptr1 / sum - *ptr0;
                }
            }
        }

        kernelList.push_back(afwMath::Kernel::Ptr(new afwMath::FixedKernel(
                                      afwImage::Image<afwMath::Kernel::Pixel>(*eigenImages[i]->getImage(),true)
                                                                          )));

        afwMath::Kernel::SpatialFunctionPtr
//            spatialFunction(new afwMath::PolynomialFunction2<double>(spatialOrder));
          spatialFunction(new afwMath::Chebyshev1Function2<double>(spatialOrder, range));
        spatialFunction->setParameter(0, 1.0); // the constant term; all others are 0
        spatialFunctionList.push_back(spatialFunction);
    }

    afwMath::LinearCombinationKernel::Ptr
        psf(new afwMath::LinearCombinationKernel(kernelList, spatialFunctionList));

    return std::make_pair(psf, eigenValues);
}

template<typename MaskedImageT >

std::vector<boost::shared_ptr<lsst::afw::detection::Footprint> > lsst::meas::algorithms::findCosmicRays	(	MaskedImageT &	mimage,
		detection::Psf const &	psf,
		double const	bkgd,
		lsst::pex::policy::Policy const &	policy,
		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.

Parameters

mimage	Image to search
psf	the Image's PSF
bkgd	unsubtracted background of frame, DN
policy	Policy directing the behavior
keep	if true, don't remove the CRs

Definition at line 341 of file CR.cc.

                {
    typedef typename MaskedImageT::Image ImageT;
    typedef typename ImageT::Pixel ImagePixel;
    typedef typename MaskedImageT::Mask::Pixel MaskPixel;

    // Parse the Policy
    double const minSigma = policy.getDouble("minSigma");    // min sigma over sky in pixel for CR candidate
    double const minDn = policy.getDouble("min_DN");         // min number of DN in an CRs
    double const cond3Fac = policy.getDouble("cond3_fac");   // fiddle factor for condition #3
    double const cond3Fac2 = policy.getDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
    int const niteration = policy.getInt("niteration");      // Number of times to look for contaminated
                                                             // pixels near CRs
    int const nCrPixelMax = policy.getInt("nCrPixelMax");    // maximum number of contaminated pixels
/*
 * thresholds for 3rd condition
 *
 * Realise PSF at center of image
 */
    lsst::afw::math::Kernel::ConstPtr kernel = psf.getLocalKernel(
        afw::geom::Point2D(mimage.getWidth() / 2.0, mimage.getHeight() / 2.0)
    );
    if (!kernel) {
        throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
    }
    detection::Psf::Image psfImage = detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
    kernel->computeImage(psfImage, true);

    int const xc = kernel->getCtrX();   // center of PSF
    int const yc = kernel->getCtrY();

    double const I0 = psfImage(xc, yc);
    double const thresH = cond3Fac2*(0.5*(psfImage(xc - 1, yc) + psfImage(xc + 1, yc)))/I0; // horizontal
    double const thresV = cond3Fac2*(0.5*(psfImage(xc, yc - 1) + psfImage(xc, yc + 1)))/I0; // vertical
    double const thresD = cond3Fac2*(0.25*(psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
                                           psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1)))/I0; // diag
/*
 * Setup desired mask planes
 */
    MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
    MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
    MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
    MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels

    MaskPixel const badMask = (badBit | interpBit | saturBit); // naughty pixels
/*
 * Go through the frame looking at each pixel (except the edge ones which we ignore)
 */
    int const ncol = mimage.getWidth();
    int const nrow = mimage.getHeight();

    std::vector<CRPixel<ImagePixel> > crpixels; // storage for detected CR-contaminated pixels
    typedef typename std::vector<CRPixel<ImagePixel> >::iterator crpixel_iter;
    typedef typename std::vector<CRPixel<ImagePixel> >::reverse_iterator crpixel_riter;

    for (int j = 1; j < nrow - 1; ++j) {
        typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data

        for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
            ImagePixel corr = 0;
            if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma,
                                           thresH, thresV, thresD, bkgd, cond3Fac)) {
                continue;
            }
/*
 * condition #4
 */
            if (loc.mask() & badMask) {
                continue;
            }
            if ((loc.mask(-1,  1) | loc.mask(0,  1) | loc.mask(1,  1) |
                 loc.mask(-1,  0) |                   loc.mask(1,  0) |
                 loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) & interpBit) {
                continue;
            }
/*
 * OK, it's a CR
 *
 * replace CR-contaminated pixels with reasonable values as we go through
 * image, which increases the detection rate
 */
            crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
            loc.image() = corr;         /* just a preliminary estimate */

            if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
                reinstateCrPixels(mimage.getImage().get(), crpixels);

                throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                                  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
            }
        }
    }
/*
 * We've found them on a pixel-by-pixel basis, now merge those pixels
 * into cosmic rays
 */
    std::vector<int> aliases;           // aliases for initially disjoint parts of CRs
    aliases.reserve(1 + crpixels.size()/2); // initial size of aliases
    
    std::vector<detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
    spans.reserve(aliases.capacity());  // initial size of spans
    
    aliases.push_back(0);               // 0 --> 0

    int ncr = 0;                        // number of detected cosmic rays
    if (!crpixels.empty()) {
        int id;                         // id number for a CR
        int x0 = -1, x1 = -1, y = -1;   // the beginning and end column, and row of this span in a CR

        // I am dummy
        CRPixel<ImagePixel> dummy(0, -1, 0, -1);
        crpixels.push_back(dummy);
        //printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id, dummy.col, dummy.row, (double)dummy.val);
        for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1 ; ++crp) {
            //printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col, crp->row, (double)crp->val);

            if (crp->id < 0) {           // not already assigned
                crp->id = ++ncr;        // a new CR
                aliases.push_back(crp->id);
                y = crp->row;
                x0 = x1 = crp->col;
                //printf("  Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
            }
            id = crp->id;
            //printf("  Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row, crp[1].col);

            if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
                //printf("  Adjoining!  Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
                crp[1].id = id;
                ++x1;
            } else {
                assert (y >= 0 && x0 >= 0 && x1 >= 0);
                spans.push_back(detection::IdSpan::Ptr(new detection::IdSpan(id, y, x0, x1)));
                //printf("  Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
            }
        }
    }

    // At the end of this loop, all crpixel entries have been assigned an ID,
    // except for the "dummy" entry at the end of the array.
    if (crpixels.size() > 0) {
        for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
            assert(cp->id >= 0);
            assert(cp->col >= 0);
            assert(cp->row >= 0);
        }
        // dummy:
        assert(crpixels[crpixels.size()-1].id == -1);
        assert(crpixels[crpixels.size()-1].col == 0);
        assert(crpixels[crpixels.size()-1].row == -1);
    }

    for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end; sp++) {
        assert((*sp)->id >= 0);
        assert((*sp)->y >= 0);
        assert((*sp)->x0 >= 0);
        assert((*sp)->x1 >= (*sp)->x0);
        for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
            assert((*sp2)->y >= (*sp)->y);
            if ((*sp2)->y == (*sp)->y) {
                assert((*sp2)->x0 > (*sp)->x1);
            }
        }
    }

/*
 * See if spans touch each other
 */
    for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end();
         sp != end; ++sp) {
        int const y = (*sp)->y;
        int const x0 = (*sp)->x0;
        int const x1 = (*sp)->x1;

        // this loop will probably run for only a few steps
        for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
            if ((*sp2)->y == y) {
                // on this row (but not adjoining columns, since it would have been merged into this span);
                // keep looking.
                continue;
            } else if ((*sp2)->y != (y + 1)) {
                // sp2 is more than one row below; can't be connected.
                break;
            } else if ((*sp2)->x0 > (x1 + 1)) {
                // sp2 is more than one column away to the right; can't be connected
                break;
            } else if ((*sp2)->x1 >= (x0 - 1)) {
                // touches
                int r1 = detection::resolve_alias(aliases, (*sp)->id);
                int r2 = detection::resolve_alias(aliases, (*sp2)->id);
                aliases[r1] = r2;
            }
        }
    }





/*
 * Resolve aliases; first alias chains, then the IDs in the spans
 */
    for (unsigned int i = 0; i != spans.size(); ++i) {
        spans[i]->id = detection::resolve_alias(aliases, spans[i]->id);
    }

/*
 * Sort spans by ID, so we can sweep through them once
 */
    if (spans.size() > 0) {
        std::sort(spans.begin(), spans.end(), detection::IdSpanCompar());
    }

/*
 * Build Footprints from spans
 */
    std::vector<detection::Footprint::Ptr> CRs; // our cosmic rays

    if (spans.size() > 0) {
        int id = spans[0]->id;
        unsigned int i0 = 0;            // initial value of i
        for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
            if (i == spans.size() || spans[i]->id != id) {
                detection::Footprint::Ptr cr(new detection::Footprint(i - i0));
                
                for (; i0 < i; ++i0) {
                    cr->addSpan(spans[i0]->y, spans[i0]->x0, spans[i0]->x1);
                }
                CRs.push_back(cr);
            }
            
            if (i < spans.size()) {
                id = spans[i]->id;
            }
        }
    }

    reinstateCrPixels(mimage.getImage().get(), crpixels);
/*
 * apply condition #1
 */
    CountsInCR<ImageT> CountDN(*mimage.getImage(), bkgd);
    for (std::vector<detection::Footprint::Ptr>::iterator cr = CRs.begin(), end = CRs.end();
         cr != end; ++cr) {
        CountDN.apply(**cr);            // find the sum of pixel values within the CR
                
        pexLogging::TTrace<10>("algorithms.CR", "CR at (%d, %d) has %g DN",
                               (*cr)->getBBox().getMinX(), (*cr)->getBBox().getMinY(), CountDN.getCounts());
        if (CountDN.getCounts() < minDn) { /* not bright enough */
            pexLogging::TTrace<11>("algorithms.CR", "Erasing CR");

            cr = CRs.erase(cr);
            --cr;                       // back up to previous CR (we're going to increment it)
            --end;
        }
    }
    ncr = CRs.size();           /* some may have been too faint */
/*
 * We've found them all, time to kill them all
 */
    bool const debias_values = true;
    bool grow = false;
    pexLogging::TTrace<2>("algorithms.CR", "Removing initial list of CRs");
    removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
#if 0                                   // Useful to see phase 2 in ds9; debugging only
    (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
                                   mimage.getMask()->getPlaneBitMask("DETECTED"));
#endif
/*
 * Now that we've removed them, go through image again, examining area around
 * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
 *
 * We iterate niteration times;  niter==1 was sufficient for SDSS data, but megacam
 * CCDs are different -- who knows for other devices?
 */
    bool too_many_crs = false;          // we've seen too many CR pixels
    int nextra = 0;                     // number of pixels added to list of CRs
    for (int i = 0; i != niteration && !too_many_crs; ++i) {
        pexLogging::TTrace<1>("algorithms.CR", "Starting iteration %d", i);
        for (std::vector<detection::Footprint::Ptr>::iterator fiter = CRs.begin();
             fiter != CRs.end(); fiter++) {
            detection::Footprint::Ptr cr = *fiter;
/*
 * Are all those `CR' pixels interpolated?  If so, don't grow it
 */
            {
                detection::Footprint::Ptr om = footprintAndMask(cr, mimage.getMask(), interpBit);
                int const npix = (om) ? om->getNpix() : 0;

                if (npix == cr->getNpix()) {
                    continue;
                }
            }
/*
 * No; some of the suspect pixels aren't interpolated
 */
            detection::Footprint extra;                     // extra pixels added to cr
            for (detection::Footprint::SpanList::const_iterator siter = cr->getSpans().begin();
                 siter != cr->getSpans().end(); siter++) {
                detection::Span::Ptr const span = *siter;

                /*
                 * Check the lines above and below the span.  We're going to check a 3x3 region around
                 * the pixels, so we need a buffer around the edge.  We check the pixels just to the
                 * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
                 * column direction, but only 1 in the row direction.
                 */
                int const y = span->getY() - mimage.getY0();
                if (y < 2 || y >= nrow - 2) {
                    continue;
                }
                int x0 = span->getX0() - mimage.getX0();
                int x1 = span->getX1() - mimage.getX0();
                x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
                x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;

                checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
                checkSpanForCRs(&extra, crpixels, y,     x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
                checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
            }

            if (extra.getSpans().size() > 0) {      // we added some pixels
                if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
                    too_many_crs = true;
                    break;
                }

                nextra += extra.getNpix();

                detection::Footprint::SpanList &espans = extra.getSpans();
                for (detection::Footprint::SpanList::const_iterator siter = espans.begin();
                     siter != espans.end(); siter++) {
                    cr->addSpan(**siter);
                }
                                cr->normalize();
            }            
        }

        if (nextra == 0) {
            break;
        }
    }
/*
 * mark those pixels as CRs
 */
    if (!too_many_crs) {
        (void)setMaskFromFootprintList(mimage.getMask().get(), CRs, crBit);
    }
/*
 * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
 * first value stored (hence the uses of rbegin/rend)
 *
 * We have to do this if we decide _not_ to remove certain CRs,
 * for example those which lie next to saturated pixels
 */
    if (keep || too_many_crs) {
        if (crpixels.size() > 0) {
            int const imageX0 = mimage.getX0();
            int const imageY0 = mimage.getY0();

            std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
        
            crpixel_riter rend = crpixels.rend();
            for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
                if (crp->row == -1)
                    // dummy; skip it.
                    continue;
                mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
            }
        }
    } else {
        if (true || nextra > 0) {
            grow = true;
            pexLogging::TTrace<2>("algorithms.CR", "Removing final list of CRs, grow = %d", grow);
            removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
        }
/*
 * we interpolated over all CR pixels, so set the interp bits too
 */
        (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
                                       static_cast<MaskPixel>(crBit | interpBit));
    }

    if (too_many_crs) {                 // we've cleaned up, so we can throw the exception
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
    }

    return CRs;
}

template<typename MaskedImageT >

std::vector<detection::Footprint::Ptr> lsst::meas::algorithms::findCosmicRays	(	MaskedImageT &	mimage,
		detection::Psf const &	psf,
		double const	bkgd,
		lsst::pex::policy::Policy const &	policy,
		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.

Parameters

mimage	Image to search
psf	the Image's PSF
bkgd	unsubtracted background of frame, DN
policy	Policy directing the behavior
keep	if true, don't remove the CRs

Definition at line 341 of file CR.cc.

                {
    typedef typename MaskedImageT::Image ImageT;
    typedef typename ImageT::Pixel ImagePixel;
    typedef typename MaskedImageT::Mask::Pixel MaskPixel;

    // Parse the Policy
    double const minSigma = policy.getDouble("minSigma");    // min sigma over sky in pixel for CR candidate
    double const minDn = policy.getDouble("min_DN");         // min number of DN in an CRs
    double const cond3Fac = policy.getDouble("cond3_fac");   // fiddle factor for condition #3
    double const cond3Fac2 = policy.getDouble("cond3_fac2"); // 2nd fiddle factor for condition #3
    int const niteration = policy.getInt("niteration");      // Number of times to look for contaminated
                                                             // pixels near CRs
    int const nCrPixelMax = policy.getInt("nCrPixelMax");    // maximum number of contaminated pixels
/*
 * thresholds for 3rd condition
 *
 * Realise PSF at center of image
 */
    lsst::afw::math::Kernel::ConstPtr kernel = psf.getLocalKernel(
        afw::geom::Point2D(mimage.getWidth() / 2.0, mimage.getHeight() / 2.0)
    );
    if (!kernel) {
        throw LSST_EXCEPT(pexExcept::NotFoundError, "Psf is unable to return a kernel");
    }
    detection::Psf::Image psfImage = detection::Psf::Image(geom::ExtentI(kernel->getWidth(), kernel->getHeight()));
    kernel->computeImage(psfImage, true);

    int const xc = kernel->getCtrX();   // center of PSF
    int const yc = kernel->getCtrY();

    double const I0 = psfImage(xc, yc);
    double const thresH = cond3Fac2*(0.5*(psfImage(xc - 1, yc) + psfImage(xc + 1, yc)))/I0; // horizontal
    double const thresV = cond3Fac2*(0.5*(psfImage(xc, yc - 1) + psfImage(xc, yc + 1)))/I0; // vertical
    double const thresD = cond3Fac2*(0.25*(psfImage(xc - 1, yc - 1) + psfImage(xc + 1, yc + 1) +
                                           psfImage(xc - 1, yc + 1) + psfImage(xc + 1, yc - 1)))/I0; // diag
/*
 * Setup desired mask planes
 */
    MaskPixel const badBit = mimage.getMask()->getPlaneBitMask("BAD"); // Generic bad pixels
    MaskPixel const crBit = mimage.getMask()->getPlaneBitMask("CR"); // CR-contaminated pixels
    MaskPixel const interpBit = mimage.getMask()->getPlaneBitMask("INTRP"); // Interpolated pixels
    MaskPixel const saturBit = mimage.getMask()->getPlaneBitMask("SAT"); // Saturated pixels

    MaskPixel const badMask = (badBit | interpBit | saturBit); // naughty pixels
/*
 * Go through the frame looking at each pixel (except the edge ones which we ignore)
 */
    int const ncol = mimage.getWidth();
    int const nrow = mimage.getHeight();

    std::vector<CRPixel<ImagePixel> > crpixels; // storage for detected CR-contaminated pixels
    typedef typename std::vector<CRPixel<ImagePixel> >::iterator crpixel_iter;
    typedef typename std::vector<CRPixel<ImagePixel> >::reverse_iterator crpixel_riter;

    for (int j = 1; j < nrow - 1; ++j) {
        typename MaskedImageT::xy_locator loc = mimage.xy_at(1, j); // locator for data

        for (int i = 1; i < ncol - 1; ++i, ++loc.x()) {
            ImagePixel corr = 0;
            if (!is_cr_pixel<MaskedImageT>(&corr, loc, minSigma,
                                           thresH, thresV, thresD, bkgd, cond3Fac)) {
                continue;
            }
/*
 * condition #4
 */
            if (loc.mask() & badMask) {
                continue;
            }
            if ((loc.mask(-1,  1) | loc.mask(0,  1) | loc.mask(1,  1) |
                 loc.mask(-1,  0) |                   loc.mask(1,  0) |
                 loc.mask(-1, -1) | loc.mask(0, -1) | loc.mask(1, -1)) & interpBit) {
                continue;
            }
/*
 * OK, it's a CR
 *
 * replace CR-contaminated pixels with reasonable values as we go through
 * image, which increases the detection rate
 */
            crpixels.push_back(CRPixel<ImagePixel>(i + mimage.getX0(), j + mimage.getY0(), loc.image()));
            loc.image() = corr;         /* just a preliminary estimate */

            if (static_cast<int>(crpixels.size()) > nCrPixelMax) {
                reinstateCrPixels(mimage.getImage().get(), crpixels);

                throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                                  (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
            }
        }
    }
/*
 * We've found them on a pixel-by-pixel basis, now merge those pixels
 * into cosmic rays
 */
    std::vector<int> aliases;           // aliases for initially disjoint parts of CRs
    aliases.reserve(1 + crpixels.size()/2); // initial size of aliases
    
    std::vector<detection::IdSpan::Ptr> spans; // y:x0,x1 for objects
    spans.reserve(aliases.capacity());  // initial size of spans
    
    aliases.push_back(0);               // 0 --> 0

    int ncr = 0;                        // number of detected cosmic rays
    if (!crpixels.empty()) {
        int id;                         // id number for a CR
        int x0 = -1, x1 = -1, y = -1;   // the beginning and end column, and row of this span in a CR

        // I am dummy
        CRPixel<ImagePixel> dummy(0, -1, 0, -1);
        crpixels.push_back(dummy);
        //printf("Created dummy CR: i %i, id %i, col %i, row %i, val %g\n", dummy.get_i(), dummy.id, dummy.col, dummy.row, (double)dummy.val);
        for (crpixel_iter crp = crpixels.begin(); crp < crpixels.end() - 1 ; ++crp) {
            //printf("Looking at CR: i %i, id %i, col %i, row %i, val %g\n", crp->get_i(), crp->id, crp->col, crp->row, (double)crp->val);

            if (crp->id < 0) {           // not already assigned
                crp->id = ++ncr;        // a new CR
                aliases.push_back(crp->id);
                y = crp->row;
                x0 = x1 = crp->col;
                //printf("  Assigned ID %i; looking at row %i, start col %i\n", crp->id, crp->row, crp->col);
            }
            id = crp->id;
            //printf("  Next CRpix has i=%i, id=%i, row %i, col %i\n", crp[1].get_i(), crp[1].id, crp[1].row, crp[1].col);

            if (crp[1].row == crp[0].row && crp[1].col == crp[0].col + 1) {
                //printf("  Adjoining!  Set next CRpix id = %i; x1=%i\n", crp[1].id, x1);
                crp[1].id = id;
                ++x1;
            } else {
                assert (y >= 0 && x0 >= 0 && x1 >= 0);
                spans.push_back(detection::IdSpan::Ptr(new detection::IdSpan(id, y, x0, x1)));
                //printf("  Not adjoining; adding span id=%i, y=%i, x = [%i, %i]\n", id, y, x0, x1);
            }
        }
    }

    // At the end of this loop, all crpixel entries have been assigned an ID,
    // except for the "dummy" entry at the end of the array.
    if (crpixels.size() > 0) {
        for (crpixel_iter cp = crpixels.begin(); cp != crpixels.end() - 1; cp++) {
            assert(cp->id >= 0);
            assert(cp->col >= 0);
            assert(cp->row >= 0);
        }
        // dummy:
        assert(crpixels[crpixels.size()-1].id == -1);
        assert(crpixels[crpixels.size()-1].col == 0);
        assert(crpixels[crpixels.size()-1].row == -1);
    }

    for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end(); sp != end; sp++) {
        assert((*sp)->id >= 0);
        assert((*sp)->y >= 0);
        assert((*sp)->x0 >= 0);
        assert((*sp)->x1 >= (*sp)->x0);
        for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; sp2++) {
            assert((*sp2)->y >= (*sp)->y);
            if ((*sp2)->y == (*sp)->y) {
                assert((*sp2)->x0 > (*sp)->x1);
            }
        }
    }

/*
 * See if spans touch each other
 */
    for (std::vector<detection::IdSpan::Ptr>::iterator sp = spans.begin(), end = spans.end();
         sp != end; ++sp) {
        int const y = (*sp)->y;
        int const x0 = (*sp)->x0;
        int const x1 = (*sp)->x1;

        // this loop will probably run for only a few steps
        for (std::vector<detection::IdSpan::Ptr>::iterator sp2 = sp + 1; sp2 != end; ++sp2) {
            if ((*sp2)->y == y) {
                // on this row (but not adjoining columns, since it would have been merged into this span);
                // keep looking.
                continue;
            } else if ((*sp2)->y != (y + 1)) {
                // sp2 is more than one row below; can't be connected.
                break;
            } else if ((*sp2)->x0 > (x1 + 1)) {
                // sp2 is more than one column away to the right; can't be connected
                break;
            } else if ((*sp2)->x1 >= (x0 - 1)) {
                // touches
                int r1 = detection::resolve_alias(aliases, (*sp)->id);
                int r2 = detection::resolve_alias(aliases, (*sp2)->id);
                aliases[r1] = r2;
            }
        }
    }





/*
 * Resolve aliases; first alias chains, then the IDs in the spans
 */
    for (unsigned int i = 0; i != spans.size(); ++i) {
        spans[i]->id = detection::resolve_alias(aliases, spans[i]->id);
    }

/*
 * Sort spans by ID, so we can sweep through them once
 */
    if (spans.size() > 0) {
        std::sort(spans.begin(), spans.end(), detection::IdSpanCompar());
    }

/*
 * Build Footprints from spans
 */
    std::vector<detection::Footprint::Ptr> CRs; // our cosmic rays

    if (spans.size() > 0) {
        int id = spans[0]->id;
        unsigned int i0 = 0;            // initial value of i
        for (unsigned int i = i0; i <= spans.size(); ++i) { // <= size to catch the last object
            if (i == spans.size() || spans[i]->id != id) {
                detection::Footprint::Ptr cr(new detection::Footprint(i - i0));
                
                for (; i0 < i; ++i0) {
                    cr->addSpan(spans[i0]->y, spans[i0]->x0, spans[i0]->x1);
                }
                CRs.push_back(cr);
            }
            
            if (i < spans.size()) {
                id = spans[i]->id;
            }
        }
    }

    reinstateCrPixels(mimage.getImage().get(), crpixels);
/*
 * apply condition #1
 */
    CountsInCR<ImageT> CountDN(*mimage.getImage(), bkgd);
    for (std::vector<detection::Footprint::Ptr>::iterator cr = CRs.begin(), end = CRs.end();
         cr != end; ++cr) {
        CountDN.apply(**cr);            // find the sum of pixel values within the CR
                
        pexLogging::TTrace<10>("algorithms.CR", "CR at (%d, %d) has %g DN",
                               (*cr)->getBBox().getMinX(), (*cr)->getBBox().getMinY(), CountDN.getCounts());
        if (CountDN.getCounts() < minDn) { /* not bright enough */
            pexLogging::TTrace<11>("algorithms.CR", "Erasing CR");

            cr = CRs.erase(cr);
            --cr;                       // back up to previous CR (we're going to increment it)
            --end;
        }
    }
    ncr = CRs.size();           /* some may have been too faint */
/*
 * We've found them all, time to kill them all
 */
    bool const debias_values = true;
    bool grow = false;
    pexLogging::TTrace<2>("algorithms.CR", "Removing initial list of CRs");
    removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
#if 0                                   // Useful to see phase 2 in ds9; debugging only
    (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
                                   mimage.getMask()->getPlaneBitMask("DETECTED"));
#endif
/*
 * Now that we've removed them, go through image again, examining area around
 * each CR for extra bad pixels. Note that we set cond3Fac = 0 for this pass
 *
 * We iterate niteration times;  niter==1 was sufficient for SDSS data, but megacam
 * CCDs are different -- who knows for other devices?
 */
    bool too_many_crs = false;          // we've seen too many CR pixels
    int nextra = 0;                     // number of pixels added to list of CRs
    for (int i = 0; i != niteration && !too_many_crs; ++i) {
        pexLogging::TTrace<1>("algorithms.CR", "Starting iteration %d", i);
        for (std::vector<detection::Footprint::Ptr>::iterator fiter = CRs.begin();
             fiter != CRs.end(); fiter++) {
            detection::Footprint::Ptr cr = *fiter;
/*
 * Are all those `CR' pixels interpolated?  If so, don't grow it
 */
            {
                detection::Footprint::Ptr om = footprintAndMask(cr, mimage.getMask(), interpBit);
                int const npix = (om) ? om->getNpix() : 0;

                if (npix == cr->getNpix()) {
                    continue;
                }
            }
/*
 * No; some of the suspect pixels aren't interpolated
 */
            detection::Footprint extra;                     // extra pixels added to cr
            for (detection::Footprint::SpanList::const_iterator siter = cr->getSpans().begin();
                 siter != cr->getSpans().end(); siter++) {
                detection::Span::Ptr const span = *siter;

                /*
                 * Check the lines above and below the span.  We're going to check a 3x3 region around
                 * the pixels, so we need a buffer around the edge.  We check the pixels just to the
                 * left/right of the span, so the buffer needs to be 2 pixels (not just 1) in the
                 * column direction, but only 1 in the row direction.
                 */
                int const y = span->getY() - mimage.getY0();
                if (y < 2 || y >= nrow - 2) {
                    continue;
                }
                int x0 = span->getX0() - mimage.getX0();
                int x1 = span->getX1() - mimage.getX0();
                x0 = (x0 < 2) ? 2 : (x0 > ncol - 3) ? ncol - 3 : x0;
                x1 = (x1 < 2) ? 2 : (x1 > ncol - 3) ? ncol - 3 : x1;

                checkSpanForCRs(&extra, crpixels, y - 1, x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
                checkSpanForCRs(&extra, crpixels, y,     x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
                checkSpanForCRs(&extra, crpixels, y + 1, x0, x1, mimage,
                                minSigma/2, thresH, thresV, thresD, bkgd, 0, keep);
            }

            if (extra.getSpans().size() > 0) {      // we added some pixels
                if (nextra + static_cast<int>(crpixels.size()) > nCrPixelMax) {
                    too_many_crs = true;
                    break;
                }

                nextra += extra.getNpix();

                detection::Footprint::SpanList &espans = extra.getSpans();
                for (detection::Footprint::SpanList::const_iterator siter = espans.begin();
                     siter != espans.end(); siter++) {
                    cr->addSpan(**siter);
                }
                                cr->normalize();
            }            
        }

        if (nextra == 0) {
            break;
        }
    }
/*
 * mark those pixels as CRs
 */
    if (!too_many_crs) {
        (void)setMaskFromFootprintList(mimage.getMask().get(), CRs, crBit);
    }
/*
 * Maybe reinstate initial values; n.b. the same pixel may appear twice, so we want the
 * first value stored (hence the uses of rbegin/rend)
 *
 * We have to do this if we decide _not_ to remove certain CRs,
 * for example those which lie next to saturated pixels
 */
    if (keep || too_many_crs) {
        if (crpixels.size() > 0) {
            int const imageX0 = mimage.getX0();
            int const imageY0 = mimage.getY0();

            std::sort(crpixels.begin(), crpixels.end()); // sort into birth order
        
            crpixel_riter rend = crpixels.rend();
            for (crpixel_riter crp = crpixels.rbegin(); crp != rend; ++crp) {
                if (crp->row == -1)
                    // dummy; skip it.
                    continue;
                mimage.at(crp->col - imageX0, crp->row - imageY0).image() = crp->val;
            }
        }
    } else {
        if (true || nextra > 0) {
            grow = true;
            pexLogging::TTrace<2>("algorithms.CR", "Removing final list of CRs, grow = %d", grow);
            removeCR(mimage, CRs, bkgd, crBit, saturBit, badMask, debias_values, grow);
        }
/*
 * we interpolated over all CR pixels, so set the interp bits too
 */
        (void)setMaskFromFootprintList(mimage.getMask().get(), CRs,
                                       static_cast<MaskPixel>(crBit | interpBit));
    }

    if (too_many_crs) {                 // we've cleaned up, so we can throw the exception
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          (boost::format("Too many CR pixels (max %d)") % nCrPixelMax).str());
    }

    return CRs;
}

template<typename Image >

std::pair<std::vector<double>, lsst::afw::math::KernelList> lsst::meas::algorithms::fitKernelParamsToImage	(	afwMath::LinearCombinationKernel const &	kernel,
		Image const &	image,
		afwGeom::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

kernel	the Kernel to fit
image	the image to be fit
pos	the position of the object

Definition at line 1061 of file SpatialModelPsf.cc.

{
    typedef afwImage::Image<afwMath::Kernel::Pixel> KernelT;

    afwMath::KernelList kernels = kernel.getKernelList();         // the Kernels that kernel adds together
    int const nKernel = kernels.size();

    if (nKernel == 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          "Your kernel must have at least one component");
    }

    /*
     * Go through all the kernels, get a copy centered at the desired sub-pixel position, and then
     * extract a subImage from the parent image at the same place
     */
    std::vector<KernelT::Ptr> kernelImages = offsetKernel<KernelT>(kernel, pos[0], pos[1]);
    afwGeom::BoxI bbox(kernelImages[0]->getBBox());
    Image const& subImage(Image(image, bbox, afwImage::PARENT, false)); // shallow copy

    /*
     * Solve the linear problem  subImage = sum x_i K_i + epsilon; we solve this for x_i by constructing the
     * normal equations, A x = b
     */
    Eigen::MatrixXd A(nKernel, nKernel);
    Eigen::VectorXd b(nKernel);

    for (int i = 0; i != nKernel; ++i) {
        b(i) = afwImage::innerProduct(*kernelImages[i], *subImage.getImage());

        for (int j = i; j != nKernel; ++j) {
            A(i, j) = A(j, i) = afwImage::innerProduct(*kernelImages[i], *kernelImages[j]);
        }
    }
    Eigen::VectorXd x(nKernel);

    if (nKernel == 1) {
        x(0) = b(0)/A(0, 0);
    } else {
        x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
    }

    // the XY0() point of the shifted Kernel basis functions
    int const x0 = kernelImages[0]->getX0(), y0 = kernelImages[0]->getY0(); 

    afwMath::KernelList newKernels(nKernel);
    std::vector<double> params(nKernel);
    for (int i = 0; i != nKernel; ++i) {
        afwMath::Kernel::Ptr newKernel(new afwMath::FixedKernel(*kernelImages[i]));
        newKernel->setCtrX(x0 + static_cast<int>(newKernel->getWidth()/2));
        newKernel->setCtrY(y0 + static_cast<int>(newKernel->getHeight()/2));

        params[i] = x[i];
        newKernels[i] = newKernel;
    }

    return std::make_pair(params, newKernels);
}

template<typename Image >

std::pair<std::vector<double>, afwMath::KernelList> lsst::meas::algorithms::fitKernelParamsToImage	(	afwMath::LinearCombinationKernel const &	kernel,
		Image const &	image,
		afwGeom::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

kernel	the Kernel to fit
image	the image to be fit
pos	the position of the object

Definition at line 1061 of file SpatialModelPsf.cc.

{
    typedef afwImage::Image<afwMath::Kernel::Pixel> KernelT;

    afwMath::KernelList kernels = kernel.getKernelList();         // the Kernels that kernel adds together
    int const nKernel = kernels.size();

    if (nKernel == 0) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                          "Your kernel must have at least one component");
    }

    /*
     * Go through all the kernels, get a copy centered at the desired sub-pixel position, and then
     * extract a subImage from the parent image at the same place
     */
    std::vector<KernelT::Ptr> kernelImages = offsetKernel<KernelT>(kernel, pos[0], pos[1]);
    afwGeom::BoxI bbox(kernelImages[0]->getBBox());
    Image const& subImage(Image(image, bbox, afwImage::PARENT, false)); // shallow copy

    /*
     * Solve the linear problem  subImage = sum x_i K_i + epsilon; we solve this for x_i by constructing the
     * normal equations, A x = b
     */
    Eigen::MatrixXd A(nKernel, nKernel);
    Eigen::VectorXd b(nKernel);

    for (int i = 0; i != nKernel; ++i) {
        b(i) = afwImage::innerProduct(*kernelImages[i], *subImage.getImage());

        for (int j = i; j != nKernel; ++j) {
            A(i, j) = A(j, i) = afwImage::innerProduct(*kernelImages[i], *kernelImages[j]);
        }
    }
    Eigen::VectorXd x(nKernel);

    if (nKernel == 1) {
        x(0) = b(0)/A(0, 0);
    } else {
        x = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
    }

    // the XY0() point of the shifted Kernel basis functions
    int const x0 = kernelImages[0]->getX0(), y0 = kernelImages[0]->getY0(); 

    afwMath::KernelList newKernels(nKernel);
    std::vector<double> params(nKernel);
    for (int i = 0; i != nKernel; ++i) {
        afwMath::Kernel::Ptr newKernel(new afwMath::FixedKernel(*kernelImages[i]));
        newKernel->setCtrX(x0 + static_cast<int>(newKernel->getWidth()/2));
        newKernel->setCtrY(y0 + static_cast<int>(newKernel->getHeight()/2));

        params[i] = x[i];
        newKernels[i] = newKernel;
    }

    return std::make_pair(params, newKernels);
}

template<typename Image >

std::pair<lsst::afw::math::Kernel::Ptr, std::pair<double, double> > lsst::meas::algorithms::fitKernelToImage	(	afwMath::LinearCombinationKernel const &	kernel,
		Image const &	image,
		afwGeom::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

kernel	the Kernel to fit
image	the image to be fit
pos	the position of the object

Definition at line 1133 of file SpatialModelPsf.cc.

{
    std::pair<std::vector<double>, afwMath::KernelList> const fit = 
        fitKernelParamsToImage(kernel, image, pos);
    std::vector<double> params = fit.first;
    afwMath::KernelList kernels = fit.second;
    int const nKernel = params.size();
    assert(kernels.size() == static_cast<unsigned int>(nKernel));

    double amp = 0.0;
    for (int i = 0; i != nKernel; ++i) {
        afwMath::Kernel::Ptr base = kernels[i];
        afwMath::FixedKernel::Ptr k = boost::static_pointer_cast<afwMath::FixedKernel>(base);
        amp += params[i] * k->getSum();
    }

    afwMath::Kernel::Ptr outputKernel(new afwMath::LinearCombinationKernel(kernels, params));
    double chisq = 0.0;
    outputKernel->setCtrX(kernels[0]->getCtrX());
    outputKernel->setCtrY(kernels[0]->getCtrY());

    return std::make_pair(outputKernel, std::make_pair(amp, chisq));
}

template<typename Image >

std::pair<afwMath::Kernel::Ptr, std::pair<double, double> > lsst::meas::algorithms::fitKernelToImage	(	afwMath::LinearCombinationKernel const &	kernel,
		Image const &	image,
		afwGeom::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

kernel	the Kernel to fit
image	the image to be fit
pos	the position of the object

Definition at line 1133 of file SpatialModelPsf.cc.

{
    std::pair<std::vector<double>, afwMath::KernelList> const fit = 
        fitKernelParamsToImage(kernel, image, pos);
    std::vector<double> params = fit.first;
    afwMath::KernelList kernels = fit.second;
    int const nKernel = params.size();
    assert(kernels.size() == static_cast<unsigned int>(nKernel));

    double amp = 0.0;
    for (int i = 0; i != nKernel; ++i) {
        afwMath::Kernel::Ptr base = kernels[i];
        afwMath::FixedKernel::Ptr k = boost::static_pointer_cast<afwMath::FixedKernel>(base);
        amp += params[i] * k->getSum();
    }

    afwMath::Kernel::Ptr outputKernel(new afwMath::LinearCombinationKernel(kernels, params));
    double chisq = 0.0;
    outputKernel->setCtrX(kernels[0]->getCtrX());
    outputKernel->setCtrY(kernels[0]->getCtrY());

    return std::make_pair(outputKernel, std::make_pair(amp, chisq));
}

template<typename PixelT >

std::pair< bool, double > lsst::meas::algorithms::fitSpatialKernelFromPsfCandidates	(	afwMath::Kernel *	kernel,
		afwMath::SpatialCellSet const &	psfCells,
		int const	nStarPerCell,
		double const	tolerance,
		double const	lambda
	)

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

Parameters

kernel	the Kernel to fit
psfCells	A SpatialCellSet containing PsfCandidates
nStarPerCell	max no. of stars per cell; <= 0 => infty
tolerance	Tolerance; how close chi^2 should be to true minimum
lambda	floor for variance is lambda*data

Definition at line 621 of file SpatialModelPsf.cc.

                                   {
    typedef typename afwImage::Image<PixelT> Image;

    int const nComponents = kernel->getNKernelParameters();
    int const nSpatialParams = kernel->getNSpatialParameters();
    //
    // visitor that evaluates the chi^2 of the current fit
    //
    evalChi2Visitor<PixelT> getChi2(*kernel, lambda);
    //
    // We have to unpack the Kernel coefficients into a linear array, coeffs
    //
    std::vector<double> coeffs;         // The coefficients we want to fit
    coeffs.assign(nComponents*nSpatialParams, 0.0);

    std::vector<double> stepSize;       // step sizes
    stepSize.assign(nComponents*nSpatialParams, 100);
    //
    // Translate that into minuit's language
    //
    ROOT::Minuit2::MnUserParameters fitPar;
    std::vector<std::string> paramNames;
    paramNames.reserve(nComponents*nSpatialParams);
    
    for (int i = 0, c = 0; c != nComponents; ++c) {
        coeffs[i] = 1;                  // the constant part of each spatial order
        for (int s = 0; s != nSpatialParams; ++s, ++i) {
            paramNames.push_back((boost::format("C%d:%d") % c % s).str());
            fitPar.Add(paramNames[i].c_str(), coeffs[i], stepSize[i]);
        }
    }
    fitPar.Fix("C0:0");
    //
    // Create the minuit object that knows how to minimise our functor
    //
    MinimizeChi2<PixelT> minimizerFunc(getChi2, kernel, psfCells, nStarPerCell, nComponents, nSpatialParams);

    double const errorDef = 1.0;       // use +- 1sigma errors
    minimizerFunc.setErrorDef(errorDef);
    //
    // tell minuit about it
    //    
    ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
    //
    // And let it loose
    //
    int maxFnCalls = 0;                 // i.e. unlimited
    ROOT::Minuit2::FunctionMinimum min =
        migrad(maxFnCalls, tolerance/(1e-4*errorDef)); // minuit uses 0.1*1e-3*tolerance*errorDef

    float minChi2 = min.Fval();
    bool const isValid = min.IsValid() && std::isfinite(minChi2);
    
    if (true || isValid) {              // calculate coeffs even in minuit is unhappy
        for (int i = 0; i != nComponents*nSpatialParams; ++i) {
            coeffs[i] = min.UserState().Value(i);
        }

        setSpatialParameters(kernel, coeffs);
    }

#if 0                                   // Estimate errors;  we don't really need this
    ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
    for (int i = 0, c = 0; c != nComponents; ++c) {
        for (int s = 0; s != nSpatialParams; ++s, ++i) {
            char const *name = paramNames[i].c_str();
            printf("%s %g", name, min.UserState().Value(name));
            if (isValid && !fitPar.Parameter(fitPar.Index(name)).IsFixed()) {
                printf(" (%g+%g)\n", minos(i).first, minos(i).second);
            }
            printf("\n");
        }
    }
#endif
    //
    // One time more through the Candidates setting their chi^2 values. We'll
    // do all the candidates this time, not just the first nStarPerCell
    //
    psfCells.visitAllCandidates(&getChi2, true);
    
    return std::make_pair(isValid, minChi2);
}

template<typename PixelT >

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

Parameters

kernel	the Kernel to fit
psfCells	A SpatialCellSet containing PsfCandidates
doNonLinearFit	Use the full-up nonlinear fitter
nStarPerCell	max no. of stars per cell; <= 0 => infty
tolerance	Tolerance; how close chi^2 should be to true minimum
lambda	floor for variance is lambda*data

Definition at line 909 of file SpatialModelPsf.cc.

{
    if (doNonLinearFit) {
        return fitSpatialKernelFromPsfCandidates<PixelT>(kernel, psfCells, nStarPerCell, tolerance);
    }

    double const tau = 0;               // softening for errors

    afwMath::LinearCombinationKernel const* lcKernel =
        dynamic_cast<afwMath::LinearCombinationKernel const*>(kernel);
    if (!lcKernel) {
        throw LSST_EXCEPT(lsst::pex::exceptions::LogicError,
                          "Failed to cast Kernel to LinearCombinationKernel while building spatial PSF model");
    }
#if 1
    //
    // Set the initial amplitudes of all our candidates
    //
    setAmplitudeVisitor<PixelT> setAmplitude;
    psfCells.visitAllCandidates(&setAmplitude, true);
#endif
    //
    // visitor that fills out the A and b matrices (we'll solve A x = b for the coeffs, x)
    //
    FillABVisitor<PixelT> getAB(*lcKernel, tau);
    //
    // Actually visit all our candidates
    //
    psfCells.visitCandidates(&getAB, nStarPerCell, true);
    //
    // Extract A and b, and solve Ax = b
    //
    Eigen::MatrixXd const& A = getAB.getA();
    Eigen::VectorXd const& b = getAB.getB();
    Eigen::VectorXd x0(b.size());       // Solution to matrix problem
    x0 = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
#if 0
    std::cout << "A " << A << std::endl;
    std::cout << "b " << b.transpose() << std::endl;
    std::cout << "x " << x.transpose() << std::endl;

    afwImage::Image<double> img(b.size(), b.size());
    for (int j = 0; j < b.size(); ++j) {
        for (int i = 0; i < b.size(); ++i) {
            img(i, j) = A(i, j);
        }
    }
    img.writeFits("a.fits");

    if (x.cols() >= 6) {
        for (int i = 0; i != 6; ++i) {
            double xcen = 25; double ycen = 35 + 35*i;
            std::cout << "x, y " << xcen << " , " << ycen << " b "
                      << (x[3] + xcen*x[4] + ycen*x[5])/(x[0] + xcen*x[1] + ycen*x[2]) << std::endl;
        }
    }
#endif

    // Generate kernel parameters (including 0th component) from matrix solution
    Eigen::VectorXd x(kernel->getNKernelParameters() * kernel->getNSpatialParameters()); // Kernel parameters
    x(0) = 1.0;
    std::fill(x.data() + 1, x.data() + kernel->getNSpatialParameters(), 0.0);
    std::copy(x0.data(), x0.data() + x0.size(), x.data() + kernel->getNSpatialParameters());

    setSpatialParameters(kernel, x);
    //
    // One time more through the Candidates setting their chi^2 values. We'll
    // do all the candidates this time, not just the first nStarPerCell
    //
    // visitor that evaluates the chi^2 of the current fit
    //
    evalChi2Visitor<PixelT> getChi2(*kernel, lambda);

    psfCells.visitAllCandidates(&getChi2, true);
    
    return std::make_pair(true, getChi2.getValue());
}

template<class ShapeletPtr >

int lsst::meas::algorithms::getImageSize ( ShapeletPtr  shapelet, int  size  )

inline

Definition at line 106 of file ShapeletKernel.cc.

    { return size == 0 ? int(ceil(shapelet->getSigma()*10.)) : size; }

Eigen::Matrix2d lsst::meas::algorithms::getJacobian	(	const lsst::afw::image::Wcs &	wcs,
		const lsst::afw::geom::PointD &	pos
	)

a helper function to deal with the fact that Wcs doesn't directly return a Jacobian.

This function calculates the local Jacobian of the distortion pattern at a given location.

The input position is given in chip coordinates, since that's what we generally have when we want to calculate the local distortion.

The return matrix is:

J = ( du/dx du/dy ) ( dv/dx dv/dy )

where (u,v) are the sky coordinates, and (x,y) are the chip coordinates.

Parameters

wcs	The wcs defining the distortion function
pos	The position in chip coordinates at which to find the Jacobian

Definition at line 209 of file Shapelet.cc.

    {
        lsst::afw::geom::AffineTransform localTransform = 
            wcs.linearizePixelToSky(pos,lsst::afw::geom::degrees);

        // J = ( du/dx  du/dy )
        //     ( dv/dx  dv/dy )
        // where (u,v) are sky coordinates, and (x,y) are chip coordinates.
        Eigen::Matrix2d J;
        // Answer comes out in degrees for u,v.  *3600 to get arcsec.
        J(0,0) = localTransform.getMatrix()(0,0) * 3600.;
        J(0,1) = localTransform.getMatrix()(0,1) * 3600.;
        J(1,0) = localTransform.getMatrix()(1,0) * 3600.;
        J(1,1) = localTransform.getMatrix()(1,1) * 3600.;
        return J;
    }

afw::geom::Box2I lsst::meas::algorithms::getOverallBBox

(

std::vector< boost::shared_ptr< afw::image::Image< double > >> const &

imgVector

)

Definition at line 183 of file CoaddPsf.cc.

                                                                                     {

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

template<typename MaskedImageT >

void lsst::meas::algorithms::interpolateOverDefects	(	MaskedImageT &	mimage,
		lsst::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

mimage	Image to patch
psf	the Image's PSF
_badList	List of Defects to patch
fallbackValue	Value to fallback to if all else fails
useFallbackValueAtEdge	Use the fallback value at the image's edge?

Definition at line 2058 of file Interp.cc.

                             {
/*
 * Allow for image's origin
 */
    int const width = mimage.getWidth();
    int const height = mimage.getHeight();

    std::vector<Defect::Ptr> badList;
    badList.reserve(_badList.size());
    for (std::vector<Defect::Ptr>::iterator ptr = _badList.begin(), end = _badList.end(); ptr != end; ++ptr) {
        geom::BoxI bbox = (*ptr)->getBBox();
        bbox.shift(geom::ExtentI(-mimage.getX0(), -mimage.getY0())); //allow for image's origin
                geom::PointI min = bbox.getMin(), max = bbox.getMax();
                if(min.getX() >= width){
            continue;
        } else if (min.getX() < 0) {
            if (max.getX() < 0) {
                continue;
            } else {
                min.setX(0);
            }
        }

        if (max.getX() < 0) {
            continue;
        } else if (max.getX() >= width) {
            max.setX(width - 1);
        }

        bbox = geom::BoxI(min, max);
        Defect::Ptr ndefect(new Defect(bbox));
        ndefect->classify((*ptr)->getPos(), (*ptr)->getType());
        badList.push_back(ndefect);
    }

    sort(badList.begin(), badList.end(), Sort_ByX0<Defect>());
/*
 * Go through the frame looking at each pixel (except the edge ones which we ignore)
 */
    typename MaskedImageT::Mask::Pixel const interpBit =
        mimage.getMask()->getPlaneBitMask("INTRP"); // interp'd pixels

    int nUseInterp = 6;                       // no. of pixels to interpolate towards edge
    assert(nUseInterp < Defect::WIDE_DEFECT); // we'd use C++11's static_assert if available

    for (int y = 0; y != height; y++) {
        std::vector<Defect::Ptr> badList1D = classify_defects(badList, y, width);

        do_defects(badList1D, y, *mimage.getImage(),
                   -std::numeric_limits<typename MaskedImageT::Image::Pixel>::max(),
                   fallbackValue, useFallbackValueAtEdge, nUseInterp);

        do_defects(badList1D, y, *mimage.getMask(), interpBit, useFallbackValueAtEdge, nUseInterp);

        do_defects(badList1D, y, *mimage.getVariance(),
                   -std::numeric_limits<typename MaskedImageT::Image::Pixel>::max(),
                   fallbackValue, useFallbackValueAtEdge, nUseInterp);
    }
}

template<typename ExposureT >

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

Factory function for ExposurePatch.

Definition at line 86 of file ExposurePatch.h.

      {
    return boost::make_shared<ExposurePatch<ExposureT> >(exp, foot, center);
}

template<typename ExposureT >

boost::shared_ptr< ExposurePatch<ExposureT> > lsst::meas::algorithms::makeExposurePatch	(	boost::shared_ptr< ExposureT const >	exp,
		afw::detection::Footprint const &	standardFoot,
		afw::geom::Point2D const &	standardCenter,
		afw::image::Wcs const &	standardWcs
	)

Definition at line 94 of file ExposurePatch.h.

      {
    return boost::make_shared<ExposurePatch<ExposureT> >(exp, standardFoot, standardCenter, standardWcs);
}

template<typename PixelT >

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

Return a PsfCandidate of the right sort

Cf. std::make_pair

Parameters

source	The detected Source
image	The image wherein lies the object

Definition at line 188 of file PsfCandidate.h.

    {
        return boost::make_shared< PsfCandidate<PixelT> >(source, image);
    }

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

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

Definition at line 520 of file SpatialModelPsf.cc.

{
    int const nComponents = kernel->getNKernelParameters();
    int const nSpatialParams = kernel->getNSpatialParameters();

    assert (nComponents*nSpatialParams == static_cast<long>(coeffs.size()));

    std::vector<std::vector<double> > kCoeffs; // coefficients rearranged for Kernel
    kCoeffs.reserve(nComponents);
    for (int i = 0; i != nComponents; ++i) {
        kCoeffs.push_back(std::vector<double>(nSpatialParams));
        std::copy(coeffs.begin() + i*nSpatialParams,
                  coeffs.begin() + (i + 1)*nSpatialParams, kCoeffs[i].begin());
    }
    
    kernel->setSpatialParameters(kCoeffs);
}

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

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

Definition at line 544 of file SpatialModelPsf.cc.

{
    int const nComponents = kernel->getNKernelParameters();
    int const nSpatialParams = kernel->getNSpatialParameters();

    assert (nComponents*nSpatialParams == vec.size());

    std::vector<std::vector<double> > kCoeffs; // coefficients rearranged for Kernel
    kCoeffs.reserve(nComponents);
    for (int i = 0; i != nComponents; ++i) {
        std::vector<double> spatialCoeffs(nSpatialParams);
        for (int j = 0; j != nSpatialParams; ++j) {
            spatialCoeffs[j] = vec[i*nSpatialParams + j];
        }
        kCoeffs.push_back(spatialCoeffs);
    }
    
    kernel->setSpatialParameters(kCoeffs);
}

template<typename ImageT >

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

template<typename MaskedImageT >

double lsst::meas::algorithms::subtractPsf	(	afwDetection::Psf const &	psf,
		MaskedImageT *	data,
		double	x,
		double	y,
		double	psfFlux
	)

Subtract a PSF from an image at a given position

Parameters

psf	the PSF to subtract
data	Image to subtract PSF from
x	column position
y	row position
psfFlux	object's PSF flux (if not NaN)

Definition at line 999 of file SpatialModelPsf.cc.

{
    if (lsst::utils::isnan(x + y)) {
        return std::numeric_limits<double>::quiet_NaN();
    }

    //
    // Get Psf candidate
    //
    afwDetection::Psf::Image::Ptr kImage = psf.computeImage(afwGeom::PointD(x, y));

    //
    // Now find the proper sub-Image
    //
    afwGeom::BoxI bbox = kImage->getBBox();
    
    typename MaskedImageT::Ptr subData(new MaskedImageT(*data, bbox, afwImage::PARENT, false)); // shallow copy
    //
    // Now we've got both; find the PSF's amplitude
    //
    double lambda = 0.0;                // floor for variance is lambda*data
    try {
        double chi2;                    // chi^2 for fit
        double amp;                     // estimate of amplitude of model at this point

        if (lsst::utils::isnan(psfFlux)) {
            std::pair<double, double> result = fitKernel(*kImage, *subData, lambda, true);
            chi2 = result.first;        // chi^2 for fit
            amp = result.second;        // estimate of amplitude of model at this point
        } else {
            chi2 = std::numeric_limits<double>::quiet_NaN();
            amp = psfFlux/afwMath::makeStatistics(*kImage, afwMath::SUM).getValue();
        }
        //
        // Convert kImage to the proper type so that I can subtract it.
        //
        typename MaskedImageT::Image::Ptr
            kImageF(new typename MaskedImageT::Image(*kImage, true)); // of data's type

        *kImageF *= amp;
        *subData->getImage() -= *kImageF;
        
        return chi2;
    } catch(lsst::pex::exceptions::RangeError &e) {
        LSST_EXCEPT_ADD(e, (boost::format("Object at (%.2f, %.2f)") % x % y).str());
        throw e;
    }
}