LSSTApplications  20.0.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 
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
y
int y
Definition: SpanSet.cc:49
lsst::afw::image::ImageBase::Array
ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:149
std::string
STL class.
std::shared_ptr
STL class.
lsst::afw::detection::GaussianPsf::~GaussianPsf
~GaussianPsf() override
lsst::afw::table::io::OutputArchiveHandle
An object passed to Persistable::write to allow it to persist itself.
Definition: OutputArchive.h:118
Persistable.cc
lsst::afw::table::io::OutputArchiveHandle::saveCatalog
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
Definition: OutputArchive.cc:211
lsst::afw
Definition: imageAlgorithm.dox:1
lsst::afw::detection::GaussianPsf::clone
std::shared_ptr< afw::detection::Psf > clone() const override
Polymorphic deep copy; should usually be unnecessary because Psfs are immutable.
Definition: GaussianPsf.cc:113
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:50
lsst::afw::detection::GaussianPsf::getPythonModule
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
astshim.keyMap.keyMapContinued.keys
def keys(self)
Definition: keyMapContinued.py:6
CatalogVector.h
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
aggregates.h
schema
afw::table::Schema schema
Definition: GaussianPsf.cc:48
lsst::afw::detection::GaussianPsf::GaussianPsf
GaussianPsf(int width, int height, double sigma)
Constructor for a GaussianPsf.
Definition: GaussianPsf.cc:99
lsst::afw::table::io::OutputArchiveHandle::makeCatalog
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
Definition: OutputArchive.cc:207
LSST_ARCHIVE_ASSERT
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
lsst::afw::detection::GaussianPsf::getSigma
double getSigma() const
Return the radius of the Gaussian.
Definition: GaussianPsf.h:77
lsst::afw::detection::Psf::Image
image::Image< Pixel > Image
Image type returned by computeImage.
Definition: Psf.h:83
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::afw::detection::GaussianPsf
A circularly symmetric Gaussian Psf class with no spatial variation, intended mostly for testing purp...
Definition: GaussianPsf.h:42
GaussianPsf.h
dimensions
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
lsst::afw::detection::GaussianPsf::resized
std::shared_ptr< afw::detection::Psf > resized(int width, int height) const override
Return a clone with specified kernel dimensions.
Definition: GaussianPsf.cc:117
lsst::afw::detection::GaussianPsf::write
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: GaussianPsf.cc:125
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
lsst.pipe.tasks.insertFakes.radius
radius
Definition: insertFakes.py:288
row
int row
Definition: CR.cc:145
std::exp
T exp(T... args)
lsst::geom::Point
A coordinate class intended to represent absolute positions.
Definition: CoordinateBase.h:39
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
InputArchive.h
lsst::afw::detection::GaussianPsf::getPersistenceName
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: GaussianPsf.cc:121
lsst::afw::detection::Psf
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
lsst::afw::table::CatalogT::addNew
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
lsst::afw::detection::Psf::computeBBox
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
lsst::afw::table::CatalogT< BaseRecord >
lsst::geom::Extent< int, 2 >
OutputArchive.h
lsst::afw::image::Color
Describe the colour of a source.
Definition: Color.h:26