28 #include <boost/any.hpp>
83 namespace SourceXtractor {
85 using namespace ModelFitting;
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) :
110 exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
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))),
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)))
142 component_list.
clear();
143 component_list.emplace_back(
std::move(moff));
162 double min_flux = 0.;
164 for (
auto pixel : pixel_coordinates) {
165 pixel -= stamp_top_left;
167 min_flux += threshold_map_stamp.getValue(pixel);
181 auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
183 double radius_guess = shape_parameters.getEllipseA() / 2.0;
186 double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.
m_y;
188 double exp_flux_guess = std::max<double>(iso_flux, min_flux);
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();
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);
197 source_model->createModels(extended_models, point_models);
198 source_model->registerParameters(manager);
204 (
size_t) source_stamp.getWidth(), (
size_t) source_stamp.getHeight(),
211 std::fill(weight->getData().begin(), weight->getData().end(), 1);
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;
219 for (
auto pixel : pixel_coordinates) {
220 pixel -= stamp_top_left;
221 weight->at(pixel.m_x, pixel.m_y) = 1;
225 SeFloat gain = frame->getGain();
226 SeFloat saturation = frame->getSaturation();
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) {
235 weight->at(
x,
y) =
sqrt(1.0 / (back_var + source_stamp.getValue(
x,
y) / gain));
237 weight->at(
x,
y) =
sqrt(1.0 / back_var);
254 auto solution = engine->solveProblem(manager, res_estimator);
264 source_model->createModels(extended_models, point_models);
269 auto final_image = frame_model_after.getImage();
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++) {
276 pixel += stamp_top_left;
279 final_stamp->setValue(
x,
y, final_stamp->getValue(
x,
y) + final_image->getValue(
x,
y));
281 total_flux += final_image->getValue(
x,
y);
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;
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(),
EngineParameter are those derived from the minimization process.
CoordinateConverter implementation using the sigmoid function.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
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
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.
Data vs model comparator which computes a modified residual, using asinh.
std::unique_ptr< T > make_unique(Args &&...args)
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)
std::shared_ptr< EngineParameter > moffat_x_scale
Provides to the LeastSquareEngine the residual values.
std::shared_ptr< EngineParameter > dx
std::shared_ptr< EngineParameter > flat_top_offset
std::shared_ptr< EngineParameter > moffat_rotation
std::shared_ptr< EngineParameter > dy