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
detect_pybind11.cc
Go to the documentation of this file.
1#include <pybind11/pybind11.h>
2#include <pybind11/numpy.h>
3#include <pybind11/stl.h>
4#include <pybind11/eigen.h>
5#include <math.h>
6#include <algorithm>
7
8namespace py = pybind11;
9using namespace pybind11::literals;
10
11typedef Eigen::Array<int, 4, 1> Bounds;
12typedef Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixB;
13
14
15// Create a boolean mask `footprint` for all of the pixels that are connected to the pixel
16// located at `i,j` and create the bounding box for the `footprint` in `image`.
17template <typename M>
19 const int i,
20 const int j,
21 Eigen::Ref<const M> image,
22 Eigen::Ref<MatrixB, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> unchecked,
23 Eigen::Ref<MatrixB, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> footprint,
24 Eigen::Ref<Bounds, 0, Eigen::Stride<4, 1>> bounds,
25 const double thresh=0
26){
27 if(not unchecked(i,j)){
28 return;
29 }
30 unchecked(i,j) = false;
31
32 if(image(i,j) > thresh){
33 footprint(i,j) = true;
34
35 if(i < bounds[0]){
36 bounds[0] = i;
37 } else if(i > bounds[1]){
38 bounds[1] = i;
39 }
40 if(j < bounds[2]){
41 bounds[2] = j;
42 } else if(j > bounds[3]){
43 bounds[3] = j;
44 }
45
46 if(i > 0){
47 get_connected_pixels(i-1, j, image, unchecked, footprint, bounds, thresh);
48 }
49 if(i < image.rows()-1){
50 get_connected_pixels(i+1, j, image, unchecked, footprint, bounds, thresh);
51 }
52 if(j > 0){
53 get_connected_pixels(i, j-1, image, unchecked, footprint, bounds, thresh);
54 }
55 if(j < image.cols()-1){
56 get_connected_pixels(i, j+1, image, unchecked, footprint, bounds, thresh);
57 }
58 }
59}
60
61
63template <typename M>
65 Eigen::Ref<const M> image,
66 const std::vector<std::vector<int>> centers,
67 const double thresh=0
68){
69 const int height = image.rows();
70 const int width = image.cols();
71 MatrixB unchecked = MatrixB::Ones(height, width);
72 MatrixB footprint = MatrixB::Zero(height, width);
73
74 for(auto center=begin(centers); center!=end(centers); ++center){
75 const int y = (*center)[0];
76 const int x = (*center)[1];
77 Bounds bounds; bounds << y, y, x, x;
78 get_connected_pixels(y, x, image, unchecked, footprint, bounds, thresh);
79 }
80
81 return footprint;
82}
83
84
88class Peak {
89public:
90 Peak(int y, int x, double flux){
91 _y = y;
92 _x = x;
93 _flux = flux;
94 }
95
96 int getY(){
97 return _y;
98 }
99
100 int getX(){
101 return _x;
102 }
103
104 double getFlux(){
105 return _flux;
106 }
107
108
109private:
110 int _y;
111 int _x;
112 double _flux;
113};
114
115
118 return a.getFlux() > b.getFlux();
119}
120
121
122// Get a list of peaks found in an image.
123// To make ut easier to cull peaks that are too close together
124// and ensure that every footprint has at least one peak,
125// this algorithm is meant to be run on a single footprint
126// created by `get_connected_pixels`.
127template <typename M>
129 M& image,
130 const double min_separation,
131 const int y0,
132 const int x0
133){
134 const int height = image.rows();
135 const int width = image.cols();
136
137 std::vector<Peak> peaks;
138
139 for(int i=0; i<height; i++){
140 for(int j=0; j<width; j++){
141 if(i > 0 && image(i, j) <= image(i-1, j)){
142 continue;
143 }
144 if(i < height-1 && image(i,j) <= image(i+1, j)){
145 continue;
146 }
147 if(j > 0 && image(i, j) <= image(i, j-1)){
148 continue;
149 }
150 if(j < width-1 && image(i,j) <= image(i, j+1)){
151 continue;
152 }
153
154 if(i > 0 && j > 0 && image(i, j) <= image(i-1, j-1)){
155 continue;
156 }
157 if(i < height-1 && j < width-1 && image(i,j) <= image(i+1, j+1)){
158 continue;
159 }
160 if(i < height-1 && j > 0 && image(i, j) <= image(i+1, j-1)){
161 continue;
162 }
163 if(i > 0 && j < width-1 && image(i,j) <= image(i-1, j+1)){
164 continue;
165 }
166
167 peaks.push_back(Peak(i+y0, j+x0, static_cast<double>(image(i, j))));
168 }
169 }
170
171 assert(peaks.size() > 0);
172
174 std::sort (peaks.begin(), peaks.end(), sortBrightness);
175
176 // Remove peaks within min_separation
177 double min_separation2 = min_separation * min_separation;
178 int i = 0;
179 while (i < peaks.size()-1){
180 int j = i+1;
181 Peak *p1 = &peaks[i];
182 while (j < peaks.size()){
183 Peak *p2 = &peaks[j];
184 double dy = p1->getY()-p2->getY();
185 double dx = p1->getX()-p2->getX();
186 double separation2 = dy*dy + dx*dx;
187 if(separation2 < min_separation2){
188 peaks.erase(peaks.begin()+j);
189 i--;
190 }
191 j++;
192 }
193 i++;
194 }
195
196 assert(peaks.size() > 0);
197
198 return peaks;
199}
200
201
202// A detected footprint
204public:
206 _data = footprint;
207 this->peaks = peaks;
208 _bounds = bounds;
209 }
210
212 return _data;
213 }
214
216
218 return _bounds;
219 }
220
221private:
222 MatrixB _data;
223 Bounds _bounds;
224};
225
226
227template <typename M>
229 Eigen::Ref<M, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> image,
230 Eigen::Ref<MatrixB, 0, Eigen::Stride<Eigen::Dynamic, Eigen::Dynamic>> footprint
231){
232 const int height = image.rows();
233 const int width = image.cols();
234
235 for(int i=0; i<height; i++){
236 for(int j=0; j<width; j++){
237 if(!footprint(i,j)){
238 image(i,j) = 0;
239 }
240 }
241 }
242}
243
244
245template <typename M, typename P>
247 Eigen::Ref<const M> image,
248 const double min_separation,
249 const int min_area,
250 const double thresh,
251 const bool find_peaks=true
252){
253 const int height = image.rows();
254 const int width = image.cols();
255
256 std::vector<Footprint> footprints;
257 MatrixB unchecked = MatrixB::Ones(height, width);
258 MatrixB footprint = MatrixB::Zero(height, width);
259
260 for(int i=0; i<height; i++){
261 for(int j=0; j<width; j++){
262 Bounds bounds; bounds << i, i, j, j;
263 get_connected_pixels(i, j, image, unchecked, footprint, bounds, thresh);
264 int subHeight = bounds[1]-bounds[0]+1;
265 int subWidth = bounds[3]-bounds[2]+1;
266 if(subHeight * subWidth > min_area){
267 MatrixB subFootprint = footprint.block(bounds[0], bounds[2], subHeight, subWidth);
268 int area = subFootprint.count();
269 if(area >= min_area){
270 M patch = image.block(bounds[0], bounds[2], subHeight, subWidth);
271 maskImage<M>(patch, subFootprint);
272 std::vector<Peak> _peaks;
273 if(find_peaks){
274 _peaks = get_peaks(
275 patch,
276 min_separation,
277 bounds[0],
278 bounds[2]
279 );
280 }
281
282 footprints.push_back(Footprint(subFootprint, _peaks, bounds));
283 }
284 }
285 footprint.block(bounds[0], bounds[2], subHeight, subWidth) = MatrixB::Zero(subHeight, subWidth);
286 }
287 }
288 return footprints;
289}
290
291
292
293PYBIND11_MODULE(detect_pybind11, mod) {
294 typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixF;
295 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixD;
296
297 mod.doc() = "Fast detection algorithms implemented in C++";
298
299 mod.def("get_connected_pixels", &get_connected_pixels<MatrixB>,
300 "Create a boolean mask for pixels that are connected");
301 mod.def("get_connected_pixels", &get_connected_pixels<MatrixF>,
302 "Create a boolean mask for pixels that are connected");
303 mod.def("get_connected_pixels", &get_connected_pixels<MatrixD>,
304 "Create a boolean mask for pixels that are connected");
305
306 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixB>,
307 "Trim pixels not conencted to a center from a list of centers");
308 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixF>,
309 "Trim pixels not conencted to a center from a list of centers");
310 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixD>,
311 "Trim pixels not conencted to a center from a list of centers");
312
313 mod.def("get_peaks", &get_peaks<MatrixF>,
314 "Get a list of peaks in a footprint created by get_connected_pixels");
315 mod.def("get_peaks", &get_peaks<MatrixD>,
316 "Get a list of peaks in a footprint created by get_connected_pixels");
317
318 mod.def("get_footprints", &get_footprints<MatrixF, float>,
319 "Create a list of all of the footprints in an image, with their peaks"
320 "image"_a, "min_separation"_a, "min_area"_a, "thresh"_a, "find_peaks"_a);
321 mod.def("get_footprints", &get_footprints<MatrixD, double>,
322 "Create a list of all of the footprints in an image, with their peaks"
323 "image"_a, "min_separation"_a, "min_area"_a, "thresh"_a, "find_peaks"_a);
324
325 py::class_<Footprint>(mod, "Footprint")
326 .def(py::init<MatrixB, std::vector<Peak>, Bounds>(),
327 "footprint"_a, "peaks"_a, "bounds"_a)
328 .def_property_readonly("data", &Footprint::getFootprint)
329 .def_readwrite("peaks", &Footprint::peaks)
330 .def_property_readonly("bounds", &Footprint::getBounds);
331
332 py::class_<Peak>(mod, "Peak")
333 .def(py::init<int, int, double>(),
334 "y"_a, "x"_a, "flux"_a)
335 .def_property_readonly("y", &Peak::getY)
336 .def_property_readonly("x", &Peak::getX)
337 .def_property_readonly("flux", &Peak::getFlux);
338}
int end
afw::table::Key< afw::table::Array< ImagePixelT > > image
int y
Definition SpanSet.cc:48
table::Key< int > b
T begin(T... args)
std::vector< Peak > peaks
Bounds getBounds()
Footprint(MatrixB footprint, std::vector< Peak > peaks, Bounds bounds)
MatrixB getFootprint()
A Peak in a Footprint This class is meant to keep track of both the position and flux at the location...
int getY()
Peak(int y, int x, double flux)
double getFlux()
bool sortBrightness(Peak a, Peak b)
Sort two peaks, placing the brightest peak first.
Eigen::Array< int, 4, 1 > Bounds
MatrixB get_connected_multipeak(Eigen::Ref< const M > image, const std::vector< std::vector< int > > centers, const double thresh=0)
Proximal operator to trim pixels not connected to one of the source centers.
PYBIND11_MODULE(detect_pybind11, mod)
void maskImage(Eigen::Ref< M, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > image, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > footprint)
std::vector< Peak > get_peaks(M &image, const double min_separation, const int y0, const int x0)
std::vector< Footprint > get_footprints(Eigen::Ref< const M > image, const double min_separation, const int min_area, const double thresh, const bool find_peaks=true)
Eigen::Matrix< bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixB
void get_connected_pixels(const int i, const int j, Eigen::Ref< const M > image, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > unchecked, Eigen::Ref< MatrixB, 0, Eigen::Stride< Eigen::Dynamic, Eigen::Dynamic > > footprint, Eigen::Ref< Bounds, 0, Eigen::Stride< 4, 1 > > bounds, const double thresh=0)
T end(T... args)
T erase(T... args)
T push_back(T... args)
T size(T... args)
T sort(T... args)