ExtendedModel.icpp

Go to the documentation of this file.

#include <iostream>
#include <cmath> // for std::sqrt

#include "ElementsKernel/Exception.h"

#include "AlexandriaKernel/memory_tools.h"

#include "ModelFitting/Image/ImageTraits.h"

using Euclid::make_unique;

namespace ModelFitting {

namespace _impl {
  
  template <typename ImageType>
  void addSharp(ImageType& image, double pixel_scale, ModelComponent& component) {
    using Traits = ImageTraits<ImageType>;
    ssize_t size_x = Traits::width(image);
    ssize_t size_y = Traits::height(image);
    for (auto& sample : component.getSharpSampling()) {
      ssize_t image_x = std::get<0>(sample) / pixel_scale + size_x / 2.;
      ssize_t image_y = std::get<1>(sample) / pixel_scale + size_y / 2.;
      if (image_x >= 0 && image_x < size_x && image_y >= 0 && image_y < size_y) {
        Traits::at(image, image_x, image_y) += std::get<2>(sample);
      }
    }
  }
  
  template <typename ImageType>
  void addSmooth(ImageType& image, double pixel_scale, ModelComponent& component) {
    using Traits = ImageTraits<ImageType>;
    auto size_x = Traits::width(image);
    auto size_y = Traits::height(image);
    for (std::size_t y = 0; y < size_y; ++y) {
      double y_model = y - (size_y - 1) / 2.;
      y_model *= pixel_scale;
      for (std::size_t x = 0; x < size_x; ++x) {
        double x_model = x - (size_x - 1) / 2.;
        x_model *= pixel_scale;

        if (!component.insideSharpRegion(x_model - pixel_scale / 2., y_model - pixel_scale / 2.) ||
            !component.insideSharpRegion(x_model - pixel_scale / 2., y_model + pixel_scale / 2.) ||
            !component.insideSharpRegion(x_model + pixel_scale / 2., y_model - pixel_scale / 2.) ||
            !component.insideSharpRegion(x_model + pixel_scale / 2., y_model + pixel_scale / 2.)) {
          Traits::at(image, x, y) = component.getValue(x_model, y_model) * pixel_scale * pixel_scale;
        }
      }
    }
  }

} // end of namespace _impl

template<typename ImageType>
ImageType ExtendedModel<ImageType>::getRasterizedImage(double pixel_scale, std::size_t size_x, std::size_t size_y) const {
  //std::cout << "]] " << getX() << " " << getY() << "\n";
  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 << ")";
  }
  ImageType image = Traits::factory(size_x, size_y);
  double r_max = std::sqrt(size_x * size_x + size_y * size_y) / 2.;
  for (auto& component : m_component_list) {
    component->updateRasterizationInfo(pixel_scale, r_max);
    ImageType comp_image = Traits::factory(size_x, size_y);
    _impl::addSharp(comp_image, pixel_scale, *component);
    _impl::addSmooth(comp_image, pixel_scale, *component);
    for (auto im_it = Traits::begin(image), comp_it = Traits::begin(comp_image);
        im_it != Traits::end(image); ++im_it, ++comp_it) {
      *im_it += *comp_it;
    }
  }
  return image;
}

template<typename ImageType>
ExtendedModel<ImageType>::ExtendedModel(std::vector<std::unique_ptr<ModelComponent>>&& component_list,
    std::shared_ptr<BasicParameter> x_scale, std::shared_ptr<BasicParameter> y_scale,
    std::shared_ptr<BasicParameter> rotation_angle, double width, double height,
    std::shared_ptr<BasicParameter> x, std::shared_ptr<BasicParameter> y)
        : PositionedModel(x, y), m_width(width), m_height(height) {
  for (auto& component : component_list) {
    auto scaled = Euclid::make_unique<ScaledModelComponent>(std::move(component), x_scale, y_scale);
    auto rotated = Euclid::make_unique<RotatedModelComponent>(std::move(scaled), rotation_angle);
    m_component_list.emplace_back(std::move(rotated));
  }
}

template<typename ImageType>
double ExtendedModel<ImageType>::getValue(double x, double y) const {
  x -= getX();
  y -= getY();
  return std::accumulate(m_component_list.begin(), m_component_list.end(), 0.,
                         [x, y](double a, const std::unique_ptr<ModelComponent>& b) {
                           return a + b->getValue(x, y);
                         });
}

} // end of namespace ModelFitting