37 #include <boost/any.hpp>
38 #include <boost/random.hpp>
40 #include <CCfits/CCfits>
67 namespace po = boost::program_options;
68 namespace fs = boost::filesystem;
91 DummyWCS(
int image_width,
int image_height,
double rotation,
double scale,
double shift_x,
double shift_y)
92 : m_image_width(image_width), m_image_height(image_height),
93 m_rotation(rotation), m_scale(1.0/scale), m_shift_x(shift_x), m_shift_y(shift_y) {}
105 auto c =
cos(m_rotation);
106 auto s =
sin(m_rotation);
109 {
"CTYPE1",
"'RA---TAN'"},
110 {
"CTYPE2",
"'DEC--TAN'"},
113 {
"RADSYS",
"'ICRS'"},
121 {
"CRPIX1",
std::to_string(m_image_width / 2.0 + 1.5 + m_shift_x)},
122 {
"CRPIX2",
std::to_string(m_image_height / 2.0 + 1.5 + m_shift_y)},
133 double m_rotation = 0.0;
134 double m_scale = 1.0;
135 double m_shift_x = 0.0;
136 double m_shift_y = 0.0;
139 template <
typename ImageType>
169 po::options_description config_options {
"TestImage options" };
172 config_options.add_options()
173 (
"output", po::value<string>()->required(),
"filename to save the created test image")
174 (
"output-weight", po::value<string>()->default_value(
""),
"filename to save the created weight map image")
175 (
"size", po::value<double>()->default_value(512.0),
"image size")
176 (
"bg-level", po::value<double>()->default_value(0.0),
"background level")
177 (
"bg-sigma", po::value<double>()->default_value(20.0),
"standard deviation of background gaussian noise")
178 (
"gain", po::value<double>()->default_value(0.0),
"gain in e-/adu, 0 for infinite gain")
179 (
"saturation", po::value<double>()->default_value(0.0),
"image saturation level, 0 for no saturation")
180 (
"psf-file", po::value<string>()->default_value(
""),
"Psf file for convolution")
181 (
"psf-fwhm", po::value<double>()->default_value(5.0),
182 "Full width half maximum for generated gaussian psf (used when no psf file is provided)")
183 (
"psf-scale", po::value<double>()->default_value(0.2),
"Pixel scale for generated gaussian psf")
184 (
"disable-psf", po::bool_switch(),
"Disable psf convolution")
185 (
"random-sources", po::value<int>()->default_value(0),
"Nb of random sources to add")
186 (
"source-list", po::value<string>()->default_value(
""),
"Use sources from file")
187 (
"source-catalog", po::value<string>()->default_value(
""),
"Use sources from file (skymaker format)")
188 (
"zero-point", po::value<double>()->default_value(0.0),
"Zero point for magnitudes in catalog")
189 (
"exposure-time", po::value<double>()->default_value(300.),
"Exposure time for objects in catalog")
190 (
"save-sources", po::value<string>()->default_value(
""),
"Filename to save final list of sources")
191 (
"bad-pixels", po::value<double>()->default_value(0.0),
"Probability for a pixel to be a bad pixel")
192 (
"bad-columns", po::value<double>()->default_value(0.0),
"Probability for a column of pixels to be bad")
193 (
"rotation", po::value<double>()->default_value(0.0),
"Rotate sources around middle point")
194 (
"scale", po::value<double>()->default_value(1.0),
"Scale factor")
195 (
"shift-x", po::value<double>()->default_value(0.0),
"Shift X")
196 (
"shift-y", po::value<double>()->default_value(0.0),
"Shift Y")
197 (
"model-size", po::value<double>()->default_value(0.0),
"Model size for the rasterization of sources, defaults to size")
198 (
"max-tile-memory", po::value<int>()->default_value(512),
"Maximum memory used for image tiles cache in megabytes")
199 (
"tile-size", po::value<int>()->default_value(256),
"Image tiles size in pixels")
200 (
"copy-coordinate-system", po::value<string>()->default_value(
""),
"Copy the coordinate system from another FITS file")
203 return config_options;
207 auto x_param = std::make_shared<ManualParameter>(source.
x);
208 auto y_param = std::make_shared<ManualParameter>(source.
y);
214 auto xs = std::make_shared<ManualParameter>(1);
215 auto ys = std::make_shared<ManualParameter>(source.
exp_aspect);
216 auto rot = std::make_shared<ManualParameter>(source.
exp_rot);
217 auto exp_n = std::make_shared<ManualParameter>(1);
219 auto exp_k = std::make_shared<ManualParameter>(1.7 / source.
exp_rad);
220 auto exp_i0 = std::make_shared<ManualParameter>(
225 auto exp = Euclid::make_unique<SersicModelComponent>(Euclid::make_unique<OldSharp>(), exp_i0, exp_n, exp_k);
226 component_list.clear();
229 std::move(component_list), xs, ys, rot, size, size, x_param, y_param, jacobian));
234 auto xs = std::make_shared<ManualParameter>(1);
235 auto ys = std::make_shared<ManualParameter>(source.
dev_aspect);
236 auto rot = std::make_shared<ManualParameter>(source.
dev_rot);
237 auto dev_n = std::make_shared<ManualParameter>(4);
239 auto dev_k = std::make_shared<ManualParameter>(
pow(3459.0 / source.
dev_rad, .25));
240 auto dev_i0 = std::make_shared<ManualParameter>(
244 auto exp = Euclid::make_unique<SersicModelComponent>(Euclid::make_unique<OldSharp>(), dev_i0, dev_n, dev_k);
245 component_list.clear();
248 std::move(component_list), xs, ys, rot, size, size, x_param, y_param, jacobian));
252 auto flux_param = std::make_shared<ManualParameter>(source.
point_flux);
253 point_models.
emplace_back(x_param, y_param, flux_param);
258 auto x_param = std::make_shared<ManualParameter>(
x);
259 auto y_param = std::make_shared<ManualParameter>(
y);
260 auto flux_param = std::make_shared<ManualParameter>(flux);
262 point_models.
emplace_back(x_param, y_param, flux_param);
268 boost::random::normal_distribution<> bg_noise_dist(background_level, background_sigma);
269 for (
int y=0;
y < image->getHeight();
y++) {
270 for (
int x=0;
x < image->getWidth();
x++) {
272 image->setValue(
x,
y, accessor.
getValue(
x,
y) + bg_noise_dist(m_rng));
281 for (
int y=0;
y < image->getHeight();
y++) {
282 for (
int x=0;
x < image->getWidth();
x++) {
284 if (pixel_value > 0.) {
285 image->setValue(
x,
y, boost::random::poisson_distribution<>(pixel_value * gain)(m_rng) / gain);
294 if (saturation_level > 0.0) {
295 for (
int y=0;
y < image->getHeight();
y++) {
296 for (
int x=0;
x < image->getWidth();
x++) {
305 for (
int y=0;
y < weight_map->getHeight();
y++) {
306 for (
int x=0;
x < weight_map->getWidth();
x++) {
307 if (boost::random::uniform_01<double>()(m_rng) < probability) {
308 weight_map->setValue(
x,
y, 0);
317 for (
int x=0;
x < weight_map->getWidth();
x++) {
318 if (boost::random::uniform_01<double>()(m_rng) < probability) {
319 for (
int y=0;
y < weight_map->getHeight();
y++) {
320 weight_map->setValue(
x,
y, 0);
328 double x_min,
double y_min,
double x_max,
double y_max) {
333 boost::random::uniform_real_distribution<> random_x(x_min, x_max);
334 boost::random::uniform_real_distribution<> random_y(y_min, y_max);
338 for (
int i = 0; i < number_of_sources; i++) {
341 random_x(m_rng), random_y(m_rng),
343 gal_exp_random_i0(m_rng),
344 boost::random::uniform_real_distribution<>(.5, 6)(m_rng),
345 boost::random::uniform_real_distribution<>(.2, .8)(m_rng),
346 boost::random::uniform_real_distribution<>(-M_PI/2, M_PI/2)(m_rng),
367 while (file.
good()) {
372 if (line.
size() == 0) {
379 linestream >> source.
x >> source.
y
384 source.
exp_rot *= -M_PI / 180.0;
385 source.
dev_rot *= -M_PI / 180.0;
397 for (
const auto& source : sources) {
398 file << source.x <<
" " << source.y <<
" "
399 << source.exp_flux <<
" " << source.exp_rad <<
" " << source.exp_aspect <<
" " << (source.exp_rot * -180.0 / M_PI) <<
" "
400 << source.dev_flux <<
" " << source.dev_rad <<
" " << source.dev_aspect <<
" " << (source.dev_rot * -180.0 / M_PI) <<
" "
401 << source.point_flux <<
"\n";
406 auto center = image_size / 2.0;
407 auto c =
cos(rot_angle);
408 auto s =
sin(rot_angle);
409 for (
auto& source : sources) {
412 double x = (source.x * c - source.y *
s) * scale + shift_x;
413 double y = (source.x *
s + source.y * c) * scale + shift_y;
414 source.x =
x + center;
415 source.y =
y + center;
417 source.exp_rot -= rot_angle;
418 source.dev_rot -= rot_angle;
430 while (file.
good()) {
435 if (line.
size() == 0) {
441 linestream >> source_type;
446 if (source_type == 200) {
447 double magnitude, bt, exp_angle_degree, dev_angle_degree;
449 linestream >> source.
x >> source.
y >> magnitude
455 source.
dev_rot = dev_angle_degree * M_PI / 180.0;
456 source.
exp_rot = exp_angle_degree * M_PI / 180.0;
458 auto total_flux =
pow(10, (magnitude - m_zero_point) / -2.5);
459 source.
exp_flux = total_flux * (1.0 - bt);
466 }
else if (source_type == 100) {
468 linestream >> source.
x >> source.
y >> magnitude;
469 source.
point_flux =
pow(10, (magnitude - m_zero_point) / -2.5) * m_exp_time;
488 auto max_tile_memory = args[
"max-tile-memory"].as<
int>();
489 auto tile_size = args[
"tile-size"].as<
int>();
492 auto image_size = args[
"size"].as<
double>();
493 auto rot_angle = args[
"rotation"].as<
double>() / 180.0 * M_PI;
494 auto scale = args[
"scale"].as<
double>();
495 auto shift_x = args[
"shift-x"].as<
double>();
496 auto shift_y = args[
"shift-y"].as<
double>();
497 auto model_size = args[
"model-size"].as<
double>();
498 if (model_size <= 0.) {
499 model_size = image_size;
502 m_zero_point = args[
"zero-point"].as<
double>();
503 m_exp_time = args[
"exposure-time"].as<
double>();
511 auto psf_filename = args[
"psf-file"].as<
std::string>();
512 if (psf_filename !=
"") {
513 logger.info() <<
"Loading psf file: " << psf_filename;
520 auto copy_coordinate_system = args[
"copy-coordinate-system"].as<
std::string>();
521 if (copy_coordinate_system !=
"") {
522 coordinate_system = std::make_shared<WCS>(copy_coordinate_system);
524 coordinate_system = std::make_shared<DummyWCS>(image_size, image_size, rot_angle, scale, shift_x, shift_y);
527 auto raster_model_size = model_size / vpsf->getPixelSampling() +
std::max(vpsf->getWidth(), vpsf->getHeight());
529 logger.fatal() <<
"The expected required memory for model rasterization exceeds the maximum size for an integer";
530 logger.fatal() <<
"Please, either reduce the model size, the image size, or increase the PSF pixel scale";
536 const auto& vpsf_components = vpsf->getComponents();
538 for (
auto i = 0u; i < psf_vals.
size(); ++i) {
539 if (vpsf_components[i] ==
"X_IMAGE" || vpsf_components[i] ==
"Y_IMAGE") {
540 psf_vals[i] = image_size / 2 - 1;
548 auto p = vpsf->getPsf(psf_vals);
549 auto psf_sum =
std::accumulate(p->getData().begin(), p->getData().end(), 0.);
551 auto psf = std::make_shared<ImagePsf>(vpsf->getPixelSampling(), p);
556 if (sources_filename !=
"") {
557 logger.info() <<
"Loading sources from < " << sources_filename <<
" >...";
558 auto loaded_sources = loadSourcesFromFile(sources_filename);
559 sources.
insert(sources.
end(), loaded_sources.begin(), loaded_sources.end());
563 if (sources_cat_filename !=
"") {
564 logger.info() <<
"Loading source catalog from < " << sources_cat_filename <<
" >...";
565 auto loaded_sources = loadSourcesFromCatalog(sources_cat_filename);
566 sources.
insert(sources.
end(), loaded_sources.begin(), loaded_sources.end());
569 int nb_of_random_sources = args[
"random-sources"].as<
int>();
570 if (nb_of_random_sources > 0) {
571 logger.info() <<
"Adding " << nb_of_random_sources <<
" random source" << (nb_of_random_sources>1 ?
"s" :
"") <<
"...";
572 auto random_sources = generateRandomSources(nb_of_random_sources, image_size * .1, image_size * .1, image_size * .9, image_size * .9);
573 sources.
insert(sources.
end(), random_sources.begin(), random_sources.end());
577 if (sources_save_filename !=
"") {
578 logger.info() <<
"Saving sources to < " << sources_save_filename <<
" >...";
579 saveSources(sources, sources_save_filename);
582 logger.info(
"Transforming sources...");
583 transformSources(sources, image_size, rot_angle, scale, shift_x, shift_y);
585 logger.info(
"Creating source models...");
586 for (
const auto& source : sources) {
587 addSource(point_models, extended_models, model_size, source,
std::make_tuple(scale, 0, 0, scale));
590 logger.info(
"Rendering...");
593 auto target_image_source = std::make_shared<FitsImageSource>(
filename, image_size, image_size,
594 ImageTile::ImageType::FloatImage, coordinate_system);
597 if (args[
"disable-psf"].as<bool>()) {
605 frame_model.rasterToImage(target_image);
615 frame_model.rasterToImage(target_image);
619 logger.info(
"Adding noise...");
621 addPoissonNoise(target_image, args[
"gain"].as<double>());
622 addBackgroundNoise(target_image, args[
"bg-level"].as<double>(), args[
"bg-sigma"].as<double>());
624 logger.info(
"Adding saturation...");
626 auto saturation_level = args[
"saturation"].as<
double>();
627 saturate(target_image, saturation_level);
629 logger.info(
"Creating weight map...");
631 auto weight_filename = args[
"output-weight"].as<
std::string>();
632 if (weight_filename !=
"") {
633 auto weight_map_source = std::make_shared<FitsImageSource>(weight_filename, image_size, image_size, ImageTile::ImageType::FloatImage);
635 for (
int y = 0;
y < image_size; ++
y) {
636 for (
int x = 0;
x < image_size; ++
x) {
637 weight_map->setValue(
x,
y, 1);
641 logger.info(
"Adding bad pixels...");
643 addBadPixels(weight_map, args[
"bad-pixels"].as<double>());
644 addBadColumns(weight_map, args[
"bad-columns"].as<double>());
647 for (
int y = 0;
y < target_image->getHeight();
y++) {
648 for (
int x = 0;
x < target_image->getWidth();
x++) {
650 target_image->setValue(
x,
y, saturation_level);
656 logger.info(
"All done ^__^");
661 boost::random::mt19937 m_rng { (
unsigned int)
time(NULL) } ;
662 double m_zero_point = 0.0, m_exp_time = 300.;
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
const double gal_dev_min_i0
const double gal_exp_max_i0
const double gal_exp_min_i0
const double gal_dev_max_i0
double getPixelScale() const
std::shared_ptr< VectorImage< SourceXtractor::SeFloat > > getScaledKernel(double) const
void convolve(ImageType &) const
std::size_t getSize() const
ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
DummyWCS(int image_width, int image_height, double rotation, double scale, double shift_x, double shift_y)
std::map< std::string, std::string > getFitsHeaders() const override
static Logging getLogger(const std::string &name="")
Elements::ExitCode mainMethod(std::map< std::string, po::variable_value > &args) override
void saveSources(const std::vector< TestImageSource > &sources, const std::string &filename)
void addPoissonNoise(std::shared_ptr< WriteableImage< SeFloat >> image, double gain)
po::options_description defineSpecificProgramOptions() override
void addSource(std::vector< PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< WriteableInterfaceTypePtr >>> &extended_models, double size, const TestImageSource &source, std::tuple< double, double, double, double > jacobian)
void addPointSource(std::vector< PointModel > &point_models, double x, double y, double flux)
void addBadPixels(std::shared_ptr< WriteableImage< SeFloat >> weight_map, double probability)
std::vector< TestImageSource > loadSourcesFromFile(const std::string &filename)
void addBackgroundNoise(std::shared_ptr< WriteableImage< SeFloat >> image, double background_level, double background_sigma)
std::vector< TestImageSource > loadSourcesFromCatalog(const std::string &filename)
void saturate(std::shared_ptr< WriteableImage< SeFloat >> image, double saturation_level)
std::vector< TestImageSource > generateRandomSources(int number_of_sources, double x_min, double y_min, double x_max, double y_max)
void transformSources(std::vector< TestImageSource > &sources, int image_size, double rot_angle, double scale, double shift_x, double shift_y)
void addBadColumns(std::shared_ptr< WriteableImage< SeFloat >> weight_map, double probability)
T emplace_back(T... args)
#define MAIN_FOR(ELEMENTS_PROGRAM_NAME)
std::unique_ptr< T > make_unique(Args &&... args)