Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
FitsReaderHelper.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2022 Euclid Science Ground Segment
3 *
4 * This library is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the Free
6 * Software Foundation; either version 3.0 of the License, or (at your option)
7 * any later version.
8 *
9 * This library is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU Lesser General Public License
15 * along with this library; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
25#include "FitsReaderHelper.h"
27#include <CCfits/CCfits>
28#include <boost/lexical_cast.hpp>
29#include <boost/tokenizer.hpp>
30
31namespace Euclid {
32namespace Table {
33
34using NdArray::NdArray;
35
36std::vector<std::string> autoDetectColumnNames(const CCfits::Table& table_hdu) {
38 for (int i = 1; i <= table_hdu.numCols(); ++i) {
39 std::string name = table_hdu.column(i).name();
40 if (name.empty()) {
41 name = "Col" + std::to_string(i);
42 }
43 names.push_back(std::move(name));
44 }
45 return names;
46}
47
49 if (format[0] == 'A') {
50 std::size_t size = std::stoi(format.substr(1));
51 return {typeid(std::string), size};
52 } else if (format[0] == 'I') {
53 return {typeid(int64_t), 0};
54 } else if (format[0] == 'F') {
55 return {typeid(double), 0};
56 } else if (format[0] == 'E') {
57 return {typeid(double), 0};
58 } else if (format[0] == 'D') {
59 return {typeid(double), 0};
60 }
61 throw Elements::Exception() << "FITS ASCII table format " << format << " is not "
62 << "yet supported";
63}
64
67 {'J', typeid(NdArray<int32_t>)}, {'B', typeid(NdArray<int32_t>)}, {'I', typeid(NdArray<int32_t>)},
68 {'K', typeid(NdArray<int64_t>)}, {'E', typeid(NdArray<float>)}, {'D', typeid(NdArray<double>)},
69 },
70 ScalarTypeMap{{'L', typeid(bool)}, {'J', typeid(int32_t)}, {'B', typeid(int32_t)}, {'I', typeid(int32_t)},
71 {'K', typeid(int64_t)}, {'E', typeid(float)}, {'D', typeid(double)}, {'A', typeid(std::string)}},
73 {'J', typeid(std::vector<int32_t>)}, {'K', typeid(std::vector<int64_t>)},
74 {'E', typeid(std::vector<float>)}, {'D', typeid(std::vector<double>)},
75 {'A', typeid(std::string)}};
76
78 const std::vector<size_t>& shape) {
79 // Get the size out of the format string
80 char ft = format.front();
81 int size = 1;
82 if (std::isdigit(format.front())) {
83 size = std::stoi(format.substr(0, format.size() - 1));
84 ft = format.back();
85 }
86
87 // If shape is set *and* it has more than one dimension, it is an NdArray
88 if (shape.size() > 1) {
89 auto i = std::find_if(NdTypeMap.begin(), NdTypeMap.end(),
90 [ft](const std::pair<char, std::type_index>& p) { return p.first == ft; });
91 if (i != NdTypeMap.end()) {
92 return {i->second, size};
93 }
94 }
95 // If the dimensionality is 1, it is a scalar
96 else if (size == 1) {
97 auto i = std::find_if(ScalarTypeMap.begin(), ScalarTypeMap.end(),
98 [ft](const std::pair<char, std::type_index>& p) { return p.first == ft; });
99 if (i != ScalarTypeMap.end()) {
100 return {i->second, size};
101 }
102 }
103 // Last, vectors
104 else {
105 auto i = std::find_if(VectorTypeMap.begin(), VectorTypeMap.end(),
106 [ft](const std::pair<char, std::type_index>& p) { return p.first == ft; });
107 if (i != VectorTypeMap.end()) {
108 return {i->second, size};
109 }
110 }
111 throw Elements::Exception() << "FITS binary table format " << format << " is not "
112 << "yet supported";
113}
114
116 std::vector<size_t> result{};
117 if (!tdim.empty() && tdim.front() == '(' && tdim.back() == ')') {
118 auto subtdim = tdim.substr(1, tdim.size() - 2);
119 boost::char_separator<char> sep{","};
120 boost::tokenizer<boost::char_separator<char>> tok{subtdim, sep};
121 std::transform(tok.begin(), tok.end(), std::back_inserter(result),
122 [](const std::string& s) { return boost::lexical_cast<size_t>(s); });
123 // Note: the shape is in fortran order, so we need to reverse
124 std::reverse(result.begin(), result.end());
125 }
126 return result;
127}
128
131 for (int i = 1; i <= table_hdu.numCols(); i++) {
132 auto& column = table_hdu.column(i);
133
134 if (typeid(table_hdu) == typeid(CCfits::BinTable)) {
135 column.setDimen();
136 types.push_back(binaryFormatToType(column.format(), parseTDIM(column.dimen())));
137 } else {
138 types.push_back(asciiFormatToType(column.format()));
139 }
140 }
141 return types;
142}
143
144std::vector<std::string> autoDetectColumnUnits(const CCfits::Table& table_hdu) {
146 for (int i = 1; i <= table_hdu.numCols(); ++i) {
147 units.push_back(table_hdu.column(i).unit());
148 }
149 return units;
150}
151
153 std::vector<std::string> descriptions{};
154 for (int i = 1; i <= table_hdu.numCols(); ++i) {
155 std::string desc;
156 auto key = table_hdu.keyWord().find("TDESC" + std::to_string(i));
157 if (key != table_hdu.keyWord().end()) {
158 key->second->value(desc);
159 }
160 descriptions.push_back(desc);
161 }
162 return descriptions;
163}
164
165template <typename T>
166std::vector<Row::cell_type> convertScalarColumn(CCfits::Column& column, long first, long last) {
167 std::vector<T> data;
168 column.read(data, first, last);
170 std::copy(data.begin(), data.end(), std::back_inserter(result));
171 return result;
172}
173
174template <typename T>
175std::vector<Row::cell_type> convertVectorColumn(CCfits::Column& column, long first, long last) {
177 column.readArrays(data, first, last);
179 result.reserve(data.size());
180 std::transform(data.begin(), data.end(), std::back_inserter(result),
181 [](const std::valarray<T>& valar) { return std::vector<T>(std::begin(valar), std::end(valar)); });
182 return result;
183}
184
185template <typename T>
186std::vector<Row::cell_type> convertNdArrayColumn(CCfits::Column& column, long first, long last) {
188 column.readArrays(data, first, last);
189 std::vector<size_t> shape = parseTDIM(column.dimen());
191 result.reserve(data.size());
192 std::transform(data.begin(), data.end(), std::back_inserter(result), [&shape](const std::valarray<T>& valar) {
193 return NdArray<T>(shape, std::move(std::vector<T>(std::begin(valar), std::end(valar))));
194 });
195 return result;
196}
197
199 return translateColumn(column, type, 1, column.rows());
200}
201
202std::vector<Row::cell_type> translateColumn(CCfits::Column& column, std::type_index type, long first, long last) {
203 if (type == typeid(bool)) {
204 return convertScalarColumn<bool>(column, first, last);
205 } else if (type == typeid(int32_t)) {
206 return convertScalarColumn<int32_t>(column, first, last);
207 } else if (type == typeid(int64_t)) {
208 return convertScalarColumn<int64_t>(column, first, last);
209 } else if (type == typeid(float)) {
210 return convertScalarColumn<float>(column, first, last);
211 } else if (type == typeid(double)) {
212 return convertScalarColumn<double>(column, first, last);
213 } else if (type == typeid(std::string)) {
214 return convertScalarColumn<std::string>(column, first, last);
215 } else if (type == typeid(std::vector<int32_t>)) {
216 return convertVectorColumn<int32_t>(column, first, last);
217 } else if (type == typeid(std::vector<int64_t>)) {
218 return convertVectorColumn<int64_t>(column, first, last);
219 } else if (type == typeid(std::vector<float>)) {
220 return convertVectorColumn<float>(column, first, last);
221 } else if (type == typeid(std::vector<double>)) {
222 return convertVectorColumn<double>(column, first, last);
223 } else if (type == typeid(NdArray<int32_t>)) {
224 return convertNdArrayColumn<int32_t>(column, first, last);
225 } else if (type == typeid(NdArray<int64_t>)) {
226 return convertNdArrayColumn<int64_t>(column, first, last);
227 } else if (type == typeid(NdArray<float>)) {
228 return convertNdArrayColumn<float>(column, first, last);
229 } else if (type == typeid(NdArray<double>)) {
230 return convertNdArrayColumn<double>(column, first, last);
231 }
232 throw Elements::Exception() << "Unsupported column type " << type.name();
233}
234
235} // namespace Table
236} // end of namespace Euclid
T back(T... args)
T back_inserter(T... args)
T begin(T... args)
T copy(T... args)
T empty(T... args)
T end(T... args)
T find_if(T... args)
T front(T... args)
T isdigit(T... args)
T move(T... args)
T name(T... args)
std::vector< Row::cell_type > convertVectorColumn(CCfits::Column &column, long first, long last)
std::map< std::string, ColumnDescription > autoDetectColumnDescriptions(std::istream &in, const std::string &comment)
Reads the column descriptions of the given stream.
std::vector< Row::cell_type > convertNdArrayColumn(CCfits::Column &column, long first, long last)
std::pair< std::type_index, std::size_t > binaryFormatToType(const std::string &format, const std::vector< size_t > &shape)
std::pair< std::type_index, std::size_t > asciiFormatToType(const std::string &format)
const std::vector< std::pair< char, std::type_index > > ScalarTypeMap
const std::vector< std::pair< char, std::type_index > > NdTypeMap
std::vector< std::string > autoDetectColumnUnits(const CCfits::Table &table_hdu)
Reads the column units based on the TUNITn keyword.
std::vector< std::string > autoDetectColumnNames(std::istream &in, const std::string &comment, size_t columns_number)
Reads the column names of the given stream.
std::vector< Row::cell_type > translateColumn(CCfits::Column &column, std::type_index type)
Returns a vector representing the given FITS table column data, converted to the requested type.
std::vector< size_t > parseTDIM(const std::string &tdim)
const std::vector< std::pair< char, std::type_index > > VectorTypeMap
std::vector< std::pair< std::type_index, std::size_t > > autoDetectColumnTypes(const CCfits::Table &table_hdu)
Reads the column types of the given table HDU.
std::vector< Row::cell_type > convertScalarColumn(CCfits::Column &column, long first, long last)
T push_back(T... args)
T reserve(T... args)
T reverse(T... args)
T size(T... args)
T stoi(T... args)
T substr(T... args)
T to_string(T... args)
T transform(T... args)