68 unsigned int max_iterations,
double modified_chi_squared_scale,
73 : m_least_squares_engine(least_squares_engine),
74 m_max_iterations(max_iterations), m_modified_chi_squared_scale(modified_chi_squared_scale),
75 m_parameters(parameters), m_frames(frames), m_priors(priors), m_scale_factor(scale_factor) {}
79 return stamp_rect.
getWidth() > 0 && stamp_rect.getHeight() > 0;
102 SeFloat saturation = frame_info.getSaturation();
107 for (
int y = 0;
y < rect.getHeight();
y++) {
108 for (
int x = 0;
x < rect.getWidth();
x++) {
109 auto back_var = variance_map->getValue(rect.getTopLeft().m_x +
x, rect.getTopLeft().m_y +
y);
110 auto pixel_val = frame_image->getValue(rect.getTopLeft().m_x +
x, rect.getTopLeft().m_y +
y);
111 if (saturation > 0 && pixel_val > saturation) {
112 weight->at(
x,
y) = 0;
114 else if (gain > 0.0 && pixel_val > 0.0) {
115 weight->at(
x,
y) =
sqrt(1.0 / (back_var + pixel_val / gain));
118 weight->at(
x,
y) =
sqrt(1.0 / back_var);
131 int frame_index = frame->getFrameNb();
133 auto frame_coordinates =
135 auto ref_coordinates =
146 auto group_psf =
ImagePsf(
pixel_scale * psf_property.getPixelSampling(), psf_property.getPsf());
152 for (
auto& source : group) {
153 for (
auto model : frame->getModels()) {
154 model->addForSource(manager, source, constant_models, point_models, extended_models, jacobian, ref_coordinates, frame_coordinates,
155 stamp_rect.getTopLeft());
161 pixel_scale, (
size_t) stamp_rect.getWidth(), (
size_t) stamp_rect.getHeight(),
172 int n_free_parameters = 0;
175 for (
auto& source : group) {
177 if (std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter)) {
181 parameter->create(parameter_manager, engine_parameter_manager, source));
192 int valid_frames = 0;
193 int n_good_pixels = 0;
195 int frame_index = frame->getFrameNb();
205 for (
int y = 0;
y < weight->getHeight(); ++
y) {
206 for (
int x = 0;
x < weight->getWidth(); ++
x) {
207 n_good_pixels += (weight->at(
x,
y) != 0.);
216 res_estimator.registerBlockProvider(
std::move(data_vs_model));
222 if (valid_frames == 0) {
225 else if (n_good_pixels < n_free_parameters) {
235 for (
auto& source : group) {
237 prior->setupPrior(parameter_manager, source, res_estimator);
246 auto solution = engine->solveProblem(engine_parameter_manager, res_estimator);
247 auto iterations = solution.iteration_no;
248 auto stop_reason = solution.engine_stop_reason;
249 switch (solution.status_flag) {
260 int total_data_points = 0;
263 int nb_of_free_parameters = 0;
264 for (
auto& source : group) {
266 bool is_free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter).get();
267 bool accessed_by_modelfitting = parameter_manager.
isParamAccessed(source, parameter);
268 if (is_free_parameter && accessed_by_modelfitting) {
269 nb_of_free_parameters++;
273 avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
276 for (
auto& source : group) {
281 bool is_dependent_parameter = std::dynamic_pointer_cast<FlexibleModelFittingDependentParameter>(parameter).get();
282 bool accessed_by_modelfitting = parameter_manager.
isParamAccessed(source, parameter);
283 auto modelfitting_parameter = parameter_manager.
getParameter(source, parameter);
285 if (is_dependent_parameter || accessed_by_modelfitting) {
286 parameter_values[parameter->getId()] = modelfitting_parameter->getValue();
287 parameter_sigmas[parameter->getId()] = parameter->getSigma(parameter_manager, source, solution.parameter_sigmas);
291 auto engine_parameter = std::dynamic_pointer_cast<EngineParameter>(modelfitting_parameter);
292 if (engine_parameter) {
301 avg_reduced_chi_squared, solution.duration,
source_flags,
302 parameter_values, parameter_sigmas,
310 logger.error() <<
"An exception occured during model fitting: " <<
e.what();
319 for (
auto& source : group) {
322 auto modelfitting_parameter = parameter_manager.
getParameter(source, parameter);
323 auto manual_parameter = std::dynamic_pointer_cast<ManualParameter>(modelfitting_parameter);
324 if (manual_parameter) {
330 dummy_values, dummy_values,
340 int frame_index = frame->getFrameNb();
344 auto final_stamp = frame_model.getImage();
351 for (
int x = 0;
x < final_stamp->getWidth();
x++) {
352 for (
int y = 0;
y < final_stamp->getHeight();
y++) {
353 auto x_coord = stamp_rect.getTopLeft().m_x +
x;
354 auto y_coord = stamp_rect.getTopLeft().m_y +
y;
355 debug_image->setValue(x_coord, y_coord,
356 debugAccessor.
getValue(x_coord, y_coord) + final_stamp->getValue(
x,
y));
367 double reduced_chi_squared = 0.0;
373 for (
int y=0;
y < image->getHeight();
y++) {
374 for (
int x=0;
x < image->getWidth();
x++) {
375 double tmp = imageAccessor.getValue(
x,
y) - modelAccessor.
getValue(
x,
y);
376 reduced_chi_squared += tmp * tmp * weightAccessor.
getValue(
x,
y) * weightAccessor.
getValue(
x,
y);
382 return reduced_chi_squared;
389 total_data_points = 0;
390 int valid_frames = 0;
392 int frame_index = frame->getFrameNb();
397 auto final_stamp = frame_model.getImage();
405 image, final_stamp, weight, data_points);
407 total_data_points += data_points;
408 total_chi_squared += chi_squared;
412 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.
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)