SourceXtractorPlusPlus  0.14
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CompactModelBase.icpp
Go to the documentation of this file.
1 /*
2  * CompactModelBase.icpp
3  *
4  * Created on: Aug 19, 2019
5  * Author: mschefer
6  */
7 
8 
9 #include <random>
10 
11 namespace ModelFitting {
12 
13 template<typename ImageType>
19  : ExtendedModel<ImageType>({}, x_scale, y_scale, rotation, width, height, x, y),
21 {
22  m_jacobian = Mat22(transform).GetTranspose();
23  m_inv_jacobian = m_jacobian.GetInverse();
24 }
25 
26 template<typename ImageType>
28  double s, c;
29 #ifdef _GNU_SOURCE
30  sincos(m_rotation->getValue(), &s, &c);
31 #else
32  s = sin(m_rotation->getValue());
33  c = cos(m_rotation->getValue());
34 #endif
35 
37  c, s,
38  -s, c);
39 
40  Mat22 scale(
41  1. / m_x_scale->getValue(), 0.0,
42  0.0, 1. / m_y_scale->getValue());
43 
44  return scale * rotation * m_inv_jacobian * pixel_scale;
45 }
46 
47 template<typename ImageType>
48 template<typename ModelEvaluator>
49 inline float CompactModelBase<ImageType>::samplePixel(const ModelEvaluator& model_eval, int x, int y, unsigned int subsampling) const {
50  double acc = 0.;
51  for (std::size_t ix=0; ix<subsampling; ++ix) {
52  float x_model = (x - 0.5 + (ix+1) * 1.0 / (subsampling+1));
53  for (std::size_t iy=0; iy<subsampling; ++iy) {
54  float y_model = (y - 0.5 + (iy+1) * 1.0 / (subsampling+1));
55  acc += model_eval.evaluateModel(x_model, y_model);
56  }
57  }
58 
59  return acc / (subsampling*subsampling);
60 }
61 
62 template<typename ImageType>
63 template<typename ModelEvaluator>
64 inline float CompactModelBase<ImageType>::sampleStochastic(const ModelEvaluator& model_eval, int x, int y, unsigned int samples) const {
65  double acc = 0.;
67  std::uniform_real_distribution<double> distribution(-.5,.5);
68 
69  for (unsigned int i=0; i<samples; i++) {
70  acc += model_eval.evaluateModel(double(x) + distribution(generator), double(y) + distribution(generator));
71  }
72 
73  return acc / samples;
74 }
75 
76 
77 
78 template<typename ImageType>
79 template<typename ModelEvaluator>
80 inline float CompactModelBase<ImageType>::adaptiveSamplePixel(const ModelEvaluator& model_eval, int x, int y, unsigned int max_subsampling, float threshold) const {
81  unsigned int steps[] = {1,3,5,7,11,15,23,31,47,63,95,127};
82  float value = samplePixel(model_eval, x,y, 1);
83  for (unsigned int i=2; i < (sizeof(steps)/sizeof(steps[0])) && steps[i] <= max_subsampling; i++) {
84  float newValue = samplePixel(model_eval, x,y, steps[i] + (max_subsampling % 2));
85 
86  double diff = fabs(newValue - value);
87  if (diff <= threshold * value) {
88  value = newValue;
89  break;
90  }
91 
92  value = newValue;
93  }
94 
95  return value;
96 }
97 
98 template<typename ImageType>
100 
101  float x_x = (size_x * 0.95f / 2.f - 1) * transform[0];
102  float x_y = (size_x * 0.95f / 2.f - 1) * transform[2];
103 
104  float y_x = (size_y * 0.95f / 2.f - 1) * transform[1];
105  float y_y = (size_y * 0.95f / 2.f - 1) * transform[3];
106 
107  float xy_x = x_x + y_x;
108  float xy_y = x_y + y_y;
109 
110  float d1 = (xy_x * x_y - xy_y * x_x) * (xy_x * x_y - xy_y * x_x) / ((xy_x - x_x) * (xy_x - x_x) + (xy_y - x_y) * (xy_y - x_y));
111  float d2 = (xy_x * y_y - xy_y * y_x) * (xy_x * y_y - xy_y * y_x) / ((xy_x - y_x) * (xy_x - y_x) + (xy_y - y_y) * (xy_y - y_y));
112 
113  return std::min(d1, d2);
114 }
115 
116 template<typename ImageType>
117 void CompactModelBase<ImageType>::renormalize(ImageType& image, double flux) const {
119 
120  int width = Traits::width(image);
121  int height = Traits::height(image);
122 
123  double acc = 0.0;
124 
125  for (int y=0; y<height; y++) {
126  for (int x=0; x<width; x++) {
127  acc += Traits::at(image, x, y);
128  }
129  }
130 
131  if (acc > 0.0) {
132  double scale = flux / acc;
133  for (int y=0; y<height; y++) {
134  for (int x=0; x<width; x++) {
135  Traits::at(image, x, y) = Traits::at(image, x, y) * scale;
136  }
137  }
138  }
139 }
140 
141 }
float samplePixel(const ModelEvaluator &model_eval, int x, int y, unsigned int subsampling) const
constexpr double s
Mat22 GetTranspose() const
Definition: Mat22.h:73
double getMaxRadiusSqr(std::size_t size_x, std::size_t size_y, const Mat22 &transform) const
T min(T...args)
T sin(T...args)
m_rotation(rotation)
Mat22 getCombinedTransform(double pixel_scale) const
T cos(T...args)
float adaptiveSamplePixel(const ModelEvaluator &model_eval, int x, int y, unsigned int max_subsampling, float threshold=1.1) const
T fabs(T...args)
m_y_scale(y_scale)
void renormalize(ImageType &image, double flux) const
ModelFitting::ImageTraits< ImageInterfaceTypePtr > Traits
Mat22 GetInverse() const
Definition: Mat22.h:58
float sampleStochastic(const ModelEvaluator &model_eval, int x, int y, unsigned int samples=100) const
m_x_scale(x_scale)
const double pixel_scale
Definition: TestImage.cpp:75
std::pair< double, double > transform(int x, int y, const std::array< double, 4 > &t)
CompactModelBase(std::shared_ptr< BasicParameter > x_scale, std::shared_ptr< BasicParameter > y_scale, std::shared_ptr< BasicParameter > rotation, double width, double height, std::shared_ptr< BasicParameter > x, std::shared_ptr< BasicParameter > y, std::tuple< double, double, double, double > transform)