Alexandria  2.25.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  for (auto& s : tok) {
122  result.push_back(boost::lexical_cast<size_t>(s));
123  }
124  // Note: the shape is in fortran order, so we need to reverse
125  std::reverse(result.begin(), result.end());
126  }
127  return result;
128 }
129 
132  for (int i = 1; i <= table_hdu.numCols(); i++) {
133  auto& column = table_hdu.column(i);
134 
135  if (typeid(table_hdu) == typeid(CCfits::BinTable)) {
136  column.setDimen();
137  types.push_back(binaryFormatToType(column.format(), parseTDIM(column.dimen())));
138  } else {
139  types.push_back(asciiFormatToType(column.format()));
140  }
141  }
142  return types;
143 }
144 
145 std::vector<std::string> autoDetectColumnUnits(const CCfits::Table& table_hdu) {
146  std::vector<std::string> units{};
147  for (int i = 1; i <= table_hdu.numCols(); ++i) {
148  units.push_back(table_hdu.column(i).unit());
149  }
150  return units;
151 }
152 
154  std::vector<std::string> descriptions{};
155  for (int i = 1; i <= table_hdu.numCols(); ++i) {
156  std::string desc;
157  auto key = table_hdu.keyWord().find("TDESC" + std::to_string(i));
158  if (key != table_hdu.keyWord().end()) {
159  key->second->value(desc);
160  }
161  descriptions.push_back(desc);
162  }
163  return descriptions;
164 }
165 
166 template <typename T>
167 std::vector<Row::cell_type> convertScalarColumn(CCfits::Column& column, long first, long last) {
168  std::vector<T> data;
169  column.read(data, first, last);
171  for (auto value : data) {
172  result.push_back(value);
173  }
174  return result;
175 }
176 
177 template <typename T>
178 std::vector<Row::cell_type> convertVectorColumn(CCfits::Column& column, long first, long last) {
180  column.readArrays(data, first, last);
182  for (auto& valar : data) {
183  result.push_back(std::vector<T>(std::begin(valar), std::end(valar)));
184  }
185  return result;
186 }
187 
188 template <typename T>
189 std::vector<Row::cell_type> convertNdArrayColumn(CCfits::Column& column, long first, long last) {
191  column.readArrays(data, first, last);
192  std::vector<size_t> shape = parseTDIM(column.dimen());
193 
195  for (auto& valar : data) {
196  result.push_back(NdArray<T>(shape, std::move(std::vector<T>(std::begin(valar), std::end(valar)))));
197  }
198  return result;
199 }
200 
202  return translateColumn(column, type, 1, column.rows());
203 }
204 
205 std::vector<Row::cell_type> translateColumn(CCfits::Column& column, std::type_index type, long first, long last) {
206  if (type == typeid(bool)) {
207  return convertScalarColumn<bool>(column, first, last);
208  } else if (type == typeid(int32_t)) {
209  return convertScalarColumn<int32_t>(column, first, last);
210  } else if (type == typeid(int64_t)) {
211  return convertScalarColumn<int64_t>(column, first, last);
212  } else if (type == typeid(float)) {
213  return convertScalarColumn<float>(column, first, last);
214  } else if (type == typeid(double)) {
215  return convertScalarColumn<double>(column, first, last);
216  } else if (type == typeid(std::string)) {
217  return convertScalarColumn<std::string>(column, first, last);
218  } else if (type == typeid(std::vector<int32_t>)) {
219  return convertVectorColumn<int32_t>(column, first, last);
220  } else if (type == typeid(std::vector<int64_t>)) {
221  return convertVectorColumn<int64_t>(column, first, last);
222  } else if (type == typeid(std::vector<float>)) {
223  return convertVectorColumn<float>(column, first, last);
224  } else if (type == typeid(std::vector<double>)) {
225  return convertVectorColumn<double>(column, first, last);
226  } else if (type == typeid(NdArray<int32_t>)) {
227  return convertNdArrayColumn<int32_t>(column, first, last);
228  } else if (type == typeid(NdArray<int64_t>)) {
229  return convertNdArrayColumn<int64_t>(column, first, last);
230  } else if (type == typeid(NdArray<float>)) {
231  return convertNdArrayColumn<float>(column, first, last);
232  } else if (type == typeid(NdArray<double>)) {
233  return convertNdArrayColumn<double>(column, first, last);
234  }
235  throw Elements::Exception() << "Unsupported column type " << type.name();
236 }
237 
238 } // namespace Table
239 } // end of namespace Euclid
const std::vector< std::pair< char, std::type_index > > ScalarTypeMap
T empty(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)
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)
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)
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)
std::vector< std::string > autoDetectColumnNames(std::istream &in, const std::string &comment, size_t columns_number)
Reads the column names of the given stream.