LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
_calcCompensatedGaussian.cc
Go to the documentation of this file.
1#include "pybind11/pybind11.h"
2#include "pybind11/numpy.h"
4#include <math.h>
5#include <vector>
6#include <stdio.h>
7
8namespace py = pybind11;
9
10namespace lsst {
11namespace meas {
12namespace base {
13
14py::tuple compensated_gaussian_filt_inner_product(py::array_t<double, py::array::c_style | py::array::forcecast> & array,
15 py::array_t<double, py::array::c_style | py::array::forcecast> & variance_array,
16 double x_mean, double y_mean, double sig, double t) {
17 //Verify the input array is conditioned in an appropriate manner.
18 py::buffer_info buffer = array.request();
19 py::buffer_info variance_buffer = variance_array.request();
20
21 if (buffer.ndim != 2) {
22 throw std::runtime_error("Number of dimensions for array must be 2");
23 }
24 if (buffer.shape[0] != buffer.shape[1]){
25 throw std::runtime_error("Array must be square");
26 }
27 if (!(buffer.shape[0] % 2)) {
28 throw std::runtime_error("Number of pixels along array side must be odd");
29 }
30
31 if (variance_buffer.ndim != buffer.ndim) {
32 throw std::runtime_error("Variance array must have same dimensions as image array");
33 }
34 if (variance_buffer.shape[0] != buffer.shape[0] || variance_buffer.shape[1] != buffer.shape[1]) {
35 throw std::runtime_error("Variance array must have same dimensions as image array");
36 }
37
38 double sig_sq = sig*sig;
39 double t_sq_inv = 1./(t*t);
40 double t_inv = 1./t;
41
42 // fast access to array since we know all the bounds are appropriate
43 auto array_unchecked = array.unchecked<2>();
44 auto variance_array_unchecked = variance_array.unchecked<2>();
45
46 // declare most variables that will be used
47 double flux = 0;
48 double x_offset_sq;
49 double y_offset_sq;
50 double y_component;
51 double y_component_out;
52 double tmp_x_div;
53 double tmp_y_div;
54 double two_sig_sq = 2*sig_sq;
55
56 int half_domain = floor(buffer.shape[0] / 2);
57 int stop = buffer.shape[0];
58
59 // adjust the x and y center to be centered about the middle of the array
60 x_mean -= (double) half_domain;
61 y_mean -= (double) half_domain;
62
63 /*
64 create containers to store the 1 dimensional x calculations, these will be
65 accessed for each loop in y, no need to do the same calculations each time
66 */
67 std::vector<float> x_container;
68 std::vector<float> x_container_out;
69 x_container.reserve(stop);
70 x_container_out.reserve(stop);
71
72 // for weighted variance calculation
73 double variance = 0;
74 double var_norm = 0;
75
76 // calculate the x profile
77 for (int j = 0; j < stop; ++j) {
78 double x_pos = j - half_domain;
79 x_offset_sq = -1*pow(x_pos - x_mean, 2);
80 tmp_x_div = x_offset_sq / two_sig_sq;
81 x_container.push_back(exp(tmp_x_div));
82 x_container_out.push_back(t_inv * exp(tmp_x_div * t_sq_inv));
83 }
84
85 // for each y row, go over the saved x vector accumulating the inner
86 // product.
87 for (int i = 0; i < stop; ++i) {
88 double y_pos = i - half_domain;
89 y_offset_sq = -1*pow(y_pos - y_mean, 2);
90 tmp_y_div = y_offset_sq / two_sig_sq;
91 y_component = exp(tmp_y_div);
92 y_component_out = t_inv * exp(tmp_y_div * t_sq_inv);
93 for (int j = 0; j < stop; ++j){
94 double weight = (y_component*x_container[j] - y_component_out*x_container_out[j]);
95 double weighted_value = weight*array_unchecked(i, j);
96 flux += weighted_value;
97
98 variance += weight*weight*variance_array_unchecked(i, j);
99 var_norm += weight*weight;
100 }
101 }
102
103 // Normalization of the normalized Gaussian filter is 4 * pi * sig**2 * (t**2 + 1)/(t**2 - 1)
104 // We have deliberately not applied the normalization factor 1. / (2 * pi * sig**2) in the
105 // inner product, leaving the normalization term as 2 * (t**2 + 1)/(t**2 - 1) for the
106 // flux, and XXX for the variance.
107 flux *= 2 * (t*t + 1) / (t*t - 1);
108
109 // The variance has been normalized and so requires the full normalization term.
110 variance /= var_norm;
111 variance *= (4 * M_PI * sig_sq * (t*t + 1) ) / (t*t);
112
113 return py::make_tuple(flux, variance);
114}
115
117 wrappers.wrap([](auto &mod) {
118 mod.def("_compensatedGaussianFiltInnerProduct", &compensated_gaussian_filt_inner_product,
119R"doc_string(Calculates flux and variance for a Compensated Gaussian filter.
120
121Parameters
122----------
123array : `np.ndarray` (N, N)
124 Array in which the Gaussian filter will be applied, must be square and
125 odd shaped.
126variance_array : `np.ndarray` (N, N)
127 Variance array in which the Gaussian filter will be applied, must be
128 same shape as array.
129x_mean : `float`
130 The x center position of the object in the reference frame of the array
131 argument.
132y_mean : `float`
133 The y center position of the object in the reference frame of the array
134 argument.
135sig : `float`
136 The std of the Gaussian kernel
137t : `float`
138 The dilating parameter used to calculate background.
139
140Returns
141-------
142result : `tuple` of `float`
143 Compensated Gaussian flux and variance.
144
145Raises
146------
147RuntimeError
148 Raised if the input array is not properly shaped
149
150)doc_string");
151 });
152}
153
154} // namespace base
155} // namespace meas
156} // namespace lsst
afw::table::Key< afw::table::Array< VariancePixelT > > variance
#define M_PI
Definition ListMatch.cc:31
A helper class for subdividing pybind11 module across multiple translation units (i....
Definition python.h:242
void wrap(WrapperCallback function)
Add a set of wrappers without defining a class.
Definition python.h:369
void wrapCalcCompensatedGaussian(WrapperCollection &)
py::tuple compensated_gaussian_filt_inner_product(py::array_t< double, py::array::c_style|py::array::forcecast > &array, py::array_t< double, py::array::c_style|py::array::forcecast > &variance_array, double x_mean, double y_mean, double sig, double t)
T push_back(T... args)
T reserve(T... args)
afw::table::Key< double > weight