2 * Copyright (C) 2012-2020 Euclid Science Ground Segment
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)
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
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
21 * @file FitsSerialize.icpp
26 #include <type_traits>
27 #include <boost/filesystem.hpp>
28 #include <CCfits/CCfits>
29 #include "XYDataset/QualifiedName.h"
30 #include "Table/Table.h"
31 #include "Table/FitsWriter.h"
32 #include "GridContainer/GridContainer.h"
33 #include "GridConstructionHelper.h"
36 namespace GridContainer {
39 struct FitsBpixTraits {
40 static_assert(!std::is_same<T,T>::value, "FITS arrays of type T are not supported");
44 struct FitsBpixTraits<std::int8_t> {
45 static constexpr int BPIX = BYTE_IMG;
49 struct FitsBpixTraits<std::int16_t> {
50 static constexpr int BPIX = SHORT_IMG;
54 struct FitsBpixTraits<std::int32_t> {
55 static constexpr int BPIX = LONG_IMG;
59 struct FitsBpixTraits<std::int64_t> {
60 static constexpr int BPIX = LONGLONG_IMG;
64 struct FitsBpixTraits<float> {
65 static constexpr int BPIX = FLOAT_IMG;
69 struct FitsBpixTraits<double> {
70 static constexpr int BPIX = DOUBLE_IMG;
74 struct GridAxisValueFitsHelper {
76 static FitsType axisToFits(const T& value) {
79 static T FitsToAxis(const FitsType& value) {
85 struct GridAxisValueFitsHelper<XYDataset::QualifiedName> {
86 using FitsType = std::string;
87 static FitsType axisToFits(const XYDataset::QualifiedName& value) {
88 return value.qualifiedName();
90 static XYDataset::QualifiedName FitsToAxis(const FitsType& value) {
95 template<typename... AxesTypes>
96 struct GridAxesToFitsHelper {
99 static void addGridAxesToFitsFile(const boost::filesystem::path& filename,
100 const std::string& array_hdu_name,
101 const std::tuple<GridAxis<AxesTypes>...>& axes_tuple,
102 const TemplateLoopCounter<I>&) {
103 addGridAxesToFitsFile(filename, array_hdu_name, axes_tuple, TemplateLoopCounter<I-1>{});
105 auto& axis = std::get<I-1>(axes_tuple);
106 using AxisType = typename std::remove_reference<decltype(axis)>::type::data_type;
107 using FitsType = typename GridAxisValueFitsHelper<AxisType>::FitsType;
109 std::vector<Table::ColumnInfo::info_type> info_list {
110 Table::ColumnInfo::info_type {"Index", typeid(int32_t)},
111 Table::ColumnInfo::info_type {"Value", typeid(FitsType)}
113 std::shared_ptr<Table::ColumnInfo> column_info {new Table::ColumnInfo {std::move(info_list)}};
115 std::vector<Table::Row> row_list {};
116 for (size_t i=0; i<axis.size(); ++i) {
117 auto fits_value = GridAxisValueFitsHelper<AxisType>::axisToFits(axis[i]);
118 row_list.push_back(Table::Row{{(int)i, fits_value}, column_info});
120 Table::Table table {row_list};
122 Table::FitsWriter{filename.string(), false}
123 .setFormat(Table::FitsWriter::Format::BINARY)
124 .setHduName(axis.name()+"_"+array_hdu_name)
128 static void addGridAxesToFitsFile(const boost::filesystem::path&,
130 const std::tuple<GridAxis<AxesTypes>...>&,
131 const TemplateLoopCounter<0>&) {
136 template<typename GridCellManager, typename... AxesTypes>
137 void gridFitsExport(const boost::filesystem::path& filename,
138 const std::string& hdu_name,
139 const GridContainer<GridCellManager, AxesTypes...>& grid) {
140 auto& axes = grid.getAxesTuple();
142 // Create the first HDU with the array. We do that in a scope so the file is
143 // created and the data are flushed into it before we continue.
145 CCfits::FITS fits (filename.string(), CCfits::Write);
147 auto ext_ax_size_t = GridConstructionHelper<AxesTypes...>::createAxesSizesVector(
148 axes, TemplateLoopCounter<sizeof...(AxesTypes)>{});
149 std::vector<long> ext_ax {ext_ax_size_t.begin(), ext_ax_size_t.end()};
151 using cell_type = typename GridCellManagerTraits<GridCellManager>::data_type;
152 auto bpix = FitsBpixTraits<cell_type>::BPIX;
153 fits.addImage(hdu_name, bpix, ext_ax);
154 std::valarray<cell_type> data (grid.size());
156 for (auto value : grid) {
160 fits.currentExtension().write(1, grid.size(), data);
163 GridAxesToFitsHelper<AxesTypes...>::addGridAxesToFitsFile(filename, hdu_name,
164 axes, TemplateLoopCounter<sizeof...(AxesTypes)>{});
167 template <typename GridType>
168 class GridAxisFitsReader {
171 using AxisType = typename std::remove_reference<decltype(std::declval<GridType>().template getAxis<I>())>::type;
173 template <int I, typename=void>
174 struct AxesTupleType {
175 using previous = typename AxesTupleType<I-1>::type;
176 using type = decltype(std::tuple_cat(std::declval<previous>(), std::declval<std::tuple<AxisType<I>>>()));
180 struct AxesTupleType<I, typename std::enable_if<I == -1>::type> {
181 using type = std::tuple<>;
185 using GridAxisType = typename std::remove_reference<decltype(std::declval<GridType>().template getAxis<I>())>::type;
188 static GridAxisType<I> readAxis(const std::string& grid_name, CCfits::ExtHDU& hdu) {
189 using KnotType = typename GridAxisType<I>::data_type;
190 using FitsType = typename GridAxisValueFitsHelper<KnotType>::FitsType;
192 auto axis_name = hdu.name().substr(0, hdu.name().size() - grid_name.size() - 1);
193 std::vector<FitsType> data {};
195 auto& column = hdu.column("Value");
196 column.read(data, 1, column.rows());
197 } catch (CCfits::FitsException e) {
198 throw Elements::Exception() << e.message();
200 std::vector<KnotType> knots {};
201 for (std::size_t i = 0; i < data.size(); ++i) {
202 knots.emplace_back(GridAxisValueFitsHelper<KnotType>::FitsToAxis(data[i]));
204 return {std::move(axis_name), std::move(knots)};
208 static typename AxesTupleType<I>::type readAxesTuple(CCfits::FITS& fits, const std::string& grid_name, int hdu_index, const TemplateLoopCounter<I>&) {
209 auto axis = readAxis<I>(grid_name, fits.extension(hdu_index));
210 auto previous = readAxesTuple(fits, grid_name, hdu_index-1, TemplateLoopCounter<I-1>{});
211 return std::tuple_cat(std::move(previous), std::tuple<decltype(axis)>{std::move(axis)});
214 static std::tuple<> readAxesTuple(CCfits::FITS&, const std::string&, int, const TemplateLoopCounter<-1>&) {
220 static typename AxesTupleType<GridType::axisNumber()-1>::type readAllAxes(CCfits::FITS& fits, int hdu_index) {
221 auto name = fits.extension(hdu_index).name();
222 return readAxesTuple(fits, name, hdu_index + GridType::axisNumber(), TemplateLoopCounter<GridType::axisNumber()-1>{});
227 template<typename GridType>
228 GridType gridFitsImport(const boost::filesystem::path& filename, int hdu_index) {
229 CCfits::FITS fits (filename.string(), CCfits::Read);
231 auto axes = GridAxisFitsReader<GridType>::readAllAxes(fits, hdu_index);
233 GridType grid {std::move(axes)};
235 std::valarray<typename GridType::cell_type> data {};
236 fits.extension(hdu_index).read(data);
239 for (auto iter = grid.begin(); iter != grid.end(); ++iter, ++i) {
246 } // end of namespace GridContainer
247 } // end of namespace Euclid