lsst::meas::base::detail Namespace Reference

Classes
class	SdssShapeImpl

Functions
template<typename ImageT >
bool	getAdaptiveMoments (ImageT const &mimage, double bkgd, double xcen, double ycen, double shiftmax, SdssShapeImpl *shape, int maxIter, float tol1, float tol2)
template<typename ImageT >
std::pair< double, double >	getFixedMomentsFlux (ImageT const &image, double bkgd, double xcen, double ycen, SdssShapeImpl const &shape_)
	Return the flux of an object, using the aperture described by the SdssShape object. More...

Variables
int const	SDSS_SHAPE_MAX_ITER = 100
float const	SDSS_SHAPE_TOL1 = 1.0e-5
float const	SDSS_SHAPE_TOL2 = 1.0e-4

Function Documentation

template<typename ImageT >

bool lsst::meas::base::detail::getAdaptiveMoments	(	ImageT const &	mimage,
		double	bkgd,
		double	xcen,
		double	ycen,
		double	shiftmax,
		SdssShapeImpl *	shape,
		int	maxIter,
		float	tol1,
		float	tol2
	)

Workhorse for adaptive moments

Calculate adaptive moments from an image

The moments are measured iteratively with a Gaussian window with width equal to the second moments from the previous iteration.

Parameters

mimage	the data to process
bkgd	background level
xcen	x-centre of object
ycen	y-centre of object
shiftmax	max allowed centroid shift
shape	a place to store desired data
maxIter	Maximum number of iterations
tol1	Convergence tolerance for e1,e2
tol2	Convergence tolerance for FWHM

Definition at line 448 of file SdssShape.cc.

{
    double I0 = 0;                      // amplitude of best-fit Gaussian
    double sum;                         // sum of intensity*weight
    double sumx, sumy;                  // sum ((int)[xy])*intensity*weight
    double sumxx, sumxy, sumyy;         // sum {x^2,xy,y^2}*intensity*weight
    double sums4;                       // sum intensity*weight*exponent^2
    float const xcen0 = xcen;           // initial centre
    float const ycen0 = ycen;           //                of object

    double sigma11W = 1.5;              // quadratic moments of the
    double sigma12W = 0.0;              //     weighting fcn;
    double sigma22W = 1.5;              //               xx, xy, and yy

    double w11 = -1, w12 = -1, w22 = -1;        // current weights for moments; always set when iter == 0
    float e1_old = 1e6, e2_old = 1e6;           // old values of shape parameters e1 and e2
    float sigma11_ow_old = 1e6;                 // previous version of sigma11_ow
    
    typename ImageAdaptor<ImageT>::Image const &image = ImageAdaptor<ImageT>().getImage(mimage);

    if (lsst::utils::isnan(xcen) || lsst::utils::isnan(ycen)) {
        // Can't do anything
        shape->setFlag(SdssShapeImpl::UNWEIGHTED_BAD);
        return false;
    }

    bool interpflag = false;            // interpolate finer than a pixel?
    lsst::afw::geom::BoxI bbox;
    int iter = 0;                       // iteration number
    for (; iter < maxIter; iter++) {
        bbox = set_amom_bbox(image.getWidth(), image.getHeight(), xcen, ycen, sigma11W, sigma12W, sigma22W);
        boost::tuple<std::pair<bool, double>, double, double, double> weights = 
            getWeights(sigma11W, sigma12W, sigma22W);
        if (!weights.get<0>().first) {
            shape->setFlag(SdssShapeImpl::UNWEIGHTED);
            break;
        }

        double const detW = weights.get<0>().second;
        
#if 0                                   // this form was numerically unstable on my G4 powerbook
        assert(detW >= 0.0);
#else
        assert(sigma11W*sigma22W >= sigma12W*sigma12W - std::numeric_limits<float>::epsilon());
#endif

        {
            const double ow11 = w11;    // old
            const double ow12 = w12;    //     values
            const double ow22 = w22;    //            of w11, w12, w22

            w11 = weights.get<1>();
            w12 = weights.get<2>();
            w22 = weights.get<3>();

            if (shouldInterp(sigma11W, sigma22W, detW)) {
                if (!interpflag) {
                    interpflag = true;       // N.b.: stays set for this object
                    if (iter > 0) {
                        sigma11_ow_old = 1.e6; // force at least one more iteration
                        w11 = ow11;
                        w12 = ow12;
                        w22 = ow22;
                        iter--;         // we didn't update wXX
                    }
                }
            }
        }

        if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22,
                           &I0, &sum, &sumx, &sumy, &sumxx, &sumxy, &sumyy, &sums4) < 0) {
            shape->setFlag(SdssShapeImpl::UNWEIGHTED);
            break;
        }
