Alexandria  2.16
Please provide a description of the project.
FitsReaderHelper.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2020 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 <CCfits/CCfits>
26 #include <boost/lexical_cast.hpp>
27 #include <boost/tokenizer.hpp>
29 #include "FitsReaderHelper.h"
30 
31 namespace Euclid {
32 namespace Table {
33 
34 using NdArray::NdArray;
35 
36 std::vector<std::string> autoDetectColumnNames(const CCfits::Table& table_hdu) {
37  std::vector<std::string> names {};
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 
48 std::type_index asciiFormatToType(const std::string& format) {
49  if (format[0] == 'A') {
50  return typeid(std::string);
51  } else if (format[0] == 'I') {
52  return typeid(int64_t);
53  } else if (format[0] == 'F') {
54  return typeid(double);
55  } else if (format[0] == 'E') {
56  return typeid(double);
57  } else if (format[0] == 'D') {
58  return typeid(double);
59  }
60  throw Elements::Exception() << "FITS ASCII table format " << format << " is not "
61  << "yet supported";
62 }
63 
64 std::type_index binaryFormatToType(const std::string& format, const std::vector<size_t>& shape) {
65  // Get the size out of the format string
66  char ft = format.front();
67  int size = 1;
68  if (std::isdigit(format.front())) {
69  size = std::stoi(format.substr(0, format.size() - 1));
70  ft = format.back();
71  }
72 
73  // If shape is set, it is an NdArray
74  if (!shape.empty()) {
75  if (ft == 'B') {
76  return typeid(NdArray<int32_t>);
77  } else if (ft == 'I') {
78  return typeid(NdArray<int32_t>);
79  } else if (ft == 'J') {
80  return typeid(NdArray<int32_t>);
81  } else if (ft == 'K') {
82  return typeid(NdArray<int64_t>);
83  } else if (ft == 'E') {
84  return typeid(NdArray<float>);
85  } else if (ft == 'D') {
86  return typeid(NdArray<double>);
87  }
88  }
89  // If the dimensionality is 1, it is a scalar
90  else if (size == 1) {
91  if (ft == 'L') {
92  return typeid(bool);
93  }
94  else if (ft == 'B') {
95  return typeid(int32_t);
96  }
97  else if (ft == 'I') {
98  return typeid(int32_t);
99  }
100  else if (ft == 'J') {
101  return typeid(int32_t);
102  }
103  else if (ft == 'K') {
104  return typeid(int64_t);
105  }
106  else if (ft == 'E') {
107  return typeid(float);
108  }
109  else if (ft == 'D') {
110  return typeid(double);
111  }
112  }
113  // Last, vectors
114  else {
115  if (ft == 'B') {
116  return typeid(std::vector<int32_t>);
117  } else if (ft == 'I') {
118  return typeid(std::vector<int32_t>);
119  } else if (ft == 'J') {
120  return typeid(std::vector<int32_t>);
121  } else if (ft == 'K') {
122  return typeid(std::vector<int64_t>);
123  } else if (ft == 'E') {
124  return typeid(std::vector<float>);
125  } else if (ft == 'D') {
126  return typeid(std::vector<double>);
127  } else if (ft == 'A') {
128  return typeid(std::string);
129  }
130  }
131  throw Elements::Exception() << "FITS binary table format " << format << " is not "
132  << "yet supported";
133 }
134 
135 std::vector<size_t> parseTDIM(const std::string &tdim) {
136  std::vector<size_t> result {};
137  if (!tdim.empty() && tdim.front() == '(' && tdim.back() == ')') {
138  auto subtdim = tdim.substr(1, tdim.size() - 2);
139  boost::char_separator<char> sep{","};
140  boost::tokenizer<boost::char_separator<char>> tok{subtdim, sep};
141  for (auto& s : tok) {
142  result.push_back(boost::lexical_cast<size_t>(s));
143  }
144  // Note: the shape is in fortran order, so we need to reverse
145  std::reverse(result.begin(), result.end());
146  }
147  return result;
148 }
149 
150 std::vector<std::type_index> autoDetectColumnTypes(const CCfits::Table& table_hdu) {
151  std::vector<std::type_index> types {};
152  for (int i=1; i<=table_hdu.numCols(); i++) {
153  auto& column = table_hdu.column(i);
154 
155  if (typeid(table_hdu) == typeid(CCfits::BinTable)) {
156  column.setDimen();
157  types.push_back(binaryFormatToType(column.format(), parseTDIM(column.dimen())));
158  } else {
159  types.push_back(asciiFormatToType(column.format()));
160  }
161  }
162  return types;
163 }
164 
165 std::vector<std::string> autoDetectColumnUnits(const CCfits::Table& table_hdu) {
166  std::vector<std::string> units {};
167  for (int i=1; i<=table_hdu.numCols(); ++i) {
168  units.push_back(table_hdu.column(i).unit());
169  }
170  return units;
171 }
172 
173 std::vector<std::string> autoDetectColumnDescriptions(const CCfits::Table& table_hdu) {
174  std::vector<std::string> descriptions {};
175  for (int i=1; i<=table_hdu.numCols(); ++i) {
176  std::string desc;
177  auto key = table_hdu.keyWord().find("TDESC" + std::to_string(i));
178  if (key != table_hdu.keyWord().end()) {
179  key->second->value(desc);
180  }
181  descriptions.push_back(desc);
182  }
183  return descriptions;
184 }
185 
186 template<typename T>
187 std::vector<Row::cell_type> convertScalarColumn(CCfits::Column& column, long first, long last) {
188  std::vector<T> data;
189  column.read(data, first, last);
190  std::vector<Row::cell_type> result;
191  for (auto value : data) {
192  result.push_back(value);
193  }
194  return result;
195 }
196 
197 template<typename T>
198 std::vector<Row::cell_type> convertVectorColumn(CCfits::Column& column, long first, long last) {
199  std::vector<std::valarray<T>> data;
200  column.readArrays(data, first, last);
201  std::vector<Row::cell_type> result;
202  for (auto& valar : data) {
203  result.push_back(std::vector<T>(std::begin(valar),std::end(valar)));
204  }
205  return result;
206 }
207 
208 template<typename T>
209 std::vector<Row::cell_type> convertNdArrayColumn(CCfits::Column& column, long first, long last) {
210  std::vector<std::valarray<T>> data;
211  column.readArrays(data, first, last);
212  std::vector<size_t> shape = parseTDIM(column.dimen());
213 
214  std::vector<Row::cell_type> result;
215  for (auto& valar : data) {
216  result.push_back(NdArray<T>(shape, std::move(std::vector<T>(std::begin(valar), std::end(valar)))));
217  }
218  return result;
219 }
220 
221 std::vector<Row::cell_type> translateColumn(CCfits::Column& column, std::type_index type) {
222  return translateColumn(column, type, 1, column.rows());
223 }
224 
225 std::vector<Row::cell_type> translateColumn(CCfits::Column& column, std::type_index type, long first, long last) {
226  if (type == typeid(bool)) {
227  return convertScalarColumn<bool>(column, first, last);
228  } else if (type == typeid(int32_t)) {
229  return convertScalarColumn<int32_t>(column, first, last);
230  } else if (type == typeid(int64_t)) {
231  return convertScalarColumn<int64_t>(column, first, last);
232  } else if (type == typeid(float)) {
233  return convertScalarColumn<float>(column, first, last);
234  } else if (type == typeid(double)) {
235  return convertScalarColumn<double>(column, first, last);
236  } else if (type == typeid(std::string)) {
237  return convertScalarColumn<std::string>(column, first, last);
238  } else if (type == typeid(std::vector<int32_t>)) {
239  return convertVectorColumn<int32_t>(column, first, last);
240  } else if (type == typeid(std::vector<int64_t>)) {
241  return convertVectorColumn<int64_t>(column, first, last);
242  } else if (type == typeid(std::vector<float>)) {
243  return convertVectorColumn<float>(column, first, last);
244  } else if (type == typeid(std::vector<double>)) {
245  return convertVectorColumn<double>(column, first, last);
246  } else if (type == typeid(NdArray<int32_t>)) {
247  return convertNdArrayColumn<int32_t>(column, first, last);
248  } else if (type == typeid(NdArray<int64_t>)) {
249  return convertNdArrayColumn<int64_t>(column, first, last);
250  } else if (type == typeid(NdArray<float>)) {
251  return convertNdArrayColumn<float>(column, first, last);
252  } else if (type == typeid(NdArray<double>)) {
253  return convertNdArrayColumn<double>(column, first, last);
254  }
255  throw Elements::Exception() << "Unsupported column type " << type.name();
256 }
257 
258 }
259 } // end of namespace Euclid
constexpr double s
std::vector< Row::cell_type > convertVectorColumn(CCfits::Column &column, long first, long last)
std::type_index binaryFormatToType(const std::string &format, const std::vector< size_t > &shape)
NdArray(const std::vector< size_t > &shape)
Definition: NdArray.h:62
std::type_index asciiFormatToType(const std::string &format)
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)
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< Row::cell_type > convertScalarColumn(CCfits::Column &column, long first, long last)
std::vector< Row::cell_type > convertNdArrayColumn(CCfits::Column &column, long first, long last)
std::vector< std::string > autoDetectColumnUnits(const CCfits::Table &table_hdu)
Reads the column units based on the TUNITn keyword.
std::vector< std::type_index > autoDetectColumnTypes(const CCfits::Table &table_hdu)
Reads the column types of the given table HDU.
std::vector< std::string > autoDetectColumnNames(std::istream &in, const std::string &comment, size_t columns_number)
Reads the column names of the given stream.