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