SourceXtractorPlusPlus  0.11
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Example4.cpp
Go to the documentation of this file.
1 
23 #include <iostream>
24 #include <tuple>
25 #include <vector>
44 #include "utils.h"
46 
47 using namespace std;
48 using namespace ModelFitting;
50 
51 int main(int argc, char **argv) {
52  std::string engine_impl("levmar");
53  if (argc > 1) {
54  engine_impl = argv[1];
55  }
56 
57  auto frames_path = Elements::pathSearchInEnvVariable("multiframe.fits", "ELEMENTS_AUX_PATH");
58  auto frames = readFrames(frames_path[0].string());
59 
60  //
61  // Model creation
62  //
63  // The frame model we will use will contain a single extended model, with a
64  // single exponential component.
65 
66  // First we define the parameters of the exponential. We are going to minimize
67  // only the I0, so it is the only EngineParameter. For the engine parameters
68  // we need to use a coordinate converter. The options are:
69  // - NeutralConverter : Does no conversion
70  // - NormalizedConverter : Normalizes the parameter so the engine value is 1
71  // for a specific world value
72  // - SigmoidConverter : Converts the parameter using the sigmoid function
73  // - ExpSigmoidConverter : Converts the parameter using the exponential sigmoid function
74  auto i0 = std::make_shared<EngineParameter>(12., make_unique<ExpSigmoidConverter>(1, 100));
75  auto n = std::make_shared<ManualParameter>(1.);
76  auto k = std::make_shared<ManualParameter>(1.);
77 
78  // We create the component list of the extended model with the single exponential
79  auto reg_man = make_unique<OnlySmooth>();
80  auto exp = make_unique<SersicModelComponent>(move(reg_man), i0, n, k);
81  vector<unique_ptr<ModelComponent>> component_list {};
82  component_list.emplace_back(move(exp));
83 
84  // We create the extended model. All of its parameters will be optimized by
85  // the minimization engine.
86  auto x = std::make_shared<EngineParameter>(10, make_unique<NormalizedConverter>(30.));
87  auto y = std::make_shared<EngineParameter>(20, make_unique<NormalizedConverter>(30.));
88  auto x_scale = std::make_shared<EngineParameter>(.5, make_unique<SigmoidConverter>(0, 1));
89  auto y_scale = std::make_shared<EngineParameter>(.5, make_unique<SigmoidConverter>(0, 1));
90  auto rot_angle = std::make_shared<EngineParameter>(2., make_unique<SigmoidConverter>(0, 2*M_PI));
91 
92  // The size of the extended model (??? from the detection step ???)
93  double width = 10;
94  double height = 10;
95 
96  // We create the extended model list with a single model
98  extended_models.emplace_back(std::make_shared<ExtendedModel<cv::Mat>>(std::move(component_list), x_scale, y_scale,
99  rot_angle, width, height, x, y));
100 
101  // We read the PSF from the file
102  auto psf_path = Elements::pathSearchInEnvVariable("psf.fits", "ELEMENTS_AUX_PATH");
103  auto psf = readPsf(psf_path[0].string());
104 
105  ResidualEstimator res_estimator {};
106 
107  for (auto& pair : frames) {
108  cv::Mat image;
109  double pixel_scale {};
110  tie(image, pixel_scale) = pair;
111 
112  FrameModel<OpenCvPsf, cv::Mat> frame_model {
113  pixel_scale, (size_t)image.cols, (size_t)image.rows, {}, {},
114  move(extended_models), move(psf)
115  };
116 
117  cv::Mat weight = cv::Mat::ones(image.rows, image.cols, CV_64F);
118  auto data_vs_model = createDataVsModelResiduals(std::move(image), std::move(frame_model),
119  std::move(weight), LogChiSquareComparator{});
120  res_estimator.registerBlockProvider(move(data_vs_model));
121  }
122 
123  //
124  // Minimization
125  //
126 
127  // We print the parameters before the minimization for comparison
128  cout << "I0 (12) = " << i0->getValue() << '\n';
129  cout << "X (14.5) = " << x->getValue() << '\n';
130  cout << "Y (15.3) = " << y->getValue() << '\n';
131  cout << "X_SCALE (.83) = " << x_scale->getValue() << '\n';
132  cout << "Y_SCALE (.25) = " << y_scale->getValue() << '\n';
133  cout << "angle (2.3) = " << rot_angle->getValue() << '\n';
134 
135  // First we need to specify which parameters are optimized by the engine
136  EngineParameterManager manager {};
137  manager.registerParameter(i0);
138  manager.registerParameter(x);
139  manager.registerParameter(y);
140  manager.registerParameter(x_scale);
141  manager.registerParameter(y_scale);
142  manager.registerParameter(rot_angle);
143 
144  // Finally we create a least square engine and we solve the problem
145  std::cout << "Registered engines: " << std::endl;
146  for (auto &e : LeastSquareEngineManager::getImplementations()) {
147  std::cout << '\t' << e << std::endl;
148  }
149  std::cout << "Using engine " << engine_impl << std::endl;
150  auto engine = LeastSquareEngineManager::create(engine_impl);
151  auto t1 = chrono::steady_clock::now();
152  auto solution = engine->solveProblem(manager, res_estimator);
153  auto t2 = chrono::steady_clock::now();
154 
155  // We print the results
156  cout << "\nTime of fitting: " << chrono::duration <double, milli> (t2-t1).count() << " ms" << endl;
157  cout << "\n";
158 
159  cout << "I0 (12) = " << i0->getValue() << '\n';
160  cout << "X (14.5) = " << x->getValue() << '\n';
161  cout << "Y (15.3) = " << y->getValue() << '\n';
162  cout << "X_SCALE (.83) = " << x_scale->getValue() << '\n';
163  cout << "Y_SCALE (.25) = " << y_scale->getValue() << '\n';
164  cout << "angle (2.3) = " << rot_angle->getValue() << '\n';
165 
166  printLevmarInfo(boost::any_cast<array<double,10>>(solution.underlying_framework_info));
167 
168 }
void printLevmarInfo(std::array< double, 10 > info)
Definition: utils.h:118
T tie(T...args)
T exp(T...args)
int main()
Definition: Example1.cpp:48
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
constexpr double e
T endl(T...args)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
STL class.
std::vector< std::pair< cv::Mat, double > > readFrames(const std::string &filename)
Definition: utils.h:95
void registerParameter(std::shared_ptr< EngineParameter > parameter)
Registers an EngineParameter to the EngineParameterManager.
T move(T...args)
T count(T...args)
STL class.
ELEMENTS_API std::vector< boost::filesystem::path > pathSearchInEnvVariable(const std::string &file_name, const std::string &path_like_env_variable, SearchType search_type=SearchType::Recursive)
T make_shared(T...args)
std::unique_ptr< T > make_unique(Args &&...args)
STL class.
Class responsible for managing the parameters the least square engine minimizes.
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)
ModelFitting::OpenCvPsf readPsf(const std::string &filename)
Definition: utils.h:53
Provides to the LeastSquareEngine the residual values.
const double pixel_scale
Definition: TestImage.cpp:75
Data vs model comparator which computes a modified residual.
T emplace_back(T...args)