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