SourceXtractorPlusPlus  0.8
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>
30 #include <mutex>
31 
35 
38 
41 
54 
55 #include "ModelFitting/utils.h"
59 
61 
62 
65 
67 
75 
77 
79 
81 
82 
83 namespace SourceXtractor {
84 
85 using namespace ModelFitting;
86 
87 namespace {
88 
89 struct SourceModel {
90  double m_size;
93 
94  double exp_i0_guess;
97 
98  SourceModel(double size, double x_guess, double y_guess, double pos_range,
99  double exp_flux_guess, double exp_radius_guess, double exp_aspect_guess, double exp_rot_guess) :
100 
101  m_size(size),
102  dx(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
103  dy(std::make_shared<EngineParameter>(0, make_unique<SigmoidConverter>(-pos_range, pos_range))),
104 
105  x(createDependentParameter([x_guess](double dx) { return dx + x_guess + 0.5; }, dx)),
106  y(createDependentParameter([y_guess](double dy) { return dy + y_guess + 0.5; }, dy)),
107 
108  // FIXME
109  exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
110  moffat_i0(std::make_shared<EngineParameter>(exp_i0_guess, make_unique<ExpSigmoidConverter>(exp_i0_guess * .00001, exp_i0_guess * 1000))),
111 
112  moffat_index(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.5, 8))),
113  minkowski_exponent(std::make_shared<EngineParameter>(2, make_unique<ExpSigmoidConverter>(0.5, 10))),
114  flat_top_offset(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.000001, 10))),
115 
116  moffat_x_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
117  moffat_y_scale(std::make_shared<EngineParameter>(1, make_unique<ExpSigmoidConverter>(0.0001, 100.0))),
118  moffat_rotation(std::make_shared<EngineParameter>(-exp_rot_guess, make_unique<SigmoidConverter>(-2*M_PI, 2*M_PI)))
119  {
120  }
121 
122  void registerParameters(EngineParameterManager& manager) {
123  manager.registerParameter(dx);
124  manager.registerParameter(dy);
125 
126  manager.registerParameter(moffat_i0);
130 
134  }
135 
136  void createModels(std::vector<TransformedModel>& extended_models, std::vector<PointModel>& /*point_models*/) {
137  // Moffat model
138  {
140  auto moff = make_unique<FlattenedMoffatComponent>(moffat_i0, moffat_index, minkowski_exponent, flat_top_offset);
141  component_list.clear();
142  component_list.emplace_back(std::move(moff));
143  extended_models.emplace_back(std::move(component_list), moffat_x_scale, moffat_y_scale, moffat_rotation, m_size, m_size, x, y);
144  }
145  }
146 };
147 
148 }
149 
150 
152  auto& source_stamp = source.getProperty<DetectionFrameSourceStamp>().getStamp();
153  auto& variance_stamp = source.getProperty<DetectionFrameSourceStamp>().getVarianceStamp();
154  auto& thresholded_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdedStamp();
155  auto& threshold_map_stamp = source.getProperty<DetectionFrameSourceStamp>().getThresholdMapStamp();
156  PixelCoordinate stamp_top_left = source.getProperty<DetectionFrameSourceStamp>().getTopLeft();
157 
158  // Computes the minimum flux that a detection should have (min. detection threshold for every pixel)
159  // This will be used instead of lower or negative fluxes that can happen for various reasons
160  double min_flux = 0.;
161  auto& pixel_coordinates = source.getProperty<PixelCoordinateList>().getCoordinateList();
162  for (auto pixel : pixel_coordinates) {
163  pixel -= stamp_top_left;
164 
165  min_flux += threshold_map_stamp.getValue(pixel);
166  }
167 
168  double pixel_scale = 1;
169 
170  EngineParameterManager manager {};
171  std::vector<ConstantModel> constant_models;
172  std::vector<TransformedModel> extended_models;
173  std::vector<PointModel> point_models;
174 
175  auto& pixel_centroid = source.getProperty<PixelCentroid>();
176  auto& shape_parameters = source.getProperty<ShapeParameters>();
177  auto iso_flux = source.getProperty<IsophotalFlux>().getFlux();
178 
179  auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
180 
181  double radius_guess = shape_parameters.getEllipseA() / 2.0;
182 
183  double guess_x = pixel_centroid.getCentroidX() - stamp_top_left.m_x;
184  double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.m_y;
185 
186  double exp_flux_guess = std::max<double>(iso_flux, min_flux);
187 
188  double exp_reff_guess = radius_guess;
189  double exp_aspect_guess = std::max<double>(shape_parameters.getEllipseB() / shape_parameters.getEllipseA(), 0.01);
190  double exp_rot_guess = shape_parameters.getEllipseTheta();
191 
192  auto source_model = make_unique<SourceModel>(size, guess_x, guess_y, radius_guess * 2,
193  exp_flux_guess, exp_reff_guess, exp_aspect_guess, exp_rot_guess);
194 
195  source_model->createModels(extended_models, point_models);
196  source_model->registerParameters(manager);
197 
198  // Full frame model with all sources
200  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model {
201  pixel_scale,
202  (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(),
203  std::move(constant_models),
204  std::move(point_models),
205  std::move(extended_models)
206  };
207 
208  auto weight = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
209  std::fill(weight->getData().begin(), weight->getData().end(), 1);
210 
211  for (int y=0; y < source_stamp.getHeight(); y++) {
212  for (int x=0; x < source_stamp.getWidth(); x++) {
213  weight->at(x, y) = (thresholded_stamp.getValue(x, y) >= 0) ? 0 : 1;
214  }
215  }
216 
217  for (auto pixel : pixel_coordinates) {
218  pixel -= stamp_top_left;
219  weight->at(pixel.m_x, pixel.m_y) = 1;
220  }
221 
222  auto frame = source.getProperty<DetectionFrame>().getFrame();
223  SeFloat gain = frame->getGain();
224  SeFloat saturation = frame->getSaturation();
225 
226  for (int y=0; y < source_stamp.getHeight(); y++) {
227  for (int x=0; x < source_stamp.getWidth(); x++) {
228  auto back_var = variance_stamp.getValue(x, y);
229  if (saturation > 0 && source_stamp.getValue(x, y) >= saturation) {
230  weight->at(x, y) = 0;
231  } else if (weight->at(x, y)>0) {
232  if (gain > 0.0) {
233  weight->at(x, y) = sqrt(1.0 / (back_var + source_stamp.getValue(x, y) / gain));
234  } else {
235  weight->at(x, y) = sqrt(1.0 / back_var); // infinite gain
236  }
237  }
238  }
239  }
240 
241  // FIXME we should be able to use the source_stamp Image interface directly
242  auto image = VectorImage<SeFloat>::create(source_stamp);
243 
244  auto data_vs_model =
245  createDataVsModelResiduals(image, std::move(frame_model), weight, AsinhChiSquareComparator{});
246 
247  ResidualEstimator res_estimator {};
248  res_estimator.registerBlockProvider(move(data_vs_model));
249 
250  // Perform the minimization
251  auto engine = LeastSquareEngineManager::create(m_least_squares_engine, m_max_iterations);
252  auto solution = engine->solveProblem(manager, res_estimator);
253  size_t iterations = (size_t) boost::any_cast<std::array<double,10>>(solution.underlying_framework_info)[5];
254 
255  auto final_stamp = VectorImage<SeFloat>::create(source_stamp.getWidth(), source_stamp.getHeight());
256 
257  {
258 
259  // renders an image of the model for a single source with the final parameters
260  std::vector<TransformedModel> extended_models {};
261  std::vector<PointModel> point_models {};
262  source_model->createModels(extended_models, point_models);
263  FrameModel<NullPsf<VectorImageType>, VectorImageType> frame_model_after {
264  1, (size_t) source_stamp.getWidth(), (size_t) source_stamp.getHeight(), std::move(constant_models), std::move(point_models),
265  std::move(extended_models)
266  };
267  auto final_image = frame_model_after.getImage();
268 
269  // integrates the flux for that source
270  double total_flux = 0;
271  for (int y=0; y < source_stamp.getHeight(); y++) {
272  for (int x=0; x < source_stamp.getWidth(); x++) {
273  PixelCoordinate pixel(x, y);
274  pixel += stamp_top_left;
275 
276  // build final stamp
277  final_stamp->setValue(x, y, final_stamp->getValue(x, y) + final_image->getValue(x, y));
278 
279  total_flux += final_image->getValue(x, y);
280  }
281  }
282 
283  auto coordinate_system = source.getProperty<DetectionFrame>().getFrame()->getCoordinateSystem();
284 
285  SeFloat x = stamp_top_left.m_x + source_model->x->getValue() - 0.5f;
286  SeFloat y = stamp_top_left.m_y + source_model->y->getValue() - 0.5f;
287 
289  x, y,
290 
291  source_model->moffat_i0->getValue(),
292  source_model->moffat_index->getValue(),
293  source_model->minkowski_exponent->getValue(),
294  source_model->flat_top_offset->getValue(),
295  source_model->m_size,
296  source_model->moffat_x_scale->getValue(),
297  source_model->moffat_y_scale->getValue(),
298  source_model->moffat_rotation->getValue(),
299 
300  iterations
301  );
302 
303  }
304 }
305 
306 }
307 
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:89
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.
std::unique_ptr< T > make_unique(Args &&...args)
Definition: utils.h:32
T move(T...args)
Data vs model comparator which computes a modified residual, using asinh.
STL class.
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:72
std::shared_ptr< EngineParameter > moffat_rotation
std::shared_ptr< EngineParameter > dy
SeFloat getCentroidX() const
X coordinate of centroid.
Definition: PixelCentroid.h:48
T emplace_back(T...args)