#if 0
/*
 * Find new centre
 *
 * This is only needed if we update the centre; if we use the input position we've already done the work
 */
        xcen = sumx/sum;
        ycen = sumy/sum;
#endif
        shape->setX(sumx/sum); // update centroid.  N.b. we're not setting errors here
        shape->setY(sumy/sum);

        if (fabs(shape->getX() - xcen0) > shiftmax || fabs(shape->getY() - ycen0) > shiftmax) {
            shape->setFlag(SdssShapeImpl::SHIFT);
        }
/*
 * OK, we have the centre. Proceed to find the second moments.
 */
        float const sigma11_ow = sumxx/sum; // quadratic moments of
        float const sigma22_ow = sumyy/sum; //          weight*object
        float const sigma12_ow = sumxy/sum; //                 xx, xy, and yy 

        if (sigma11_ow <= 0 || sigma22_ow <= 0) {
            shape->setFlag(SdssShapeImpl::UNWEIGHTED);
            break;
        }

        float const d = sigma11_ow + sigma22_ow; // current values of shape parameters
        float const e1 = (sigma11_ow - sigma22_ow)/d;
        float const e2 = 2.0*sigma12_ow/d;
/*
 * Did we converge?
 */
        if (iter > 0 &&
           fabs(e1 - e1_old) < tol1 && fabs(e2 - e2_old) < tol1 &&
           fabs(sigma11_ow/sigma11_ow_old - 1.0) < tol2 ) {
            break;                              // yes; we converged
        }

        e1_old = e1;
        e2_old = e2;
        sigma11_ow_old = sigma11_ow;
/*
 * Didn't converge, calculate new values for weighting function
 *
 * The product of two Gaussians is a Gaussian:
 * <x^2 exp(-a x^2 - 2bxy - cy^2) exp(-Ax^2 - 2Bxy - Cy^2)> = 
 *                            <x^2 exp(-(a + A) x^2 - 2(b + B)xy - (c + C)y^2)>
 * i.e. the inverses of the covariances matrices add.
 *
 * We know sigmaXX_ow and sigmaXXW, the covariances of the weighted object
 * and of the weights themselves.  We can estimate the object's covariance as
 *   sigmaXX_ow^-1 - sigmaXXW^-1
 * and, as we want to find a set of weights with the _same_ covariance as the
 * object we take this to be the an estimate of our correct weights.
 *
 * N.b. This assumes that the object is roughly Gaussian.
 * Consider the object:
 *   O == delta(x + p) + delta(x - p)
 * the covariance of the weighted object is equal to that of the unweighted
 * object, and this prescription fails badly.  If we detect this, we set
 * the UNWEIGHTED flag, and calculate the UNweighted moments
 * instead.
 */
        {
            float n11, n12, n22;                // elements of inverse of next guess at weighting function
            float ow11, ow12, ow22;             // elements of inverse of sigmaXX_ow

            boost::tuple<std::pair<bool, double>, double, double, double> weights = 
                getWeights(sigma11_ow, sigma12_ow, sigma22_ow);
            if (!weights.get<0>().first) {
                shape->setFlag(SdssShapeImpl::UNWEIGHTED);
                break;
            }
         
            ow11 = weights.get<1>();
            ow12 = weights.get<2>();
            ow22 = weights.get<3>();

            n11 = ow11 - w11;
            n12 = ow12 - w12;
            n22 = ow22 - w22;

            weights = getWeights(n11, n12, n22, false);
            if (!weights.get<0>().first) {
                // product-of-Gaussians assumption failed
                shape->setFlag(SdssShapeImpl::UNWEIGHTED);
                break;
            }
      
            sigma11W = weights.get<1>();
            sigma12W = weights.get<2>();
            sigma22W = weights.get<3>();
        }

        if (sigma11W <= 0 || sigma22W <= 0) {
            shape->setFlag(SdssShapeImpl::UNWEIGHTED);
            break;
        }
    }

    if (iter == maxIter) {
        shape->setFlag(SdssShapeImpl::UNWEIGHTED);
        shape->setFlag(SdssShapeImpl::MAXITER);
    }

    if (sumxx + sumyy == 0.0) {
        shape->setFlag(SdssShapeImpl::UNWEIGHTED);
    }
