LSST Applications g0f08755f38+c89d42e150,g1635faa6d4+b6cf076a36,g1653933729+a8ce1bb630,g1a0ca8cf93+4c08b13bf7,g28da252d5a+f33f8200ef,g29321ee8c0+0187be18b1,g2bbee38e9b+9634bc57db,g2bc492864f+9634bc57db,g2cdde0e794+c2c89b37c4,g3156d2b45e+41e33cbcdc,g347aa1857d+9634bc57db,g35bb328faa+a8ce1bb630,g3a166c0a6a+9634bc57db,g3e281a1b8c+9f2c4e2fc3,g414038480c+077ccc18e7,g41af890bb2+e740673f1a,g5fbc88fb19+17cd334064,g7642f7d749+c89d42e150,g781aacb6e4+a8ce1bb630,g80478fca09+f8b2ab54e1,g82479be7b0+e2bd23ab8b,g858d7b2824+c89d42e150,g9125e01d80+a8ce1bb630,g9726552aa6+10f999ec6a,ga5288a1d22+065360aec4,gacf8899fa4+9553554aa7,gae0086650b+a8ce1bb630,gb58c049af0+d64f4d3760,gbd46683f8f+ac57cbb13d,gc28159a63d+9634bc57db,gcf0d15dbbd+e37acf7834,gda3e153d99+c89d42e150,gda6a2b7d83+e37acf7834,gdaeeff99f8+1711a396fd,ge2409df99d+cb1e6652d6,ge79ae78c31+9634bc57db,gf0baf85859+147a0692ba,gf3967379c6+02b11634a5,w.2024.45
LSST Data Management Base Package
Loading...
Searching...
No Matches
Typedefs | Functions
operators_pybind11.cc File Reference
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
#include <math.h>
#include <algorithm>

Go to the source code of this file.

Typedefs

typedef Eigen::Array< int, Eigen::Dynamic, 1 > IndexVector
 
typedef Eigen::Array< int, 4, 1 > Bounds
 
typedef Eigen::Matrix< bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixB
 

Functions

template<typename T , typename M >
void new_monotonicity (Eigen::Ref< const IndexVector > coord_y, Eigen::Ref< const IndexVector > coord_x, std::vector< M > weights, Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image)
 
template<typename T , typename M , typename V >
void prox_weighted_monotonic (Eigen::Ref< V > flat_img, Eigen::Ref< const M > weights, Eigen::Ref< const IndexVector > offsets, Eigen::Ref< const IndexVector > dist_idx, T const &min_gradient)
 
template<typename M , typename V >
void apply_filter (Eigen::Ref< const M > image, Eigen::Ref< const V > values, Eigen::Ref< const IndexVector > y_start, Eigen::Ref< const IndexVector > y_end, Eigen::Ref< const IndexVector > x_start, Eigen::Ref< const IndexVector > x_end, Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > result)
 
template<typename M >
void get_valid_monotonic_pixels (const int i, const int j, Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > orphans, const double variance, Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds, const double thresh=0)
 
template<typename M , typename P >
void linear_interpolate_invalid_pixels (Eigen::Ref< const IndexVector > row_indices, Eigen::Ref< const IndexVector > column_indices, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked, Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > model, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > orphans, const double variance, bool recursive, Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds)
 
 PYBIND11_MODULE (operators_pybind11, mod)
 

Typedef Documentation

◆ Bounds

typedef Eigen::Array<int, 4, 1> Bounds

Definition at line 11 of file operators_pybind11.cc.

◆ IndexVector

typedef Eigen::Array<int, Eigen::Dynamic, 1> IndexVector

Definition at line 10 of file operators_pybind11.cc.

◆ MatrixB

typedef Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixB

Definition at line 12 of file operators_pybind11.cc.

Function Documentation

◆ apply_filter()

template<typename M , typename V >
void apply_filter ( Eigen::Ref< const M > image,
Eigen::Ref< const V > values,
Eigen::Ref< const IndexVector > y_start,
Eigen::Ref< const IndexVector > y_end,
Eigen::Ref< const IndexVector > x_start,
Eigen::Ref< const IndexVector > x_end,
Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > result )

Definition at line 62 of file operators_pybind11.cc.

70 {
71 result.setZero();
72 for(int n=0; n<values.size(); n++){
73 int rows = image.rows()-y_start(n)-y_end(n);
74 int cols = image.cols()-x_start(n)-x_end(n);
75 result.block(y_start(n), x_start(n), rows, cols) +=
76 values(n) * image.block(y_end(n), x_end(n), rows, cols);
77 }
78}
py::object result
Definition _schema.cc:429

◆ get_valid_monotonic_pixels()

template<typename M >
void get_valid_monotonic_pixels ( const int i,
const int j,
Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image,
Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked,
Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > orphans,
const double variance,
Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds,
const double thresh = 0 )

Definition at line 84 of file operators_pybind11.cc.

