LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
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#include <stack>
8#include <queue>
9#include <vector>
10#include <utility> // For std::pair
11#include <stdexcept>
12#include <iostream>
13
14namespace py = pybind11;
15using namespace pybind11::literals;
16
17typedef Eigen::Array<int, 4, 1> Bounds;
18typedef Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixB;
19
20
21// Create a boolean mask `footprint` for all of the pixels that are connected to the pixel
22// located at `i,j` and create the bounding box for the `footprint` in `image`.
23template <typename M>
25 const int start_i,
26 const int start_j,
27 py::EigenDRef<const M> image,
28 py::EigenDRef<MatrixB> unchecked,
29 py::EigenDRef<MatrixB> footprint,
30 Eigen::Ref<Bounds> bounds,
31 const double thresh=0
32){
34 stack.push(std::make_pair(start_i, start_j));
35
36 while (!stack.empty()) {
37 int i, j;
38 std::tie(i, j) = stack.top();
39 stack.pop();
40
41 if (!unchecked(i, j)) {
42 continue;
43 }
44 unchecked(i, j) = false;
45
46 if (image(i, j) > thresh) {
47 footprint(i, j) = true;
48
49 if (i < bounds[0]) {
50 bounds[0] = i;
51 } else if (i > bounds[1]) {
52 bounds[1] = i;
53 }
54 if (j < bounds[2]) {
55 bounds[2] = j;
56 } else if (j > bounds[3]) {
57 bounds[3] = j;
58 }
59
60 if (i > 0 && unchecked(i-1, j)) {
61 stack.push(std::make_pair(i-1, j));
62 }
63 if (i < image.rows() - 1 && unchecked(i+1, j)) {
64 stack.push(std::make_pair(i+1, j));
65 }
66 if (j > 0 && unchecked(i, j-1)) {
67 stack.push(std::make_pair(i, j-1));
68 }
69 if (j < image.cols() - 1 && unchecked(i, j+1)) {
70 stack.push(std::make_pair(i, j+1));
71 }
72 }
73 }
74}
75
76
78template <typename M>
80 py::EigenDRef<const M> image,
81 const std::vector<std::vector<int>>& centers,
82 const double thresh=0
83){
84 const int height = image.rows();
85 const int width = image.cols();
86 MatrixB footprint = MatrixB::Zero(height, width);
88
89 // Seed the queue with peaks
90 for(const auto& center : centers){
91 const int y = center[0];
92 const int x = center[1];
93
94 // Validate center coordinates
95 if (y < 0 || y >= height || x < 0 || x >= width) {
96 throw std::out_of_range("Center coordinates (" + std::to_string(y) + ", " +
97 std::to_string(x) + ") are out of image bounds [0, " +
98 std::to_string(height) + ") x [0, " + std::to_string(width) + ")");
99 }
100
101 if (!footprint(y, x) && image(y, x) > thresh) {
102 footprint(y, x) = true;
103 pixel_queue.emplace(y, x);
104 }
105 }
106
107 // 4-connectivity offsets
108 const std::vector<std::pair<int, int>> offsets = {{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
109
110 // Flood fill
111 while (!pixel_queue.empty()) {
112 auto [i, j] = pixel_queue.front();
113 pixel_queue.pop();
114
115 for (const auto& [di, dj] : offsets) {
116 int ni = i + di;
117 int nj = j + dj;
118 if (ni >= 0 && ni < height && nj >= 0 && nj < width) {
119 if (!footprint(ni, nj) && image(ni, nj) > thresh) {
120 footprint(ni, nj) = true;
121 pixel_queue.emplace(ni, nj);
122 }
123 }
124 }
125 }
126
127 return footprint;
128}
129
130
134class Peak {
135public:
136 Peak(int y, int x, double flux){
137 _y = y;
138 _x = x;
139 _flux = flux;
140 }
141
142 int getY(){
143 return _y;
144 }
145
146 int getX(){
147 return _x;
148 }
149
150 double getFlux(){
151 return _flux;
152 }
153
154
155private:
156 int _y;
157 int _x;
158 double _flux;
159};
160
161
164 return a.getFlux() > b.getFlux();
165}
166
167
168// Get a list of peaks found in an image.
169// To make ut easier to cull peaks that are too close together
170// and ensure that every footprint has at least one peak,
171// this algorithm is meant to be run on a single footprint
172// created by `get_connected_pixels`.
173template <typename M>
175 M& image,
176 const double min_separation,
177 const double peak_thresh,
178 const int y0,
179 const int x0
180){
181 const int height = image.rows();
182 const int width = image.cols();
183
184 std::vector<Peak> peaks;
185
186 for(int i=0; i<height; i++){
187 for(int j=0; j<width; j++){
188 if(image(i, j) < peak_thresh){
189 continue;
190 }
191 if(i > 0 && image(i, j) <= image(i-1, j)){
192 continue;
193 }
194 if(i < height-1 && image(i,j) <= image(i+1, j)){
195 continue;
196 }
197 if(j > 0 && image(i, j) <= image(i, j-1)){
198 continue;
199 }
200 if(j < width-1 && image(i,j) <= image(i, j+1)){
201 continue;
202 }
203
204 if(i > 0 && j > 0 && image(i, j) <= image(i-1, j-1)){
205 continue;
206 }
207 if(i < height-1 && j < width-1 && image(i,j) <= image(i+1, j+1)){
208 continue;
209 }
210 if(i < height-1 && j > 0 && image(i, j) <= image(i+1, j-1)){
211 continue;
212 }
213 if(i > 0 && j < width-1 && image(i,j) <= image(i-1, j+1)){
214 continue;
215 }
216
217 peaks.push_back(Peak(i+y0, j+x0, static_cast<double>(image(i, j))));
218 }
219 }
220
221 if(peaks.empty()){
222 return peaks;
223 }
224
226 std::sort (peaks.begin(), peaks.end(), sortBrightness);
227
228 // Remove peaks within min_separation
229 double min_separation2 = min_separation * min_separation;
230 for (size_t i = 0; i < peaks.size() - 1; ++i) {
231 for (size_t j = i + 1; j < peaks.size();) {
232 Peak *p1 = &peaks[i];
233 Peak *p2 = &peaks[j];
234 double dy = p1->getY() - p2->getY();
235 double dx = p1->getX() - p2->getX();
236 double separation2 = dy*dy + dx*dx;
237 if (separation2 < min_separation2) {
238 peaks.erase(peaks.begin() + j);
239 } else {
240 ++j;
241 }
242 }
243 }
244 return peaks;
245}
246
247
248// A detected footprint
250public:
252 _data = footprint;
253 this->peaks = peaks;
254 _bounds = bounds;
255 }
256
258 return _data;
259 }
260
262
264 return _bounds;
265 }
266
267 void addPeak(Peak peak){
268 peaks.push_back(peak);
269 }
270
271private:
272 MatrixB _data;
273 Bounds _bounds;
274};
275
276
277template <typename M>
279 py::EigenDRef<M> image,
280 py::EigenDRef<MatrixB> footprint
281){
282 const int height = image.rows();
283 const int width = image.cols();
284
285 for(int i=0; i<height; i++){
286 for(int j=0; j<width; j++){
287 if(!footprint(i,j)){
288 image(i,j) = 0;
289 }
290 }
291 }
292}
293
308template <typename M, typename P>
310 py::EigenDRef<const M> image,
311 const double min_separation,
312 const int min_area,
313 const double peak_thresh,
314 const double footprint_thresh,
315 const bool find_peaks=true,
316 const int y0=0,
317 const int x0=0
318){
319 const int height = image.rows();
320 const int width = image.cols();
321
322 std::vector<Footprint> footprints;
323 MatrixB unchecked = MatrixB::Ones(height, width);
324 MatrixB footprint = MatrixB::Zero(height, width);
325
326 for(int i=0; i<height; i++){
327 for(int j=0; j<width; j++){
328 Bounds bounds; bounds << i, i, j, j;
329 get_connected_pixels(i, j, image, unchecked, footprint, bounds, footprint_thresh);
330 int subHeight = bounds[1]-bounds[0]+1;
331 int subWidth = bounds[3]-bounds[2]+1;
332 if(subHeight * subWidth > min_area){
333 MatrixB subFootprint = footprint.block(bounds[0], bounds[2], subHeight, subWidth);
334 int area = subFootprint.count();
335 if(area >= min_area){
336 M patch = image.block(bounds[0], bounds[2], subHeight, subWidth);
337 maskImage<M>(patch, subFootprint);
338 std::vector<Peak> _peaks;
339 if(find_peaks){
340 _peaks = get_peaks(
341 patch,
342 min_separation,
343 peak_thresh,
344 bounds[0] + y0,
345 bounds[2] + x0
346 );
347 }
348 // Only add footprints that have at least one peak above the
349 // minimum peak_thresh.
350 if(!_peaks.empty() || !find_peaks){
351 Bounds trueBounds; trueBounds << bounds[0] + y0,
352 bounds[1] + y0, bounds[2] + x0, bounds[3] + x0;
353 footprints.push_back(Footprint(subFootprint, _peaks, trueBounds));
354 }
355 }
356 }
357 footprint.block(bounds[0], bounds[2], subHeight, subWidth) = MatrixB::Zero(subHeight, subWidth);
358 }
359 }
360 return footprints;
361}
362
363
364
365PYBIND11_MODULE(detect_pybind11, mod) {
366 typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixF;
367 typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> MatrixD;
368
369 mod.doc() = "Fast detection algorithms implemented in C++";
370
371 mod.def("get_connected_pixels", &get_connected_pixels<MatrixB>,
372 "Create a boolean mask for pixels that are connected");
373 mod.def("get_connected_pixels", &get_connected_pixels<MatrixF>,
374 "Create a boolean mask for pixels that are connected");
375 mod.def("get_connected_pixels", &get_connected_pixels<MatrixD>,
376 "Create a boolean mask for pixels that are connected");
377
378 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixB>,
379 "Trim pixels not conencted to a center from a list of centers");
380 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixF>,
381 "Trim pixels not conencted to a center from a list of centers");
382 mod.def("get_connected_multipeak", &get_connected_multipeak<MatrixD>,
383 "Trim pixels not conencted to a center from a list of centers");
384
385 mod.def("get_peaks", &get_peaks<MatrixF>,
386 "Get a list of peaks in a footprint created by get_connected_pixels");
387 mod.def("get_peaks", &get_peaks<MatrixD>,
388 "Get a list of peaks in a footprint created by get_connected_pixels");
389
390 mod.def("get_footprints", &get_footprints<MatrixF, float>,
391 "Create a list of all of the footprints in an image, with their peaks",
392 "image"_a, "min_separation"_a, "min_area"_a, "peak_thresh"_a, "footprint_thresh"_a,
393 "find_peaks"_a=true, "y0"_a=0, "x0"_a=0);
394 mod.def("get_footprints", &get_footprints<MatrixD, double>,
395 "Create a list of all of the footprints in an image, with their peaks",
396 "image"_a, "min_separation"_a, "min_area"_a, "peak_thresh"_a, "footprint_thresh"_a,
397 "find_peaks"_a=true, "y0"_a=0, "x0"_a=0);
398
399 py::class_<Footprint>(mod, "Footprint")
400 .def(py::init<MatrixB, std::vector<Peak>, Bounds>(),
401 "footprint"_a, "peaks"_a, "bounds"_a)
402 .def_property_readonly("data", &Footprint::getFootprint)
403 .def_readwrite("peaks", &Footprint::peaks)
404 .def_property_readonly("bounds", &Footprint::getBounds)
405 .def("add_peak", &Footprint::addPeak);
406
407 py::class_<Peak>(mod, "Peak")
408 .def(py::init<int, int, double>(),
409 "y"_a, "x"_a, "flux"_a)
410 .def_property_readonly("y", &Peak::getY)
411 .def_property_readonly("x", &Peak::getX)
412 .def_property_readonly("flux", &Peak::getFlux);
413}
T begin(T... args)
std::vector< Peak > peaks
Bounds getBounds()
void addPeak(Peak peak)
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...
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
PYBIND11_MODULE(detect_pybind11, mod)
void get_connected_pixels(const int start_i, const int start_j, py::EigenDRef< const M > image, py::EigenDRef< MatrixB > unchecked, py::EigenDRef< MatrixB > footprint, Eigen::Ref< Bounds > bounds, const double thresh=0)
std::vector< Peak > get_peaks(M &image, const double min_separation, const double peak_thresh, const int y0, const int x0)
MatrixB get_connected_multipeak(py::EigenDRef< 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.
Eigen::Matrix< bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixB
std::vector< Footprint > get_footprints(py::EigenDRef< const M > image, const double min_separation, const int min_area, const double peak_thresh, const double footprint_thresh, const bool find_peaks=true, const int y0=0, const int x0=0)
Get all footprints in an image.
void maskImage(py::EigenDRef< M > image, py::EigenDRef< MatrixB > footprint)
T emplace(T... args)
T empty(T... args)
T end(T... args)
T erase(T... args)
T front(T... args)
T make_pair(T... args)
T pop(T... args)
T push_back(T... args)
T push(T... args)
T size(T... args)
T sort(T... args)
T tie(T... args)
T to_string(T... args)
T top(T... args)