.. py:currentmodule:: lsst.afw.math

.. _lsst.afw.math-SpatialCellSetExample:

#######################################
Example of lsst.afw.math.SpatialCellSet
#######################################

Demonstrate the use of `SpatialCellSet`\ s; the code's in `spatialCellExample.cc`_.

.. _spatialCellExample.cc: https://github.com/lsst/afw/blob/master/examples/spatialCellExample.cc
.. _spatialCellExample.py: https://github.com/lsst/afw/blob/master/examples/spatialCellExample.py
.. _testSpatialCell.h: https://github.com/lsst/afw/blob/master/examples/testSpatialCell.h
.. _testSpatialCell.cc: https://github.com/lsst/afw/blob/master/examples/testSpatialCell.cc

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

.. code-block:: cpp

   #include <string>
   #include "lsst/geom.h"
   #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 "testSpatialCell.h"

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

   std::pair<std::shared_ptr<afwImage::MaskedImage<PixelT>>, std::shared_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).

.. code-block:: cpp

   void SpatialCellSetDemo() {
       std::shared_ptr<afwImage::MaskedImage<PixelT>> im;
       std::shared_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 `SpatialCell`\ s, 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 `SpatialCell`\ s.

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.)

.. code-block:: cpp
   :name: dummy-to-force-indent1

       /*
        * 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 `SpatialCellCandidate` or `SpatialCellImageCandidate` (the latter adds some extra virtual methods).
I chose to save a pointer to the parent image, and the object's bounding box.

.. code-block:: cpp
   :name: dummy-to-force-indent2

       /*
        * 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) {
           lsst::geom::Box2I const bbox = (*ptr)->getBBox();
           float const xc = (bbox.getMinX() + bbox.getMaxX()) / 2.0;
           float const yc = (bbox.getMinY() + bbox.getMaxY()) / 2.0;
           std::shared_ptr<ExampleCandidate> 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.

.. code-block:: cpp
   :name: dummy-to-force-indent3

       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 `~SpatialCellCandidate.BAD`.

.. code-block:: cpp
   :name: dummy-to-force-indent4

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

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


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

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

.. code-block:: cpp
   :name: dummy-to-force-indent5

       for (unsigned int i = 0; i != cellSet.getCellList().size(); ++i) {
           std::shared_ptr<afwMath::SpatialCell> 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

.. code-block:: cpp
   :name: dummy-to-force-indent6

       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

.. code-block:: none

   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 `SpatialCell`\ s, but for completeness...

.. code-block:: cpp

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

First read a part of the FITS file.
We use `lsst.utils.getPackageDir` 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.

.. code-block:: cpp
   :name: dummy-to-force-indent7

       std::shared_ptr<afwImage::MaskedImage<PixelT>> mi;

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

           std::string filename = dataDir + "/CFHT/D4/cal-53535-i-797722_1.fits";

           lsst::geom::Box2I bbox =
                   lsst::geom::Box2I(lsst::geom::Point2I(270, 2530), lsst::geom::Extent2I(512, 512));

           std::shared_ptr<lsst::daf::base::PropertySet> 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.

.. code-block:: cpp
   :name: dummy-to-force-indent8

       /*
        * 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);

       std::shared_ptr<afwImage::Image<PixelT>> 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

.. code-block:: cpp
   :name: dummy-to-force-indent9

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

And return the desired data

.. code-block:: cpp

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