Tickets/dm 2961 - psfex port

bin/psfex.py

@@ -2,12 +2,12 @@
 import argparse
 import os
 import sys
-from lsst.meas.extensions.psfex import *
+from lsst.meas.extensions.psfex import psfex, readPrefs, makeitLsst, makeit, showPsf, dafBase
 
 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 
 if __name__ == "__main__":
-    parser = argparse.ArgumentParser(description="PSFEX")
+    parser = argparse.ArgumentParser(description="Wrapper for Point Spread Function Extractor (PSFEX)")
 
     parser.add_argument('catalogs', type=str, nargs="+", help="Input catalogues from SExtractor")
     parser.add_argument('-c', type=str, dest="defaultsFile",
@@ -17,12 +17,12 @@
                         help="Overrides for default parameters", default=[])
     parser.add_argument('--plot', type=str, nargs="+",
                         help="Desired plots", default=[])
-    parser.add_argument('--ds9', type=int, 
+    parser.add_argument('--ds9', type=int,
                         help="Show the PSF on ds9", default=None)
     parser.add_argument('--diagnostics', action="store_true",
                         help="Write output diagnostic plots to outDir")
     parser.add_argument('--verbose', action="store_true", help="How chatty should I be?", default=False)
-    
+
     argv = sys.argv[:]                  # argparse will mess with sys.argv
     args = parser.parse_args()
 

config/default-lsst.psfex

@@ -10,13 +10,13 @@ PSF_ACCURACY    0.01            # Accuracy to expect from PSF "pixel" values
 
 #------------------------- Point source measurements -------------------------
 
-CENTER_KEYS     initial.centroid.sdss.x,initial.centroid.sdss.y # Catalogue parameters for source pre-centering
-PHOTFLUX_KEY    initial.flux.sinc    # Catalogue parameter for photometric norm.
-PHOTFLUXERR_KEY initial.flux.sinc.err # Catalogue parameter for photometric error
+CENTER_KEYS     slot_Centroid_x,slot_Centroid_y # Catalogue parameters for source pre-centering
+PHOTFLUX_KEY    slot_ApFlux_flux    # Catalogue parameter for photometric norm.
+PHOTFLUXERR_KEY slot_ApFlux_fluxErr # Catalogue parameter for photometric error
 
 #----------------------------- PSF variability -------------------------------
 
-PSFVAR_KEYS     initial.centroid.sdss.x,initial.centroid.sdss.y # Catalogue or FITS (preceded by :) params
+PSFVAR_KEYS     slot_Centroid_x,slot_Centroid_y # Catalogue or FITS (preceded by :) params
 PSFVAR_GROUPS   1,1             # Group tag for each context key
 PSFVAR_DEGREES  2               # Polynom degree for each group
 

include/lsst/meas/extensions/psfex/Field.hh

@@ -1,4 +1,5 @@
 // -*- lsst-C++ -*-
+// Wrapper and extension of the psfex field.h header
 #if !defined(ASTROMATIC_PSFEX_FIELD_HH)
 #define ASTROMATIC_PSFEX_FIELD_HH
 
@@ -21,6 +22,7 @@ struct structtab; typedef struct structtab tabstruct;
 
 #include "define.h"
 #include "field.h"
+#undef VERSION
 }
 
 namespace lsst { namespace meas { namespace extensions { namespace psfex {

include/lsst/meas/extensions/psfex/prefs.hh

@@ -1,4 +1,5 @@
 // -*- lsst-C++ -*-
+// Wrapper and extension of the psfex prefs.h header
 #if !defined(ASTROMATIC_PSFEX_PREFS_HH)
 #define ASTROMATIC_PSFEX_PREFS_HH
 
@@ -15,6 +16,7 @@ struct structtab; typedef struct structtab tabstruct;
 
 #include "define.h"
 #include "prefs.h"
+#undef VERSION
 }
 
 namespace lsst { namespace meas { namespace extensions { namespace psfex {

python/lsst/meas/extensions/psfex/psfex.py

@@ -378,7 +378,6 @@ def load_samples(prefs, context, ext=psfex.Prefs.ALL_EXTENSIONS, next=1, plot=di
     ncat = len(filenames)
     fwhmmin = np.empty(ncat)
     fwhmmax = np.empty(ncat)
-    fwhmode = np.empty(ncat)
 
     if not prefs.getAutoselectFlag():
         fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
@@ -747,7 +746,6 @@ def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, co
     shape = tab.getShapeDefinition()
     ixx = tab.get("%s.xx" % shape)
     iyy = tab.get("%s.yy" % shape)
-    ixy = tab.get("%s.xy" % shape)
 
     rmsSize = np.sqrt(0.5*(ixx + iyy))
     elong = 0.5*(ixx - iyy)/(ixx + iyy)
@@ -756,7 +754,6 @@ def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, co
     fluxErr = tab.get(prefs.getPhotfluxerrRkey())
     flags = getLsstFlags(tab)
 
-    nobj = len(xm)
     #
     # Now the VIGNET data
     #
@@ -783,7 +780,7 @@ def read_samplesLsst(prefs, set, filename, frmin, frmax, ext, next, catindex, co
         try:
             pstamp = mi[x - vigw//2:x + vigw//2 + 1, y - vigh//2:y + vigh//2 + 1]
             vignet[i] = pstamp.getImage().getArray().transpose()
-        except Exception, e:
+        except Exception:
             flags[i] |= edgeBit         # mark star as bad
 
     # Try to load the set of context keys
@@ -871,7 +868,6 @@ def load_samplesLsst(prefs, context, ext=psfex.Prefs.ALL_EXTENSIONS, next=1, plo
     ncat = len(filenames)
     fwhmmin = np.empty(ncat)
     fwhmmax = np.empty(ncat)
-    fwhmode = np.empty(ncat)
 
     if not prefs.getAutoselectFlag():
         fwhmmin = np.zeros(ncat) + prefs.getFwhmrange()[0]
@@ -896,7 +892,6 @@ def load_samplesLsst(prefs, context, ext=psfex.Prefs.ALL_EXTENSIONS, next=1, plo
             shape = tab.getShapeDefinition()
             ixx = tab.get("%s.xx" % shape)
             iyy = tab.get("%s.yy" % shape)
-            ixy = tab.get("%s.xy" % shape)
 
             rmsSize = np.sqrt(0.5*(ixx + iyy))
             elong = 0.5*(ixx - iyy)/(ixx + iyy)
@@ -1020,9 +1015,7 @@ def makeit(prefs, context, saveWcs=False, plot=dict()):
         field.finalize()
         fields.append(field)
 
-    next = fields[0].getNext()          # number of extensions
 
-    psfstep = prefs.getPsfStep()
 
     sets = psfex.vectorSet()
     for set in load_samples(prefs, context, plot=plot):

python/lsst/meas/extensions/psfex/psfexPsfDeterminer.py

@@ -253,18 +253,18 @@ def determinePsf(self, exposure, psfCandidateList, metadata=None, flagKey=None):
         prefs.addCatalog("psfexPsfDeterminer")
 
         prefs.use()
-
+        principalComponentExclusionFlag = bool(bool(psfex.Context.REMOVEHIDDEN)
+                if False else psfex.Context.KEEPHIDDEN)
         context = psfex.Context(prefs.getContextName(), prefs.getContextGroup(),
-                                prefs.getGroupDeg(),
-                                psfex.Context.REMOVEHIDDEN if False else psfex.Context.KEEPHIDDEN)
+                                prefs.getGroupDeg(),principalComponentExclusionFlag)
         set = psfex.Set(context)
         set.setVigSize(pixKernelSize, pixKernelSize)
         set.setFwhm(2*np.sqrt(2*np.log(2))*np.median(sizes))
         set.setRecentroid(self.config.recentroid)
 
         catindex, ext = 0, 0
         backnoise2 = afwMath.makeStatistics(mi.getImage(), afwMath.VARIANCECLIP).getValue()
-        ccd = afwCG.cast_Ccd(exposure.getDetector())
+        ccd = exposure.getDetector()
         if ccd:
             gain = np.mean(np.array([a.getElectronicParams().getGain() for a in ccd]))
         else:

src/PsfexPsf.cc

@@ -43,6 +43,7 @@ extern "C" {
 }
 #include "lsst/meas/extensions/psfex/PsfexPsf.h"
 #include "lsst/meas/algorithms/KernelPsfFactory.h"
+#include "lsst/afw/table/aggregates.h"
 
 namespace lsst { namespace meas { namespace extensions { namespace psfex {
 
@@ -98,7 +99,7 @@ PsfexPsf::getKernel(afw::geom::Point2D position) const
     double pos[MAXCONTEXT];
     int const ndim = _context.size();
     if (ndim != 2) {                    // we're only handling spatial variation for now
-        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterException,
+        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                           str(boost::format("Only spatial variation (ndim == 2) is supported; saw %d")
                               % ndim));
 
@@ -188,7 +189,7 @@ PsfexPsf::_doComputeImage(afw::geom::Point2D const& position,
     double pos[MAXCONTEXT];
     int const ndim = _context.size();
     if (ndim != 2) {                    // we're only handling spatial variation for now
-        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterException,
+        throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
                           str(boost::format("Only spatial variation (ndim == 2) is supported; saw %d")
                               % ndim));
 
@@ -279,8 +280,7 @@ class PsfexPsfSchema1 {
         _comp_size(schema.addField<int>("_comp_size", "Size of _comp array")),
         _context_size(schema.addField<int>("_context_size", "Size of _context array")),
     // Other scalars
-        averagePosition(schema.addField< afw::table::Point<double> >(
-                        "averagePosition", "average position of stars used to make the PSF")),
+        averagePosition(afw::table::PointKey<double>::addFields(schema,"averagePosition","average position of stars used to make the PSF","pixels")),
         _pixstep(schema.addField<float>("_pixstep", "oversampling", "pixel"))
     {
         ;
@@ -297,7 +297,7 @@ class PsfexPsfSchema1 {
     table::Key<int> _comp_size;
     table::Key<int> _context_size;
     // Other scalars
-    table::Key<table::Point<double> > averagePosition;
+    table::PointKey<double> averagePosition;
     table::Key<float> _pixstep;
 };
 

src/prefsImpl.cc

@@ -14,11 +14,18 @@ Prefs::Prefs(std::string const& filename, lsst::daf::base::PropertySet const* va
         std::vector<char *> argkey(narg);
         std::vector<char *> argval(narg);
         std::vector<std::string> const names = values->paramNames();
+        // Next, decalre a vector to store string values. These strings must
+        // live long enough such that pointers to them (returned by c_str) stay
+        // valid until the readprefs function can read their values. Once this
+        // happens, the srings and pointers are no longer needed, and will be
+        // cleaned up when the vector goes out of scope
+        std::vector<std::string> paramvalues(narg);
 
         for (int i = 0; i != narg; ++i) {
             std::string const& name = names[i];
             argkey[i] = const_cast<char *>(name.c_str());
-            argval[i] = const_cast<char *>(values->getAsString(name).c_str());
+            paramvalues[i] = values->getAsString(name);
+            argval[i] = const_cast<char *>(paramvalues[i].c_str());
         }
 
         readprefs(cfilename, &argkey[0], &argval[0], narg);
@@ -56,7 +63,7 @@ Prefs::setCommandLine(std::vector<std::string> const& argv)
 void
 Prefs::addCatalog(std::string const& filename) {
     if (prefs.ncat >= MAXFILE) {
-        throw LSST_EXCEPT(lsst::pex::exceptions::LengthErrorException, "Too many input catalogues");
+        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError, "Too many input catalogues");
     }
     _catalogs.push_back(filename);
     prefs.incat_name[prefs.ncat++] = const_cast<char *>((_catalogs.end() - 1)->c_str());

src/psfImpl.cc

@@ -66,7 +66,7 @@ void Sample::setVig(ndarray::Array<float,2,2> const& img)
     int const w = img.getSize<0>(), h = img.getSize<1>();
 
     if (_vigsize[0]*_vigsize[1] != w*h) {
-        throw LSST_EXCEPT(lsst::pex::exceptions::LengthErrorException,
+        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                           str(boost::format("Expected %dx%d array, but was passed %dx%d")
                               % _vigsize[0] % _vigsize[1] % w % h));
     }
@@ -146,7 +146,7 @@ Psf::build(double x, double y,
     }
 
     psf_build(impl, &pos[0]);
-}		
+}	
 
 RETURN_IMAGE_FIELD(Psf::getLoc,  loc,  impl->size)
 RETURN_IMAGE_FIELD(Psf::getResi, resi, impl->size)

src/psfexAdaptors.cc

@@ -74,7 +74,7 @@ makeit(std::vector<boost::shared_ptr<Field> > &fields_,
       )
 {
     if (sets.size() > MAXFILE) {
-        throw LSST_EXCEPT(lsst::pex::exceptions::LengthErrorException,
+        throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
                           (boost::format("Too many sets: %d") % sets.size()).str());
     }
 

tests/testPsfexPsf.py

@@ -55,7 +55,8 @@
 import lsst.meas.algorithms.utils as maUtils
 import lsst.afw.cameraGeom as cameraGeom
 import lsst.meas.extensions.psfex.psfexPsfDeterminer as psfexPsfDeterminer
-
+from lsst.meas.base import SingleFrameMeasurementTask
+from lsst.meas.algorithms.detection import SourceDetectionTask
 try:
     type(verbose)
 except NameError:
@@ -84,32 +85,21 @@ class SpatialModelPsfTestCase(unittest.TestCase):
     @staticmethod
     def measure(footprintSet, exposure):
         """Measure a set of Footprints, returning a SourceCatalog"""
-        config = measAlg.SourceMeasurementConfig()
-        config.prefix = "initial."
-        config.algorithms.names = ["flags.pixel", "flux.psf", "flux.sinc", "flux.gaussian", "shape.sdss"]
-        config.centroider.name = "centroid.sdss"
-        config.algorithms["flux.naive"].radius = 3.0
-        config.slots.centroid = "initial.centroid.sdss"
-        config.slots.psfFlux = "initial.flux.psf"
-        config.slots.apFlux = "initial.flux.sinc"
-        config.slots.modelFlux = None
-        config.slots.instFlux = None
-        config.slots.calibFlux = None
-        config.slots.shape = "initial.shape.sdss"
-
+        config = SingleFrameMeasurementTask.ConfigClass()
+        config.slots.apFlux = 'base_CircularApertureFlux_12_0'
         schema = afwTable.SourceTable.makeMinimalSchema()
-        measureSources = config.makeMeasureSources(schema)
+
+        measureSources = SingleFrameMeasurementTask(schema,config=config)
+
+        tab = afwTable.SourceTable.make(schema)
         catalog = afwTable.SourceCatalog(schema)
-        config.slots.setupTable(catalog.table)
 
         if display:
             ds9.mtv(exposure, title="Original", frame=0)
 
         footprintSet.makeSources(catalog)
 
-        for i, source in enumerate(catalog):
-            measureSources.applyWithPeak(source, exposure)
-
+        measureSources.run(catalog,exposure)
         return catalog
 
     def setUp(self):
@@ -134,14 +124,6 @@ def setUp(self):
         wcs = afwImage.makeWcs(crval, afwGeom.PointD(0, 0), 1.0, 0, 0, 1.0)
         self.exposure.setWcs(wcs)
 
-        ccd = cameraGeom.Ccd(cameraGeom.Id(1))
-        ccd.addAmp(cameraGeom.Amp(cameraGeom.Id(0),
-                                  afwGeom.BoxI(afwGeom.PointI(0,0), self.exposure.getDimensions()),
-                                  afwGeom.BoxI(afwGeom.PointI(0,0), afwGeom.ExtentI(0,0)),
-                                  afwGeom.BoxI(afwGeom.PointI(0,0), self.exposure.getDimensions()),
-                                  cameraGeom.ElectronicParams(1.0, 100.0, 65535)))
-        self.exposure.setDetector(ccd)
-        self.exposure.getDetector().setDistortion(None)        
         #
         # Make a kernel with the exactly correct basis functions.  Useful for debugging
         #
@@ -241,13 +223,12 @@ def setupDeterminer(exposure):
         starSelectorFactory = measAlg.starSelectorRegistry["objectSize"]
         starSelectorConfig = starSelectorFactory.ConfigClass()
 
-        starSelectorConfig.sourceFluxField = "initial.flux.gaussian"
-        starSelectorConfig.badFlags = ["initial.flags.pixel.edge",
-                                       "initial.flags.pixel.interpolated.center",
-                                       "initial.flags.pixel.saturated.center",
-                                       "initial.flags.pixel.cr.center",
+        starSelectorConfig.sourceFluxField = "base_GaussianFlux_flux"
+        starSelectorConfig.badFlags = ["base_PixelFlags_flag_edge",
+                                       "base_PixelFlags_flag_interpolatedCenter",
+                                       "base_PixelFlags_flag_saturatedCenter",
+                                       "base_PixelFlags_flag_crCenter",
                                        ]
-        starSelectorConfig.widthStdAllowed = 0.5
 
         starSelector = starSelectorFactory(starSelectorConfig)
 
@@ -270,7 +251,6 @@ def subtractStars(self, exposure, catalog, chi_lim=-1):
         mi, psf = exposure.getMaskedImage(), exposure.getPsf()
 
         subtracted =  mi.Factory(mi, True)
-
         for s in catalog:
             xc, yc = s.getX(), s.getY()
             bbox = subtracted.getBBox(afwImage.PARENT)
@@ -279,7 +259,6 @@ def subtractStars(self, exposure, catalog, chi_lim=-1):
                     measAlg.subtractPsf(psf, subtracted, xc, yc)
                 except:
                     pass
-
         chi = subtracted.Factory(subtracted, True)
         var = subtracted.getVariance()
         np.sqrt(var.getArray(), var.getArray()) # inplace sqrt

ups/meas_extensions_psfex.cfg

@@ -4,8 +4,8 @@ import lsst.sconsUtils
 
 dependencies = {
     "required": [],
-    "buildRequired": ["afw", "boost_test", "daf_base", "meas_algorithms",
-                      "swig", "utils", "psfex"],
+    "buildRequired": ["afw", "boost_test", "daf_base", "fftw",
+                      "meas_base","meas_algorithms","swig", "utils", "psfex"],
 }
 
 config = lsst.sconsUtils.Configuration(

ups/meas_extensions_psfex.table

@@ -1,6 +1,12 @@
+setupRequired(afw)
+setupRequired(boost)
+setupRequired(daf_base)
+setupRequired(fftw)
+setupRequired(swig)
+setupRequired(utils)
+setupRequired(meas_base)
 setupRequired(meas_algorithms)
 setupRequired(psfex)
-#setupRequired(pyfits)
 
 envPrepend(PATH, ${PRODUCT_DIR}/bin)