SourceXtractorPlusPlus  0.13
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 
SourceXtractor::PixelBoundaries
The bounding box of all the pixels in the source. Both min and max coordinate are inclusive.
Definition: PixelBoundaries.h:37
AutoSharp.h
SourceXtractor::PixelCoordinate
A pixel coordinate made of two integers m_x and m_y.
Definition: PixelCoordinate.h:37
PixelBoundaries.h
std::shared_ptr
STL class.
SourceXtractor::FlexibleModelFittingExponentialModel::addForSource
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
Definition: FlexibleModelFittingModel.cpp:109
CompactExponentialModel.h
EngineParameter.h
SourceXtractor::FlexibleModelFittingConstantModel::addForSource
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
Definition: FlexibleModelFittingModel.cpp:234
SourceXtractor::FlexibleModelFittingSersicModel::addForSource
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
Definition: FlexibleModelFittingModel.cpp:196
std::make_shared
T make_shared(T... args)
std::vector< ModelFitting::ConstantModel >
CompactSersicModel.h
SourceXtractor::FlexibleModelFittingPointModel::addForSource
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
Definition: FlexibleModelFittingModel.cpp:84
FlexibleModelFittingModel.h
OnlySmooth.h
std::tuple< double, double, double, double >
Euclid::Configuration
ModelFitting::CompactExponentialModel
Definition: CompactExponentialModel.h:16
SourceXtractor::FlexibleModelFittingParameterManager
Definition: FlexibleModelFittingParameterManager.h:43
SourceXtractor
Definition: Aperture.h:30
CircularlySymmetricModelComponent.h
SourceXtractor::PixelCoordinate::m_x
int m_x
Definition: PixelCoordinate.h:38
SourceXtractor::FlexibleModelFittingParameterManager::getParameter
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
Definition: FlexibleModelFittingParameterManager.h:51
SourceXtractor::PixelCoordinate::m_y
int m_y
Definition: PixelCoordinate.h:38
ManualParameter.h
ImageInterfaceTraits.h
SourceXtractor::FlexibleModelFittingDevaucouleursModel::addForSource
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
Definition: FlexibleModelFittingModel.cpp:150
SourceXtractor::ImageCoordinate
Definition: CoordinateSystem.h:42
Euclid::make_unique
std::unique_ptr< T > make_unique(Args &&... args)
std::tgamma
T tgamma(T... args)
SigmoidConverter.h
FlexibleModelFittingParameter.h
OldSharp.h
std::vector::emplace_back
T emplace_back(T... args)
DependentParameter.h
ModelFitting::createDependentParameter
std::shared_ptr< DependentParameter< Parameters... > > createDependentParameter(typename DependentParameter< Parameters... >::ValueCalculator value_calculator, Parameters... parameters)
Definition: DependentParameter.h:131
ExpSigmoidConverter.h
x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
Definition: MoffatModelFittingTask.cpp:94
ModelFitting::ExtendedModel
Definition: ExtendedModel.h:39
SourceXtractor::computeBn
static double computeBn(double n)
Definition: FlexibleModelFittingModel.cpp:190
SourceXtractor::SourceInterface::getProperty
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
Definition: SourceInterface.h:57
SourceXtractor::MODEL_SIZE_FACTOR
static const double MODEL_SIZE_FACTOR
Definition: FlexibleModelFittingModel.cpp:64
TransformedModel.h
SourceXtractor::MODEL_MIN_SIZE
static const double MODEL_MIN_SIZE
Definition: FlexibleModelFittingModel.cpp:63
ModelFitting::CompactSersicModel
Definition: CompactSersicModel.h:16
SourceXtractor::SourceInterface
The SourceInterface is an abstract "source" that has properties attached to it.
Definition: SourceInterface.h:46
std::max
T max(T... args)
y
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
Definition: MoffatModelFittingTask.cpp:94
ConfigManager.h
ModelFitting
Definition: AsinhChiSquareComparator.h:30
memory_tools.h
SamplingConfig.h
EngineParameterManager.h
FlexibleModelFittingParameterManager.h
std::pow
T pow(T... args)