48 unsigned int max_iterations,
double modified_chi_squared_scale,
54 double deblend_factor,
55 double meta_iteration_stop,
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) {
71 if (fitting_rect.getWidth() <= 0 || fitting_rect.getHeight() <= 0) {
74 const auto& frame_info = source.
getProperty<MeasurementFrameInfo>(frame_index);
76 auto min = fitting_rect.getTopLeft();
77 auto max = fitting_rect.getBottomRight();
80 PixelCoordinate border = (max - min) * .8 + PixelCoordinate(2, 2);
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);
91 return PixelRectangle(min, max);
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;
101 const auto& frame_images = source.getProperty<MeasurementFrameImages>(frame_index);
102 auto rect = getFittingRect(source, frame_index);
105 LayerSubtractedImage, rect.getTopLeft().m_x, rect.getTopLeft().m_y, rect.getWidth(), rect.getHeight()));
111 const auto& frame_images = source.getProperty<MeasurementFrameImages>(frame_index);
116 const auto& frame_info = source.getProperty<MeasurementFrameInfo>(frame_index);
117 SeFloat gain = frame_info.getGain();
118 SeFloat saturation = frame_info.getSaturation();
120 auto rect = getFittingRect(source, frame_index);
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;
130 else if (gain > 0.0 && pixel_val > 0.0) {
131 weight->at(
x,
y) =
sqrt(1.0 / (back_var + pixel_val / gain));
134 weight->at(
x,
y) =
sqrt(1.0 / back_var);
144 SourceInterface& source,
double pixel_scale, FlexibleModelFittingParameterManager& manager,
147 int frame_index = frame->getFrameNb();
149 auto frame_coordinates = source.getProperty<MeasurementFrameCoordinates>(frame_index).getCoordinateSystem();
150 auto ref_coordinates = source.getProperty<DetectionFrameCoordinates>().getCoordinateSystem();
152 auto psf_property = source.getProperty<SourcePsfProperty>(frame_index);
153 auto jacobian = source.getProperty<JacobianSource>(frame_index).asTuple();
159 auto source_psf = DownSampledImagePsf(psf_property.getPixelSampling(), psf_property.getPsf(), down_scaling);
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());
172 pixel_scale, (
size_t) stamp_rect.getWidth(), (
size_t) stamp_rect.getHeight(),
185 for (
auto& source : group) {
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);
208 double prev_chi_squared = 999999.9;
211 for (
auto& source : group) {
212 fitSource(group, source, index, fitting_state);
218 double chi_squared = 0.0;
220 chi_squared += source_state.reduced_chi_squared;
228 prev_chi_squared = chi_squared;
234 for (
size_t index = 0; index < group.
size(); index++){
237 auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
239 if (free_parameter !=
nullptr && !source_state.parameters_fitted[parameter->getId()]) {
248 for (
auto& source : group) {
251 int meta_iterations = source_state.chi_squared_per_meta.size();
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,
275 dummy_values, dummy_values,
283 int frame_index = frame->getFrameNb();
284 auto rect = getFittingRect(source, frame_index);
289 int n_free_parameters = 0;
292 for (
auto& src : group) {
293 if (index != source_index) {
295 auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
297 if (free_parameter !=
nullptr) {
302 free_parameter->create(parameter_manager, engine_parameter_manager, src,
303 state.
source_states[index].parameters_values.at(free_parameter->getId())));
307 parameter->create(parameter_manager, engine_parameter_manager, src));
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();
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);
330 return deblend_image;
337 int free_parameters_nb = 0;
339 auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
341 if (free_parameter !=
nullptr) {
342 ++free_parameters_nb;
346 free_parameter->create(parameter_manager, engine_parameter_manager, source,
347 state.
source_states[index].parameters_values.at(free_parameter->getId())));
350 parameter->create(parameter_manager, engine_parameter_manager, source));
358 return free_parameters_nb;
367 int valid_frames = 0;
369 int frame_index = frame->getFrameNb();
371 if (isFrameValid(source, frame_index)) {
374 auto stamp_rect = getFittingRect(source, frame_index);
375 auto frame_model = createFrameModel(source,
pixel_scale, parameter_manager, frame, stamp_rect, down_scaling);
377 auto image = createImageCopy(source, frame_index);
380 for (
int y = 0;
y < image->getHeight(); ++
y) {
381 for (
int x = 0;
x < image->getWidth(); ++
x) {
386 auto weight = createWeightImage(source, frame_index);
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.);
412 int total_data_points = 0;
415 int nb_of_free_parameters = 0;
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++;
423 avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
425 return avg_reduced_chi_squared;
430 SeFloat avg_reduced_chi_squared,
SeFloat duration,
unsigned int iterations,
unsigned int stop_reason,
Flags flags,
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);
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;
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;
453 auto engine_parameter = std::dynamic_pointer_cast<EngineParameter>(modelfitting_parameter);
454 if (engine_parameter) {
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);
469 state.
source_states[index].iterations_per_meta.emplace_back(iterations);
481 int frame_index = frame->getFrameNb();
483 if (isFrameValid(source, 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()));
495 <<
" scaling factor: " << down_scaling;
504 parameter_manager, engine_parameter_manager, source, index, state);
509 int n_good_pixels = 0;
511 parameter_manager, res_estimator, n_good_pixels, group, source, index, state, down_scaling);
518 if (valid_frames == 0) {
521 else if (n_good_pixels < n_free_parameters) {
531 if (down_scaling < 1.0) {
539 prior->setupPrior(parameter_manager, source, res_estimator);
546 auto solution = engine->solveProblem(engine_parameter_manager, res_estimator);
548 auto iterations = solution.iteration_no;
549 auto stop_reason = solution.engine_stop_reason;
553 auto duration = solution.duration;
562 fitSourceUpdateState(parameter_manager, source, avg_reduced_chi_squared, duration, iterations, stop_reason, flags, solution,
575 for (
auto& src : group) {
577 auto free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter);
579 if (free_parameter !=
nullptr) {
582 free_parameter->create(parameter_manager, engine_parameter_manager, src,
583 state.
source_states[index].parameters_values.at(free_parameter->getId())));
587 parameter->create(parameter_manager, engine_parameter_manager, src));
593 for (
auto& src : group) {
595 int frame_index = frame->getFrameNb();
597 if (isFrameValid(src, frame_index)) {
598 auto stamp_rect = getFittingRect(src, frame_index);
600 auto frame_model = createFrameModel(src,
pixel_scale, parameter_manager, frame, stamp_rect);
601 auto final_stamp = frame_model.getImage();
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));
623 double reduced_chi_squared = 0.0;
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);
638 return reduced_chi_squared;
644 total_data_points = 0;
645 int valid_frames = 0;
647 int frame_index = frame->getFrameNb();
649 if (isFrameValid(source, frame_index)) {
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();
655 auto image = createImageCopy(source, frame_index);
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);
663 auto weight = createWeightImage(source, frame_index);
668 total_data_points += data_points;
669 total_chi_squared += chi_squared;
673 return total_chi_squared;
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
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.
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)
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