37 #include <boost/any.hpp>
38 #include <boost/random.hpp>
40 #include <CCfits/CCfits>
68 namespace po = boost::program_options;
69 namespace fs = boost::filesystem;
84 double exp_flux, exp_rad, exp_aspect,
exp_rot;
85 double dev_flux, dev_rad, dev_aspect,
dev_rot;
92 DummyWCS(
int image_width,
int image_height,
double rotation,
double scale,
double shift_x,
double shift_y)
93 : m_image_width(image_width), m_image_height(image_height),
94 m_rotation(rotation), m_scale(1.0/scale), m_shift_x(shift_x), m_shift_y(shift_y) {}
106 auto c =
cos(m_rotation);
107 auto s =
sin(m_rotation);
110 {
"CTYPE1",
"'RA---TAN'"},
111 {
"CTYPE2",
"'DEC--TAN'"},
114 {
"RADSYS",
"'ICRS'"},
122 {
"CRPIX1",
std::to_string(m_image_width / 2.0 + 1.5 + m_shift_x)},
123 {
"CRPIX2",
std::to_string(m_image_height / 2.0 + 1.5 + m_shift_y)},
134 double m_rotation = 0.0;
135 double m_scale = 1.0;
136 double m_shift_x = 0.0;
137 double m_shift_y = 0.0;
140 template <
typename ImageType>
170 po::options_description config_options {
"TestImage options" };
173 config_options.add_options()
174 (
"output", po::value<string>()->required(),
"filename to save the created test image")
175 (
"output-weight", po::value<string>()->default_value(
""),
"filename to save the created weight map image")
176 (
"size", po::value<double>()->default_value(512.0),
"image size")
177 (
"bg-level", po::value<double>()->default_value(0.0),
"background level")
178 (
"bg-sigma", po::value<double>()->default_value(20.0),
"standard deviation of background gaussian noise")
179 (
"gain", po::value<double>()->default_value(0.0),
"gain in e-/adu, 0 for infinite gain")
180 (
"saturation", po::value<double>()->default_value(0.0),
"image saturation level, 0 for no saturation")
181 (
"psf-file", po::value<string>()->default_value(
""),
"Psf file for convolution")
182 (
"psf-fwhm", po::value<double>()->default_value(5.0),
183 "Full width half maximum for generated gaussian psf (used when no psf file is provided)")
184 (
"psf-scale", po::value<double>()->default_value(0.2),
"Pixel scale for generated gaussian psf")
185 (
"disable-psf", po::bool_switch(),
"Disable psf convolution")
186 (
"random-sources", po::value<int>()->default_value(0),
"Nb of random sources to add")
187 (
"source-list", po::value<string>()->default_value(
""),
"Use sources from file")
188 (
"source-catalog", po::value<string>()->default_value(
""),
"Use sources from file (skymaker format)")
189 (
"zero-point", po::value<double>()->default_value(0.0),
"Zero point for magnitudes in catalog")
190 (
"exposure-time", po::value<double>()->default_value(300.),
"Exposure time for objects in catalog")
191 (
"save-sources", po::value<string>()->default_value(
""),
"Filename to save final list of sources")
192 (
"bad-pixels", po::value<double>()->default_value(0.0),
"Probability for a pixel to be a bad pixel")
193 (
"bad-columns", po::value<double>()->default_value(0.0),
"Probability for a column of pixels to be bad")
194 (
"rotation", po::value<double>()->default_value(0.0),
"Rotate sources around middle point")
195 (
"scale", po::value<double>()->default_value(1.0),
"Scale factor")
196 (
"shift-x", po::value<double>()->default_value(0.0),
"Shift X")
197 (
"shift-y", po::value<double>()->default_value(0.0),
"Shift Y")
198 (
"model-size", po::value<double>()->default_value(0.0),
"Model size for the rasterization of sources, defaults to size")
199 (
"max-tile-memory", po::value<int>()->default_value(512),
"Maximum memory used for image tiles cache in megabytes")
200 (
"tile-size", po::value<int>()->default_value(256),
"Image tiles size in pixels")
201 (
"copy-coordinate-system", po::value<string>()->default_value(
""),
"Copy the coordinate system from another FITS file")
204 return config_options;
208 auto x_param = std::make_shared<ManualParameter>(source.
x);
209 auto y_param = std::make_shared<ManualParameter>(source.
y);
215 auto xs = std::make_shared<ManualParameter>(1);
216 auto ys = std::make_shared<ManualParameter>(source.
exp_aspect);
217 auto rot = std::make_shared<ManualParameter>(source.
exp_rot);
218 auto exp_n = std::make_shared<ManualParameter>(1);
220 auto exp_k = std::make_shared<ManualParameter>(1.7 / source.
exp_rad);
221 auto exp_i0 = std::make_shared<ManualParameter>(
226 auto exp = make_unique<SersicModelComponent>(make_unique<OldSharp>(), exp_i0, exp_n, exp_k);
227 component_list.clear();
230 std::move(component_list), xs, ys, rot, size, size, x_param, y_param, jacobian));
235 auto xs = std::make_shared<ManualParameter>(1);
236 auto ys = std::make_shared<ManualParameter>(source.
dev_aspect);
237 auto rot = std::make_shared<ManualParameter>(source.
dev_rot);
238 auto dev_n = std::make_shared<ManualParameter>(4);
240 auto dev_k = std::make_shared<ManualParameter>(
pow(3459.0 / source.
dev_rad, .25));
241 auto dev_i0 = std::make_shared<ManualParameter>(
245 auto exp = make_unique<SersicModelComponent>(make_unique<OldSharp>(), dev_i0, dev_n, dev_k);
246 component_list.clear();
249 std::move(component_list), xs, ys, rot, size, size, x_param, y_param, jacobian));
253 auto flux_param = std::make_shared<ManualParameter>(source.
point_flux);
254 point_models.
emplace_back(x_param, y_param, flux_param);
259 auto x_param = std::make_shared<ManualParameter>(
x);
260 auto y_param = std::make_shared<ManualParameter>(
y);
261 auto flux_param = std::make_shared<ManualParameter>(flux);
263 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, image->getValue(
x,
y) + bg_noise_dist(m_rng));
280 for (
int y=0;
y < image->getHeight();
y++) {
281 for (
int x=0;
x < image->getWidth();
x++) {
282 auto pixel_value = image->getValue(
x,
y);
283 if (pixel_value > 0.) {
284 image->setValue(
x,
y, boost::random::poisson_distribution<>(pixel_value * gain)(m_rng) / gain);
292 if (saturation_level > 0.0) {
293 for (
int y=0;
y < image->getHeight();
y++) {
294 for (
int x=0;
x < image->getWidth();
x++) {
295 image->setValue(
x,
y,
std::min(image->getValue(
x,
y), (
float) saturation_level));
303 for (
int y=0;
y < weight_map->getHeight();
y++) {
304 for (
int x=0;
x < weight_map->getWidth();
x++) {
305 if (boost::random::uniform_01<double>()(m_rng) < probability) {
306 weight_map->setValue(
x,
y, 0);
315 for (
int x=0;
x < weight_map->getWidth();
x++) {
316 if (boost::random::uniform_01<double>()(m_rng) < probability) {
317 for (
int y=0;
y < weight_map->getHeight();
y++) {
318 weight_map->setValue(
x,
y, 0);
326 double x_min,
double y_min,
double x_max,
double y_max) {
331 boost::random::uniform_real_distribution<> random_x(x_min, x_max);
332 boost::random::uniform_real_distribution<> random_y(y_min, y_max);
336 for (
int i = 0; i < number_of_sources; i++) {
339 random_x(m_rng), random_y(m_rng),
341 gal_exp_random_i0(m_rng),
342 boost::random::uniform_real_distribution<>(.5, 6)(m_rng),
343 boost::random::uniform_real_distribution<>(.2, .8)(m_rng),
344 boost::random::uniform_real_distribution<>(-M_PI/2, M_PI/2)(m_rng),
365 while (file.
good()) {
370 if (line.
size() == 0) {
377 linestream >> source.
x >> source.
y
382 source.
exp_rot *= -M_PI / 180.0;
383 source.
dev_rot *= -M_PI / 180.0;
395 for (
const auto& source : sources) {
396 file << source.x <<
" " << source.y <<
" "
397 << source.exp_flux <<
" " << source.exp_rad <<
" " << source.exp_aspect <<
" " << (source.exp_rot * -180.0 / M_PI) <<
" "
398 << source.dev_flux <<
" " << source.dev_rad <<
" " << source.dev_aspect <<
" " << (source.dev_rot * -180.0 / M_PI) <<
" "
399 << source.point_flux <<
"\n";
404 auto center = image_size / 2.0;
405 auto c =
cos(rot_angle);
406 auto s =
sin(rot_angle);
407 for (
auto& source : sources) {
410 double x = (source.x * c - source.y *
s) * scale + shift_x;
411 double y = (source.x *
s + source.y * c) * scale + shift_y;
412 source.x =
x + center;
413 source.y =
y + center;
415 source.exp_rot -= rot_angle;
416 source.dev_rot -= rot_angle;
428 while (file.
good()) {
433 if (line.
size() == 0) {
439 linestream >> source_type;
444 if (source_type == 200) {
445 double magnitude, bt, exp_angle_degree, dev_angle_degree;
447 linestream >> source.
x >> source.
y >> magnitude
453 source.
dev_rot = dev_angle_degree * M_PI / 180.0;
454 source.
exp_rot = exp_angle_degree * M_PI / 180.0;
456 auto total_flux =
pow(10, (magnitude - m_zero_point) / -2.5);
457 source.
exp_flux = total_flux * (1.0 - bt);
464 }
else if (source_type == 100) {
466 linestream >> source.
x >> source.
y >> magnitude;
467 source.
point_flux =
pow(10, (magnitude - m_zero_point) / -2.5) * m_exp_time;
486 auto max_tile_memory = args[
"max-tile-memory"].as<
int>();
487 auto tile_size = args[
"tile-size"].as<
int>();
490 auto image_size = args[
"size"].as<
double>();
491 auto rot_angle = args[
"rotation"].as<
double>() / 180.0 * M_PI;
492 auto scale = args[
"scale"].as<
double>();
493 auto shift_x = args[
"shift-x"].as<
double>();
494 auto shift_y = args[
"shift-y"].as<
double>();
495 auto model_size = args[
"model-size"].as<
double>();
496 if (model_size <= 0.) {
497 model_size = image_size;
500 m_zero_point = args[
"zero-point"].as<
double>();
501 m_exp_time = args[
"exposure-time"].as<
double>();
509 auto psf_filename = args[
"psf-file"].as<
std::string>();
510 if (psf_filename !=
"") {
511 logger.
info() <<
"Loading psf file: " << psf_filename;
518 auto copy_coordinate_system = args[
"copy-coordinate-system"].as<
std::string>();
519 if (copy_coordinate_system !=
"") {
520 coordinate_system = std::make_shared<WCS>(copy_coordinate_system);
522 coordinate_system = std::make_shared<DummyWCS>(image_size, image_size, rot_angle, scale, shift_x, shift_y);
527 logger.
fatal() <<
"The expected required memory for model rasterization exceeds the maximum size for an integer";
528 logger.
fatal() <<
"Please, either reduce the model size, the image size, or increase the PSF pixel scale";
536 for (
auto i = 0u; i < psf_vals.
size(); ++i) {
537 if (vpsf_components[i].name ==
"X_IMAGE" || vpsf_components[i].name ==
"Y_IMAGE") {
538 psf_vals[i] = image_size / 2 - 1;
546 auto p = vpsf->
getPsf(psf_vals);
547 auto psf_sum =
std::accumulate(p->getData().begin(), p->getData().end(), 0.);
554 if (sources_filename !=
"") {
555 logger.
info() <<
"Loading sources from < " << sources_filename <<
" >...";
556 auto loaded_sources = loadSourcesFromFile(sources_filename);
557 sources.
insert(sources.
end(), loaded_sources.begin(), loaded_sources.end());
561 if (sources_cat_filename !=
"") {
562 logger.
info() <<
"Loading source catalog from < " << sources_cat_filename <<
" >...";
563 auto loaded_sources = loadSourcesFromCatalog(sources_cat_filename);
564 sources.
insert(sources.
end(), loaded_sources.begin(), loaded_sources.end());
567 int nb_of_random_sources = args[
"random-sources"].as<
int>();
568 if (nb_of_random_sources > 0) {
569 logger.
info() <<
"Adding " << nb_of_random_sources <<
" random source" << (nb_of_random_sources>1 ?
"s" :
"") <<
"...";
570 auto random_sources = generateRandomSources(nb_of_random_sources, image_size * .1, image_size * .1, image_size * .9, image_size * .9);
571 sources.
insert(sources.
end(), random_sources.begin(), random_sources.end());
575 if (sources_save_filename !=
"") {
576 logger.
info() <<
"Saving sources to < " << sources_save_filename <<
" >...";
577 saveSources(sources, sources_save_filename);
581 transformSources(sources, image_size, rot_angle, scale, shift_x, shift_y);
584 for (
const auto& source : sources) {
585 addSource(point_models, extended_models, model_size, source,
std::make_tuple(scale, 0, 0, scale));
591 auto target_image_source = std::make_shared<FitsImageSource<SeFloat>>(
filename, image_size, image_size, coordinate_system);
594 if (args[
"disable-psf"].as<bool>()) {
602 frame_model.rasterToImage(target_image);
612 frame_model.rasterToImage(target_image);
618 addPoissonNoise(target_image, args[
"gain"].as<double>());
619 addBackgroundNoise(target_image, args[
"bg-level"].as<double>(), args[
"bg-sigma"].as<double>());
623 auto saturation_level = args[
"saturation"].as<
double>();
624 saturate(target_image, saturation_level);
628 auto weight_filename = args[
"output-weight"].as<
std::string>();
629 if (weight_filename !=
"") {
630 auto weight_map_source = std::make_shared<FitsImageSource<SeFloat>>(weight_filename, image_size, image_size);
632 for (
int y = 0;
y < image_size; ++
y) {
633 for (
int x = 0;
x < image_size; ++
x) {
634 weight_map->setValue(
x,
y, 1);
640 addBadPixels(weight_map, args[
"bad-pixels"].as<double>());
641 addBadColumns(weight_map, args[
"bad-columns"].as<double>());
643 for (
int y = 0;
y < target_image->getHeight();
y++) {
644 for (
int x = 0;
x < target_image->getWidth();
x++) {
645 if (weight_map->getValue(
x,
y) == 0) {
646 target_image->setValue(
x,
y, saturation_level);
657 boost::random::mt19937 m_rng { (
unsigned int)
time(NULL) } ;
658 double m_zero_point = 0.0, m_exp_time = 300.;