28 #include <CCfits/CCfits>
32 namespace po = boost::program_options;
47 logger.debug() <<
"Loading a PSF stack file.";
54 CCfits::ExtHDU& psf_data = pFits->extension(
"PSF_DATA");
57 psf_data.readKey(
"POLNAXIS", n_components);
59 double pixel_sampling;
60 psf_data.readKey(
"PSF_SAMP", pixel_sampling);
64 for (
int i = 0; i < n_components; ++i) {
67 psf_data.readKey(
"POLGRP" + index_str, components[i].group_id);
68 psf_data.readKey(
"POLNAME" + index_str, components[i].name);
69 psf_data.readKey(
"POLZERO" + index_str, components[i].offset);
70 psf_data.readKey(
"POLSCAL" + index_str, components[i].scale);
73 --components[i].group_id;
77 psf_data.readKey(
"POLNGRP", n_groups);
81 for (
int i = 0; i < n_groups; ++i) {
84 psf_data.readKey(
"POLDEG" + index_str, group_degrees[i]);
87 int width, height, n_coeffs, n_pixels;
88 psf_data.readKey(
"PSFAXIS1", width);
89 psf_data.readKey(
"PSFAXIS2", height);
90 psf_data.readKey(
"PSFAXIS3", n_coeffs);
92 if (width != height) {
93 throw Elements::Exception() <<
"Non square PSFEX format PSF (" << width <<
'X' << height <<
')';
99 n_pixels = width * height;
102 psf_data.column(
"PSF_MASK").readArrays(all_data, 0, 0);
103 auto& raw_coeff_data = all_data[0];
107 for (
int i = 0; i < n_coeffs; ++i) {
108 auto offset =
std::begin(raw_coeff_data) + i * n_pixels;
112 logger.debug() <<
"Loaded variable PSF " << pFits->name() <<
" (" << width <<
", " << height <<
") with " << n_coeffs
115 ll <<
"Components: ";
116 for (
auto c : components) {
123 return std::make_shared<VariablePsf>(pixel_sampling, components, group_degrees, coefficients);
124 }
catch (CCfits::FITS::NoSuchHDU&) {
126 }
catch (CCfits::FitsException &
e) {
134 double pixel_sampling;
136 image_hdu.readKey(
"SAMPLING", pixel_sampling);
138 catch (CCfits::HDU::NoSuchKeyword&) {
142 if (image_hdu.axis(0) != image_hdu.axis(1)) {
144 << image_hdu.axis(1) <<
')';
147 auto size = image_hdu.axis(0);
150 image_hdu.read(data);
154 logger.debug() <<
"Loaded image PSF(" << size <<
", " << size <<
") with sampling step " << pixel_sampling;
156 return std::make_shared<VariablePsf>(pixel_sampling, kernel);
165 auto& image_hdu = pFits->pHDU();
167 auto axes = image_hdu.axes();
170 if (hdu_number == 1) {
173 auto& extension = pFits->extension(hdu_number - 1);
180 }
catch (CCfits::FITS::NoSuchHDU&
e) {
181 logger.debug() <<
"Failed Reading a PsfEx file!";
185 }
catch (CCfits::FitsException&
e) {
191 int size = int(fwhm / pixel_sampling + 1) * 6 + 1;
195 double sigma_squared = (fwhm / (pixel_sampling * 2.355)) * (fwhm / (pixel_sampling * 2.355));
198 for (
int y = 0;
y < size;
y++) {
199 for (
int x = 0;
x < size;
x++) {
200 double sx = (
x - size / 2);
201 double sy = (
y - size / 2);
202 double pixel_value =
exp(-(sx * sx + sy * sy) / (2 * sigma_squared));
203 total += pixel_value;
204 kernel->setValue(
x,
y, pixel_value);
207 for (
int y = 0;
y < size;
y++) {
208 for (
int x = 0;
x < size;
x++) {
209 kernel->setValue(
x,
y, kernel->getValue(
x,
y) / total);
213 logger.info() <<
"Using gaussian PSF(" << fwhm <<
") with pixel sampling step size " << pixel_sampling <<
" (size " << size <<
")";
215 return std::make_shared<VariablePsf>(pixel_sampling, kernel);
219 return {{
"Variable PSF", {
221 "Psf image file (FITS format)."},
223 "Generate a gaussian PSF with the given full-width half-maximum (in pixels)"},
225 "Generate a gaussian PSF with the given pixel sampling step size"}
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
static Logging getLogger(const std::string &name="")