EleFits  4.0.0
A modern C++ API on top of CFitsIO
Tutorial

Introduction

In this tutorial, we will show how to read and write multi-extension Fits (MEF) files. This means understanding the usage of the following service classes:

as well as Record, Raster and Column data classes.

We strongly recommend reading first What's a Fits file?, even if you're familiar with the Fits format, because it introduces EleFits-specific wording.

At the end of the tutorial, you will be able to create a MEF file from scratch with unlimited number of various extensions, and to read the metadata and data back!

The tutorial is built together with an example program: EleFitsTutorial.cpp. We've embedded the calls to the logger in the code snippets below, so that you can easily map the execution log to the following explanations. You can already run the program and watch the resulting file:

EleFitsTutorial --output tuto.fits
fv tuto.fits

We'll first discover the data classes, then use them to create a MEF file, and finally read the file and values back.

See also
To go further, at the end of each section, a link points to the reference page for the introduced topic.

Setup

First things first, we have to declare the dependency to EleFits and to use the right headers and namespaces. For the first part, head to the install_guide. For the headers and namespace, here's the only thing you'll have to do:

#include "EleFits/MefFile.h"
using namespace Euclid;
// EleFits API is in the Euclid::Fits namespace.
// We could have be using namespace Euclid::Fits instead,
// but things would have been less obvious in the snippets.

Data classes

Data classes are the classes which hold the pieces of data read from and written to a Fits file. There are three main classes:

  • Record for header units,
  • Raster for data units of image HDUs,
  • Column for data units of binary table HDUs.

Record

A keyword record as defined in the Fits standard is a triplet of keyword, value and optional comment rendered in the file as:

KEYWORD = VALUE / comment

A unit can be included in the comment as follows:

KEYWORD = VALUE / [unit] comment

The value can be a Boolean, an integer, real or complex number, or a string.

Record is a simple template class which merely stores those fields. For the purpose of the tutorial, we define a structure to hold our records:

struct TutoRecords {
Fits::Record<std::string> stringRecord;
Fits::Record<int> intRecord;
Fits::Record<float> floatRecord;
Fits::Record<std::complex<double>> complexRecord;
};

Here's a how they can be built:

/* Create a record with unit and comment */
Fits::Record<std::string> stringRecord("STRING", "VALUE", "unit", "comment");
/* Create a record with keyword and value only */
Fits::Record<int> intRecord("INT", 0);
/* Create a record from an initialization list */
Fits::Record<float> floatRecord { "FLOAT", 3.14F, "", "A piece of Pi" };
// This is often used as a shortcut to create records as function parameters.
/* Generate a random record */
auto complexRecord = Fits::Test::generateRandomRecord<std::complex<double>>("COMPLEX");
See also
Header data classes

Raster

Images in Fits are n-dimensional arrays. Raster is a class which represents this kind of data and provides constant time pixel accessors. Template arguments are the pixel type and number of axes. There are two kinds of rasters:

  • VecRaster owns the data as an std::vector;
  • PtrRaster points to the data owned by another class, as a raw pointer (be careful not to destroy the data while the raster is alive).

The rasters of this tutorial are stored in a dedicated structure:

struct TutoRasters {
Fits::VecRaster<std::int16_t, 2> int16Raster2D;
Fits::VecRaster<std::int32_t, 3> int32Raster3D;
Fits::VecRaster<std::int64_t, 4> int64Raster4D;
};

Again, let's show an example of how to create rasters:

/* Initialize and later fill a raster */
Fits::VecRaster<std::int16_t, 2> int16Raster2D({ 4, 3 });
for (const auto& position : int16Raster2D.domain()) {
int16Raster2D[position] = position[0] + position[1];
}
// This demonstrates the iteration over positions;
// It is possible to use two nested loops instead.
/* Create a raster from a vector */
std::vector<std::int32_t> int32Vec(16 * 9 * 3, 0);
// ... do what you have to do with the vector, and then move it to the raster ...
Fits::VecRaster<std::int32_t, 3> int32Raster3D({ 16, 9, 3 }, std::move(int32Vec));
// Instead of moving a vector, it's also possible to work with
// a raw pointer with the PtrRaster class.
/* Generate a random raster */
auto int64Raster4D = Fits::Test::RandomRaster<std::int64_t, 4>({ 17, 9, 3, 24 });
See also
Image data classes

