LSSTApplications  18.1.0
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 
33 #include "lsst/afw/table/io/Persistable.cc"
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  schema.getCitizen().markPersistent();
72  }
73 };
74 
75 class GaussianPsfFactory : public afw::table::io::PersistableFactory {
76 public:
77  std::shared_ptr<afw::table::io::Persistable> read(InputArchive const& archive,
78  CatalogVector const& catalogs) const override {
79  static GaussianPsfPersistenceHelper const& keys = GaussianPsfPersistenceHelper::get();
80  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
81  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
82  afw::table::BaseRecord const& record = catalogs.front().front();
83  LSST_ARCHIVE_ASSERT(record.getSchema() == keys.schema);
84  return std::make_shared<GaussianPsf>(record.get(keys.dimensions.getX()),
85  record.get(keys.dimensions.getY()), record.get(keys.sigma));
86  }
87 
88  GaussianPsfFactory(std::string const& name) : afw::table::io::PersistableFactory(name) {}
89 };
90 
91 GaussianPsfFactory registration("GaussianPsf");
92 
93 void checkDimensions(lsst::geom::Extent2I const& dimensions) {
94  if (dimensions.getX() % 2 == 0 || dimensions.getY() % 2 == 2) {
95  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "GaussianPsf dimensions must be odd");
96  }
97 }
98 
99 } // namespace
100 
101 GaussianPsf::GaussianPsf(int width, int height, double sigma)
102  : Psf(true), _dimensions(width, height), _sigma(sigma) {
103  checkDimensions(_dimensions);
104 }
105 
107  : Psf(true), _dimensions(dimensions), _sigma(sigma) {
108  checkDimensions(_dimensions);
109 }
110 
111 GaussianPsf::GaussianPsf(GaussianPsf const&) = default;
113 GaussianPsf::~GaussianPsf() = default;
114 
116  return std::make_shared<GaussianPsf>(_dimensions, _sigma);
117 }
118 
120  return std::make_shared<GaussianPsf>(lsst::geom::Extent2I(width, height), _sigma);
121 }
122 
123 std::string GaussianPsf::getPersistenceName() const { return "GaussianPsf"; }
124 
125 std::string GaussianPsf::getPythonModule() const { return "lsst.afw.detection"; }
126 
128  static GaussianPsfPersistenceHelper const& keys = GaussianPsfPersistenceHelper::get();
129  afw::table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
131  (*record)[keys.dimensions.getX()] = _dimensions.getX();
132  (*record)[keys.dimensions.getY()] = _dimensions.getY();
133  (*record)[keys.sigma] = getSigma();
134  handle.saveCatalog(catalog);
135 }
136 
137 std::shared_ptr<GaussianPsf::Image> GaussianPsf::doComputeKernelImage(lsst::geom::Point2D const&,
138  image::Color const&) const {
140  Image::Array array = r->getArray();
141  double sum = 0.0;
142  for (int yIndex = 0, y = r->getY0(); yIndex < _dimensions.getY(); ++yIndex, ++y) {
143  Image::Array::Reference row = array[yIndex];
144  for (int xIndex = 0, x = r->getX0(); xIndex < _dimensions.getX(); ++xIndex, ++x) {
145  sum += row[xIndex] = std::exp(-0.5 * (x * x + y * y) / (_sigma * _sigma));
146  }
147  }
148  ndarray::asEigenMatrix(array) /= sum;
149  return r;
150 }
151 
152 double GaussianPsf::doComputeApertureFlux(double radius, lsst::geom::Point2D const& position,
153  image::Color const& color) const {
154  return 1.0 - std::exp(-0.5 * radius * radius / (_sigma * _sigma));
155 }
156 
157 geom::ellipses::Quadrupole GaussianPsf::doComputeShape(lsst::geom::Point2D const& position,
158  image::Color const& color) const {
159  return geom::ellipses::Quadrupole(_sigma * _sigma, _sigma * _sigma, 0.0);
160 }
161 
162 lsst::geom::Box2I GaussianPsf::doComputeBBox(lsst::geom::Point2D const& position,
163  image::Color const& color) const {
164  return lsst::geom::Box2I(lsst::geom::Point2I(-_dimensions / 2),
165  _dimensions); // integer truncation intentional
166 }
167 } // namespace detection
168 } // namespace afw
169 } // namespace lsst
An ellipse core with quadrupole moments as parameters.
Definition: Quadrupole.h:47
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
An object passed to Persistable::write to allow it to persist itself.
T exp(T... args)
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:125
GaussianPsf(int width, int height, double sigma)
Constructor for a GaussianPsf.
Definition: GaussianPsf.cc:101
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:119
STL class.
A base class for image defects.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy; should usually be unnecessary because Psfs are immutable.
Definition: GaussianPsf.cc:115
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
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
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:76
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:129
ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:150
An integer coordinate rectangle.
Definition: Box.h:54
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:85
afw::table::Schema schema
Definition: GaussianPsf.cc:48
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: GaussianPsf.cc:127
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:123