SourceXtractorPlusPlus  0.15
Please provide a description of the project.
VariablePsfStack.cpp
Go to the documentation of this file.
1 
17 /*
18  * VariablePsfStack.cpp
19  *
20  * Created on: July 10, 2019
21  * Author: Martin Kuemmel
22  */
23 #include <algorithm>
24 #include <ElementsKernel/Logging.h>
27 
28 static auto stack_logger = Elements::Logging::getLogger("VarStackPsf");
29 
30 namespace SourceXtractor {
31 
33  try {
34  // basic check: load the primary HDU
35  CCfits::PHDU& phdu = pFits->pHDU();
36  if (phdu.axes() != 0) {
37  throw Elements::Exception() << "The primary HDU is not empty! File: " << pFits->name();
38  }
39 
40  // basic check: load the first extension
41  CCfits::ExtHDU& psf_data = pFits->extension(1);
42  if (psf_data.axes() != 2) {
43  throw Elements::Exception() << "The first HDU is not an image! File: " << pFits->name();
44  }
45 
46  // basic check: load the second extension
47  CCfits::ExtHDU& position_data = pFits->extension(2);
48  if (position_data.axes() != 2) {
49  throw Elements::Exception() << "The second HDU has not the right dimension! File: " << pFits->name();
50  }
51 
52  // give some feedback
53  stack_logger.debug() << "Checked the file: " << pFits->name();
54 
55  // get and store the stamp size;
56  // define the offset from GRIDX and GRIDY
57  psf_data.readKey("STMPSIZE", m_psf_size);
58  m_grid_offset = m_psf_size / 2; // this is for CCfits where indices start with 1!
59 
60  // make sure the PSF size is odd
61  if (m_psf_size % 2 == 0) {
62  throw Elements::Exception() << "PSF kernel must have odd size, but has: " << m_psf_size;
63  }
64 
65  try {
66  // try to get the sampling
67  psf_data.readKey("SAMPLING", mm_pixel_sampling);
68  } catch (CCfits::HDU::NoSuchKeyword&) {
69  // use a default value
70  mm_pixel_sampling = 1.;
71  }
72 
73  // read the nrows value
74  m_nrows = position_data.rows();
75 
76  try {
77  // read in all the EXT specific columns
78  position_data.column("GRIDX", false).read(m_gridx_values, 0, m_nrows);
79  position_data.column("GRIDY", false).read(m_gridy_values, 0, m_nrows);
80  } catch (CCfits::Table::NoSuchColumn) {
81  position_data.column("X_CENTER", false).read(m_gridx_values, 0, m_nrows);
82  position_data.column("Y_CENTER", false).read(m_gridy_values, 0, m_nrows);
83  }
84  position_data.column("RA", false).read(m_ra_values, 0, m_nrows);
85  position_data.column("DEC", false).read(m_dec_values, 0, m_nrows);
86  position_data.column("X", false).read(m_x_values, 0, m_nrows);
87  position_data.column("Y", false).read(m_y_values, 0, m_nrows);
88 
89  } catch (CCfits::FitsException& e) {
90  throw Elements::Exception() << "Error loading stacked PSF file: " << e.message();
91  }
92 }
93 
95  int naxis1, naxis2;
96 
97  // read in the min/max grid values in x/y
98  const auto x_grid_minmax = std::minmax_element(begin(m_gridx_values), end(m_gridx_values));
99  const auto y_grid_minmax = std::minmax_element(begin(m_gridy_values), end(m_gridy_values));
100 
101  // read the image size
102  m_pFits->extension(1).readKey("NAXIS1", naxis1);
103  m_pFits->extension(1).readKey("NAXIS2", naxis2);
104 
105  // make sure all PSF in the grid are there
106  if (*x_grid_minmax.first - m_grid_offset < 1)
107  throw Elements::Exception() << "The PSF at the smallest x-grid starts at: " << *x_grid_minmax.first - m_grid_offset;
108  if (*y_grid_minmax.first - m_grid_offset < 1)
109  throw Elements::Exception() << "The PSF at the smallest y-grid starts at: " << *y_grid_minmax.first - m_grid_offset;
110  if (*x_grid_minmax.second + m_grid_offset > naxis1)
111  throw Elements::Exception() << "The PSF at the largest x-grid is too large: " << *x_grid_minmax.second + m_grid_offset
112  << " NAXIS1: " << naxis1;
113  if (*y_grid_minmax.second + m_grid_offset > naxis2)
114  throw Elements::Exception() << "The PSF at the largest y-grid is too large: " << *y_grid_minmax.second + m_grid_offset
115  << " NAXIS2: " << naxis1;
116 }
117 
119  long index_min_distance = 0;
120  double min_distance = 1.0e+32;
121 
122  // make sure there are only two positions
123  if (values.size() != 2)
124  throw Elements::Exception() << "There can be only two positional value for the stacked PSF!";
125 
126  // find the position of minimal distance
127  for (int act_index = 0; act_index < m_nrows; act_index++) {
128  double act_distance = (values[0] - m_x_values[act_index]) * (values[0] - m_x_values[act_index]) +
129  (values[1] - m_y_values[act_index]) * (values[1] - m_y_values[act_index]);
130  if (act_distance < min_distance) {
131  index_min_distance = act_index;
132  min_distance = act_distance;
133  }
134  }
135  // give some feedback
136  stack_logger.debug() << "Distance: " << sqrt(min_distance) << " (" << values[0] << "," << values[1] << ")<-->("
137  << m_x_values[index_min_distance] << "," << m_y_values[index_min_distance]
138  << ") index: " << index_min_distance;
139 
140  // get the first and last pixels for the PSF to be extracted
141  // NOTE: CCfits has 1-based indices, also the last index is *included* in the reading
142  std::vector<long> first_vertex = {m_gridx_values[index_min_distance] - m_grid_offset,
143  m_gridy_values[index_min_distance] - m_grid_offset};
144  std::vector<long> last_vertex = {first_vertex[0] + m_psf_size - 1, first_vertex[1] + m_psf_size - 1};
145  std::vector<long> stride = {1, 1};
146 
147  // read out the image
148  std::valarray<SeFloat> stamp_data;
149  m_pFits->extension(1).read(stamp_data, first_vertex, last_vertex, stride);
150 
151  // create and return the psf image
152  return VectorImage<SeFloat>::create(m_psf_size, m_psf_size, std::begin(stamp_data), std::end(stamp_data));
153 }
154 
155 } // end SExtractor
static auto stack_logger
T begin(T... args)
static Logging getLogger(const std::string &name="")
std::vector< SeFloat > m_ra_values
virtual std::shared_ptr< VectorImage< SeFloat > > getPsf(const std::vector< double > &values) const
void setup(std::shared_ptr< CCfits::FITS > pFits)
std::shared_ptr< CCfits::FITS > m_pFits
std::vector< SeFloat > m_x_values
std::vector< SeFloat > m_dec_values
std::vector< SeFloat > m_y_values
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
T end(T... args)
T minmax_element(T... args)
constexpr double e
T size(T... args)
T sqrt(T... args)