SourceXtractorPlusPlus  0.8
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Example_gal2.cpp
Go to the documentation of this file.
1 
23 #include <iostream>
24 #include <tuple>
25 #include <vector>
46 #include "utils.h"
48 
49 using namespace std;
50 using namespace ModelFitting;
51 
52 int main(int argc, char **argv) {
53  std::string engine_impl("levmar");
54  if (argc > 1) {
55  engine_impl = argv[1];
56  }
57 
58  // We read the image from the aux dir. Note that we will use a cv:Mat type,
59  // so the ModelFitting/Image/OpenCvMatImageTraits.h must be included.
60  cv::Mat image;
61  double pixel_scale {};
62  auto image_path = Elements::pathSearchInEnvVariable("gal.fits", "ELEMENTS_AUX_PATH");
63  tie(image, pixel_scale) = readImage(image_path[0].string());
64  size_t image_cols = image.cols;
65  size_t image_rows = image.rows;
66 
67  //
68  // Model creation
69  //
70  // The frame model we will use will contain a single extended model, with a
71  // single exponential component.
72 
73  // First we define the parameters of the exponential. We are going to minimize
74  // only the I0, so it is the only EngineParameter. For the engine parameters
75  // we need to use a coordinate converter. The options are:
76  // - NeutralConverter : Does no conversion
77  // - NormalizedConverter : Normalizes the parameter so the engine value is 1
78  // for a specific world value
79  // - SigmoidConverter : Converts the parameter using the sigmoid function
80  // - ExpSigmoidConverter : Converts the parameter using the exponential sigmoid function
81  auto i0 = std::make_shared<EngineParameter>(50000., make_unique<ExpSigmoidConverter>(1, 1000000.));
82  auto n = std::make_shared<ManualParameter>(1.);
83  auto k = std::make_shared<ManualParameter>(1.);
84 
85  // We create the component list of the extended model with the single exponential
86  auto reg_man = make_unique<OnlySmooth>();
87  auto exp = make_unique<SersicModelComponent>(move(reg_man), i0, n, k);
88  vector<unique_ptr<ModelComponent>> component_list {};
89  component_list.emplace_back(move(exp));
90 
91  // We create the extended model. All of its parameters will be optimized by
92  // the minimization engine.
93  auto x = std::make_shared<EngineParameter>(120, make_unique<NormalizedConverter>(1500.));
94  auto y = std::make_shared<EngineParameter>(140, make_unique<NormalizedConverter>(1500.));
95  auto x_scale = std::make_shared<EngineParameter>(1.0, make_unique<SigmoidConverter>(0, 10.));
96  auto y_scale = std::make_shared<EngineParameter>(1.0, make_unique<SigmoidConverter>(0, 10.));
97  auto rot_angle = std::make_shared<EngineParameter>(20.0 * M_PI/180.0, make_unique<SigmoidConverter>(0, 2*M_PI));
98 
99  // The size of the extended model (??? from the detection step ???)
100  double width = 128;
101  double height = 128;
102 
103  // We create the extended model list with a single model
104  vector<TransformedModel> extended_models {};
105  extended_models.emplace_back(std::move(component_list), x_scale, y_scale,
106  rot_angle, width, height, x, y);
107 
108  // We add a constant background
109  auto back = std::make_shared<EngineParameter>(100., make_unique<ExpSigmoidConverter>(1, 1000000.));
110  vector<ConstantModel> constant_models {};
111  constant_models.emplace_back(back);
112 
113  // We read the PSF from the file
114  auto psf_path = Elements::pathSearchInEnvVariable("psf_gal.fits", "ELEMENTS_AUX_PATH");
115  auto psf = readPsf(psf_path[0].string());
116 
117  // Finally we create the FrameModel with same pixel scale and size as the
118  // input image
119  FrameModel<OpenCvPsf, cv::Mat> frame_model {
120  pixel_scale, image_cols, image_rows, move(constant_models), {},
121  move(extended_models), move(psf)
122  };
123 
124  writeToFits(frame_model.getImage(), "example_gal2_before.fits");
125 
126  //
127  // Minimization
128  //
129 
130  // First we need to specify which parameters are optimized by the engine
131  EngineParameterManager manager {};
132  manager.registerParameter(i0);
133  manager.registerParameter(x);
134  manager.registerParameter(y);
135  manager.registerParameter(x_scale);
136  manager.registerParameter(y_scale);
137  manager.registerParameter(rot_angle);
138  manager.registerParameter(back);
139 
140  // Now we need to create the DataVsModelResiduals. We will set all the weights
141  // as ones and we will use the LogChiSquareComparator.
142  // Note that because we use cv::Mat as input we have to include the file
143  // ModelFitting/Engine/OpenCvDataVsModelInputTraits.h
144  cv::Mat weight = cv::Mat::ones(image.rows, image.cols, CV_64F);
145  auto data_vs_model = createDataVsModelResiduals(std::move(image), std::move(frame_model),
146  std::move(weight), LogChiSquareComparator{});
147 
148  // We create a residual estimator and we add our block provider
149  ResidualEstimator res_estimator {};
150  res_estimator.registerBlockProvider(move(data_vs_model));
151 
152  // Here we want to add a prior to the aspect ratio. To do that we first need
153  // A dependent parameter which computes the aspect ratio from the x and y scales
154  auto aspect_ratio = createDependentParameter(
155  [](double x, double y){return x/y;}, // This is a lambda expression computing the aspect ratio
156  x_scale, y_scale // These are the parameters the input of the lambda expression are taken
157  );
158  // Now we can create a prior to the newly created parameter
159  res_estimator.registerBlockProvider(make_unique<WorldValueResidual>(
160  aspect_ratio, // The parameter to apply the prior for
161  0.3, // The expected value. Note that this will invert the X and Y
162  1000. // The weight (optional). We give a high value to drive the minimization
163  ));
164 
165  // We print the parameters before the minimization for comparison
166  cout << "I0 = " << i0->getValue() << '\n';
167  cout << "X = " << x->getValue() << '\n';
168  cout << "Y = " << y->getValue() << '\n';
169  cout << "X_SCALE = " << x_scale->getValue() << '\n';
170  cout << "Y_SCALE = " << y_scale->getValue() << '\n';
171  cout << "angle = " << rot_angle->getValue() << '\n';
172  cout << "Background = " << back->getValue() << '\n';
173 
174  // Finally we create a levmar engine and we solve the problem
175  auto engine = LeastSquareEngineManager::create(engine_impl);
176  auto t1 = chrono::steady_clock::now();
177  auto solution = engine->solveProblem(manager, res_estimator);
178  auto t2 = chrono::steady_clock::now();
179 
180  // We print the results
181  cout << "\nTime of fitting: " << chrono::duration <double, milli> (t2-t1).count() << " ms" << endl;
182  cout << "\n";
183 
184  cout << "I0 = " << i0->getValue() << '\n';
185  cout << "X = " << x->getValue() << '\n';
186  cout << "Y = " << y->getValue() << '\n';
187  cout << "X_SCALE = " << x_scale->getValue() << '\n';
188  cout << "Y_SCALE = " << y_scale->getValue() << '\n';
189  cout << "angle = " << rot_angle->getValue() << '\n';
190  cout << "Background = " << back->getValue() << '\n';
191 
192  printLevmarInfo(boost::any_cast<array<double,10>>(solution.underlying_framework_info));
193 
194  // We create the component list of the extended model with the single exponential
195  reg_man = make_unique<OnlySmooth>();
196  exp = make_unique<SersicModelComponent>(move(reg_man), i0, n, k);
197  component_list.clear();
198  component_list.emplace_back(move(exp));
199  extended_models.clear();
200  extended_models.emplace_back(move(component_list),x_scale, y_scale,
201  rot_angle, width, height, x, y);
202  constant_models.clear();
203  constant_models.emplace_back(back);
204  FrameModel<OpenCvPsf, cv::Mat> frame_model_after {
205  pixel_scale, image_cols, image_rows, move(constant_models), {},
206  move(extended_models), readPsf(psf_path[0].string())
207  };
208  writeToFits(frame_model_after.getImage(), "example_gal2_after.fits");
209 
210 }
std::pair< cv::Mat, double > readImage(const std::string &filename)
Definition: utils.h:77
T tie(T...args)
T exp(T...args)
int main()
Definition: Example1.cpp:47
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
T endl(T...args)
ModelFitting::OpenCvPsf readPsf(const std::string &filename)
Definition: utils.h:53
std::shared_ptr< DependentParameter< Parameters...> > createDependentParameter(typename DependentParameter< Parameters...>::ValueCalculator value_calculator, Parameters...parameters)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
STL class.
void registerBlockProvider(std::unique_ptr< ResidualBlockProvider > provider)
Registers a ResidualBlockProvider to the ResidualEstimator.
void registerParameter(std::shared_ptr< EngineParameter > parameter)
Registers an EngineParameter to the EngineParameterManager.
T move(T...args)
T count(T...args)
void writeToFits(const cv::Mat &image, const std::string &filename)
Definition: utils.h:40
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)
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)
void printLevmarInfo(std::array< double, 10 > info)
Definition: utils.h:118
Provides to the LeastSquareEngine the residual values.
const double pixel_scale
Definition: TestImage.cpp:72
Data vs model comparator which computes a modified residual.
T emplace_back(T...args)