93 {
94 // Check the pixel below this one
95 if(i>0 && unchecked(i-1,j)){
96 if(image(i-1,j) <= image(i,j)+variance && image(i-1,j) > thresh){
97 unchecked(i-1,j) = false;
98 orphans(i-1,j) = false;
99 if(i-1 < bounds(0)){
100 bounds(0) = i-1;
101 }
102 get_valid_monotonic_pixels(i-1, j, image, unchecked, orphans, variance, bounds);
103 } else {
104 orphans(i-1,j) = true;
105 }
106 }
107 // Check the pixel above this one
108 if(i < image.rows()-1 && unchecked(i+1,j)){
109 if(image(i+1,j) <= image(i,j)+variance && image(i+1,j) > thresh){
110 unchecked(i+1,j) = false;
111 orphans(i+1,j) = false;
112 if(i+1 > bounds(1)){
113 bounds(1) = i+1;
114 }
115 get_valid_monotonic_pixels(i+1, j, image, unchecked, orphans, variance, bounds);
116 } else {
117 orphans(i+1,j) = true;
118 }
119 }
120 // Check the pixel to the left of this one
121 if(j>0 && unchecked(i,j-1)){
122 if(image(i,j-1) <= image(i,j)+variance && image(i,j-1) > thresh){
123 unchecked(i,j-1) = false;
124 orphans(i,j-1) = false;
125 if(j-1 < bounds(2)){
126 bounds(2) = j-1;
127 }
128 get_valid_monotonic_pixels(i, j-1, image, unchecked, orphans, variance, bounds);
129 } else {
130 orphans(i,j-1) = true;
131 }
132 }
133 // Check the pixel to the right of this one
134 if(j < image.cols()-1 && unchecked(i,j+1)){
135 if(image(i,j+1) <= image(i,j)+variance && image(i,j+1) > thresh){
136 unchecked(i,j+1) = false;
137 orphans(i,j+1) = false;
138 if(j+1 > bounds(3)){
139 bounds(3) = j+1;
140 }
141 get_valid_monotonic_pixels(i, j+1, image, unchecked, orphans, variance, bounds);
142 } else {
143 orphans(i,j+1) = true;
144 }
145 }
146}
afw::table::Key< afw::table::Array< ImagePixelT > > image
afw::table::Key< afw::table::Array< VariancePixelT > > variance
void get_valid_monotonic_pixels(const int i, const int j, Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > orphans, const double variance, Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds, const double thresh=0)

◆ linear_interpolate_invalid_pixels()

template<typename M , typename P >
void linear_interpolate_invalid_pixels ( Eigen::Ref< const IndexVector > row_indices,
Eigen::Ref< const IndexVector > column_indices,
Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked,
Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > model,
Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > orphans,
const double variance,
bool recursive,
Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds )

Definition at line 150 of file operators_pybind11.cc.

159 {
160 for(int n=0; n<row_indices.size(); n++){
161 int i = row_indices(n);
162 int j = column_indices(n);
163 P neighborTotal = 0.0;
164 int validNeighbors = 0;
165 bool uncheckedNeighbors = false;
166
167 if(!unchecked(i,j)){
168 // This pixel has already been updated
169 continue;
170 }
171 // Even if this orphan cannot be updated, we remove it from unchecked
172 // so that it isn't attempted again in future iterations.
173 unchecked(i,j) = false;
174
175 // Check all of the neighboring pixels with negative gradients and
176 // use the maximum value for the interpolation
177 if(i < model.rows()-2 && model(i+2,j) > model(i+1,j)){
178 if(unchecked(i+2, j) || unchecked(i+1, j)){
179 uncheckedNeighbors = true;
180 } else {
181 P grad = model(i+2,j) - model(i+1,j);
182 neighborTotal += model(i+1,j)-grad;
183 validNeighbors += 1;
184 }
185 }
186 if(i > 2 && model(i-2,j) > model(i-1,j)){
187 if(unchecked(i-2, j) || unchecked(i-1, j)){
188 uncheckedNeighbors = true;
189 } else {
190 P grad = model(i-2,j) - model(i-1,j);
191 neighborTotal += model(i-1,j)-grad;
192 validNeighbors += 1;
193 }
194 }
195 if(j < model.cols()-2 && model(i,j+2) > model(i,j+1)){
196 if(unchecked(i,j+2), unchecked(i,j+1)){
197 uncheckedNeighbors = true;
198 } else {
199 P grad = model(i,j+2) - model(i,j+1);
200 neighborTotal += model(i,j+1)-grad;
201 validNeighbors += 1;
202 }
203 }
204 if(j > 2 && model(i,j-2) > model(i,j-1)){
205 if(unchecked(i, j-2), unchecked(i,j-1)){
206 uncheckedNeighbors = true;
207 } else {
208 P grad = model(i,j-2) - model(i,j-1);
209 neighborTotal += model(i,j-1)-grad;
210 validNeighbors += 1;
211 }
212 }
213 // If the non-monotonic pixel was updated then update the
214 // model with the interpolated value and search for more monotonic pixels
215 if(neighborTotal > 0){
216 // Update the model and orphan status _before_ checking neighbors
217 model(i,j) = neighborTotal / validNeighbors;
218 orphans(i,j) = false;
219
220 if(i < bounds(0)){
221 bounds(0) = i;
222 } else if(i > bounds(1)){
223 bounds(1) = i;
224 }
225 if(j < bounds(2)){
226 bounds(2) = j;
227 } else if(j > bounds(3)){
228 bounds(3) = j;
229 }
230 if(recursive){
231 get_valid_monotonic_pixels(i, j, model, unchecked, orphans, variance, bounds);
232 } else {
233 if(unchecked(i-1,j)){
234 orphans(i-1,j) = true;
235 }
236 if(unchecked(i+1,j)){
237 orphans(i+1,j) = true;
238 }
239 if(unchecked(i,j-1)){
240 orphans(i,j-1) = true;
241 }
242 if(unchecked(i,j+1)){
243 orphans(i,j+1) = true;
244 }
245 }
246 } else if(uncheckedNeighbors){
247 unchecked(i,j) = false;
248 } else {
249 // This is still an orphan, but we checked it so it won't be iterated on again.
250 orphans(i,j) = true;
251 model(i,j) = 0;
252 }
253 }
254}
model(data, psfmodels, sources)

