SourceXtractorPlusPlus  0.15
Please provide a description of the project.
CompactSersicModel.icpp
Go to the documentation of this file.
1 /*
2  * CompactSersicModel.icpp
3  *
4  * Created on: Jul 25, 2019
5  * Author: mschefer
6  */
7 
8 #include <math.h>
9 
10 namespace ModelFitting {
11 
12 template<typename ImageType>
13 CompactSersicModel<ImageType>::CompactSersicModel(double sharp_radius,
14  std::shared_ptr<BasicParameter> i0, std::shared_ptr<BasicParameter> k, std::shared_ptr<BasicParameter> n,
15  std::shared_ptr<BasicParameter> x_scale, std::shared_ptr<BasicParameter> y_scale,
16  std::shared_ptr<BasicParameter> rotation, double width, double height,
17  std::shared_ptr<BasicParameter> x, std::shared_ptr<BasicParameter> y,
18  std::shared_ptr<BasicParameter> flux,
19  std::tuple<double, double, double, double> transform
20 )
21  : CompactModelBase<ImageType>(x_scale, y_scale, rotation, width, height, x, y, transform),
22  m_sharp_radius_squared(float(sharp_radius * sharp_radius)),
23  m_i0(i0), m_k(k), m_n(n), m_flux(flux)
24 {}
25 
26 template<typename ImageType>
27 double CompactSersicModel<ImageType>::getValue(double x, double y) const {
28  SersicModelEvaluator model_eval;
29  model_eval.transform = getCombinedTransform(1.0);
30  model_eval.i0 = m_i0->getValue();
31  model_eval.k = m_k->getValue();
32  model_eval.n = m_n->getValue();
33  model_eval.max_r_sqr = 1e30;
34 
35  auto area_correction = (1.0 / fabs(m_jacobian[0] * m_jacobian[3] - m_jacobian[1] * m_jacobian[2]));
36  return model_eval.evaluateModel(x, y) * area_correction;
37 }
38 
39 
40 template<typename ImageType>
41 ImageType CompactSersicModel<ImageType>::getRasterizedImage(double pixel_scale, std::size_t size_x, std::size_t size_y) const {
42  using Traits = ImageTraits<ImageType>;
43 
44  if (size_x % 2 == 0 || size_y % 2 == 0) {
45  throw Elements::Exception() << "Rasterized image dimensions must be odd numbers "
46  << "but got (" << size_x << ',' << size_y << ")";
47  }
48 
49  ImageType image = Traits::factory(size_x, size_y);
50 
51  auto combined_tranform = getCombinedTransform(pixel_scale);
52 
53  SersicModelEvaluator model_eval;
54  model_eval.transform = combined_tranform;
55  model_eval.i0 = m_i0->getValue();
56  model_eval.k = m_k->getValue();
57  model_eval.n = m_n->getValue();
58  model_eval.max_r_sqr = getMaxRadiusSqr(size_x, size_y, combined_tranform);
59 
60  float area_correction = (1.0 / fabs(m_jacobian[0] * m_jacobian[3] - m_jacobian[1] * m_jacobian[2])) * pixel_scale * pixel_scale;
61 
62 
63  for (int x=0; x<(int)size_x; ++x) {
64  int dx = x - size_x / 2;
65  for (int y=0; y<(int)size_y; ++y) {
66  int dy = y - size_y / 2;
67  if (dx*dx + dy*dy < m_sharp_radius_squared) {
68  Traits::at(image, x, y) = adaptiveSamplePixel(model_eval, dx, dy, 7, 0.01) * area_correction;
69  } else {
70  Traits::at(image, x, y) = model_eval.evaluateModel(dx, dy) * area_correction;
71  }
72  }
73  }
74 
75  renormalize(image, m_flux->getValue());
76 
77  return image;
78 }
79 
80 }
81