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 auto size_x = Traits::width(image);
42 auto size_y = Traits::height(image);
43 for (auto& sample : component.getSharpSampling()) {
44 std::size_t image_x = std::get<0>(sample) / pixel_scale + size_x/2.;
45 std::size_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 x=0; x<size_x; ++x) {
58 double x_model = x - (size_x-1) / 2.;
59 x_model *= pixel_scale;
60 for (std::size_t y=0; y<size_y; ++y) {
61 double y_model = y - (size_y-1) / 2.;
62 y_model *= pixel_scale;
63 if (!component.insideSharpRegion(x_model-pixel_scale/2., y_model-pixel_scale/2.) ||
64 !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 Traits::at(image, x, y) = component.getValue(x_model, y_model) * pixel_scale*pixel_scale;
73 } // end of namespace _impl
75 template<typename ImageType>
76 ImageType ExtendedModel<ImageType>::getRasterizedImage(double pixel_scale, std::size_t size_x, std::size_t size_y) const {
77 //std::cout << "]] " << getX() << " " << getY() << "\n";
78 using Traits = ImageTraits<ImageType>;
79 if (size_x % 2 == 0 || size_y % 2 == 0) {
80 throw Elements::Exception() << "Rasterized image dimensions must be odd numbers "
81 << "but got (" << size_x << ',' << size_y << ")";
83 ImageType image = Traits::factory(size_x, size_y);
84 double r_max = std::sqrt(size_x * size_x + size_y * size_y) / 2.;
85 for (auto& component : m_component_list) {
86 component->updateRasterizationInfo(pixel_scale, r_max);
87 ImageType comp_image = Traits::factory(size_x, size_y);
88 _impl::addSharp(comp_image, pixel_scale, *component);
89 _impl::addSmooth(comp_image, pixel_scale, *component);
90 for (auto im_it = Traits::begin(image), comp_it = Traits::begin(comp_image);
91 im_it != Traits::end(image); ++im_it, ++comp_it) {
98 template<typename ImageType>
99 ExtendedModel<ImageType>::ExtendedModel(std::vector<std::unique_ptr<ModelComponent>>&& component_list,
100 std::shared_ptr<BasicParameter> x_scale, std::shared_ptr<BasicParameter> y_scale,
101 std::shared_ptr<BasicParameter> rotation_angle, double width, double height,
102 std::shared_ptr<BasicParameter> x, std::shared_ptr<BasicParameter> y)
103 : PositionedModel(x, y), m_width(width), m_height(height) {
104 for (auto& component : component_list) {
105 auto scaled = Euclid::make_unique<ScaledModelComponent>(std::move(component), x_scale, y_scale);
106 auto rotated = Euclid::make_unique<RotatedModelComponent>(std::move(scaled), rotation_angle);
107 m_component_list.emplace_back(std::move(rotated));
111 template<typename ImageType>
112 double ExtendedModel<ImageType>::getValue(double x, double y) const {
115 return std::accumulate(m_component_list.begin(), m_component_list.end(), 0.,
116 [x, y](double a, const std::unique_ptr<ModelComponent>& b) {
117 return a + b->getValue(x, y);
121 } // end of namespace ModelFitting