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 File Reference
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
#include <math.h>
#include <algorithm>
#include <stack>
#include <queue>
#include <vector>
#include <utility>
#include <stdexcept>
#include <iostream>

Go to the source code of this file.

Classes

class  Peak
 A Peak in a Footprint This class is meant to keep track of both the position and flux at the location of a maximum in a Footprint. More...
 
class  Footprint
 

Typedefs

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

Functions

template<typename M>
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)
 
template<typename M>
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.
 
bool sortBrightness (Peak a, Peak b)
 Sort two peaks, placing the brightest peak first.
 
template<typename M>
std::vector< Peakget_peaks (M &image, const double min_separation, const double peak_thresh, const int y0, const int x0)
 
template<typename M>
void maskImage (py::EigenDRef< M > image, py::EigenDRef< MatrixB > footprint)
 
template<typename M, typename P>
std::vector< Footprintget_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.
 
 PYBIND11_MODULE (detect_pybind11, mod)
 

Typedef Documentation

◆ Bounds

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

Definition at line 17 of file detect_pybind11.cc.

◆ MatrixB

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

Definition at line 18 of file detect_pybind11.cc.

Function Documentation

◆ get_connected_multipeak()

template<typename M>
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.

Definition at line 79 of file detect_pybind11.cc.

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}
Eigen::Matrix< bool, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixB
T emplace(T... args)
T empty(T... args)
T front(T... args)
T pop(T... args)
T to_string(T... args)

◆ get_connected_pixels()

template<typename M>
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 )

Definition at line 24 of file detect_pybind11.cc.

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}
T make_pair(T... args)
T push(T... args)
T tie(T... args)
T top(T... args)

◆ get_footprints()

template<typename M, typename P>
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.

Parameters
imageThe image to search for footprints
min_separationThe minimum separation (in pixels) between peaks in a footprint
min_areaThe minimum area of a footprint in pixels
peak_threshThe minimum flux of a peak to be detected.
footprint_threshThe minimum flux of a pixel to be included in a footprint
find_peaksIf True, find peaks in each footprint
y0The y-coordinate of the top-left corner of the image
x0The x-coordinate of the top-left corner of the image
Returns
: A list of Footprints

Definition at line 309 of file detect_pybind11.cc.

318 {
319 const int height = image.rows();
320 const int width = image.cols();
321
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}
Eigen::Array< int, 4, 1 > Bounds
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)
void maskImage(py::EigenDRef< M > image, py::EigenDRef< MatrixB > footprint)

◆ get_peaks()

template<typename M>
std::vector< Peak > get_peaks ( M & image,
const double min_separation,
const double peak_thresh,
const int y0,
const int x0 )

Sort the peaks in the footprint so that the brightest are first

Definition at line 174 of file detect_pybind11.cc.

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}
T begin(T... args)
A Peak in a Footprint This class is meant to keep track of both the position and flux at the location...
bool sortBrightness(Peak a, Peak b)
Sort two peaks, placing the brightest peak first.
T end(T... args)
T erase(T... args)
T push_back(T... args)
T size(T... args)
T sort(T... args)

◆ maskImage()

template<typename M>
void maskImage ( py::EigenDRef< M > image,
py::EigenDRef< MatrixB > footprint )

Definition at line 278 of file detect_pybind11.cc.

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}

◆ PYBIND11_MODULE()

PYBIND11_MODULE ( detect_pybind11 ,
mod  )

Definition at line 365 of file detect_pybind11.cc.

365 {
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}
std::vector< Peak > peaks
Bounds getBounds()
void addPeak(Peak peak)
MatrixB getFootprint()
double getFlux()
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.
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.

◆ sortBrightness()

bool sortBrightness ( Peak a,
Peak b )

Sort two peaks, placing the brightest peak first.

Definition at line 163 of file detect_pybind11.cc.

163 {
164 return a.getFlux() > b.getFlux();
165}