/*
 * Problems; try calculating the un-weighted moments
 */
    if (shape->getFlag(SdssShapeImpl::UNWEIGHTED)) {
        w11 = w22 = w12 = 0;
        if (calcmom<false>(image, xcen, ycen, bbox, bkgd, interpflag, w11, w12, w22,
                           &I0, &sum, &sumx, &sumy, &sumxx, &sumxy, &sumyy, NULL) < 0 || sum <= 0) {
            shape->resetFlag(SdssShapeImpl::UNWEIGHTED);
            shape->setFlag(SdssShapeImpl::UNWEIGHTED_BAD);

            if (sum > 0) {
                shape->setIxx(1/12.0);      // a single pixel
                shape->setIxy(0.0);
                shape->setIyy(1/12.0);
            }
            
            return false;
        }

        sigma11W = sumxx/sum;          // estimate of object moments
        sigma12W = sumxy/sum;          //   usually, object == weight
        sigma22W = sumyy/sum;          //      at this point
    }

    shape->setI0(I0);
    shape->setIxx(sigma11W);
    shape->setIxy(sigma12W);
    shape->setIyy(sigma22W);
    shape->setIxy4(sums4/sum);

    if (shape->getIxx() + shape->getIyy() != 0.0) {
        int const ix = lsst::afw::image::positionToIndex(xcen);
        int const iy = lsst::afw::image::positionToIndex(ycen);
        
        if (ix >= 0 && ix < mimage.getWidth() && iy >= 0 && iy < mimage.getHeight()) {
            float const bkgd_var =
                ImageAdaptor<ImageT>().getVariance(mimage, ix, iy); // XXX Overestimate as it includes object

            if (bkgd_var > 0.0) {                                   // NaN is not > 0.0
                if (!(shape->getFlag(SdssShapeImpl::UNWEIGHTED))) {
                    SdssShapeImpl::Matrix4 fisher = calc_fisher(*shape, bkgd_var); // Fisher matrix 
                    shape->setCovar(fisher.inverse());
                }
            }
        }
    }

    return true;
}

template<typename ImageT >

std::pair< double, double > lsst::meas::base::detail::getFixedMomentsFlux	(	ImageT const &	image,
		double	bkgd,
		double	xcen,
		double	ycen,
		SdssShapeImpl const &	shape_
	)

Return the flux of an object, using the aperture described by the SdssShape object.

The SdssShape algorithm calculates an elliptical Gaussian fit to an object, so the "aperture" is an elliptical Gaussian

Returns

A std::pair of the flux and its error

Parameters

image	the data to process
bkgd	background level
xcen	x-centre of object (PARENT coordinates)
ycen	y-centre of object (PARENT coordinates)
shape_	The SdssShape of the object

Definition at line 692 of file SdssShape.cc.

{
    SdssShapeImpl shape = shape_; // we need a mutable copy

    afwGeom::BoxI const& bbox = set_amom_bbox(image.getWidth(), image.getHeight(), xcen, ycen,
                                              shape.getIxx(), shape.getIxy(), shape.getIyy());

    boost::tuple<std::pair<bool, double>, double, double, double> weights =
        getWeights(shape.getIxx(), shape.getIxy(), shape.getIyy());
    double const NaN = std::numeric_limits<double>::quiet_NaN();
    if (!weights.get<0>().first) {
        return std::make_pair(NaN, NaN);
    }

    double const w11 = weights.get<1>();
    double const w12 = weights.get<2>();
    double const w22 = weights.get<3>();
    bool const interp = shouldInterp(shape.getIxx(), shape.getIyy(), weights.get<0>().second);

    double I0 = 0;                      // amplitude of Gaussian
    calcmom<true>(ImageAdaptor<ImageT>().getImage(image), xcen - image.getX0(), ycen - image.getY0(),
                  bbox, bkgd, interp, w11, w12, w22,
                  &I0, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
    /*
     * We have enerything we need, but it isn't quite packaged right; we need an initialised SdssShapeImpl
     */
    shape.setI0(I0);

    {
        int ix = static_cast<int>(xcen - image.getX0());
        int iy = static_cast<int>(ycen - image.getY0());
        float bkgd_var = 
            ImageAdaptor<ImageT>().getVariance(image, ix, iy); // XXX Overestimate as it includes object

        SdssShapeImpl::Matrix4 fisher = calc_fisher(shape, bkgd_var); // Fisher matrix 

        shape.setCovar(fisher.inverse());
    }
    
    double const scale = shape.getFluxScale();
    return std::make_pair(shape.getI0()*scale, shape.getI0Err()*scale);
}

Variable Documentation

int const lsst::meas::base::detail::SDSS_SHAPE_MAX_ITER = 100

Definition at line 12 of file SdssShapeImpl.h.

float const lsst::meas::base::detail::SDSS_SHAPE_TOL1 = 1.0e-5

Definition at line 13 of file SdssShapeImpl.h.

float const lsst::meas::base::detail::SDSS_SHAPE_TOL2 = 1.0e-4

Definition at line 14 of file SdssShapeImpl.h.