Column

A binary table is made of columns which can be either scalar (each cell contains a value) or vector (each cell contains several values). In the latter case, the number of values per cell is named repeat count.

For string columns, the repeat count must be greater than the longest string value.

A Column is a simple structure which holds the column name, repeat count (=1 for scalar columns, >1 for vector columns), unit, and owns or references some data, respectively as a VecColumn or as a PtrColumn.

The following structure stores the tutorial columns:

struct TutoColumns {
Fits::VecColumn<std::string> stringColumn;
Fits::VecColumn<std::int32_t> int32Column;
Fits::VecColumn<float> float32Column;
};

Here's a set of examples on how create them:

/* Initialize and later fill a column */
Fits::VecColumn<std::string> stringColumn({ "STRING", "unit", 3 }, 100);
// String columns must be wide-enough to hold each character.
for (long i = 0; i < stringColumn.rowCount(); ++i) {
stringColumn(i) = std::to_string(i);
}
// operator() takes two parameters: the row index, and repeat index (=0 by default)
/* Create a column from a vector */
// ... do what you have to do with the vector, and then move it to the column ...
Fits::VecColumn<std::int32_t> int32Column({ "INT32", "", 1 }, std::move(int32Vec));
// Analogously to rasters, columns can be managed with the lightweight PtrColumn classe.
/* Generate a random column */
auto float32Column = Fits::Test::RandomVectorColumn<float>(8, 100);
See also
Binary table data classes

Open (and close) a MEF file

Now that we have learned about the data classes, we can move to the service classes, or handlers. The first thing we can do is to create a new file.

The MefFile class represents multi-extension Fits (MEF) files. It extends the FitsFile class, which is also the base class of SifFile for single-image Fits (SIF) files.

Creating (or opening) a file is simply done with the constructor of MefFile (or SifFile for SIF files!):

Fits::MefFile f(filename, Fits::FileMode::Create);

A newly created Fits file consists in an empty Primary, which can then be accessed and modified, but is never created by hand.

The mode parameter controls the access rights on the file. For example, an existing file can be opened with read-only mode:

Fits::MefFile f(filename, Fits::FileMode::Read);

The file is closed when the destructor of MefFile is called (although a FitsFile::close() method is provided for convenience).

See also
File handlers

Write a MEF file

Create extensions

MefFile provides services to create and access extensions. Conceptually, MefFile is a FitsFile with a vector of HDUs (in contrast, SifFile is a FitsFile with a single HDU).

HDUs are represented by two classes:

Both derive from parent class Hdu.

To sum up, a MefFile is a vector of Hdus, which can be a mix of ImageHdus and BintableHdus.

There are two kinds of services in MefFile for creating extensions: They can be either initialized with header only or assigned directly with data. Here's an example of creating image HDUs:

/* Fill the header and data units */
const auto& image1 = f.assignImageExt("IMAGE1", rasters.int32Raster3D);
/* Fill the header only (for now) */
const auto& image2 = f.initImageExt<std::int16_t>("IMAGE2", rasters.int16Raster2D.shape());

And here's one of creating binary table HDUs:

/* Fill the header and data units */
const auto& table1 = f.assignBintableExt("TABLE1", columns.stringColumn, columns.int32Column, columns.float32Column);
/* Fill the header unit only (for now) */
const auto& table2 = f.initBintableExt(
"TABLE2",
columns.stringColumn.info(),
columns.int32Column.info(),
columns.float32Column.info());
See also
File handlers

Write records

