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