3D NeXus writing

main.cpp

//---------------------------------------------------------------------------
// Includes
//---------------------------------------------------------------------------
#include <cmath>
#include <iostream>
#include <unistd.h>
#include <vector>

#include <nexus/NeXusFile.hpp>

using namespace NeXus;
using std::size_t;

//---------------------------------------------------------------------------
// Typedefs
//---------------------------------------------------------------------------
typedef std::vector<double> WorkspaceRow;
typedef std::vector<WorkspaceRow> Workspace;

//---------------------------------------------------------------------------
// Forward declarations
//---------------------------------------------------------------------------
Workspace * createWorkspace(const size_t nfiles,
                            const size_t npixel_rows,
                            const size_t npixel_cols);
void saveToFile(const std::string & filename,
                const Workspace & workspace);
File * createFile(const std::string & filename);

//---------------------------------------------------------------------------
// Definitions
//---------------------------------------------------------------------------

int main(int, char**)
{
  const int nfiles(6), npixel_rows(64), npixel_cols(64);
  auto * workspace = createWorkspace(nfiles, npixel_rows,
                                     npixel_cols);
  saveToFile("test.nxs", *workspace);
  delete workspace;
}

/**
 * Create a 2D data matrix with dimensions
 * NX=nfiles
 * NY=npixel_cols*npixel_rows
 */
Workspace * createWorkspace(const size_t nfiles,
                            const size_t npixel_rows,
                            const size_t npixel_cols)
{
  size_t nrows = npixel_rows*npixel_cols;
  auto * workspace = new Workspace(nrows, WorkspaceRow(nfiles));
  // Pattern
  // 1       nrows+1
  // 2       nrows+2
  // ..      ..
  // nrows   2*rows
  for(size_t i = 0; i < nrows; ++i)
  {
    auto & row = (*workspace)[i];
    for(size_t j = 0; j < nfiles; ++j)
    {
      row[j] = j*nrows + i + 1;
    }
  }
  return workspace;
}

/**
 * Save to a NeXus file with a 3D structure:
 *   nfiles x (npixel_rows x npixel_cols)
 */
void saveToFile(const std::string & filename,
                const Workspace & workspace)
{
  auto *nxfile = createFile(filename);

  // Each column in the workspace is a separate file of (npixel_rows x npixel_cols) values
  // Write to the file as a 3D array where the first index gives a
  // 2D (npixel_rows x npixel_cols) array for each image
  size_t npixels(workspace.size());
  size_t nfiles(workspace[0].size());
  size_t npix_row(static_cast<size_t>(sqrt(npixels))); // we know it's square
  size_t npix_col(npix_row);
  std::vector<int> dims(3);
  dims[0] = nfiles;
  dims[1] = npix_row;
  dims[2] = npix_col;
  bool autoOpen(true);
  nxfile->makeData("data1", FLOAT64, dims, autoOpen);

  // Buffer for writing a row from the workspace at a time
  // A workspace row corresponds to a point in each file
  // so a 3D buffer that is 1D in the block dimensions will do
  std::vector<int> slabSize(3, 0), slabStart(3, 0);
  slabSize[0] = nfiles;
  slabSize[1] = 1; // fixed to 1 while reading a single row at a time
  slabSize[2] = 1; // fixed to 1 while reading a single row at a time
  double *buffer = new double[nfiles];

  // Write to file
  for(size_t detrow = 0; detrow < npixels; ++detrow) // workspace rows (detectors)
  {
    const auto & workspaceRow = workspace[detrow]; // workspace row reference
    // map workspace row to 2D subblock row/col indices
    const size_t blockRow = detrow / npix_row;
    const size_t blockCol = detrow % npix_col;
    // slab offsets
    slabStart[0] = 0; // each write puts whole length in file dimension
    slabStart[1] = blockRow; // start in 2d row
    slabStart[2] = blockCol; // start in 2d column
    // fill buffer
    for(size_t i = 0; i < nfiles; ++i)
    {
      buffer[i] = workspaceRow[i];
    }
    // push to file
    nxfile->putSlab(buffer, slabStart, slabSize);
  }
  nxfile->closeData();
  delete [] buffer;

  nxfile->closeGroup();
  nxfile->closeGroup();
  delete nxfile;
}

/**
 * Create new file hierchy (possibly deleting the old one)
 * Leaves the data open
 */
File * createFile(const std::string & filename)
{
   // remove existing
  unlink(filename.c_str());
  // Create file and structure
  auto *nxfile = new File(filename, NXACC_CREATE5);
  const bool autoOpen(true);
  nxfile->makeGroup("workspace", "NXentry", autoOpen);
  nxfile->makeGroup("data", "NXdata", autoOpen);

  return nxfile;
}