Example of SpatialCellSet

Demonstrate the use of SpatialCellSets; the code's in spatialCellExample.cc.

Start by including needed headers, and declaring namespace aliases and a routine readImage

/* 
 * LSST Data Management System
 * Copyright 2008, 2009, 2010 LSST Corporation.
 * 
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the LSST License Statement and 
 * the GNU General Public License along with this program.  If not, 
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */
 
/*
 * This C++ example does the same thing as SpatialCellExample.py.  The latter of the python version
 * is that you can set display == True and see what's going on
 */
#include <string>
#include "lsst/utils/Utils.h"
#include "lsst/pex/exceptions.h"
#include "lsst/daf/base.h"
#include "lsst/afw/detection.h"
#include "lsst/afw/image.h"
#include "lsst/afw/math.h"
#include "lsst/afw/geom.h"


#include "testSpatialCell.h"

namespace afwDetect = lsst::afw::detection;
namespace afwImage = lsst::afw::image;
namespace afwMath = lsst::afw::math;
namespace afwGeom = lsst::afw::geom;
typedef float PixelT;

std::pair<afwImage::MaskedImage<PixelT>::Ptr, PTR(afwDetect::FootprintSet)> readImage();

We start by calling readImage, and use boost::tie to unpack the std::pair. The tie call does what you think, unpacking a pair into a couple of variables (it works for boost::tuple too, and is in TR1's <tuple> header).

void SpatialCellSetDemo() {
    afwImage::MaskedImage<PixelT>::Ptr im;
    PTR(afwDetect::FootprintSet) fs;
    boost::tie(im, fs) = readImage();

We want to learn something about the objects in this image, and would like to ensure that the ones we study are spread reasonably uniformly. We accordingly create a SpatialCellSet; a collection of SpatialCells, each of which will maintain its one list of candidate objects for study. For example, if we were estimating the PSF we'd want a set of isolated bright stars, but we wouldn't want them to all come from the top right corner of the image. A SpatialCellSet allows us to take the best n candidates from each SpatialCell, ensuring a reasonable distribution across the image.

The constructor's first argument is the image's bounding box — it'd be nice to simply pass the image, wouldn't it, but that's not currently supported. The second and third arguments 260, 200 define the size (in pixels) of the SpatialCells.

If you run the python version of this example, spatialCellExample.py, with display = True the 6 cells will be shown in green (why 6? Because the image is 512x512 and you can fit 260x200 into 512x512 6 times.)

    /*
     * Create an (empty) SpatialCellSet
     */
    afwMath::SpatialCellSet cellSet(
        im->getBBox(),
        260, 200
    );

Our SpatialCellSet is empty, so let's insert all the objects in the frame into it. We have a list of detections in the FootprintSet fs, so this is easy. We package each object into an ExampleCandidate, and insert it into the set. The SpatialCellSet is responsible for putting it into the correct cell, and SpatialCell for maintaining an order within each cell; this ordering is defined by a virtual function double ExampleCandidate::getCandidateRating() const. The ExampleCandidate class is implemented in testSpatialCell.h and testSpatialCell.cc

You can store anything you like in your candidate class, the only requirement is that it inherit from lsst::afw::math::SpatialCellCandidate or lsst::afw::math::SpatialCellImageCandidate (the latter adds some extra virtual methods). I chose to save a pointer to the parent image, and the object's bounding box.

    /*
     * Populate the cellSet using the detected object in the FootprintSet
     */
    for (afwDetect::FootprintSet::FootprintList::iterator ptr = fs->getFootprints()->begin(),
             end = fs->getFootprints()->end(); ptr != end; ++ptr) {
        afwGeom::Box2I const bbox = (*ptr)->getBBox();
        float const xc = (bbox.getMinX() + bbox.getMaxX())/2.0;
        float const yc = (bbox.getMinY() + bbox.getMaxY())/2.0;
        ExampleCandidate::Ptr tc(new ExampleCandidate(xc, yc, im, bbox));
        cellSet.insertCandidate(tc);
    }

It's possible to iterate over all the objects in a SpatialCellSet (we'll do so in a moment), but the simplest way to visit all cells is to pass in a visitor object. The ExampleCandidateVisitor object (defined in testSpatialCell.h) counts the candidates and the number of pixels contained in their bounding boxes.

    ExampleCandidateVisitor visitor;

    cellSet.visitCandidates(&visitor);
    std::cout << boost::format("There are %d candidates\n") % visitor.getN();

Now we'll visit each of our objects by explicit iteration. The iterator returns a base-class pointer so we need a dynamic_cast (this cast is also available from python via a little swiggery). We decided that we don't like small objects, defined as those with less than 75 pixels in their bounding boxes, so we'll label them as BAD.

    for (unsigned int i = 0; i != cellSet.getCellList().size(); ++i) {
        afwMath::SpatialCell::Ptr cell = cellSet.getCellList()[i];

        for (afwMath::SpatialCell::iterator candidate = cell->begin(), candidateEnd = cell->end();
             candidate != candidateEnd; ++candidate) {
            afwGeom::Box2I box =
                dynamic_cast<ExampleCandidate *>((*candidate).get())->getBBox();
            if (box.getArea() < 75) {
                (*candidate)->setStatus(afwMath::SpatialCellCandidate::BAD);
            }
        }
    }

What does BAD mean (other options are UNKNOWN and GOOD)? Basically that that object is to be ignored. It no longer appears in the size of the SpatialCells, it is skipped by the iterators, and the visitors pass it by. You can turn this behaviour off with setIgnoreBad.

Note that we pass the visitor before we decide to ignore BAD so getN() and getNPix() return the number of good objects/pixels.

    for (unsigned int i = 0; i != cellSet.getCellList().size(); ++i) {
        afwMath::SpatialCell::Ptr cell = cellSet.getCellList()[i];
        cell->visitCandidates(&visitor);

        cell->setIgnoreBad(false);       // include BAD in cell.size()
        std::cout << boost::format("%s nobj=%d N_good=%d NPix_good=%d\n") %
            cell->getLabel() % cell->size() % visitor.getN() % visitor.getNPix();
    }

And count the good candidate again

    cellSet.setIgnoreBad(true);           // don't visit BAD candidates
    cellSet.visitCandidates(&visitor);
    std::cout << boost::format("There are %d good candidates\n") % visitor.getN();
}

