2 * CompactModelBase.icpp
4 * Created on: Aug 19, 2019
11 namespace ModelFitting {
13 template<typename ImageType>
14 CompactModelBase<ImageType>::CompactModelBase(
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::tuple<double, double, double, double> transform)
19 : ExtendedModel<ImageType>({}, x_scale, y_scale, rotation, width, height, x, y),
20 m_x_scale(x_scale), m_y_scale(y_scale), m_rotation(rotation)
22 m_jacobian = Mat22(transform).GetTranspose();
23 m_inv_jacobian = m_jacobian.GetInverse();
26 template<typename ImageType>
27 Mat22 CompactModelBase<ImageType>::getCombinedTransform(double pixel_scale) const {
30 sincos(m_rotation->getValue(), &s, &c);
32 s = sin(m_rotation->getValue());
33 c = cos(m_rotation->getValue());
41 1. / m_x_scale->getValue(), 0.0,
42 0.0, 1. / m_y_scale->getValue());
44 return scale * rotation * m_inv_jacobian * pixel_scale;
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 {
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);
59 return acc / (subsampling*subsampling);
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 {
66 std::default_random_engine generator;
67 std::uniform_real_distribution<double> distribution(-.5,.5);
69 for (unsigned int i=0; i<samples; i++) {
70 acc += model_eval.evaluateModel(double(x) + distribution(generator), double(y) + distribution(generator));
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));
86 double diff = fabs(newValue - value);
87 if (diff <= threshold * value) {
98 template<typename ImageType>
99 double CompactModelBase<ImageType>::getMaxRadiusSqr(std::size_t size_x, std::size_t size_y, const Mat22& transform) const {
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];
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];
107 float xy_x = x_x + y_x;
108 float xy_y = x_y + y_y;
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));
113 return std::min(d1, d2);
116 template<typename ImageType>
117 void CompactModelBase<ImageType>::renormalize(ImageType& image, double flux) const {
118 using Traits = ImageTraits<ImageType>;
120 int width = Traits::width(image);
121 int height = Traits::height(image);
125 for (int y=0; y<height; y++) {
126 for (int x=0; x<width; x++) {
127 acc += Traits::at(image, x, y);
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;