Alexandria  2.27.0
SDC-CH common library for the Euclid project
FitsSerialize.icpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012-2021 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 
19 /*
20  * @file FitsSerialize.icpp
21  * @author nikoapos
22  */
23 
24 #include "GridConstructionHelper.h"
25 #include "GridContainer/GridContainer.h"
26 #include "Table/FitsWriter.h"
27 #include "Table/Table.h"
28 #include "XYDataset/QualifiedName.h"
29 #include <CCfits/CCfits>
30 #include <boost/filesystem.hpp>
31 #include <type_traits>
32 #include <valarray>
33 
34 namespace Euclid {
35 namespace GridContainer {
36 
37 template <typename T>
38 struct FitsBpixTraits {
39  static_assert(!std::is_same<T, T>::value, "FITS arrays of type T are not supported");
40 };
41 
42 template <>
43 struct FitsBpixTraits<std::int8_t> {
44  static constexpr int BPIX = BYTE_IMG;
45 };
46 
47 template <>
48 struct FitsBpixTraits<std::int16_t> {
49  static constexpr int BPIX = SHORT_IMG;
50 };
51 
52 template <>
53 struct FitsBpixTraits<std::int32_t> {
54  static constexpr int BPIX = LONG_IMG;
55 };
56 
57 template <>
58 struct FitsBpixTraits<std::int64_t> {
59  static constexpr int BPIX = LONGLONG_IMG;
60 };
61 
62 template <>
63 struct FitsBpixTraits<float> {
64  static constexpr int BPIX = FLOAT_IMG;
65 };
66 
67 template <>
68 struct FitsBpixTraits<double> {
69  static constexpr int BPIX = DOUBLE_IMG;
70 };
71 
72 template <typename T>
73 struct GridAxisValueFitsHelper {
74  using FitsType = T;
75  static FitsType axisToFits(const T& value) {
76  return value;
77  }
78  static T FitsToAxis(const FitsType& value) {
79  return value;
80  }
81 };
82 
83 template <>
84 struct GridAxisValueFitsHelper<XYDataset::QualifiedName> {
85  using FitsType = std::string;
86  static FitsType axisToFits(const XYDataset::QualifiedName& value) {
87  return value.qualifiedName();
88  }
89  static XYDataset::QualifiedName FitsToAxis(const FitsType& value) {
90  return value;
91  }
92 };
93 
94 template <typename... AxesTypes>
95 struct GridAxesToFitsHelper {
96 
97  template <int I>
98  static void addGridAxesToFitsFile(const boost::filesystem::path& filename, const std::string& array_hdu_name,
99  const std::tuple<GridAxis<AxesTypes>...>& axes_tuple, const TemplateLoopCounter<I>&) {
100  addGridAxesToFitsFile(filename, array_hdu_name, axes_tuple, TemplateLoopCounter<I - 1>{});
101 
102  auto& axis = std::get<I - 1>(axes_tuple);
103  using AxisType = typename std::remove_reference<decltype(axis)>::type::data_type;
104  using FitsType = typename GridAxisValueFitsHelper<AxisType>::FitsType;
105 
106  std::vector<Table::ColumnInfo::info_type> info_list{Table::ColumnInfo::info_type{"Index", typeid(int32_t)},
107  Table::ColumnInfo::info_type{"Value", typeid(FitsType)}};
108  std::shared_ptr<Table::ColumnInfo> column_info{new Table::ColumnInfo{std::move(info_list)}};
109 
110  std::vector<Table::Row> row_list{};
111  for (size_t i = 0; i < axis.size(); ++i) {
112  auto fits_value = GridAxisValueFitsHelper<AxisType>::axisToFits(axis[i]);
113  row_list.push_back(Table::Row{{(int)i, fits_value}, column_info});
114  }
115  Table::Table table{row_list};
116 
117  Table::FitsWriter{filename.string(), false}
118  .setFormat(Table::FitsWriter::Format::BINARY)
119  .setHduName(axis.name() + "_" + array_hdu_name)
120  .addData(table);
121  }
122 
123  static void addGridAxesToFitsFile(const boost::filesystem::path&, const std::string&, const std::tuple<GridAxis<AxesTypes>...>&,
124  const TemplateLoopCounter<0>&) {}
125 };
126 
127 template <typename GridCellManager, typename... AxesTypes>
128 void gridFitsExport(const boost::filesystem::path& filename, const std::string& hdu_name,
129  const GridContainer<GridCellManager, AxesTypes...>& grid) {
130  auto& axes = grid.getAxesTuple();
131 
132  // Create the first HDU with the array. We do that in a scope so the file is
133  // created and the data are flushed into it before we continue.
134  {
135  CCfits::FITS fits(filename.string(), CCfits::Write);
136 
137  auto ext_ax_size_t =
138  GridConstructionHelper<AxesTypes...>::createAxesSizesVector(axes, TemplateLoopCounter<sizeof...(AxesTypes)>{});
139  std::vector<long> ext_ax{ext_ax_size_t.begin(), ext_ax_size_t.end()};
140 
141  using cell_type = typename GridCellManagerTraits<GridCellManager>::data_type;
142  auto bpix = FitsBpixTraits<cell_type>::BPIX;
143  fits.addImage(hdu_name, bpix, ext_ax);
144  std::valarray<cell_type> data(grid.size());
145  int i = 0;
146  for (auto value : grid) {
147  data[i] = value;
148  ++i;
149  }
150  fits.currentExtension().write(1, grid.size(), data);
151  }
152 
153  GridAxesToFitsHelper<AxesTypes...>::addGridAxesToFitsFile(filename, hdu_name, axes, TemplateLoopCounter<sizeof...(AxesTypes)>{});
154 }
155 
156 template <typename GridType>
157 class GridAxisFitsReader {
158 
159  template <int I>
160  using AxisType = typename std::remove_reference<decltype(std::declval<GridType>().template getAxis<I>())>::type;
161 
162  template <int I, typename = void>
163  struct AxesTupleType {
164  using previous = typename AxesTupleType<I - 1>::type;
165  using type = decltype(std::tuple_cat(std::declval<previous>(), std::declval<std::tuple<AxisType<I>>>()));
166  };
167 
168  template <int I>
169  struct AxesTupleType<I, typename std::enable_if<I == -1>::type> {
170  using type = std::tuple<>;
171  };
172 
173  template <int I>
174  using GridAxisType = typename std::remove_reference<decltype(std::declval<GridType>().template getAxis<I>())>::type;
175 
176  template <int I>
177  static GridAxisType<I> readAxis(const std::string& grid_name, CCfits::ExtHDU& hdu) {
178  using KnotType = typename GridAxisType<I>::data_type;
179  using FitsType = typename GridAxisValueFitsHelper<KnotType>::FitsType;
180 
181  auto axis_name = hdu.name().substr(0, hdu.name().size() - grid_name.size() - 1);
182  std::vector<FitsType> data{};
183  try {
184  auto& column = hdu.column("Value");
185  column.read(data, 1, column.rows());
186  } catch (CCfits::FitsException e) {
187  throw Elements::Exception() << e.message();
188  }
189  std::vector<KnotType> knots{};
190  for (std::size_t i = 0; i < data.size(); ++i) {
191  knots.emplace_back(GridAxisValueFitsHelper<KnotType>::FitsToAxis(data[i]));
192  }
193  return {std::move(axis_name), std::move(knots)};
194  }
195 
196  template <int I>
197  static typename AxesTupleType<I>::type readAxesTuple(CCfits::FITS& fits, const std::string& grid_name, int hdu_index,
198  const TemplateLoopCounter<I>&) {
199  auto axis = readAxis<I>(grid_name, fits.extension(hdu_index));
200  auto previous = readAxesTuple(fits, grid_name, hdu_index - 1, TemplateLoopCounter<I - 1>{});
201  return std::tuple_cat(std::move(previous), std::tuple<decltype(axis)>{std::move(axis)});
202  }
203 
204  static std::tuple<> readAxesTuple(CCfits::FITS&, const std::string&, int, const TemplateLoopCounter<-1>&) {
205  return {};
206  }
207 
208 public:
209  static typename AxesTupleType<GridType::axisNumber() - 1>::type readAllAxes(CCfits::FITS& fits, int hdu_index) {
210  auto name = fits.extension(hdu_index).name();
211  return readAxesTuple(fits, name, hdu_index + GridType::axisNumber(), TemplateLoopCounter<GridType::axisNumber() - 1>{});
212  }
213 };
214 
215 template <typename GridType>
216 GridType gridFitsImport(const boost::filesystem::path& filename, int hdu_index) {
217  CCfits::FITS fits(filename.string(), CCfits::Read);
218 
219  auto axes = GridAxisFitsReader<GridType>::readAllAxes(fits, hdu_index);
220 
221  GridType grid{std::move(axes)};
222 
223  std::valarray<typename GridType::cell_type> data{};
224  fits.extension(hdu_index).read(data);
225 
226  int i = 0;
227  for (auto iter = grid.begin(); iter != grid.end(); ++iter, ++i) {
228  *iter = data[i];
229  }
230 
231  return grid;
232 }
233 
234 } // end of namespace GridContainer
235 } // end of namespace Euclid