Running the example should print

There are 22 candidates
Cell 0x0 nobj=2 N_good=2 NPix_good=1858
Cell 1x0 nobj=2 N_good=1 NPix_good=210
Cell 0x1 nobj=4 N_good=4 NPix_good=1305
Cell 1x1 nobj=4 N_good=1 NPix_good=360
Cell 0x2 nobj=3 N_good=1 NPix_good=99
Cell 1x2 nobj=7 N_good=2 NPix_good=288
There are 11 good candidates

---

Here's the function that reads a FITS file and finds a set of object in it. It isn't really anything to do with SpatialCells, but for completeness...

std::pair<afwImage::MaskedImage<PixelT>::Ptr, PTR(afwDetect::FootprintSet)>
readImage() {

First read a part of the FITS file. We use eups::productDir to find the directory, and only read a part of the image (that's the BBox). The use of a boost::shared_ptr<MaskedImage> (written as MaskedImage::Ptr) is because I want to call the actual constructor in the scope of the try block, but I want to use the image at function scope.

    afwImage::MaskedImage<PixelT>::Ptr mi;

    try {
        std::string dataDir = lsst::utils::getPackageDir("afwdata");

        std::string filename = dataDir + "/CFHT/D4/cal-53535-i-797722_1.fits";
        
        afwGeom::Box2I bbox = afwGeom::Box2I(
            afwGeom::Point2I(270, 2530), 
            afwGeom::Extent2I(512, 512)
        );
        
        lsst::daf::base::PropertySet::Ptr md;
        mi.reset(new afwImage::MaskedImage<PixelT>(filename, md, bbox));
        
    } catch (lsst::pex::exceptions::NotFoundError &e) {
        std::cerr << e << std::endl;
        exit(1);
    }

Subtract the background; the try block is in case the image is too small for a spline fit.

    /*
     * Subtract the background.  We can't fix those pesky cosmic rays, as that's in a dependent product
     * (meas/algorithms)
     */
    afwMath::BackgroundControl bctrl(afwMath::Interpolate::NATURAL_SPLINE);
    bctrl.setNxSample(mi->getWidth()/256 + 1);
    bctrl.setNySample(mi->getHeight()/256 + 1);
    bctrl.getStatisticsControl()->setNumSigmaClip(3.0);
    bctrl.getStatisticsControl()->setNumIter(2);

    afwImage::Image<PixelT>::Ptr im = mi->getImage();
    try {
        *mi->getImage() -= *afwMath::makeBackground(*im, bctrl)->getImage<PixelT>();
    } catch(std::exception &) {
        bctrl.setInterpStyle(afwMath::Interpolate::CONSTANT);
        *mi->getImage() -= *afwMath::makeBackground(*im, bctrl)->getImage<PixelT>();
    }

Run an object detector

    /*
     * Find sources
     */
    afwDetect::Threshold threshold(5, afwDetect::Threshold::STDEV);
    int npixMin = 5;                    // we didn't smooth
    PTR(afwDetect::FootprintSet) fs(new afwDetect::FootprintSet(*mi, threshold, "DETECTED", npixMin));
    int const grow = 1;
    bool const isotropic = false;
    PTR(afwDetect::FootprintSet) grownFs(new afwDetect::FootprintSet(*fs, grow, isotropic));
    grownFs->setMask(mi->getMask(), "DETECTED");

And return the desired data

    return std::make_pair(mi, grownFs);
}