LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
GaussianPsf.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * LSST Data Management System
4  * Copyright 2008-2014 LSST Corporation.
5  *
6  * This product includes software developed by the
7  * LSST Project (http://www.lsst.org/).
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the LSST License Statement and
20  * the GNU General Public License along with this program. If not,
21  * see <http://www.lsstcorp.org/LegalNotices/>.
22  */
23 
24 #include <memory>
25 
26 #include "ndarray/eigen.h"
27 
34 
35 namespace lsst {
36 namespace afw {
37 
38 template std::shared_ptr<detection::GaussianPsf> table::io::PersistableFacade<
39  detection::GaussianPsf>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
40 
41 namespace detection {
42 
43 namespace {
44 
45 // Read-only singleton struct containing the schema and keys that a GaussianPsf is mapped
46 // to in record persistence.
47 struct GaussianPsfPersistenceHelper {
48  afw::table::Schema schema;
49  afw::table::PointKey<int> dimensions;
50  afw::table::Key<double> sigma;
51 
52  static GaussianPsfPersistenceHelper const& get() {
53  static GaussianPsfPersistenceHelper instance;
54  return instance;
55  }
56 
57  // No copying
58  GaussianPsfPersistenceHelper(const GaussianPsfPersistenceHelper&) = delete;
59  GaussianPsfPersistenceHelper& operator=(const GaussianPsfPersistenceHelper&) = delete;
60 
61  // No moving
62  GaussianPsfPersistenceHelper(GaussianPsfPersistenceHelper&&) = delete;
63  GaussianPsfPersistenceHelper& operator=(GaussianPsfPersistenceHelper&&) = delete;
64 
65 private:
66  GaussianPsfPersistenceHelper()
67  : schema(),
69  "width/height of realization of Psf", "pixel")),
70  sigma(schema.addField<double>("sigma", "radius of Gaussian", "pixel")) {}
71 };
72 
73 class GaussianPsfFactory : public afw::table::io::PersistableFactory {
74 public:
75  std::shared_ptr<afw::table::io::Persistable> read(InputArchive const& archive,
76  CatalogVector const& catalogs) const override {
77  static GaussianPsfPersistenceHelper const& keys = GaussianPsfPersistenceHelper::get();
78  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
79  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
80  afw::table::BaseRecord const& record = catalogs.front().front();
81  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
82  return std::make_shared<GaussianPsf>(record.get(keys.dimensions.getX()),
83  record.get(keys.dimensions.getY()), record.get(keys.sigma));
84  }
85 
86  GaussianPsfFactory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
87 };
88 
89 GaussianPsfFactory registration("GaussianPsf");
90 
91 void checkDimensions(lsst::geom::Extent2I const& dimensions) {
92  if (dimensions.getX() % 2 == 0 || dimensions.getY() % 2 == 2) {
93  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "GaussianPsf dimensions must be odd");
94  }
95 }
96 
97 } // namespace
98 
99 GaussianPsf::GaussianPsf(int width, int height, double sigma)
100  : Psf(true), _dimensions(width, height), _sigma(sigma) {
101  checkDimensions(_dimensions);
102 }
103 
105  : Psf(true), _dimensions(dimensions), _sigma(sigma) {
106  checkDimensions(_dimensions);
107 }
108 
109 GaussianPsf::GaussianPsf(GaussianPsf const&) = default;
111 GaussianPsf::~GaussianPsf() = default;
112 
114  return std::make_shared<GaussianPsf>(_dimensions, _sigma);
115 }
116 
118  return std::make_shared<GaussianPsf>(lsst::geom::Extent2I(width, height), _sigma);
119 }
120 
121 std::string GaussianPsf::getPersistenceName() const { return "GaussianPsf"; }
122 
123 std::string GaussianPsf::getPythonModule() const { return "lsst.afw.detection"; }
124 
126  static GaussianPsfPersistenceHelper const& keys = GaussianPsfPersistenceHelper::get();
127  afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
129  (*record)[keys.dimensions.getX()] = _dimensions.getX();
130  (*record)[keys.dimensions.getY()] = _dimensions.getY();
131  (*record)[keys.sigma] = getSigma();
132  handle.saveCatalog(catalog);
133 }
134 
135 std::shared_ptr<GaussianPsf::Image> GaussianPsf::doComputeKernelImage(lsst::geom::Point2D const&,
136  image::Color const&) const {
138  Image::Array array = r->getArray();
139  double sum = 0.0;
140  for (int yIndex = 0, y = r->getY0(); yIndex < _dimensions.getY(); ++yIndex, ++y) {
141  Image::Array::Reference row = array[yIndex];
142  for (int xIndex = 0, x = r->getX0(); xIndex < _dimensions.getX(); ++xIndex, ++x) {
143  sum += row[xIndex] = std::exp(-0.5 * (x * x + y * y) / (_sigma * _sigma));
144  }
145  }
146  ndarray::asEigenMatrix(array) /= sum;
147  return r;
148 }
149 
150 double GaussianPsf::doComputeApertureFlux(double radius, lsst::geom::Point2D const& position,
151  image::Color const& color) const {
152  return 1.0 - std::exp(-0.5 * radius * radius / (_sigma * _sigma));
153 }
154 
155 geom::ellipses::Quadrupole GaussianPsf::doComputeShape(lsst::geom::Point2D const& position,
156  image::Color const& color) const {
157  return geom::ellipses::Quadrupole(_sigma * _sigma, _sigma * _sigma, 0.0);
158 }
159 
160 lsst::geom::Box2I GaussianPsf::doComputeBBox(lsst::geom::Point2D const& position,
161  image::Color const& color) const {
162  return lsst::geom::Box2I(lsst::geom::Point2I(-_dimensions / 2),
163  _dimensions); // integer truncation intentional
164 }
165 } // namespace detection
166 } // namespace afw
167 } // namespace lsst
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
An object passed to Persistable::write to allow it to persist itself.
T exp(T... args)
afw::table::Schema schema
Definition: GaussianPsf.cc:48
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
Definition: GaussianPsf.cc:123
GaussianPsf(int width, int height, double sigma)
Constructor for a GaussianPsf.
Definition: GaussianPsf.cc:99
int y
Definition: SpanSet.cc:49
static PointKey addFields(Schema &schema, std::string const &name, std::string const &doc, std::string const &unit)
Add a pair of _x, _y fields to a Schema, and return a PointKey that points to them.
Definition: aggregates.cc:36
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
Definition: GaussianPsf.cc:117
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
STL class.
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
A base class for image defects.
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy; should usually be unnecessary because Psfs are immutable.
Definition: GaussianPsf.cc:113
double x
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Extent< int, 2 > Extent2I
Definition: Extent.h:397
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
double getSigma() const
Return the radius of the Gaussian.
Definition: GaussianPsf.h:77
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
Describe the colour of a source.
Definition: Color.h:26
A polymorphic base class for representing an image&#39;s Point Spread Function.
Definition: Psf.h:75
lsst::geom::Box2I computeBBox(lsst::geom::Point2D position=makeNullPoint(), image::Color color=image::Color()) const
Return the bounding box of the image returned by computeKernelImage()
Definition: Psf.cc:128
ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:149
An integer coordinate rectangle.
Definition: Box.h:55
A circularly symmetric Gaussian Psf class with no spatial variation, intended mostly for testing purp...
Definition: GaussianPsf.h:42
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: GaussianPsf.cc:125
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:485
int row
Definition: CR.cc:145
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: GaussianPsf.cc:121