SourceXtractorPlusPlus  0.16
Please provide a description of the project.
FlexibleModelFittingModel.cpp
Go to the documentation of this file.
1 
17 /*
18  * FlexibleModelFittingModel.cpp
19  *
20  * Created on: Oct 9, 2018
21  * Author: mschefer
22  */
23 
25 
27 
32 
36 
38 
41 
44 
46 
48 
51 
54 
55 #ifdef WITH_ONNX_MODELS
57 #endif
58 
60 
61 namespace SourceXtractor {
62 
63 using namespace ModelFitting;
65 using namespace Euclid::Configuration;
66 
67 static const double MODEL_MIN_SIZE = 4.0;
68 static const double MODEL_SIZE_FACTOR = 1.2;
69 
70 // Reference for Sersic related quantities:
71 // See https://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html
72 
73 
74 // Note about pixel coordinates:
75 
76 // The model fitting is made in pixel coordinates of the detection image
77 
78 // Internally we use a coordinate system with (0,0) at the center of the first pixel. But for compatibility with
79 // SExtractor 2, all pixel coordinates visible to the end user need to follow the FITS convention of (1,1) being the
80 // center of the coordinate system.
81 
82 // The ModelFitting module uses the more common standard of (0, 0) being the corner of the first pixel.
83 
84 // So we first convert the Python parameter to our internal coordinates, then do the transformation of coordinate,
85 // subtract the offset to the image cut-out and shift the result by 0.5 pixels
86 
87 
89  const SourceInterface& source,
90  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
94  std::shared_ptr<CoordinateSystem> reference_coordinates,
95  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
96 
97  //auto pixel_x = std::make_shared<DependentParameter<std::shared_ptr<BasicParameter>, std::shared_ptr<BasicParameter>>>(
98 
99  auto pixel_x = createDependentParameter(
100  [reference_coordinates, coordinates, offset](double x, double y) {
101  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
102  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
103 
104 
105  auto pixel_y = createDependentParameter(
106  [reference_coordinates, coordinates, offset](double x, double y) {
107  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
108  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
109 
110  point_models.emplace_back(pixel_x, pixel_y, manager.getParameter(source, m_flux));
111 }
112 
114  const SourceInterface& source,
115  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
116  std::vector<ModelFitting::PointModel>& /*point_models*/,
119  std::shared_ptr<CoordinateSystem> reference_coordinates,
120  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
121 
122  auto pixel_x = createDependentParameter(
123  [reference_coordinates, coordinates, offset](double x, double y) {
124  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
125  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
126 
127 
128  auto pixel_y = createDependentParameter(
129  [reference_coordinates, coordinates, offset](double x, double y) {
130  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
131  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
132 
133  //auto n = std::make_shared<ManualParameter>(1); // Sersic index for exponential
134  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
135 
136 // ManualParameter x_scale(1); // we don't scale the x coordinate
137  auto i0 = createDependentParameter(
138  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.35513 * radius * radius * aspect); },
139  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
140  manager.getParameter(source, m_aspect_ratio));
141 
142  auto k = createDependentParameter(
143  [](double eff_radius) { return 1.678 / eff_radius; },
144  manager.getParameter(source, m_effective_radius));
145 
146  auto& boundaries = source.getProperty<PixelBoundaries>();
147  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
148 
150  2.0, i0, k, x_scale, manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_angle),
151  size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
152 }
153 
155  const SourceInterface& source,
156  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
157  std::vector<ModelFitting::PointModel>& /*point_models*/,
160  std::shared_ptr<CoordinateSystem> reference_coordinates,
161  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
162 
163  auto pixel_x = createDependentParameter(
164  [reference_coordinates, coordinates, offset](double x, double y) {
165  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
166  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
167 
168 
169  auto pixel_y = createDependentParameter(
170  [reference_coordinates, coordinates, offset](double x, double y) {
171  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
172  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
173 
174  auto n = std::make_shared<ManualParameter>(4); // Sersic index for Devaucouleurs
175  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
176 
177  auto i0 = createDependentParameter(
178  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.001684925 * radius * radius * aspect); },
179  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
180  manager.getParameter(source, m_aspect_ratio));
181 
182  auto k = createDependentParameter(
183  [](double eff_radius) { return 7.669 / pow(eff_radius, .25); },
184  manager.getParameter(source, m_effective_radius));
185 
186  auto& boundaries = source.getProperty<PixelBoundaries>();
187  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
188 
189  extended_models.emplace_back(std::make_shared<CompactSersicModel<ImageInterfaceTypePtr>>(
190  3.0, i0, k, n, x_scale, manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_angle),
191  size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
192 }
193 
194 static double computeBn(double n) {
195  // Using approximation from MacArthur, L.A., Courteau, S., & Holtzman, J.A. 2003, ApJ, 582, 689
196 
197  // The approximation only works for n >= 0.36, so we clamp the value to avoid numerical problems
198  n = std::max(0.36, n);
199 
200  return 2 * n - 1.0 / 3.0 + 4 / (405 * n)
201  + 46 / (25515 * n * n) + 131 / (1148175 * n * n * n) - 2194697 / (30690717750 * n * n * n * n);
202 }
203 
205  const SourceInterface& source,
206  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
207  std::vector<ModelFitting::PointModel>& /*point_models*/,
210  std::shared_ptr<CoordinateSystem> reference_coordinates,
211  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
212  auto pixel_x = createDependentParameter(
213  [reference_coordinates, coordinates, offset](double x, double y) {
214  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
215  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
216 
217  auto pixel_y = createDependentParameter(
218  [reference_coordinates, coordinates, offset](double x, double y) {
219  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
220  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
221 
222  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
223 
224  auto i0 = createDependentParameter(
225  [](double flux, double radius, double aspect, double n) { return flux / (2 * M_PI * pow(computeBn(n), -2*n) * n * std::tgamma(2*n) * radius * radius * aspect); },
226  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
227  manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_sersic_index));
228 
229  auto k = createDependentParameter(
230  [](double eff_radius, double n) { return computeBn(n) / pow(eff_radius, 1.0 / n); },
231  manager.getParameter(source, m_effective_radius), manager.getParameter(source, m_sersic_index));
232 
233  auto& boundaries = source.getProperty<PixelBoundaries>();
234  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
235 
236  extended_models.emplace_back(std::make_shared<CompactSersicModel<ImageInterfaceTypePtr>>(
237  3.0, i0, k, manager.getParameter(source, m_sersic_index), x_scale, manager.getParameter(source, m_aspect_ratio),
238  manager.getParameter(source, m_angle), size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
239 }
240 
242  const SourceInterface& source,
244  std::vector<ModelFitting::PointModel>& /* point_models */,
247  std::shared_ptr<CoordinateSystem> /* reference_coordinates */,
248  std::shared_ptr<CoordinateSystem> /* coordinates */, PixelCoordinate /* offset */) const {
249  constant_models.emplace_back(manager.getParameter(source, m_value));
250 }
251 
252 #ifdef WITH_ONNX_MODELS
253 
254 void FlexibleModelFittingOnnxModel::addForSource(FlexibleModelFittingParameterManager& manager,
255  const SourceInterface& source,
256  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
257  std::vector<ModelFitting::PointModel>& /*point_models*/,
260  std::shared_ptr<CoordinateSystem> reference_coordinates,
261  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
262 
263  auto pixel_x = createDependentParameter(
264  [reference_coordinates, coordinates, offset](double x, double y) {
265  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
266  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
267 
268 
269  auto pixel_y = createDependentParameter(
270  [reference_coordinates, coordinates, offset](double x, double y) {
271  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
272  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
273 
274  auto y_scale = createDependentParameter(
275  [](double scale, double ratio) {
276  return scale * ratio;
277  }, manager.getParameter(source, m_scale), manager.getParameter(source, m_aspect_ratio));
278 
279  auto& boundaries = source.getProperty<PixelBoundaries>();
280  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
281 
283  for (auto it : m_params) {
284  params[it.first] = manager.getParameter(source, it.second);
285  }
286 
287  extended_models.emplace_back(std::make_shared<OnnxCompactModel<ImageInterfaceTypePtr>>(m_models,
288  manager.getParameter(source, m_scale), y_scale, manager.getParameter(source, m_angle),
289  size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), params, jacobian));
290 }
291 
292 FlexibleModelFittingOnnxModel::FlexibleModelFittingOnnxModel(
301  : m_x(x),
302  m_y(y),
303  m_flux(flux),
304  m_aspect_ratio(aspect_ratio),
305  m_angle(angle),
306  m_scale(scale),
307  m_params(params),
308  m_models(models){
309 
310  std::sort(m_models.begin(), m_models.end(),
311  [](const std::shared_ptr<OnnxModel>& a, const std::shared_ptr<OnnxModel>& b) -> bool {
312  return a->getOutputShape()[2] < b->getOutputShape()[2];
313  });
314 }
315 
316 #endif
317 
318 }
319 
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive.
The SourceInterface is an abstract "source" that has properties attached to it.
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
T emplace_back(T... args)
T make_shared(T... args)
T max(T... args)
std::unique_ptr< T > make_unique(Args &&... args)
std::shared_ptr< DependentParameter< Parameters... > > createDependentParameter(typename DependentParameter< Parameters... >::ValueCalculator value_calculator, Parameters... parameters)
static const double MODEL_SIZE_FACTOR
static double computeBn(double n)
static const double MODEL_MIN_SIZE
T pow(T... args)
T sort(T... args)
A pixel coordinate made of two integers m_x and m_y.
T tgamma(T... args)