◆ new_monotonicity()

template<typename T , typename M >
void new_monotonicity ( Eigen::Ref< const IndexVector > coord_y,
Eigen::Ref< const IndexVector > coord_x,
std::vector< M > weights,
Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image )

Definition at line 15 of file operators_pybind11.cc.

20 {
21 for(int n=0; n<coord_x.size(); n++){
22 int px = coord_x[n];
23 int py = coord_y[n];
24
25 T ref_flux = 0;
26 for(int i=0; i<3; i++){
27 for(int j=0; j<3; j++){
28 int weight_index = 3*i + j;
29 ref_flux += image(py + i, px + j) * weights[weight_index](py, px);
30 }
31 }
32 image(py + 1, px + 1) = std::min(image(py+1, px+1), ref_flux);
33 }
34}
T min(T... args)

◆ prox_weighted_monotonic()

template<typename T , typename M , typename V >
void prox_weighted_monotonic ( Eigen::Ref< V > flat_img,
Eigen::Ref< const M > weights,
Eigen::Ref< const IndexVector > offsets,
Eigen::Ref< const IndexVector > dist_idx,
T const & min_gradient )

Definition at line 37 of file operators_pybind11.cc.

44 {
45 // Start at the center of the image and set each pixel to the minimum
46 // between itself and its reference pixel (which is closer to the peak)
47 for(int d=0; d<dist_idx.size(); d++){
48 int didx = dist_idx(d);
49 T ref_flux = 0;
50 for(int i=0; i<offsets.size(); i++){
51 if(weights(i,didx)>0){
52 int nidx = offsets[i] + didx;
53 ref_flux += flat_img(nidx) * weights(i, didx);
54 }
55 }
56 flat_img(didx) = std::min(flat_img(didx), ref_flux*(1-min_gradient));
57 }
58}

◆ PYBIND11_MODULE()

PYBIND11_MODULE ( operators_pybind11 ,
mod  )

Definition at line 256 of file operators_pybind11.cc.

257{
258 mod.doc() = "operators_pybind11", "Fast proximal operators";
259
260 typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixF;
261 typedef Eigen::Matrix<float, Eigen::Dynamic, 1> VectorF;
262 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixD;
263 typedef Eigen::Matrix<double, Eigen::Dynamic, 1> VectorD;
264
265 mod.def("new_monotonicity", &new_monotonicity<float, MatrixF>, "Weighted Monotonic Proximal Operator");
266 mod.def("new_monotonicity", &new_monotonicity<double, MatrixD>, "Weighted Monotonic Proximal Operator");
267
268 mod.def("prox_weighted_monotonic", &prox_weighted_monotonic<float, MatrixF, VectorF>,
269 "Weighted Monotonic Proximal Operator");
270 mod.def("prox_weighted_monotonic", &prox_weighted_monotonic<double, MatrixD, VectorD>,
271 "Weighted Monotonic Proximal Operator");
272
273 mod.def("apply_filter", &apply_filter<MatrixF, VectorF>, "Apply a filter to a 2D image");
274 mod.def("apply_filter", &apply_filter<MatrixD, VectorD>, "Apply a filter to a 2D image");
275
276 mod.def("get_valid_monotonic_pixels", &get_valid_monotonic_pixels<MatrixF>,
277 "Create a boolean mask for pixels that are monotonic from the center along some path");
278 mod.def("get_valid_monotonic_pixels", &get_valid_monotonic_pixels<MatrixD>,
279 "Create a boolean mask for pixels that are monotonic from the center along some path");
280
281 mod.def("linear_interpolate_invalid_pixels", &linear_interpolate_invalid_pixels<MatrixF, float>,
282 "Fill in non-monotonic pixels by interpolating based on the gradients of its neighbors");
283 mod.def("linear_interpolate_invalid_pixels", &linear_interpolate_invalid_pixels<MatrixD, double>,
284 "Fill in non-monotonic pixels by interpolating based on the gradients of its neighbors");
285}