SourceXtractorPlusPlus  0.10
Please provide a description of the project.
ExtendedModel.icpp
Go to the documentation of this file.
1 /** Copyright © 2019 Université de Genève, LMU Munich - Faculty of Physics, IAP-CNRS/Sorbonne Université
2  *
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)
6  * any later version.
7  *
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
11  * details.
12  *
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
16  */
17 /**
18  * @file ExtendedModel.icpp
19  * @date September 1, 2015
20  * @author Nikolaos Apostolakos
21  */
22 
23 #include <iostream>
24 #include <cmath> // for std::sqrt
25 
26 #include "ElementsKernel/Exception.h"
27 
28 #include "AlexandriaKernel/memory_tools.h"
29 
30 #include "ModelFitting/Image/ImageTraits.h"
31 
32 using Euclid::make_unique;
33 
34 namespace ModelFitting {
35 
36 namespace _impl {
37 
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);
48  }
49  }
50  }
51 
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;
68  }
69  }
70  }
71  }
72 
73 } // end of namespace _impl
74 
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 << ")";
82  }
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) {
92  *im_it += *comp_it;
93  }
94  }
95  return image;
96 }
97 
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 = make_unique<ScaledModelComponent>(std::move(component), x_scale, y_scale);
106  auto rotated = make_unique<RotatedModelComponent>(std::move(scaled), rotation_angle);
107  m_component_list.emplace_back(std::move(rotated));
108  }
109 }
110 
111 template<typename ImageType>
112 double ExtendedModel<ImageType>::getValue(double x, double y) const {
113  x -= getX();
114  y -= getY();
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);
118  });
119 }
120 
121 } // end of namespace ModelFitting