From a579ebbe1b596ca30759911633197be9f96c8c3c Mon Sep 17 00:00:00 2001 From: Jeremy Kubica <104161096+jeremykubica@users.noreply.github.com> Date: Thu, 12 Oct 2023 09:59:42 -0400 Subject: [PATCH] Mechanical application of clang-format --- src/kbmod/search/bindings.cpp | 47 ++- src/kbmod/search/common.h | 217 +++++------ src/kbmod/search/filtering.cpp | 39 +- src/kbmod/search/filtering.h | 11 +- src/kbmod/search/image_stack.cpp | 143 ++++--- src/kbmod/search/image_stack.h | 9 +- src/kbmod/search/layered_image.cpp | 220 +++++------ src/kbmod/search/layered_image.h | 20 +- src/kbmod/search/psf.cpp | 186 +++++---- src/kbmod/search/psf.h | 9 +- src/kbmod/search/raw_image.cpp | 596 ++++++++++++++--------------- src/kbmod/search/raw_image.h | 25 +- src/kbmod/search/stack_search.cpp | 564 +++++++++++++-------------- src/kbmod/search/stack_search.h | 13 +- 14 files changed, 1027 insertions(+), 1072 deletions(-) diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index cec9e510..46f5e342 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -13,29 +13,28 @@ using pp = search::PixelPos; using std::to_string; - PYBIND11_MODULE(search, m) { - m.attr("KB_NO_DATA") = pybind11::float_(search::NO_DATA); - m.attr("HAS_GPU") = pybind11::bool_(search::HAVE_GPU); - py::enum_(m, "StampType") - .value("STAMP_SUM", search::StampType::STAMP_SUM) - .value("STAMP_MEAN", search::StampType::STAMP_MEAN) - .value("STAMP_MEDIAN", search::StampType::STAMP_MEDIAN) - .export_values(); - search::psf_bindings(m); - search::raw_image_bindings(m); - search::layered_image_bindings(m); - search::image_stack_bindings(m); - search::stack_search_bindings(m); - search::trajectory_bindings(m); - search::pixel_pos_bindings(m); - search::image_moments_bindings(m); - search::stamp_parameters_bindings(m); - // Functions from raw_image.cpp - m.def("create_median_image", &search::create_median_image); - m.def("create_summed_image", &search::create_summed_image); - m.def("create_mean_image", &search::create_mean_image); - // Functions from filtering.cpp - m.def("sigmag_filtered_indices", &search::sigmaGFilteredIndices); - m.def("calculate_likelihood_psi_phi", &search::calculateLikelihoodFromPsiPhi); + m.attr("KB_NO_DATA") = pybind11::float_(search::NO_DATA); + m.attr("HAS_GPU") = pybind11::bool_(search::HAVE_GPU); + py::enum_(m, "StampType") + .value("STAMP_SUM", search::StampType::STAMP_SUM) + .value("STAMP_MEAN", search::StampType::STAMP_MEAN) + .value("STAMP_MEDIAN", search::StampType::STAMP_MEDIAN) + .export_values(); + search::psf_bindings(m); + search::raw_image_bindings(m); + search::layered_image_bindings(m); + search::image_stack_bindings(m); + search::stack_search_bindings(m); + search::trajectory_bindings(m); + search::pixel_pos_bindings(m); + search::image_moments_bindings(m); + search::stamp_parameters_bindings(m); + // Functions from raw_image.cpp + m.def("create_median_image", &search::create_median_image); + m.def("create_summed_image", &search::create_summed_image); + m.def("create_mean_image", &search::create_mean_image); + // Functions from filtering.cpp + m.def("sigmag_filtered_indices", &search::sigmaGFilteredIndices); + m.def("calculate_likelihood_psi_phi", &search::calculateLikelihoodFromPsiPhi); } diff --git a/src/kbmod/search/common.h b/src/kbmod/search/common.h index 453edf86..8412af72 100644 --- a/src/kbmod/search/common.h +++ b/src/kbmod/search/common.h @@ -1,33 +1,31 @@ #ifndef COMMON_H_ #define COMMON_H_ - #include #include "pydocs/common_docs.h" - namespace search { #ifdef HAVE_CUDA - constexpr bool HAVE_GPU = true; +constexpr bool HAVE_GPU = true; #else - constexpr bool HAVE_GPU = false; +constexpr bool HAVE_GPU = false; #endif - constexpr unsigned int MAX_KERNEL_RADIUS = 15; - constexpr unsigned short MAX_STAMP_EDGE = 64; - constexpr unsigned short CONV_THREAD_DIM = 32; - constexpr unsigned short THREAD_DIM_X = 128; - constexpr unsigned short THREAD_DIM_Y = 2; - constexpr unsigned short RESULTS_PER_PIXEL = 8; - constexpr float NO_DATA = -9999.0; - - enum StampType { STAMP_SUM = 0, STAMP_MEAN, STAMP_MEDIAN }; - - /* - * Data structure to represent an objects trajectory - * through a stack of images - */ - struct Trajectory { +constexpr unsigned int MAX_KERNEL_RADIUS = 15; +constexpr unsigned short MAX_STAMP_EDGE = 64; +constexpr unsigned short CONV_THREAD_DIM = 32; +constexpr unsigned short THREAD_DIM_X = 128; +constexpr unsigned short THREAD_DIM_Y = 2; +constexpr unsigned short RESULTS_PER_PIXEL = 8; +constexpr float NO_DATA = -9999.0; + +enum StampType { STAMP_SUM = 0, STAMP_MEAN, STAMP_MEDIAN }; + +/* + * Data structure to represent an objects trajectory + * through a stack of images + */ +struct Trajectory { // Trajectory velocities float vx; float vy; @@ -43,45 +41,34 @@ namespace search { // I can't believe string::format is not a thing until C++ 20 const std::string to_string() const { - return "lh: " + std::to_string(lh) + - " flux: " + std::to_string(flux) + - " x: " + std::to_string(x) + - " y: " + std::to_string(y) + - " vx: " + std::to_string(vx) + - " vy: " + std::to_string(vy) + - " obs_count: " + std::to_string(obs_count); + return "lh: " + std::to_string(lh) + " flux: " + std::to_string(flux) + " x: " + std::to_string(x) + + " y: " + std::to_string(y) + " vx: " + std::to_string(vx) + " vy: " + std::to_string(vy) + + " obs_count: " + std::to_string(obs_count); } // returns a yaml-compliant string const std::string to_yaml() const { - return "{lh: " + std::to_string(lh) + - ", flux: " + std::to_string(flux) + - ", x: " + std::to_string(x) + - ", y: " + std::to_string(y) + - ", vx: " + std::to_string(vx) + - ", vy: " + std::to_string(vy) + - ", obs_count: " + std::to_string(obs_count) - +"}"; + return "{lh: " + std::to_string(lh) + ", flux: " + std::to_string(flux) + + ", x: " + std::to_string(x) + ", y: " + std::to_string(y) + ", vx: " + std::to_string(vx) + + ", vy: " + std::to_string(vy) + ", obs_count: " + std::to_string(obs_count) + "}"; } - }; +}; - // The position (in pixels) of a trajectory. - struct PixelPos { +// The position (in pixels) of a trajectory. +struct PixelPos { float x; float y; - const std::string to_string() const { - return "x: " + std::to_string(x) + " y: " + std::to_string(y); - } + const std::string to_string() const { return "x: " + std::to_string(x) + " y: " + std::to_string(y); } const std::string to_yaml() const { - return "{x: " + std::to_string(x) + " y: " + std::to_string(y) + "}"; + return "{x: " + std::to_string(x) + " y: " + std::to_string(y) + "}"; } - }; +}; - /* The parameters to use for the on device search. */ +/* The parameters to use for the on device search. */ - struct SearchParameters { +struct SearchParameters { // Basic filtering paramets. int min_observations; float min_lh; @@ -105,24 +92,24 @@ namespace search { // Provide debugging output. bool debug; - }; +}; - struct scaleParameters { +struct scaleParameters { float min_val; float max_val; float scale; - }; +}; - // Search data on a per-image basis. - struct PerImageData { +// Search data on a per-image basis. +struct PerImageData { int num_images = 0; - float* image_times = nullptr; - scaleParameters* psi_params = nullptr; - scaleParameters* phi_params = nullptr; - }; + float *image_times = nullptr; + scaleParameters *psi_params = nullptr; + scaleParameters *phi_params = nullptr; +}; - struct StampParameters { +struct StampParameters { int radius = 10; StampType stamp_type = STAMP_SUM; bool do_filtering = false; @@ -138,89 +125,83 @@ namespace search { float m11_limit; float m02_limit; float m20_limit; - }; +}; - // Basic image moments use for analysis. - struct ImageMoments { +// Basic image moments use for analysis. +struct ImageMoments { float m00; float m01; float m10; float m11; float m02; float m20; - }; +}; #ifdef Py_PYTHON_H - namespace py = pybind11; +namespace py = pybind11; - static void trajectory_bindings(py::module &m) { +static void trajectory_bindings(py::module &m) { using tj = Trajectory; py::class_(m, "Trajectory", pydocs::DOC_Trajectory) - .def(py::init<>()) - .def_readwrite("vx", &tj::vx) - .def_readwrite("vy", &tj::vy) - .def_readwrite("lh", &tj::lh) - .def_readwrite("flux", &tj::flux) - .def_readwrite("x", &tj::x) - .def_readwrite("y", &tj::y) - .def_readwrite("obs_count", &tj::obs_count) - .def("__repr__", [](const tj &t) { return "Trajectory(" + t.to_string() + ")"; }) - .def("__str__", &tj::to_string) - .def(py::pickle( - [] (const tj &p) { // __getstate__ + .def(py::init<>()) + .def_readwrite("vx", &tj::vx) + .def_readwrite("vy", &tj::vy) + .def_readwrite("lh", &tj::lh) + .def_readwrite("flux", &tj::flux) + .def_readwrite("x", &tj::x) + .def_readwrite("y", &tj::y) + .def_readwrite("obs_count", &tj::obs_count) + .def("__repr__", [](const tj &t) { return "Trajectory(" + t.to_string() + ")"; }) + .def("__str__", &tj::to_string) + .def(py::pickle( + [](const tj &p) { // __getstate__ return py::make_tuple(p.vx, p.vy, p.lh, p.flux, p.x, p.y, p.obs_count); - }, - [] (py::tuple t) { // __setstate__ - if (t.size() != 7) - throw std::runtime_error("Invalid state!"); - tj trj = { - t[0].cast(), t[1].cast(), t[2].cast(), - t[3].cast(), t[4].cast(), t[5].cast(), - t[6].cast() - }; + }, + [](py::tuple t) { // __setstate__ + if (t.size() != 7) throw std::runtime_error("Invalid state!"); + tj trj = {t[0].cast(), t[1].cast(), t[2].cast(), + t[3].cast(), t[4].cast(), t[5].cast(), + t[6].cast()}; return trj; - }) - ); - } + })); +} - static void pixel_pos_bindings(py::module &m) { +static void pixel_pos_bindings(py::module &m) { py::class_(m, "PixelPos", pydocs::DOC_PixelPos) - .def(py::init<>()) - .def_readwrite("x", &PixelPos::x) - .def_readwrite("y", &PixelPos::y) - .def("__repr__", [] (const PixelPos &p) { - return "PixelPos(" + p.to_string() + ")"; - }) - .def("__str__", &PixelPos::to_string); - } - - static void image_moments_bindings(py::module &m) { + .def(py::init<>()) + .def_readwrite("x", &PixelPos::x) + .def_readwrite("y", &PixelPos::y) + .def("__repr__", [](const PixelPos &p) { return "PixelPos(" + p.to_string() + ")"; }) + .def("__str__", &PixelPos::to_string); +} + +static void image_moments_bindings(py::module &m) { py::class_(m, "ImageMoments", pydocs::DOC_ImageMoments) - .def(py::init<>()) - .def_readwrite("m00", &ImageMoments::m00) - .def_readwrite("m01", &ImageMoments::m01) - .def_readwrite("m10", &ImageMoments::m10) - .def_readwrite("m11", &ImageMoments::m11) - .def_readwrite("m02", &ImageMoments::m02) - .def_readwrite("m20", &ImageMoments::m20); - } - - static void stamp_parameters_bindings(py::module &m) { + .def(py::init<>()) + .def_readwrite("m00", &ImageMoments::m00) + .def_readwrite("m01", &ImageMoments::m01) + .def_readwrite("m10", &ImageMoments::m10) + .def_readwrite("m11", &ImageMoments::m11) + .def_readwrite("m02", &ImageMoments::m02) + .def_readwrite("m20", &ImageMoments::m20); +} + +static void stamp_parameters_bindings(py::module &m) { py::class_(m, "StampParameters", pydocs::DOC_StampParameters) - .def(py::init<>()) - .def_readwrite("radius", &StampParameters::radius) - .def_readwrite("stamp_type", &StampParameters::stamp_type) - .def_readwrite("do_filtering", &StampParameters::do_filtering) - .def_readwrite("center_thresh", &StampParameters::center_thresh) - .def_readwrite("peak_offset_x", &StampParameters::peak_offset_x) - .def_readwrite("peak_offset_y", &StampParameters::peak_offset_y) - .def_readwrite("m01_limit", &StampParameters::m01_limit) - .def_readwrite("m10_limit", &StampParameters::m10_limit) - .def_readwrite("m11_limit", &StampParameters::m11_limit) - .def_readwrite("m02_limit", &StampParameters::m02_limit) - .def_readwrite("m20_limit", &StampParameters::m20_limit); - } + .def(py::init<>()) + .def_readwrite("radius", &StampParameters::radius) + .def_readwrite("stamp_type", &StampParameters::stamp_type) + .def_readwrite("do_filtering", &StampParameters::do_filtering) + .def_readwrite("center_thresh", &StampParameters::center_thresh) + .def_readwrite("peak_offset_x", &StampParameters::peak_offset_x) + .def_readwrite("peak_offset_y", &StampParameters::peak_offset_y) + .def_readwrite("m01_limit", &StampParameters::m01_limit) + .def_readwrite("m10_limit", &StampParameters::m10_limit) + .def_readwrite("m11_limit", &StampParameters::m11_limit) + .def_readwrite("m02_limit", &StampParameters::m02_limit) + .def_readwrite("m20_limit", &StampParameters::m20_limit); +} #endif /* Py_PYTHON_H */ diff --git a/src/kbmod/search/filtering.cpp b/src/kbmod/search/filtering.cpp index 7727cfa6..912a82aa 100644 --- a/src/kbmod/search/filtering.cpp +++ b/src/kbmod/search/filtering.cpp @@ -1,19 +1,18 @@ #include "filtering.h" #include - namespace search { #ifdef HAVE_CUDA - /* The filter_kenerls.cu functions. */ - extern "C" void SigmaGFilteredIndicesCU(float *values, int num_values, float sgl0, float sgl1, float sg_coeff, - float width, int *idx_array, int *min_keep_idx, int *max_keep_idx); +/* The filter_kenerls.cu functions. */ +extern "C" void SigmaGFilteredIndicesCU(float *values, int num_values, float sgl0, float sgl1, float sg_coeff, + float width, int *idx_array, int *min_keep_idx, int *max_keep_idx); #endif - /* Return the list of indices from the values array such that those elements - pass the sigmaG filtering defined by percentiles [sgl0, sgl1] with coefficient - sigma_g_coeff and a multiplicative factor of width. */ - std::vector sigmaGFilteredIndices(const std::vector &values, float sgl0, float sgl1, - float sigma_g_coeff, float width) { +/* Return the list of indices from the values array such that those elements + pass the sigmaG filtering defined by percentiles [sgl0, sgl1] with coefficient + sigma_g_coeff and a multiplicative factor of width. */ +std::vector sigmaGFilteredIndices(const std::vector &values, float sgl0, float sgl1, + float sigma_g_coeff, float width) { // Bounds check the percentile values. assert(sgl0 > 0.0); assert(sgl1 < 1.0); @@ -23,7 +22,7 @@ namespace search { float values_arr[num_values]; int idx_array[num_values]; for (int i = 0; i < num_values; ++i) { - values_arr[i] = values[i]; + values_arr[i] = values[i]; } int min_keep_idx = 0; @@ -39,28 +38,28 @@ namespace search { // Copy the result into a vector and return it. std::vector result; for (int i = min_keep_idx; i <= max_keep_idx; ++i) { - result.push_back(idx_array[i]); + result.push_back(idx_array[i]); } return result; - } +} - /* Given a set of psi and phi values, - return a likelihood value */ - double calculateLikelihoodFromPsiPhi(std::vector psi_values, std::vector phi_values) { +/* Given a set of psi and phi values, + return a likelihood value */ +double calculateLikelihoodFromPsiPhi(std::vector psi_values, std::vector phi_values) { assert(psi_values.size() == phi_values.size()); double psi_sum = 0.0; double phi_sum = 0.0; for (int i = 0; i < psi_values.size(); i++) { - psi_sum += psi_values[i]; - phi_sum += phi_values[i]; + psi_sum += psi_values[i]; + phi_sum += phi_values[i]; } if (psi_sum == 0.0 || phi_sum <= 0.0) { - return 0.0; + return 0.0; } return psi_sum / sqrt(phi_sum); - } - +} + } /* namespace search */ diff --git a/src/kbmod/search/filtering.h b/src/kbmod/search/filtering.h index bcb19f4c..96c658c0 100644 --- a/src/kbmod/search/filtering.h +++ b/src/kbmod/search/filtering.h @@ -3,13 +3,12 @@ #include - namespace search { - /* Return the list of indices from the values array such that those elements - pass the sigmaG filtering defined by percentiles [sGL0, sGL1] with coefficient - sigmag_coeff and a multiplicative factor of width. */ - std::vector sigmaGFilteredIndices(const std::vector& values, float sgl0, float sgl1, - float sigma_g_coeff, float width); +/* Return the list of indices from the values array such that those elements + pass the sigmaG filtering defined by percentiles [sGL0, sGL1] with coefficient + sigmag_coeff and a multiplicative factor of width. */ +std::vector sigmaGFilteredIndices(const std::vector& values, float sgl0, float sgl1, + float sigma_g_coeff, float width); } /* namespace search */ diff --git a/src/kbmod/search/image_stack.cpp b/src/kbmod/search/image_stack.cpp index bcc9e07f..241a5881 100644 --- a/src/kbmod/search/image_stack.cpp +++ b/src/kbmod/search/image_stack.cpp @@ -1,153 +1,150 @@ #include "image_stack.h" - namespace py = pybind11; - namespace search { - ImageStack::ImageStack(const std::vector& filenames, const std::vector& psfs) { +ImageStack::ImageStack(const std::vector& filenames, const std::vector& psfs) { verbose = true; images = std::vector(); load_images(filenames, psfs); global_mask = RawImage(get_width(), get_height()); global_mask.set_all_pix(0.0); - } +} - ImageStack::ImageStack(const std::vector& imgs) { +ImageStack::ImageStack(const std::vector& imgs) { verbose = true; images = imgs; global_mask = RawImage(get_width(), get_height()); global_mask.set_all_pix(0.0); - } +} - void ImageStack::load_images(const std::vector& filenames, - const std::vector& psfs) { +void ImageStack::load_images(const std::vector& filenames, const std::vector& psfs) { const int num_files = filenames.size(); if (num_files == 0) { - std::cout << "No files provided" - << "\n"; + std::cout << "No files provided" + << "\n"; } if (psfs.size() != num_files) throw std::runtime_error("Mismatched PSF array in ImageStack creation."); // Load images from file for (int i = 0; i < num_files; ++i) { - images.push_back(LayeredImage(filenames[i], psfs[i])); - if (verbose) std::cout << "." << std::flush; + images.push_back(LayeredImage(filenames[i], psfs[i])); + if (verbose) std::cout << "." << std::flush; } if (verbose) std::cout << "\n"; - } +} - LayeredImage& ImageStack::get_single_image(int index) { +LayeredImage& ImageStack::get_single_image(int index) { if (index < 0 || index > images.size()) throw std::out_of_range("ImageStack index out of bounds."); return images[index]; - } +} - float ImageStack::get_obstime(int index) const { +float ImageStack::get_obstime(int index) const { if (index < 0 || index > images.size()) throw std::out_of_range("ImageStack index out of bounds."); return images[index].get_obstime(); - } +} - float ImageStack::get_zeroed_time(int index) const { +float ImageStack::get_zeroed_time(int index) const { if (index < 0 || index > images.size()) throw std::out_of_range("ImageStack index out of bounds."); return images[index].get_obstime() - images[0].get_obstime(); - } +} - std::vector ImageStack::build_zeroed_times() const { +std::vector ImageStack::build_zeroed_times() const { std::vector zeroed_times = std::vector(); if (images.size() > 0) { - float t0 = images[0].get_obstime(); - for (auto& i : images) { - zeroed_times.push_back(i.get_obstime() - t0); - } + float t0 = images[0].get_obstime(); + for (auto& i : images) { + zeroed_times.push_back(i.get_obstime() - t0); + } } return zeroed_times; - } - - void ImageStack::convolve_psf() { +} + +void ImageStack::convolve_psf() { for (auto& i : images) i.convolve_psf(); - } +} - void ImageStack::save_global_mask(const std::string& path) { global_mask.save_to_file(path); } +void ImageStack::save_global_mask(const std::string& path) { global_mask.save_to_file(path); } - void ImageStack::save_images(const std::string& path) { +void ImageStack::save_images(const std::string& path) { for (auto& i : images) i.save_layers(path); - } +} - const RawImage& ImageStack::get_global_mask() const { return global_mask; } +const RawImage& ImageStack::get_global_mask() const { return global_mask; } - void ImageStack::apply_mask_flags(int flags, const std::vector& exceptions) { +void ImageStack::apply_mask_flags(int flags, const std::vector& exceptions) { for (auto& i : images) { - i.apply_mask_flags(flags, exceptions); + i.apply_mask_flags(flags, exceptions); } - } +} - void ImageStack::apply_global_mask(int flags, int threshold) { +void ImageStack::apply_global_mask(int flags, int threshold) { create_global_mask(flags, threshold); for (auto& i : images) { - i.apply_global_mask(global_mask); + i.apply_global_mask(global_mask); } - } +} - void ImageStack::apply_mask_threshold(float thresh) { +void ImageStack::apply_mask_threshold(float thresh) { for (auto& i : images) i.apply_mask_threshold(thresh); - } +} - void ImageStack::grow_mask(int steps) { +void ImageStack::grow_mask(int steps) { for (auto& i : images) i.grow_mask(steps); - } +} - void ImageStack::create_global_mask(int flags, int threshold) { +void ImageStack::create_global_mask(int flags, int threshold) { int npixels = get_npixels(); // For each pixel count the number of images where it is masked. std::vector counts(npixels, 0); for (unsigned int img = 0; img < images.size(); ++img) { - float* imgMask = images[img].get_mask().data(); - // Count the number of times a pixel has any of the flags - for (unsigned int pixel = 0; pixel < npixels; ++pixel) { - if ((flags & static_cast(imgMask[pixel])) != 0) counts[pixel]++; - } + float* imgMask = images[img].get_mask().data(); + // Count the number of times a pixel has any of the flags + for (unsigned int pixel = 0; pixel < npixels; ++pixel) { + if ((flags & static_cast(imgMask[pixel])) != 0) counts[pixel]++; + } } // Set all pixels below threshold to 0 and all above to 1 float* global_m = global_mask.data(); for (unsigned int p = 0; p < npixels; ++p) { - global_m[p] = counts[p] < threshold ? 0.0 : 1.0; + global_m[p] = counts[p] < threshold ? 0.0 : 1.0; } - } +} #ifdef Py_PYTHON_H - static void image_stack_bindings(py::module &m) { +static void image_stack_bindings(py::module& m) { using is = search::ImageStack; using li = search::LayeredImage; using pf = search::PSF; py::class_(m, "ImageStack", pydocs::DOC_ImageStack) - .def(py::init, std::vector>()) - .def(py::init>()) - .def("get_images", &is::get_images, pydocs::DOC_ImageStack_get_images) - .def("get_single_image", &is::get_single_image, - py::return_value_policy::reference_internal, - pydocs::DOC_ImageStack_get_single_image) - .def("get_obstime", &is::get_obstime, pydocs::DOC_ImageStack_get_obstime) - .def("get_zeroed_time", &is::get_zeroed_time, pydocs::DOC_ImageStack_get_zeroed_time) - .def("build_zeroed_times", &is::build_zeroed_times, pydocs::DOC_ImageStack_build_zeroed_times) - .def("img_count", &is::img_count, pydocs::DOC_ImageStack_img_count) - .def("apply_mask_flags", &is::apply_mask_flags, pydocs::DOC_ImageStack_apply_mask_flags) - .def("apply_mask_threshold", &is::apply_mask_threshold, pydocs::DOC_ImageStack_apply_mask_threshold) - .def("apply_global_mask", &is::apply_global_mask, pydocs::DOC_ImageStack_apply_global_mask) - .def("grow_mask", &is::grow_mask, pydocs::DOC_ImageStack_grow_mask) - .def("save_global_mask", &is::save_global_mask, pydocs::DOC_ImageStack_save_global_mask) - .def("save_images", &is::save_images, pydocs::DOC_ImageStack_save_images) - .def("get_global_mask", &is::get_global_mask, pydocs::DOC_ImageStack_get_global_mask) - .def("convolve_psf", &is::convolve_psf, pydocs::DOC_ImageStack_convolve_psf) - .def("get_width", &is::get_width, pydocs::DOC_ImageStack_get_width) - .def("get_height", &is::get_height, pydocs::DOC_ImageStack_get_height) - .def("get_npixels", &is::get_npixels, pydocs::DOC_ImageStack_get_npixels); - } + .def(py::init, std::vector>()) + .def(py::init>()) + .def("get_images", &is::get_images, pydocs::DOC_ImageStack_get_images) + .def("get_single_image", &is::get_single_image, py::return_value_policy::reference_internal, + pydocs::DOC_ImageStack_get_single_image) + .def("get_obstime", &is::get_obstime, pydocs::DOC_ImageStack_get_obstime) + .def("get_zeroed_time", &is::get_zeroed_time, pydocs::DOC_ImageStack_get_zeroed_time) + .def("build_zeroed_times", &is::build_zeroed_times, pydocs::DOC_ImageStack_build_zeroed_times) + .def("img_count", &is::img_count, pydocs::DOC_ImageStack_img_count) + .def("apply_mask_flags", &is::apply_mask_flags, pydocs::DOC_ImageStack_apply_mask_flags) + .def("apply_mask_threshold", &is::apply_mask_threshold, + pydocs::DOC_ImageStack_apply_mask_threshold) + .def("apply_global_mask", &is::apply_global_mask, pydocs::DOC_ImageStack_apply_global_mask) + .def("grow_mask", &is::grow_mask, pydocs::DOC_ImageStack_grow_mask) + .def("save_global_mask", &is::save_global_mask, pydocs::DOC_ImageStack_save_global_mask) + .def("save_images", &is::save_images, pydocs::DOC_ImageStack_save_images) + .def("get_global_mask", &is::get_global_mask, pydocs::DOC_ImageStack_get_global_mask) + .def("convolve_psf", &is::convolve_psf, pydocs::DOC_ImageStack_convolve_psf) + .def("get_width", &is::get_width, pydocs::DOC_ImageStack_get_width) + .def("get_height", &is::get_height, pydocs::DOC_ImageStack_get_height) + .def("get_npixels", &is::get_npixels, pydocs::DOC_ImageStack_get_npixels); +} #endif /* Py_PYTHON_H */ diff --git a/src/kbmod/search/image_stack.h b/src/kbmod/search/image_stack.h index f9c5acac..f1758d50 100644 --- a/src/kbmod/search/image_stack.h +++ b/src/kbmod/search/image_stack.h @@ -10,10 +10,9 @@ #include "layered_image.h" #include "pydocs/image_stack_docs.h" - namespace search { - class ImageStack { - public: +class ImageStack { +public: ImageStack(const std::vector& filenames, const std::vector& psfs); ImageStack(const std::vector& imgs); @@ -45,13 +44,13 @@ namespace search { virtual ~ImageStack(){}; - private: +private: void load_images(const std::vector& filenames, const std::vector& psfs); void create_global_mask(int flags, int threshold); std::vector images; RawImage global_mask; bool verbose; - }; +}; } /* namespace search */ diff --git a/src/kbmod/search/layered_image.cpp b/src/kbmod/search/layered_image.cpp index f536d978..05e210ce 100644 --- a/src/kbmod/search/layered_image.cpp +++ b/src/kbmod/search/layered_image.cpp @@ -1,8 +1,7 @@ #include "layered_image.h" - namespace search { - LayeredImage::LayeredImage(std::string path, const PSF& psf) : psf(psf) { +LayeredImage::LayeredImage(std::string path, const PSF& psf) : psf(psf) { int f_begin = path.find_last_of("/"); int f_end = path.find_last_of(".fits") - 4; filename = path.substr(f_begin, f_end - f_begin); @@ -19,36 +18,35 @@ namespace search { variance.load_from_file(path, 3); if (width != variance.get_width() or height != variance.get_height()) - throw std::runtime_error("Science and Variance layers are not the same size."); + throw std::runtime_error("Science and Variance layers are not the same size."); if (width != mask.get_width() or height != mask.get_height()) - throw std::runtime_error("Science and Mask layers are not the same size."); - } + throw std::runtime_error("Science and Mask layers are not the same size."); +} - LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, - const PSF& psf) - : psf(psf) { +LayeredImage::LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf) + : psf(psf) { // Get the dimensions of the science layer and check for consistency with // the other two layers. width = sci.get_width(); height = sci.get_height(); if (width != var.get_width() or height != var.get_height()) - throw std::runtime_error("Science and Variance layers are not the same size."); + throw std::runtime_error("Science and Variance layers are not the same size."); if (width != msk.get_width() or height != msk.get_height()) - throw std::runtime_error("Science and Mask layers are not the same size."); + throw std::runtime_error("Science and Mask layers are not the same size."); // Copy the image layers. science = sci; mask = msk; variance = var; - } +} - LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, double time, - const PSF& psf) - : LayeredImage(name, w, h, noise_stdev, pixel_variance, time, psf, -1) {} +LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, + double time, const PSF& psf) + : LayeredImage(name, w, h, noise_stdev, pixel_variance, time, psf, -1) {} - LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, double time, - const PSF& psf, int seed) - : psf(psf) { +LayeredImage::LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, + double time, const PSF& psf, int seed) + : psf(psf) { filename = name; width = w; height = h; @@ -57,7 +55,7 @@ namespace search { std::random_device r; std::default_random_engine generator(r()); if (seed >= 0) { - generator.seed(seed); + generator.seed(seed); } std::normal_distribution distrib(0.0, noise_stdev); for (float& p : raw_sci) p = distrib(generator); @@ -67,67 +65,63 @@ namespace search { mask = RawImage(w, h, std::vector(w * h, 0.0)); variance = RawImage(w, h, std::vector(w * h, pixel_variance)); - } +} - void LayeredImage::set_psf(const PSF& new_psf) { - psf = new_psf; - } +void LayeredImage::set_psf(const PSF& new_psf) { psf = new_psf; } - void LayeredImage::grow_mask(int steps) { +void LayeredImage::grow_mask(int steps) { science.grow_mask(steps); variance.grow_mask(steps); - } +} - void LayeredImage::convolve_given_psf(const PSF& given_psf) { +void LayeredImage::convolve_given_psf(const PSF& given_psf) { science.convolve(given_psf); // Square the PSF use that on the variance image. PSF psfsq = PSF(given_psf); // Copy psfsq.square_psf(); variance.convolve(psfsq); - } +} - void LayeredImage::convolve_psf() { - convolve_given_psf(psf); - } +void LayeredImage::convolve_psf() { convolve_given_psf(psf); } - void LayeredImage::apply_mask_flags(int flags, const std::vector& exceptions) { +void LayeredImage::apply_mask_flags(int flags, const std::vector& exceptions) { science.apply_mask(flags, exceptions, mask); variance.apply_mask(flags, exceptions, mask); - } +} - /* Mask all pixels that are not 0 in global mask */ - void LayeredImage::apply_global_mask(const RawImage& global_mask) { +/* Mask all pixels that are not 0 in global mask */ +void LayeredImage::apply_global_mask(const RawImage& global_mask) { science.apply_mask(0xFFFFFF, {}, global_mask); variance.apply_mask(0xFFFFFF, {}, global_mask); - } +} - void LayeredImage::apply_mask_threshold(float thresh) { +void LayeredImage::apply_mask_threshold(float thresh) { const int num_pixels = get_npixels(); float* sci_pixels = science.data(); float* var_pix = variance.data(); for (int i = 0; i < num_pixels; ++i) { - if (sci_pixels[i] > thresh) { - sci_pixels[i] = NO_DATA; - var_pix[i] = NO_DATA; - } + if (sci_pixels[i] > thresh) { + sci_pixels[i] = NO_DATA; + var_pix[i] = NO_DATA; + } } - } +} - void LayeredImage::subtract_template(const RawImage& sub_template) { +void LayeredImage::subtract_template(const RawImage& sub_template) { assert(get_height() == sub_template.get_height() && get_width() == sub_template.get_width()); const int num_pixels = get_npixels(); float* sci_pixels = science.data(); const std::vector& tem_pixels = sub_template.get_pixels(); for (unsigned i = 0; i < num_pixels; ++i) { - if ((sci_pixels[i] != NO_DATA) && (tem_pixels[i] != NO_DATA)) { - sci_pixels[i] -= tem_pixels[i]; - } + if ((sci_pixels[i] != NO_DATA) && (tem_pixels[i] != NO_DATA)) { + sci_pixels[i] -= tem_pixels[i]; + } } - } +} - void LayeredImage::save_layers(const std::string& path) { +void LayeredImage::save_layers(const std::string& path) { fitsfile* fptr; int status = 0; long naxes[2] = {0, 0}; @@ -138,12 +132,12 @@ namespace search { // If we are unable to create the file, check if it already exists // and, if so, delete it and retry the create. if (status == 105) { - status = 0; - fits_open_file(&fptr, (path + filename + ".fits").c_str(), READWRITE, &status); - if (status == 0) { - fits_delete_file(fptr, &status); - fits_create_file(&fptr, (path + filename + ".fits").c_str(), &status); - } + status = 0; + fits_open_file(&fptr, (path + filename + ".fits").c_str(), READWRITE, &status); + if (status == 0) { + fits_delete_file(fptr, &status); + fits_create_file(&fptr, (path + filename + ".fits").c_str(), &status); + } } fits_create_img(fptr, SHORT_IMG, 0, naxes, &status); @@ -154,29 +148,29 @@ namespace search { science.append_layer_to_file(path + filename + ".fits"); mask.append_layer_to_file(path + filename + ".fits"); variance.append_layer_to_file(path + filename + ".fits"); - } +} - void LayeredImage::set_science(RawImage& im) { +void LayeredImage::set_science(RawImage& im) { check_dims(im); science = im; - } +} - void LayeredImage::set_mask(RawImage& im) { +void LayeredImage::set_mask(RawImage& im) { check_dims(im); mask = im; - } +} - void LayeredImage::set_variance(RawImage& im) { +void LayeredImage::set_variance(RawImage& im) { check_dims(im); variance = im; - } +} - void LayeredImage::check_dims(RawImage& im) { +void LayeredImage::check_dims(RawImage& im) { if (im.get_width() != get_width()) throw std::runtime_error("Image width does not match"); if (im.get_height() != get_height()) throw std::runtime_error("Image height does not match"); - } +} - RawImage LayeredImage::generate_psi_image() { +RawImage LayeredImage::generate_psi_image() { RawImage result(width, height); float* result_arr = result.data(); float* sci_array = science.data(); @@ -185,21 +179,21 @@ namespace search { // Set each of the result pixels. const int num_pixels = get_npixels(); for (int p = 0; p < num_pixels; ++p) { - float var_pix = var_array[p]; - if (var_pix != NO_DATA) { - result_arr[p] = sci_array[p] / var_pix; - } else { - result_arr[p] = NO_DATA; - } + float var_pix = var_array[p]; + if (var_pix != NO_DATA) { + result_arr[p] = sci_array[p] / var_pix; + } else { + result_arr[p] = NO_DATA; + } } // Convolve with the PSF. result.convolve(psf); return result; - } +} - RawImage LayeredImage::generate_phi_image() { +RawImage LayeredImage::generate_phi_image() { RawImage result(width, height); float* result_arr = result.data(); float* var_array = variance.data(); @@ -207,12 +201,12 @@ namespace search { // Set each of the result pixels. const int num_pixels = get_npixels(); for (int p = 0; p < num_pixels; ++p) { - float var_pix = var_array[p]; - if (var_pix != NO_DATA) { - result_arr[p] = 1.0 / var_pix; - } else { - result_arr[p] = NO_DATA; - } + float var_pix = var_array[p]; + if (var_pix != NO_DATA) { + result_arr[p] = 1.0 / var_pix; + } else { + result_arr[p] = NO_DATA; + } } // Convolve with the PSF squared. @@ -221,53 +215,49 @@ namespace search { result.convolve(psfsq); return result; - } +} #ifdef Py_PYTHON_H - static void layered_image_bindings(py::module &m) { +static void layered_image_bindings(py::module& m) { using li = search::LayeredImage; using ri = search::RawImage; using pf = search::PSF; py::class_
  • (m, "LayeredImage", pydocs::DOC_LayeredImage) - .def(py::init()) - .def(py::init()) - .def(py::init()) - .def(py::init()) - .def("set_psf", &li::set_psf, pydocs::DOC_LayeredImage_set_psf) - .def("get_psf", &li::get_psf, - py::return_value_policy::reference_internal, - pydocs::DOC_LayeredImage_get_psf) - .def("apply_mask_flags", &li::apply_mask_flags, pydocs::DOC_LayeredImage_apply_mask_flags) - .def("apply_mask_threshold", &li::apply_mask_threshold, pydocs::DOC_LayeredImage_apply_mask_threshold) - .def("sub_template", &li::subtract_template, pydocs::DOC_LayeredImage_sub_template) - .def("save_layers", &li::save_layers, pydocs::DOC_LayeredImage_save_layers) - .def("get_science", &li::get_science, - py::return_value_policy::reference_internal, - pydocs::DOC_LayeredImage_get_science) - .def("get_mask", &li::get_mask, - py::return_value_policy::reference_internal, - pydocs::DOC_LayeredImage_get_mask) - .def("get_variance", &li::get_variance, - py::return_value_policy::reference_internal, - pydocs::DOC_LayeredImage_get_variance) - .def("set_science", &li::set_science, pydocs::DOC_LayeredImage_set_science) - .def("set_mask", &li::set_mask, pydocs::DOC_LayeredImage_set_mask) - .def("set_variance", &li::set_variance, pydocs::DOC_LayeredImage_set_variance) - .def("convolve_psf", &li::convolve_psf, pydocs::DOC_LayeredImage_convolve_psf) - .def("convolve_given_psf", &li::convolve_given_psf, pydocs::DOC_LayeredImage_convolve_given_psf) - .def("grow_mask", &li::grow_mask, pydocs::DOC_LayeredImage_grow_mask) - .def("get_name", &li::get_name, pydocs::DOC_LayeredImage_get_name) - .def("get_width", &li::get_width, pydocs::DOC_LayeredImage_get_width) - .def("get_height", &li::get_height, pydocs::DOC_LayeredImage_get_height) - .def("get_npixels", &li::get_npixels, pydocs::DOC_LayeredImage_get_npixels) - .def("get_obstime", &li::get_obstime, pydocs::DOC_LayeredImage_get_obstime) - .def("set_obstime", &li::set_obstime, pydocs::DOC_LayeredImage_set_obstime) - .def("generate_psi_image", &li::generate_psi_image, pydocs::DOC_LayeredImage_generate_psi_image) - .def("generate_phi_image", &li::generate_phi_image, pydocs::DOC_LayeredImage_generate_phi_image); - } + .def(py::init()) + .def(py::init()) + .def(py::init()) + .def(py::init()) + .def("set_psf", &li::set_psf, pydocs::DOC_LayeredImage_set_psf) + .def("get_psf", &li::get_psf, py::return_value_policy::reference_internal, + pydocs::DOC_LayeredImage_get_psf) + .def("apply_mask_flags", &li::apply_mask_flags, pydocs::DOC_LayeredImage_apply_mask_flags) + .def("apply_mask_threshold", &li::apply_mask_threshold, + pydocs::DOC_LayeredImage_apply_mask_threshold) + .def("sub_template", &li::subtract_template, pydocs::DOC_LayeredImage_sub_template) + .def("save_layers", &li::save_layers, pydocs::DOC_LayeredImage_save_layers) + .def("get_science", &li::get_science, py::return_value_policy::reference_internal, + pydocs::DOC_LayeredImage_get_science) + .def("get_mask", &li::get_mask, py::return_value_policy::reference_internal, + pydocs::DOC_LayeredImage_get_mask) + .def("get_variance", &li::get_variance, py::return_value_policy::reference_internal, + pydocs::DOC_LayeredImage_get_variance) + .def("set_science", &li::set_science, pydocs::DOC_LayeredImage_set_science) + .def("set_mask", &li::set_mask, pydocs::DOC_LayeredImage_set_mask) + .def("set_variance", &li::set_variance, pydocs::DOC_LayeredImage_set_variance) + .def("convolve_psf", &li::convolve_psf, pydocs::DOC_LayeredImage_convolve_psf) + .def("convolve_given_psf", &li::convolve_given_psf, pydocs::DOC_LayeredImage_convolve_given_psf) + .def("grow_mask", &li::grow_mask, pydocs::DOC_LayeredImage_grow_mask) + .def("get_name", &li::get_name, pydocs::DOC_LayeredImage_get_name) + .def("get_width", &li::get_width, pydocs::DOC_LayeredImage_get_width) + .def("get_height", &li::get_height, pydocs::DOC_LayeredImage_get_height) + .def("get_npixels", &li::get_npixels, pydocs::DOC_LayeredImage_get_npixels) + .def("get_obstime", &li::get_obstime, pydocs::DOC_LayeredImage_get_obstime) + .def("set_obstime", &li::set_obstime, pydocs::DOC_LayeredImage_set_obstime) + .def("generate_psi_image", &li::generate_psi_image, pydocs::DOC_LayeredImage_generate_psi_image) + .def("generate_phi_image", &li::generate_phi_image, pydocs::DOC_LayeredImage_generate_phi_image); +} #endif /* Py_PYTHON_H */ } /* namespace search */ - diff --git a/src/kbmod/search/layered_image.h b/src/kbmod/search/layered_image.h index b4b3abb4..8a0d061f 100644 --- a/src/kbmod/search/layered_image.h +++ b/src/kbmod/search/layered_image.h @@ -12,17 +12,15 @@ #include "common.h" #include "pydocs/layered_image_docs.h" - namespace search { - class LayeredImage { - public: +class LayeredImage { +public: explicit LayeredImage(std::string path, const PSF& psf); - explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, - const PSF& psf); - explicit LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, double time, - const PSF& psf); - explicit LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, double time, - const PSF& psf, int seed); + explicit LayeredImage(const RawImage& sci, const RawImage& var, const RawImage& msk, const PSF& psf); + explicit LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, + double time, const PSF& psf); + explicit LayeredImage(std::string name, int w, int h, float noise_stdev, float pixel_variance, + double time, const PSF& psf, int seed); // Set an image specific point spread function. void set_psf(const PSF& psf); @@ -68,7 +66,7 @@ namespace search { RawImage generate_psi_image(); RawImage generate_phi_image(); - private: +private: void check_dims(RawImage& im); std::string filename; @@ -79,7 +77,7 @@ namespace search { RawImage science; RawImage mask; RawImage variance; - }; +}; } /* namespace search */ diff --git a/src/kbmod/search/psf.cpp b/src/kbmod/search/psf.cpp index c75107e4..696d66dd 100644 --- a/src/kbmod/search/psf.cpp +++ b/src/kbmod/search/psf.cpp @@ -1,19 +1,16 @@ #include "psf.h" - namespace py = pybind11; - - namespace search { - PSF::PSF() : kernel(1, 1.0) { +PSF::PSF() : kernel(1, 1.0) { dim = 1; radius = 0; width = 1e-12; - sum = 1.0; - } + sum = 1.0; +} - PSF::PSF(float stdev) { +PSF::PSF(float stdev) { width = stdev; float simple_gauss[MAX_KERNEL_RADIUS]; double psf_coverage = 0.0; @@ -22,15 +19,15 @@ namespace search { // Create 1D gaussian array while (psf_coverage < 0.98 && i < MAX_KERNEL_RADIUS) { - float current_bin = - 0.5 * (std::erf((float(i) + 0.5) / norm_factor) - std::erf((float(i) - 0.5) / norm_factor)); - simple_gauss[i] = current_bin; - if (i == 0) { - psf_coverage += current_bin; - } else { - psf_coverage += 2.0 * current_bin; - } - i++; + float current_bin = + 0.5 * (std::erf((float(i) + 0.5) / norm_factor) - std::erf((float(i) - 0.5) / norm_factor)); + simple_gauss[i] = current_bin; + if (i == 0) { + psf_coverage += current_bin; + } else { + psf_coverage += 2.0 * current_bin; + } + i++; } radius = i - 1; // This value is good for @@ -39,137 +36,138 @@ namespace search { // Create 2D gaussain by multiplying with itself kernel = std::vector(); for (int ii = 0; ii < dim; ++ii) { - for (int jj = 0; jj < dim; ++jj) { - float current = simple_gauss[abs(radius - ii)] * simple_gauss[abs(radius - jj)]; - kernel.push_back(current); - } + for (int jj = 0; jj < dim; ++jj) { + float current = simple_gauss[abs(radius - ii)] * simple_gauss[abs(radius - jj)]; + kernel.push_back(current); + } } calc_sum(); - } +} - // Copy constructor. - PSF::PSF(const PSF& other) { +// Copy constructor. +PSF::PSF(const PSF& other) { kernel = other.kernel; dim = other.dim; radius = other.radius; width = other.width; sum = other.sum; - } +} - // Copy assignment. - PSF& PSF::operator=(const PSF& other) { +// Copy assignment. +PSF& PSF::operator=(const PSF& other) { kernel = other.kernel; dim = other.dim; radius = other.radius; width = other.width; sum = other.sum; return *this; - } - - // Move constructor. - PSF::PSF(PSF&& other) - : kernel(std::move(other.kernel)), - dim(other.dim), - radius(other.radius), - width(other.width), - sum(other.sum) {} - - // Move assignment. - PSF& PSF::operator=(PSF&& other) { +} + +// Move constructor. +PSF::PSF(PSF&& other) + : kernel(std::move(other.kernel)), + dim(other.dim), + radius(other.radius), + width(other.width), + sum(other.sum) {} + +// Move assignment. +PSF& PSF::operator=(PSF&& other) { if (this != &other) { - kernel = std::move(other.kernel); - dim = other.dim; - radius = other.radius; - width = other.width; - sum = other.sum; + kernel = std::move(other.kernel); + dim = other.dim; + radius = other.radius; + width = other.width; + sum = other.sum; } return *this; - } +} - void PSF::calc_sum() { +void PSF::calc_sum() { sum = 0.0; for (auto& i : kernel) sum += i; - } +} - void PSF::square_psf() { +void PSF::square_psf() { for (float& i : kernel) { - i = i * i; + i = i * i; } calc_sum(); - } +} - std::string PSF::print() { +std::string PSF::print() { std::stringstream ss; ss.setf(std::ios::fixed, std::ios::floatfield); ss.precision(3); for (int row = 0; row < dim; ++row) { - ss << "| "; - for (int col = 0; col < dim; ++col) { - ss << kernel[row * dim + col] << " | "; - } - ss << "\n "; - for (int space = 0; space < dim * 8 - 1; ++space) ss << "-"; - ss << "\n"; + ss << "| "; + for (int col = 0; col < dim; ++col) { + ss << kernel[row * dim + col] << " | "; + } + ss << "\n "; + for (int space = 0; space < dim * 8 - 1; ++space) ss << "-"; + ss << "\n"; } ss << 100.0 * sum << "% of PSF contained within kernel\n"; return ss.str(); - } +} #ifdef Py_PYTHON_H - PSF::PSF(pybind11::array_t arr) { set_array(arr); } +PSF::PSF(pybind11::array_t arr) { set_array(arr); } - void PSF::set_array(pybind11::array_t arr) { +void PSF::set_array(pybind11::array_t arr) { pybind11::buffer_info info = arr.request(); if (info.ndim != 2) - throw std::runtime_error( - "Array must have 2 dimensions. (It " - " must also be a square with odd dimensions)"); + throw std::runtime_error( + "Array must have 2 dimensions. (It " + " must also be a square with odd dimensions)"); if (info.shape[0] != info.shape[1]) - throw std::runtime_error( - "Array must be square (x-dimension == y-dimension)." - "It also must have odd dimensions."); + throw std::runtime_error( + "Array must be square (x-dimension == y-dimension)." + "It also must have odd dimensions."); float* pix = static_cast(info.ptr); dim = info.shape[0]; if (dim % 2 == 0) - throw std::runtime_error( - "Array dimension must be odd. The " - "middle of an even numbered array is ambiguous."); + throw std::runtime_error( + "Array dimension must be odd. The " + "middle of an even numbered array is ambiguous."); radius = dim / 2; // Rounds down sum = 0.0; kernel = std::vector(pix, pix + dim * dim); calc_sum(); width = 0.0; - } +} - static void psf_bindings(py::module &m) { +static void psf_bindings(py::module& m) { using psf = search::PSF; py::class_(m, "PSF", py::buffer_protocol(), pydocs::DOC_PSF) - .def_buffer([](psf &m) -> py::buffer_info { - return py::buffer_info(m.data(), // void *ptr; - sizeof(float), // py::ssize_t itemsize; - py::format_descriptor::format(), // std::string format; - 2, // py::ssize_t ndim; - {m.get_dim(), m.get_dim()}, // std::vector shape; - {sizeof(float) * m.get_dim(), sizeof(float)}); // std::vector strides; - }) - .def(py::init<>()) - .def(py::init()) - .def(py::init>()) - .def(py::init()) - .def("set_array", &psf::set_array, pydocs::DOC_PSF_set_array) - .def("get_std", &psf::get_std, pydocs::DOC_PSF_get_std) - .def("get_sum", &psf::get_sum, pydocs::DOC_PSF_get_sum) - .def("get_dim", &psf::get_dim, pydocs::DOC_PSF_get_dim) - .def("get_radius", &psf::get_radius, pydocs::DOC_PSF_get_radius) - .def("get_size", &psf::get_size, pydocs::DOC_PSF_get_size) - .def("get_kernel", &psf::get_kernel, pydocs::DOC_PSF_get_kernel) - .def("get_value", &psf::get_value, pydocs::DOC_PSF_get_value) - .def("square_psf", &psf::square_psf, pydocs::DOC_PSF_square_psf) - .def("print", &psf::print, pydocs::DOC_PSF_print); - } + .def_buffer([](psf& m) -> py::buffer_info { + return py::buffer_info( + m.data(), // void *ptr; + sizeof(float), // py::ssize_t itemsize; + py::format_descriptor::format(), // std::string format; + 2, // py::ssize_t ndim; + {m.get_dim(), m.get_dim()}, // std::vector shape; + {sizeof(float) * m.get_dim(), sizeof(float)}); // std::vector strides; + }) + .def(py::init<>()) + .def(py::init()) + .def(py::init>()) + .def(py::init()) + .def("set_array", &psf::set_array, pydocs::DOC_PSF_set_array) + .def("get_std", &psf::get_std, pydocs::DOC_PSF_get_std) + .def("get_sum", &psf::get_sum, pydocs::DOC_PSF_get_sum) + .def("get_dim", &psf::get_dim, pydocs::DOC_PSF_get_dim) + .def("get_radius", &psf::get_radius, pydocs::DOC_PSF_get_radius) + .def("get_size", &psf::get_size, pydocs::DOC_PSF_get_size) + .def("get_kernel", &psf::get_kernel, pydocs::DOC_PSF_get_kernel) + .def("get_value", &psf::get_value, pydocs::DOC_PSF_get_value) + .def("square_psf", &psf::square_psf, pydocs::DOC_PSF_square_psf) + .def("print", &psf::print, pydocs::DOC_PSF_print); +} #endif } /* namespace search */ diff --git a/src/kbmod/search/psf.h b/src/kbmod/search/psf.h index 0a38c873..37b7a4d2 100644 --- a/src/kbmod/search/psf.h +++ b/src/kbmod/search/psf.h @@ -9,10 +9,9 @@ #include "common.h" #include "pydocs/psf_docs.h" - namespace search { - class PSF { - public: +class PSF { +public: PSF(); // Create a no-op PSF. PSF(float stdev); PSF(const PSF& other); // Copy constructor @@ -43,13 +42,13 @@ namespace search { void square_psf(); std::string print(); - private: +private: std::vector kernel; float width; float sum; int dim; int radius; - }; +}; } /* namespace search */ diff --git a/src/kbmod/search/raw_image.cpp b/src/kbmod/search/raw_image.cpp index e7f4ed70..a1db1b11 100644 --- a/src/kbmod/search/raw_image.cpp +++ b/src/kbmod/search/raw_image.cpp @@ -1,84 +1,82 @@ #include "raw_image.h" - namespace py = pybind11; - namespace search { #ifdef HAVE_CUDA - // Performs convolution between an image represented as an array of floats - // and a PSF on a GPU device. - extern "C" void deviceConvolve(float* source_img, float* result_img, int width, int height, float* psf_kernel, - int psf_size, int psf_dim, int psf_radius, float psf_sum); +// Performs convolution between an image represented as an array of floats +// and a PSF on a GPU device. +extern "C" void deviceConvolve(float* source_img, float* result_img, int width, int height, float* psf_kernel, + int psf_size, int psf_dim, int psf_radius, float psf_sum); #endif - RawImage::RawImage() : width(0), height(0), obstime(-1.0) { pixels = std::vector(); } +RawImage::RawImage() : width(0), height(0), obstime(-1.0) { pixels = std::vector(); } - // Copy constructor - RawImage::RawImage(const RawImage& old) { +// Copy constructor +RawImage::RawImage(const RawImage& old) { width = old.get_width(); height = old.get_height(); pixels = old.get_pixels(); obstime = old.get_obstime(); - } +} - // Copy assignment - RawImage& RawImage::operator=(const RawImage& source) { +// Copy assignment +RawImage& RawImage::operator=(const RawImage& source) { width = source.width; height = source.height; pixels = source.pixels; obstime = source.obstime; return *this; - } +} - // Move constructor - RawImage::RawImage(RawImage&& source) - : width(source.width), - height(source.height), - obstime(source.obstime), - pixels(std::move(source.pixels)) {} +// Move constructor +RawImage::RawImage(RawImage&& source) + : width(source.width), + height(source.height), + obstime(source.obstime), + pixels(std::move(source.pixels)) {} - // Move assignment - RawImage& RawImage::operator=(RawImage&& source) { +// Move assignment +RawImage& RawImage::operator=(RawImage&& source) { if (this != &source) { - width = source.width; - height = source.height; - pixels = std::move(source.pixels); - obstime = source.obstime; + width = source.width; + height = source.height; + pixels = std::move(source.pixels); + obstime = source.obstime; } return *this; - } +} - RawImage::RawImage(unsigned w, unsigned h) : height(h), width(w), obstime(-1.0), pixels(w * h) {} +RawImage::RawImage(unsigned w, unsigned h) : height(h), width(w), obstime(-1.0), pixels(w * h) {} - RawImage::RawImage(unsigned w, unsigned h, const std::vector& pix) - : width(w), height(h), obstime(-1.0), pixels(pix) { +RawImage::RawImage(unsigned w, unsigned h, const std::vector& pix) + : width(w), height(h), obstime(-1.0), pixels(pix) { assert(w * h == pix.size()); - } +} - bool RawImage::approx_equal(const RawImage& img_b, float atol) const { +bool RawImage::approx_equal(const RawImage& img_b, float atol) const { if ((width != img_b.width) || (height != img_b.height)) return false; for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - float p1 = get_pixel(x, y); - float p2 = img_b.get_pixel(x, y); + for (int x = 0; x < width; ++x) { + float p1 = get_pixel(x, y); + float p2 = img_b.get_pixel(x, y); - // NO_DATA values must match exactly. - if ((p1 == NO_DATA) && (p2 != NO_DATA)) return false; - if ((p1 != NO_DATA) && (p2 == NO_DATA)) return false; + // NO_DATA values must match exactly. + if ((p1 == NO_DATA) && (p2 != NO_DATA)) return false; + if ((p1 != NO_DATA) && (p2 == NO_DATA)) return false; - // Other values match within an absolute tolerance. - if (fabs(p1 - p2) > atol) return false; - } + // Other values match within an absolute tolerance. + if (fabs(p1 - p2) > atol) return false; + } } return true; - } +} - void RawImage::load_time_from_file(fitsfile* fptr) { +void RawImage::load_time_from_file(fitsfile* fptr) { int mjd_status = 0; obstime = -1.0; - + // Read image observation time, trying the MJD field first and DATE-AVG second. // Ignore error if does not exist. if (fits_read_key(fptr, TDOUBLE, "MJD", &obstime, NULL, &mjd_status)) { @@ -86,10 +84,10 @@ namespace search { obstime = -1.0; } } - } - - // Load the image data from a specific layer of a FITS file. - void RawImage::load_from_file(const std::string& file_path, int layer_num) { +} + +// Load the image data from a specific layer of a FITS file. +void RawImage::load_from_file(const std::string& file_path, int layer_num) { // Open the file's header and read in the obstime and the dimensions. fitsfile* fptr; int status = 0; @@ -100,34 +98,34 @@ namespace search { // Open the correct layer to extract the RawImage. std::string layerPath = file_path + "[" + std::to_string(layer_num) + "]"; if (fits_open_file(&fptr, layerPath.c_str(), READONLY, &status)) { - fits_report_error(stderr, status); - throw std::runtime_error("Could not open FITS file to read RawImage"); + fits_report_error(stderr, status); + throw std::runtime_error("Could not open FITS file to read RawImage"); } // Read image dimensions. long dimensions[2]; if (fits_read_keys_lng(fptr, "NAXIS", 1, 2, dimensions, &file_not_found, &status)) - fits_report_error(stderr, status); + fits_report_error(stderr, status); width = dimensions[0]; height = dimensions[1]; // Read in the image. pixels = std::vector(width * height); if (fits_read_img(fptr, TFLOAT, 1, get_npixels(), &nullval, pixels.data(), &anynull, &status)) - fits_report_error(stderr, status); + fits_report_error(stderr, status); load_time_from_file(fptr); if (fits_close_file(fptr, &status)) fits_report_error(stderr, status); - + // If we are reading from a sublayer and did not find a time, try the overall header. if (obstime < 0.0) { - if (fits_open_file(&fptr, file_path.c_str(), READONLY, &status)) - throw std::runtime_error("Could not open FITS file to read RawImage"); - load_time_from_file(fptr); + if (fits_open_file(&fptr, file_path.c_str(), READONLY, &status)) + throw std::runtime_error("Could not open FITS file to read RawImage"); + load_time_from_file(fptr); } - } +} - void RawImage::save_to_file(const std::string& filename) { +void RawImage::save_to_file(const std::string& filename) { fitsfile* fptr; int status = 0; long naxes[2] = {0, 0}; @@ -137,12 +135,12 @@ namespace search { // If we are unable to create the file, check if it already exists // and, if so, delete it and retry the create. if (status == 105) { - status = 0; - fits_open_file(&fptr, filename.c_str(), READWRITE, &status); - if (status == 0) { - fits_delete_file(fptr, &status); - fits_create_file(&fptr, filename.c_str(), &status); - } + status = 0; + fits_open_file(&fptr, filename.c_str(), READWRITE, &status); + if (status == 0) { + fits_delete_file(fptr, &status); + fits_create_file(&fptr, filename.c_str(), &status); + } } // Create the primary array image (32-bit float pixels) @@ -162,16 +160,16 @@ namespace search { fits_close_file(fptr, &status); fits_report_error(stderr, status); - } +} - void RawImage::append_layer_to_file(const std::string& filename) { +void RawImage::append_layer_to_file(const std::string& filename) { int status = 0; fitsfile* f; // Check that we can open the file. if (fits_open_file(&f, filename.c_str(), READWRITE, &status)) { - fits_report_error(stderr, status); - throw std::runtime_error("Unable to open FITS file for appending."); + fits_report_error(stderr, status); + throw std::runtime_error("Unable to open FITS file for appending."); } // This appends a layer (extension) if the file exists) @@ -192,131 +190,131 @@ namespace search { fits_close_file(f, &status); fits_report_error(stderr, status); - } +} - RawImage RawImage::create_stamp(float x, float y, int radius, bool interpolate, bool keep_no_data) const { +RawImage RawImage::create_stamp(float x, float y, int radius, bool interpolate, bool keep_no_data) const { if (radius < 0) throw std::runtime_error("stamp radius must be at least 0"); int dim = radius * 2 + 1; RawImage stamp(dim, dim); for (int yoff = 0; yoff < dim; ++yoff) { - for (int xoff = 0; xoff < dim; ++xoff) { - float pix_val; - if (interpolate) - pix_val = get_pixel_interp(x + static_cast(xoff - radius), - y + static_cast(yoff - radius)); - else - pix_val = get_pixel(static_cast(x) + xoff - radius, static_cast(y) + yoff - radius); - if ((pix_val == NO_DATA) && !keep_no_data) pix_val = 0.0; - stamp.set_pixel(xoff, yoff, pix_val); - } + for (int xoff = 0; xoff < dim; ++xoff) { + float pix_val; + if (interpolate) + pix_val = get_pixel_interp(x + static_cast(xoff - radius), + y + static_cast(yoff - radius)); + else + pix_val = get_pixel(static_cast(x) + xoff - radius, static_cast(y) + yoff - radius); + if ((pix_val == NO_DATA) && !keep_no_data) pix_val = 0.0; + stamp.set_pixel(xoff, yoff, pix_val); + } } stamp.set_obstime(obstime); return stamp; - } +} - void RawImage::convolve_cpu(const PSF& psf) { +void RawImage::convolve_cpu(const PSF& psf) { std::vector result(width * height, 0.0); const int psf_rad = psf.get_radius(); const float psf_total = psf.get_sum(); for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - // Pixels with NO_DATA remain NO_DATA. - if (pixels[y * width + x] == NO_DATA) { - result[y * width + x] = NO_DATA; - continue; - } + for (int x = 0; x < width; ++x) { + // Pixels with NO_DATA remain NO_DATA. + if (pixels[y * width + x] == NO_DATA) { + result[y * width + x] = NO_DATA; + continue; + } - float sum = 0.0; - float psf_portion = 0.0; - for (int j = -psf_rad; j <= psf_rad; j++) { - for (int i = -psf_rad; i <= psf_rad; i++) { - if ((x + i >= 0) && (x + i < width) && (y + j >= 0) && (y + j < height)) { - float current_pixel = pixels[(y + j) * width + (x + i)]; - if (current_pixel != NO_DATA) { - float current_psf = psf.get_value(i + psf_rad, j + psf_rad); - psf_portion += current_psf; - sum += current_pixel * current_psf; - } + float sum = 0.0; + float psf_portion = 0.0; + for (int j = -psf_rad; j <= psf_rad; j++) { + for (int i = -psf_rad; i <= psf_rad; i++) { + if ((x + i >= 0) && (x + i < width) && (y + j >= 0) && (y + j < height)) { + float current_pixel = pixels[(y + j) * width + (x + i)]; + if (current_pixel != NO_DATA) { + float current_psf = psf.get_value(i + psf_rad, j + psf_rad); + psf_portion += current_psf; + sum += current_pixel * current_psf; + } + } + } } - } + result[y * width + x] = (sum * psf_total) / psf_portion; } - result[y * width + x] = (sum * psf_total) / psf_portion; - } } // Copy the data into the pixels vector. const int npixels = get_npixels(); for (int i = 0; i < npixels; ++i) { - pixels[i] = result[i]; + pixels[i] = result[i]; } - } +} - void RawImage::convolve(PSF psf) { +void RawImage::convolve(PSF psf) { #ifdef HAVE_CUDA deviceConvolve(pixels.data(), pixels.data(), get_width(), get_height(), psf.data(), psf.get_size(), psf.get_dim(), psf.get_radius(), psf.get_sum()); #else convolve_cpu(psf); #endif - } +} - void RawImage::apply_mask(int flags, const std::vector& exceptions, const RawImage& mask) { +void RawImage::apply_mask(int flags, const std::vector& exceptions, const RawImage& mask) { const std::vector& mask_pix = mask.get_pixels(); const int num_pixels = get_npixels(); assert(num_pixels == mask.get_npixels()); for (unsigned int p = 0; p < num_pixels; ++p) { - int pix_flags = static_cast(mask_pix[p]); - bool is_exception = false; - for (auto& e : exceptions) is_exception = is_exception || e == pix_flags; - if (!is_exception && ((flags & pix_flags) != 0)) pixels[p] = NO_DATA; + int pix_flags = static_cast(mask_pix[p]); + bool is_exception = false; + for (auto& e : exceptions) is_exception = is_exception || e == pix_flags; + if (!is_exception && ((flags & pix_flags) != 0)) pixels[p] = NO_DATA; } - } +} - /* This implementation of grow_mask is optimized for steps > 1 - (which is how the code is generally used. If you are only - growing the mask by 1, the extra copy will be a little slower. - */ - void RawImage::grow_mask(int steps) { +/* This implementation of grow_mask is optimized for steps > 1 + (which is how the code is generally used. If you are only + growing the mask by 1, the extra copy will be a little slower. +*/ +void RawImage::grow_mask(int steps) { const int num_pixels = width * height; // Set up the initial masked vector that stores the number of steps // each pixel is from a masked pixel in the original image. std::vector masked(num_pixels, -1); for (int i = 0; i < num_pixels; ++i) { - if (pixels[i] == NO_DATA) masked[i] = 0; + if (pixels[i] == NO_DATA) masked[i] = 0; } // Grow out the mask one for each step. for (int itr = 1; itr <= steps; ++itr) { - for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - int center = width * y + x; - if (masked[center] == -1) { - // Mask pixels that are adjacent to a pixel masked during - // the last iteration only. - if ((x + 1 < width && masked[center + 1] == itr - 1) || - (x - 1 >= 0 && masked[center - 1] == itr - 1) || - (y + 1 < height && masked[center + width] == itr - 1) || - (y - 1 >= 0 && masked[center - width] == itr - 1)) { - masked[center] = itr; + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + int center = width * y + x; + if (masked[center] == -1) { + // Mask pixels that are adjacent to a pixel masked during + // the last iteration only. + if ((x + 1 < width && masked[center + 1] == itr - 1) || + (x - 1 >= 0 && masked[center - 1] == itr - 1) || + (y + 1 < height && masked[center + width] == itr - 1) || + (y - 1 >= 0 && masked[center - width] == itr - 1)) { + masked[center] = itr; + } + } } - } } - } } // Mask the pixels in the image. for (std::size_t i = 0; i < num_pixels; ++i) { - if (masked[i] > -1) { - pixels[i] = NO_DATA; - } + if (masked[i] > -1) { + pixels[i] = NO_DATA; + } } - } +} - std::vector RawImage::bilinear_interp(float x, float y) const { +std::vector RawImage::bilinear_interp(float x, float y) const { // Linear interpolation // Find the 4 pixels (aPix, bPix, cPix, dPix) // that the corners (a, b, c, d) of the @@ -357,9 +355,9 @@ namespace search { float diff = std::abs(a_amt + b_amt + c_amt + d_amt - 1.0); if (diff > 0.01) std::cout << "warning: bilinear_interpSum == " << diff << "\n"; return {a_px, a_py, a_amt, b_px, b_py, b_amt, c_px, c_py, c_amt, d_px, d_py, d_amt}; - } +} - void RawImage::add_pixel_interp(float x, float y, float value) { +void RawImage::add_pixel_interp(float x, float y, float value) { // Interpolation values std::vector iv = bilinear_interp(x, y); @@ -367,18 +365,18 @@ namespace search { add_to_pixel(iv[3], iv[4], value * iv[5]); add_to_pixel(iv[6], iv[7], value * iv[8]); add_to_pixel(iv[9], iv[10], value * iv[11]); - } +} - void RawImage::add_to_pixel(float fx, float fy, float value) { +void RawImage::add_to_pixel(float fx, float fy, float value) { assert(fx - floor(fx) == 0.0 && fy - floor(fy) == 0.0); int x = static_cast(fx); int y = static_cast(fy); if (x >= 0 && x < width && y >= 0 && y < height) pixels[y * width + x] += value; - } +} - float RawImage::get_pixel_interp(float x, float y) const { +float RawImage::get_pixel_interp(float x, float y) const { if ((x < 0.0 || y < 0.0) || (x > static_cast(width) || y > static_cast(height))) - return NO_DATA; + return NO_DATA; std::vector iv = bilinear_interp(x, y); float a = get_pixel(iv[0], iv[1]); float b = get_pixel(iv[3], iv[4]); @@ -387,41 +385,41 @@ namespace search { float interpSum = 0.0; float total = 0.0; if (a != NO_DATA) { - interpSum += iv[2]; - total += a * iv[2]; + interpSum += iv[2]; + total += a * iv[2]; } if (b != NO_DATA) { - interpSum += iv[5]; - total += b * iv[5]; + interpSum += iv[5]; + total += b * iv[5]; } if (c != NO_DATA) { - interpSum += iv[8]; - total += c * iv[8]; + interpSum += iv[8]; + total += c * iv[8]; } if (d != NO_DATA) { - interpSum += iv[11]; - total += d * iv[11]; + interpSum += iv[11]; + total += d * iv[11]; } if (interpSum == 0.0) { - return NO_DATA; + return NO_DATA; } else { - return total / interpSum; + return total / interpSum; } - } +} - void RawImage::set_all_pix(float value) { +void RawImage::set_all_pix(float value) { for (auto& p : pixels) p = value; - } +} - std::array RawImage::compute_bounds() const { +std::array RawImage::compute_bounds() const { const int num_pixels = get_npixels(); float min_val = FLT_MAX; float max_val = -FLT_MAX; for (unsigned p = 0; p < num_pixels; ++p) { - if (pixels[p] != NO_DATA) { - min_val = std::min(min_val, pixels[p]); - max_val = std::max(max_val, pixels[p]); - } + if (pixels[p] != NO_DATA) { + min_val = std::min(min_val, pixels[p]); + max_val = std::max(max_val, pixels[p]); + } } // Assert that we have seen at least some valid data. @@ -433,10 +431,10 @@ namespace search { res[0] = min_val; res[1] = max_val; return res; - } +} - // The maximum value of the image and return the coordinates. - PixelPos RawImage::find_peak(bool furthest_from_center) const { +// The maximum value of the image and return the coordinates. +PixelPos RawImage::find_peak(bool furthest_from_center) const { int c_x = width / 2; int c_y = height / 2; @@ -447,34 +445,34 @@ namespace search { // Search each pixel for the peak. for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - float pix_val = pixels[y * width + x]; - if (pix_val > max_val) { - max_val = pix_val; - result.x = x; - result.y = y; - dist2 = (c_x - x) * (c_x - x) + (c_y - y) * (c_y - y); - } else if (pix_val == max_val) { - int new_dist2 = (c_x - x) * (c_x - x) + (c_y - y) * (c_y - y); - if ((furthest_from_center && (new_dist2 > dist2)) || - (!furthest_from_center && (new_dist2 < dist2))) { - max_val = pix_val; - result.x = x; - result.y = y; - dist2 = new_dist2; - } + for (int x = 0; x < width; ++x) { + float pix_val = pixels[y * width + x]; + if (pix_val > max_val) { + max_val = pix_val; + result.x = x; + result.y = y; + dist2 = (c_x - x) * (c_x - x) + (c_y - y) * (c_y - y); + } else if (pix_val == max_val) { + int new_dist2 = (c_x - x) * (c_x - x) + (c_y - y) * (c_y - y); + if ((furthest_from_center && (new_dist2 > dist2)) || + (!furthest_from_center && (new_dist2 < dist2))) { + max_val = pix_val; + result.x = x; + result.y = y; + dist2 = new_dist2; + } + } } - } } return result; - } +} - // Find the basic image moments in order to test if stamps have a gaussian shape. - // It computes the moments on the "normalized" image where the minimum - // value has been shifted to zero and the sum of all elements is 1.0. - // Elements with NO_DATA are treated as zero. - ImageMoments RawImage::find_central_moments() const { +// Find the basic image moments in order to test if stamps have a gaussian shape. +// It computes the moments on the "normalized" image where the minimum +// value has been shifted to zero and the sum of all elements is 1.0. +// Elements with NO_DATA are treated as zero. +ImageMoments RawImage::find_central_moments() const { const int num_pixels = width * height; const int c_x = width / 2; const int c_y = height / 2; @@ -485,35 +483,35 @@ namespace search { // Find the min (non-NO_DATA) value to subtract off. float min_val = FLT_MAX; for (int p = 0; p < num_pixels; ++p) { - min_val = ((pixels[p] != NO_DATA) && (pixels[p] < min_val)) ? pixels[p] : min_val; + min_val = ((pixels[p] != NO_DATA) && (pixels[p] < min_val)) ? pixels[p] : min_val; } // Find the sum of the zero-shifted (non-NO_DATA) pixels. double sum = 0.0; for (int p = 0; p < num_pixels; ++p) { - sum += (pixels[p] != NO_DATA) ? (pixels[p] - min_val) : 0.0; + sum += (pixels[p] != NO_DATA) ? (pixels[p] - min_val) : 0.0; } if (sum == 0.0) return res; // Compute the rest of the moments. for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - int ind = y * width + x; - float pix_val = (pixels[ind] != NO_DATA) ? (pixels[ind] - min_val) / sum : 0.0; - - res.m00 += pix_val; - res.m10 += (x - c_x) * pix_val; - res.m20 += (x - c_x) * (x - c_x) * pix_val; - res.m01 += (y - c_y) * pix_val; - res.m02 += (y - c_y) * (y - c_y) * pix_val; - res.m11 += (x - c_x) * (y - c_y) * pix_val; - } + for (int x = 0; x < width; ++x) { + int ind = y * width + x; + float pix_val = (pixels[ind] != NO_DATA) ? (pixels[ind] - min_val) / sum : 0.0; + + res.m00 += pix_val; + res.m10 += (x - c_x) * pix_val; + res.m20 += (x - c_x) * (x - c_x) * pix_val; + res.m01 += (y - c_y) * pix_val; + res.m02 += (y - c_y) * (y - c_y) * pix_val; + res.m11 += (x - c_x) * (y - c_y) * pix_val; + } } return res; - } +} - RawImage create_median_image(const std::vector& images) { +RawImage create_median_image(const std::vector& images) { int num_images = images.size(); assert(num_images > 0); @@ -524,41 +522,41 @@ namespace search { RawImage result = RawImage(width, height); std::vector pix_array(num_images); for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - int num_unmasked = 0; - for (int i = 0; i < num_images; ++i) { - // Only used the unmasked pixels. - float pix_val = images[i].get_pixel(x, y); - if ((pix_val != NO_DATA) && (!std::isnan(pix_val))) { - pix_array[num_unmasked] = pix_val; - num_unmasked += 1; - } - } + for (int x = 0; x < width; ++x) { + int num_unmasked = 0; + for (int i = 0; i < num_images; ++i) { + // Only used the unmasked pixels. + float pix_val = images[i].get_pixel(x, y); + if ((pix_val != NO_DATA) && (!std::isnan(pix_val))) { + pix_array[num_unmasked] = pix_val; + num_unmasked += 1; + } + } - if (num_unmasked > 0) { - std::sort(pix_array.begin(), pix_array.begin() + num_unmasked); - - // If we have an even number of elements, take the mean of the two - // middle ones. - int median_ind = num_unmasked / 2; - if (num_unmasked % 2 == 0) { - float ave_middle = (pix_array[median_ind] + pix_array[median_ind - 1]) / 2.0; - result.set_pixel(x, y, ave_middle); - } else { - result.set_pixel(x, y, pix_array[median_ind]); - } - } else { - // We use a 0.0 value if there is no data to allow for visualization - // and value based filtering. - result.set_pixel(x, y, 0.0); + if (num_unmasked > 0) { + std::sort(pix_array.begin(), pix_array.begin() + num_unmasked); + + // If we have an even number of elements, take the mean of the two + // middle ones. + int median_ind = num_unmasked / 2; + if (num_unmasked % 2 == 0) { + float ave_middle = (pix_array[median_ind] + pix_array[median_ind - 1]) / 2.0; + result.set_pixel(x, y, ave_middle); + } else { + result.set_pixel(x, y, pix_array[median_ind]); + } + } else { + // We use a 0.0 value if there is no data to allow for visualization + // and value based filtering. + result.set_pixel(x, y, 0.0); + } } - } } return result; - } +} - RawImage create_summed_image(const std::vector& images) { +RawImage create_summed_image(const std::vector& images) { int num_images = images.size(); assert(num_images > 0); @@ -568,21 +566,21 @@ namespace search { RawImage result = RawImage(width, height); for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - float sum = 0.0; - for (int i = 0; i < num_images; ++i) { - float pix_val = images[i].get_pixel(x, y); - if ((pix_val == NO_DATA) || (std::isnan(pix_val))) pix_val = 0.0; - sum += pix_val; + for (int x = 0; x < width; ++x) { + float sum = 0.0; + for (int i = 0; i < num_images; ++i) { + float pix_val = images[i].get_pixel(x, y); + if ((pix_val == NO_DATA) || (std::isnan(pix_val))) pix_val = 0.0; + sum += pix_val; + } + result.set_pixel(x, y, sum); } - result.set_pixel(x, y, sum); - } } return result; - } +} - RawImage create_mean_image(const std::vector& images) { +RawImage create_mean_image(const std::vector& images) { int num_images = images.size(); assert(num_images > 0); @@ -592,38 +590,37 @@ namespace search { RawImage result = RawImage(width, height); for (int y = 0; y < height; ++y) { - for (int x = 0; x < width; ++x) { - float sum = 0.0; - float count = 0.0; - for (int i = 0; i < num_images; ++i) { - float pix_val = images[i].get_pixel(x, y); - if ((pix_val != NO_DATA) && (!std::isnan(pix_val))) { - count += 1.0; - sum += pix_val; - } - } + for (int x = 0; x < width; ++x) { + float sum = 0.0; + float count = 0.0; + for (int i = 0; i < num_images; ++i) { + float pix_val = images[i].get_pixel(x, y); + if ((pix_val != NO_DATA) && (!std::isnan(pix_val))) { + count += 1.0; + sum += pix_val; + } + } - if (count > 0.0) { - result.set_pixel(x, y, sum / count); - } else { - // We use a 0.0 value if there is no data to allow for visualization - // and value based filtering. - result.set_pixel(x, y, 0.0); + if (count > 0.0) { + result.set_pixel(x, y, sum / count); + } else { + // We use a 0.0 value if there is no data to allow for visualization + // and value based filtering. + result.set_pixel(x, y, 0.0); + } } - } } return result; - } - +} #ifdef Py_PYTHON_H - RawImage::RawImage(pybind11::array_t arr) { +RawImage::RawImage(pybind11::array_t arr) { obstime = -1.0; set_array(arr); - } +} - void RawImage::set_array(pybind11::array_t& arr) { +void RawImage::set_array(pybind11::array_t& arr) { pybind11::buffer_info info = arr.request(); if (info.ndim != 2) throw std::runtime_error("Array must have 2 dimensions."); @@ -633,48 +630,47 @@ namespace search { float* pix = static_cast(info.ptr); pixels = std::vector(pix, pix + get_npixels()); - } - +} - static void raw_image_bindings(py::module &m) { +static void raw_image_bindings(py::module& m) { using ri = search::RawImage; py::class_(m, "RawImage", py::buffer_protocol(), pydocs::DOC_RawImage) - .def_buffer([](ri &m) -> py::buffer_info { - return py::buffer_info(m.data(), sizeof(float), py::format_descriptor::format(), - 2, {m.get_height(), m.get_width()}, - {sizeof(float) * m.get_width(), sizeof(float)}); - }) - .def(py::init()) - .def(py::init()) - .def(py::init>()) - .def("get_height", &ri::get_height, pydocs::DOC_RawImage_get_height ) - .def("get_width", &ri::get_width, pydocs:: DOC_RawImage_get_width) - .def("get_npixels", &ri::get_npixels, pydocs:: DOC_RawImage_get_npixels) - .def("get_all_pixels", &ri::get_pixels, pydocs:: DOC_RawImage_get_all_pixels) - .def("set_array", &ri::set_array, pydocs:: DOC_RawImage_set_array) - .def("get_obstime", &ri::get_obstime, pydocs:: DOC_RawImage_get_obstime) - .def("set_obstime", &ri::set_obstime, pydocs:: DOC_RawImage_set_obstime) - .def("approx_equal", &ri::approx_equal, pydocs:: DOC_RawImage_approx_equal) - .def("compute_bounds", &ri::compute_bounds, pydocs:: DOC_RawImage_compute_bounds) - .def("find_peak", &ri::find_peak, pydocs:: DOC_RawImage_find_peak) - .def("find_central_moments", &ri::find_central_moments, pydocs:: DOC_RawImage_find_central_moments) - .def("create_stamp", &ri::create_stamp, pydocs:: DOC_RawImage_create_stamp) - .def("set_pixel", &ri::set_pixel, pydocs:: DOC_RawImage_set_pixel) - .def("add_pixel", &ri::add_to_pixel, pydocs:: DOC_RawImage_add_pixel) - .def("add_pixel_interp", &ri::add_pixel_interp, pydocs:: DOC_RawImage_add_pixel_interp) - .def("apply_mask", &ri::apply_mask, pydocs:: DOC_RawImage_apply_mask) - .def("grow_mask", &ri::grow_mask, pydocs:: DOC_RawImage_grow_mask) - .def("pixel_has_data", &ri::pixel_has_data, pydocs:: DOC_RawImage_pixel_has_data) - .def("set_all", &ri::set_all_pix, pydocs:: DOC_RawImage_set_all) - .def("get_pixel", &ri::get_pixel, pydocs:: DOC_RawImage_get_pixel) - .def("get_pixel_interp", &ri::get_pixel_interp, pydocs:: DOC_RawImage_get_pixel_interp) - .def("convolve", &ri::convolve, pydocs:: DOC_RawImage_convolve) - .def("convolve_cpu", &ri::convolve_cpu, pydocs:: DOC_RawImage_convolve_cpu) - .def("load_fits", &ri::load_from_file, pydocs:: DOC_RawImage_load_fits) - .def("save_fits", &ri::save_to_file, pydocs:: DOC_RawImage_save_fits) - .def("append_fits_layer", &ri::append_layer_to_file, pydocs:: DOC_RawImage_append_fits_layer); - } + .def_buffer([](ri& m) -> py::buffer_info { + return py::buffer_info(m.data(), sizeof(float), py::format_descriptor::format(), 2, + {m.get_height(), m.get_width()}, + {sizeof(float) * m.get_width(), sizeof(float)}); + }) + .def(py::init()) + .def(py::init()) + .def(py::init>()) + .def("get_height", &ri::get_height, pydocs::DOC_RawImage_get_height) + .def("get_width", &ri::get_width, pydocs::DOC_RawImage_get_width) + .def("get_npixels", &ri::get_npixels, pydocs::DOC_RawImage_get_npixels) + .def("get_all_pixels", &ri::get_pixels, pydocs::DOC_RawImage_get_all_pixels) + .def("set_array", &ri::set_array, pydocs::DOC_RawImage_set_array) + .def("get_obstime", &ri::get_obstime, pydocs::DOC_RawImage_get_obstime) + .def("set_obstime", &ri::set_obstime, pydocs::DOC_RawImage_set_obstime) + .def("approx_equal", &ri::approx_equal, pydocs::DOC_RawImage_approx_equal) + .def("compute_bounds", &ri::compute_bounds, pydocs::DOC_RawImage_compute_bounds) + .def("find_peak", &ri::find_peak, pydocs::DOC_RawImage_find_peak) + .def("find_central_moments", &ri::find_central_moments, pydocs::DOC_RawImage_find_central_moments) + .def("create_stamp", &ri::create_stamp, pydocs::DOC_RawImage_create_stamp) + .def("set_pixel", &ri::set_pixel, pydocs::DOC_RawImage_set_pixel) + .def("add_pixel", &ri::add_to_pixel, pydocs::DOC_RawImage_add_pixel) + .def("add_pixel_interp", &ri::add_pixel_interp, pydocs::DOC_RawImage_add_pixel_interp) + .def("apply_mask", &ri::apply_mask, pydocs::DOC_RawImage_apply_mask) + .def("grow_mask", &ri::grow_mask, pydocs::DOC_RawImage_grow_mask) + .def("pixel_has_data", &ri::pixel_has_data, pydocs::DOC_RawImage_pixel_has_data) + .def("set_all", &ri::set_all_pix, pydocs::DOC_RawImage_set_all) + .def("get_pixel", &ri::get_pixel, pydocs::DOC_RawImage_get_pixel) + .def("get_pixel_interp", &ri::get_pixel_interp, pydocs::DOC_RawImage_get_pixel_interp) + .def("convolve", &ri::convolve, pydocs::DOC_RawImage_convolve) + .def("convolve_cpu", &ri::convolve_cpu, pydocs::DOC_RawImage_convolve_cpu) + .def("load_fits", &ri::load_from_file, pydocs::DOC_RawImage_load_fits) + .def("save_fits", &ri::save_to_file, pydocs::DOC_RawImage_save_fits) + .def("append_fits_layer", &ri::append_layer_to_file, pydocs::DOC_RawImage_append_fits_layer); +} #endif } /* namespace search */ diff --git a/src/kbmod/search/raw_image.h b/src/kbmod/search/raw_image.h index a8b4e66c..0742eba5 100644 --- a/src/kbmod/search/raw_image.h +++ b/src/kbmod/search/raw_image.h @@ -14,10 +14,9 @@ #include "psf.h" #include "pydocs/raw_image_docs.h" - namespace search { - class RawImage { - public: +class RawImage { +public: RawImage(); RawImage(const RawImage& old); // Copy constructor RawImage(RawImage&& source); // Move constructor @@ -39,15 +38,15 @@ namespace search { // Inline pixel functions. float get_pixel(int x, int y) const { - return (x >= 0 && x < width && y >= 0 && y < height) ? pixels[y * width + x] : NO_DATA; + return (x >= 0 && x < width && y >= 0 && y < height) ? pixels[y * width + x] : NO_DATA; } bool pixel_has_data(int x, int y) const { - return (x >= 0 && x < width && y >= 0 && y < height) ? pixels[y * width + x] != NO_DATA : false; + return (x >= 0 && x < width && y >= 0 && y < height) ? pixels[y * width + x] != NO_DATA : false; } void set_pixel(int x, int y, float value) { - if (x >= 0 && x < width && y >= 0 && y < height) pixels[y * width + x] = value; + if (x >= 0 && x < width && y >= 0 && y < height) pixels[y * width + x] = value; } const std::vector& get_pixels() const { return pixels; } float* data() { return pixels.data(); } // Get pointer to pixels @@ -111,19 +110,19 @@ namespace search { virtual ~RawImage(){}; - private: +private: void load_time_from_file(fitsfile* fptr); - + unsigned width; unsigned height; std::vector pixels; double obstime; - }; +}; - // Helper functions for creating composite images. - RawImage create_median_image(const std::vector& images); - RawImage create_summed_image(const std::vector& images); - RawImage create_mean_image(const std::vector& images); +// Helper functions for creating composite images. +RawImage create_median_image(const std::vector& images); +RawImage create_summed_image(const std::vector& images); +RawImage create_mean_image(const std::vector& images); } /* namespace search */ diff --git a/src/kbmod/search/stack_search.cpp b/src/kbmod/search/stack_search.cpp index 6a8e9168..59fc38cd 100644 --- a/src/kbmod/search/stack_search.cpp +++ b/src/kbmod/search/stack_search.cpp @@ -1,18 +1,17 @@ #include "stack_search.h" - namespace search { #ifdef HAVE_CUDA - extern "C" void deviceSearchFilter(int num_images, int width, int height, float* psi_vect, float* phi_vect, - PerImageData img_data, SearchParameters params, int num_trajectories, - Trajectory* trj_to_search, int num_results, Trajectory* best_results); +extern "C" void deviceSearchFilter(int num_images, int width, int height, float* psi_vect, float* phi_vect, + PerImageData img_data, SearchParameters params, int num_trajectories, + Trajectory* trj_to_search, int num_results, Trajectory* best_results); - void deviceGetCoadds(ImageStack& stack, PerImageData image_data, int num_trajectories, - Trajectory* trajectories, StampParameters params, - std::vector >& use_index_vect, float* results); +void deviceGetCoadds(ImageStack& stack, PerImageData image_data, int num_trajectories, + Trajectory* trajectories, StampParameters params, + std::vector>& use_index_vect, float* results); #endif - StackSearch::StackSearch(ImageStack& imstack) : stack(imstack) { +StackSearch::StackSearch(ImageStack& imstack) : stack(imstack) { debug_info = false; psi_phi_generated = false; @@ -37,49 +36,48 @@ namespace search { params.y_start_max = stack.get_height(); params.debug = false; - } +} - void StackSearch::set_debug(bool d) { +void StackSearch::set_debug(bool d) { debug_info = d; params.debug = d; - } +} - void StackSearch::enable_gpu_sigmag_filter(std::vector percentiles, float sigmag_coeff, - float min_lh) { +void StackSearch::enable_gpu_sigmag_filter(std::vector percentiles, float sigmag_coeff, float min_lh) { params.do_sigmag_filter = true; params.sgl_L = percentiles[0]; params.sgl_H = percentiles[1]; params.sigmag_coeff = sigmag_coeff; params.min_lh = min_lh; - } +} - void StackSearch::enable_gpu_encoding(int psi_num_bytes, int phi_num_bytes) { +void StackSearch::enable_gpu_encoding(int psi_num_bytes, int phi_num_bytes) { // Make sure the encoding is one of the supported options. // Otherwise use default float (aka no encoding). if (psi_num_bytes == 1 || psi_num_bytes == 2) { - params.psi_num_bytes = psi_num_bytes; + params.psi_num_bytes = psi_num_bytes; } else { - params.psi_num_bytes = -1; + params.psi_num_bytes = -1; } if (phi_num_bytes == 1 || phi_num_bytes == 2) { - params.phi_num_bytes = phi_num_bytes; + params.phi_num_bytes = phi_num_bytes; } else { - params.phi_num_bytes = -1; + params.phi_num_bytes = -1; } - } +} - void StackSearch::set_start_bounds_x(int x_min, int x_max) { +void StackSearch::set_start_bounds_x(int x_min, int x_max) { params.x_start_min = x_min; params.x_start_max = x_max; - } +} - void StackSearch::set_start_bounds_y(int y_min, int y_max) { +void StackSearch::set_start_bounds_y(int y_min, int y_max) { params.y_start_min = y_min; params.y_start_max = y_max; - } +} - void StackSearch::search(int ang_steps, int vel_steps, float min_ang, float max_ang, float min_vel, - float mavx, int min_observations) { +void StackSearch::search(int ang_steps, int vel_steps, float min_ang, float max_ang, float min_vel, + float mavx, int min_observations) { prepare_psi_phi(); create_search_list(ang_steps, vel_steps, min_ang, max_ang, min_vel, mavx); @@ -100,22 +98,22 @@ namespace search { std::vector psi_scale_vect; std::vector phi_scale_vect; if (params.psi_num_bytes > 0) { - psi_scale_vect = compute_image_scaling(psi_images, params.psi_num_bytes); - img_data.psi_params = psi_scale_vect.data(); + psi_scale_vect = compute_image_scaling(psi_images, params.psi_num_bytes); + img_data.psi_params = psi_scale_vect.data(); } if (params.phi_num_bytes > 0) { - phi_scale_vect = compute_image_scaling(phi_images, params.phi_num_bytes); - img_data.phi_params = phi_scale_vect.data(); + phi_scale_vect = compute_image_scaling(phi_images, params.phi_num_bytes); + img_data.phi_params = phi_scale_vect.data(); } // Allocate a vector for the results. int num_search_pixels = - ((params.x_start_max - params.x_start_min) * (params.y_start_max - params.y_start_min)); + ((params.x_start_max - params.x_start_min) * (params.y_start_max - params.y_start_min)); int max_results = num_search_pixels * RESULTS_PER_PIXEL; if (debug_info) { - std::cout << "Searching X=[" << params.x_start_min << ", " << params.x_start_max << "]" - << " Y=[" << params.y_start_min << ", " << params.y_start_max << "]\n"; - std::cout << "Allocating space for " << max_results << " results.\n"; + std::cout << "Searching X=[" << params.x_start_min << ", " << params.x_start_max << "]" + << " Y=[" << params.y_start_min << ", " << params.y_start_max << "]\n"; + std::cout << "Allocating space for " << max_results << " results.\n"; } results = std::vector(max_results); if (debug_info) std::cout << search_list.size() << " trajectories... \n" << std::flush; @@ -126,8 +124,9 @@ namespace search { // Do the actual search on the GPU. start_timer("Searching"); #ifdef HAVE_CUDA - deviceSearchFilter(stack.img_count(), stack.get_width(), stack.get_height(), psi_vect.data(), phi_vect.data(), - img_data, params, search_list.size(), search_list.data(), max_results, results.data()); + deviceSearchFilter(stack.img_count(), stack.get_width(), stack.get_height(), psi_vect.data(), + phi_vect.data(), img_data, params, search_list.size(), search_list.data(), max_results, + results.data()); #else throw std::runtime_error("Non-GPU search is not implemented."); #endif @@ -136,98 +135,97 @@ namespace search { start_timer("Sorting results"); sort_results(); end_timer(); - } +} - void StackSearch::save_psiphi(const std::string& path) { +void StackSearch::save_psiphi(const std::string& path) { prepare_psi_phi(); save_images(path); - } +} - void StackSearch::prepare_psi_phi() { +void StackSearch::prepare_psi_phi() { if (!psi_phi_generated) { - psi_images.clear(); - phi_images.clear(); - - // Compute Phi and Psi from convolved images - // while leaving masked pixels alone - // Reinsert 0s for NO_DATA? - const int num_images = stack.img_count(); - for (int i = 0; i < num_images; ++i) { - LayeredImage& img = stack.get_single_image(i); - psi_images.push_back(img.generate_psi_image()); - phi_images.push_back(img.generate_phi_image()); - } - - psi_phi_generated = true; + psi_images.clear(); + phi_images.clear(); + + // Compute Phi and Psi from convolved images + // while leaving masked pixels alone + // Reinsert 0s for NO_DATA? + const int num_images = stack.img_count(); + for (int i = 0; i < num_images; ++i) { + LayeredImage& img = stack.get_single_image(i); + psi_images.push_back(img.generate_psi_image()); + phi_images.push_back(img.generate_phi_image()); + } + + psi_phi_generated = true; } - } +} - std::vector StackSearch::compute_image_scaling(const std::vector& vect, - int encoding_bytes) const { +std::vector StackSearch::compute_image_scaling(const std::vector& vect, + int encoding_bytes) const { std::vector result; const int num_images = vect.size(); for (int i = 0; i < num_images; ++i) { - scaleParameters params; - params.scale = 1.0; + scaleParameters params; + params.scale = 1.0; - std::array bnds = vect[i].compute_bounds(); - params.min_val = bnds[0]; - params.max_val = bnds[1]; + std::array bnds = vect[i].compute_bounds(); + params.min_val = bnds[0]; + params.max_val = bnds[1]; - // Increase width to avoid divide by zero. - float width = (params.max_val - params.min_val); - if (width < 1e-6) width = 1e-6; + // Increase width to avoid divide by zero. + float width = (params.max_val - params.min_val); + if (width < 1e-6) width = 1e-6; - // Set the scale if we are encoding the values. - if (encoding_bytes == 1 || encoding_bytes == 2) { - long int num_values = (1 << (8 * encoding_bytes)) - 1; - params.scale = width / (double)num_values; - } + // Set the scale if we are encoding the values. + if (encoding_bytes == 1 || encoding_bytes == 2) { + long int num_values = (1 << (8 * encoding_bytes)) - 1; + params.scale = width / (double)num_values; + } - result.push_back(params); + result.push_back(params); } return result; - } +} - void StackSearch::save_images(const std::string& path) { +void StackSearch::save_images(const std::string& path) { for (int i = 0; i < stack.img_count(); ++i) { - std::string number = std::to_string(i); - // Add leading zeros - number = std::string(4 - number.length(), '0') + number; - psi_images[i].save_to_file(path + "/psi/PSI" + number + ".fits"); - phi_images[i].save_to_file(path + "/phi/PHI" + number + ".fits"); + std::string number = std::to_string(i); + // Add leading zeros + number = std::string(4 - number.length(), '0') + number; + psi_images[i].save_to_file(path + "/psi/PSI" + number + ".fits"); + phi_images[i].save_to_file(path + "/phi/PHI" + number + ".fits"); } - } +} - void StackSearch::create_search_list(int angle_steps, int velocity_steps, float min_ang, float max_ang, - float min_vel, float mavx) { +void StackSearch::create_search_list(int angle_steps, int velocity_steps, float min_ang, float max_ang, + float min_vel, float mavx) { std::vector angles(angle_steps); float ang_stepsize = (max_ang - min_ang) / float(angle_steps); for (int i = 0; i < angle_steps; ++i) { - angles[i] = min_ang + float(i) * ang_stepsize; + angles[i] = min_ang + float(i) * ang_stepsize; } std::vector velocities(velocity_steps); float vel_stepsize = (mavx - min_vel) / float(velocity_steps); for (int i = 0; i < velocity_steps; ++i) { - velocities[i] = min_vel + float(i) * vel_stepsize; + velocities[i] = min_vel + float(i) * vel_stepsize; } int trajCount = angle_steps * velocity_steps; search_list = std::vector(trajCount); for (int a = 0; a < angle_steps; ++a) { - for (int v = 0; v < velocity_steps; ++v) { - search_list[a * velocity_steps + v].vx = cos(angles[a]) * velocities[v]; - search_list[a * velocity_steps + v].vy = sin(angles[a]) * velocities[v]; - } + for (int v = 0; v < velocity_steps; ++v) { + search_list[a * velocity_steps + v].vx = cos(angles[a]) * velocities[v]; + search_list[a * velocity_steps + v].vy = sin(angles[a]) * velocities[v]; + } } - } +} - void StackSearch::fill_psi_phi(const std::vector& psi_imgs, - const std::vector& phi_imgs, std::vector* psi_vect, - std::vector* phi_vect) { +void StackSearch::fill_psi_phi(const std::vector& psi_imgs, const std::vector& phi_imgs, + std::vector* psi_vect, std::vector* phi_vect) { assert(psi_vect != NULL); assert(phi_vect != NULL); @@ -237,8 +235,8 @@ namespace search { int num_pixels = psi_imgs[0].get_npixels(); for (int i = 0; i < num_images; ++i) { - assert(psi_imgs[i].get_npixels() == num_pixels); - assert(phi_imgs[i].get_npixels() == num_pixels); + assert(psi_imgs[i].get_npixels() == num_pixels); + assert(phi_imgs[i].get_npixels() == num_pixels); } psi_vect->clear(); @@ -247,66 +245,66 @@ namespace search { phi_vect->reserve(num_images * num_pixels); for (int i = 0; i < num_images; ++i) { - const std::vector& psi_ref = psi_imgs[i].get_pixels(); - const std::vector& phi_ref = phi_imgs[i].get_pixels(); - for (unsigned p = 0; p < num_pixels; ++p) { - psi_vect->push_back(psi_ref[p]); - phi_vect->push_back(phi_ref[p]); - } + const std::vector& psi_ref = psi_imgs[i].get_pixels(); + const std::vector& phi_ref = phi_imgs[i].get_pixels(); + for (unsigned p = 0; p < num_pixels; ++p) { + psi_vect->push_back(psi_ref[p]); + phi_vect->push_back(phi_ref[p]); + } } - } +} - std::vector StackSearch::create_stamps(const Trajectory& trj, int radius, bool interpolate, - bool keep_no_data, const std::vector& use_index) { +std::vector StackSearch::create_stamps(const Trajectory& trj, int radius, bool interpolate, + bool keep_no_data, const std::vector& use_index) { if (use_index.size() > 0 && use_index.size() != stack.img_count()) { - throw std::runtime_error("Wrong size use_index passed into create_stamps()"); + throw std::runtime_error("Wrong size use_index passed into create_stamps()"); } bool use_all_stamps = use_index.size() == 0; std::vector stamps; int num_times = stack.img_count(); for (int i = 0; i < num_times; ++i) { - if (use_all_stamps || use_index[i]) { - PixelPos pos = get_trajectory_position(trj, i); - RawImage& img = stack.get_single_image(i).get_science(); - stamps.push_back(img.create_stamp(pos.x, pos.y, radius, interpolate, keep_no_data)); - } + if (use_all_stamps || use_index[i]) { + PixelPos pos = get_trajectory_position(trj, i); + RawImage& img = stack.get_single_image(i).get_science(); + stamps.push_back(img.create_stamp(pos.x, pos.y, radius, interpolate, keep_no_data)); + } } return stamps; - } +} - // For stamps used for visualization we interpolate the pixel values, replace - // NO_DATA tages with zeros, and return all the stamps (regardless of whether - // individual timesteps have been filtered). - std::vector StackSearch::get_stamps(const Trajectory& t, int radius) { +// For stamps used for visualization we interpolate the pixel values, replace +// NO_DATA tages with zeros, and return all the stamps (regardless of whether +// individual timesteps have been filtered). +std::vector StackSearch::get_stamps(const Trajectory& t, int radius) { std::vector empty_vect; return create_stamps(t, radius, true /*=interpolate*/, false /*=keep_no_data*/, empty_vect); - } +} - // For creating coadded stamps, we do not interpolate the pixel values and keep - // NO_DATA tagged (so we can filter it out of mean/median). - RawImage StackSearch::get_median_stamp(const Trajectory& trj, int radius, - const std::vector& use_index) { +// For creating coadded stamps, we do not interpolate the pixel values and keep +// NO_DATA tagged (so we can filter it out of mean/median). +RawImage StackSearch::get_median_stamp(const Trajectory& trj, int radius, + const std::vector& use_index) { return create_median_image( - create_stamps(trj, radius, false /*=interpolate*/, true /*=keep_no_data*/, use_index)); - } + create_stamps(trj, radius, false /*=interpolate*/, true /*=keep_no_data*/, use_index)); +} - // For creating coadded stamps, we do not interpolate the pixel values and keep - // NO_DATA tagged (so we can filter it out of mean/median). - RawImage StackSearch::get_mean_stamp(const Trajectory& trj, int radius, const std::vector& use_index) { +// For creating coadded stamps, we do not interpolate the pixel values and keep +// NO_DATA tagged (so we can filter it out of mean/median). +RawImage StackSearch::get_mean_stamp(const Trajectory& trj, int radius, const std::vector& use_index) { return create_mean_image( - create_stamps(trj, radius, false /*=interpolate*/, true /*=keep_no_data*/, use_index)); - } + create_stamps(trj, radius, false /*=interpolate*/, true /*=keep_no_data*/, use_index)); +} - // For creating summed stamps, we do not interpolate the pixel values and replace NO_DATA - // with zero (which is the same as filtering it out for the sum). - RawImage StackSearch::get_summed_stamp(const Trajectory& trj, int radius, - const std::vector& use_index) { +// For creating summed stamps, we do not interpolate the pixel values and replace NO_DATA +// with zero (which is the same as filtering it out for the sum). +RawImage StackSearch::get_summed_stamp(const Trajectory& trj, int radius, + const std::vector& use_index) { return create_summed_image( - create_stamps(trj, radius, false /*=interpolate*/, false /*=keep_no_data*/, use_index)); - } + create_stamps(trj, radius, false /*=interpolate*/, false /*=keep_no_data*/, use_index)); +} - bool StackSearch::filter_stamp(const RawImage& img, const StampParameters& params) { +bool StackSearch::filter_stamp(const RawImage& img, const StampParameters& params) { // Allocate space for the coadd information and initialize to zero. const int stamp_width = 2 * params.radius + 1; const int stamp_ppi = stamp_width * stamp_width; @@ -316,21 +314,21 @@ namespace search { PixelPos pos = img.find_peak(true); if ((abs(pos.x - params.radius) >= params.peak_offset_x) || (abs(pos.y - params.radius) >= params.peak_offset_y)) { - return true; + return true; } // Filter on the percentage of flux in the central pixel. if (params.center_thresh > 0.0) { - const std::vector& pixels = img.get_pixels(); - float center_val = pixels[(int)pos.y * stamp_width + (int)pos.x]; - float pixel_sum = 0.0; - for (int p = 0; p < stamp_ppi; ++p) { - pixel_sum += pixels[p]; - } - - if (center_val / pixel_sum < params.center_thresh) { - return true; - } + const std::vector& pixels = img.get_pixels(); + float center_val = pixels[(int)pos.y * stamp_width + (int)pos.x]; + float pixel_sum = 0.0; + for (int p = 0; p < stamp_ppi; ++p) { + pixel_sum += pixels[p]; + } + + if (center_val / pixel_sum < params.center_thresh) { + return true; + } } // Filter on the image moments. @@ -338,69 +336,69 @@ namespace search { if ((fabs(moments.m01) >= params.m01_limit) || (fabs(moments.m10) >= params.m10_limit) || (fabs(moments.m11) >= params.m11_limit) || (moments.m02 >= params.m02_limit) || (moments.m20 >= params.m20_limit)) { - return true; + return true; } return false; - } +} - std::vector StackSearch::get_coadded_stamps(std::vector& t_array, - std::vector >& use_index_vect, - const StampParameters& params, bool use_gpu) { +std::vector StackSearch::get_coadded_stamps(std::vector& t_array, + std::vector>& use_index_vect, + const StampParameters& params, bool use_gpu) { if (use_gpu) { #ifdef HAVE_CUDA - return get_coadded_stamps_gpu(t_array, use_index_vect, params); + return get_coadded_stamps_gpu(t_array, use_index_vect, params); #else - std::cout << "WARNING: GPU is not enabled. Performing co-adds on the CPU."; - //py::print("WARNING: GPU is not enabled. Performing co-adds on the CPU."); + std::cout << "WARNING: GPU is not enabled. Performing co-adds on the CPU."; + // py::print("WARNING: GPU is not enabled. Performing co-adds on the CPU."); #endif } return get_coadded_stamps_cpu(t_array, use_index_vect, params); - } +} - std::vector StackSearch::get_coadded_stamps_cpu(std::vector& t_array, - std::vector >& use_index_vect, - const StampParameters& params) { +std::vector StackSearch::get_coadded_stamps_cpu(std::vector& t_array, + std::vector>& use_index_vect, + const StampParameters& params) { const int num_trajectories = t_array.size(); std::vector results(num_trajectories); std::vector empty_pixels(1, NO_DATA); for (int i = 0; i < num_trajectories; ++i) { - std::vector stamps = - create_stamps(t_array[i], params.radius, false, true, use_index_vect[i]); - - RawImage coadd(1, 1); - switch (params.stamp_type) { - case STAMP_MEDIAN: - coadd = create_median_image(stamps); - break; - case STAMP_MEAN: - coadd = create_mean_image(stamps); - break; - case STAMP_SUM: - coadd = create_summed_image(stamps); - break; - default: - throw std::runtime_error("Invalid stamp coadd type."); - } - - // Do the filtering if needed. - if (params.do_filtering && filter_stamp(coadd, params)) { - results[i] = RawImage(1, 1, empty_pixels); - } else { - results[i] = coadd; - } + std::vector stamps = + create_stamps(t_array[i], params.radius, false, true, use_index_vect[i]); + + RawImage coadd(1, 1); + switch (params.stamp_type) { + case STAMP_MEDIAN: + coadd = create_median_image(stamps); + break; + case STAMP_MEAN: + coadd = create_mean_image(stamps); + break; + case STAMP_SUM: + coadd = create_summed_image(stamps); + break; + default: + throw std::runtime_error("Invalid stamp coadd type."); + } + + // Do the filtering if needed. + if (params.do_filtering && filter_stamp(coadd, params)) { + results[i] = RawImage(1, 1, empty_pixels); + } else { + results[i] = coadd; + } } return results; - } +} - std::vector StackSearch::get_coadded_stamps_gpu(std::vector& t_array, - std::vector >& use_index_vect, - const StampParameters& params) { +std::vector StackSearch::get_coadded_stamps_gpu(std::vector& t_array, + std::vector>& use_index_vect, + const StampParameters& params) { // Right now only limited stamp sizes are allowed. if (2 * params.radius + 1 > MAX_STAMP_EDGE || params.radius <= 0) { - throw std::runtime_error("Invalid Radius."); + throw std::runtime_error("Invalid Radius."); } const int num_images = stack.img_count(); @@ -432,49 +430,49 @@ namespace search { std::vector current_pixels(stamp_ppi, 0.0); std::vector empty_pixels(1, NO_DATA); for (int t = 0; t < num_trajectories; ++t) { - // Copy the data into a single RawImage. - int offset = t * stamp_ppi; - for (unsigned p = 0; p < stamp_ppi; ++p) { - current_pixels[p] = stamp_data[offset + p]; - } - RawImage current_image = RawImage(stamp_width, stamp_width, current_pixels); - - if (params.do_filtering && filter_stamp(current_image, params)) { - results[t] = RawImage(1, 1, empty_pixels); - } else { - results[t] = RawImage(stamp_width, stamp_width, current_pixels); - } + // Copy the data into a single RawImage. + int offset = t * stamp_ppi; + for (unsigned p = 0; p < stamp_ppi; ++p) { + current_pixels[p] = stamp_data[offset + p]; + } + RawImage current_image = RawImage(stamp_width, stamp_width, current_pixels); + + if (params.do_filtering && filter_stamp(current_image, params)) { + results[t] = RawImage(1, 1, empty_pixels); + } else { + results[t] = RawImage(stamp_width, stamp_width, current_pixels); + } } return results; - } +} - std::vector StackSearch::create_stamps(Trajectory t, int radius, const std::vector& imgs, - bool interpolate) { +std::vector StackSearch::create_stamps(Trajectory t, int radius, const std::vector& imgs, + bool interpolate) { if (radius < 0) throw std::runtime_error("stamp radius must be at least 0"); std::vector stamps; for (int i = 0; i < imgs.size(); ++i) { - PixelPos pos = get_trajectory_position(t, i); - stamps.push_back(imgs[i]->create_stamp(pos.x, pos.y, radius, interpolate, false)); + PixelPos pos = get_trajectory_position(t, i); + stamps.push_back(imgs[i]->create_stamp(pos.x, pos.y, radius, interpolate, false)); } return stamps; - } +} - PixelPos StackSearch::get_trajectory_position(const Trajectory& t, int i) const { +PixelPos StackSearch::get_trajectory_position(const Trajectory& t, int i) const { float time = stack.get_zeroed_time(i); return {t.x + time * t.vx, t.y + time * t.vy}; - } +} - std::vector StackSearch::get_trajectory_positions(Trajectory& t) const { +std::vector StackSearch::get_trajectory_positions(Trajectory& t) const { std::vector results; int num_times = stack.img_count(); for (int i = 0; i < num_times; ++i) { - PixelPos pos = get_trajectory_position(t, i); - results.push_back(pos); + PixelPos pos = get_trajectory_position(t, i); + results.push_back(pos); } return results; - } +} - std::vector StackSearch::create_curves(Trajectory t, const std::vector& imgs) { +std::vector StackSearch::create_curves(Trajectory t, const std::vector& imgs) { /*Create a lightcurve from an image along a trajectory * * INPUT- @@ -490,18 +488,17 @@ namespace search { lightcurve.reserve(img_size); std::vector times = stack.build_zeroed_times(); for (int i = 0; i < img_size; ++i) { - /* Do not use get_pixel_interp(), because results from create_curves must - * be able to recover the same likelihoods as the ones reported by the - * gpu search.*/ - float pix_val = imgs[i].get_pixel(t.x + int(times[i] * t.vx + 0.5), - t.y + int(times[i] * t.vy + 0.5)); - if (pix_val == NO_DATA) pix_val = 0.0; - lightcurve.push_back(pix_val); + /* Do not use get_pixel_interp(), because results from create_curves must + * be able to recover the same likelihoods as the ones reported by the + * gpu search.*/ + float pix_val = imgs[i].get_pixel(t.x + int(times[i] * t.vx + 0.5), t.y + int(times[i] * t.vy + 0.5)); + if (pix_val == NO_DATA) pix_val = 0.0; + lightcurve.push_back(pix_val); } return lightcurve; - } +} - std::vector StackSearch::get_psi_curves(Trajectory& t) { +std::vector StackSearch::get_psi_curves(Trajectory& t) { /*Generate a psi lightcurve for further analysis * INPUT- * Trajectory& t - The trajectory along which to find the lightcurve @@ -510,9 +507,9 @@ namespace search { */ prepare_psi_phi(); return create_curves(t, psi_images); - } +} - std::vector StackSearch::get_phi_curves(Trajectory& t) { +std::vector StackSearch::get_phi_curves(Trajectory& t) { /*Generate a phi lightcurve for further analysis * INPUT- * Trajectory& t - The trajectory along which to find the lightcurve @@ -521,59 +518,59 @@ namespace search { */ prepare_psi_phi(); return create_curves(t, phi_images); - } +} - std::vector& StackSearch::get_psi_images() { return psi_images; } +std::vector& StackSearch::get_psi_images() { return psi_images; } - std::vector& StackSearch::get_phi_images() { return phi_images; } +std::vector& StackSearch::get_phi_images() { return phi_images; } - void StackSearch::sort_results() { +void StackSearch::sort_results() { __gnu_parallel::sort(results.begin(), results.end(), [](Trajectory a, Trajectory b) { return b.lh < a.lh; }); - } +} - void StackSearch::filter_results(int min_observations) { +void StackSearch::filter_results(int min_observations) { results.erase(std::remove_if(results.begin(), results.end(), std::bind([](Trajectory t, int cutoff) { return t.obs_count < cutoff; }, std::placeholders::_1, min_observations)), results.end()); - } +} - void StackSearch::filter_results_lh(float min_lh) { +void StackSearch::filter_results_lh(float min_lh) { results.erase(std::remove_if(results.begin(), results.end(), std::bind([](Trajectory t, float cutoff) { return t.lh < cutoff; }, std::placeholders::_1, min_lh)), results.end()); - } +} - std::vector StackSearch::get_results(int start, int count) { +std::vector StackSearch::get_results(int start, int count) { if (start + count >= results.size()) { - count = results.size() - start; + count = results.size() - start; } if (start < 0) throw std::runtime_error("start must be 0 or greater"); return std::vector(results.begin() + start, results.begin() + start + count); - } +} - // This function is used only for testing by injecting known result trajectories. - void StackSearch::set_results(const std::vector& new_results) { results = new_results; } +// This function is used only for testing by injecting known result trajectories. +void StackSearch::set_results(const std::vector& new_results) { results = new_results; } - void StackSearch::start_timer(const std::string& message) { +void StackSearch::start_timer(const std::string& message) { if (debug_info) { - std::cout << message << "... " << std::flush; - t_start = std::chrono::system_clock::now(); + std::cout << message << "... " << std::flush; + t_start = std::chrono::system_clock::now(); } - } +} - void StackSearch::end_timer() { +void StackSearch::end_timer() { if (debug_info) { - t_end = std::chrono::system_clock::now(); - t_delta = t_end - t_start; - std::cout << " Took " << t_delta.count() << " seconds.\n" << std::flush; + t_end = std::chrono::system_clock::now(); + t_delta = t_end - t_start; + std::cout << " Took " << t_delta.count() << " seconds.\n" << std::flush; } - } +} #ifdef Py_PYTHON_H - static void stack_search_bindings(py::module &m) { +static void stack_search_bindings(py::module& m) { using tj = search::Trajectory; using pf = search::PSF; using ri = search::RawImage; @@ -581,43 +578,48 @@ namespace search { using ks = search::StackSearch; py::class_(m, "StackSearch", pydocs::DOC_StackSearch) - .def(py::init()) - .def("save_psi_phi", &ks::save_psiphi, pydocs::DOC_StackSearch_save_psi_phi) - .def("search", &ks::search, pydocs::DOC_StackSearch_search) - .def("enable_gpu_sigmag_filter", &ks::enable_gpu_sigmag_filter, pydocs::DOC_StackSearch_enable_gpu_sigmag_filter) - .def("enable_gpu_encoding", &ks::enable_gpu_encoding, pydocs::DOC_StackSearch_enable_gpu_encoding) - .def("set_start_bounds_x", &ks::set_start_bounds_x, pydocs::DOC_StackSearch_set_start_bounds_x) - .def("set_start_bounds_y", &ks::set_start_bounds_y, pydocs::DOC_StackSearch_set_start_bounds_y) - .def("set_debug", &ks::set_debug, pydocs::DOC_StackSearch_set_debug) - .def("filter_min_obs", &ks::filter_results, pydocs::DOC_StackSearch_filter_min_obs) - .def("get_num_images", &ks::num_images, pydocs::DOC_StackSearch_get_num_images) - .def("get_image_width", &ks::get_image_width, pydocs::DOC_StackSearch_get_image_width) - .def("get_image_height", &ks::get_image_height, pydocs::DOC_StackSearch_get_image_height) - .def("get_image_npixels", &ks::get_image_npixels, pydocs::DOC_StackSearch_get_image_npixels) - .def("get_imagestack", &ks::get_imagestack, - py::return_value_policy::reference_internal, - pydocs::DOC_StackSearch_get_imagestack) - // Science Stamp Functions - .def("get_stamps", &ks::get_stamps, pydocs::DOC_StackSearch_get_stamps) - .def("get_median_stamp", &ks::get_median_stamp, pydocs::DOC_StackSearch_get_median_stamp) - .def("get_mean_stamp", &ks::get_mean_stamp, pydocs::DOC_StackSearch_get_mean_stamp) - .def("get_summed_stamp", &ks::get_summed_stamp, pydocs::DOC_StackSearch_get_summed_stamp) - .def("get_coadded_stamps", //wth is happening here - (std::vector(ks::*)(std::vector &, std::vector> &, - const search::StampParameters &, bool)) & - ks::get_coadded_stamps, pydocs::DOC_StackSearch_get_coadded_stamps) - // For testing - .def("filter_stamp", &ks::filter_stamp, pydocs::DOC_StackSearch_filter_stamp) - .def("get_trajectory_position", &ks::get_trajectory_position, pydocs::DOC_StackSearch_get_trajectory_position) - .def("get_trajectory_positions", &ks::get_trajectory_positions, pydocs::DOC_StackSearch_get_trajectory_positions) - .def("get_psi_curves", (std::vector(ks::*)(tj &)) & ks::get_psi_curves, pydocs::DOC_StackSearch_get_psi_curves) - .def("get_phi_curves", (std::vector(ks::*)(tj &)) & ks::get_phi_curves, pydocs::DOC_StackSearch_get_phi_curves) - .def("prepare_psi_phi", &ks::prepare_psi_phi, pydocs::DOC_StackSearch_prepare_psi_phi) - .def("get_psi_images", &ks::get_psi_images, pydocs::DOC_StackSearch_get_psi_images) - .def("get_phi_images", &ks::get_phi_images, pydocs::DOC_StackSearch_get_phi_images) - .def("get_results", &ks::get_results, pydocs::DOC_StackSearch_get_results) - .def("set_results", &ks::set_results, pydocs::DOC_StackSearch_set_results); - } + .def(py::init()) + .def("save_psi_phi", &ks::save_psiphi, pydocs::DOC_StackSearch_save_psi_phi) + .def("search", &ks::search, pydocs::DOC_StackSearch_search) + .def("enable_gpu_sigmag_filter", &ks::enable_gpu_sigmag_filter, + pydocs::DOC_StackSearch_enable_gpu_sigmag_filter) + .def("enable_gpu_encoding", &ks::enable_gpu_encoding, pydocs::DOC_StackSearch_enable_gpu_encoding) + .def("set_start_bounds_x", &ks::set_start_bounds_x, pydocs::DOC_StackSearch_set_start_bounds_x) + .def("set_start_bounds_y", &ks::set_start_bounds_y, pydocs::DOC_StackSearch_set_start_bounds_y) + .def("set_debug", &ks::set_debug, pydocs::DOC_StackSearch_set_debug) + .def("filter_min_obs", &ks::filter_results, pydocs::DOC_StackSearch_filter_min_obs) + .def("get_num_images", &ks::num_images, pydocs::DOC_StackSearch_get_num_images) + .def("get_image_width", &ks::get_image_width, pydocs::DOC_StackSearch_get_image_width) + .def("get_image_height", &ks::get_image_height, pydocs::DOC_StackSearch_get_image_height) + .def("get_image_npixels", &ks::get_image_npixels, pydocs::DOC_StackSearch_get_image_npixels) + .def("get_imagestack", &ks::get_imagestack, py::return_value_policy::reference_internal, + pydocs::DOC_StackSearch_get_imagestack) + // Science Stamp Functions + .def("get_stamps", &ks::get_stamps, pydocs::DOC_StackSearch_get_stamps) + .def("get_median_stamp", &ks::get_median_stamp, pydocs::DOC_StackSearch_get_median_stamp) + .def("get_mean_stamp", &ks::get_mean_stamp, pydocs::DOC_StackSearch_get_mean_stamp) + .def("get_summed_stamp", &ks::get_summed_stamp, pydocs::DOC_StackSearch_get_summed_stamp) + .def("get_coadded_stamps", // wth is happening here + (std::vector(ks::*)(std::vector&, std::vector>&, + const search::StampParameters&, bool)) & + ks::get_coadded_stamps, + pydocs::DOC_StackSearch_get_coadded_stamps) + // For testing + .def("filter_stamp", &ks::filter_stamp, pydocs::DOC_StackSearch_filter_stamp) + .def("get_trajectory_position", &ks::get_trajectory_position, + pydocs::DOC_StackSearch_get_trajectory_position) + .def("get_trajectory_positions", &ks::get_trajectory_positions, + pydocs::DOC_StackSearch_get_trajectory_positions) + .def("get_psi_curves", (std::vector(ks::*)(tj&)) & ks::get_psi_curves, + pydocs::DOC_StackSearch_get_psi_curves) + .def("get_phi_curves", (std::vector(ks::*)(tj&)) & ks::get_phi_curves, + pydocs::DOC_StackSearch_get_phi_curves) + .def("prepare_psi_phi", &ks::prepare_psi_phi, pydocs::DOC_StackSearch_prepare_psi_phi) + .def("get_psi_images", &ks::get_psi_images, pydocs::DOC_StackSearch_get_psi_images) + .def("get_phi_images", &ks::get_phi_images, pydocs::DOC_StackSearch_get_phi_images) + .def("get_results", &ks::get_results, pydocs::DOC_StackSearch_get_results) + .def("set_results", &ks::set_results, pydocs::DOC_StackSearch_set_results); +} #endif /* Py_PYTHON_H */ diff --git a/src/kbmod/search/stack_search.h b/src/kbmod/search/stack_search.h index e58ac5de..bdeec4aa 100644 --- a/src/kbmod/search/stack_search.h +++ b/src/kbmod/search/stack_search.h @@ -15,10 +15,9 @@ #include "psf.h" #include "pydocs/stack_search_docs.h" - namespace search { - class StackSearch { - public: +class StackSearch { +public: StackSearch(ImageStack& imstack); int num_images() const { return stack.img_count(); } @@ -90,7 +89,7 @@ namespace search { virtual ~StackSearch(){}; - protected: +protected: void save_images(const std::string& path); void sort_results(); std::vector create_curves(Trajectory t, const std::vector& imgs); @@ -110,8 +109,8 @@ namespace search { bool interpolate); // Creates list of trajectories to search. - void create_search_list(int angle_steps, int velocity_steps, float min_ang, float max_ang, - float min_vel, float max_vel); + void create_search_list(int angle_steps, int velocity_steps, float min_ang, float max_ang, float min_vel, + float max_vel); std::vector get_coadded_stamps_gpu(std::vector& t_array, std::vector >& use_index_vect, @@ -138,7 +137,7 @@ namespace search { // Parameters for the GPU search. SearchParameters params; - }; +}; } /* namespace search */