LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
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(data.getPsf()->getAveragePosition()).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
afw::table::Key< afw::table::Array< ImagePixelT > > image
int y
Definition SpanSet.cc:48
typename ndarray::Array< PixelT, 2, 1 > Array
A mutable ndarray representation of the image.
Definition ImageBase.h:149
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.
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)
double constexpr PI
The ratio of a circle's circumference to diameter.
Definition Angle.h:40
T sin(T... args)
T sqrt(T... args)
int row
Definition CR.cc:145