28 #include <boost/any.hpp>
83 namespace SourceXtractor {
85 using namespace ModelFitting;
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) :
109 exp_i0_guess(exp_flux_guess / (M_PI * 2.0 * 0.346 * exp_radius_guess * exp_radius_guess * exp_aspect_guess)),
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))),
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)))
141 component_list.
clear();
142 component_list.emplace_back(
std::move(moff));
160 double min_flux = 0.;
162 for (
auto pixel : pixel_coordinates) {
163 pixel -= stamp_top_left;
165 min_flux += threshold_map_stamp.getValue(pixel);
179 auto size = std::max<double>(source_stamp.getWidth(), source_stamp.getHeight());
181 double radius_guess = shape_parameters.getEllipseA() / 2.0;
184 double guess_y = pixel_centroid.getCentroidY() - stamp_top_left.
m_y;
186 double exp_flux_guess = std::max<double>(iso_flux, min_flux);
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();
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);
195 source_model->createModels(extended_models, point_models);
196 source_model->registerParameters(manager);
202 (
size_t) source_stamp.getWidth(), (
size_t) source_stamp.getHeight(),
209 std::fill(weight->getData().begin(), weight->getData().end(), 1);
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;
217 for (
auto pixel : pixel_coordinates) {
218 pixel -= stamp_top_left;
219 weight->at(pixel.m_x, pixel.m_y) = 1;
223 SeFloat gain = frame->getGain();
224 SeFloat saturation = frame->getSaturation();
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) {
233 weight->at(
x,
y) =
sqrt(1.0 / (back_var + source_stamp.getValue(
x,
y) / gain));
235 weight->at(
x,
y) =
sqrt(1.0 / back_var);
252 auto solution = engine->solveProblem(manager, res_estimator);
262 source_model->createModels(extended_models, point_models);
267 auto final_image = frame_model_after.getImage();
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++) {
274 pixel += stamp_top_left;
277 final_stamp->setValue(
x,
y, final_stamp->getValue(
x,
y) + final_image->getValue(
x,
y));
279 total_flux += final_image->getValue(
x,
y);
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;
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(),
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.
std::unique_ptr< T > make_unique(Args &&...args)
Data vs model comparator which computes a modified residual, using asinh.
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