SourceXtractorPlusPlus  0.12
Please provide a description of the project.
MoffatModelFittingTask.cpp
Go to the documentation of this file.
1 
17 /*
18  * MoffatModelTask.cpp
19  *
20  * Created on: May 2, 2017
21  * Author: mschefer
22  */
23 
24 #include <iostream>
25 #include <tuple>
26 #include <vector>
27 #include <valarray>
28 #include <boost/any.hpp>
29 #include <mutex>
30 
36 
39 
42 
55 
59 
61 
62 
65 
67 
76 
78 
80 
82 
83 
84 namespace SourceXtractor {
85 
86 using namespace ModelFitting;
88 
89 namespace {
90 
91 struct SourceModel {
92  double m_size;
95 
96  double exp_i0_guess;
99 
100  SourceModel(double size, double x_guess, double y_guess, double pos_range,
101  double exp_flux_guess, double exp_radius_guess, double exp_aspect_guess, double exp_rot_guess) :
102 
103  m_size(size),
104  dx(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
105  dy(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
106 
107  x(createDependentParameter([x_guess](double dx) { return dx + x_guess + 0.5; }, dx)),
108  y(createDependentParameter([y_guess](double dy) { return dy + y_guess + 0.5; }, dy)),
109 
110  // FIXME
111  exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
112  moffat_i0(std::make_shared<EngineParameter>(exp_i0_guess, make_unique<ExpSigmoidConverter>(exp_i0_guess * .00001, exp_i0_guess * 1000))),
113 
114  moffat_index(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.5, 8))),
115  minkowski_exponent(std::make_shared<EngineParameter>(2, make_unique<ExpSigmoidConverter>(0.5, 10))),
116  flat_top_offset(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.000001, 10))),
117 
118  moffat_x_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
119  moffat_y_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
120  moffat_rotation(std::make_shared<EngineParameter>(-exp_rot_guess, make_unique<SigmoidConverter>(-2*M_PI, 2*M_PI)))
121  {
122  }
123 
124  void registerParameters(EngineParameterManager& manager) {
125  manager.registerParameter(dx);
126  manager.registerParameter(dy);
127 
128  manager.registerParameter(moffat_i0);
132 
136  }
137 
138  void createModels(std::vector<std::shared_ptr<ModelFitting::ExtendedModel<ImageInterfaceTypePtr>>>& extended_models, std::vector<PointModel>& /*point_models*/) {
139  // Moffat model
140  {
142  auto moff = Euclid::make_unique<FlattenedMoffatComponent>(moffat_i0, moffat_index, minkowski_exponent, flat_top_offset);
143  component_list.clear();
144  component_list.emplace_back(std::move(moff));
147  }
148  }
149 };
150 
151 }
152 
153 
155  auto& source_stamp = source.getProperty<DetectionFrameSourceStamp>().getStamp();
156  auto& variance_stamp = source.getProperty<DetectionFrameSourceStamp>().getVarianceStamp();
157  auto& thresholded_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdedStamp();
158  auto& threshold_map_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdMapStamp();
159  PixelCoordinate stamp_top_left = source.getProperty<DetectionFrameSourceStamp>().getTopLeft();
160 
161  // Computes the minimum flux that a detection should have (min. detection threshold for every pixel)
162  // This will be used instead of lower or negative fluxes that can happen for various reasons
163  double min_flux = 0.;
164  auto& pixel_coordinates = source.getProperty<PixelCoordinateList>().getCoordinateList();
165  for (auto pixel : pixel_coordinates) {
166  pixel -= stamp_top_left;
167 
168  min_flux += threshold_map_stamp.getValue(pixel);
169  }
170 
171  double pixel_scale = 1;
172 
173  EngineParameterManager manager {};
174  std::vector<ConstantModel> constant_models;
176  std::vector<PointModel> point_models;
177 
178  auto& pixel_centroid = source.getProperty<PixelCentroid>();
179  auto& shape_parameters = source.getProperty<ShapeParameters>();
180  auto iso_flux = source.getProperty<IsophotalFlux>().getFlux();
181 
182  auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
183 
184  double radius_guess = shape_parameters.getEllipseA() / 2.0;
185 
186  double guess_x = pixel_centroid.getCentroidX() - stamp_top_left.m_x;
187  double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.m_y;
188 
189  double exp_flux_guess = std::max<double>(iso_flux, min_flux);
190 
191  double exp_reff_guess = radius_guess;
192  double exp_aspect_guess = std::max<double>(shape_parameters.getEllipseB() / shape_parameters.getEllipseA(), 0.01);
193  double exp_rot_guess = shape_parameters.getEllipseTheta();
194 
195  auto source_model = make_unique<SourceModel>(size, guess_x, guess_y, radius_guess * 2,
196  exp_flux_guess, exp_reff_guess, exp_aspect_guess, exp_rot_guess);
197 
198  source_model->createModels(extended_models, point_models);
199  source_model->registerParameters(manager);
200 
201  // Full frame model with all sources
203  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model {
204  pixel_scale,
205  (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(),
206  std::move(constant_models),
207  std::move(point_models),
208  std::move(extended_models)
209  };
210 
211  auto weight = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
212  std::fill(weight->getData().begin(), weight->getData().end(), 1);
213 
214  for (int y=0; y < source_stamp.getHeight(); y++) {
215  for (int x=0; x < source_stamp.getWidth(); x++) {
216  weight->at(x, y) = (thresholded_stamp.getValue(x, y) >= 0) ? 0 : 1;
217  }
218  }
219 
220  for (auto pixel : pixel_coordinates) {
221  pixel -= stamp_top_left;
222  weight->at(pixel.m_x, pixel.m_y) = 1;
223  }
224 
225  const auto& detection_frame_info = source.getProperty<DetectionFrameInfo>();
226  SeFloat gain = detection_frame_info.getGain();
227  SeFloat saturation = detection_frame_info.getSaturation();
228 
229  for (int y=0; y < source_stamp.getHeight(); y++) {
230  for (int x=0; x < source_stamp.getWidth(); x++) {
231  auto back_var = variance_stamp.getValue(x, y);
232  if (saturation > 0 && source_stamp.getValue(x, y) >= saturation) {
233  weight->at(x, y) = 0;
234  } else if (weight->at(x, y)>0) {
235  if (gain > 0.0) {
236  weight->at(x, y) = sqrt(1.0 / (back_var + source_stamp.getValue(x, y) / gain));
237  } else {
238  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
239  }
240  }
241  }
242  }
243 
244  // FIXME we should be able to use the source_stamp Image interface directly
245  auto image = VectorImage<SeFloat>::create(source_stamp);
246 
247  auto data_vs_model =
248  createDataVsModelResiduals(image, std::move(frame_model), weight, AsinhChiSquareComparator{});
249 
250  ResidualEstimator res_estimator {};
251  res_estimator.registerBlockProvider(move(data_vs_model));
252 
253  // Perform the minimization
254  auto engine = LeastSquareEngineManager::create(m_least_squares_engine, m_max_iterations);
255  auto solution = engine->solveProblem(manager, res_estimator);
256  size_t iterations = (size_t) boost::any_cast<std::array<double,10>>(solution.underlying_framework_info)[5];
257 
258  auto final_stamp = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
259 
260  {
261 
262  // renders an image of the model for a single source with the final parameters
264  std::vector<PointModel> point_models {};
265  source_model->createModels(extended_models, point_models);
266  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model_after {
267  1, (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(), std::move(constant_models), std::move(point_models),
268  std::move(extended_models)
269  };
270  auto final_image = frame_model_after.getImage();
271 
272  // integrates the flux for that source
273  double total_flux = 0;
274  for (int y=0; y < source_stamp.getHeight(); y++) {
275  for (int x=0; x < source_stamp.getWidth(); x++) {
276  PixelCoordinate pixel(x, y);
277  pixel += stamp_top_left;
278 
279  // build final stamp
280  final_stamp->setValue(x, y, final_stamp->getValue(x, y) + final_image->getValue(x, y));
281 
282  total_flux += final_image->getValue(x, y);
283  }
284  }
285 
286  auto coordinate_system = source.getProperty<DetectionFrameCoordinates>().getCoordinateSystem();
287 
288  SeFloat x = stamp_top_left.m_x + source_model->x->getValue() - 0.5f;
289  SeFloat y = stamp_top_left.m_y + source_model->y->getValue() - 0.5f;
290 
292  x, y,
293 
294  source_model->moffat_i0->getValue(),
295  source_model->moffat_index->getValue(),
296  source_model->minkowski_exponent->getValue(),
297  source_model->flat_top_offset->getValue(),
298  source_model->m_size,
299  source_model->moffat_x_scale->getValue(),
300  source_model->moffat_y_scale->getValue(),
301  source_model->moffat_rotation->getValue(),
302 
303  iterations
304  );
305 
306  }
307 }
308 
309 }
310 
virtual void computeProperties(SourceInterface &source) const override
Computes one or more properties for the Source.
EngineParameter are those derived from the minimization process.
CoordinateConverter implementation using the sigmoid function.
Computes the isophotal flux and magnitude.
Definition: IsophotalFlux.h:36
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
STL namespace.
SeFloat32 SeFloat
Definition: Types.h:32
double exp_i0_guess
The centroid of all the pixels in the source, weighted by their DetectionImage pixel values...
Definition: PixelCentroid.h:37
std::shared_ptr< EngineParameter > moffat_y_scale
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
std::shared_ptr< EngineParameter > minkowski_exponent
SeFloat getCentroidX() const
X coordinate of centroid.
Definition: PixelCentroid.h:48
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > moffat_index
A copy of the rectangular region of the detection image just large enough to include the whole Source...
void registerBlockProvider(std::unique_ptr< ResidualBlockProvider > provider)
Registers a ResidualBlockProvider to the ResidualEstimator.
static std::shared_ptr< VectorImage< T > > create(Args &&... args)
Definition: VectorImage.h:89
std::shared_ptr< EngineParameter > moffat_i0
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
void registerParameter(std::shared_ptr< EngineParameter > parameter)
Registers an EngineParameter to the EngineParameterManager.
T clear(T... args)
A pixel coordinate made of two integers m_x and m_y.
T move(T... args)
Data vs model comparator which computes a modified residual, using asinh.
STL class.
T make_shared(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)
STL class.
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)
T sqrt(T... args)
std::shared_ptr< EngineParameter > moffat_x_scale
T fill(T... args)
Provides to the LeastSquareEngine the residual values.
std::shared_ptr< DependentParameter< Parameters... > > createDependentParameter(typename DependentParameter< Parameters... >::ValueCalculator value_calculator, Parameters... parameters)
double m_size
The SourceInterface is an abstract "source" that has properties attached to it.
std::shared_ptr< EngineParameter > flat_top_offset
const double pixel_scale
Definition: TestImage.cpp:75
std::shared_ptr< EngineParameter > moffat_rotation
std::shared_ptr< EngineParameter > dy
std::unique_ptr< T > make_unique(Args &&... args)