SourceXtractorPlusPlus  0.15
Please provide a description of the project.
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 
28 #include <CCfits/CCfits>
30 #include <ElementsKernel/Logging.h>
31 
32 namespace po = boost::program_options;
34 
35 static auto logger = Elements::Logging::getLogger("PsfPlugin");
36 
37 namespace SourceXtractor {
38 
39 static const std::string PSF_FILE{"psf-filename"};
40 static const std::string PSF_FWHM {"psf-fwhm" };
41 static const std::string PSF_PIXEL_SAMPLING {"psf-pixel-sampling" };
42 
43 /*
44  * Reading in a stacked PSF as it is being developed for co-added images in Euclid
45  */
47  logger.debug() << "Loading a PSF stack file.";
48  std::shared_ptr<VariablePsfStack> act_stack = std::make_shared<VariablePsfStack>(std::move(pFits));
49  return act_stack;
50 }
51 
53  try {
54  CCfits::ExtHDU& psf_data = pFits->extension("PSF_DATA");
55 
56  int n_components;
57  psf_data.readKey("POLNAXIS", n_components);
58 
59  double pixel_sampling;
60  psf_data.readKey("PSF_SAMP", pixel_sampling);
61 
62  std::vector<VariablePsf::Component> components(n_components);
63 
64  for (int i = 0; i < n_components; ++i) {
65  auto index_str = std::to_string(i + 1);
66 
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);
71 
72  // Groups in the FITS file are indexed starting at 1, but we use a zero-based index
73  --components[i].group_id;
74  }
75 
76  int n_groups;
77  psf_data.readKey("POLNGRP", n_groups);
78 
79  std::vector<int> group_degrees(n_groups);
80 
81  for (int i = 0; i < n_groups; ++i) {
82  auto index_str = std::to_string(i + 1);
83 
84  psf_data.readKey("POLDEG" + index_str, group_degrees[i]);
85  }
86 
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);
91 
92  if (width != height) {
93  throw Elements::Exception() << "Non square PSFEX format PSF (" << width << 'X' << height << ')';
94  }
95  if (width % 2 == 0) {
96  throw Elements::Exception() << "PSF kernel must have odd size";
97  }
98 
99  n_pixels = width * height;
100 
102  psf_data.column("PSF_MASK").readArrays(all_data, 0, 0);
103  auto& raw_coeff_data = all_data[0];
104 
105  std::vector<std::shared_ptr<VectorImage<SeFloat>>> coefficients(n_coeffs);
106 
107  for (int i = 0; i < n_coeffs; ++i) {
108  auto offset = std::begin(raw_coeff_data) + i * n_pixels;
109  coefficients[i] = VectorImage<SeFloat>::create(width, height, offset, offset + n_pixels);
110  }
111 
112  logger.debug() << "Loaded variable PSF " << pFits->name() << " (" << width << ", " << height << ") with " << n_coeffs
113  << " coefficients";
114  auto ll = logger.debug();
115  ll << "Components: ";
116  for (auto c : components) {
117  ll << c.name << " ";
118  if (component_value_getters.find(c.name) == component_value_getters.end()) {
119  throw Elements::Exception() << "Can not find a getter for the component " << c.name;
120  }
121  }
122 
123  return std::make_shared<VariablePsf>(pixel_sampling, components, group_degrees, coefficients);
124  } catch (CCfits::FITS::NoSuchHDU&) {
125  throw;
126  } catch (CCfits::FitsException &e) {
127  throw Elements::Exception() << "Error loading PSFEx file: " << e.message();
128  }
129 }
130 
131 // templated to work around CCfits limitation that primary and extension HDUs are different classes
132 template<typename T>
134  double pixel_sampling;
135  try {
136  image_hdu.readKey("SAMPLING", pixel_sampling);
137  }
138  catch (CCfits::HDU::NoSuchKeyword&) {
139  pixel_sampling = 1.;
140  }
141 
142  if (image_hdu.axis(0) != image_hdu.axis(1)) {
143  throw Elements::Exception() << "Non square image PSF (" << image_hdu.axis(0) << 'X'
144  << image_hdu.axis(1) << ')';
145  }
146 
147  auto size = image_hdu.axis(0);
148 
149  std::valarray<double> data{};
150  image_hdu.read(data);
151  auto kernel = VectorImage<SeFloat>::create(size, size);
152  std::copy(begin(data), end(data), kernel->getData().begin());
153 
154  logger.debug() << "Loaded image PSF(" << size << ", " << size << ") with sampling step " << pixel_sampling;
155 
156  return std::make_shared<VariablePsf>(pixel_sampling, kernel);
157 }
158 
162  try {
163  // Read the HDU from the file
164  std::unique_ptr<CCfits::FITS> pFits{new CCfits::FITS(filename, CCfits::Read)};
165  auto& image_hdu = pFits->pHDU();
166 
167  auto axes = image_hdu.axes();
168  // PSF as image
169  if (axes == 2) {
170  if (hdu_number == 1) {
171  return readImage(image_hdu);
172  } else {
173  auto& extension = pFits->extension(hdu_number - 1);
174  return readImage(extension);
175  }
176  } else {
177  try {
178  // PSFEx format
179  return readPsfEx(pFits);
180  } catch (CCfits::FITS::NoSuchHDU& e) {
181  logger.debug() << "Failed Reading a PsfEx file!";
182  return readStackedPsf(pFits);
183  }
184  }
185  } catch (CCfits::FitsException& e) {
186  throw Elements::Exception() << "Error loading PSF file: " << e.message();
187  }
188 }
189 
191  int size = int(fwhm / pixel_sampling + 1) * 6 + 1;
192  auto kernel = VectorImage<SeFloat>::create(size, size);
193 
194  // convert full width half maximum to standard deviation
195  double sigma_squared = (fwhm / (pixel_sampling * 2.355)) * (fwhm / (pixel_sampling * 2.355));
196 
197  double total = 0;
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);
205  }
206  }
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);
210  }
211  }
212 
213  logger.info() << "Using gaussian PSF(" << fwhm << ") with pixel sampling step size " << pixel_sampling << " (size " << size << ")";
214 
215  return std::make_shared<VariablePsf>(pixel_sampling, kernel);
216 }
217 
219  return {{"Variable PSF", {
220  {PSF_FILE.c_str(), po::value<std::string>(),
221  "Psf image file (FITS format)."},
222  {PSF_FWHM.c_str(), po::value<double>(),
223  "Generate a gaussian PSF with the given full-width half-maximum (in pixels)"},
224  {PSF_PIXEL_SAMPLING.c_str(), po::value<double>(),
225  "Generate a gaussian PSF with the given pixel sampling step size"}
226  }}};
227 }
228 
230  // Fail if there is no PSF file, but PSF_FWHM is set and PSF_PIXEL_SAMPLING is not
231  if (args.find(PSF_FILE) == args.end() && args.find(PSF_FWHM) != args.end() &&
232  args.find(PSF_PIXEL_SAMPLING) == args.end()) {
233  throw Elements::Exception(PSF_PIXEL_SAMPLING + " is required when using " + PSF_FWHM);
234  }
235 }
236 
238  if (args.find(PSF_FILE) != args.end()) {
239  m_vpsf = readPsf(args.find(PSF_FILE)->second.as<std::string>());
240  } else if (args.find(PSF_FWHM) != args.end()) {
241  m_vpsf = generateGaussianPsf(args.find(PSF_FWHM)->second.as<double>(),
242  args.find(PSF_PIXEL_SAMPLING)->second.as<double>());
243  }
244 }
245 
247  return m_vpsf;
248 }
249 
250 } // end SourceXtractor
PsfPluginConfig.h
SourceXtractor::PsfPluginConfig::m_vpsf
std::shared_ptr< Psf > m_vpsf
Definition: PsfPluginConfig.h:49
std::string
STL class.
std::shared_ptr
STL class.
std::move
T move(T... args)
conf.filename
string filename
Definition: conf.py:63
SourceXtractor::SeFloat
SeFloat32 SeFloat
Definition: Types.h:32
std::vector
STL class.
SourceXtractor::readImage
static std::shared_ptr< VariablePsf > readImage(T &image_hdu)
Definition: PsfPluginConfig.cpp:133
std::map::find
T find(T... args)
SourceXtractor::PSF_FILE
static const std::string PSF_FILE
Definition: PsfPluginConfig.cpp:39
SourceXtractor::PsfPluginConfig::generateGaussianPsf
static std::shared_ptr< Psf > generateGaussianPsf(SeFloat fwhm, SeFloat pixel_sampling)
Definition: PsfPluginConfig.cpp:190
SourceXtractor::PsfPluginConfig::getPsf
const std::shared_ptr< Psf > & getPsf() const
Definition: PsfPluginConfig.cpp:246
SourceXtractor::readStackedPsf
static std::shared_ptr< VariablePsfStack > readStackedPsf(std::unique_ptr< CCfits::FITS > &pFits)
Definition: PsfPluginConfig.cpp:46
SourceXtractor::readPsfEx
static std::shared_ptr< VariablePsf > readPsfEx(std::unique_ptr< CCfits::FITS > &pFits)
Definition: PsfPluginConfig.cpp:52
SourceXtractor
Definition: Aperture.h:30
logger
static auto logger
Definition: PsfPluginConfig.cpp:35
SourceXtractor::PsfPluginConfig::getProgramOptions
std::map< std::string, OptionDescriptionList > getProgramOptions() override
Definition: PsfPluginConfig.cpp:218
std::string::c_str
T c_str(T... args)
std::to_string
T to_string(T... args)
VariablePsfStack.h
SourceXtractor::VectorImage::create
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
std::valarray
STL class.
std::copy
T copy(T... args)
Elements::Exception
std::map
STL class.
SourceXtractor::logger
static auto logger
Definition: WCS.cpp:44
e
constexpr double e
Elements::Logging::getLogger
static Logging getLogger(const std::string &name="")
SourceXtractor::PsfPluginConfig::readPsf
static std::shared_ptr< Psf > readPsf(const std::string &filename, int hdu_number=1)
Definition: PsfPluginConfig.cpp:161
x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
Definition: MoffatModelFittingTask.cpp:94
PsfTask.h
std::exp
T exp(T... args)
std::begin
T begin(T... args)
SourceXtractor::PSF_FWHM
static const std::string PSF_FWHM
Definition: PsfPluginConfig.cpp:40
VariablePsf.h
SourceXtractor::PsfPluginConfig::initialize
void initialize(const UserValues &args) override
Definition: PsfPluginConfig.cpp:237
Euclid::Configuration::Configuration
Logging.h
std::end
T end(T... args)
SourceXtractor::PSF_PIXEL_SAMPLING
static const std::string PSF_PIXEL_SAMPLING
Definition: PsfPluginConfig.cpp:41
y
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
Definition: MoffatModelFittingTask.cpp:94
std::unique_ptr< CCfits::FITS >
SourceXtractor::component_value_getters
std::map< std::string, ValueGetter > component_value_getters
Definition: PsfTask.cpp:45
ProgramOptionsHelper.h
SourceXtractor::PsfPluginConfig::preInitialize
void preInitialize(const UserValues &args) override
Definition: PsfPluginConfig.cpp:229