diff --git a/src/kbmod/fake_data_creator.py b/src/kbmod/fake_data_creator.py index 7b0dfc14..2665eb0d 100644 --- a/src/kbmod/fake_data_creator.py +++ b/src/kbmod/fake_data_creator.py @@ -110,8 +110,8 @@ def make_fake_ImageStack(self): for i in range(self.num_times): img = LayeredImage( ("%06i" % i), - self.width, self.height, + self.width, self.noise_level, self.noise_level**2, self.times[i], @@ -142,7 +142,7 @@ def insert_object(self, trj): # re-set the image. This last step needs to be done # explicitly because of how pybind handles references. current = self.stack.get_single_image(i) - add_fake_object(current, px, py, trj.flux, current.get_psf()) + add_fake_object(current, py, px, trj.flux, current.get_psf()) # Save the trajectory into the internal list. self.trajectories.append(trj) diff --git a/src/kbmod/search/bindings.cpp b/src/kbmod/search/bindings.cpp index 61f41a67..5faa4c80 100644 --- a/src/kbmod/search/bindings.cpp +++ b/src/kbmod/search/bindings.cpp @@ -32,6 +32,7 @@ PYBIND11_MODULE(search, m) { search::layered_image_bindings(m); search::image_stack_bindings(m); search::stack_search_bindings(m); + search::stamp_creator_bindings(m); search::trajectory_bindings(m); search::pixel_pos_bindings(m); search::image_moments_bindings(m); diff --git a/src/kbmod/search/raw_image.cpp b/src/kbmod/search/raw_image.cpp index eda51bb6..73937432 100644 --- a/src/kbmod/search/raw_image.cpp +++ b/src/kbmod/search/raw_image.cpp @@ -297,12 +297,12 @@ RawImage& RawImage::operator=(RawImage&& source) { for(int j = 0; j < height; ++j) { for(int i = 0; i < width; ++i){ if (bitmask(j, i) == -1){ - if (((j-1 > 0) && (bitmask(j-1, i) == itr-1)) || - ((i-1 > 0) && (bitmask(j, i-1) == itr-1)) || + if (((j-1 >= 0) && (bitmask(j-1, i) == itr-1)) || + ((i-1 >= 0) && (bitmask(j, i-1) == itr-1)) || ((j+1 < height) && (bitmask(j+1, i) == itr-1)) || ((i+1 < width) && (bitmask(j, i+1) == itr-1))){ bitmask(j, i) = itr; - } + } } } // for i } // for j @@ -351,10 +351,10 @@ RawImage& RawImage::operator=(RawImage&& source) { dist2 = new_dist2; } } - } - + } // for x + } // for y return result; -} + } // Find the basic image moments in order to test if stamps have a gaussian shape. @@ -581,8 +581,8 @@ RawImage& RawImage::operator=(RawImage&& source) { // and value based filtering. result(y, x) = 0.0; } - } - + } // for x + } // for y return RawImage(result); } diff --git a/src/kbmod/search/stack_search.cpp b/src/kbmod/search/stack_search.cpp index 70649a27..69749c82 100644 --- a/src/kbmod/search/stack_search.cpp +++ b/src/kbmod/search/stack_search.cpp @@ -2,12 +2,11 @@ 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); #endif - StackSearch::StackSearch(ImageStack& imstack) : stack(imstack) { debug_info = false; psi_phi_generated = false; @@ -216,73 +215,68 @@ extern "C" void deviceSearchFilter(int num_images, int width, int height, float* } } - void StackSearch::create_search_list(int angle_steps, int velocity_steps, float min_ang, float max_ang, + void StackSearch::create_search_list(int angle_steps, int velocity_steps, + float min_ang, float max_ang, float min_vel, float mavx) { DebugTimer timer = DebugTimer("Creating search candidate list", debug_info); + 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; + } - 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; - } - - 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; - } + 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; + } - 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]; - } + 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]; } - - timer.stop(); } + timer.stop(); + } + // I'm not a 100% sure what we're copying here and I can't see this being + // called somewhere + 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); - // I'm not a 100% sure what we're copying here and I can't see this being - // called somewhere - 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); - - int num_images = psi_imgs.size(); - assert(num_images > 0); - assert(phi_imgs.size() == num_images); + int num_images = psi_imgs.size(); + assert(num_images > 0); + assert(phi_imgs.size() == num_images); - 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); - } + 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); + } - psi_vect->clear(); - psi_vect->reserve(num_images * num_pixels); + psi_vect->clear(); + psi_vect->reserve(num_images * num_pixels); - for (int i = 0; i < num_images; ++i) { - const auto& psi_ref = psi_imgs[i].get_image().reshaped(); - const auto& phi_ref = phi_imgs[i].get_image().reshaped(); - for (unsigned p = 0; p < num_pixels; ++p) { - psi_vect->push_back(psi_ref[p]); - phi_vect->push_back(phi_ref[p]); - } + for (int i = 0; i < num_images; ++i) { + const auto& psi_ref = psi_imgs[i].get_image().reshaped(); + const auto& phi_ref = phi_imgs[i].get_image().reshaped(); + 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_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- @@ -292,110 +286,94 @@ extern "C" void deviceSearchFilter(int num_images, int width, int height, float* * Output- * std::vector lightcurve - The computed trajectory */ + int img_size = imgs.size(); std::vector lightcurve; 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.*/ - Point p({ - t.y + times[i] * t.vy + 0.5f, - t.x + times[i] * t.vx + 0.5f - }); - float pix_val = imgs[i].get_pixel(p.to_index()); - if (pix_val == NO_DATA) pix_val = 0.0; - lightcurve.push_back(pix_val); - } - return lightcurve; - } - - - 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 - * OUTPUT- - * std::vector - A vector of the lightcurve values - */ - prepare_psi_phi(); - return create_curves(t, psi_images); - } - - - 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 - * OUTPUT- - * std::vector - A vector of the lightcurve values - */ - prepare_psi_phi(); - return create_curves(t, phi_images); + /* 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.*/ + Point p({ + t.y + times[i] * t.vy + 0.5f, + t.x + times[i] * t.vx + 0.5f + }); + float pix_val = imgs[i].get_pixel(p.to_index()); + if (pix_val == NO_DATA) pix_val = 0.0; + lightcurve.push_back(pix_val); } + return lightcurve; + } - std::vector& StackSearch::get_psi_images() { return psi_images; } + 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 + * OUTPUT- + * std::vector - A vector of the lightcurve values + */ + prepare_psi_phi(); + return create_curves(t, psi_images); + } - std::vector& StackSearch::get_phi_images() { return phi_images; } + 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 + * OUTPUT- + * std::vector - A vector of the lightcurve values + */ + prepare_psi_phi(); + return create_curves(t, phi_images); + } - void StackSearch::sort_results() { - __gnu_parallel::sort(results.begin(), results.end(), - [](Trajectory a, Trajectory b) { return b.lh < a.lh; }); - } + std::vector& StackSearch::get_psi_images() { return psi_images; } - 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()); - } + std::vector& StackSearch::get_phi_images() { return phi_images; } - 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()); - } + void StackSearch::sort_results() { + __gnu_parallel::sort(results.begin(), results.end(), + [](Trajectory a, Trajectory b) { return b.lh < a.lh; }); + } - std::vector StackSearch::get_results(int start, int count) { - if (start + count >= results.size()) { - 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); - } + 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()); + } - // 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::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()); + } - void StackSearch::start_timer(const std::string& message) { - if (debug_info) { - std::cout << message << "... " << std::flush; - t_start = std::chrono::system_clock::now(); - } + std::vector StackSearch::get_results(int start, int count) { + if (start + count >= results.size()) { + 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); + } - 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; - } - } + // This function is used only for testing by injecting known result trajectories. + void StackSearch::set_results(const std::vector& new_results) { results = new_results; } #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; @@ -403,36 +381,35 @@ static void stack_search_bindings(py::module& m) { 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) - // For testings - .def("get_trajectory_position", &ks::get_trajectory_position, - pydocs::DOC_StackSearch_get_trajectory_position) - .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) + // For testings + .def("get_trajectory_position", &ks::get_trajectory_position, + pydocs::DOC_StackSearch_get_trajectory_position) + .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 */ - } /* namespace search */ +} /* namespace search */ diff --git a/src/kbmod/search/stack_search.h b/src/kbmod/search/stack_search.h index 6efbe67e..21916ad9 100644 --- a/src/kbmod/search/stack_search.h +++ b/src/kbmod/search/stack_search.h @@ -74,7 +74,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); @@ -102,7 +102,7 @@ namespace search { // Parameters for the GPU search. SearchParameters params; -}; + }; } /* namespace search */ diff --git a/src/kbmod/search/stamp_creator.cpp b/src/kbmod/search/stamp_creator.cpp index b012045b..b6fe335e 100644 --- a/src/kbmod/search/stamp_creator.cpp +++ b/src/kbmod/search/stamp_creator.cpp @@ -32,9 +32,9 @@ std::vector StampCreator::create_stamps(ImageStack& stack, const Traje if (use_all_stamps || use_index[i]) { // Calculate the trajectory position. float time = stack.get_zeroed_time(i); - PixelPos pos = {trj.x + time * trj.vx, trj.y + time * trj.vy}; + Point pos{trj.y + time * trj.vy, trj.x + time * trj.vx}; RawImage& img = stack.get_single_image(i).get_science(); - stamps.push_back(img.create_stamp(pos.x, pos.y, radius, interpolate, keep_no_data)); + stamps.push_back(img.create_stamp(pos, radius, interpolate, keep_no_data)); } } return stamps; diff --git a/src/kbmod/work_unit.py b/src/kbmod/work_unit.py index 27644c69..f261b1ab 100644 --- a/src/kbmod/work_unit.py +++ b/src/kbmod/work_unit.py @@ -153,11 +153,8 @@ def raw_image_to_hdu(img): hdu : `astropy.io.fits.hdu.image.ImageHDU` The image extension. """ - # Expensive copy. To be removed with RawImage refactor. - np_pixels = np.array(img.get_all_pixels(), dtype=np.single) - np_array = np_pixels.reshape((img.get_height(), img.get_width())) - hdu = fits.hdu.image.ImageHDU(np_array) - hdu.header["MJD"] = img.get_obstime() + hdu = fits.hdu.image.ImageHDU(img.image) + hdu.header["MJD"] = img.obstime return hdu @@ -176,8 +173,9 @@ def hdu_to_raw_image(hdu): """ img = None if isinstance(hdu, fits.hdu.image.ImageHDU): - # Expensive copy. To be removed with RawImage refactor. - img = RawImage(hdu.data) + # This will be a copy whenever dtype != np.single including when + # endianness doesn't match the native float. + img = RawImage(hdu.data.astype(np.single)) if "MJD" in hdu.header: - img.set_obstime(hdu.header["MJD"]) + imgobstime = hdu.header["MJD"] return img diff --git a/tests/test_fake_data_creator.py b/tests/test_fake_data_creator.py index 7c38bfce..6c0d8462 100644 --- a/tests/test_fake_data_creator.py +++ b/tests/test_fake_data_creator.py @@ -45,7 +45,7 @@ def test_insert_object(self): self.assertLess(py, 256) # Check that there is a bright spot at the predicted position. - pix_val = ds.stack.get_single_image(i).get_science().get_pixel(px, py) + pix_val = ds.stack.get_single_image(i).get_science().get_pixel(py, px) self.assertGreaterEqual(pix_val, 50.0) def test_save_and_clean(self): @@ -98,8 +98,8 @@ def test_save_work_unit(self): self.assertEqual(work2.im_stack.img_count(), num_images) for i in range(num_images): li = work2.im_stack.get_single_image(i) - self.assertEqual(li.get_width(), 15) - self.assertEqual(li.get_height(), 10) + self.assertEqual(li.get_width(), 10) + self.assertEqual(li.get_height(), 15) if __name__ == "__main__": diff --git a/tests/test_image_stack.py b/tests/test_image_stack.py index b9a22487..eb43d6dc 100644 --- a/tests/test_image_stack.py +++ b/tests/test_image_stack.py @@ -66,6 +66,7 @@ def test_apply_mask(self): self.assertTrue(sci.pixel_has_data(y, x)) self.im_stack.apply_mask_flags(1, []) + # Check that one pixel is masked in each time. for i in range(self.num_images): sci = self.im_stack.get_single_image(i).get_science() diff --git a/tests/test_masking.py b/tests/test_masking.py index e5c94149..b1331274 100644 --- a/tests/test_masking.py +++ b/tests/test_masking.py @@ -59,8 +59,8 @@ def test_threshold_masker(self): for i in range(self.img_count): img = self.stack.get_single_image(i) sci = img.get_science() - sci.set_pixel(2 + i, 8, 501.0) - sci.set_pixel(1 + i, 9, 499.0) + sci.set_pixel(8, 2 + i, 501.0) + sci.set_pixel(9, 1 + i, 499.0) # With a threshold of 500 one pixel per image should be masked. mask = ThresholdMask(500) @@ -70,9 +70,9 @@ def test_threshold_masker(self): for x in range(self.dim_x): for y in range(self.dim_y): if x == 2 + i and y == 8: - self.assertFalse(sci.pixel_has_data(x, y)) + self.assertFalse(sci.pixel_has_data(y, x)) else: - self.assertTrue(sci.pixel_has_data(x, y)) + self.assertTrue(sci.pixel_has_data(y, x)) def test_per_image_dictionary_mask(self): # Set each mask pixel in a row to one masking reason. @@ -80,7 +80,7 @@ def test_per_image_dictionary_mask(self): img = self.stack.get_single_image(i) msk = img.get_mask() for x in range(self.dim_x): - msk.set_pixel(x, 3, 2**x) + msk.set_pixel(3, x, 2**x) # Mask with two keys. mask = DictionaryMasker(self.mask_bits_dict, ["BAD", "EDGE"]) @@ -90,9 +90,9 @@ def test_per_image_dictionary_mask(self): for x in range(self.dim_x): for y in range(self.dim_y): if y == 3 and (x == 0 or x == 4): - self.assertFalse(sci.pixel_has_data(x, y)) + self.assertFalse(sci.pixel_has_data(y, x)) else: - self.assertTrue(sci.pixel_has_data(x, y)) + self.assertTrue(sci.pixel_has_data(y, x)) # Mask with all the default keys. mask = DictionaryMasker(self.mask_bits_dict, self.default_flag_keys) @@ -102,9 +102,9 @@ def test_per_image_dictionary_mask(self): for x in range(self.dim_x): for y in range(self.dim_y): if y == 3 and (x == 0 or x == 4 or x == 7 or x == 8 or x == 15): - self.assertFalse(sci.pixel_has_data(x, y)) + self.assertFalse(sci.pixel_has_data(y, x)) else: - self.assertTrue(sci.pixel_has_data(x, y)) + self.assertTrue(sci.pixel_has_data(y, x)) def test_mask_grow(self): # Mask one pixel per image. @@ -112,7 +112,7 @@ def test_mask_grow(self): img = self.stack.get_single_image(i) msk = img.get_mask() for x in range(self.dim_x): - msk.set_pixel(2 + i, 8, 1) + msk.set_pixel(8, 2 + i, 1) # Apply the bit vector based mask and check that one pixel per image is masked. self.stack = BitVectorMasker(1, []).apply_mask(self.stack) @@ -120,7 +120,7 @@ def test_mask_grow(self): sci = self.stack.get_single_image(i).get_science() for x in range(self.dim_x): for y in range(self.dim_y): - self.assertEqual(sci.pixel_has_data(x, y), x != (2 + i) or y != 8) + self.assertEqual(sci.pixel_has_data(y, x), x != (2 + i) or y != 8) # Grow the mask by two pixels and recheck. self.stack = GrowMask(2).apply_mask(self.stack) @@ -129,7 +129,7 @@ def test_mask_grow(self): for x in range(self.dim_x): for y in range(self.dim_y): dist = abs(2 + i - x) + abs(y - 8) - self.assertEqual(sci.pixel_has_data(x, y), dist > 2) + self.assertEqual(sci.pixel_has_data(y, x), dist > 2) def test_global_mask(self): # Set each mask pixel in a single row depending on the image number. @@ -151,7 +151,7 @@ def test_global_mask(self): sci = self.stack.get_single_image(i).get_science() for x in range(self.dim_x): for y in range(self.dim_y): - self.assertEqual(sci.pixel_has_data(x, y), x != 1 or y != 1) + self.assertEqual(sci.pixel_has_data(y, x), x != 1 or y != 1) class test_run_search_masking(unittest.TestCase): @@ -185,33 +185,33 @@ def test_apply_masks(self): bad_pixels = [] for i in range(self.img_count): - # Make the pixel (i, 1) in each image above threshold and (i, 10) + # Make the pixel (1, i) in each image above threshold and (10, i) # right below the threshold. img = self.stack.get_single_image(i) sci = img.get_science() - sci.set_pixel(i, 1, 1000.0) - bad_pixels.append((i, i, 1)) - sci.set_pixel(i, 10, 895.0) + sci.set_pixel(1, i, 1000.0) + bad_pixels.append((i, 1, i)) + sci.set_pixel(10, i, 895.0) - # Set the "BAD" key on pixel (15 + i, 2) in each image. + # Set the "BAD" key on pixel (2, 15 + i) in each image. msk = img.get_mask() - msk.set_pixel(15 + i, 2, 1) - bad_pixels.append((i, 15 + i, 2)) + msk.set_pixel(2, 15 + i, 1) + bad_pixels.append((i, 2, 15 + i)) - # Set the "INTRP" key on pixel (30 + i, 3) in each image. + # Set the "INTRP" key on pixel (3, 30 + i) in each image. # But we won't mark this as a pixel to filter. - msk.set_pixel(30 + i, 3, 4) + msk.set_pixel(3, 30 + i, 4) # Set the "EDGE" key on pixel (2 * i, 4) in every third image. if i % 3 == 1: - msk.set_pixel(2 * i, 4, 16) - bad_pixels.append((i, 2 * i, 4)) + msk.set_pixel(4, 2 * i, 16) + bad_pixels.append((i, 4, 2 * i)) - # Set the "CR" key on pixel (6, 5) in every other image. + # Set the "CR" key on pixel (5, 6) in every other image. # It will be bad in every image because of global masking. if i % 2 == 0: - msk.set_pixel(6, 5, 8) - bad_pixels.append((i, 6, 5)) + msk.set_pixel(5, 6, 8) + bad_pixels.append((i, 5, 6)) bad_set = set(bad_pixels) @@ -225,18 +225,18 @@ def test_apply_masks(self): # Test the the correct pixels have been masked. for i in range(self.img_count): sci = self.stack.get_single_image(i).get_science() - for x in range(self.dim_x): - for y in range(self.dim_y): + for y in range(self.dim_y): + for x in range(self.dim_x): # Extend is_bad to check for pixels marked bad in our list OR # pixels adjacent to a bad pixel (from growing the mask 1). is_bad = ( - (i, x, y) in bad_set - or (i, x + 1, y) in bad_set - or (i, x - 1, y) in bad_set - or (i, x, y + 1) in bad_set - or (i, x, y - 1) in bad_set + (i, y, x) in bad_set + or (i, y + 1, x) in bad_set + or (i, y - 1, x) in bad_set + or (i, y, x + 1) in bad_set + or (i, y, x - 1) in bad_set ) - self.assertEqual(sci.pixel_has_data(x, y), not is_bad) + self.assertEqual(sci.pixel_has_data(y, x), not is_bad) if __name__ == "__main__": diff --git a/tests/test_search.py b/tests/test_search.py index 769a7906..2da6e8f6 100644 --- a/tests/test_search.py +++ b/tests/test_search.py @@ -153,7 +153,6 @@ def test_results(self): results = self.search.get_results(0, 10) best = results[0] - breakpoint() self.assertAlmostEqual(best.x, self.start_x, delta=self.pixel_error) self.assertAlmostEqual(best.y, self.start_y, delta=self.pixel_error) self.assertAlmostEqual(best.vx / self.vxel, 1, delta=self.velocity_error) @@ -751,137 +750,137 @@ def test_coadd_cpu_use_inds(self): self.assertAlmostEqual(sum_0 / count_0, meanStamps[0].get_pixel(stamp_y, stamp_x), delta=1e-3) self.assertAlmostEqual(sum_1 / count_1, meanStamps[1].get_pixel(stamp_y, stamp_x), delta=1e-3) -# @unittest.skipIf(not HAS_GPU, "Skipping test (no GPU detected)") -# def test_coadd_gpu_use_inds(self): -# params = StampParameters() -# params.radius = 1 -# params.do_filtering = False -# params.stamp_type = StampType.STAMP_MEAN -# -# # Mark a few of the observations as "do not use" -# inds = [[True] * self.imCount, [True] * self.imCount] -# inds[0][5] = False -# inds[1][3] = False -# inds[1][6] = False -# inds[1][7] = False -# inds[1][11] = False -# -# # Compute the stacked science (summed and mean) from a single Trajectory. -# meanStamps = StampCreator.get_coadded_stamps( -# self.search.get_imagestack(), [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, False -# ) -# -# # Compute the true summed and mean pixels for all of the pixels in the stamp. -# times = self.stack.build_zeroed_times() -# for stamp_x in range(2 * params.radius + 1): -# for stamp_y in range(2 * params.radius + 1): -# x_offset = stamp_x - params.radius -# y_offset = stamp_y - params.radius -# -# sum_0 = 0.0 -# sum_1 = 0.0 -# count_0 = 0.0 -# count_1 = 0.0 -# for i in range(self.imCount): -# t = times[i] -# x = int(self.trj.x + self.trj.vx * t) + x_offset -# y = int(self.trj.y + self.trj.vy * t) + y_offset -# pixVal = self.imlist[i].get_science().get_pixel(y, x) -# -# if pixVal != KB_NO_DATA and inds[0][i] > 0: -# sum_0 += pixVal -# count_0 += 1.0 -# -# if pixVal != KB_NO_DATA and inds[1][i] > 0: -# sum_1 += pixVal -# count_1 += 1.0 -# -# # Check that we get the correct answers. -# self.assertAlmostEqual(count_0, 19.0) -# self.assertAlmostEqual(count_1, 16.0) -# self.assertAlmostEqual(sum_0 / count_0, meanStamps[0].get_pixel(stamp_y, stamp_x), delta=1e-3) -# self.assertAlmostEqual(sum_1 / count_1, meanStamps[1].get_pixel(stamp_y, stamp_x), delta=1e-3) -# -# def test_coadd_filter_cpu(self): -# # Create a second Trajectory that isn't any good. -# trj2 = Trajectory() -# trj2.x = 1 -# trj2.y = 1 -# trj2.vx = 0 -# trj2.vy = 0 -# -# # Create a third Trajectory that is close to good, but offset. -# trj3 = Trajectory() -# trj3.x = self.trj.x + 2 -# trj3.y = self.trj.y + 2 -# trj3.vx = self.trj.vx -# trj3.vy = self.trj.vy -# -# # Create a fourth Trajectory that is close enough -# trj4 = Trajectory() -# trj4.x = self.trj.x + 1 -# trj4.y = self.trj.y + 1 -# trj4.vx = self.trj.vx -# trj4.vy = self.trj.vy -# -# # Compute the stacked science from a single Trajectory. -# all_valid_vect = [(self.all_valid) for i in range(4)] -# meanStamps = StampCreator.get_coadded_stamps( -# self.search.get_imagestack(), [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, True -# ) -# -# # The first and last are unfiltered -# self.assertEqual(meanStamps[0].width, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[0].height, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[3].width, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[3].height, 2 * self.params.radius + 1) -# -# # The second and third are filtered. -# self.assertEqual(meanStamps[1].width, 1) -# self.assertEqual(meanStamps[1].height, 1) -# self.assertEqual(meanStamps[2].width, 1) -# self.assertEqual(meanStamps[2].height, 1) -# -# @unittest.skipIf(not HAS_GPU, "Skipping test (no GPU detected)") -# def test_coadd_filter_gpu(self): -# # Create a second Trajectory that isn't any good. -# trj2 = Trajectory() -# trj2.x = 1 -# trj2.y = 1 -# trj2.vx = 0 -# trj2.vy = 0 -# -# # Create a third Trajectory that is close to good, but offset. -# trj3 = Trajectory() -# trj3.x = self.trj.x + 2 -# trj3.y = self.trj.y + 2 -# trj3.vx = self.trj.vx -# trj3.vy = self.trj.vy -# -# # Create a fourth Trajectory that is close enough -# trj4 = Trajectory() -# trj4.x = self.trj.x + 1 -# trj4.y = self.trj.y + 1 -# trj4.vx = self.trj.vx -# trj4.vy = self.trj.vy -# -# # Compute the stacked science from a single Trajectory. -# all_valid_vect = [(self.all_valid) for i in range(4)] -# meanStamps = self.search.get_coadded_stamps( -# [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, True -# ) -# -# # The first and last are unfiltered -# self.assertEqual(meanStamps[0].width, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[0].height, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[3].width, 2 * self.params.radius + 1) -# self.assertEqual(meanStamps[3].height, 2 * self.params.radius + 1) -# -# # The second and third are filtered. -# self.assertEqual(meanStamps[1].width, 1) -# self.assertEqual(meanStamps[1].height, 1) -# self.assertEqual(meanStamps[2].width, 1) -# self.assertEqual(meanStamps[2].height, 1) + @unittest.skipIf(not HAS_GPU, "Skipping test (no GPU detected)") + def test_coadd_gpu_use_inds(self): + params = StampParameters() + params.radius = 1 + params.do_filtering = False + params.stamp_type = StampType.STAMP_MEAN + + # Mark a few of the observations as "do not use" + inds = [[True] * self.imCount, [True] * self.imCount] + inds[0][5] = False + inds[1][3] = False + inds[1][6] = False + inds[1][7] = False + inds[1][11] = False + + # Compute the stacked science (summed and mean) from a single Trajectory. + meanStamps = StampCreator.get_coadded_stamps( + self.search.get_imagestack(), [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, False + ) + + # Compute the true summed and mean pixels for all of the pixels in the stamp. + times = self.stack.build_zeroed_times() + for stamp_x in range(2 * params.radius + 1): + for stamp_y in range(2 * params.radius + 1): + x_offset = stamp_x - params.radius + y_offset = stamp_y - params.radius + + sum_0 = 0.0 + sum_1 = 0.0 + count_0 = 0.0 + count_1 = 0.0 + for i in range(self.imCount): + t = times[i] + x = int(self.trj.x + self.trj.vx * t) + x_offset + y = int(self.trj.y + self.trj.vy * t) + y_offset + pixVal = self.imlist[i].get_science().get_pixel(y, x) + + if pixVal != KB_NO_DATA and inds[0][i] > 0: + sum_0 += pixVal + count_0 += 1.0 + + if pixVal != KB_NO_DATA and inds[1][i] > 0: + sum_1 += pixVal + count_1 += 1.0 + + # Check that we get the correct answers. + self.assertAlmostEqual(count_0, 19.0) + self.assertAlmostEqual(count_1, 16.0) + self.assertAlmostEqual(sum_0 / count_0, meanStamps[0].get_pixel(stamp_y, stamp_x), delta=1e-3) + self.assertAlmostEqual(sum_1 / count_1, meanStamps[1].get_pixel(stamp_y, stamp_x), delta=1e-3) + + def test_coadd_filter_cpu(self): + # Create a second Trajectory that isn't any good. + trj2 = Trajectory() + trj2.x = 1 + trj2.y = 1 + trj2.vx = 0 + trj2.vy = 0 + + # Create a third Trajectory that is close to good, but offset. + trj3 = Trajectory() + trj3.x = self.trj.x + 2 + trj3.y = self.trj.y + 2 + trj3.vx = self.trj.vx + trj3.vy = self.trj.vy + + # Create a fourth Trajectory that is close enough + trj4 = Trajectory() + trj4.x = self.trj.x + 1 + trj4.y = self.trj.y + 1 + trj4.vx = self.trj.vx + trj4.vy = self.trj.vy + + # Compute the stacked science from a single Trajectory. + all_valid_vect = [(self.all_valid) for i in range(4)] + meanStamps = StampCreator.get_coadded_stamps( + self.search.get_imagestack(), [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, True + ) + + # The first and last are unfiltered + self.assertEqual(meanStamps[0].width, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[0].height, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[3].width, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[3].height, 2 * self.params.radius + 1) + + # The second and third are filtered. + self.assertEqual(meanStamps[1].width, 1) + self.assertEqual(meanStamps[1].height, 1) + self.assertEqual(meanStamps[2].width, 1) + self.assertEqual(meanStamps[2].height, 1) + + @unittest.skipIf(not HAS_GPU, "Skipping test (no GPU detected)") + def test_coadd_filter_gpu(self): + # Create a second Trajectory that isn't any good. + trj2 = Trajectory() + trj2.x = 1 + trj2.y = 1 + trj2.vx = 0 + trj2.vy = 0 + + # Create a third Trajectory that is close to good, but offset. + trj3 = Trajectory() + trj3.x = self.trj.x + 2 + trj3.y = self.trj.y + 2 + trj3.vx = self.trj.vx + trj3.vy = self.trj.vy + + # Create a fourth Trajectory that is close enough + trj4 = Trajectory() + trj4.x = self.trj.x + 1 + trj4.y = self.trj.y + 1 + trj4.vx = self.trj.vx + trj4.vy = self.trj.vy + + # Compute the stacked science from a single Trajectory. + all_valid_vect = [(self.all_valid) for i in range(4)] + meanStamps = self.search.get_coadded_stamps( + [self.trj, trj2, trj3, trj4], all_valid_vect, self.params, True + ) + + # The first and last are unfiltered + self.assertEqual(meanStamps[0].width, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[0].height, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[3].width, 2 * self.params.radius + 1) + self.assertEqual(meanStamps[3].height, 2 * self.params.radius + 1) + + # The second and third are filtered. + self.assertEqual(meanStamps[1].width, 1) + self.assertEqual(meanStamps[1].height, 1) + self.assertEqual(meanStamps[2].width, 1) + self.assertEqual(meanStamps[2].height, 1) if __name__ == "__main__": diff --git a/tests/test_stamp_parity.py b/tests/test_stamp_parity.py index 21676af5..549daafe 100644 --- a/tests/test_stamp_parity.py +++ b/tests/test_stamp_parity.py @@ -58,12 +58,12 @@ def setUp(self): for i in range(self.imCount): time = i / self.imCount im = LayeredImage( - str(i), self.dim_x, self.dim_y, self.noise_level, self.variance, time, self.p, i + str(i), self.dim_y, self.dim_x, self.noise_level, self.variance, time, self.p, i ) add_fake_object( im, - self.start_x + time * self.vxel + 0.5, self.start_y + time * self.vyel + 0.5, + self.start_x + time * self.vxel + 0.5, self.object_flux, self.p, ) @@ -71,7 +71,7 @@ def setUp(self): # Mask a pixel in half the images. if i % 2 == 0: mask = im.get_mask() - mask.set_pixel(self.masked_x, self.masked_y, 1) + mask.set_pixel(self.masked_y, self.masked_x, 1) im.apply_mask_flags(1, []) self.imlist.append(im) diff --git a/tests/test_work_unit.py b/tests/test_work_unit.py index ca6e6e8e..2b6ad6d0 100644 --- a/tests/test_work_unit.py +++ b/tests/test_work_unit.py @@ -20,8 +20,8 @@ def setUp(self): self.p[i] = kb.PSF(5.0 / float(2 * i + 1)) self.images[i] = kb.LayeredImage( ("layered_test_%i" % i), - self.width, self.height, + self.width, 2.0, # noise_level 4.0, # variance 2.0 * i + 1.0, # time @@ -79,9 +79,9 @@ def test_save_and_load_fits(self): for y in range(self.height): for x in range(self.width): - self.assertAlmostEqual(sci1.get_pixel(x, y), sci2.get_pixel(x, y)) - self.assertAlmostEqual(var1.get_pixel(x, y), var2.get_pixel(x, y)) - self.assertAlmostEqual(msk1.get_pixel(x, y), msk2.get_pixel(x, y)) + self.assertAlmostEqual(sci1.get_pixel(y, x), sci2.get_pixel(y, x)) + self.assertAlmostEqual(var1.get_pixel(y, x), var2.get_pixel(y, x)) + self.assertAlmostEqual(msk1.get_pixel(y, x), msk2.get_pixel(y, x)) # Check the PSF layer matches. p1 = self.p[i] @@ -90,7 +90,7 @@ def test_save_and_load_fits(self): for y in range(p1.get_dim()): for x in range(p1.get_dim()): - self.assertAlmostEqual(p1.get_value(x, y), p2.get_value(x, y)) + self.assertAlmostEqual(p1.get_value(y, x), p2.get_value(y, x)) # Check that we read in the configuration values correctly. self.assertEqual(work2.config["im_filepath"], "Here")