LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
LSST Data Management Base Package
VeresModel.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2/*
3 * This file is part of meas_extensions_trailedSources.
4 *
5 * Developed for the LSST Data Management System.
6 * This product includes software developed by the LSST Project
7 * (https://www.lsst.org).
8 * See the COPYRIGHT file at the top-level directory of this distribution
9 * for details of code ownership.
10 *
11 * This program is free software: you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * This program is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with this program. If not, see <https://www.gnu.org/licenses/>.
23 */
24
25#include "lsst/geom.h"
26#include "lsst/afw/image.h"
27#include "lsst/afw/detection.h"
29
30namespace lsst {
31namespace meas {
32namespace extensions {
33namespace trailedSources {
34
38
40 ExposureF const& data
41) : _sigma(data.getPsf()->computeShape().getTraceRadius()),
42 _bbox(data.getBBox()),
43 _data(data.getMaskedImage().getImage()->getArray()),
44 _variance(data.getMaskedImage().getVariance()->getArray()) {}
45
46double VeresModel::operator()(std::vector<double> const& params) const {
47
48 double xc = params[0]; // Centroid x
49 double yc = params[1]; // Centroid y
50 double flux = params[2]; // Flux
51 double length = params[3]; // Trail length
52 double theta = params[4]; // Angle from +x-axis
53
54 // Compute model image and chi-squared
55 double chiSq = 0.0;
56 // Loop is adapted from lsst::afw::detection::Psf::computeKernelImage()
57 for (int yIndex = 0, yp = _bbox.getBeginY(); yIndex < _bbox.getHeight(); ++yIndex, ++yp) {
58 ImageF::Array::Reference dataRow = _data[yIndex];
59 ImageF::Array::Reference varRow = _variance[yIndex];
60 for (int xIndex = 0, xp = _bbox.getBeginX(); xIndex < _bbox.getWidth(); ++xIndex, ++xp) {
61 double model = _computeModel(xp,yp,xc,yc,flux,length,theta);
62 double diff = dataRow[xIndex] - model;
63 chiSq += diff*diff/varRow[xIndex];
64 }
65 }
66
67 return chiSq;
68}
69
71
72 double xc = params[0]; // Centroid x
73 double yc = params[1]; // Centroid y
74 double flux = params[2]; // Flux
75 double length = params[3]; // Trail length
76 double theta = params[4]; // Angle from +x-axis
77
78 // Compute gradients of the model and of chi-squared
79 std::vector<double> gradChiSq = {0.0,0.0,0.0,0.0,0.0};
80 for (int yIndex = 0, yp = _bbox.getBeginY(); yIndex < _bbox.getHeight(); ++yIndex, ++yp) {
81 ImageF::Array::Reference dataRow = _data[yIndex];
82 ImageF::Array::Reference varRow = _variance[yIndex];
83 for (int xIndex = 0, xp = _bbox.getBeginX(); xIndex < _bbox.getWidth(); ++xIndex, ++xp) {
84 double model = _computeModel(xp,yp,xc,yc,flux,length,theta);
85 double gradDiff = -2.0 * (dataRow[xIndex] - model) / varRow[xIndex];
86 std::array<double, 5> gradModel = _computeGradient(xp,yp,xc,yc,flux,length,theta);
87 for (int k=0; k<5; ++k) {
88 gradChiSq[k] += gradModel[k] * gradDiff;
89 }
90 }
91 }
92 return gradChiSq;
93}
94
96
97 double xc = params[0]; // Centroid x
98 double yc = params[1]; // Centroid y
99 // double flux = params[2]; // Flux (in this case, always 1)
100 double length = params[3]; // Trail length
101 double theta = params[4]; // Angle from +x-axis
102
103 // Compute the flux and gradient wrt the other model parameters
104 double m2 = 0.0; // sum_i model_i*model_i
105 double md = 0.0; // sum_i model_i*data_i
106 std::vector<double> gradmd = {0.0, 0.0, 0.0, 0.0, 0.0}; // sum_i (gradModel_(i,k)*data_i)
107 std::vector<double> gradmm = {0.0, 0.0, 0.0, 0.0, 0.0}; // sum_i (gradModel_(i,k)*model_i)
108 for (int yIndex = 0, yp = _bbox.getBeginY(); yIndex < _bbox.getHeight(); ++yIndex, ++yp) {
109 ImageF::Array::Reference dataRow = _data[yIndex];
110 for (int xIndex = 0, xp = _bbox.getBeginX(); xIndex < _bbox.getWidth(); ++xIndex, ++xp) {
111 double data = dataRow[xIndex];
112 double model = _computeModel(xp, yp, xc, yc, 1.0, length, theta);
113 std::array<double, 5> gradModel = _computeGradient(xp, yp, xc, yc, 1.0, length, theta);
114 m2 += model*model;
115 md += model*dataRow[xIndex];
116 for (int k=0; k<5; ++k) {
117 gradmd[k] += gradModel[k] * data;
118 gradmm[k] += gradModel[k] * model;
119 }
120 }
121 }
122 double flux = md / m2;
123 std::vector<double> gradFlux = {0.0, 0.0, 0.0, 0.0, 0.0};
124 for (int k=0; k<5; ++k) {
125 gradFlux[k] = (gradmd[k] - 2.0*flux*gradmm[k]) / m2;
126 }
127 gradFlux[2] = 0.0; // Make dfluxdflux = 0
129 return results;
130}
131
133 double xc = params[0]; // Centroid x
134 double yc = params[1]; // Centroid y
135 double flux = params[2]; // Flux
136 double length = params[3]; // Trail length
137 double theta = params[4]; // Angle from +x-axis
138
139 // Loop is adapted from lsst::afw::detection::GaussianPsf::doComputeKernelImage()
141 ImageF::Array array = image->getArray();
142 for (int yIndex = 0, yp = _bbox.getBeginY(); yIndex < _bbox.getHeight(); ++yIndex, ++yp) {
143 ImageF::Array::Reference row = array[yIndex];
144 for (int xIndex = 0, xp = _bbox.getBeginX(); xIndex < _bbox.getWidth(); ++xIndex, ++xp) {
145 row[xIndex] = _computeModel(xp,yp,xc,yc,flux,length,theta);
146 }
147 }
148 return image;
149}
150
151double VeresModel::_computeModel(double x, double y, double xc, double yc,
152 double flux, double length, double theta) const noexcept {
153 double xp = (x-xc)*cos(theta) + (y-yc)*sin(theta);
154 double yp = (x-xc)*sin(theta) - (y-yc)*cos(theta);
155 double A = exp(-0.5 * yp*yp / (_sigma*_sigma));
156 double B = erf((xp+length/2) / (sqrt(2.0) * _sigma));
157 double C = erf((xp-length/2) / (sqrt(2.0) * _sigma));
158 return flux * A * (B - C) / (length * 2 * sqrt(2.0 * geom::PI) * _sigma);
159}
160
161std::array<double, 5> VeresModel::_computeGradient(double x, double y, double xc, double yc,
162 double flux, double length, double theta) const noexcept {
163 double xp = (x-xc)*cos(theta) + (y-yc)*sin(theta);
164 double yp = (x-xc)*sin(theta) - (y-yc)*cos(theta);
165
166 // Duplicated quantities
167 double flux2L = flux/(2.0*length);
168 double ypSq = yp*yp;
169 double sqrt2 = sqrt(2.0);
170 double sqrt2Pi = sqrt(2.0*geom::PI);
171 double sigmaSq = _sigma*_sigma;
172 double sigmaSq8 = sigmaSq * 8.0;
173 double eypSq = exp(-ypSq/(2.0*sigmaSq));
174 double lengthPlus = length+2.0*xp;
175 double lengthMinus= length-2.0*xp;
176 double erfPlus = erf(lengthPlus/(2.0*sqrt2*_sigma));
177 double erfMinus = erf(lengthMinus/(2.0*sqrt2*_sigma));
178 double expPlus = exp(-lengthPlus*lengthPlus/sigmaSq8);
179
180 // Compute partials wrt the transformed coordinates
181 double dfdxp = flux2L/(geom::PI*sigmaSq)*exp(-4.0*ypSq/sigmaSq8)*expPlus*
182 (1.0 - exp(length*xp/sigmaSq));
183 double dfdyp = -flux2L*yp/(sqrt2Pi*_sigma*sigmaSq)*eypSq*(erfMinus+erfPlus);
184
185 // Use the chain rule to get partials wrt the centroid and rotation angle
186 double dxpdxc = -cos(theta);
187 double dxpdyc = -sin(theta);
188 double dxpdTheta = -yp;
189 double dypdxc = -sin(theta);
190 double dypdyc = cos(theta);
191 double dypdTheta = xp;
192 double dfdxc = dfdxp*dxpdxc + dfdyp*dypdxc;
193 double dfdyc = dfdxp*dxpdyc + dfdyp*dypdyc;
194 double dfdTheta = dfdxp*dxpdTheta + dfdyp*dypdTheta;
195
196 double dfdFlux = _computeModel(x,y,xc,yc,1.0,length,theta); // dfdFlux = f / flux
197
198 double dfdLength = flux2L/(length*sqrt2Pi*_sigma)*eypSq*(length/(sqrt2Pi*_sigma)*
199 (exp(-lengthMinus*lengthMinus/sigmaSq8)+expPlus) - erfMinus - erfPlus);
200
201 std::array<double, 5> gradModel = {dfdxc, dfdyc, dfdFlux, dfdLength, dfdTheta};
202 return gradModel;
203}
204
205}}}} // lsst::meas::extensions::trailedSources
char * data
Definition: BaseRecord.cc:61
double x
afw::table::Key< afw::table::Array< ImagePixelT > > image
int y
Definition: SpanSet.cc:48
A class to contain the data, WCS, and other information needed to describe an image of the sky.
Definition: Exposure.h:72
typename ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition: ImageBase.h:149
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:51
int getBeginX() const noexcept
Definition: Box.h:172
int getHeight() const noexcept
Definition: Box.h:188
int getWidth() const noexcept
Definition: Box.h:187
int getBeginY() const noexcept
Definition: Box.h:173
VeresModel(ExposureF const &data)
Constructor for VeresModel.
Definition: VeresModel.cc:39
std::vector< double > gradient(std::vector< double > const &params) const
Compute the gradient of chi-squared of the model given the data.
Definition: VeresModel.cc:70
std::shared_ptr< ImageF > computeModelImage(std::vector< double > const &params) const
Compute an image for a trail generated from the Veres model.
Definition: VeresModel.cc:132
std::tuple< double, std::vector< double > > computeFluxWithGradient(std::vector< double > const &params) const
Computes the flux and the gradient with respect to the other model parameters.
Definition: VeresModel.cc:95
double operator()(std::vector< double > const &params) const
Compute chi-squared of the model given the data.
Definition: VeresModel.cc:46
T cos(T... args)
T erf(T... args)
T exp(T... args)
T make_tuple(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
constexpr double PI
The ratio of a circle's circumference to diameter.
Definition: Angle.h:40
double sin(Angle const &a)
Definition: Angle.h:102
double cos(Angle const &a)
Definition: Angle.h:103
A base class for image defects.
T sin(T... args)
T sqrt(T... args)
int row
Definition: CR.cc:145