SourceXtractorPlusPlus  0.15
Please provide a description of the project.
FlexibleModelFittingModel.cpp
Go to the documentation of this file.
1 
17 /*
18  * FlexibleModelFittingModel.cpp
19  *
20  * Created on: Oct 9, 2018
21  * Author: mschefer
22  */
23 
25 
27 
32 
36 
38 
41 
44 
46 
48 
51 
54 
56 
57 namespace SourceXtractor {
58 
59 using namespace ModelFitting;
61 using namespace Euclid::Configuration;
62 
63 static const double MODEL_MIN_SIZE = 4.0;
64 static const double MODEL_SIZE_FACTOR = 1.2;
65 
66 // Reference for Sersic related quantities:
67 // See https://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html
68 
69 
70 // Note about pixel coordinates:
71 
72 // The model fitting is made in pixel coordinates of the detection image
73 
74 // Internally we use a coordinate system with (0,0) at the center of the first pixel. But for compatibility with
75 // SExtractor 2, all pixel coordinates visible to the end user need to follow the FITS convention of (1,1) being the
76 // center of the coordinate system.
77 
78 // The ModelFitting module uses the more common standard of (0, 0) being the corner of the first pixel.
79 
80 // So we first convert the Python parameter to our internal coordinates, then do the transformation of coordinate,
81 // subtract the offset to the image cut-out and shift the result by 0.5 pixels
82 
83 
85  const SourceInterface& source,
86  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
90  std::shared_ptr<CoordinateSystem> reference_coordinates,
91  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
92 
93  //auto pixel_x = std::make_shared<DependentParameter<std::shared_ptr<BasicParameter>, std::shared_ptr<BasicParameter>>>(
94 
95  auto pixel_x = createDependentParameter(
96  [reference_coordinates, coordinates, offset](double x, double y) {
97  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
98  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
99 
100 
101  auto pixel_y = createDependentParameter(
102  [reference_coordinates, coordinates, offset](double x, double y) {
103  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
104  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
105 
106  point_models.emplace_back(pixel_x, pixel_y, manager.getParameter(source, m_flux));
107 }
108 
110  const SourceInterface& source,
111  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
112  std::vector<ModelFitting::PointModel>& /*point_models*/,
115  std::shared_ptr<CoordinateSystem> reference_coordinates,
116  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
117 
118  auto pixel_x = createDependentParameter(
119  [reference_coordinates, coordinates, offset](double x, double y) {
120  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
121  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
122 
123 
124  auto pixel_y = createDependentParameter(
125  [reference_coordinates, coordinates, offset](double x, double y) {
126  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
127  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
128 
129  //auto n = std::make_shared<ManualParameter>(1); // Sersic index for exponential
130  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
131 
132 // ManualParameter x_scale(1); // we don't scale the x coordinate
133  auto i0 = createDependentParameter(
134  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.35513 * radius * radius * aspect); },
135  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
136  manager.getParameter(source, m_aspect_ratio));
137 
138  auto k = createDependentParameter(
139  [](double eff_radius) { return 1.678 / eff_radius; },
140  manager.getParameter(source, m_effective_radius));
141 
142  auto& boundaries = source.getProperty<PixelBoundaries>();
143  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
144 
146  2.0, i0, k, x_scale, manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_angle),
147  size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
148 }
149 
151  const SourceInterface& source,
152  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
153  std::vector<ModelFitting::PointModel>& /*point_models*/,
156  std::shared_ptr<CoordinateSystem> reference_coordinates,
157  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
158 
159  auto pixel_x = createDependentParameter(
160  [reference_coordinates, coordinates, offset](double x, double y) {
161  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
162  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
163 
164 
165  auto pixel_y = createDependentParameter(
166  [reference_coordinates, coordinates, offset](double x, double y) {
167  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
168  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
169 
170  auto n = std::make_shared<ManualParameter>(4); // Sersic index for Devaucouleurs
171  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
172 
173  auto i0 = createDependentParameter(
174  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.001684925 * radius * radius * aspect); },
175  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
176  manager.getParameter(source, m_aspect_ratio));
177 
178  auto k = createDependentParameter(
179  [](double eff_radius) { return 7.669 / pow(eff_radius, .25); },
180  manager.getParameter(source, m_effective_radius));
181 
182  auto& boundaries = source.getProperty<PixelBoundaries>();
183  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
184 
185  extended_models.emplace_back(std::make_shared<CompactSersicModel<ImageInterfaceTypePtr>>(
186  3.0, i0, k, n, x_scale, manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_angle),
187  size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
188 }
189 
190 static double computeBn(double n) {
191  // Using approximation from MacArthur, L.A., Courteau, S., & Holtzman, J.A. 2003, ApJ, 582, 689
192  return 2 * n - 1.0 / 3.0 + 4 / (405 * n)
193  + 46 / (25515 * n * n) + 131 / (1148175 * n * n * n) - 2194697 / (30690717750 * n * n * n * n);
194 }
195 
197  const SourceInterface& source,
198  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
199  std::vector<ModelFitting::PointModel>& /*point_models*/,
202  std::shared_ptr<CoordinateSystem> reference_coordinates,
203  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
204  auto pixel_x = createDependentParameter(
205  [reference_coordinates, coordinates, offset](double x, double y) {
206  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
207  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
208 
209 
210  auto pixel_y = createDependentParameter(
211  [reference_coordinates, coordinates, offset](double x, double y) {
212  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
213  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
214 
215  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
216 
217  auto i0 = createDependentParameter(
218  [](double flux, double radius, double aspect, double n) { return flux / (2 * M_PI * pow(computeBn(n), -2*n) * n * std::tgamma(2*n) * radius * radius * aspect); },
219  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
220  manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_sersic_index));
221 
222  auto k = createDependentParameter(
223  [](double eff_radius, double n) { return computeBn(n) / pow(eff_radius, 1.0 / n); },
224  manager.getParameter(source, m_effective_radius), manager.getParameter(source, m_sersic_index));
225 
226  auto& boundaries = source.getProperty<PixelBoundaries>();
227  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
228 
229  extended_models.emplace_back(std::make_shared<CompactSersicModel<ImageInterfaceTypePtr>>(
230  3.0, i0, k, manager.getParameter(source, m_sersic_index), x_scale, manager.getParameter(source, m_aspect_ratio),
231  manager.getParameter(source, m_angle), size, size, pixel_x, pixel_y, manager.getParameter(source, m_flux), jacobian));
232 }
233 
235  const SourceInterface& source,
237  std::vector<ModelFitting::PointModel>& /* point_models */,
240  std::shared_ptr<CoordinateSystem> /* reference_coordinates */,
241  std::shared_ptr<CoordinateSystem> /* coordinates */, PixelCoordinate /* offset */) const {
242  constant_models.emplace_back(manager.getParameter(source, m_value));
243 }
244 
245 }
246 
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< std::shared_ptr< ModelFitting::ExtendedModel< ImageInterfaceTypePtr >>> &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive.
The SourceInterface is an abstract "source" that has properties attached to it.
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
T emplace_back(T... args)
T make_shared(T... args)
T max(T... args)
std::unique_ptr< T > make_unique(Args &&... args)
std::shared_ptr< DependentParameter< Parameters... > > createDependentParameter(typename DependentParameter< Parameters... >::ValueCalculator value_calculator, Parameters... parameters)
static const double MODEL_SIZE_FACTOR
static double computeBn(double n)
static const double MODEL_MIN_SIZE
T pow(T... args)
A pixel coordinate made of two integers m_x and m_y.
T tgamma(T... args)