LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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(),
68  dimensions(afw::table::PointKey<int>::addFields(schema, "dimensions",
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
table::Key< std::string > name
Definition: Amplifier.cc:116
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
afw::table::Schema schema
Definition: GaussianPsf.cc:48
int y
Definition: SpanSet.cc:49
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
A circularly symmetric Gaussian Psf class with no spatial variation, intended mostly for testing purp...
Definition: GaussianPsf.h:42
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
Definition: GaussianPsf.cc:117
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: GaussianPsf.cc:121
double getSigma() const
Return the radius of the Gaussian.
Definition: GaussianPsf.h:77
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: GaussianPsf.cc:125
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
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy; should usually be unnecessary because Psfs are immutable.
Definition: GaussianPsf.cc:113
GaussianPsf(int width, int height, double sigma)
Constructor for a GaussianPsf.
Definition: GaussianPsf.cc:99
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
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
Describe the colour of a source.
Definition: Color.h:26
ndarray::Array< Pixel, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:149
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
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
An integer coordinate rectangle.
Definition: Box.h:55
T exp(T... args)
FilterProperty & operator=(FilterProperty const &)=default
A base class for image defects.
int row
Definition: CR.cc:145