SourceXtractorPlusPlus  0.8
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
24 #include "ModelFitting/utils.h"
25 
27 
32 
36 
38 
41 
45 
47 
48 namespace SourceXtractor {
49 
50 using namespace ModelFitting;
51 
52 static const double MODEL_MIN_SIZE = 4.0;
53 static const double MODEL_SIZE_FACTOR = 1.2;
54 
55 // Reference for Sersic related quantities:
56 // See https://ned.ipac.caltech.edu/level5/March05/Graham/Graham2.html
57 
58 
59 // Note about pixel coordinates:
60 
61 // The model fitting is made in pixel coordinates of the detection image
62 
63 // Internally we use a coordinate system with (0,0) at the center of the first pixel. But for compatibility with
64 // SExtractor 2, all pixel coordinates visible to the end user need to follow the FITS convention of (1,1) being the
65 // center of the coordinate system.
66 
67 // The ModelFitting module uses the more common standard of (0, 0) being the corner of the first pixel.
68 
69 // So we first convert the Python parameter to our internal coordinates, then do the transformation of coordinate,
70 // subtract the offset to the image cut-out and shift the result by 0.5 pixels
71 
72 
74  const SourceInterface& source,
75  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
79  std::shared_ptr<CoordinateSystem> reference_coordinates,
80  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
81 
82  //auto pixel_x = std::make_shared<DependentParameter<std::shared_ptr<BasicParameter>, std::shared_ptr<BasicParameter>>>(
83 
84  auto pixel_x = createDependentParameter(
85  [reference_coordinates, coordinates, offset](double x, double y) {
86  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
87  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
88 
89 
90  auto pixel_y = createDependentParameter(
91  [reference_coordinates, coordinates, offset](double x, double y) {
92  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
93  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
94 
95  point_models.emplace_back(pixel_x, pixel_y, manager.getParameter(source, m_flux));
96 }
97 
99  const SourceInterface& source,
100  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
101  std::vector<ModelFitting::PointModel>& /*point_models*/,
104  std::shared_ptr<CoordinateSystem> reference_coordinates,
105  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
106 
107  auto pixel_x = createDependentParameter(
108  [reference_coordinates, coordinates, offset](double x, double y) {
109  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
110  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
111 
112 
113  auto pixel_y = createDependentParameter(
114  [reference_coordinates, coordinates, offset](double x, double y) {
115  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
116  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
117 
118  auto n = std::make_shared<ManualParameter>(1); // Sersic index for exponential
119  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
120 
121  auto i0 = createDependentParameter(
122  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.35513 * radius * radius * aspect); },
123  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
124  manager.getParameter(source, m_aspect_ratio));
125 
126  auto k = createDependentParameter(
127  [](double eff_radius) { return 1.678 / eff_radius; },
128  manager.getParameter(source, m_effective_radius));
129 
131  sersic_component.emplace_back(new SersicModelComponent(make_unique<OldSharp>(), i0, n, k));
132 
133  auto& boundaries = source.getProperty<PixelBoundaries>();
134  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
135 
136  auto minus_angle = createDependentParameter(
137  [](double angle) { return -angle; },
138  manager.getParameter(source, m_angle));
139 
140  extended_models.emplace_back(
141  std::move(sersic_component), x_scale, manager.getParameter(source, m_aspect_ratio), minus_angle,
142  size, size, pixel_x, pixel_y, jacobian);
143 }
144 
146  const SourceInterface& source,
147  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
148  std::vector<ModelFitting::PointModel>& /*point_models*/,
151  std::shared_ptr<CoordinateSystem> reference_coordinates,
152  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
153 
154  auto pixel_x = createDependentParameter(
155  [reference_coordinates, coordinates, offset](double x, double y) {
156  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
157  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
158 
159 
160  auto pixel_y = createDependentParameter(
161  [reference_coordinates, coordinates, offset](double x, double y) {
162  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
163  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
164 
165  auto n = std::make_shared<ManualParameter>(4); // Sersic index for Devaucouleurs
166  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
167 
168  auto i0 = createDependentParameter(
169  [](double flux, double radius, double aspect) { return flux / (2 * M_PI * 0.001684925 * radius * radius * aspect); },
170  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
171  manager.getParameter(source, m_aspect_ratio));
172 
173  auto k = createDependentParameter(
174  [](double eff_radius) { return 7.669 / pow(eff_radius, .25); },
175  manager.getParameter(source, m_effective_radius));
176 
178  sersic_component.emplace_back(new SersicModelComponent(make_unique<OldSharp>(), i0, n, k));
179 
180  auto& boundaries = source.getProperty<PixelBoundaries>();
181  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
182 
183  auto minus_angle = createDependentParameter(
184  [](double angle) { return -angle; },
185  manager.getParameter(source, m_angle));
186 
187  extended_models.emplace_back(
188  std::move(sersic_component), x_scale, manager.getParameter(source, m_aspect_ratio), minus_angle,
189  size, size, pixel_x, pixel_y, jacobian);
190 }
191 
192 static double computeBn(double n) {
193  // Using approximation from MacArthur, L.A., Courteau, S., & Holtzman, J.A. 2003, ApJ, 582, 689
194  return 2 * n - 1.0 / 3.0 + 4 / (405 * n)
195  + 46 / (25515 * n * n) + 131 / (1148175 * n * n * n) - 2194697 / (30690717750 * n * n * n * n);
196 }
197 
199  const SourceInterface& source,
200  std::vector<ModelFitting::ConstantModel>& /* constant_models */,
201  std::vector<ModelFitting::PointModel>& /*point_models*/,
204  std::shared_ptr<CoordinateSystem> reference_coordinates,
205  std::shared_ptr<CoordinateSystem> coordinates, PixelCoordinate offset) const {
206 
207  auto pixel_x = createDependentParameter(
208  [reference_coordinates, coordinates, offset](double x, double y) {
209  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_x - offset.m_x + 0.5;
210  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
211 
212 
213  auto pixel_y = createDependentParameter(
214  [reference_coordinates, coordinates, offset](double x, double y) {
215  return coordinates->worldToImage(reference_coordinates->imageToWorld(ImageCoordinate(x-1, y-1))).m_y - offset.m_y + 0.5;
216  }, manager.getParameter(source, m_x), manager.getParameter(source, m_y));
217 
218  auto x_scale = std::make_shared<ManualParameter>(1); // we don't scale the x coordinate
219 
220  auto i0 = createDependentParameter(
221  [](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); },
222  manager.getParameter(source, m_flux), manager.getParameter(source, m_effective_radius),
223  manager.getParameter(source, m_aspect_ratio), manager.getParameter(source, m_sersic_index));
224 
225  auto k = createDependentParameter(
226  [](double eff_radius, double n) { return computeBn(n) / pow(eff_radius, 1.0 / n); },
227  manager.getParameter(source, m_effective_radius), manager.getParameter(source, m_sersic_index));
228 
230  sersic_component.emplace_back(new SersicModelComponent(make_unique<OldSharp>(), i0,
231  manager.getParameter(source, m_sersic_index), k));
232 
233  auto& boundaries = source.getProperty<PixelBoundaries>();
234  int size = std::max(MODEL_MIN_SIZE, MODEL_SIZE_FACTOR * std::max(boundaries.getWidth(), boundaries.getHeight()));
235 
236  extended_models.emplace_back(
237  std::move(sersic_component), x_scale, manager.getParameter(source, m_aspect_ratio),
238  manager.getParameter(source, m_angle), size, size, pixel_x, pixel_y, jacobian);
239 }
240 
242  const SourceInterface& source,
244  std::vector<ModelFitting::PointModel>& /* point_models */,
245  std::vector<ModelFitting::TransformedModel>& /* extended_models */,
247  std::shared_ptr<CoordinateSystem> /* reference_coordinates */,
248  std::shared_ptr<CoordinateSystem> /* coordinates */, PixelCoordinate /* offset */) const {
249  constant_models.emplace_back(manager.getParameter(source, m_value));
250 }
251 
252 }
253 
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< ModelFitting::TransformedModel > &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
static double computeBn(double n)
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
static const double MODEL_SIZE_FACTOR
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive...
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
CircularlySymmetricModelComponent< SersicProfile > SersicModelComponent
std::shared_ptr< DependentParameter< Parameters...> > createDependentParameter(typename DependentParameter< Parameters...>::ValueCalculator value_calculator, Parameters...parameters)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T tgamma(T...args)
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< ModelFitting::TransformedModel > &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< ModelFitting::TransformedModel > &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
T max(T...args)
A pixel coordinate made of two integers m_x and m_y.
T move(T...args)
static const double MODEL_MIN_SIZE
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
T pow(T...args)
The SourceInterface is an abstract &quot;source&quot; that has properties attached to it.
virtual void addForSource(FlexibleModelFittingParameterManager &manager, const SourceInterface &source, std::vector< ModelFitting::ConstantModel > &constant_models, std::vector< ModelFitting::PointModel > &point_models, std::vector< ModelFitting::TransformedModel > &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< ModelFitting::TransformedModel > &extended_models, std::tuple< double, double, double, double > jacobian, std::shared_ptr< CoordinateSystem > reference_coordinates, std::shared_ptr< CoordinateSystem > coordinates, PixelCoordinate offset) const
T emplace_back(T...args)