Psf: use new lsst::utils::Cache

This allows caching many more realisations (they're relatively cheap:
50x50 pixels per realisation @ 8 bytes per pixel); the number of
realisations is configurable.

This gains about 22% in the runtime of measureCoaddSources.py using the
default 100 cache slots for all Psf instances. Allowing CoaddPsf to use
100000 cache slots increases the gain to about 25%.

include/lsst/afw/detection/Psf.h

@@ -29,6 +29,7 @@
 #include <memory>
 
 #include "lsst/daf/base.h"
+#include "lsst/utils/Cache.h"
 #include "lsst/afw/geom/ellipses/Quadrupole.h"
 #include "lsst/afw/math/Kernel.h"
 #include "lsst/afw/image/Color.h"
@@ -254,14 +255,28 @@ class Psf : public daf::base::Citizen,
                                                       std::string const& warpAlgorithm = "lanczos5",
                                                       unsigned int warpBuffer = 5);
 
+    /// Return the capacity of the caches
+    ///
+    /// Both the image and kernel image caches have the same capacity.
+    std::size_t getCacheCapacity() const { return _kernelImageCache.capacity(); }
+
+    /// Set the capacity of the caches
+    ///
+    /// Both the image and kernel image caches will be set to this capacity.
+    void setCacheCapacity(std::size_t capacity) {
+        _imageCache.reserve(capacity);
+        _kernelImageCache.reserve(capacity);
+    }
+
 protected:
     /**
      *  Main constructor for subclasses.
      *
      *  @param[in] isFixed  Should be true for Psf for which doComputeKernelImage always returns
      *                      the same image, regardless of color or position arguments.
+     *  @param[in] capacity  Capacity of the caches.
      */
-    explicit Psf(bool isFixed = false);
+    explicit Psf(bool isFixed = false, std::size_t capacity=100);
 
 private:
     //@{
@@ -283,13 +298,29 @@ class Psf : public daf::base::Citizen,
     virtual geom::Box2I doComputeBBox(geom::Point2D const& position, image::Color const& color) const = 0;
     //@}
 
+    struct CacheKey {
+        geom::Point2D position;
+        image::Color color;
+
+        CacheKey(geom::Point2D const& position_, image::Color color_=image::Color())
+          : position(position_), color(color_) {}
+
+        friend std::size_t hash_value(CacheKey const& key) {
+            return boost::hash<geom::Point2D>()(key.position);  // Currently don't care about color
+        }
+
+        bool operator==(CacheKey const& other) const {
+            return position == other.position;  // Currently don't care about color
+        }
+
+        friend std::ostream & operator<<(std::ostream &os, CacheKey const &key) {
+            return os << key.position;
+        }
+    };
+
     bool const _isFixed;
-    mutable std::shared_ptr<Image> _cachedImage;
-    mutable std::shared_ptr<Image> _cachedKernelImage;
-    mutable image::Color _cachedImageColor;
-    mutable image::Color _cachedKernelImageColor;
-    mutable geom::Point2D _cachedImagePosition;
-    mutable geom::Point2D _cachedKernelImagePosition;
+    mutable utils::Cache<CacheKey, std::shared_ptr<Image>> _imageCache;
+    mutable utils::Cache<CacheKey, std::shared_ptr<Image>> _kernelImageCache;
 
     LSST_PERSIST_FORMATTER(PsfFormatter)
 };

python/lsst/afw/detection/psf.cc

@@ -74,6 +74,8 @@ PYBIND11_PLUGIN(_psf) {
     cls.def("getAveragePosition", &Psf::getAveragePosition);
     cls.def_static("recenterKernelImage", &Psf::recenterKernelImage, "im"_a, "position"_a,
                    "warpAlgorithm"_a = "lanczos5", "warpBuffer"_a = 5);
+    cls.def("getCacheCapacity", &Psf::getCacheCapacity);
+    cls.def("setCacheCapacity", &Psf::setCacheCapacity);
 
     return mod.ptr();
 }

src/detection/Psf.cc

@@ -12,19 +12,16 @@ namespace detection {
 
 namespace {
 
-// Comparison function that determines when we used the cached image instead of recomputing it.
-// We'll probably want a tolerance for colors someday too, but they're just a placeholder right now
-// so it's not worth the effort.
-bool comparePsfEvalPoints(geom::Point2D const &a, geom::Point2D const &b) {
-    // n.b. desired tolerance is actually sqrt(eps), so tolerance squared is eps.
-    return (a - b).computeSquaredNorm() < std::numeric_limits<double>::epsilon();
-}
-
 bool isPointNull(geom::Point2D const &p) { return std::isnan(p.getX()) && std::isnan(p.getY()); }
 
 }  // anonymous
 
-Psf::Psf(bool isFixed) : daf::base::Citizen(typeid(this)), _isFixed(isFixed) {}
+Psf::Psf(bool isFixed, std::size_t capacity)
+  : daf::base::Citizen(typeid(this)),
+  _isFixed(isFixed),
+  _imageCache(capacity),
+  _kernelImageCache(capacity)
+{}
 
 std::shared_ptr<image::Image<double>> Psf::recenterKernelImage(std::shared_ptr<Image> im,
                                                                geom::Point2D const &position,
@@ -46,15 +43,10 @@ std::shared_ptr<Psf::Image> Psf::computeImage(geom::Point2D position, image::Col
                                               ImageOwnerEnum owner) const {
     if (isPointNull(position)) position = getAveragePosition();
     if (color.isIndeterminate()) color = getAverageColor();
-    std::shared_ptr<Psf::Image> result;
-    if (_cachedImage && color == _cachedImageColor && comparePsfEvalPoints(position, _cachedImagePosition)) {
-        result = _cachedImage;
-    } else {
-        result = doComputeImage(position, color);
-        _cachedImage = result;
-        _cachedImageColor = color;
-        _cachedImagePosition = position;
-    }
+    std::shared_ptr<Psf::Image> result = _imageCache(
+        CacheKey(position, color),
+        [this](CacheKey const& key) { return doComputeImage(key.position, key.color); }
+    );
     if (owner == COPY) {
         result = std::make_shared<Image>(*result, true);
     }
@@ -63,18 +55,12 @@ std::shared_ptr<Psf::Image> Psf::computeImage(geom::Point2D position, image::Col
 
 std::shared_ptr<Psf::Image> Psf::computeKernelImage(geom::Point2D position, image::Color color,
                                                     ImageOwnerEnum owner) const {
-    if (isPointNull(position)) position = getAveragePosition();
-    if (color.isIndeterminate()) color = getAverageColor();
-    std::shared_ptr<Psf::Image> result;
-    if (_cachedKernelImage && (_isFixed || (color == _cachedKernelImageColor &&
-                                            comparePsfEvalPoints(position, _cachedKernelImagePosition)))) {
-        result = _cachedKernelImage;
-    } else {
-        result = doComputeKernelImage(position, color);
-        _cachedKernelImage = result;
-        _cachedKernelImageColor = color;
-        _cachedKernelImagePosition = position;
-    }
+    if (_isFixed || isPointNull(position)) position = getAveragePosition();
+    if (_isFixed || color.isIndeterminate()) color = getAverageColor();
+    std::shared_ptr<Psf::Image> result = _kernelImageCache(
+        CacheKey(position, color),
+        [this](CacheKey const& key) { return doComputeKernelImage(key.position, key.color); }
+    );
     if (owner == COPY) {
         result = std::make_shared<Image>(*result, true);
     }