1 /** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
3 * This library is free software; you can redistribute it and/or modify it under
4 * the terms of the GNU Lesser General Public License as published by the Free
5 * Software Foundation; either version 3.0 of the License, or (at your option)
8 * This library is distributed in the hope that it will be useful, but WITHOUT
9 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
10 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 * You should have received a copy of the GNU Lesser General Public License
14 * along with this library; if not, write to the Free Software Foundation, Inc.,
15 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 * @file ExtendedModel.icpp
19 * @date September 1, 2015
20 * @author Nikolaos Apostolakos
24 #include <cmath> // for std::sqrt
26 #include "ElementsKernel/Exception.h"
28 #include "AlexandriaKernel/memory_tools.h"
30 #include "ModelFitting/Image/ImageTraits.h"
32 using Euclid::make_unique;
34 namespace ModelFitting {
38 template <typename ImageType>
39 void addSharp(ImageType& image, double pixel_scale, ModelComponent& component) {
40 using Traits = ImageTraits<ImageType>;
41 ssize_t size_x = Traits::width(image);
42 ssize_t size_y = Traits::height(image);
43 for (auto& sample : component.getSharpSampling()) {
44 ssize_t image_x = std::get<0>(sample) / pixel_scale + size_x / 2.;
45 ssize_t image_y = std::get<1>(sample) / pixel_scale + size_y / 2.;
46 if (image_x >= 0 && image_x < size_x && image_y >= 0 && image_y < size_y) {
47 Traits::at(image, image_x, image_y) += std::get<2>(sample);
52 template <typename ImageType>
53 void addSmooth(ImageType& image, double pixel_scale, ModelComponent& component) {
54 using Traits = ImageTraits<ImageType>;
55 auto size_x = Traits::width(image);
56 auto size_y = Traits::height(image);
57 for (std::size_t y = 0; y < size_y; ++y) {
58 double y_model = y - (size_y - 1) / 2.;
59 y_model *= pixel_scale;
60 for (std::size_t x = 0; x < size_x; ++x) {
61 double x_model = x - (size_x - 1) / 2.;
62 x_model *= pixel_scale;
64 if (!component.insideSharpRegion(x_model - pixel_scale / 2., y_model - pixel_scale / 2.) ||
65 !component.insideSharpRegion(x_model - pixel_scale / 2., y_model + pixel_scale / 2.) ||
66 !component.insideSharpRegion(x_model + pixel_scale / 2., y_model - pixel_scale / 2.) ||
67 !component.insideSharpRegion(x_model + pixel_scale / 2., y_model + pixel_scale / 2.)) {
68 Traits::at(image, x, y) = component.getValue(x_model, y_model) * pixel_scale * pixel_scale;
74 } // end of namespace _impl
76 template<typename ImageType>
77 ImageType ExtendedModel<ImageType>::getRasterizedImage(double pixel_scale, std::size_t size_x, std::size_t size_y) const {
78 //std::cout << "]] " << getX() << " " << getY() << "\n";
79 using Traits = ImageTraits<ImageType>;
80 if (size_x % 2 == 0 || size_y % 2 == 0) {
81 throw Elements::Exception() << "Rasterized image dimensions must be odd numbers "
82 << "but got (" << size_x << ',' << size_y << ")";
84 ImageType image = Traits::factory(size_x, size_y);
85 double r_max = std::sqrt(size_x * size_x + size_y * size_y) / 2.;
86 for (auto& component : m_component_list) {
87 component->updateRasterizationInfo(pixel_scale, r_max);
88 ImageType comp_image = Traits::factory(size_x, size_y);
89 _impl::addSharp(comp_image, pixel_scale, *component);
90 _impl::addSmooth(comp_image, pixel_scale, *component);
91 for (auto im_it = Traits::begin(image), comp_it = Traits::begin(comp_image);
92 im_it != Traits::end(image); ++im_it, ++comp_it) {
99 template<typename ImageType>
100 ExtendedModel<ImageType>::ExtendedModel(std::vector<std::unique_ptr<ModelComponent>>&& component_list,
101 std::shared_ptr<BasicParameter> x_scale, std::shared_ptr<BasicParameter> y_scale,
102 std::shared_ptr<BasicParameter> rotation_angle, double width, double height,
103 std::shared_ptr<BasicParameter> x, std::shared_ptr<BasicParameter> y)
104 : PositionedModel(x, y), m_width(width), m_height(height) {
105 for (auto& component : component_list) {
106 auto scaled = Euclid::make_unique<ScaledModelComponent>(std::move(component), x_scale, y_scale);
107 auto rotated = Euclid::make_unique<RotatedModelComponent>(std::move(scaled), rotation_angle);
108 m_component_list.emplace_back(std::move(rotated));
112 template<typename ImageType>
113 double ExtendedModel<ImageType>::getValue(double x, double y) const {
116 return std::accumulate(m_component_list.begin(), m_component_list.end(), 0.,
117 [x, y](double a, const std::unique_ptr<ModelComponent>& b) {
118 return a + b->getValue(x, y);
122 } // end of namespace ModelFitting