CompactSersicModel.icpp

Go to the documentation of this file.

/*
 * CompactSersicModel.icpp
 *
 *  Created on: Jul 25, 2019
 *      Author: mschefer
 */
 
#include <math.h>
 
namespace ModelFitting {
 
template<typename ImageType>
CompactSersicModel<ImageType>::CompactSersicModel(double sharp_radius,
    std::shared_ptr<BasicParameter> i0, std::shared_ptr<BasicParameter> k, std::shared_ptr<BasicParameter> n,
    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::shared_ptr<BasicParameter> flux,
    std::tuple<double, double, double, double> transform
)
        : CompactModelBase<ImageType>(x_scale, y_scale, rotation, width, height, x, y, transform),
          m_sharp_radius_squared(float(sharp_radius * sharp_radius)),
          m_i0(i0), m_k(k), m_n(n), m_flux(flux)
{}
 
template<typename ImageType>
double CompactSersicModel<ImageType>::getValue(double x, double y) const {
  SersicModelEvaluator model_eval;
  model_eval.transform = getCombinedTransform(1.0);
  model_eval.i0 = m_i0->getValue();
  model_eval.k = m_k->getValue();
  model_eval.n = m_n->getValue();
  model_eval.max_r_sqr = 1e30;
 
  auto area_correction = (1.0 / fabs(m_jacobian[0] * m_jacobian[3] - m_jacobian[1] * m_jacobian[2]));
  return model_eval.evaluateModel(x, y) * area_correction;
}
 
 
template<typename ImageType>
ImageType CompactSersicModel<ImageType>::getRasterizedImage(double pixel_scale, std::size_t size_x, std::size_t size_y) const {
  using Traits = ImageTraits<ImageType>;
 
  if (size_x % 2 == 0 || size_y % 2 == 0) {
    throw Elements::Exception() << "Rasterized image dimensions must be odd numbers "
        << "but got (" << size_x << ',' << size_y << ")";
  }
 
  // Add a 4 pixels border to prevent cropping from Lanczos resize
  ImageType image = Traits::factory(size_x+8, size_y+8);
 
  auto combined_tranform = getCombinedTransform(pixel_scale);
 
  SersicModelEvaluator model_eval;
  model_eval.transform = combined_tranform;
  model_eval.i0 = m_i0->getValue();
  model_eval.k = m_k->getValue();
  model_eval.n = m_n->getValue();
  model_eval.max_r_sqr = getMaxRadiusSqr(size_x, size_y, combined_tranform);
 
  float area_correction = (1.0 / fabs(m_jacobian[0] * m_jacobian[3] - m_jacobian[1] * m_jacobian[2])) * pixel_scale * pixel_scale;
 
  for (int y = 0; y < (int)size_y; ++y) {
    int dy = y - size_y / 2;
    for (int x = 0; x < (int)size_x; ++x) {
      int dx = x - size_x / 2;
      if (dx * dx + dy * dy < m_sharp_radius_squared) {
        Traits::at(image, x+4, y+4) = adaptiveSamplePixel(model_eval, dx, dy, 7, 0.01) * area_correction;
      } else {
        Traits::at(image, x+4, y+4) = model_eval.evaluateModel(dx, dy) * area_correction;
      }
    }
  }
 
  renormalize(image, m_flux->getValue());
 
  return image;
}
 
}