SourceXtractorPlusPlus  0.10
Please provide a description of the project.
MoffatModelFittingTask.cpp
Go to the documentation of this file.
1 
17 /*
18  * MoffatModelTask.cpp
19  *
20  * Created on: May 2, 2017
21  * Author: mschefer
22  */
23 
24 #include <iostream>
25 #include <tuple>
26 #include <vector>
27 #include <valarray>
28 #include <boost/any.hpp>
29 #include <mutex>
30 
36 
39 
42 
55 
59 
61 
62 
65 
67 
75 
77 
79 
81 
82 
83 namespace SourceXtractor {
84 
85 using namespace ModelFitting;
87 
88 namespace {
89 
90 struct SourceModel {
91  double m_size;
94 
95  double exp_i0_guess;
98 
99  SourceModel(double size, double x_guess, double y_guess, double pos_range,
100  double exp_flux_guess, double exp_radius_guess, double exp_aspect_guess, double exp_rot_guess) :
101 
102  m_size(size),
103  dx(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
104  dy(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
105 
106  x(createDependentParameter([x_guess](double dx) { return dx + x_guess + 0.5; }, dx)),
107  y(createDependentParameter([y_guess](double dy) { return dy + y_guess + 0.5; }, dy)),
108 
109  // FIXME
110  exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
111  moffat_i0(std::make_shared<EngineParameter>(exp_i0_guess, make_unique<ExpSigmoidConverter>(exp_i0_guess * .00001, exp_i0_guess * 1000))),
112 
113  moffat_index(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.5, 8))),
114  minkowski_exponent(std::make_shared<EngineParameter>(2, make_unique<ExpSigmoidConverter>(0.5, 10))),
115  flat_top_offset(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.000001, 10))),
116 
117  moffat_x_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
118  moffat_y_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
119  moffat_rotation(std::make_shared<EngineParameter>(-exp_rot_guess, make_unique<SigmoidConverter>(-2*M_PI, 2*M_PI)))
120  {
121  }
122 
123  void registerParameters(EngineParameterManager& manager) {
124  manager.registerParameter(dx);
125  manager.registerParameter(dy);
126 
127  manager.registerParameter(moffat_i0);
131 
135  }
136 
137  void createModels(std::vector<std::shared_ptr<ModelFitting::ExtendedModel<ImageInterfaceTypePtr>>>& extended_models, std::vector<PointModel>& /*point_models*/) {
138  // Moffat model
139  {
141  auto moff = make_unique<FlattenedMoffatComponent>(moffat_i0, moffat_index, minkowski_exponent, flat_top_offset);
142  component_list.clear();
143  component_list.emplace_back(std::move(moff));
146  }
147  }
148 };
149 
150 }
151 
152 
154  auto& source_stamp = source.getProperty<DetectionFrameSourceStamp>().getStamp();
155  auto& variance_stamp = source.getProperty<DetectionFrameSourceStamp>().getVarianceStamp();
156  auto& thresholded_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdedStamp();
157  auto& threshold_map_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdMapStamp();
158  PixelCoordinate stamp_top_left = source.getProperty<DetectionFrameSourceStamp>().getTopLeft();
159 
160  // Computes the minimum flux that a detection should have (min. detection threshold for every pixel)
161  // This will be used instead of lower or negative fluxes that can happen for various reasons
162  double min_flux = 0.;
163  auto& pixel_coordinates = source.getProperty<PixelCoordinateList>().getCoordinateList();
164  for (auto pixel : pixel_coordinates) {
165  pixel -= stamp_top_left;
166 
167  min_flux += threshold_map_stamp.getValue(pixel);
168  }
169 
170  double pixel_scale = 1;
171 
172  EngineParameterManager manager {};
173  std::vector<ConstantModel> constant_models;
175  std::vector<PointModel> point_models;
176 
177  auto& pixel_centroid = source.getProperty<PixelCentroid>();
178  auto& shape_parameters = source.getProperty<ShapeParameters>();
179  auto iso_flux = source.getProperty<IsophotalFlux>().getFlux();
180 
181  auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
182 
183  double radius_guess = shape_parameters.getEllipseA() / 2.0;
184 
185  double guess_x = pixel_centroid.getCentroidX() - stamp_top_left.m_x;
186  double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.m_y;
187 
188  double exp_flux_guess = std::max<double>(iso_flux, min_flux);
189 
190  double exp_reff_guess = radius_guess;
191  double exp_aspect_guess = std::max<double>(shape_parameters.getEllipseB() / shape_parameters.getEllipseA(), 0.01);
192  double exp_rot_guess = shape_parameters.getEllipseTheta();
193 
194  auto source_model = make_unique<SourceModel>(size, guess_x, guess_y, radius_guess * 2,
195  exp_flux_guess, exp_reff_guess, exp_aspect_guess, exp_rot_guess);
196 
197  source_model->createModels(extended_models, point_models);
198  source_model->registerParameters(manager);
199 
200  // Full frame model with all sources
202  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model {
203  pixel_scale,
204  (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(),
205  std::move(constant_models),
206  std::move(point_models),
207  std::move(extended_models)
208  };
209 
210  auto weight = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
211  std::fill(weight->getData().begin(), weight->getData().end(), 1);
212 
213  for (int y=0; y < source_stamp.getHeight(); y++) {
214  for (int x=0; x < source_stamp.getWidth(); x++) {
215  weight->at(x, y) = (thresholded_stamp.getValue(x, y) >= 0) ? 0 : 1;
216  }
217  }
218 
219  for (auto pixel : pixel_coordinates) {
220  pixel -= stamp_top_left;
221  weight->at(pixel.m_x, pixel.m_y) = 1;
222  }
223 
224  auto frame = source.getProperty<DetectionFrame>().getFrame();
225  SeFloat gain = frame->getGain();
226  SeFloat saturation = frame->getSaturation();
227 
228  for (int y=0; y < source_stamp.getHeight(); y++) {
229  for (int x=0; x < source_stamp.getWidth(); x++) {
230  auto back_var = variance_stamp.getValue(x, y);
231  if (saturation > 0 && source_stamp.getValue(x, y) >= saturation) {
232  weight->at(x, y) = 0;
233  } else if (weight->at(x, y)>0) {
234  if (gain > 0.0) {
235  weight->at(x, y) = sqrt(1.0 / (back_var + source_stamp.getValue(x, y) / gain));
236  } else {
237  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
238  }
239  }
240  }
241  }
242 
243  // FIXME we should be able to use the source_stamp Image interface directly
244  auto image = VectorImage<SeFloat>::create(source_stamp);
245 
246  auto data_vs_model =
247  createDataVsModelResiduals(image, std::move(frame_model), weight, AsinhChiSquareComparator{});
248 
249  ResidualEstimator res_estimator {};
250  res_estimator.registerBlockProvider(move(data_vs_model));
251 
252  // Perform the minimization
253  auto engine = LeastSquareEngineManager::create(m_least_squares_engine, m_max_iterations);
254  auto solution = engine->solveProblem(manager, res_estimator);
255  size_t iterations = (size_t) boost::any_cast<std::array<double,10>>(solution.underlying_framework_info)[5];
256 
257  auto final_stamp = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
258 
259  {
260 
261  // renders an image of the model for a single source with the final parameters
263  std::vector<PointModel> point_models {};
264  source_model->createModels(extended_models, point_models);
265  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model_after {
266  1, (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(), std::move(constant_models), std::move(point_models),
267  std::move(extended_models)
268  };
269  auto final_image = frame_model_after.getImage();
270 
271  // integrates the flux for that source
272  double total_flux = 0;
273  for (int y=0; y < source_stamp.getHeight(); y++) {
274  for (int x=0; x < source_stamp.getWidth(); x++) {
275  PixelCoordinate pixel(x, y);
276  pixel += stamp_top_left;
277 
278  // build final stamp
279  final_stamp->setValue(x, y, final_stamp->getValue(x, y) + final_image->getValue(x, y));
280 
281  total_flux += final_image->getValue(x, y);
282  }
283  }
284 
285  auto coordinate_system = source.getProperty<DetectionFrame>().getFrame()->getCoordinateSystem();
286 
287  SeFloat x = stamp_top_left.m_x + source_model->x->getValue() - 0.5f;
288  SeFloat y = stamp_top_left.m_y + source_model->y->getValue() - 0.5f;
289 
291  x, y,
292 
293  source_model->moffat_i0->getValue(),
294  source_model->moffat_index->getValue(),
295  source_model->minkowski_exponent->getValue(),
296  source_model->flat_top_offset->getValue(),
297  source_model->m_size,
298  source_model->moffat_x_scale->getValue(),
299  source_model->moffat_y_scale->getValue(),
300  source_model->moffat_rotation->getValue(),
301 
302  iterations
303  );
304 
305  }
306 }
307 
308 }
309 
virtual void computeProperties(SourceInterface &source) const override
Computes one or more properties for the Source.
STL class.
EngineParameter are those derived from the minimization process.
CoordinateConverter implementation using the sigmoid function.
Computes the isophotal flux and magnitude.
Definition: IsophotalFlux.h:36
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
STL namespace.
SeFloat32 SeFloat
Definition: Types.h:32
double exp_i0_guess
The centroid of all the pixels in the source, weighted by their DetectionImage pixel values.
Definition: PixelCentroid.h:37
std::shared_ptr< EngineParameter > moffat_y_scale
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
std::shared_ptr< EngineParameter > minkowski_exponent
SeFloat getCentroidX() const
X coordinate of centroid.
Definition: PixelCentroid.h:48
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > moffat_index
A copy of the rectangular region of the detection image just large enough to include the whole Source...
void registerBlockProvider(std::unique_ptr< ResidualBlockProvider > provider)
Registers a ResidualBlockProvider to the ResidualEstimator.
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:89
std::shared_ptr< EngineParameter > moffat_i0
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
void registerParameter(std::shared_ptr< EngineParameter > parameter)
Registers an EngineParameter to the EngineParameterManager.
A pixel coordinate made of two integers m_x and m_y.
T move(T... args)
Data vs model comparator which computes a modified residual, using asinh.
STL class.
T make_shared(T... args)
std::unique_ptr< DataVsModelResiduals< typename std::remove_reference< DataType >::type, typename std::remove_reference< ModelType >::type, typename std::remove_reference< WeightType >::type, typename std::remove_reference< Comparator >::type > > createDataVsModelResiduals(DataType &&data, ModelType &&model, WeightType &&weight, Comparator &&comparator)
STL class.
Class responsible for managing the parameters the least square engine minimizes.
static std::shared_ptr< LeastSquareEngine > create(const std::string &name, unsigned max_iterations=1000)
T sqrt(T... args)
std::shared_ptr< EngineParameter > moffat_x_scale
T fill(T... args)
Provides to the LeastSquareEngine the residual values.
std::shared_ptr< DependentParameter< Parameters... > > createDependentParameter(typename DependentParameter< Parameters... >::ValueCalculator value_calculator, Parameters... parameters)
double m_size
The SourceInterface is an abstract "source" that has properties attached to it.
std::shared_ptr< EngineParameter > flat_top_offset
const double pixel_scale
Definition: TestImage.cpp:75
std::shared_ptr< EngineParameter > moffat_rotation
std::shared_ptr< EngineParameter > dy
std::unique_ptr< T > make_unique(Args &&... args)