26 #include "ndarray/eigen.h"
47 struct GaussianPsfPersistenceHelper {
52 static GaussianPsfPersistenceHelper
const& get() {
53 static GaussianPsfPersistenceHelper instance;
58 GaussianPsfPersistenceHelper(
const GaussianPsfPersistenceHelper&) =
delete;
59 GaussianPsfPersistenceHelper& operator=(
const GaussianPsfPersistenceHelper&) =
delete;
62 GaussianPsfPersistenceHelper(GaussianPsfPersistenceHelper&&) =
delete;
63 GaussianPsfPersistenceHelper& operator=(GaussianPsfPersistenceHelper&&) =
delete;
66 GaussianPsfPersistenceHelper()
69 "width/height of realization of Psf",
"pixel")),
70 sigma(
schema.addField<double>(
"sigma",
"radius of Gaussian",
"pixel")) {}
73 class GaussianPsfFactory :
public afw::table::io::PersistableFactory {
76 CatalogVector
const& catalogs)
const override {
77 static GaussianPsfPersistenceHelper
const&
keys = GaussianPsfPersistenceHelper::get();
80 afw::table::BaseRecord
const& record = catalogs.front().front();
82 return std::make_shared<GaussianPsf>(record.get(
keys.dimensions.getX()),
83 record.get(
keys.dimensions.getY()), record.get(
keys.sigma));
89 GaussianPsfFactory registration(
"GaussianPsf");
93 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
"GaussianPsf dimensions must be odd");
100 :
Psf(true), _dimensions(width, height), _sigma(
sigma) {
101 checkDimensions(_dimensions);
106 checkDimensions(_dimensions);
114 return std::make_shared<GaussianPsf>(_dimensions, _sigma);
126 static GaussianPsfPersistenceHelper
const&
keys = GaussianPsfPersistenceHelper::get();
129 (*record)[
keys.dimensions.getX()] = _dimensions.getX();
130 (*record)[
keys.dimensions.getY()] = _dimensions.getY();
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));
146 ndarray::asEigenMatrix(array) /= sum;
155 geom::ellipses::Quadrupole GaussianPsf::doComputeShape(
lsst::geom::Point2D const& position,
157 return geom::ellipses::Quadrupole(_sigma * _sigma, _sigma * _sigma, 0.0);