SourceXtractorPlusPlus  0.10
Please provide a description of the project.
FlexibleModelFittingTask.cpp
Go to the documentation of this file.
1 
17 /*
18  * FlexibleModelFittingTask.cpp
19  *
20  * Created on: Sep 17, 2018
21  * Author: mschefer
22  */
23 
24 #include <mutex>
25 
35 
37 
39 
42 
49 
53 
55 
56 namespace SourceXtractor {
57 
58 using namespace ModelFitting;
59 
60 
61 namespace {
62 
63 __attribute__((unused))
65  std::cerr << "\nMinimization info:\n";
66  std::cerr << " ||e||_2 at initial p: " << info[0] << '\n';
67  std::cerr << " ||e||_2: " << info[1] << '\n';
68  std::cerr << " ||J^T e||_inf: " << info[2] << '\n';
69  std::cerr << " ||Dp||_2: " << info[3] << '\n';
70  std::cerr << " mu/max[J^T J]_ii: " << info[4] << '\n';
71  std::cerr << " # iterations: " << info[5] << '\n';
72  switch ((int) info[6]) {
73  case 1:
74  std::cerr << " stopped by small gradient J^T e\n";
75  break;
76  case 2:
77  std::cerr << " stopped by small Dp\n";
78  break;
79  case 3:
80  std::cerr << " stopped by itmax\n";
81  break;
82  case 4:
83  std::cerr << " singular matrix. Restart from current p with increased mu\n";
84  break;
85  case 5:
86  std::cerr << " no further error reduction is possible. Restart with increased mu\n";
87  break;
88  case 6:
89  std::cerr << " stopped by small ||e||_2\n";
90  break;
91  case 7:
92  std::cerr << " stopped by invalid (i.e. NaN or Inf) func values; a user error\n";
93  break;
94  }
95  std::cerr << " # function evaluations: " << info[7] << '\n';
96  std::cerr << " # Jacobian evaluations: " << info[8] << '\n';
97  std::cerr << " # linear systems solved: " << info[9] << "\n\n";
98 }
99 
100 }
101 
103  unsigned int max_iterations, double modified_chi_squared_scale,
107  : m_least_squares_engine(least_squares_engine),
108  m_max_iterations(max_iterations), m_modified_chi_squared_scale(modified_chi_squared_scale),
109  m_parameters(parameters), m_frames(frames), m_priors(priors) {}
110 
112  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
113  return stamp_rect.getWidth() > 0 && stamp_rect.getHeight() > 0;
114 }
115 
117  SourceGroupInterface& group, int frame_index) const {
119 
120  auto frame = group.begin()->getProperty<MeasurementFrame>(frame_index).getFrame();
121  auto frame_image = frame->getSubtractedImage();
122 
123  auto rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
124  auto image = VectorImage<SeFloat>::create(rect.getWidth(), rect.getHeight());
125  for (int y = 0; y < rect.getHeight(); y++) {
126  for (int x = 0; x < rect.getWidth(); x++) {
127  image->at(x, y) = frame_image->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
128  }
129  }
130 
131  return image;
132 }
133 
135  SourceGroupInterface& group, int frame_index) const {
137 
138  auto frame = group.begin()->getProperty<MeasurementFrame>(frame_index).getFrame();
139  auto frame_image = frame->getSubtractedImage();
140  auto frame_image_thresholded = frame->getThresholdedImage();
141  auto variance_map = frame->getVarianceMap();
142 
143  auto rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
144  auto weight = VectorImage<SeFloat>::create(rect.getWidth(), rect.getHeight());
145  std::fill(weight->getData().begin(), weight->getData().end(), 1);
146 
147  auto measurement_frame = group.begin()->getProperty<MeasurementFrame>(frame_index).getFrame();
148  SeFloat gain = measurement_frame->getGain();
149  SeFloat saturation = measurement_frame->getSaturation();
150 
151  for (int y = 0; y < rect.getHeight(); y++) {
152  for (int x = 0; x < rect.getWidth(); x++) {
153  auto back_var = variance_map->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
154  if (saturation > 0 && frame_image->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y) > saturation) {
155  weight->at(x, y) = 0;
156  } else if (weight->at(x, y) > 0) {
157  if (gain > 0.0) {
158  weight->at(x, y) = sqrt(
159  1.0 / (back_var + frame_image->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y) / gain));
160  } else {
161  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
162  }
163  }
164  }
165  }
166 
167  return weight;
168 }
169 
171  SourceGroupInterface& group,
174 
175  int frame_index = frame->getFrameNb();
176 
177  auto frame_coordinates =
178  group.begin()->getProperty<MeasurementFrame>(frame_index).getFrame()->getCoordinateSystem();
179  auto ref_coordinates =
180  group.begin()->getProperty<DetectionFrame>().getFrame()->getCoordinateSystem();
181 
182  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
183  auto psf_property = group.getProperty<PsfProperty>(frame_index);
184  auto jacobian = group.getProperty<JacobianGroup>(frame_index).asTuple();
185 
186  // The model fitting module expects to get a PSF with a pixel scale, but we have the pixel sampling step size
187  // It will be used to compute the rastering grid size, and after convolving with the PSF the result will be
188  // downscaled before copied into the frame image.
189  // We can multiply here then, as the unit is pixel/pixel, rather than "/pixel or similar
190  auto group_psf = ImagePsf(pixel_scale * psf_property.getPixelSampling(), psf_property.getPsf());
191 
192  std::vector<ConstantModel> constant_models;
193  std::vector<PointModel> point_models;
195 
196  for (auto& source : group) {
197  for (auto model : frame->getModels()) {
198  model->addForSource(manager, source, constant_models, point_models, extended_models, jacobian, ref_coordinates, frame_coordinates,
199  stamp_rect.getTopLeft());
200  }
201  }
202 
203  // Full frame model with all sources
205  pixel_scale, (size_t) stamp_rect.getWidth(), (size_t) stamp_rect.getHeight(),
206  std::move(constant_models), std::move(point_models), std::move(extended_models), group_psf);
207 
208  return frame_model;
209 }
210 
211 
213  double pixel_scale = 1;
214  FlexibleModelFittingParameterManager parameter_manager;
215  ModelFitting::EngineParameterManager engine_parameter_manager{};
216  int n_free_parameters = 0;
217 
218  {
220 
221  // Prepare parameters
222  for (auto& source : group) {
223  for (auto parameter : m_parameters) {
224  if (std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter)) {
225  ++n_free_parameters;
226  }
227  parameter_manager.addParameter(source, parameter,
228  parameter->create(parameter_manager, engine_parameter_manager, source));
229  }
230  }
231  }
232 
233  // Reset access checks, as a dependent parameter could have triggered it
234  parameter_manager.clearAccessCheck();
235 
236  // Add models for all frames
237  ResidualEstimator res_estimator{};
238 
239  int valid_frames = 0;
240  int n_good_pixels = 0;
241  for (auto frame : m_frames) {
242  int frame_index = frame->getFrameNb();
243  // Validate that each frame covers the model fitting region
244  if (isFrameValid(group, frame_index)) {
245  valid_frames++;
246 
247  auto frame_model = createFrameModel(group, pixel_scale, parameter_manager, frame);
248 
249  auto image = createImageCopy(group, frame_index);
250  auto weight = createWeightImage(group, frame_index);
251 
252  for (int y = 0; y < weight->getHeight(); ++y) {
253  for (int x = 0; x < weight->getWidth(); ++x) {
254  n_good_pixels += (weight->at(x, y) != 0.);
255  }
256  }
257 
258  // Setup residuals
259  auto data_vs_model =
260  createDataVsModelResiduals(image, std::move(frame_model), weight,
261  //LogChiSquareComparator(m_modified_chi_squared_scale));
263  res_estimator.registerBlockProvider(std::move(data_vs_model));
264  }
265  }
266 
267  // Check that we had enough data for the fit
268  Flags group_flags = Flags::NONE;
269  if (valid_frames == 0) {
270  group_flags = Flags::OUTSIDE;
271  }
272  else if (n_good_pixels < n_free_parameters) {
273  group_flags = Flags::INSUFFICIENT_DATA;
274  }
275 
276  if (group_flags != Flags::NONE) {
277  // Can't do model fitting as no measurement frame overlaps the detected source, or there are not enough pixels
278  // We still need to provide a property
279  for (auto& source : group) {
280  std::unordered_map<int, double> dummy_values;
281  for (auto parameter : m_parameters) {
282  auto modelfitting_parameter = parameter_manager.getParameter(source, parameter);
283  auto manual_parameter = std::dynamic_pointer_cast<ManualParameter>(modelfitting_parameter);
284  if (manual_parameter) {
285  manual_parameter->setValue(std::numeric_limits<double>::quiet_NaN());
286  }
287  dummy_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
288  }
289  source.setProperty<FlexibleModelFitting>(0, std::numeric_limits<double>::quiet_NaN(), group_flags,
290  dummy_values, dummy_values);
291  }
292  return;
293  }
294 
295  // Add priors
296  for (auto& source : group) {
297  for (auto prior : m_priors) {
298  prior->setupPrior(parameter_manager, source, res_estimator);
299  }
300  }
301 
302  // Model fitting
303 
304  // FIXME we can no longer specify different settings with LeastSquareEngineManager!!
305  // LevmarEngine engine{m_max_iterations, 1E-3, 1E-6, 1E-6, 1E-6, 1E-4};
307  auto solution = engine->solveProblem(engine_parameter_manager, res_estimator);
308 
309  size_t iterations = (size_t) boost::any_cast<std::array<double, 10>>(solution.underlying_framework_info)[5];
310 
311  int total_data_points = 0;
312  SeFloat avg_reduced_chi_squared = computeChiSquared(group, pixel_scale, parameter_manager, total_data_points);
313 
314  int nb_of_free_parameters = 0;
315  for (auto& source : group) {
316  for (auto parameter : m_parameters) {
317  bool is_free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter).get();
318  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
319  if (is_free_parameter && accessed_by_modelfitting) {
320  nb_of_free_parameters++;
321  }
322  }
323  }
324  avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
325 
326  // Collect parameters for output
327  for (auto& source : group) {
328  std::unordered_map<int, double> parameter_values, parameter_sigmas;
329  auto source_flags = Flags::NONE;
330 
331  for (auto parameter : m_parameters) {
332  bool is_dependent_parameter = std::dynamic_pointer_cast<FlexibleModelFittingDependentParameter>(parameter).get();
333  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
334  auto modelfitting_parameter = parameter_manager.getParameter(source, parameter);
335 
336  if (is_dependent_parameter || accessed_by_modelfitting) {
337  parameter_values[parameter->getId()] = modelfitting_parameter->getValue();
338  parameter_sigmas[parameter->getId()] = parameter->getSigma(parameter_manager, source, solution.parameter_sigmas);
339  }
340  else {
341  // Need to cascade the NaN to any potential dependent parameter
342  auto engine_parameter = std::dynamic_pointer_cast<EngineParameter>(modelfitting_parameter);
343  if (engine_parameter) {
344  engine_parameter->setEngineValue(std::numeric_limits<double>::quiet_NaN());
345  }
346  parameter_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
347  parameter_sigmas[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
349  }
350  }
351 
352  source.setProperty<FlexibleModelFitting>(iterations, avg_reduced_chi_squared, source_flags, parameter_values,
353  parameter_sigmas);
354  }
355 
356  updateCheckImages(group, pixel_scale, parameter_manager);
357 }
358 
360  double pixel_scale, FlexibleModelFittingParameterManager& manager) const {
361 
362  int frame_id = 0;
363  for (auto frame : m_frames) {
364  int frame_index = frame->getFrameNb();
365  // Validate that each frame covers the model fitting region
366  if (isFrameValid(group, frame_index)) {
367  auto frame_model = createFrameModel(group, pixel_scale, manager, frame);
368  auto final_stamp = frame_model.getImage();
369 
370  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
371  auto frame = group.begin()->getProperty<MeasurementFrame>(frame_index).getFrame();
372 
373  auto debug_image = CheckImages::getInstance().getModelFittingImage(frame);
374  if (debug_image) {
376  for (int x = 0; x < final_stamp->getWidth(); x++) {
377  for (int y = 0; y < final_stamp->getHeight(); y++) {
378  auto x_coord = stamp_rect.getTopLeft().m_x + x;
379  auto y_coord = stamp_rect.getTopLeft().m_y + y;
380  debug_image->setValue(x_coord, y_coord,
381  debug_image->getValue(x_coord, y_coord) + final_stamp->getValue(x, y));
382  }
383  }
384  }
385  }
386  frame_id++;
387  }
388 }
389 
391  std::shared_ptr<const Image<SeFloat>> model, std::shared_ptr<const Image<SeFloat>> weights, int& data_points) const {
392  double reduced_chi_squared = 0.0;
393  data_points = 0;
394  for (int y=0; y < image->getHeight(); y++) {
395  for (int x=0; x < image->getWidth(); x++) {
396  double tmp = image->getValue(x, y) - model->getValue(x, y);
397  reduced_chi_squared += tmp * tmp * weights->getValue(x, y) * weights->getValue(x, y);
398  if (weights->getValue(x, y) > 0) {
399  data_points++;
400  }
401  }
402  }
403  return reduced_chi_squared;
404 }
405 
407  double pixel_scale, FlexibleModelFittingParameterManager& manager, int& total_data_points) const {
408 
409  SeFloat total_chi_squared = 0;
410  total_data_points = 0;
411  int valid_frames = 0;
412  for (auto frame : m_frames) {
413  int frame_index = frame->getFrameNb();
414  // Validate that each frame covers the model fitting region
415  if (isFrameValid(group, frame_index)) {
416  valid_frames++;
417  auto frame_model = createFrameModel(group, pixel_scale, manager, frame);
418  auto final_stamp = frame_model.getImage();
419 
420  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
421  auto image = createImageCopy(group, frame_index);
422  auto weight = createWeightImage(group, frame_index);
423 
424  int data_points = 0;
425  SeFloat chi_squared = computeChiSquaredForFrame(
426  image, final_stamp, weight, data_points);
427 
428  total_data_points += data_points;
429  total_chi_squared += chi_squared;
430  }
431  }
432 
433  return total_chi_squared;
434 }
435 
437 }
438 
439 }
static StaticPlugin< SourceFlagsPlugin > source_flags
bool isParamAccessed(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
STL class.
void printLevmarInfo(std::array< double, 10 > info)
Definition: utils.h:118
void addParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter, std::shared_ptr< ModelFitting::BasicParameter > engine_parameter)
SeFloat computeChiSquared(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager, int &total_data_points) const
The object is completely outside of the measurement frame.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
ModelFitting::FrameModel< ImagePsf, std::shared_ptr< VectorImage< SourceXtractor::SeFloat > > > createFrameModel(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager, std::shared_ptr< FlexibleModelFittingFrame > frame) const
virtual void computeProperties(SourceGroupInterface &group) const override
Computes one or more properties for the SourceGroup and/or the Sources it contains.
SeFloat32 SeFloat
Definition: Types.h:32
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
STL class.
FlexibleModelFittingTask(const std::string &least_squares_engine, unsigned int max_iterations, double modified_chi_squared_scale, std::vector< std::shared_ptr< FlexibleModelFittingParameter >> parameters, std::vector< std::shared_ptr< FlexibleModelFittingFrame >> frames, std::vector< std::shared_ptr< FlexibleModelFittingPrior >> priors)
std::shared_ptr< WriteableImage< MeasurementImage::PixelType > > getModelFittingImage(std::shared_ptr< const MeasurementImageFrame > frame)
SeFloat computeChiSquaredForFrame(std::shared_ptr< const Image< SeFloat >> image, std::shared_ptr< const Image< SeFloat >> model, std::shared_ptr< const Image< SeFloat >> weights, int &data_points) const
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:89
static CheckImages & getInstance()
Definition: CheckImages.h:114
std::shared_ptr< VectorImage< SeFloat > > createImageCopy(SourceGroupInterface &group, int frame_index) const
T lock(T... args)
#define __attribute__(xyz)
T move(T... args)
Data vs model comparator which computes a modified residual, using asinh.
bool isFrameValid(SourceGroupInterface &group, int frame_index) const
Defines the interface used to group sources.
std::vector< std::shared_ptr< FlexibleModelFittingPrior > > m_priors
std::vector< std::shared_ptr< FlexibleModelFittingFrame > > m_frames
STL class.
There are not enough good pixels to fit the parameters.
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
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)
void updateCheckImages(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager) const
Some/all of the model parameters could not be fitted.
Flags
Flagging of bad sources.
Definition: SourceFlags.h:34
STL class.
Class responsible for managing the parameters the least square engine minimizes.
std::shared_ptr< VectorImage< SeFloat > > createWeightImage(SourceGroupInterface &group, int frame_index) const
static std::shared_ptr< LeastSquareEngine > create(const std::string &name, unsigned max_iterations=1000)
std::vector< std::shared_ptr< FlexibleModelFittingParameter > > m_parameters
T quiet_NaN(T... args)
T sqrt(T... args)
T fill(T... args)
Provides to the LeastSquareEngine the residual values.
T getValue(int x, int y) const override
Returns the value of the pixel with the coordinates (x,y)
Definition: VectorImage.h:106
const double pixel_scale
Definition: TestImage.cpp:75
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const