Records are read and written through a Header object, which is accessible from an ImageHdu or a BintableHdu as header(). They can be written and updated one-by-one or by sequences.

An optional template parameter controls the write policy, e.g. what to do if the given keyword already exists. By default, a record is created if the keyword doesn't exists, or updated if the keyword already exists.

Record writing functions have different overloards to allow calling them in different ways. Here are a few examples:

/* Write a single record */
hdu.header().write(records.stringRecord);
/* Write several records */
hdu.header().writeSeq(records.intRecord, records.floatRecord, records.complexRecord);
/* Update using initialization lists */
Fits::Record<int>("INT", 1),
Fits::Record<float>("FLOAT", 3.14159F, "", "A larger piece of Pi"),
Fits::Record<std::complex<double>>("COMPLEX", { 180., 90. }));
See also
Header unit handlers

Write image data

Remember that we have left IMAGE2 extension without data? Filling it is done through the ImageRaster class, accessed as ImageHdu::raster(). For example, here is how to write all the pixels at once:

image2.raster().write(rasters.int16Raster2D);

Many more options are available to write data region-wise. They are described in details in Image data unit handlers.

See also
Image data unit handlers

Write binary table data

Analogously to ImageRaster for image HDUs, BintableColumns provides read/write services for the data unit of the binary table extension.

As the name implies, data is stored, read and written column-wise, as opposed to the row-major ordering of the Fits file. This is to avoid working with very complex structures like tuples with holes, and rather rely on Column, std::vector, arrays... Yet, this also means that I/Os could be very inefficient! To workaround this, several columns can be written (and read) at a time. In this case, I/Os are internally performed chunk-wise, using some buffer, and performances drammatically improve.

/* Write a single column */
table2.columns().write(columns.stringColumn);
/* Write several columns */
table2.columns().writeSeq(columns.int32Column, columns.float32Column);

Like for image HDUs, it is possible to write the data unit only partially, e.g. to append rows, using richer methods of BintableColumns.

See also
Binary table data unit handlers

Read a MEF file

Access HDUs

HDUs can be accessed with a set of methods, templated with the type of HDU: BintableHdu, ImageHdu, or Hdu. For metadata work, we don't need to know the type of HDU: whether this is an image or binary table HDU has no impact, and a Hdu will be returned by default.

HDUs are accessed either by their name (first HDU whose name matches argument is returned) or by their index, and a shortcut is provided for the primary HDU (which has index 0):

/* Access the Primary HDU */
const auto& primary = f.primary();
const auto primaryIndex = primary.index();
// Indices are 0-based.
/* Access an HDU by its index */
const auto& image2 = f.access<Fits::ImageHdu>(2);
const auto imageName = image2.readName();
/* Access an HDU by its name */
const auto& table1 = f.accessFirst<Fits::BintableHdu>("TABLE1");
const auto tableIndex = table1.index();
// If several HDUs have the same name, the first one is returned.

You've probably noticed that we use references to constant HDU handlers. Indeed, HDU handlers are not modified by reading and writing services, only the MefFile is. Actually, HDU-level and data unit-level handlers are very light classes which aim at providing a clear interface.

See also
File handlers

Parse records

Records are parsed using Header services. Like for writing, you can parse several records at once:

/* Read a single record */
const auto intRecord = hdu.header().parse<int>("INT");
// Records can be sliced as their value for immediate use:
const int intValue = hdu.header().parse<int>("INT");
/* Read several records */
const auto someRecords = hdu.header().parseSeq(
Fits::Named<std::string>("STRING"),
Fits::Named<int>("INT"),
Fits::Named<float>("FLOAT"),
Fits::Named<std::complex<double>>("COMPLEX"));
const auto& thirdRecord = std::get<2>(someRecords);
/* Read as VariantValue */
const auto variantRecords = hdu.header().parseSeq<>({ "INT", "COMPLEX" });
const auto complexRecord = variantRecords.as<std::complex<double>>("COMPLEX");
/* Read as a user-defined structure */
const auto tutoRecords = hdu.header().parseStruct<TutoRecords>(
Fits::Named<std::string>("STRING"),
Fits::Named<int>("INT"),
Fits::Named<float>("FLOAT"),
Fits::Named<std::complex<double>>("COMPLEX"));
const auto& stringRecord = tutoRecords.stringRecord;
See also
Header unit handlers

