CompactModelBase.icpp

Go to the documentation of this file.

/*
 * CompactModelBase.icpp
 *
 *  Created on: Aug 19, 2019
 *      Author: mschefer
 */


#include <random>

namespace ModelFitting {

template<typename ImageType>
CompactModelBase<ImageType>::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)
        : ExtendedModel<ImageType>({}, x_scale, y_scale, rotation, width, height, x, y),
          m_x_scale(x_scale), m_y_scale(y_scale), m_rotation(rotation)
{
  m_jacobian = Mat22(transform).GetTranspose();
  m_inv_jacobian = m_jacobian.GetInverse();
}

template<typename ImageType>
Mat22 CompactModelBase<ImageType>::getCombinedTransform(double pixel_scale) const {
  double s, c;
#if defined(__APPLE__)
  __sincos(m_rotation->getValue(), &s, &c);
#elif defined(_GNU_SOURCE)
  sincos(m_rotation->getValue(), &s, &c);
#else
  s = sin(m_rotation->getValue());
  c = cos(m_rotation->getValue());
#endif

  Mat22 rotation(
      c, s,
      -s, c);

  Mat22 scale(
      1. / m_x_scale->getValue(), 0.0,
      0.0, 1. / m_y_scale->getValue());

  return scale * rotation * m_inv_jacobian * pixel_scale;
}

template<typename ImageType>
template<typename ModelEvaluator>
inline float CompactModelBase<ImageType>::samplePixel(const ModelEvaluator& model_eval, int x, int y, unsigned int subsampling) const {
  double acc = 0.;
  for (std::size_t ix=0; ix<subsampling; ++ix) {
    float x_model = (x - 0.5 + (ix+1) * 1.0 / (subsampling+1));
    for (std::size_t iy=0; iy<subsampling; ++iy) {
      float y_model = (y - 0.5 + (iy+1) * 1.0 / (subsampling+1));
      acc += model_eval.evaluateModel(x_model, y_model);
    }
  }

  return acc / (subsampling*subsampling);
}

template<typename ImageType>
template<typename ModelEvaluator>
inline float CompactModelBase<ImageType>::sampleStochastic(const ModelEvaluator& model_eval, int x, int y, unsigned int samples) const {
  double acc = 0.;
  std::default_random_engine generator;
  std::uniform_real_distribution<double> distribution(-.5,.5);

  for (unsigned int i=0; i<samples; i++) {
    acc += model_eval.evaluateModel(double(x) + distribution(generator), double(y) + distribution(generator));
  }

  return acc / samples;
}



template<typename ImageType>
template<typename ModelEvaluator>
inline float CompactModelBase<ImageType>::adaptiveSamplePixel(const ModelEvaluator& model_eval, int x, int y, unsigned int max_subsampling, float threshold) const {
  unsigned int steps[] = {1,3,5,7,11,15,23,31,47,63,95,127};
  float value = samplePixel(model_eval, x,y, 1);
  for (unsigned int i=2; i < (sizeof(steps)/sizeof(steps[0])) && steps[i] <= max_subsampling; i++) {
    float newValue = samplePixel(model_eval, x,y, steps[i] + (max_subsampling % 2));

    double diff = fabs(newValue - value);
    if (diff <= threshold * value) {
      value = newValue;
      break;
    }

    value = newValue;
  }

  return value;
}

// computes the square of distance from origin to a line defined by 2 points
template<typename ImageType>
double CompactModelBase<ImageType>::computeSqrDistanceLineToOrigin(double x1, double y1, double x2, double y2) const {
  double a = (x2-x1) * y1 - (y2-y1) * x1;
  return a * a / ((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
}

template<typename ImageType>
double CompactModelBase<ImageType>::getMaxRadiusSqr(std::size_t size_x, std::size_t size_y, const Mat22& transform) const {

  // transform coordinates of 3 corner points of the rendering rectangle, to be in the same space the model is evaluated

  double x = size_x * 0.95f / 2.f;
  double y = size_y * 0.95f / 2.f;

  double p0_x =  x * transform[0] + y * transform[1];
  double p0_y =  x * transform[2] + y * transform[3];
  double p1_x = -x * transform[0] + y * transform[1];
  double p1_y = -x * transform[2] + y * transform[3];
  double p2_x =  x * transform[0] - y * transform[1];
  double p2_y =  x * transform[2] - y * transform[3];

  // Now computes distance from border lines in that space and use that as maximum radius to be evaluated

  return std::min(computeSqrDistanceLineToOrigin(p0_x, p0_y, p1_x, p1_y),
                  computeSqrDistanceLineToOrigin(p0_x, p0_y, p2_x, p2_y));
}

template<typename ImageType>
void CompactModelBase<ImageType>::renormalize(ImageType& image, double flux) const {
  using Traits = ImageTraits<ImageType>;

  int width = Traits::width(image);
  int height = Traits::height(image);

  double acc = 0.0;

  for (int y=0; y<height; y++) {
    for (int x=0; x<width; x++) {
      acc += Traits::at(image, x, y);
    }
  }

  if (acc > 0.0) {
    double scale = flux / acc;
    for (int y=0; y<height; y++) {
      for (int x=0; x<width; x++) {
        Traits::at(image, x, y) = Traits::at(image, x, y) * scale;
      }
    }
  }
}

}