56 namespace SourceXtractor {
58 using namespace ModelFitting;
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]) {
74 std::cerr <<
" stopped by small gradient J^T e\n";
83 std::cerr <<
" singular matrix. Restart from current p with increased mu\n";
86 std::cerr <<
" no further error reduction is possible. Restart with increased mu\n";
89 std::cerr <<
" stopped by small ||e||_2\n";
92 std::cerr <<
" stopped by invalid (i.e. NaN or Inf) func values; a user error\n";
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";
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) {}
113 return stamp_rect.
getWidth() > 0 && stamp_rect.getHeight() > 0;
121 auto frame_image = frame->getSubtractedImage();
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);
139 auto frame_image = frame->getSubtractedImage();
140 auto frame_image_thresholded = frame->getThresholdedImage();
141 auto variance_map = frame->getVarianceMap();
145 std::fill(weight->getData().begin(), weight->getData().end(), 1);
148 SeFloat gain = measurement_frame->getGain();
149 SeFloat saturation = measurement_frame->getSaturation();
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) {
159 1.0 / (back_var + frame_image->getValue(rect.getTopLeft().m_x +
x, rect.getTopLeft().m_y +
y) / gain));
161 weight->at(
x,
y) =
sqrt(1.0 / back_var);
175 int frame_index = frame->getFrameNb();
177 auto frame_coordinates =
179 auto ref_coordinates =
190 auto group_psf =
ImagePsf(pixel_scale * psf_property.getPixelSampling(), psf_property.getPsf());
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());
205 pixel_scale, (
size_t) stamp_rect.getWidth(), (
size_t) stamp_rect.getHeight(),
216 int n_free_parameters = 0;
222 for (
auto& source : group) {
224 if (std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter)) {
228 parameter->create(parameter_manager, engine_parameter_manager, source));
239 int valid_frames = 0;
240 int n_good_pixels = 0;
242 int frame_index = frame->getFrameNb();
247 auto frame_model =
createFrameModel(group, pixel_scale, parameter_manager, frame);
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.);
262 res_estimator.registerBlockProvider(
std::move(data_vs_model));
268 if (valid_frames == 0) {
271 else if (n_good_pixels < n_free_parameters) {
278 for (
auto& source : group) {
281 auto modelfitting_parameter = parameter_manager.
getParameter(source, parameter);
283 if (manual_parameter) {
289 dummy_values, dummy_values);
295 for (
auto& source : group) {
297 prior->setupPrior(parameter_manager, source, res_estimator);
303 auto solution = engine->solveProblem(engine_parameter_manager, res_estimator);
306 int total_data_points = 0;
309 int nb_of_free_parameters = 0;
310 for (
auto& source : group) {
313 bool accessed_by_modelfitting = parameter_manager.
isParamAccessed(source, parameter);
314 if (is_free_parameter && accessed_by_modelfitting) {
315 nb_of_free_parameters++;
319 avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
322 for (
auto& source : group) {
328 bool accessed_by_modelfitting = parameter_manager.
isParamAccessed(source, parameter);
329 auto modelfitting_parameter = parameter_manager.
getParameter(source, parameter);
331 if (is_dependent_parameter || accessed_by_modelfitting) {
332 parameter_values[parameter->getId()] = modelfitting_parameter->getValue();
333 parameter_sigmas[parameter->getId()] = parameter->getSigma(parameter_manager, source, solution.parameter_sigmas);
338 if (engine_parameter) {
359 int frame_index = frame->getFrameNb();
363 auto final_stamp = frame_model.getImage();
371 for (
int x = 0;
x < final_stamp->getWidth();
x++) {
372 for (
int y = 0;
y < final_stamp->getHeight();
y++) {
373 auto x_coord = stamp_rect.getTopLeft().m_x +
x;
374 auto y_coord = stamp_rect.getTopLeft().m_y +
y;
375 debug_image->setValue(x_coord, y_coord,
376 debug_image->getValue(x_coord, y_coord) + final_stamp->getValue(
x,
y));
389 double reduced_chi_squared = 0.0;
391 for (
int y=0;
y < image->getHeight();
y++) {
392 for (
int x=0;
x < image->getWidth();
x++) {
393 double tmp = image->getValue(
x,
y) - model->getValue(
x,
y);
394 reduced_chi_squared += tmp * tmp * weights->getValue(
x,
y) * weights->getValue(
x,
y);
395 if (weights->getValue(
x,
y) > 0) {
400 return reduced_chi_squared;
407 total_data_points = 0;
408 int valid_frames = 0;
410 int frame_index = frame->getFrameNb();
415 auto final_stamp = frame_model.getImage();
423 image, final_stamp, weight, data_points);
425 total_data_points += data_points;
426 total_chi_squared += chi_squared;
430 return total_chi_squared;
EngineParameter are those derived from the minimization process.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
void setEngineValue(const double engine_value)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T dynamic_pointer_cast(T...args)
Data vs model comparator which computes a modified residual, using asinh.
void setValue(const double new_value)
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)
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 printLevmarInfo(std::array< double, 10 > info)
Provides to the LeastSquareEngine the residual values.