SourceXtractorPlusPlus  0.15
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
40 
43 
47 
48 
54 
58 
60 
61 namespace SourceXtractor {
62 
63 using namespace ModelFitting;
64 
65 static auto logger = Elements::Logging::getLogger("FlexibleModelFitting");
66 
68  unsigned int max_iterations, double modified_chi_squared_scale,
72  double scale_factor)
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) {}
76 
77 bool FlexibleModelFittingTask::isFrameValid(SourceGroupInterface& group, int frame_index) const {
78  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
79  return stamp_rect.getWidth() > 0 && stamp_rect.getHeight() > 0;
80 }
81 
83  SourceGroupInterface& group, int frame_index) const {
84  const auto& frame_images = group.begin()->getProperty<MeasurementFrameImages>(frame_index);
85  auto rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
86  auto image = VectorImage<SeFloat>::create(frame_images.getImageChunk(
87  LayerSubtractedImage, rect.getTopLeft().m_x, rect.getTopLeft().m_y, rect.getWidth(), rect.getHeight()));
88 
89  return image;
90 }
91 
93  SourceGroupInterface& group, int frame_index) const {
94  const auto& frame_images = group.begin()->getProperty<MeasurementFrameImages>(frame_index);
95 
96  auto frame_image = frame_images.getLockedImage(LayerSubtractedImage);
97  auto frame_image_thresholded = frame_images.getLockedImage(LayerThresholdedImage);
98  auto variance_map = frame_images.getLockedImage(LayerVarianceMap);
99 
100  const auto& frame_info = group.begin()->getProperty<MeasurementFrameInfo>(frame_index);
101  SeFloat gain = frame_info.getGain();
102  SeFloat saturation = frame_info.getSaturation();
103 
104  auto rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
105  auto weight = VectorImage<SeFloat>::create(rect.getWidth(), rect.getHeight());
106 
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;
113  }
114  else if (gain > 0.0 && pixel_val > 0.0) {
115  weight->at(x, y) = sqrt(1.0 / (back_var + pixel_val / gain));
116  }
117  else {
118  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
119  }
120  }
121  }
122 
123  return weight;
124 }
125 
127  SourceGroupInterface& group,
130 
131  int frame_index = frame->getFrameNb();
132 
133  auto frame_coordinates =
134  group.begin()->getProperty<MeasurementFrameCoordinates>(frame_index).getCoordinateSystem();
135  auto ref_coordinates =
136  group.begin()->getProperty<DetectionFrameCoordinates>().getCoordinateSystem();
137 
138  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
139  auto psf_property = group.getProperty<PsfProperty>(frame_index);
140  auto jacobian = group.getProperty<JacobianGroup>(frame_index).asTuple();
141 
142  // The model fitting module expects to get a PSF with a pixel scale, but we have the pixel sampling step size
143  // It will be used to compute the rastering grid size, and after convolving with the PSF the result will be
144  // downscaled before copied into the frame image.
145  // We can multiply here then, as the unit is pixel/pixel, rather than "/pixel or similar
146  auto group_psf = ImagePsf(pixel_scale * psf_property.getPixelSampling(), psf_property.getPsf());
147 
148  std::vector<ConstantModel> constant_models;
149  std::vector<PointModel> point_models;
151 
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());
156  }
157  }
158 
159  // Full frame model with all sources
161  pixel_scale, (size_t) stamp_rect.getWidth(), (size_t) stamp_rect.getHeight(),
162  std::move(constant_models), std::move(point_models), std::move(extended_models), group_psf);
163 
164  return frame_model;
165 }
166 
167 
169  double pixel_scale = 1 / m_scale_factor;
170  FlexibleModelFittingParameterManager parameter_manager;
171  ModelFitting::EngineParameterManager engine_parameter_manager{};
172  int n_free_parameters = 0;
173 
174  // Prepare parameters
175  for (auto& source : group) {
176  for (auto parameter : m_parameters) {
177  if (std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter)) {
178  ++n_free_parameters;
179  }
180  parameter_manager.addParameter(source, parameter,
181  parameter->create(parameter_manager, engine_parameter_manager, source));
182  }
183  }
184 
185  try {
186  // Reset access checks, as a dependent parameter could have triggered it
187  parameter_manager.clearAccessCheck();
188 
189  // Add models for all frames
190  ResidualEstimator res_estimator{};
191 
192  int valid_frames = 0;
193  int n_good_pixels = 0;
194  for (auto frame : m_frames) {
195  int frame_index = frame->getFrameNb();
196  // Validate that each frame covers the model fitting region
197  if (isFrameValid(group, frame_index)) {
198  valid_frames++;
199 
200  auto frame_model = createFrameModel(group, pixel_scale, parameter_manager, frame);
201 
202  auto image = createImageCopy(group, frame_index);
203  auto weight = createWeightImage(group, frame_index);
204 
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.);
208  }
209  }
210 
211  // Setup residuals
212  auto data_vs_model =
213  createDataVsModelResiduals(image, std::move(frame_model), weight,
214  //LogChiSquareComparator(m_modified_chi_squared_scale));
216  res_estimator.registerBlockProvider(std::move(data_vs_model));
217  }
218  }
219 
220  // Check that we had enough data for the fit
221  Flags group_flags = Flags::NONE;
222  if (valid_frames == 0) {
223  group_flags = Flags::OUTSIDE;
224  }
225  else if (n_good_pixels < n_free_parameters) {
226  group_flags = Flags::INSUFFICIENT_DATA;
227  }
228 
229  if (group_flags != Flags::NONE) {
230  setDummyProperty(group, parameter_manager, group_flags);
231  return;
232  }
233 
234  // Add priors
235  for (auto& source : group) {
236  for (auto prior : m_priors) {
237  prior->setupPrior(parameter_manager, source, res_estimator);
238  }
239  }
240 
241  // Model fitting
242 
243  // FIXME we can no longer specify different settings with LeastSquareEngineManager!!
244  // LevmarEngine engine{m_max_iterations, 1E-3, 1E-6, 1E-6, 1E-6, 1E-4};
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  if (solution.status_flag == LeastSquareSummary::ERROR) {
250  group_flags |= Flags::ERROR;
251  }
252 
253  int total_data_points = 0;
254  SeFloat avg_reduced_chi_squared = computeChiSquared(group, pixel_scale, parameter_manager, total_data_points);
255 
256  int nb_of_free_parameters = 0;
257  for (auto& source : group) {
258  for (auto parameter : m_parameters) {
259  bool is_free_parameter = std::dynamic_pointer_cast<FlexibleModelFittingFreeParameter>(parameter).get();
260  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
261  if (is_free_parameter && accessed_by_modelfitting) {
262  nb_of_free_parameters++;
263  }
264  }
265  }
266  avg_reduced_chi_squared /= (total_data_points - nb_of_free_parameters);
267 
268  // Collect parameters for output
269  for (auto& source : group) {
270  std::unordered_map<int, double> parameter_values, parameter_sigmas;
271  auto source_flags = group_flags;
272 
273  for (auto parameter : m_parameters) {
274  bool is_dependent_parameter = std::dynamic_pointer_cast<FlexibleModelFittingDependentParameter>(parameter).get();
275  bool accessed_by_modelfitting = parameter_manager.isParamAccessed(source, parameter);
276  auto modelfitting_parameter = parameter_manager.getParameter(source, parameter);
277 
278  if (is_dependent_parameter || accessed_by_modelfitting) {
279  parameter_values[parameter->getId()] = modelfitting_parameter->getValue();
280  parameter_sigmas[parameter->getId()] = parameter->getSigma(parameter_manager, source, solution.parameter_sigmas);
281  }
282  else {
283  // Need to cascade the NaN to any potential dependent parameter
284  auto engine_parameter = std::dynamic_pointer_cast<EngineParameter>(modelfitting_parameter);
285  if (engine_parameter) {
287  }
288  parameter_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
289  parameter_sigmas[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
291  }
292  }
293  source.setProperty<FlexibleModelFitting>(iterations, stop_reason,
294  avg_reduced_chi_squared, source_flags,
295  parameter_values, parameter_sigmas);
296  }
297  updateCheckImages(group, pixel_scale, parameter_manager);
298 
299  }
300  catch (const Elements::Exception& e) {
301  logger.error() << "An exception occured during model fitting: " << e.what();
302  setDummyProperty(group, parameter_manager, Flags::ERROR);
303  }
304 }
305 
306 // Used to set a dummy property in case of error that contains no result but just an error flag
308  FlexibleModelFittingParameterManager& parameter_manager,
309  Flags flags) const {
310  for (auto& source : group) {
311  std::unordered_map<int, double> dummy_values;
312  for (auto parameter : m_parameters) {
313  auto modelfitting_parameter = parameter_manager.getParameter(source, parameter);
314  auto manual_parameter = std::dynamic_pointer_cast<ManualParameter>(modelfitting_parameter);
315  if (manual_parameter) {
317  }
318  dummy_values[parameter->getId()] = std::numeric_limits<double>::quiet_NaN();
319  }
320  source.setProperty<FlexibleModelFitting>(0, 0, std::numeric_limits<double>::quiet_NaN(), flags,
321  dummy_values, dummy_values);
322  }
323 }
324 
326  double pixel_scale, FlexibleModelFittingParameterManager& manager) const {
327 
328  int frame_id = 0;
329  for (auto frame : m_frames) {
330  int frame_index = frame->getFrameNb();
331  // Validate that each frame covers the model fitting region
332  if (isFrameValid(group, frame_index)) {
333  auto frame_model = createFrameModel(group, pixel_scale, manager, frame);
334  auto final_stamp = frame_model.getImage();
335 
336  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
337 
338  auto debug_image = CheckImages::getInstance().getModelFittingImage(frame_index);
339  if (debug_image) {
340  ImageAccessor<SeFloat> debugAccessor(debug_image);
341  for (int x = 0; x < final_stamp->getWidth(); x++) {
342  for (int y = 0; y < final_stamp->getHeight(); y++) {
343  auto x_coord = stamp_rect.getTopLeft().m_x + x;
344  auto y_coord = stamp_rect.getTopLeft().m_y + y;
345  debug_image->setValue(x_coord, y_coord,
346  debugAccessor.getValue(x_coord, y_coord) + final_stamp->getValue(x, y));
347  }
348  }
349  }
350  }
351  frame_id++;
352  }
353 }
354 
356  std::shared_ptr<const Image<SeFloat>> model, std::shared_ptr<const Image<SeFloat>> weights, int& data_points) const {
357  double reduced_chi_squared = 0.0;
358  data_points = 0;
359 
360  ImageAccessor<SeFloat> imageAccessor(image), modelAccessor(model);
361  ImageAccessor<SeFloat> weightAccessor(weights);
362 
363  for (int y=0; y < image->getHeight(); y++) {
364  for (int x=0; x < image->getWidth(); x++) {
365  double tmp = imageAccessor.getValue(x, y) - modelAccessor.getValue(x, y);
366  reduced_chi_squared += tmp * tmp * weightAccessor.getValue(x, y) * weightAccessor.getValue(x, y);
367  if (weightAccessor.getValue(x, y) > 0) {
368  data_points++;
369  }
370  }
371  }
372  return reduced_chi_squared;
373 }
374 
376  double pixel_scale, FlexibleModelFittingParameterManager& manager, int& total_data_points) const {
377 
378  SeFloat total_chi_squared = 0;
379  total_data_points = 0;
380  int valid_frames = 0;
381  for (auto frame : m_frames) {
382  int frame_index = frame->getFrameNb();
383  // Validate that each frame covers the model fitting region
384  if (isFrameValid(group, frame_index)) {
385  valid_frames++;
386  auto frame_model = createFrameModel(group, pixel_scale, manager, frame);
387  auto final_stamp = frame_model.getImage();
388 
389  auto stamp_rect = group.getProperty<MeasurementFrameGroupRectangle>(frame_index);
390  auto image = createImageCopy(group, frame_index);
391  auto weight = createWeightImage(group, frame_index);
392 
393  int data_points = 0;
394  SeFloat chi_squared = computeChiSquaredForFrame(
395  image, final_stamp, weight, data_points);
396 
397  total_data_points += data_points;
398  total_chi_squared += chi_squared;
399  }
400  }
401 
402  return total_chi_squared;
403 }
404 
406 }
407 
408 }
static StaticPlugin< SourceFlagsPlugin > source_flags
void addParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter, std::shared_ptr< ModelFitting::BasicParameter > engine_parameter)
static auto logger
Definition: WCS.cpp:46
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
void updateCheckImages(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager) const
EngineParameter are those derived from the minimization process.
The object is completely outside of the measurement frame.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
constexpr double e
void setEngineValue(const double engine_value)
virtual void computeProperties(SourceGroupInterface &group) const override
Computes one or more properties for the SourceGroup and/or the Sources it contains.
std::shared_ptr< VectorImage< SeFloat > > createWeightImage(SourceGroupInterface &group, int frame_index) const
SeFloat32 SeFloat
Definition: Types.h:32
static std::shared_ptr< VectorImage< T > > create(Args &&...args)
Definition: VectorImage.h:100
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
STL class.
static CheckImages & getInstance()
Definition: CheckImages.h:136
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
T dynamic_pointer_cast(T...args)
T move(T...args)
Data vs model comparator which computes a modified residual, using asinh.
bool isParamAccessed(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
void setValue(const double new_value)
Defines the interface used to group sources.
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, double scale_factor=1.0)
std::vector< std::shared_ptr< FlexibleModelFittingPrior > > m_priors
std::vector< std::shared_ptr< FlexibleModelFittingFrame > > m_frames
bool isFrameValid(SourceGroupInterface &group, int frame_index) const
STL class.
There are not enough good pixels to fit the parameters.
std::shared_ptr< ModelFitting::BasicParameter > getParameter(const SourceInterface &source, std::shared_ptr< const FlexibleModelFittingParameter > parameter) const
SeFloat computeChiSquared(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager, int &total_data_points) const
Some/all of the model parameters could not be fitted.
Flags
Flagging of bad sources.
Definition: SourceFlags.h:37
Interface representing an image.
Definition: Image.h:43
std::shared_ptr< VectorImage< SeFloat > > createImageCopy(SourceGroupInterface &group, int frame_index) const
Class responsible for managing the parameters the least square engine minimizes.
void setDummyProperty(SourceGroupInterface &group, FlexibleModelFittingParameterManager &parameter_manager, Flags flags) const
static std::shared_ptr< LeastSquareEngine > create(const std::string &name, unsigned max_iterations=1000)
std::vector< std::shared_ptr< FlexibleModelFittingParameter > > m_parameters
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)
T quiet_NaN(T...args)
T sqrt(T...args)
Provides to the LeastSquareEngine the residual values.
const char * what() const noexceptoverride
static Logging getLogger(const std::string &name="")
const double pixel_scale
Definition: TestImage.cpp:74
ModelFitting::FrameModel< ImagePsf, std::shared_ptr< VectorImage< SourceXtractor::SeFloat > > > createFrameModel(SourceGroupInterface &group, double pixel_scale, FlexibleModelFittingParameterManager &manager, std::shared_ptr< FlexibleModelFittingFrame > frame) const
Error flag: something bad happened during the measurement, model fitting, etc.
std::shared_ptr< WriteableImage< MeasurementImage::PixelType > > getModelFittingImage(unsigned int frame_number)
std::shared_ptr< ImageAccessor< SeFloat > > getLockedImage(FrameImageLayer layer) const