26 #include <CCfits/CCfits>
31 namespace po = boost::program_options;
36 namespace SourceXtractor {
45 CCfits::ExtHDU &psf_data = pFits->extension(hdu_number);
46 if (psf_data.name() !=
"PSF_DATA") {
47 throw Elements::Exception() <<
"No PSFEX data in file " << pFits->name() <<
" HDU " << hdu_number;
51 psf_data.readKey(
"POLNAXIS", n_components);
53 double pixel_sampling;
54 psf_data.readKey(
"PSF_SAMP", pixel_sampling);
58 for (
int i = 0; i < n_components; ++i) {
61 psf_data.readKey(
"POLGRP" + index_str, components[i].group_id);
62 psf_data.readKey(
"POLNAME" + index_str, components[i].name);
63 psf_data.readKey(
"POLZERO" + index_str, components[i].offset);
64 psf_data.readKey(
"POLSCAL" + index_str, components[i].scale);
67 --components[i].group_id;
71 psf_data.readKey(
"POLNGRP", n_groups);
75 for (
int i = 0; i < n_groups; ++i) {
78 psf_data.readKey(
"POLDEG" + index_str, group_degrees[i]);
82 psf_data.readKey(
"PSFAXIS1", width);
83 psf_data.readKey(
"PSFAXIS2", height);
84 psf_data.readKey(
"PSFAXIS3", n_coeffs);
86 if (width != height) {
87 throw Elements::Exception() <<
"Non square PSFEX format PSF (" << width <<
'X' << height <<
')';
96 psf_data.column(
"PSF_MASK").readArrays(all_data, 0, 0);
97 auto& raw_coeff_data = all_data[0];
101 for (
int i = 0; i < n_coeffs; ++i) {
102 auto offset =
std::begin(raw_coeff_data) + i * n_pixels;
106 logger.
debug() <<
"Loaded variable PSF " << pFits->name() <<
" (" << width <<
", " << height <<
") with "
107 << n_coeffs <<
" coefficients";
109 ll <<
"Components: ";
110 for (
auto c : components) {
117 return std::make_shared<VariablePsf>(pixel_sampling, components, group_degrees, coefficients);
119 }
catch (CCfits::FitsException &
e) {
127 double pixel_sampling;
129 image_hdu.readKey(
"SAMPLING", pixel_sampling);
131 catch (CCfits::HDU::NoSuchKeyword&) {
135 if (image_hdu.axis(0) != image_hdu.axis(1)) {
137 << image_hdu.axis(1) <<
')';
140 auto size = image_hdu.axis(0);
143 image_hdu.read(data);
147 logger.
debug() <<
"Loaded image PSF(" << size <<
", " << size <<
") with sampling step " << pixel_sampling;
149 return std::make_shared<VariablePsf>(pixel_sampling, kernel);
158 auto& image_hdu = pFits->pHDU();
160 auto axes = image_hdu.axes();
163 if (hdu_number == 1) {
166 auto& extension = pFits->extension(hdu_number - 1);
174 }
catch (CCfits::FitsException &
e) {
180 int size = int(fwhm / pixel_sampling + 1) * 6 + 1;
184 double sigma_squared = (fwhm / (pixel_sampling * 2.355)) * (fwhm / (pixel_sampling * 2.355));
187 for (
int y = 0;
y < size;
y++) {
188 for (
int x = 0;
x < size;
x++) {
189 double sx = (
x - size / 2);
190 double sy = (
y - size / 2);
191 double pixel_value =
exp(-(sx * sx + sy * sy) / (2 * sigma_squared));
192 total += pixel_value;
193 kernel->setValue(
x,
y, pixel_value);
196 for (
int y = 0;
y < size;
y++) {
197 for (
int x = 0;
x < size;
x++) {
198 kernel->setValue(
x,
y, kernel->getValue(
x,
y) / total);
202 logger.
info() <<
"Using gaussian PSF(" << fwhm <<
") with pixel sampling step size " << pixel_sampling <<
" (size " << size <<
")";
204 return std::make_shared<VariablePsf>(pixel_sampling, kernel);
208 return {{
"Variable PSF", {
210 "Psf image file (FITS format)."},
212 "Generate a gaussian PSF with the given full-width half-maximum (in pixels)"},
214 "Generate a gaussian PSF with the given pixel sampling step size"}
220 if (args.find(
PSF_FILE) == args.end() && args.find(
PSF_FWHM) != args.end() &&
227 if (args.find(
PSF_FILE) != args.end()) {
229 }
else if (args.find(
PSF_FWHM) != args.end()) {
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void info(const std::string &logMessage)
void debug(const std::string &logMessage)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
static Logging getLogger(const std::string &name="")