SourceXtractorPlusPlus  0.15
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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.
const PropertyType & getProperty(unsigned int index=0) const
Convenience template method to call getProperty() with a more user-friendly syntax.
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
SeFloat32 SeFloat
Definition: Types.h:32
double exp_i0_guess
static std::shared_ptr< VectorImage< T > > create(Args &&...args)
Definition: VectorImage.h:100
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< Parameters...> > createDependentParameter(typename DependentParameter< Parameters...>::ValueCalculator value_calculator, Parameters...parameters)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
std::shared_ptr< EngineParameter > minkowski_exponent
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.
std::shared_ptr< EngineParameter > moffat_i0
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< T > make_unique(Args &&...args)
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)
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 sqrt(T...args)
std::shared_ptr< EngineParameter > moffat_x_scale
T fill(T...args)
Provides to the LeastSquareEngine the residual values.
double m_size
std::shared_ptr< EngineParameter > dx
The SourceInterface is an abstract &quot;source&quot; that has properties attached to it.
std::shared_ptr< EngineParameter > flat_top_offset
const double pixel_scale
Definition: TestImage.cpp:74
std::shared_ptr< EngineParameter > moffat_rotation
std::shared_ptr< EngineParameter > dy
SeFloat getCentroidX() const
X coordinate of centroid.
Definition: PixelCentroid.h:48