Alexandria  2.27.0
SDC-CH common library for the Euclid project
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
31 namespace Euclid {
32 namespace Table {
33 
34 using NdArray::NdArray;
35 
36 std::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 
66  NdTypeMap{
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)}},
72  VectorTypeMap{{'B', typeid(std::vector<int32_t>)}, {'I', typeid(std::vector<int32_t>)},
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 
144 std::vector<std::string> autoDetectColumnUnits(const CCfits::Table& table_hdu) {
145  std::vector<std::string> units{};
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 
165 template <typename T>
166 std::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 
174 template <typename T>
175 std::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 
185 template <typename T>
186 std::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 
202 std::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
const std::vector< std::pair< char, std::type_index > > ScalarTypeMap
T empty(T...args)
T copy(T...args)
T front(T...args)
T to_string(T...args)
std::vector< Row::cell_type > convertVectorColumn(CCfits::Column &column, long first, long last)
T end(T...args)
STL class.
T reverse(T...args)
std::pair< std::type_index, std::size_t > asciiFormatToType(const std::string &format)
T push_back(T...args)
std::pair< std::type_index, std::size_t > binaryFormatToType(const std::string &format, const std::vector< size_t > &shape)
constexpr double s
T isdigit(T...args)
STL class.
T move(T...args)
std::map< std::string, ColumnDescription > autoDetectColumnDescriptions(std::istream &in, const std::string &comment)
Reads the column descriptions of the given stream.
std::vector< size_t > parseTDIM(const std::string &tdim)
T find_if(T...args)
T size(T...args)
T name(T...args)
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...
T begin(T...args)
T back_inserter(T...args)
const std::vector< std::pair< char, std::type_index > > NdTypeMap
std::vector< Row::cell_type > convertScalarColumn(CCfits::Column &column, long first, long last)
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 > convertNdArrayColumn(CCfits::Column &column, long first, long last)
T back(T...args)
T substr(T...args)
T transform(T...args)
const std::vector< std::pair< char, std::type_index > > VectorTypeMap
std::vector< std::string > autoDetectColumnUnits(const CCfits::Table &table_hdu)
Reads the column units based on the TUNITn keyword.
T stoi(T...args)
T reserve(T...args)
std::vector< std::string > autoDetectColumnNames(std::istream &in, const std::string &comment, size_t columns_number)
Reads the column names of the given stream.