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