SourceXtractorPlusPlus 1.0.3
SourceXtractor++, the next generation SExtractor
Loading...
Searching...
No Matches
AssocModeConfig.cpp
Go to the documentation of this file.
1
17
18#include <map>
19#include <boost/algorithm/string.hpp>
20#include <fstream>
21
22#include <CCfits/CCfits>
23
24#include "ElementsKernel/Logging.h"
25
26#include "Table/AsciiReader.h"
27#include "Table/FitsReader.h"
28#include "Table/CastVisitor.h"
29
33
35
37
38using namespace Euclid::Configuration;
39namespace po = boost::program_options;
40
41namespace SourceXtractor {
42
44
45static const std::string ASSOC_CATALOG { "assoc-catalog" };
46static const std::string ASSOC_MODE { "assoc-mode" };
47static const std::string ASSOC_RADIUS { "assoc-radius" };
48static const std::string ASSOC_FILTER { "assoc-filter" };
49static const std::string ASSOC_COPY { "assoc-copy" };
50static const std::string ASSOC_COLUMNS { "assoc-columns" };
51static const std::string ASSOC_COORD_TYPE { "assoc-coord-type" };
52static const std::string ASSOC_SOURCE_SIZES { "assoc-source-sizes" };
53static const std::string ASSOC_SOURCE_HALF_WIDTHS { "assoc-source-half-widths" };
54static const std::string ASSOC_SOURCE_HALF_HEIGHTS { "assoc-source-half-heights" };
55static const std::string ASSOC_DEFAULT_PIXEL_SIZE { "assoc-default-pixel-size" };
56static const std::string ASSOC_GROUP_ID { "assoc-group-id" };
57static const std::string ASSOC_CONFIG { "assoc-config" };
58static const std::string ASSOC_TEST { "assoc-test" };
59
60namespace {
61
72};
73
78};
79
83};
84
85const std::string COLUMN_X = "x";
86const std::string COLUMN_Y = "y";
87const std::string COLUMN_RA = "ra";
88const std::string COLUMN_DEC = "dec";
89const std::string COLUMN_WEIGHT = "weight";
90const std::string COLUMN_GROUP_ID = "group_id";
91const std::string COLUMN_PIXEL_SIZE = "pixel_size";
92const std::string COLUMN_PIXEL_HALF_WIDTH = "pixel_half_width";
93const std::string COLUMN_PIXEL_HALF_HEIGHT = "pixel_half_height";
94
95std::vector<int> parseColumnList(const std::string& arg) {
96 if (arg.size() > 0) {
97 try {
99 boost::split(parts, arg, boost::is_any_of(","));
100
101 std::vector<int> column_list;
102 for (auto& part : parts) {
103 // the input is a 1-based index, the internal index is 0-based
104 column_list.emplace_back(boost::lexical_cast<int>(part)-1);
105 }
106 return column_list;
107 } catch(...) {
108 throw Elements::Exception() << "Can't parse column list to int: " << arg;
109 }
110 } else {
111 return {};
112 }
113}
114
115}
116
126
128 return { {"Assoc config", {
129 {ASSOC_CATALOG.c_str(), po::value<std::string>(),
130 "Assoc catalog file"},
131 {ASSOC_COLUMNS.c_str(), po::value<std::string>()->default_value("1,2"),
132 "Assoc columns to specify x/ra,y/dec[,weight] (the index of the first column is 1)"},
133 {ASSOC_MODE.c_str(), po::value<std::string>()->default_value("NEAREST"),
134 "Assoc mode [FIRST, NEAREST, MEAN, MAG_MEAN, SUM, MAG_SUM, MIN, MAX]"},
135 {ASSOC_RADIUS.c_str(), po::value<double>()->default_value(2.0),
136 "Assoc radius (Always in pixels of the detection image)"},
137 {ASSOC_FILTER.c_str(), po::value<std::string>()->default_value("ALL"),
138 "Assoc catalog filter setting: ALL, MATCHED, UNMATCHED"},
139 {ASSOC_COPY.c_str(), po::value<std::string>()->default_value(""),
140 "List of columns indices in the assoc catalog to copy on match (the index of the first column is 1). "},
141 {ASSOC_COORD_TYPE.c_str(), po::value<std::string>()->default_value("PIXEL"),
142 "Assoc coordinates type: PIXEL, WORLD"},
143 {ASSOC_SOURCE_SIZES.c_str(), po::value<int>()->default_value(-1),
144 "Column containing the source radius (in reference frame pixels)"},
145 {ASSOC_SOURCE_HALF_WIDTHS.c_str(), po::value<int>()->default_value(-1),
146 "Column containing the source half width (in reference frame pixels)"},
147 {ASSOC_SOURCE_HALF_HEIGHTS.c_str(), po::value<int>()->default_value(-1),
148 "Column containing the source half height (in reference frame pixels)"},
149 {ASSOC_DEFAULT_PIXEL_SIZE.c_str(), po::value<double>()->default_value(5.0),
150 "Default source size (in reference frame pixels)"},
151 {ASSOC_GROUP_ID.c_str(), po::value<int>()->default_value(-1),
152 "Column containing the group id"},
153 {ASSOC_CONFIG.c_str(), po::value<std::string>(),
154 "Text file containing the assoc columns configuration"},
155 {ASSOC_TEST.c_str(), po::bool_switch(),
156 "Prints the assoc configuration and quits"},
157 }}};
158}
159
161 // read configuration from command line arguments
162 readCommonConfig(args);
163
164 // read columns from columns config file
165 if (args.find(ASSOC_CONFIG) != args.end()) {
166 auto assoc_config_filename = args.at(ASSOC_CONFIG).as<std::string>();
167 readConfigFromFile(assoc_config_filename);
168 } else {
170 }
171
172 // sanity check that the configuration is coherent
173
174 // read the catalogs
175 if (m_filename != "") {
176 checkConfig();
178
179 if (args.at(ASSOC_TEST).as<bool>()) {
180 printConfig();
181 throw Elements::Exception() << "Exiting by user request";
182 }
183 }
184}
185
187 auto filter = boost::to_upper_copy(args.at(ASSOC_FILTER).as<std::string>());
188 if (assoc_filter_table.find(filter) != assoc_filter_table.end()) {
189 auto assoc_filter = assoc_filter_table.at(filter);
190 if (assoc_filter == AssocFilter::MATCHED) {
191 getDependency<PartitionStepConfig>().addPartitionStepCreator(
194 }
195 );
196 } else if (assoc_filter == AssocFilter::UNMATCHED) {
197 getDependency<PartitionStepConfig>().addPartitionStepCreator(
200 }
201 );
202 }
203 } else {
204 throw Elements::Exception() << "Invalid assoc filter: " << filter;
205 }
206
207 if (args.find(ASSOC_MODE) != args.end()) {
208 auto assoc_mode = boost::to_upper_copy(args.at(ASSOC_MODE).as<std::string>());
209 if (assoc_mode_table.find(assoc_mode) != assoc_mode_table.end()) {
210 m_assoc_mode = assoc_mode_table.at(assoc_mode);
211 } else {
212 throw Elements::Exception() << "Invalid association mode: " << assoc_mode;
213 }
214 }
215
216 m_assoc_radius = args.at(ASSOC_RADIUS).as<double>();
218
219 if (args.find(ASSOC_CATALOG) != args.end()) {
221 }
222}
223
225 m_columns = parseColumnList(args.at(ASSOC_COLUMNS).as<std::string>());
226 m_columns_idx = parseColumnList(args.at(ASSOC_COPY).as<std::string>());
227
228 m_pixel_width_column = args.at(ASSOC_SOURCE_SIZES).as<int>() - 1; // config uses 1 as first column
229 m_pixel_height_column = args.at(ASSOC_SOURCE_SIZES).as<int>() - 1; // config uses 1 as first column
230
231 if (args.find(ASSOC_SOURCE_HALF_WIDTHS) != args.end() && args.at(ASSOC_SOURCE_HALF_WIDTHS).as<int>() > 0) {
232 m_pixel_width_column = args.at(ASSOC_SOURCE_HALF_WIDTHS).as<int>() - 1; // config uses 1 as first column
233 }
234 if (args.find(ASSOC_SOURCE_HALF_HEIGHTS) != args.end() && args.at(ASSOC_SOURCE_HALF_HEIGHTS).as<int>() > 0) {
235 m_pixel_height_column = args.at(ASSOC_SOURCE_HALF_HEIGHTS).as<int>() - 1; // config uses 1 as first column
236 }
237
238 m_group_id_column = args.at(ASSOC_GROUP_ID).as<int>() - 1; // config uses 1 as first column
239
241}
242
245
246 if (m_assoc_columns.find(COLUMN_X) != m_assoc_columns.end() && m_assoc_columns.find(COLUMN_Y) != m_assoc_columns.end()) {
248
249 m_columns.push_back(m_assoc_columns.at(COLUMN_X));
250 m_assoc_columns.erase(COLUMN_X);
251 m_columns.push_back(m_assoc_columns.at(COLUMN_Y));
252 m_assoc_columns.erase(COLUMN_Y);
253
254 if (m_assoc_columns.find(COLUMN_RA) != m_assoc_columns.end() ||
255 m_assoc_columns.find(COLUMN_DEC) != m_assoc_columns.end()) {
256 throw Elements::Exception() << "Use either X/Y or RA/DEC coordinates in assoc config file but not both";
257 }
258 } else if (m_assoc_columns.find(COLUMN_RA) != m_assoc_columns.end() &&
259 m_assoc_columns.find(COLUMN_DEC) != m_assoc_columns.end()) {
260
262
263 m_columns.push_back(m_assoc_columns.at(COLUMN_RA));
264 m_assoc_columns.erase(COLUMN_RA);
265 m_columns.push_back(m_assoc_columns.at(COLUMN_DEC));
266 m_assoc_columns.erase(COLUMN_DEC);
267
268 if (m_assoc_columns.find(COLUMN_X) != m_assoc_columns.end() ||
269 m_assoc_columns.find(COLUMN_Y) != m_assoc_columns.end()) {
270 throw Elements::Exception() << "Use either X/Y or RA/DEC coordinates in assoc config file but not both";
271 }
272 } else {
273 throw Elements::Exception() << "Missing X/Y or RA/DEC coordinates in assoc config file";
274 }
275
276 if (m_assoc_columns.find(COLUMN_WEIGHT) != m_assoc_columns.end()) {
277 m_columns.push_back(m_assoc_columns.at(COLUMN_WEIGHT));
278 m_assoc_columns.erase(COLUMN_WEIGHT);
279 }
280
281 if (m_assoc_columns.find(COLUMN_PIXEL_SIZE) != m_assoc_columns.end()) {
282 m_pixel_width_column = m_assoc_columns.at(COLUMN_PIXEL_SIZE);
283 m_pixel_height_column = m_assoc_columns.at(COLUMN_PIXEL_SIZE);
284 m_assoc_columns.erase(COLUMN_PIXEL_SIZE);
285 }
286
287 if (m_assoc_columns.find(COLUMN_PIXEL_HALF_WIDTH) != m_assoc_columns.end()) {
288 m_pixel_width_column = m_assoc_columns.at(COLUMN_PIXEL_HALF_WIDTH);
289 m_assoc_columns.erase(COLUMN_PIXEL_HALF_WIDTH);
290 }
291
292 if (m_assoc_columns.find(COLUMN_PIXEL_HALF_HEIGHT) != m_assoc_columns.end()) {
293 m_pixel_height_column = m_assoc_columns.at(COLUMN_PIXEL_HALF_HEIGHT);
294 m_assoc_columns.erase(COLUMN_PIXEL_HALF_HEIGHT);
295 }
296
297 if (m_assoc_columns.find(COLUMN_GROUP_ID) != m_assoc_columns.end()) {
298 m_group_id_column = m_assoc_columns.at(COLUMN_GROUP_ID);
299 m_assoc_columns.erase(COLUMN_GROUP_ID);
300 }
301
302 for (auto& column_info : m_assoc_columns) {
303 m_custom_column_names.push_back(column_info.first);
304 m_columns_idx.push_back(column_info.second);
305 }
306}
307
309 if (m_columns.size() < 2) {
310 throw Elements::Exception() << "At least 2 columns must be specified for x,y coordinates in the assoc catalog";
311 }
312 if (m_columns.size() > 3) {
313 throw Elements::Exception() << "Maximum 3 columns for x, y and weight must be specified in the assoc catalog";
314 }
315
317 logger.warn() <<
318 "Using Assoc catalog matching in pixel coordinates with multiple detection images";
319 }
320
322 throw Elements::Exception() << "Using Assoc catalog matching in pixel coordinates without a detection images";
323 }
324
325}
326
328 AssocCoordType assoc_coord_type = AssocCoordType::PIXEL;
329 if (args.find(ASSOC_COORD_TYPE) != args.end()) {
330 auto assoc_coord_type_str = boost::to_upper_copy(args.at(ASSOC_COORD_TYPE).as<std::string>());
331 if (assoc_coord_type_table.find(assoc_coord_type_str) != assoc_coord_type_table.end()) {
332 assoc_coord_type = assoc_coord_type_table.at(assoc_coord_type_str);
333 } else {
334 throw Elements::Exception() << "Invalid association coordinate type: " << assoc_coord_type_str;
335 }
336 }
337 return assoc_coord_type;
338}
339
341 const std::vector<int>& columns, AssocCoordType assoc_coord_type) {
342 try {
344 try {
346 } catch (...) {
347 // If FITS not successful try reading as ascii
349 }
350 auto table = reader->read();
351
352 size_t exts_nb = getDependency<DetectionImageConfig>().getExtensionsNb();
353 if (exts_nb == 0) {
354 // No detection image
355 m_catalogs.emplace_back(readTable(table, columns, m_columns_idx, true));
356 } else {
357 for (size_t i = 0; i < exts_nb; i++) {
358 auto coordinate_system = getDependency<DetectionImageConfig>().getCoordinateSystem(i);
359 if (assoc_coord_type == AssocCoordType::WORLD) {
360 m_catalogs.emplace_back(readTable(table, columns, m_columns_idx, true, coordinate_system));
361 } else {
362 m_catalogs.emplace_back(readTable(table, columns, m_columns_idx, false, coordinate_system));
363 }
364 }
365 }
366 } catch (const std::exception& e) {
367 throw Elements::Exception() << "Can't either open or read assoc catalog: " << filename << " (" << e.what() << ")";
368 } catch(...) {
369 throw Elements::Exception() << "Can't either open or read assoc catalog: " << filename;
370 }
371}
372
374 const Euclid::Table::Table& table, const std::vector<int>& columns,
375 const std::vector<int>& copy_columns, bool use_world, std::shared_ptr<CoordinateSystem> coordinate_system) {
377
379 for (auto& row : table) {
380
381 ImageCoordinate coord;
382 WorldCoordinate world_coord;
383 if (use_world) {
384 world_coord = WorldCoordinate {
385 boost::apply_visitor(CastVisitor<double>{}, row[columns.at(0)]),
386 boost::apply_visitor(CastVisitor<double>{}, row[columns.at(1)]),
387 };
388 if (coordinate_system != nullptr) {
389 coord = coordinate_system->worldToImage(world_coord);
390 }
391 } else {
392 coord = ImageCoordinate {
393 // our internal pixel coordinates are zero-based
394 boost::apply_visitor(CastVisitor<double>{}, row[columns.at(0)]) - 1.0,
395 boost::apply_visitor(CastVisitor<double>{}, row[columns.at(1)]) - 1.0,
396 };
397 if (coordinate_system != nullptr) {
398 world_coord = coordinate_system->imageToWorld(coord);
399 }
400 }
401 catalog.emplace_back(CatalogEntry { coord, world_coord, 1.0, {}, 1.0, 1.0, 0 });
402 if (columns.size() == 3 && columns.at(2) >= 0) {
403 catalog.back().weight = boost::apply_visitor(CastVisitor<double>{}, row[columns.at(2)]);
404 }
405 for (auto column : copy_columns) {
406 if (column >= static_cast<int>(row.size())) {
407 throw Elements::Exception() << "Column index " << column << " is out of bounds";
408 }
409 if (row[column].type() == typeid(int)) {
410 catalog.back().assoc_columns.emplace_back(boost::get<int>(row[column]));
411 } else if (row[column].type() == typeid(double)) {
412 catalog.back().assoc_columns.emplace_back(boost::get<double>(row[column]));
413 } else if (row[column].type() == typeid(int64_t)) {
414 catalog.back().assoc_columns.emplace_back(boost::get<int64_t>(row[column]));
415 } else if (row[column].type() == typeid(SeFloat)) {
416 catalog.back().assoc_columns.emplace_back(boost::get<SeFloat>(row[column]));
417 } else {
418 throw Elements::Exception() << "Wrong type in assoc column (must be a numeric type)";
419 }
420 }
421
422 if (m_group_id_column >= 0) {
423 catalog.back().group_id = boost::apply_visitor(CastVisitor<int64_t>{}, row[m_group_id_column]);
424 }
425
427 //catalog.back().source_radius_pixels = boost::apply_visitor(CastVisitor<double>{}, row[m_pixel_size_column]);
428 catalog.back().source_pixel_width = boost::apply_visitor(CastVisitor<double>{}, row[m_pixel_width_column]);
429 catalog.back().source_pixel_height = boost::apply_visitor(CastVisitor<double>{}, row[m_pixel_height_column]);
430 } else {
431 catalog.back().source_pixel_width = m_default_pixel_size;
432 catalog.back().source_pixel_height = m_default_pixel_size;
433 }
434 }
435 return catalog;
436}
437
440
441 const std::vector<std::string> reserved_names {
442 COLUMN_X, COLUMN_Y, COLUMN_RA, COLUMN_DEC, COLUMN_WEIGHT,
443 COLUMN_GROUP_ID, COLUMN_PIXEL_SIZE, COLUMN_PIXEL_HALF_WIDTH, COLUMN_PIXEL_HALF_HEIGHT
444 };
445
446 std::ifstream config_file(filename);
447 if (!config_file.is_open()) {
448 throw Elements::Exception() << "Can't either open or read assoc config file: " << filename;
449 }
450
451 std::string line;
452 int current_column_nb = 1; // column indices start at 1
453 while (std::getline(config_file, line)) {
454 boost::trim(line);
455 // Skip lines starting with '#'
456 if (line.empty() || line[0] == '#') {
457 continue;
458 }
459
460 std::string name;
461 int number;
462 // Find the position of the '=' character
463 auto equal_sign_pos = line.find('=');
464 if (equal_sign_pos != std::string::npos) {
465 name = line.substr(0, equal_sign_pos);
466 std::string number_string = line.substr(equal_sign_pos + 1);
467 std::istringstream iss(number_string);
468 iss >> number;
469 number--;
470 current_column_nb = number + 1;
471 } else {
472 name = line;
473 number = current_column_nb++;
474 }
475 boost::trim(name);
476
477 if (std::find(reserved_names.begin(), reserved_names.end(), boost::to_lower_copy(name)) != reserved_names.end()) {
478 boost::to_lower(name);
479 }
480
481
482 // Store the parsed information into the vector
483 columns[name] = number;
484 }
485
486 config_file.close();
487 return columns;
488}
489
491 std::cout << "Assoc catalog configuration" << std::endl;
492
493 auto& catalog = m_catalogs.at(0);
494
496 std::cout << "X" << "\t";
497 std::cout << "Y" << "\t";
498 } else {
499 std::cout << "RA" << "\t";
500 std::cout << "DEC" << "\t";
501 }
502 if (m_columns.size() >= 3) {
503 std::cout << "WEIGHT" << "\t";
504 }
505
506 if (m_pixel_width_column >= 0) {
507 std::cout << "PIXEL_HALF_WIDTH" << "\t";
508 }
509 if (m_pixel_height_column >= 0) {
510 std::cout << "PIXEL_HALF_HEIGHT" << "\t";
511 }
512 if (m_group_id_column >= 0) {
513 std::cout << "GROUP_ID" << "\t";
514 }
515
516 for (const auto& name : m_custom_column_names) {
517 std::cout << name << "\t";
518 }
520
521 int lines = 0;
522 for (const auto& entry : catalog) {
524 std::cout << entry.coord.m_x << "\t";
525 std::cout << entry.coord.m_y << "\t";
526 } else {
527 std::cout << entry.world_coord.m_alpha << "\t";
528 std::cout << entry.world_coord.m_delta << "\t";
529 }
530 if (m_columns.size() >= 3) {
531 std::cout << entry.weight << "\t";
532 }
533 if (m_pixel_width_column >= 0) {
534 std::cout << entry.source_pixel_width << "\t";
535 }
536 if (m_pixel_height_column >= 0) {
537 std::cout << entry.source_pixel_height << "\t";
538 }
539 if (m_group_id_column >= 0) {
540 std::cout << entry.group_id << "\t";
541 }
542 for (auto& column : entry.assoc_columns) {
543 std::cout << column << "\t";
544 }
546
547 if (++lines >= 10) {
548 break;
549 }
550 }
551}
552
553
554}
T at(T... args)
T back(T... args)
T begin(T... args)
static Logging getLogger(const std::string &name="")
static ConfigManager & getInstance(long id)
std::map< std::string, boost::program_options::variable_value > UserValues
std::map< std::string, OptionDescriptionList > getProgramOptions() override
std::vector< CatalogEntry > readTable(const Euclid::Table::Table &table, const std::vector< int > &columns, const std::vector< int > &copy_columns, bool use_world, std::shared_ptr< CoordinateSystem > coordinate_system=nullptr)
void readConfigFromFile(const std::string &filename)
void readConfigFromParams(const UserValues &args)
void initialize(const UserValues &args) override
std::map< std::string, unsigned int > parseConfigFile(const std::string &filename)
void readCommonConfig(const UserValues &args)
std::vector< std::vector< CatalogEntry > > m_catalogs
std::map< std::string, unsigned int > m_assoc_columns
void readCatalogs(const std::string &filename, const std::vector< int > &columns, AssocCoordType assoc_coord_type)
AssocCoordType getCoordinateType(const UserValues &args) const
std::vector< std::string > m_custom_column_names
T close(T... args)
T emplace_back(T... args)
T empty(T... args)
T end(T... args)
T endl(T... args)
T find(T... args)
T getline(T... args)
T is_open(T... args)
T make_pair(T... args)
T make_shared(T... args)
static Elements::Logging logger
static const std::string ASSOC_COORD_TYPE
static const std::string ASSOC_GROUP_ID
static const std::string ASSOC_COPY
static const std::string ASSOC_CONFIG
static const std::string ASSOC_SOURCE_SIZES
static const std::string ASSOC_SOURCE_HALF_WIDTHS
static const std::string ASSOC_FILTER
static const std::string ASSOC_CATALOG
static const std::string ASSOC_SOURCE_HALF_HEIGHTS
SeFloat32 SeFloat
Definition Types.h:32
static const std::string ASSOC_MODE
static const std::string ASSOC_RADIUS
static const std::string ASSOC_TEST
static const std::string ASSOC_DEFAULT_PIXEL_SIZE
static const std::string ASSOC_COLUMNS
T size(T... args)
T substr(T... args)