SourceXtractorPlusPlus  0.13
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PsfPluginConfig.cpp
Go to the documentation of this file.
1 
17 /*
18  * PsfPluginConfig.cpp
19  *
20  * Created on: Jun 25, 2018
21  * Author: Alejandro Álvarez Ayllón (greatly adapted from mschefer's code)
22  */
23 
26 #include <CCfits/CCfits>
27 #include <ElementsKernel/Logging.h>
30 
31 namespace po = boost::program_options;
33 
34 static auto logger = Elements::Logging::getLogger("PsfPlugin");
35 
36 namespace SourceXtractor {
37 
38 static const std::string PSF_FILE{"psf-filename"};
39 static const std::string PSF_FWHM {"psf-fwhm" };
40 static const std::string PSF_PIXEL_SAMPLING {"psf-pixel-sampling" };
41 
42 
44  try {
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;
48  }
49 
50  int n_components;
51  psf_data.readKey("POLNAXIS", n_components);
52 
53  double pixel_sampling;
54  psf_data.readKey("PSF_SAMP", pixel_sampling);
55 
56  std::vector<VariablePsf::Component> components(n_components);
57 
58  for (int i = 0; i < n_components; ++i) {
59  auto index_str = std::to_string(i + 1);
60 
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);
65 
66  // Groups in the FITS file are indexed starting at 1, but we use a zero-based index
67  --components[i].group_id;
68  }
69 
70  int n_groups;
71  psf_data.readKey("POLNGRP", n_groups);
72 
73  std::vector<int> group_degrees(n_groups);
74 
75  for (int i = 0; i < n_groups; ++i) {
76  auto index_str = std::to_string(i + 1);
77 
78  psf_data.readKey("POLDEG" + index_str, group_degrees[i]);
79  }
80 
81  int width, height, n_coeffs, n_pixels;
82  psf_data.readKey("PSFAXIS1", width);
83  psf_data.readKey("PSFAXIS2", height);
84  psf_data.readKey("PSFAXIS3", n_coeffs);
85 
86  if (width != height) {
87  throw Elements::Exception() << "Non square PSFEX format PSF (" << width << 'X' << height << ')';
88  }
89  if (width % 2 == 0) {
90  throw Elements::Exception() << "PSF kernel must have odd size";
91  }
92 
93  n_pixels = width * height;
94 
96  psf_data.column("PSF_MASK").readArrays(all_data, 0, 0);
97  auto& raw_coeff_data = all_data[0];
98 
99  std::vector<std::shared_ptr<VectorImage<SeFloat>>> coefficients(n_coeffs);
100 
101  for (int i = 0; i < n_coeffs; ++i) {
102  auto offset = std::begin(raw_coeff_data) + i * n_pixels;
103  coefficients[i] = VectorImage<SeFloat>::create(width, height, offset, offset + n_pixels);
104  }
105 
106  logger.debug() << "Loaded variable PSF " << pFits->name() << " (" << width << ", " << height << ") with "
107  << n_coeffs << " coefficients";
108  auto ll = logger.debug();
109  ll << "Components: ";
110  for (auto c : components) {
111  ll << c.name << " ";
112  if (component_value_getters.find(c.name) == component_value_getters.end()) {
113  throw Elements::Exception() << "Can not find a getter for the component " << c.name;
114  }
115  }
116 
117  return std::make_shared<VariablePsf>(pixel_sampling, components, group_degrees, coefficients);
118 
119  } catch (CCfits::FitsException &e) {
120  throw Elements::Exception() << "Error loading PSFEx file: " << e.message();
121  }
122 }
123 
124 // templated to work around CCfits limitation that primary and extension HDUs are different classes
125 template<typename T>
127  double pixel_sampling;
128  try {
129  image_hdu.readKey("SAMPLING", pixel_sampling);
130  }
131  catch (CCfits::HDU::NoSuchKeyword&) {
132  pixel_sampling = 1.;
133  }
134 
135  if (image_hdu.axis(0) != image_hdu.axis(1)) {
136  throw Elements::Exception() << "Non square image PSF (" << image_hdu.axis(0) << 'X'
137  << image_hdu.axis(1) << ')';
138  }
139 
140  auto size = image_hdu.axis(0);
141 
142  std::valarray<double> data{};
143  image_hdu.read(data);
144  auto kernel = VectorImage<SeFloat>::create(size, size);
145  std::copy(begin(data), end(data), kernel->getData().begin());
146 
147  logger.debug() << "Loaded image PSF(" << size << ", " << size << ") with sampling step " << pixel_sampling;
148 
149  return std::make_shared<VariablePsf>(pixel_sampling, kernel);
150 }
151 
155  try {
156  // Read the HDU from the file
157  std::unique_ptr<CCfits::FITS> pFits{new CCfits::FITS(filename, CCfits::Read)};
158  auto& image_hdu = pFits->pHDU();
159 
160  auto axes = image_hdu.axes();
161  // PSF as image
162  if (axes == 2) {
163  if (hdu_number == 1) {
164  return readImage(image_hdu);
165  } else {
166  auto& extension = pFits->extension(hdu_number - 1);
167  return readImage(extension);
168  }
169  }
170  // PSFEx format
171  else {
172  return readPsfEx(pFits, hdu_number);
173  }
174  } catch (CCfits::FitsException &e) {
175  throw Elements::Exception() << "Error loading PSF file: " << e.message();
176  }
177 }
178 
180  int size = int(fwhm / pixel_sampling + 1) * 6 + 1;
181  auto kernel = VectorImage<SeFloat>::create(size, size);
182 
183  // convert full width half maximum to standard deviation
184  double sigma_squared = (fwhm / (pixel_sampling * 2.355)) * (fwhm / (pixel_sampling * 2.355));
185 
186  double total = 0;
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);
194  }
195  }
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);
199  }
200  }
201 
202  logger.info() << "Using gaussian PSF(" << fwhm << ") with pixel sampling step size " << pixel_sampling << " (size " << size << ")";
203 
204  return std::make_shared<VariablePsf>(pixel_sampling, kernel);
205 }
206 
208  return {{"Variable PSF", {
209  {PSF_FILE.c_str(), po::value<std::string>(),
210  "Psf image file (FITS format)."},
211  {PSF_FWHM.c_str(), po::value<double>(),
212  "Generate a gaussian PSF with the given full-width half-maximum (in pixels)"},
213  {PSF_PIXEL_SAMPLING.c_str(), po::value<double>(),
214  "Generate a gaussian PSF with the given pixel sampling step size"}
215  }}};
216 }
217 
218 void PsfPluginConfig::preInitialize(const Euclid::Configuration::Configuration::UserValues &args) {
219  // Fail if there is no PSF file, but PSF_FWHM is set and PSF_PIXEL_SAMPLING is not
220  if (args.find(PSF_FILE) == args.end() && args.find(PSF_FWHM) != args.end() &&
221  args.find(PSF_PIXEL_SAMPLING) == args.end()) {
222  throw Elements::Exception(PSF_PIXEL_SAMPLING + " is required when using " + PSF_FWHM);
223  }
224 }
225 
226 void PsfPluginConfig::initialize(const UserValues &args) {
227  if (args.find(PSF_FILE) != args.end()) {
228  m_vpsf = readPsf(args.find(PSF_FILE)->second.as<std::string>());
229  } else if (args.find(PSF_FWHM) != args.end()) {
230  m_vpsf = generateGaussianPsf(args.find(PSF_FWHM)->second.as<double>(),
231  args.find(PSF_PIXEL_SAMPLING)->second.as<double>());
232  }
233 }
234 
236  return m_vpsf;
237 }
238 
239 } // end SourceXtractor
static const std::string PSF_FWHM
std::map< std::string, OptionDescriptionList > getProgramOptions() override
T copy(T...args)
static std::shared_ptr< VariablePsf > readImage(T &image_hdu)
static auto logger
Definition: WCS.cpp:46
T exp(T...args)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
constexpr double e
T to_string(T...args)
void initialize(const UserValues &args) override
static std::shared_ptr< VariablePsf > readPsfEx(std::unique_ptr< CCfits::FITS > &pFits, int hdu_number=1)
const std::shared_ptr< VariablePsf > & getPsf() const
SeFloat32 SeFloat
Definition: Types.h:32
T end(T...args)
static std::shared_ptr< VectorImage< T > > create(Args &&...args)
Definition: VectorImage.h:89
STL class.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
STL class.
void preInitialize(const UserValues &args) override
static const std::string PSF_FILE
string filename
Definition: conf.py:63
STL class.
static const std::string PSF_PIXEL_SAMPLING
std::map< std::string, ValueGetter > component_value_getters
Definition: PsfTask.cpp:45
STL class.
static std::shared_ptr< VariablePsf > generateGaussianPsf(SeFloat fwhm, SeFloat pixel_sampling)
T begin(T...args)
static std::shared_ptr< VariablePsf > readPsf(const std::string &filename, int hdu_number=1)
T c_str(T...args)
std::shared_ptr< VariablePsf > m_vpsf
static auto logger
static Logging getLogger(const std::string &name="")