Read image data

Like for writing, reading images is ensured by ImageRaster. A VecRaster is generally instantiated, and pixel values can be accessed with the subscript operator []:

const auto image = hdu.raster().read<std::int16_t, 2>();
const auto& firstPixel = image[{ 0, 0 }];
const auto& lastPixel = image.at({ -1, -1 });
// `operator[]` performs no bound checking, while `at` does and enables backward indexing.

Again, regions can be read.

See also
Image data unit handlers

Read binary table data

Column reading is provided by BintableColumns. Columns are generally read as VecColumns, values of which are accessed with the call operator ():

/* Read a single column */
const auto vectorColumn = hdu.readColumn<double>("VECTOR");
/* Read several columns by their name */
const auto byName = hdu.columns().readSeq(Fits::Named<std::string>("STRING"), Fits::Named<std::int32_t>("INT32"));
const auto& stringColumn = std::get<0>(byName);
/* Read several columns by their index */
const auto byIndex = hdu.columns().readSeq(Fits::Indexed<std::string>(0), Fits::Indexed<std::int32_t>(1));
const auto& intColumn = std::get<1>(byIndex);
/* Use values */
const auto& firstString = stringColumn(0);
const auto& firstInt = intColumn(0);
const auto& lastFloat = vectorColumn.at(-1, -1);
// There is no operator[]() for columns, because vector columns require 2 indices (row and repeat).
// operator()() performs no bound checking, while at() does and enables backward indexing.
See also
Binary table data unit handlers

Wrap up

A picture is worth a thousand words; below is the class diagram of what we've just learned (i.e. the main contents of the Euclid::Fits namespace).

To go further or out of curiosity, head to the module pages or related pages. Specifically, data classes and service classes are described in more details in Data classes and File and HDU handlers. To practice or test your own code, you can also browse the Euclid::Fits::Test namespace, which provides ready-to-use functions and classes, such as the random records, rasters and columns we've seen above.

std::move
T move(T... args)
std::vector
STL class.
Euclid::Fits::BintableHdu::columns
const BintableColumns & columns() const
Access the data unit column-wise.
Euclid::Fits::RecordMode::CreateOrUpdate
@ CreateOrUpdate
Modify a record if keyword already exists, create a record otherwise.
Euclid::Fits::Header::parse
Record< T > parse(const std::string &keyword) const
Parse a record.
Euclid::Fits::Header::parseSeq
RecordVec< T > parseSeq(const std::vector< std::string > &keywords) const
Parse a sequence of homogeneous records.
Euclid::Fits::Hdu::index
long index() const
Get the 0-based index of the HDU.
Euclid::Fits::Hdu::readName
std::string readName() const
Read the extension name.
std::to_string
T to_string(T... args)
std::int16_t
Euclid::Fits::Header::writeSeq
void writeSeq(const Record< Ts > &... records) const
Write a homogeneous sequence of records.
Euclid::Fits::Header::write
void write(const Record< T > &record) const
Write a record.
Euclid::Fits::Header::parseStruct
TOut parseStruct(const Named< Ts > &... keywords) const
Parse a sequence of records.
std::complex< double >
Euclid::Fits::BintableHdu::readColumn
VecColumn< T > readColumn(long index) const
Read a column with given index.
Euclid::Fits::Hdu::header
const Header & header() const
Access the header unit to read and write records.
Euclid
Euclid SGS namespace.
Definition: cfitsio_toc.dox:1
Euclid::Fits::BintableColumns::readSeq
std::tuple< VecColumn< Ts >... > readSeq(const Named< Ts > &... names) const
Read the columns with given names.