28 #include <boost/any.hpp>
84 namespace SourceXtractor {
86 using namespace ModelFitting;
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) :
111 exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
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))),
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)))
143 component_list.
clear();
144 component_list.emplace_back(
std::move(moff));
163 double min_flux = 0.;
165 for (
auto pixel : pixel_coordinates) {
166 pixel -= stamp_top_left;
168 min_flux += threshold_map_stamp.getValue(pixel);
182 auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
184 double radius_guess = shape_parameters.getEllipseA() / 2.0;
187 double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.
m_y;
189 double exp_flux_guess = std::max<double>(iso_flux, min_flux);
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();
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);
198 source_model->createModels(extended_models, point_models);
199 source_model->registerParameters(manager);
205 (
size_t) source_stamp.getWidth(), (
size_t) source_stamp.getHeight(),
212 std::fill(weight->getData().begin(), weight->getData().end(), 1);
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;
220 for (
auto pixel : pixel_coordinates) {
221 pixel -= stamp_top_left;
222 weight->at(pixel.m_x, pixel.m_y) = 1;
227 SeFloat saturation = detection_frame_info.getSaturation();
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) {
236 weight->at(
x,
y) =
sqrt(1.0 / (back_var + source_stamp.getValue(
x,
y) / gain));
238 weight->at(
x,
y) =
sqrt(1.0 / back_var);
255 auto solution = engine->solveProblem(manager, res_estimator);
265 source_model->createModels(extended_models, point_models);
270 auto final_image = frame_model_after.getImage();
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++) {
277 pixel += stamp_top_left;
280 final_stamp->setValue(
x,
y, final_stamp->getValue(
x,
y) + final_image->getValue(
x,
y));
282 total_flux += final_image->getValue(
x,
y);
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;
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(),
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