SourceXtractorPlusPlus  0.16
Please provide a description of the project.
FlexibleModelFittingIterativeTask.cpp
Go to the documentation of this file.
1 
22 
26 
28 
36 
40 
41 namespace SourceXtractor {
42 
43 using namespace ModelFitting;
44 
45 static auto logger = Elements::Logging::getLogger("FlexibleModelFitting");
46 
48  unsigned int max_iterations, double modified_chi_squared_scale,
52  double scale_factor,
53  int meta_iterations,
54  double deblend_factor,
55  double meta_iteration_stop,
56  size_t max_fit_size)
57  : m_least_squares_engine(least_squares_engine), m_max_iterations(max_iterations),
58  m_modified_chi_squared_scale(modified_chi_squared_scale), m_scale_factor(scale_factor),
59  m_meta_iterations(meta_iterations), m_deblend_factor(deblend_factor), m_meta_iteration_stop(meta_iteration_stop),
60  m_max_fit_size(max_fit_size * max_fit_size), m_parameters(parameters), m_frames(frames), m_priors(priors) {
61 }
62 
64 }
65 
66 namespace {
67 
68 PixelRectangle getFittingRect(SourceInterface& source, int frame_index) {
69  auto fitting_rect = source.getProperty<MeasurementFrameRectangle>(frame_index).getRect();
70 
71  if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) {
72  return PixelRectangle();
73  } else {
74  const auto& frame_info = source.getProperty<MeasurementFrameInfo>(frame_index);
75 
76  auto min = fitting_rect.getTopLeft();
77  auto max = fitting_rect.getBottomRight();
78 
79  // FIXME temporary, for now just enlarge the area by a fixed amount of pixels
80  PixelCoordinate border = (max - min) * .8 + PixelCoordinate(2, 2);
81 
82  min -= border;
83  max += border;
84 
85  // clip to image size
86  min.m_x = std::max(min.m_x, 0);
87  min.m_y = std::max(min.m_y, 0);
88  max.m_x = std::min(max.m_x, frame_info.getWidth() - 1);
89  max.m_y = std::min(max.m_y, frame_info.getHeight() - 1);
90 
91  return PixelRectangle(min, max);
92  }
93 }
94 
95 bool isFrameValid(SourceInterface& source, int frame_index) {
96  auto stamp_rect = getFittingRect(source, frame_index);
97  return stamp_rect.getWidth() > 0 && stamp_rect.getHeight() > 0;
98 }
99 
100 std::shared_ptr<VectorImage<SeFloat>> createImageCopy(SourceInterface& source, int frame_index) {
101  const auto& frame_images = source.getProperty<MeasurementFrameImages>(frame_index);
102  auto rect = getFittingRect(source, frame_index);
103 
104  auto image = VectorImage<SeFloat>::create(frame_images.getImageChunk(
105  LayerSubtractedImage, rect.getTopLeft().m_x, rect.getTopLeft().m_y, rect.getWidth(), rect.getHeight()));
106 
107  return image;
108 }
109 
110 std::shared_ptr<VectorImage<SeFloat>> createWeightImage(SourceInterface& source, int frame_index) {
111  const auto& frame_images = source.getProperty<MeasurementFrameImages>(frame_index);
112  auto frame_image = frame_images.getLockedImage(LayerSubtractedImage);
113  auto frame_image_thresholded = frame_images.getLockedImage(LayerThresholdedImage);
114  auto variance_map = frame_images.getLockedImage(LayerVarianceMap);
115 
116  const auto& frame_info = source.getProperty<MeasurementFrameInfo>(frame_index);
117  SeFloat gain = frame_info.getGain();
118  SeFloat saturation = frame_info.getSaturation();
119 
120  auto rect = getFittingRect(source, frame_index);
121  auto weight = VectorImage<SeFloat>::create(rect.getWidth(), rect.getHeight());
122 
123  for (int y = 0; y < rect.getHeight(); y++) {
124  for (int x = 0; x < rect.getWidth(); x++) {
125  auto back_var = variance_map->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
126  auto pixel_val = frame_image->getValue(rect.getTopLeft().m_x + x, rect.getTopLeft().m_y + y);
127  if (saturation > 0 && pixel_val > saturation) {
128  weight->at(x, y) = 0;
129  }
130  else if (gain > 0.0 && pixel_val > 0.0) {
131  weight->at(x, y) = sqrt(1.0 / (back_var + pixel_val / gain));
132  }
133  else {
134  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
135  }
136  }
137  }
138 
139 
140  return weight;
141 }
142 
144  SourceInterface& source, double pixel_scale, FlexibleModelFittingParameterManager& manager,
145  std::shared_ptr<FlexibleModelFittingFrame> frame, PixelRectangle stamp_rect, double down_scaling=1.0) {
146 
147  int frame_index = frame->getFrameNb();
148 
149  auto frame_coordinates = source.getProperty<MeasurementFrameCoordinates>(frame_index).getCoordinateSystem();
150  auto ref_coordinates = source.getProperty<DetectionFrameCoordinates>().getCoordinateSystem();
151 
152  auto psf_property = source.getProperty<SourcePsfProperty>(frame_index);
153  auto jacobian = source.getProperty<JacobianSource>(frame_index).asTuple();
154 
155  // The model fitting module expects to get a PSF with a pixel scale, but we have the pixel sampling step size
156  // It will be used to compute the rastering grid size, and after convolving with the PSF the result will be
157  // downscaled before copied into the frame image.
158  // We can multiply here then, as the unit is pixel/pixel, rather than "/pixel or similar
159  auto source_psf = DownSampledImagePsf(psf_property.getPixelSampling(), psf_property.getPsf(), down_scaling);
160 
161  std::vector<ConstantModel> constant_models;
162  std::vector<PointModel> point_models;
164 
165  for (auto model : frame->getModels()) {
166  model->addForSource(manager, source, constant_models, point_models, extended_models, jacobian, ref_coordinates,
167  frame_coordinates, stamp_rect.getTopLeft());
168  }
169 
170  // Full frame model with all sources
172  pixel_scale, (size_t) stamp_rect.getWidth(), (size_t) stamp_rect.getHeight(),
173  std::move(constant_models), std::move(point_models), std::move(extended_models), source_psf);
174 
175  return frame_model;
176 }
177 
178 }
179 
180 
182  FittingState fitting_state;
183 
184  //std::cout << "gs " << group.size() << "\n";
185  for (auto& source : group) {
186  SourceState initial_state;
187  initial_state.flags = Flags::NONE;
188  initial_state.iterations = 0;
189  initial_state.stop_reason = 0;
190  initial_state.reduced_chi_squared = 0.0;
191  initial_state.duration = 0.0;
192 
193  for (auto parameter : m_parameters) {
194  auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
195  if (free_parameter != nullptr) {
196  initial_state.parameters_values[free_parameter->getId()] = free_parameter->getInitialValue(source);
197  } else {
198  initial_state.parameters_values[parameter->getId()] = 0.0;
199  }
200  }
201  fitting_state.source_states.emplace_back(std::move(initial_state));
202  }
203 
204  // TODO Sort sources by flux to fit brightest sources first?
205 
206  // iterate over the whole group, fitting sources one at a time
207 
208  double prev_chi_squared = 999999.9;
209  for (int iteration = 0; iteration < m_meta_iterations; iteration++) {
210  int index = 0;
211  for (auto& source : group) {
212  fitSource(group, source, index, fitting_state);
213  index++;
214  }
215 
216  // evaluate reduced chi squared to bail out of meta iterations if no longer improving the fit
217 
218  double chi_squared = 0.0;
219  for (auto& source_state : fitting_state.source_states) {
220  chi_squared += source_state.reduced_chi_squared;
221  }
222  chi_squared /= fitting_state.source_states.size();
223 
224  if (fabs(chi_squared - prev_chi_squared) / chi_squared < m_meta_iteration_stop) {
225  break;
226  }
227 
228  prev_chi_squared = chi_squared;
229  }
230 
231 
232  // Remove parameters that couldn't be fit from the output
233 
234  for (size_t index = 0; index < group.size(); index++){
235  auto& source_state = fitting_state.source_states.at(index);
236  for (auto parameter : m_parameters) {
237  auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
238 
239  if (free_parameter != nullptr && !source_state.parameters_fitted[parameter->getId()]) {
240  source_state.parameters_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
241  source_state.parameters_sigmas[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
242  }
243  }
244  }
245 
246  // output a property for every source
247  size_t index = 0;
248  for (auto& source : group) {
249  auto& source_state = fitting_state.source_states.at(index);
250 
251  int meta_iterations = source_state.chi_squared_per_meta.size();
252  source_state.chi_squared_per_meta.resize(m_meta_iterations);
253  source_state.iterations_per_meta.resize(m_meta_iterations);
254 
255  source.setProperty<FlexibleModelFitting>(source_state.iterations, source_state.stop_reason,
256  source_state.reduced_chi_squared, source_state.duration, source_state.flags,
257  source_state.parameters_values, source_state.parameters_sigmas,
258  source_state.chi_squared_per_meta, source_state.iterations_per_meta,
259  meta_iterations);
260 
261  index++;
262  }
263 
264  updateCheckImages(group, 1.0, fitting_state);
265 }
266 
267 
268 // Used to set a dummy property in case of error that contains no result but just an error flag
270  std::unordered_map<int, double> dummy_values;
271  for (auto parameter : m_parameters) {
272  dummy_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
273  }
275  dummy_values, dummy_values,
278 }
279 
281  SourceGroupInterface& group, SourceInterface& source, int source_index,
283  int frame_index = frame->getFrameNb();
284  auto rect = getFittingRect(source, frame_index);
285 
286  double pixel_scale = 1.0;
287  FlexibleModelFittingParameterManager parameter_manager;
288  ModelFitting::EngineParameterManager engine_parameter_manager {};
289  int n_free_parameters = 0;
290 
291  int index = 0;
292  for (auto& src : group) {
293  if (index != source_index) {
294  for (auto parameter : m_parameters) {
295  auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
296 
297  if (free_parameter != nullptr) {
298  ++n_free_parameters;
299 
300  // Initial with the values from the current iteration run
301  parameter_manager.addParameter(src, parameter,
302  free_parameter->create(parameter_manager, engine_parameter_manager, src,
303  state.source_states[index].parameters_values.at(free_parameter->getId())));
304 
305  } else {
306  parameter_manager.addParameter(src, parameter,
307  parameter->create(parameter_manager, engine_parameter_manager, src));
308  }
309  }
310  }
311  index++;
312  }
313 
314  auto deblend_image = VectorImage<SeFloat>::create(rect.getWidth(), rect.getHeight());
315  index = 0;
316  for (auto& src : group) {
317  if (index != source_index) {
318  auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, rect);
319  auto final_stamp = frame_model.getImage();
320 
321  for (int y = 0; y < final_stamp->getHeight(); ++y) {
322  for (int x = 0; x < final_stamp->getWidth(); ++x) {
323  deblend_image->at(x, y) += final_stamp->at(x, y);
324  }
325  }
326  }
327  index++;
328  }
329 
330  return deblend_image;
331 }
332 
334  FlexibleModelFittingParameterManager& parameter_manager,
335  ModelFitting::EngineParameterManager& engine_parameter_manager,
336  SourceInterface& source, int index, FittingState& state) const {
337  int free_parameters_nb = 0;
338  for (auto parameter : m_parameters) {
339  auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
340 
341  if (free_parameter != nullptr) {
342  ++free_parameters_nb;
343 
344  // Initial with the values from the current iteration run
345  parameter_manager.addParameter(source, parameter,
346  free_parameter->create(parameter_manager, engine_parameter_manager, source,
347  state.source_states[index].parameters_values.at(free_parameter->getId())));
348  } else {
349  parameter_manager.addParameter(source, parameter,
350  parameter->create(parameter_manager, engine_parameter_manager, source));
351  }
352 
353  }
354 
355  // Reset access checks, as a dependent parameter could have triggered it
356  parameter_manager.clearAccessCheck();
357 
358  return free_parameters_nb;
359 }
360 
362  ResidualEstimator& res_estimator, int& good_pixels,
363  SourceGroupInterface& group, SourceInterface& source, int index, FittingState& state, double down_scaling) const {
364 
365  double pixel_scale = 1.0;
366 
367  int valid_frames = 0;
368  for (auto frame : m_frames) {
369  int frame_index = frame->getFrameNb();
370  // Validate that each frame covers the model fitting region
371  if (isFrameValid(source, frame_index)) {
372  valid_frames++;
373 
374  auto stamp_rect = getFittingRect(source, frame_index);
375  auto frame_model = createFrameModel(source, pixel_scale, parameter_manager, frame, stamp_rect, down_scaling);
376 
377  auto image = createImageCopy(source, frame_index);
378 
379  auto deblend_image = createDeblendImage(group, source, index, frame, state);
380  for (int y = 0; y < image->getHeight(); ++y) {
381  for (int x = 0; x < image->getWidth(); ++x) {
382  image->at(x, y) -= m_deblend_factor * deblend_image->at(x, y);
383  }
384  }
385 
386  auto weight = createWeightImage(source, frame_index);
387 
388  // count number of pixels that can be used for fitting
389  for (int y = 0; y < weight->getHeight(); ++y) {
390  for (int x = 0; x < weight->getWidth(); ++x) {
391  good_pixels += (weight->at(x, y) != 0.);
392  }
393  }
394 
395  // Setup residuals
396  auto data_vs_model =
397  createDataVsModelResiduals(image, std::move(frame_model), weight,
398  //LogChiSquareComparator(m_modified_chi_squared_scale));
400  res_estimator.registerBlockProvider(std::move(data_vs_model));
401  }
402  }
403 
404  return valid_frames;
405 }
406 
408  SourceGroupInterface& group, SourceInterface& source, int index, FittingState& state) const {
409 
410  double pixel_scale = 1.0;
411 
412  int total_data_points = 0;
413  SeFloat avg_reduced_chi_squared = computeChiSquared(group, source, index,pixel_scale, parameter_manager, total_data_points, state);
414 
415  int nb_of_free_parameters = 0;
416  for (auto parameter : m_parameters) {
417  bool is_free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter).get();
418  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
419  if (is_free_parameter && accessed_by_modelfitting) {
420  nb_of_free_parameters++;
421  }
422  }
423  avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
424 
425  return avg_reduced_chi_squared;
426 }
427 
429  FlexibleModelFittingParameterManager& parameter_manager, SourceInterface& source,
430  SeFloat avg_reduced_chi_squared, SeFloat duration, unsigned int iterations, unsigned int stop_reason, Flags flags,
432  int index, FittingState& state) const {
434  // Collect parameters for output
435  std::unordered_map<int, double> parameters_values, parameters_sigmas;
436  std::unordered_map<int, bool> parameters_fitted;
437 
438  for (auto parameter : m_parameters) {
439  bool is_dependent_parameter = std::dynamic_pointer_cast<FlexibleModelFittingDependentParameter>(parameter).get();
440  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
441  auto modelfitting_parameter = parameter_manager.getParameter(source, parameter);
442 
443  if (is_dependent_parameter || accessed_by_modelfitting) {
444  parameters_values[parameter->getId()] = modelfitting_parameter->getValue();
445  parameters_sigmas[parameter->getId()] = parameter->getSigma(parameter_manager, source, solution.parameter_sigmas);
446  parameters_fitted[parameter->getId()] = true;
447  } else {
448  parameters_values[parameter->getId()] = state.source_states[index].parameters_values[parameter->getId()];
449  parameters_sigmas[parameter->getId()] = state.source_states[index].parameters_sigmas[parameter->getId()];
450  parameters_fitted[parameter->getId()] = false;
451 
452  // Need to cascade the NaN to any potential dependent parameter
453  auto engine_parameter = std::dynamic_pointer_cast<EngineParameter>(modelfitting_parameter);
454  if (engine_parameter) {
455  engine_parameter->setEngineValue(std::numeric_limits<double>::quiet_NaN());
456  }
457 
458  flags |= Flags::PARTIAL_FIT;
459  }
460  }
461 
462  state.source_states[index].parameters_values = parameters_values;
463  state.source_states[index].parameters_sigmas = parameters_sigmas;
464  state.source_states[index].parameters_fitted = parameters_fitted;
465  state.source_states[index].reduced_chi_squared = avg_reduced_chi_squared;
466  state.source_states[index].chi_squared_per_meta.emplace_back(avg_reduced_chi_squared);
467  state.source_states[index].duration += duration;
468  state.source_states[index].iterations += iterations;
469  state.source_states[index].iterations_per_meta.emplace_back(iterations);
470  state.source_states[index].stop_reason = stop_reason;
471  state.source_states[index].flags = flags;
472 }
473 
475 
477  // Determine size of fitted area and if needed downsize factor
478 
479  double fit_size = 0;
480  for (auto frame : m_frames) {
481  int frame_index = frame->getFrameNb();
482  // Validate that each frame covers the model fitting region
483  if (isFrameValid(source, frame_index)) {
484  auto psf_property = source.getProperty<SourcePsfProperty>(frame_index);
485  auto stamp_rect = getFittingRect(source, frame_index);
486  fit_size = std::max(fit_size, stamp_rect.getWidth() * stamp_rect.getHeight() /
487  (psf_property.getPixelSampling() * psf_property.getPixelSampling()));
488  }
489  }
490 
491  double down_scaling = m_scale_factor;
492  if (fit_size > m_max_fit_size * 2.0) {
493  down_scaling *= sqrt(double(m_max_fit_size) / fit_size);
494  logger.warn() << "Exceeding max fit size: " << fit_size << " / " << m_max_fit_size
495  << " scaling factor: " << down_scaling;
496  }
497 
499  // Prepare parameters
500 
501  FlexibleModelFittingParameterManager parameter_manager;
502  ModelFitting::EngineParameterManager engine_parameter_manager{};
503  int n_free_parameters = fitSourcePrepareParameters(
504  parameter_manager, engine_parameter_manager, source, index, state);
505 
507  // Add models for all frames
508  ResidualEstimator res_estimator {};
509  int n_good_pixels = 0;
510  int valid_frames = fitSourcePrepareModels(
511  parameter_manager, res_estimator, n_good_pixels, group, source, index, state, down_scaling);
512 
514  // Check that we had enough data for the fit
515 
516  Flags flags = Flags::NONE;
517 
518  if (valid_frames == 0) {
519  flags = Flags::OUTSIDE;
520  }
521  else if (n_good_pixels < n_free_parameters) {
522  flags = Flags::INSUFFICIENT_DATA;
523  }
524 
525  // Do not run the model fitting for the flags above
526  if (flags != Flags::NONE) {
527  setDummyProperty(source, flags);
528  return;
529  }
530 
531  if (down_scaling < 1.0) {
532  flags |= Flags::DOWNSAMPLED;
533  }
534 
535 
537  // Add priors
538  for (auto prior : m_priors) {
539  prior->setupPrior(parameter_manager, source, res_estimator);
540  }
541 
543  // Model fitting
544 
546  auto solution = engine->solveProblem(engine_parameter_manager, res_estimator);
547 
548  auto iterations = solution.iteration_no;
549  auto stop_reason = solution.engine_stop_reason;
550  if (solution.status_flag == LeastSquareSummary::ERROR) {
551  flags |= Flags::ERROR;
552  }
553  auto duration = solution.duration;
554 
556  // compute chi squared
557 
558  SeFloat avg_reduced_chi_squared = fitSourceComputeChiSquared(parameter_manager, group, source, index, state);
559 
561  // update state with results
562  fitSourceUpdateState(parameter_manager, source, avg_reduced_chi_squared, duration, iterations, stop_reason, flags, solution,
563  index, state);
564 }
565 
567  double pixel_scale, FittingState& state) const {
568 
569  // recreate parameters
570 
571  FlexibleModelFittingParameterManager parameter_manager;
572  ModelFitting::EngineParameterManager engine_parameter_manager {};
573 
574  int index = 0;
575  for (auto& src : group) {
576  for (auto parameter : m_parameters) {
577  auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
578 
579  if (free_parameter != nullptr) {
580  // Initialize with the values from the current iteration run
581  parameter_manager.addParameter(src, parameter,
582  free_parameter->create(parameter_manager, engine_parameter_manager, src,
583  state.source_states[index].parameters_values.at(free_parameter->getId())));
584 
585  } else {
586  parameter_manager.addParameter(src, parameter,
587  parameter->create(parameter_manager, engine_parameter_manager, src));
588  }
589  }
590  index++;
591  }
592 
593  for (auto& src : group) {
594  for (auto frame : m_frames) {
595  int frame_index = frame->getFrameNb();
596 
597  if (isFrameValid(src, frame_index)) {
598  auto stamp_rect = getFittingRect(src, frame_index);
599 
600  auto frame_model = createFrameModel(src, pixel_scale, parameter_manager, frame, stamp_rect);
601  auto final_stamp = frame_model.getImage();
602 
603  auto debug_image = CheckImages::getInstance().getModelFittingImage(frame_index);
604  if (debug_image) {
605  ImageAccessor<SeFloat> debugAccessor(debug_image);
606  for (int x = 0; x < final_stamp->getWidth(); x++) {
607  for (int y = 0; y < final_stamp->getHeight(); y++) {
608  auto x_coord = stamp_rect.getTopLeft().m_x + x;
609  auto y_coord = stamp_rect.getTopLeft().m_y + y;
610  debug_image->setValue(x_coord, y_coord,
611  debugAccessor.getValue(x_coord, y_coord) + final_stamp->getValue(x, y));
612  }
613  }
614  }
615 
616  }
617  }
618  }
619 }
620 
622  std::shared_ptr<const Image<SeFloat>> model, std::shared_ptr<const Image<SeFloat>> weights, int& data_points) const {
623  double reduced_chi_squared = 0.0;
624  data_points = 0;
625 
626  ImageAccessor<SeFloat> imageAccessor(image), modelAccessor(model);
627  ImageAccessor<SeFloat> weightAccessor(weights);
628 
629  for (int y=0; y < image->getHeight(); y++) {
630  for (int x=0; x < image->getWidth(); x++) {
631  double tmp = imageAccessor.getValue(x, y) - modelAccessor.getValue(x, y);
632  reduced_chi_squared += tmp * tmp * weightAccessor.getValue(x, y) * weightAccessor.getValue(x, y);
633  if (weightAccessor.getValue(x, y) > 0) {
634  data_points++;
635  }
636  }
637  }
638  return reduced_chi_squared;
639 }
640 
642  double pixel_scale, FlexibleModelFittingParameterManager& manager, int& total_data_points, FittingState& state) const {
643  SeFloat total_chi_squared = 0;
644  total_data_points = 0;
645  int valid_frames = 0;
646  for (auto frame : m_frames) {
647  int frame_index = frame->getFrameNb();
648  // Validate that each frame covers the model fitting region
649  if (isFrameValid(source, frame_index)) {
650  valid_frames++;
651  auto stamp_rect = getFittingRect(source, frame_index);
652  auto frame_model = createFrameModel(source, pixel_scale, manager, frame, stamp_rect);
653  auto final_stamp = frame_model.getImage();
654 
655  auto image = createImageCopy(source, frame_index);
656  auto deblend_image = createDeblendImage(group, source, index, frame, state);
657  for (int y = 0; y < image->getHeight(); ++y) {
658  for (int x = 0; x < image->getWidth(); ++x) {
659  image->at(x, y) -= deblend_image->at(x, y);
660  }
661  }
662 
663  auto weight = createWeightImage(source, frame_index);
664 
665  int data_points = 0;
666  SeFloat chi_squared = computeChiSquaredForFrame(image, final_stamp, weight, data_points);
667 
668  total_data_points += data_points;
669  total_chi_squared += chi_squared;
670  }
671  }
672 
673  return total_chi_squared;
674 }
675 
676 }
677 
678 
679 
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
const double pixel_scale
Definition: TestImage.cpp:74
static Logging getLogger(const std::string &name="")
Data vs model comparator which computes a modified residual, using asinh.
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)
Provides to the LeastSquareEngine the residual values.
void registerBlockProvider(std::unique_ptr< ResidualBlockProvider > provider)
Registers a ResidualBlockProvider to the ResidualEstimator.
static CheckImages & getInstance()
Definition: CheckImages.h:138
std::shared_ptr< WriteableImage< MeasurementImage::PixelType > > getModelFittingImage(unsigned int frame_number)
void fitSourceUpdateState(FlexibleModelFittingParameterManager &parameter_manager, SourceInterface &source, SeFloat avg_reduced_chi_squared, SeFloat duration, unsigned int iterations, unsigned int stop_reason, Flags flags, ModelFitting::LeastSquareSummary solution, int index, FittingState &state) const
void fitSource(SourceGroupInterface &group, SourceInterface &source, int index, FittingState &state) const
int fitSourcePrepareParameters(FlexibleModelFittingParameterManager &parameter_manager, ModelFitting::EngineParameterManager &engine_parameter_manager, SourceInterface &source, int index, FittingState &state) const
virtual void computeProperties(SourceGroupInterface &group) const override
Computes one or more properties for the SourceGroup and/or the Sources it contains.
SeFloat fitSourceComputeChiSquared(FlexibleModelFittingParameterManager &parameter_manager, SourceGroupInterface &group, SourceInterface &source, int index, FittingState &state) const
std::vector< std::shared_ptr< FlexibleModelFittingFrame > > m_frames
FlexibleModelFittingIterativeTask(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, double scale_factor=1.0, int meta_iterations=3, double deblend_factor=1.0, double meta_iteration_stop=0.0001, size_t max_fit_size=100)
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
std::shared_ptr< VectorImage< SeFloat > > createDeblendImage(SourceGroupInterface &group, SourceInterface &source, int source_index, std::shared_ptr< FlexibleModelFittingFrame > frame, FittingState &state) const
int fitSourcePrepareModels(FlexibleModelFittingParameterManager &parameter_manager, ModelFitting::ResidualEstimator &res_estimator, int &good_pixels, SourceGroupInterface &group, SourceInterface &source, int index, FittingState &state, double downscaling) const
void setDummyProperty(SourceInterface &source, Flags flags) const
std::vector< std::shared_ptr< FlexibleModelFittingParameter > > m_parameters
SeFloat computeChiSquared(SourceGroupInterface &group, SourceInterface &source, int index, double pixel_scale, FlexibleModelFittingParameterManager &manager, int &total_data_points, FittingState &state) const
void updateCheckImages(SourceGroupInterface &group, double pixel_scale, FittingState &state) const
std::vector< std::shared_ptr< FlexibleModelFittingPrior > > m_priors
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
void addParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter, std::shared_ptr< ModelFitting::BasicParameter > engine_parameter)
bool isParamAccessed(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
Interface representing an image.
Definition: Image.h:43
Defines the interface used to group sources.
virtual unsigned int size() const =0
The SourceInterface is an abstract "source" that has properties attached to it.
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:100
T fabs(T... args)
T max(T... args)
T min(T... args)
T move(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)
Flags
Flagging of bad sources.
Definition: SourceFlags.h:37
@ DOWNSAMPLED
The fit was done on a downsampled image due to exceeding max size.
@ OUTSIDE
The object is completely outside of the measurement frame.
@ NONE
No flag is set.
@ ERROR
Error flag: something bad happened during the measurement, model fitting, etc.
@ INSUFFICIENT_DATA
There are not enough good pixels to fit the parameters.
@ PARTIAL_FIT
Some/all of the model parameters could not be fitted.
SeFloat32 SeFloat
Definition: Types.h:32
static auto logger
Definition: WCS.cpp:44
@ LayerVarianceMap
Definition: Frame.h:44
@ LayerThresholdedImage
Definition: Frame.h:40
@ LayerSubtractedImage
Definition: Frame.h:38
T quiet_NaN(T... args)
T sqrt(T... args)
Class containing the summary information of solving a least square minimization problem.
std::vector< double > parameter_sigmas
1-sigma margin of error for all the parameters