SourceXtractorPlusPlus  0.15
Please provide a description of the project.
VariablePsf.cpp
Go to the documentation of this file.
1 
17 /*
18  * VariablePsf.cpp
19  *
20  * Created on: Jun 25, 2018
21  * Author: Alejandro Álvarez Ayllón
22  */
23 
25 #include <algorithm>
27 
28 
29 namespace SourceXtractor {
30 
31 VariablePsf::VariablePsf(double pixel_sampling, const std::vector<Component> &components,
32  const std::vector<int> &group_degrees,
33  const std::vector<std::shared_ptr<VectorImage<SeFloat>>> &coefficients):
34  m_pixel_sampling(pixel_sampling), m_components(components), m_group_degrees(group_degrees), m_coefficients(coefficients)
35 {
36  selfTest();
39  [](const Component& c) { return c.name; });
40 }
41 
42 VariablePsf::VariablePsf(double pixel_sampling, const std::shared_ptr<VectorImage<SeFloat>> &constant):
43  m_pixel_sampling(pixel_sampling), m_coefficients{constant}
44 {
45  selfTest();
47 }
48 
49 int VariablePsf::getWidth() const {
50  return m_coefficients[0]->getWidth();
51 }
52 
54  return m_coefficients[0]->getHeight();
55 }
56 
58  return m_pixel_sampling;
59 }
60 
62  return m_component_names;
63 }
64 
66 {
67  auto scaled_props = scaleProperties(values);
68 
69  // Initialize with the constant component
71 
72  // Add the rest of the components
73  for (auto i = 1u; i < m_coefficients.size(); ++i) {
74  const auto& exp = m_exponents[i];
75  const auto& coef = m_coefficients[i];
76 
77  double acc = 1.;
78  for (auto j = 0u; j < scaled_props.size(); ++j) {
79  acc *= std::pow(scaled_props[j], exp[j]);
80  }
81 
82  for (auto x = 0; x < getWidth(); ++x) {
83  for (auto y = 0; y < getHeight(); ++y) {
84  result->at(x, y) += acc * coef->at(x, y);
85  }
86  }
87  }
88 
89  return result;
90 }
91 
93  // Pre-condition: There is at least a constant component
94  if (m_coefficients.size() == 0) {
95  throw Elements::Exception() << "A variable PSF needs at least one set of coefficients";
96  }
97 
98  // Pre-condition: There is a degree value per unique group
99  std::vector<int> n_component_per_group(m_group_degrees.size());
100  for (auto& component : m_components) {
101  if (component.group_id >= (int) m_group_degrees.size()) {
102  throw Elements::Exception() << "Component group out of range for " << component.name;
103  }
104  ++n_component_per_group[component.group_id];
105  }
106 
107  // Pre-condition: There are enough coefficients - (n+d)!/(n!d!) per group
108  unsigned int n_coefficients = 1;
109  for (unsigned int g = 0; g < n_component_per_group.size(); ++g) {
110  int dmax = m_group_degrees[g];
111  int n = n_component_per_group[g];
112  int d = std::min<int>(dmax, n);
113  int num, den;
114 
115  for (num = 1, den = 1; d > 0; num *= (n+dmax--), den*= d--);
116 
117  n_coefficients *= num / den;
118  }
119 
120  if (n_coefficients != m_coefficients.size()) {
121  throw Elements::Exception() << "Invalid number of coefficients. Got " << m_coefficients.size()
122  << " expected " << n_coefficients;
123  }
124 
125  // Pre-condition: All components have the same size
126  auto psf_width = m_coefficients[0]->getWidth();
127  auto psf_height = m_coefficients[0]->getHeight();
128 
129  for (auto coeff : m_coefficients) {
130  if (coeff->getWidth() != psf_width || coeff->getHeight() != psf_height) {
131  throw Elements::Exception() << "Malformed variable PSF, coefficient matrices do not have the same dimensions";
132  }
133  }
134 }
135 
137 {
138  if (values.size() != m_components.size()) {
139  throw Elements::Exception()
140  << "Expecting " << m_components.size() << " values, got " << values.size();
141  }
142  std::vector<double> scaled(values.size());
143  for (auto i = 0u; i < values.size(); ++i) {
144  scaled[i] = (values[i] - m_components[i].offset) / m_components[i].scale;
145  }
146  return scaled;
147 }
148 
150  std::vector<int> group_exponents(m_group_degrees);
152 
154  if (m_components.empty()) {
155  return;
156  }
157 
158  // Constant
159  m_exponents[0].resize(m_components.size());
160  --group_exponents[m_components[0].group_id];
161 
162  // Polynomial
163  exp[0] = 1;
164  for (auto e = m_exponents.begin() + 1; e != m_exponents.end(); ++e) {
165  *e = exp;
166 
167  size_t ei = 0;
168  for (auto component : m_components) {
169  if (group_exponents[component.group_id] > 0) {
170  --group_exponents[component.group_id];
171  ++exp[ei];
172  break;
173  }
174  else {
175  group_exponents[component.group_id] = exp[ei];
176  exp[ei++] = 0;
177  }
178  }
179  }
180 }
181 
182 } // end SourceXtractor
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T back_inserter(T... args)
T begin(T... args)
std::shared_ptr< VectorImage< SeFloat > > getPsf(const std::vector< double > &values) const override
Definition: VariablePsf.cpp:65
int getHeight() const override
Definition: VariablePsf.cpp:53
std::vector< std::string > m_component_names
Definition: VariablePsf.h:123
double getPixelSampling() const override
Definition: VariablePsf.cpp:57
VariablePsf(double pixel_sampling, const std::vector< Component > &components, const std::vector< int > &group_degrees, const std::vector< std::shared_ptr< VectorImage< SeFloat >>> &coefficients)
Definition: VariablePsf.cpp:31
std::vector< double > scaleProperties(const std::vector< double > &values) const
Normalizes the values.
std::vector< std::vector< int > > m_exponents
Definition: VariablePsf.h:126
std::vector< std::shared_ptr< VectorImage< SeFloat > > > m_coefficients
Definition: VariablePsf.h:125
int getWidth() const override
Definition: VariablePsf.cpp:49
std::vector< Component > m_components
Definition: VariablePsf.h:122
const std::vector< std::string > & getComponents() const override
Definition: VariablePsf.cpp:61
std::vector< int > m_group_degrees
Definition: VariablePsf.h:124
void selfTest()
Verify that the preconditions of getPsf are met at construction time.
Definition: VariablePsf.cpp:92
Image implementation which keeps the pixel values in memory.
Definition: VectorImage.h:52
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
T end(T... args)
T exp(T... args)
constexpr double e
constexpr double g
T pow(T... args)
T resize(T... args)
T size(T... args)
T transform(T... args)